MC_segment.cpp

Go to the documentation of this file.

#include <MC_segment.hpp>

#include <MC_v3d_vector.hpp>
#include <MC_double_vector.hpp>
#include <MC_int_pair.hpp>

namespace mesh_conv
{

    MC_segment::MC_segment(){}
    MC_segment::MC_segment(const MC_v3d& start_point,const MC_v3d& end_point)
    {x0=start_point;x1=end_point;}
    MC_segment::MC_segment(const MC_segment& s){x0=s.x0;x1=s.x1;}

    MC_segment::~MC_segment(){}


    MC_v3d MC_segment::unit_vector() const {return vector().normalized();}
    MC_v3d MC_segment::vector() const {return x1-x0;}
    MC_segment MC_segment::fliped() const {return MC_segment(x1,x0);}


    double MC_segment::length() const {return vector().norm();}
    MC_v3d MC_segment::closest_point(const MC_v3d& x) const
    {
        MC_v3d u=unit_vector();
        double proj=0.0;
        proj = (x-x0).dot(u);
        if(proj<0)
            return x0;
        else if(proj>length())
            return x1;
        else
            return x0+proj*u;
    }

    double MC_segment::distance_to_point(const MC_v3d& x) const {return (x-closest_point(x)).norm();}
    MC_segment MC_segment::closest_to_segment(const MC_segment& s) const
    {
        //****************************************************//
        // 1. Search the closest position between the two lines
        //    = line perpendicular to both segments
        //****************************************************//


        // useful vectors
        MC_v3d u1,u2;
        MC_v3d A12;

        u1 = vector();
        u2 = s.vector();
        A12 = s[0]-x0;

        double det=0.0;
        det = powf(float(u1.dot(u2)),2.0) - u1.dot(u1)*u2.dot(u2);

        // check if the lines are parallel
        double epsilon=0.00001;
        if(fabs(det)<epsilon)
        {std::cout<<"Error in MC_segment::closest_to_segment(), lines are almost parallel"<<std::endl;exit(-1);}

        // distance of the intersection
        double s1=0.0,s2=0.0;
        s1 = 1/det*(-(A12.dot(u1))*(u2.dot(u2)) + (A12.dot(u2))*(u1.dot(u2)) );
        s2 = 1/det*(-(A12.dot(u1))*(u2.dot(u1)) + (A12.dot(u2))*(u1.dot(u1)) );

        MC_v3d A1=x0,A2=s[0];
        MC_v3d P1;P1 = A1 + s1*u1;
        MC_v3d P2;P2 = A2 + s2*u2;

        // check if the points (P1,P2) are inside the two segments
        double dot_p1=0.0,dot_p2=0.0;
        dot_p1 = (P1-A1).dot(u1);
        dot_p2 = (P2-A2).dot(u2);

        MC_segment seg;
        int is_valid_p1=0,is_valid_p2=0;
        if(dot_p1>=0 && dot_p1<=u1.dot(u1))
            is_valid_p1=1;
        if(dot_p2>=0 && dot_p2<=u2.dot(u2))
            is_valid_p2=1;

        if(is_valid_p1==1 && is_valid_p2==1)
            seg=MC_segment(P1,P2);
        else
        {
            // else the closest points are at the extremities
            // 9 possibilities in this cases: check all of these
            MC_v3d y0[3]={A1,A1+u1,P1};
            MC_v3d y1[3]={A2,A2+u2,P2};
            double min_dist=99999.99;
            double current_dist=0.0;
            int min_k1=-1,min_k2=-1;
            for(int k1=0;k1<3;k1++)
                for(int k2=0;k2<3;k2++)
                {
                current_dist=(y0[k1]-y1[k2]).norm();
                if(current_dist<min_dist)
                    if(k1!=2 || (k1==2 && is_valid_p1==1))
                        if(k2!=2 || (k2==2 && is_valid_p2==1))
                        {min_dist=current_dist;min_k1=k1;min_k2=k2;}
            }

            // 9 cases
            if(min_k1==0 && min_k2==0) {seg=MC_segment(A1,A2);}
            else if(min_k1==0 && min_k2==1) {seg=MC_segment(A1,A2+u2);}
            else if(min_k1==0 && min_k2==2) {seg=MC_segment(A1,P2);}

            else if(min_k1==1 && min_k2==0) {seg=MC_segment(A1+u1,A2);}
            else if(min_k1==1 && min_k2==1) {seg=MC_segment(A1+u1,A2+u2);}
            else if(min_k1==1 && min_k2==2) {seg=MC_segment(A1+u1,P2);}

            else if(min_k1==2 && min_k2==0) {seg=MC_segment(P1,A2);}
            else if(min_k1==2 && min_k2==1) {seg=MC_segment(P1,A2+u2);}
            else if(min_k1==2 && min_k2==2) {seg=MC_segment(P1,P2);}

        }

        return seg;

    }


    double MC_segment::distance_to_segment(const MC_segment& s) const
    {return closest_to_segment(s).length();}

