MC_triangle.cpp

Go to the documentation of this file.

#include <MC_segment.hpp>
#include <MC_double_vector.hpp>


#include <MC_triangle.hpp>


namespace mesh_conv
{

    MC_v3d MC_triangle::normal() const
    {

        MC_v3d v0=v[1]-v[0];
        MC_v3d v1=v[2]-v[0];

        MC_v3d n=v0.cross(v1);
        double nn=n.norm();
        double epsilon=0.00000001;
        if(nn<epsilon)
        {std::cout<<"Warning in MC_triangle::normal(), from triangle "<<*this<<", normal is "<<n<<" of norm()="<<nn<<std::endl;return MC_v3d(0,0,1);}
        n/=nn;
        return n;
    }

    MC_v3d MC_triangle::barycenter() const
    {return 1.0/3.0*(v[0]+v[1]+v[2]);}

    MC_triangle::MC_triangle(){resize(3);}
    MC_triangle::MC_triangle(const MC_v3d& v0,const MC_v3d& v1,const MC_v3d& v2)
            :MC_v3d_vector(v0,v1,v2)
    {}
    MC_triangle::MC_triangle(const MC_v3d_vector& vec)
            :MC_v3d_vector()
    {
        if(vec.size()!=3)
        {std::cout<<"Error in MC_triangle::MC_triangle("<<vec<<"), size should be 3"<<std::endl;exit(-1);}
        resize(3);
        v[0]=vec[0];v[1]=vec[1];v[2]=vec[2];
    }
    MC_triangle::MC_triangle(const std::vector <MC_v3d>& vec)
            :MC_v3d_vector()
    {
        if(vec.size()!=3)
        {std::cout<<"Error in MC_triangle::MC_triangle(std::vector <MC_v3d>), size should be 3 and not "<<vec.size()<<std::endl;exit(-1);}
        resize(3);
        v[0]=vec[0];v[1]=vec[1];v[2]=vec[2];
    }
    MC_triangle::MC_triangle(const MC_triangle& p)
            :MC_v3d_vector(p)
    {v=p.v;}

    MC_v3d MC_triangle::closest_point(const MC_v3d& x,int *type) const
    {

        double epsilon=0.00001;

        // First project and look if the projection is inside the triangle (barycentric coordinates)
        MC_v3d n=normal();
        MC_v3d proj = x-((x-v[0]).dot(n))*n;

        // if the proj is inside the triangle, this is the closest point
        if(is_inside(proj)==1)
        {
            if(type!=0)
                *type=0;
            return proj;
        }


        // if not, take the closest line
        if(type!=0)
            *type=1;
        MC_segment s0=segment(0),s1=segment(1),s2=segment(2);
        MC_v3d y0=s0.closest_point(x),y1=s1.closest_point(x),y2=s2.closest_point(x);
        double dist_1=(y0-x).norm(),dist_2=(y1-x).norm(),dist_3=(y2-x).norm();
        if(dist_1<=dist_2 && dist_1<=dist_3)
        {
            //check if its a vertex point
            for(int k=0;k<3;k++){if((y0-v[k]).norm()<epsilon) *type=2;}
            return y0;
        }
        else if(dist_2<=dist_1 && dist_2<=dist_3)
        {
            for(int k=0;k<3;k++){if((y1-v[k]).norm()<epsilon) *type=2;}
            return y1;
        }
        else
        {
            for(int k=0;k<3;k++){if((y2-v[k]).norm()<epsilon) *type=2;}
            return y2;
        }

    }

    bool MC_triangle::is_inside(const MC_v3d& _x) const
    {
        std::pair <MC_double_vector,bool> bar=barycentric_coordinates(_x);
        if(bar.second==false) //not in plane
            return false;
        double epsilon=0.0001;
        if(bar.first[0]>1.0 || bar.first[1]>1.0 || bar.first[2]>1.0 || (bar.first[0]+bar.first[1]+bar.first[2]-1.0)>epsilon) // outside the polygon
            return false;

        // in the polygon
        return true;
    }

    std::pair<MC_double_vector,bool> MC_triangle::barycentric_coordinates(const MC_v3d& x) const
    {
        // check if _x is in the plane
        double epsilon=0.0001;
        double L = (x-v[0]).dot(normal());
        if(L>epsilon)//not in the plane
            return std::pair <MC_double_vector,bool> (MC_double_vector(-1,-1,-1),false);

        // calculate the barycentric coordinates
        double area_1=MC_v3d::area(v[1]-v[0],x-v[0]);
        double area_2=MC_v3d::area(v[2]-v[1],x-v[1]);
        double area_3=MC_v3d::area(v[0]-v[2],x-v[2]);
        double area_tot = area();

        MC_double_vector b(area_1/area_tot,area_2/area_tot,area_3/area_tot);
        return std::pair <MC_double_vector,bool> (b,true);
    }

