a_curve.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.h"
#include "a_segment.h"
#include <iomanip>
#include <stdlib.h>
 
// .PORTABILITY : ansi C++
 
//***************************************************************************
const a_point horiz0(0,1,0);
//---------------------------------------------------------------------------
a_curve::a_curve(double da) : da_(da)
{
        initialised_ = false;
        orient_ = 1;
        n_pts_ = 0;
        horiz_ = horiz0;
}
//---------------------------------------------------------------------------
a_curve::a_curve(const a_curve& ifn)
{
        initialised_ = false;
        da_ = ifn.da();
        orient_ = ifn.orient_;
        horiz_ = ifn.horiz_;
        for (int i = 0; i < ifn.size(); i++)
        {
                auto p = new a_point;
                *p = *ifn[i];
                x_.push_back(p);
        }
        this->init();
}
//---------------------------------------------------------------------------
a_curve::~a_curve()
{
        for (auto c:x_)
                delete c;
        x_.clear();
}
//---------------------------------------------------------------------------
a_curve& a_curve::init()
{
        if (!initialised_)
        {
                initialised_ = true;
                n_pts_ = x_.size();
                this->construct_axis();
        }
        return *this;
}
//---------------------------------------------------------------------------
a_point a_curve::dz_axis(double a) const
{
        //tangent along the path
        a_point xi,xs;
        if (a-da_<0)
        {
                if (this->is_closed())
                {
                        xi = (*this)(a-da_+1.);
                        xs = (*this)(a+da_);
                        return (.5*(xs-xi)/da_);
                }
                else
                {
                        xi = (*this)(a);
                        xs = (*this)(a+da_);
                        return ((xs-xi)/da_);
                }
        }
        else if (a+da_>1)
        {
                if (this->is_closed())
                {
                        xi = (*this)(a-da_);
                        xs = (*this)(a+da_-1.);
                        return (.5*(xs-xi)/da_);
                }
                else
                {
                        xi = (*this)(a-da_);
                        xs = (*this)(a);
                        return ((xs-xi)/da_);
                }
        }
        else
        {
                xi = (*this)(a-da_);
                xs = (*this)(a+da_);
                return (.5*(xs-xi)/da_);
        }
}
//---------------------------------------------------------------------------
bool a_curve::is_closed() const
{
        a_point * p0 = (*this)[0];
        a_point * pn = (*this)[this->size()-1];
        return ((*p0)==(*pn));
}
//---------------------------------------------------------------------------
a_point a_curve::dy_axis(double a) const
{
        a_point xi,xs;
        if (a-da_<0)
        {
                xi = this->dz_axis(a);
                xs = this->dz_axis(a+da_);
                return ((xs-xi)/da_);
        }
        else if (a-da_>1)
        {
                xi = this->dz_axis(a-da_);
                xs = this->dz_axis(a);
                return ((xs-xi)/da_);
        }
        else
        {
                xi = this->dz_axis(a-da_);
                xs = this->dz_axis(a+da_);
                return (.5*(xs-xi)/da_);
        }
}
//---------------------------------------------------------------------------
a_point a_curve::dx_axis(double a) const
{
        a_point z = this->dz_axis(a);
        a_point y = this->dy_axis(a);
        return cross(y,z);
}
//---------------------------------------------------------------------------
a_curve& a_curve::origin(a_point & p)
{
        a_point pd = p-*x_[0];
        return this->translate(pd);
}
//---------------------------------------------------------------------------
a_curve& a_curve::origin(double x, double y, double z)
{
        a_point p(x,y,z);
        return this->origin(p);
}
//---------------------------------------------------------------------------
a_point a_curve::centre() const
{
        a_point sum;
        int n_pts = x_.size();
        for (int i = 0; i < n_pts; i++)
                sum += *x_[i];
        return sum/n_pts;
}
//---------------------------------------------------------------------------
/*a_curve& a_curve::invert()
        {
        int n_pts = x_.size();
        std::vector<a_point *>  x(n_pts);
        for (int i = 0; i < n_pts; i++)
        x[i] = x_[n_pts-i-1];
        for (int i = 0; i < n_pts; i++)
        x_[i] = x[i];
        x.clear();
        return *this;
        }*/
