a_curve_lin.cxx

Go to the documentation of this file.

/*
         Copyright 2002-2019 Pierre SMARS (smars@yuntech.edu.tw)
         This program is free software: you can redistribute it and/or modify
         it under the terms of the GNU General Public License as published by
         the Free Software Foundation, either version 2 of the License, or
         (at your option) any later version.
 
         This program is distributed in the hope that it will be useful,
         but WITHOUT ANY WARRANTY; without even the implied warranty of
         MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
         GNU General Public License for more details.
 
         You should have received a copy of the GNU General Public License
         along with this program.  If not, see <http://www.gnu.org/licenses/>.
         */
#include "a_curve_lin.h"
#include <math.h>
 
// .PORTABILITY : ansi C++
 
//---------------------------------------------------------------------------
a_curve_lin::a_curve_lin(double da) : a_curve(da)
{
}
//---------------------------------------------------------------------------
a_curve_lin::a_curve_lin(const a_curve_lin& ifn) :      a_curve(ifn)
{
        ltot_ = ifn.ltot();
        for (int i = 0; i < ifn.size(); i++)
                a_.push_back(ifn.a()[i]);
}
//---------------------------------------------------------------------------
a_curve_lin::~a_curve_lin()
{
        a_.clear();
}
//---------------------------------------------------------------------------
a_curve_lin& a_curve_lin::flatten(a_curve_lin& path)
{
        int n_pt_profile = this->size();
        //find the position of a mean plane through the profile along the path
 
        a_point xg;
        for (int i = 0; i < n_pt_profile; i++)
                xg += *x_[i];
        xg /= n_pt_profile;
 
        int ref;
        double m;
        double dist = path.dist2(xg,ref,m);//find the closest reference along the path
 
        a_point x0 = *path[ref]+ m*(*path[ref+1]-*path[ref]);//reference point on the path
        double a = path.a()[ref]+m*(path.a()[ref+1]-path.a()[ref]);//relative position on the path
 
 
        /*      int reft = 0;
                        for (int i = 0; i < n_pt_profile; i++)
                        {
                        a_point * x = x_[i];//for each pt of the profile
                        int ref;
                        double dist = path.dist(*x,ref);//find the closest reference along the path
                        reft += ref;
        //              std::cerr << ref << " ";
        }
        int refm = double(reft)/n_pt_profile+.5;//take the mean
        //      std::cerr << std::endl << refm << std::endl;
 
        a_point x0 = *path[refm];//reference point on the path
        double a = path.a()[refm];//relative position on the path
        */
 
 
        a_point centre = path.curvature_centre();
 
        //      a_point centre = path.curvature_centre(a);
        //define the profile plane
 
        a_point y_axis = (x0-centre).normalise();
        a_point z_dir = path.dz_axis(a);//not perpendicular to y_axis, used to define x_axis
        a_point x_axis = cross(y_axis,z_dir).normalise();//the axis of rotation
        a_point z_axis = cross(x_axis,y_axis).normalise();
        //project the profile's points
        std::vector<a_point> points_n;
        for (int i = 0; i < n_pt_profile; i++)
        {
                a_point p = *x_[i]-centre;//the vector to rotate
                a_point px = (p*x_axis)*x_axis;//component along axis do not change
                a_point pyz = ((p*y_axis)*y_axis)+((p*z_axis)*z_axis);//component perpendicular will change
                a_point pyz_r = (pyz.norm())*y_axis;//component rotated
                a_point xn = centre+pyz_r+px;//new vector
                points_n.push_back(xn);
        }
        for (int i = 0; i < n_pt_profile; i++)
                *x_[i] = points_n[i];
        points_n.clear();
        this->init();
        return *this;
}
//---------------------------------------------------------------------------
a_curve& a_curve_lin::init()
{
        if (!initialised_)
        {
                a_curve::init();
                double l = 0;
                a_.resize(n_pts_);
                a_[0] = 0.;
                a_point * p0 = x_[0];
                for (int i = 1; i < n_pts_; i++)
                {
                        a_point * p = x_[i];
                        l += (*p-*p0).norm();
                        a_[i] = l;
                        p0 = p;
                }
                ltot_ = l;
                for (int i = 1; i < n_pts_; i++)
                        a_[i] /= ltot_;
        }
        return *this;
}
//---------------------------------------------------------------------------
int a_curve_lin::ref(double a) const
{
        int i = -1;
        double a2;
        do
                a2 = a_[++i];
        while (a2 <= a);
        return i-1;
}
//---------------------------------------------------------------------------
a_point a_curve_lin::operator()(double a) const
{
        if (a<=0)
                return *x_[0];
        if (a>=1)
                return **(x_.end()-1);
        int ref = this->ref(a);
        a_point x1 = *x_[ref];
        a_point x2 = *x_[ref+1];
        double a1 = a_[ref];
        double a2 = a_[ref+1];
        double f = (a-a1)/(a2-a1);
        return (x1 + (x2-x1)*f);
}
//---------------------------------------------------------------------------
a_point a_curve_lin::dx_axis(double a) const
{
        a_point z = this->dz_axis(a);
        a_point y = this->dy_axis(a);
        return cross(y,z).normalise();
}
//---------------------------------------------------------------------------
a_point a_curve_lin::dy_axis(double a) const
{
        //      a_point horizontal(0,1,0);
        a_point y;
        a_point z = this->dz_axis(a);
        if ((z-horiz_).norm()<1e-5)//the plane is not vertical // xz
                y = a_point(0,0,1);
        else if ((z+horiz_).norm()<1e-5)//the plane is not vertical // xz
                y = a_point(0,0,-1);
        else
                y = (cross(z,horiz_)).normalise();
        return y;
}
//---------------------------------------------------------------------------
a_point a_curve_lin::dz_axis(double a) const
{
        if (a>1)
                return this->dz_axis(a-1);
        else if (a==1)//at the extremity, take the direction of the last segment
        {
                int n = x_.size();
                a_point x1 = *x_[n-2];
                a_point x2 = *x_[n-1];
                return (x2-x1).normalise();
        }
        int ref = this->ref(a);
        //interpolate the tangent
        double a1 = a_[ref];
        double a2 = a_[ref+1];
        double f = (a-a1)/(a2-a1);
        a_point x1 = a_curve::dz_axis(a1);
        a_point x2 = a_curve::dz_axis(a2);
        return (x1 + (x2-x1)*f).normalise();
}
//---------------------------------------------------------------------------
a_point a_curve_lin::tangent(double a) const
{
        if (a>=1)
                return a_curve::dz_axis(a);
        int ref = this->ref(a);
        //interpolate the tangent
        double a1 = a_[ref];
        double a2 = a_[ref+1];
        double f = (a-a1)/(a2-a1);
        a_point x1 = a_curve::dz_axis(a1);
        a_point x2 = a_curve::dz_axis(a2);
        return x1 + (x2-x1)*f;
}
//---------------------------------------------------------------------------
a_point a_curve_lin::curvature_centre() const
{
        int n = x_.size();
        a_point p1,p2,p3;
        if (this->is_closed())
        {
                p1 = *x_[0];
                p2 = *x_[int((n-1)/3.+.5)];
                p3 = *x_[int((n-1)/1.5+.5)];
        }
        else
        {
                p1 = *x_[0];
                p2 = *x_[int((n-1)/2.+.5)];
                p3 = *x_[n-1];
        }
        return circle_centre(p1,p2,p3);
}
//---------------------------------------------------------------------------
a_point a_curve_lin::curvature_centre(double a) const
{
        a_point p1,p2,p3;
        if (a>0.99)
        {
                p1 = (*this)(0.98);
                p2 = (*this)(0.99);
                p3 = (*this)(1.);
        }
        else if (a>0.01)
        {
                p1 = (*this)(a-0.01);
                p2 = (*this)(a);
                p3 = (*this)(a+0.01);
        }
        else
        {
                p1 = (*this)(0.);
                p2 = (*this)(0.01);
                p3 = (*this)(0.02);
        }
        return circle_centre(p1,p2,p3);
}
//---------------------------------------------------------------------------
double a_curve_lin::curvature(double a) const
{
        a_point xm = this->tangent(a-da_);
        a_point xM = this->tangent(a+da_);
        a_point dx = (xM-xm)/da_;
        //crash if 0
        return 1/dx.norm();
}
//---------------------------------------------------------------------------
a_curve_lin&  a_curve_lin::extrapolate(double f1, double f2)
{
        if (f1>0) return *this;
        if (f2<1) return *this;
        a_point centre = this->curvature_centre();
        a_point p1 = *x_[0];
        a_point p2 = *x_[n_pts_-1];
        a_point r1 = p1-centre;
        a_point r2 = p2-centre;
        double rn1 = r1.norm();
        double rn2 = r2.norm();
        r1.normalise();
        r2.normalise();
        a_point d = p1-p2;
        a_point normal = cross(r1,d).normalise();
        double da = acos(r1*r2)/(2.*n_pts_);
        if (f1<0)
        {
                int n = (int)fabs(f1)*n_pts_;
                double ds = 2*rn1*sin(da);
                a_point x = p1;
                for (int i=0; i < n; i++)
                {
                        a_point uy = (centre-x).normalise();
                        a_point ux = cross(normal,uy).normalise();
                        x -= ds*(ux*cos(da)+uy*sin(da));
                        auto xp = new a_point(x);
                        x_.push_front(xp);
                }
        }
        if (f2>1)
        {
                int n = (int)(f2-1)*n_pts_;
                double ds = 2*rn2*sin(da);
                a_point x = p2;
                for (int i=0; i < n; i++)
                {
                        a_point uy = (centre-x).normalise();
                        a_point ux = cross(normal,uy).normalise();
                        x += ds*(ux*cos(da)+uy*sin(da));
                        auto xp = new a_point(x);
                        x_.push_back(xp);
                }
        }
        a_.clear();
        this->init();
        return *this;
}
//***************************************************************************
//---------------------------------------------------------------------------
std::istream& operator>> (std::istream& in, a_curve_lin & f)
{
        in >> *((a_curve *)(&f));
        f.init();
        return in;
}
//***************************************************************************
//---------------------------------------------------------------------------
a_curve_lin interpolate(const a_curve_lin& pat1,const a_curve_lin& pat2, double a)
{
        std::vector<double> merge;
        int n1 = pat1.size();
        int n2 = pat2.size();
        if (n1!=n2)
                throw a_curve_lin::pattern_incompatible_error();
        a_curve_lin patn;
        for (int i = 0; i < n1; i++)
        {
                a_point p1 = *pat1[i];
                a_point p2 = *pat2[i];
                auto p = new a_point(p1+a*(p2-p1));
                patn.x_.push_back(p);
        }
        patn.init();
        return patn;
}