    MC_segment MC_triangle::segment(const int& edge_number) const
    {
        if(edge_number<0 || edge_number>2)
        {std::cout<<"MC_triangle::segment("<<edge_number<<"),_k_edge has to be between 0 and 2"<<std::endl;exit(-1);}

        return MC_segment(v[edge_number],v[(edge_number+1)%3]);
    }

    double MC_triangle::area() const
    {
        return MC_v3d::area(v[1]-v[0],v[2]-v[0]);
    }

    MC_triangle::MC_triangle(const MC_polygon& p)
    {
        if(p.size()!=3)
        {std::cout<<"Error in MC_triangle::MC_triangle("<<p<<"), size must be 3 and not "<<p.size()<<std::endl;exit(-1);}
        resize(3);
        v[0]=p[0];v[1]=p[1];v[2]=p[2];
    }
    std::vector<MC_triangle> MC_triangle::triangulate(const MC_polygon& p)
    {
        int N=p.size();
        std::vector <MC_triangle> vtri(N-2);
        for(int k_triangle=0;k_triangle<N-2;k_triangle++)
            vtri[k_triangle]=MC_triangle(p[0],p[k_triangle+1],p[k_triangle+2]);
        return vtri;
    }

    MC_v3d MC_triangle::segment_intersection(const MC_segment& s,int *type) const
    {
        int temp_type=0;
        MC_v3d i=s.plane_intersection(normal(),v[0],&temp_type);

        if(type!=0)
            *type=temp_type;

        if(temp_type==0 || temp_type==3)
            return i;

        if(is_inside(i)==false)
            if(type!=0)
                *type=0;

        return i;
    }

    MC_v3d MC_triangle::pn_triangle(const MC_v3d_vector& X,const MC_v3d_vector& N,const MC_double_vector& x)
    {
        MC_double_vector w=MC_double_vector::zeros(9);
        for(int i=0;i<3;++i)
            for(int j=0;j<3;++j)
                w[3*i+j]=(X[j]-X[i]).dot(N[i]);

        MC_v3d b300=X[0];
        MC_v3d b030=X[1];
        MC_v3d b003=X[2];

        MC_v3d b210=(2*X[0]+X[1]-w[0*3+1]*N[0])/3.0;
        MC_v3d b120=(2*X[1]+X[0]-w[1*3+0]*N[1])/3.0;
        MC_v3d b021=(2*X[1]+X[2]-w[1*3+2]*N[1])/3.0;
        MC_v3d b012=(2*X[2]+X[1]-w[2*3+1]*N[2])/3.0;
        MC_v3d b102=(2*X[2]+X[0]-w[2*3+0]*N[2])/3.0;
        MC_v3d b201=(2*X[0]+X[2]-w[0*3+2]*N[0])/3.0;
        MC_v3d E=(b210+b120+b021+b012+b102+b201)/6.0;
        MC_v3d V=(X[0]+X[1]+X[2])/3.0;
        MC_v3d b111=E+(E-V)/2.0;

        MC_v3d y=b300*x[0]*x[0]*x[0]+b030*x[1]*x[1]*x[1]+b003*x[2]*x[2]*x[2]+
                 b210*3.0*x[0]*x[0]*x[1]+b120*3.0*x[0]*x[1]*x[1]+b201*3.0*x[0]*x[0]*x[2]+
                 b021*3.0*x[1]*x[1]*x[2]+b102*3.0*x[0]*x[2]*x[2]+b012*3.0*x[1]*x[2]*x[2]+
                 b111*6.0*x[0]*x[1]*x[2];

        return y;
    }

    MC_matrix MC_triangle::inertia() const
    {
        MC_v3d x0=(*this)[0];
        MC_v3d x1=(*this)[1];
        MC_v3d x2=(*this)[2];

        double a=((x1-x0).cross(x2-x0)).norm();

        MC_matrix S(3,3);
        S(0,0)=2;S(0,1)=1;S(0,2)=1;
        S(1,0)=1;S(1,1)=2;S(1,2)=1;
        S(2,0)=1;S(2,1)=1;S(2,2)=2;
        S*=1.0/24.0;

        MC_matrix V(x0,x1,x2);

        MC_matrix C=a*V*S*V.transposed();
        MC_matrix J=C.trace()*MC_matrix(3)-C;

        return J;
    }

}