//---------------------------------------------------------------------------
a_curve& a_curve::translate(const a_point & p)
{
        int n_pts = x_.size();
        for (int i = 0; i < n_pts; i++)
                *x_[i] += p;
        return *this;
}
//---------------------------------------------------------------------------
a_curve& a_curve::translate(double x, double y, double z)
{
        a_point p(x,y,z);
        return this->translate(p);
}
//---------------------------------------------------------------------------
a_curve& a_curve::flatten()
{
        //project the set of point on the plan passing by the centre of the set and having 
        double f0 = (*x_[0])*z_axis_;
        for (int i = 0; i < n_pts_; i++)
        {
                double f = (*x_[i])*z_axis_;
                a_point p = *x_[i] - z_axis_*(f-f0);
                *x_[i] = p;
        }
        initialised_ = false;
        return this->init();
}
//---------------------------------------------------------------------------
a_curve& a_curve::flatten2D()
{
        a_point x0 =  *x_[0];
        for (int i = 0; i < n_pts_; i++)
        {
                a_point x1 = *x_[i]-x0;
                x_[i]->set(x1*x_axis_,x1*y_axis_,0);
        }
        initialised_ = false;
        return this->init();
}
//---------------------------------------------------------------------------
a_curve& a_curve::place3D(a_point& origin, a_point& x_axis, a_point& y_axis)
{
        for (int i = 0; i < n_pts_; i++)
        {
                a_point po = *x_[i];
                *x_[i] = origin+po.x()*x_axis+po.y()*y_axis;
        }
        initialised_ = false;
        return this->init();
}
//---------------------------------------------------------------------------
a_curve& a_curve::construct_x_axis()
{
        x_axis_ = cross(y_axis_,z_axis_).normalise();
        return *this;
}
//---------------------------------------------------------------------------
a_curve& a_curve::construct_y_axis()
{
        //      a_point horizontal(1,0,0);
        if ((z_axis_-horiz_).norm()<1e-5)//the plane is vertical // xz
                y_axis_ = a_point(0,0,1);
        else if ((z_axis_+horiz_).norm()<1e-5)//the plane is vertical // xz
                y_axis_ = a_point(0,0,-1);
        else
                y_axis_ = (cross(z_axis_,horiz_)).normalise();
        return *this;
}
//---------------------------------------------------------------------------
//find the average perpendicular to the set of points
a_curve& a_curve::construct_z_axis()
{
        a_point centre = this->centre();
        a_point perp;
        for (int i = 0; i < n_pts_-1; i++)
        {
                a_point pt = *x_[i+1]-*x_[i];
                a_point pr = *x_[i]-centre;
                perp += cross(pr,pt);
        }
        a_point pt = *x_[0]-**(x_.end()-1);
        a_point pr = **(x_.end()-1)-centre;
        perp += cross(pr,pt);
        perp.normalise();
        z_axis_ = orient_*perp;
        return *this;
}
//---------------------------------------------------------------------------
a_curve& a_curve::construct_axis()
{
        a_point centre = this->centre();
        a_point perp;
        for (int i = 0; i < n_pts_-1; i++)
        {
                a_point pt = *x_[i+1]-*x_[i];
                a_point pr = *x_[i]-centre;
                perp += cross(pr,pt);
        }
        a_point pt = *x_[0]-**(x_.end()-1);
        a_point pr = **(x_.end()-1)-centre;
        perp += cross(pr,pt);
        if (perp.norm()<1e-8) //only a line
        {
                x_axis_ = (*x_[1]-*x_[0]).normalise();
                a_point p(0,0,1);
                if (cross(p,x_axis_).norm()<1e-8)
                {
                        std::cerr << "this case is not yet supported (see code a_curve)" << std::endl;
                        exit(-1);
                }
                p += x_axis_;
                perp = cross(p,x_axis_).normalise();
                y_axis_ = orient_*perp;
                z_axis_ = cross(x_axis_,y_axis_).normalise();
                return *this;
        }
        else
                return this->construct_z_axis().construct_y_axis().construct_x_axis();
}
//---------------------------------------------------------------------------
double a_curve::dist(const a_point & p, int & ref) const
{
        double d0 = 1e30;
        int i0;
        int n_pts = x_.size();
        for (int i = 0; i < n_pts; i++)
        {
                double d = ((*x_[i])-p).norm();
                if (d<d0)
                {
                        d0 = d;
                        ref = i;
                }
        }
        return d0;
}
//---------------------------------------------------------------------------
double a_curve::dist2(const a_point & p, int & ref, double & m) const
{
        double d0 = 1e30;
        a_point a0;
        int i0;
        int n_pts = x_.size();
        for (int i = 0; i < n_pts-1; i++)
        {
                a_segment s(*x_[i],*x_[i+1]);
                a_point a = s.closest(p);
                double d = (a-p).norm();
                if (d<d0)
                {
                        d0 = d;
                        a0 = a;
                        ref = i;
                }
        }
        m = (a0-*x_[ref]).norm()/(*x_[ref+1]-*x_[ref]).norm();
        return d0;
}
//---------------------------------------------------------------------------
double a_curve::dist_last(const a_point & p, int & ref) const
{
        double d0 = 1e30;
        int i0;
        int n_pts = x_.size();
        for (int i = 0; i < n_pts; i++)
        {
                double d = ((*x_[i])-p).norm();
                if (d<=d0)
                {
                        d0 = d;
                        ref = i;
                }
        }
        return d0;
}
//---------------------------------------------------------------------------
double a_curve::dist(const a_curve & fn, int & ref1, int & ref2) const
{
        double d0 = 1e30;
        int i0;
        int n_pts = fn.size();
        for (int i = 0; i < n_pts; i++)
        {
                int ref1t;
                double d = this->dist(*fn[i],ref1t);
                if (d<d0)
                {
                        d0 = d;
                        ref1 = ref1t;
                        ref2 = i;
                }
        }
        return d0;
}
//---------------------------------------------------------------------------
double a_curve::dist_last(const a_curve & fn, int & ref1, int & ref2) const
{
        double d0 = 1e30;
        int i0;
        int n_pts = fn.size();
        for (int i = 0; i < n_pts; i++)
        {
                int ref1t;
                double d = this->dist_last(*fn[i],ref1t);
                if (d<=d0)
                {
                        d0 = d;
                        ref1 = ref1t;
                        ref2 = i;
                }
        }
        return d0;
}
//---------------------------------------------------------------------------
a_segment a_curve::shortest(const a_curve & fn, int & ref1, int & ref2, double& m1, double& m2) const
{
        a_segment sshort;
        const double eps = 1e-10;
        double d0 = 1e30;
        double mm1,mm2;
        int n_pt2 = fn.size();
        for (int i2 = 0; i2 < n_pt2-1; i2++)
        {
                a_segment seg2(*fn[i2],*fn[i2+1]);
                if (seg2.length()<eps)
                {
                        std::cerr << "error: |seg2(" << i2+1 << ")| ~= 0" << std::endl;
                        exit(-1);
                }
                int n_pt1 = x_.size();
                for (int i1 = 0; i1 < n_pt1-1; i1++)
                {
                        a_segment seg1(*x_[i1],*x_[i1+1]);
                        if (seg1.length()<eps)
                        {
                                std::cerr << "error: |seg1(" << i1+1 << ")| ~= 0" << std::endl;
                                exit(-1);
                        }
                        a_segment inter = seg1.shortest(seg2,mm1,mm2);
                        double d = inter.length();
                        if (d<d0)
                        {
                                sshort = inter;
                                d0 = d;
                                ref1 = i1;
                                ref2 = i2;
                                m1 = mm1;
                                m2 = mm2;
                                if (mm2<0)
                                {
                                        std::cerr << "-----------------" << std::endl;
                                        std::cerr << seg1 << std::endl;
                                        std::cerr << seg2 << std::endl;
                                        exit(-1);
                                }
                        }
                }
        }
        return sshort;
}
//---------------------------------------------------------------------------
double a_curve::dist2(const a_curve & fn, int & ref1, int & ref2, double& m1, double& m2) const
{
        const double eps = 1e-10;
        double d0 = 1e30;
        double mm1,mm2;
        int n_pt2 = fn.size();
        for (int i2 = 0; i2 < n_pt2-1; i2++)
        {
                a_segment seg2(*fn[i2],*fn[i2+1]);
                if (seg2.length()<eps)
                {
                        std::cerr << "error: |seg2(" << i2+1 << ")| ~= 0" << std::endl;
                        exit(-1);
                }
                int n_pt1 = x_.size();
                for (int i1 = 0; i1 < n_pt1-1; i1++)
                {
                        a_segment seg1(*x_[i1],*x_[i1+1]);
                        if (seg1.length()<eps)
                        {
                                std::cerr << "error: |seg1(" << i1+1 << ")| ~= 0" << std::endl;
                                exit(-1);
                        }
                        a_segment inter = seg1.shortest(seg2,mm1,mm2);
                        double d = inter.length();
                        if (d<d0)
                        {
                                d0 = d;
                                ref1 = i1;
                                ref2 = i2;
                                m1 = mm1;
                                m2 = mm2;
                        }
                }
        }
        return d0;
}
//***************************************************************************
//---------------------------------------------------------------------------
std::istream& operator>> (std::istream& in, a_curve & f)
{
        int n_pts;
        in >> n_pts;
        for (int i = 0; i < n_pts; i++)
        {
                double x,y,z;
                in >> x >> y >> z;
                auto p = new a_point(x,y,z);
                f.x_.push_back(p);
        }
        return in;
}
//---------------------------------------------------------------------------
std::ostream& operator<< (std::ostream& o, const a_curve & f)
{
        int n_pts = f.x_.size();
        o << n_pts << std::endl;
        for (int i = 0; i < n_pts; i++)
                o << f.x_[i] << std::endl;
        return o;
}