    MC_v3d MC_segment::plane_intersection(const MC_v3d& n,const MC_v3d& _x0,int *type) const
    {
        double epsilon=0.00001;
        //first check if the segment is parallel to the plane
        MC_v3d AB=x1-x0;
        if(fabs(AB.dot(n))<epsilon)
        {
            //plane is \\ to the segment
            //now test if the plane pass through the line
            if(fabs(n.dot(_x0-x0))<epsilon)//pass through
            {
                if(type!=0)
                    *type=3;
                return MC_v3d(-1,-1,-1);
            }
            else//no intersection
            {
                if(type!=0)
                    *type=0;
                return MC_v3d(-1,-1,-1);
            }
        }

        //We are sure that the plane is not \\ to the segment
        double t=n.dot(_x0-x0)/AB.dot(n);

        MC_v3d intersection=x0+t*AB;



        //check if the intersection is inside the segment
        if(t>=0 && t<=1)
        {
            //check if its close to a vertex or not
            if( ((intersection-x0).norm()<epsilon) || ((intersection-x1).norm()<epsilon) )
            {if(type!=0) *type=2;}
            else
            {if(type!=0) *type=1;}
            return intersection;
        }

        //else outside the segment => no intersection
        if(type!=0) *type=0;
        return intersection;
    }

    bool MC_segment::is_belonging(const MC_v3d& x,int *type) const
    {
        bool belongs=0;
        double t = relative_position(x,&belongs);
        if(belongs!=1)
        {
            if(type!=0) *type=-2;
            return false;
        }
        else if(t<0 || t>1)
        {
            if(type!=0) *type=-1;
            return false;
        }
        if(type!=0) *type=1;
        return true;
    }
    bool MC_segment::is_aligned(const MC_v3d& x) const
    {
        double epsilon=0.0001;
        return ((x-x0)-(x-x0).dot(unit_vector())*unit_vector()).norm()<epsilon;
    }
    double MC_segment::relative_position(const MC_v3d& x,bool *_is_aligned) const
    {
        // check if it is aligned
        if(is_aligned(x)==false)
        {
            if(_is_aligned!=0)
                *_is_aligned=false;
            return 0.0;}

        if(_is_aligned!=0)
            *_is_aligned=true;
        return (x-x0).dot(unit_vector())/length();
    }
    void MC_segment::assert_bounds(const int& u) const
    {
        if(u!=0 && u!=1)
        {std::cout<<"Error in MC_segment::assert_bounds("<<u<<"), value should be 0 or 1"<<std::endl;exit(-1);}
    }
    MC_v3d  MC_segment::operator()(const int& index) const
    {assert_bounds(index); if(index==0)return x0;if(index==1) return x1;exit(-1);}
    MC_v3d& MC_segment::operator()(const int& index)
    {assert_bounds(index); if(index==0)return x0;if(index==1) return x1;exit(-1);}
    MC_v3d  MC_segment::operator[](const int& index) const
    {assert_bounds(index); if(index==0)return x0;if(index==1) return x1;exit(-1);}
    MC_v3d& MC_segment::operator[](const int& index)
    {assert_bounds(index); if(index==0)return x0;if(index==1) return x1;exit(-1);}

    MC_v3d MC_segment::value(const double& t) const
    {return (1-t)*x0+t*x1;}
    MC_v3d_vector MC_segment::value(const MC_double_vector& t) const
    {return (1-t)*x0+t*x1;}



    std::pair<MC_v3d_vector,MC_double_vector> MC_segment::linear_sampling_intervals(const double& d_L) const
    {
        double epsilon=0.00001;
        if(d_L<=epsilon)
        {std::cout<<"Error in MC_segment::linear_sampling_intervals("<<d_L<<"), interval is too low"<<std::endl;exit(-1);}

        double total_length=length();
        int N_subdiv=int(total_length/d_L+1);

        return linear_sampling(N_subdiv);
    }
    std::pair<MC_v3d_vector,MC_double_vector> MC_segment::linear_sampling(const int& N) const
    {
        MC_double_vector t_value=MC_double_vector::zeros(N);
        MC_v3d_vector sampling=MC_v3d_vector::zeros(N);
        for(int k=0;k<N;k++)
        {
            double u=double(k)/double(N-1);
            sampling(k) = (1-u)*x0 + u*x1;
            t_value(k)=u;
        }
        return std::pair<MC_v3d_vector,MC_double_vector> (sampling,t_value);
    }
    MC_segment operator+(const MC_segment& s,const MC_v3d& x)
    {return MC_segment(s[0]+x,s[1]+x);}
    MC_segment operator+(const MC_v3d& x,const MC_segment& s)
    {return s+x;}
    MC_segment operator*(const MC_segment& s,const double& a)
    {return MC_segment(s[0]*a,s[1]*a);}
    MC_segment operator*(const double& a,const MC_segment& s)
    {return s*a;}

    std::ostream& operator << (std::ostream& stream, const MC_segment& s)
    {stream<<"["<<s[0]<<" ; "<<s[1]<<"]";return stream;}

    std::vector<MC_segment> MC_segment::build_segment_vector(const MC_v3d_vector& vertices,const std::vector<MC_int_pair>& constraints)
    {
        unsigned int N_vertices=vertices.size();
        std::vector<MC_segment> v_seg;
        for(int k=0,N=constraints.size();k<N;++k)
        {
            const MC_int_pair& index=constraints[k];
            if(index[0]>=N_vertices || index[1]>=N_vertices)
            {std::cout<<"Error in MC_segment::build_segment_vector(), index "<<index<<" is too large compared to N_vertex="<<N_vertices<<" at k="<<k<<std::endl;exit(-1);}
            v_seg.push_back(MC_segment(vertices[index[0]],vertices[index[1]]));
        }
        return v_seg;
    }


}