MC_polygon.cpp

Go to the documentation of this file.

#include <MC_triangle.hpp>
#include <MC_segment.hpp>
#include <MC_int_vector.hpp>
#include <MC_matrix.hpp>
#include <MC_int_pair.hpp>
#include <MC_mesh_index_vector.hpp>

#include <MC_polygon.hpp>


namespace mesh_conv
{
    bool MC_polygon::is_planar() const
    {

        int N=size();
        if(N<3)
        {std::cout<<"Error in MC_polygon::is_planar(), size is "<<N<<"<3"<<std::endl;exit(-1);}
        MC_v3d n=normal();

        double tol=0.0001;
        for(int k=1;k<N;k++)
            if( fabs((v[k]-v[0]).dot(n))>tol )
                return false;
        return true;
    }

    MC_v3d MC_polygon::normal() const
    {
        int N=size();
        if(N<3)
        {std::cout<<"Error in MC_polygon::normal(), size is "<<N<<"<3"<<std::endl;exit(-1);}

        MC_v3d n;
        for(int k=0;k<N-2;k++)
        {
            MC_v3d v0,v1;
            v0=(v[k+1]-v[k]).normalized();
            v1=(v[k+2]-v[k]).normalized();
            n+=v0.cross(v1);
        }

        double nn=n.norm();

        double epsilon=0.000000001;
        if(nn<epsilon)
        {
            //std::cout<<"Warning in MC_polygon::normal(), from polygon "<<*this<<", normal is "<<n<<" of norm()="<<nn<<std::endl;
            return MC_v3d(0,0,1);}
        n/=nn;
        return n;
    }

    MC_v3d MC_polygon::barycenter() const
    {return MC_v3d_vector::sum(v)/size();}

    MC_polygon::MC_polygon()
            :MC_v3d_vector()
    {}
    MC_polygon::MC_polygon(const MC_v3d& v0,const MC_v3d& v1,const MC_v3d& v2)
            :MC_v3d_vector(v0,v1,v2)
    {}
    MC_polygon::MC_polygon(const MC_v3d& v0,const MC_v3d& v1,const MC_v3d& v2,const MC_v3d& v3)
            :MC_v3d_vector(v0,v1,v2,v3)
    {}
    MC_polygon::MC_polygon(const MC_v3d& v0,const MC_v3d& v1,const MC_v3d& v2,const MC_v3d& v3,const MC_v3d& v4)
            :MC_v3d_vector(v0,v1,v2,v3,v4)
    {}

    MC_polygon::MC_polygon(const MC_v3d_vector& vec)
            :MC_v3d_vector(vec)
    {}
    MC_polygon::MC_polygon(const std::vector <MC_v3d>& vec)
            :MC_v3d_vector(vec)
    {}
    MC_polygon::MC_polygon(const MC_polygon& p)
            :MC_v3d_vector(p)
    {v=p.v;}

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


        if(size()<=2)
        {std::cout<<"MC_polygon::closest_point(), polygon size is not ok ["<<size()<<"]"<<std::endl;exit(-1);}

        double epsilon=0.0001;
        MC_v3d n = normal();

        // get the projection on the plane
        MC_v3d projection=x-((x-v[0]).dot(n))*n;

        // get triangulation
        std::vector <MC_triangle> triangles=MC_triangle::triangulate(*this);

        // first, try to find if the projection is inside the polygon (in one of the triangle)
        for(unsigned int k=0;k<triangles.size();k++)
        {
            if(triangles[k].is_inside(projection)==1)
            {
                if(type!=0)
                    *type=0;
                return projection;
            }
        }

        // then, the projection is outside the polygon
        // find the closest edge
        MC_v3d x0,x1,closest;MC_segment s;double L=0.0,L_min=9999.9;
        int closest_edge=-1;
        for(int k=0;k<size();k++)
        {

            x0 = v[k];
            x1 = v[(k+1)%size()];
            s=MC_segment(x0,x1);

            L = s.distance_to_point(projection);
            if(L<L_min)
            {
                closest = s.closest_point(projection);
                closest_edge = k;
                L_min=L;
            }
        }

        if(type!=0)
            *type = 1;
        //look if the closest point is inside an edge, or a vertex
        for(int k=0;k<size();k++)
            if( (closest-v[k]).norm() < epsilon) //then it is close to a vertex
                if(type!=0)
                    *type=2;
        return closest;
    }

    MC_v3d_vector MC_polygon::plane_intersection(const MC_v3d& n,const MC_v3d& x0,int *type,MC_int_vector *type_edge) const
    {
        std::cout<<"enter plane intersection"<<std::endl;

        MC_v3d_vector intersections;

        // First, take every intersection with every segment

        std::vector <MC_segment> s;
        MC_v3d_vector inter;
        if(type_edge!=0)
            type_edge->resize(size());
        MC_int_vector temp_type_edge=MC_int_vector::zeros(size());

        for(int k=0;k<size();k++)
        {
            s.push_back(segment(k));
            inter.add(s[k].plane_intersection(n,x0,&temp_type_edge[k]));

            std::cout<<temp_type_edge[k]<<" :"<<inter.last()<<std::endl;
            if(type_edge!=0)
                (*type_edge)[k]=temp_type_edge[k];
        }
//std::cout<<"============"<<std::endl;

        // look if there is no intersection that all
        int is_intersection=0;
        for(int k=0;is_intersection==0 && k<size();k++)
            if(temp_type_edge[k]!=0)
                is_intersection=1;
        if(is_intersection==0)//no intersection
        {
            if(type!=0)
                *type=0;
            return intersections;
        }


        // look if there is only one vertex in common
        int is_two_intersection=0;
        int k_vertex=-1;
        for(int k=0;is_two_intersection==0 && k<size();k++)
        {
            if( temp_type_edge[k]==1 || temp_type_edge[k]==2) is_two_intersection=1;
            if( temp_type_edge[k]==2) k_vertex=k;
        }

        if(is_two_intersection==0)
        {
            if(k_vertex==-1)
            {
                //printf("Error in Polygon in plane_intersection, types are not corrects\n");
                if(type!=0)
                    *type=0;
                return intersections;exit(-1);}
            if(type!=0)
                *type=2;
            intersections.add(v[k_vertex]); return intersections;
        }



        // there is two intersections
        int k_1=0,k_2=0;
        for(k_1=0;k_1<size();k_1++){
            for(k_2=k_1+1;k_2<size();k_2++){
                if( (temp_type_edge[k_1]==1 || temp_type_edge[k_1]==2) &&
                    (temp_type_edge[k_2]==1 || temp_type_edge[k_2]==2) &&
                    (inter[k_1]!=inter[k_2])
                    )
                {
                    if(type!=0)
                        *type=1;
                    intersections.add(inter[k_1]);
                    intersections.add(inter[k_2]);

                    //std::cout<<intersections<<std::endl;

                    return intersections;
                }
            }
        }


        //else probleme
        if(type!=0)
            *type=0;
        return intersections;

    }


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

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

    MC_polygon MC_polygon::half_space_intersection(const MC_v3d& n,const MC_v3d& x0,int *type) const
    {
        //first take the intersection of the polygon and the plane
        MC_v3d_vector intersection;
        int type_edges=0;
        intersection = plane_intersection(n,x0,&type_edges);



        //check if we need to cut or not the polygon
        MC_polygon new_polygon;
        if(type_edges!=1) // no change or just destroy the polygon
        {
            if(((*this)[0]-x0).dot(n)<0)//check the orientation of the polygon to the plane
            {
                //negative so there is no more polygons
                if(type!=0)
                    *type=2;
                new_polygon.add(MC_v3d(-1.0,-1.0,-1.0));
                return new_polygon;
            }
            else
            {
                // the polygon is unchanged
                if(type!=0)
                    *type=0;
                return *this;
            }
        }

        //else there is a modification to bring to the polygon

        if(type!=0)
            *type=1;
        MC_segment s(intersection[0],intersection[1]);
        MC_v3d old_normal=normal();

        //first add every points in the correct half space
        //(just to be sure that we know the vertex before the line)

        //find the first vertex in half space
        int k=0;
        while(((*this)[k]-x0).dot(n)<0)
            k++;
        int k_2=k;

        int added_line=0;
        do
        {
            if( ((*this)[k_2]-x0).dot(n)>0 )
                new_polygon.add((*this)[k_2]);
            else if(added_line==0)
            {
                // check if the line is in the right sense (we are sure that (*this[k2-1] already exists))
                if( ((s[1]-s[0]).cross(s[0]-(*this)[(k_2-1)<0?size()-1:k_2-1])).dot(old_normal)<0)
                {new_polygon.add(s[0]);new_polygon.add(s[1]);}
                else{new_polygon.add(s[1]);new_polygon.add(s[0]);}
                added_line=1;
            }
            k_2=(k_2+1)%size();
        }while(k_2!=k);

        return new_polygon;
    }

    MC_v3d MC_polygon::segment_intersection(const MC_segment& s,int *type) const
    {
        std::vector <MC_triangle> v_t = MC_triangle::triangulate(*this);
        int N=v_t.size();
        for(int k=0;k<N;k++)
        {
            int temp_type=0;
            MC_v3d i=v_t[k].segment_intersection(s,&temp_type);

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

            if(temp_type==1 || temp_type==2 || temp_type==3)
                if(is_inside(i)==1)
                    return i;
        }

        if(type!=0)
            *type=0;
        return MC_v3d(0,0,0);
    }

    bool MC_polygon::is_inside(const MC_v3d& x) const
    {
        std::vector <MC_triangle> t = MC_triangle::triangulate(*this);
        int N=t.size();
        for(int k=0;k<N;k++)
            if(t[k].is_inside(x)==true)
                return true;
        return false;
    }
    bool MC_polygon::is_outside(const MC_v3d& _x) const
    {return !is_inside(_x);}

    std::pair <MC_v3d,std::pair<int,std::pair<int,double> > > MC_polygon::projected_direction(const MC_v3d& c0,const MC_v3d& _d,const MC_int_vector& forbidden_edge) const
    {


        int type=0,sub_type=0;double relative=0.0;

        MC_v3d d=_d.normalized();
        double epsilon=0.00001;

        //first get the polygon in 2D
        MC_matrix R=MC_matrix::rotation_axis_to_axis(normal(),MC_v3d(0,0,1));

        bool is_intersected=false;
        MC_v3d intersection;

        MC_polygon poly_2d=R*(*this);
        MC_v3d c0_2d=MC_v3d(R*c0);
        MC_v3d d_2d=MC_v3d(R*d);

        // direction equation
        // c0+t*d => ax+by=c
        double a=d_2d[1],b=-d_2d[0],c=a*c0_2d[0]+b*c0_2d[1];


        int N=size();
        for(int k_seg=0;k_seg<N;k_seg++)
        {
            if(forbidden_edge.find(k_seg)==-1)
            {

                MC_v3d c_x0=poly_2d[k_seg];
                MC_v3d c_x1=poly_2d[(k_seg+1)%N];

                double n_segment=(c_x0-c_x1).norm();

                MC_v3d d_c=c_x1-c_x0;

                double ra=d_c[1],rb=-d_c[0],rc=ra*c_x0[0]+rb*c_x0[1];
                double det=rb*a-b*ra;

                if(abs(det)>epsilon)//parallel
                {
                    double xi=(rb*c-rc*b)/det;
                    double yi=(rc*a-ra*c)/det;

                    //if (xi,yi) in the positive sens
                    if( (xi-c0_2d[0])*d_2d[0]+(yi-c0_2d[1])*d_2d[1] >0 )
                    {

                        // check if inside
                        double proj=(xi-c_x0[0])*(c_x1[0]-c_x0[0])+(yi-c_x0[1])*(c_x1[1]-c_x0[1]);
                        double n_xi=pow((xi-c_x0[0])*(xi-c_x0[0])+(yi-c_x0[1])*(yi-c_x0[1]),0.5);

                        if(proj>=epsilon && n_xi<=n_segment)
                        {


                            if(is_intersected==false || (is_intersected==true && (intersection-c0).norm()>
                                                         pow((xi-c0_2d[0])*(xi-c0_2d[0])+(yi-c0_2d[1])*(yi-c0_2d[0]),0.5) ) )
                            {
                                is_intersected = true;
                                intersection = MC_v3d(R.inverted()*MC_v3d(xi,yi,c0_2d[2]));


                                if(proj<epsilon) //vertex 0
                                {type=0; sub_type=k_seg; relative=0.0;}
                                else if(abs(proj-1)<epsilon) //vertex 1
                                {type=0; sub_type=(k_seg+1)%N; relative=0.0;}
                                else //edge
                                {type=1; sub_type=k_seg;
                                    if(n_segment<epsilon)
                                        relative=0.0;
                                    else
                                        relative=n_xi/n_segment;
                                }
                            }
                        }
                    }
                }
            }
        }

        if(is_intersected==false)//no intersection (point outside)
        {intersection=MC_v3d(-1,-1,-1);type=-1;sub_type=-1;relative=0.0;}


        std::pair <MC_v3d,std::pair<int,std::pair<int,double> > > res;
        res.first=intersection;
        res.second.first=type;
        res.second.second.first=sub_type;
        res.second.second.second=relative;
        return res;

    }

    MC_double_vector MC_polygon::barycentric_coordinates(const MC_v3d& p) const
    {
        int N=size();
        MC_double_vector bar=MC_double_vector::zeros(N);

        //check the position
        MC_int_vector pos=position(p);

        if(pos[0]==3)
        {std::cout<<"Warning in Polygon::barycentric_coordinates("<<p<<") is not in polygon"<<std::endl;return bar;}
        if(pos[0]==2) // vertex
        {bar[pos[1]]=1.0;return bar;}
        if(pos[0]==1) // edge
        {
            double alpha=segment(pos[1]).relative_position(p);
            bar[pos[1]] = 1-alpha;
            bar[(pos[1]+1)%N] = alpha;

            return bar;
        }

        //else inside the polygon
        for(int k=0;k<N;k++)
        {
            int kp=k-1<0?N-1:k-1;
            int kn=(k+1)%N;

            double dist2=(p-v[k]).dot(p-v[k]);
            bar[k]=(MC_v3d::cotan(p-v[k],v[kp]-v[k])+MC_v3d::cotan(p-v[k],v[kn]-v[k]))/dist2;
        }

        bar/=MC_double_vector::sum(bar);
        return bar;
    }


    MC_int_vector MC_polygon::position(const MC_v3d& x) const
    {
        double epsilon=0.0001;
        int N=size();
        for(int k=0;k<N;k++)
            if( (x-v[k]).norm()<epsilon )
                return MC_int_vector (2,k);
        if(is_outside(x)==1)
            return MC_int_vector (3,-1);
        for(int k=0;k<N;k++)
            if(segment(k).is_belonging(x)==1)
                return MC_int_vector (1,k);

        return MC_int_vector (0,-1);
    }

    std::vector <MC_polygon> MC_polygon::subdivide_mid_edge_no_constraint() const
    {
        // the new midpoints
        MC_v3d_vector mid_points;
        // as much mid point as there is vertices in the polygon
        mid_points.resize(size());

        int k=0,N_mid_points=size();
        for(k=0;k<N_mid_points;k++)
            mid_points[k]=0.5*(v[k]+v[(k+1)%N_mid_points]);

        // the new polygons
        std::vector <MC_polygon> new_polygon;
        // size of new_polygon = (old)N_vertex + 1
        new_polygon.resize(N_mid_points+1);

        // all triangles excepted the last on linking the mid_points
        for(k=0;k<N_mid_points;k++)
        {
            new_polygon[k]=MC_polygon(mid_points[k],v[k],mid_points[k-1>=0?k-1:N_mid_points-1]);
        }
        //last one linking every mid_points
        for(k=0;k<N_mid_points;k++)
            new_polygon[N_mid_points].add(mid_points[k]);


        return new_polygon;
    }
    std::vector <MC_polygon> MC_polygon::subdivide_barycenter_mid_edge_no_constraint() const
    {

        MC_v3d barycenter;

        // the new midpoints
        MC_v3d_vector mid_points;

        // as much mid point as there is vertices in the polygon
        mid_points.resize(size());

        int k=0,N_mid_points=size();
        for(k=0;k<N_mid_points;k++)
        {
            mid_points[k]=0.5*(v[k]+v[(k+1)%N_mid_points]);
            barycenter+=v[k];
        }
        barycenter/=double(N_mid_points);

        std::vector <MC_polygon> new_polygon(N_mid_points);

        for(k=0;k<N_mid_points;k++)
        {
            new_polygon[k].add(v[k]);
            new_polygon[k].add(mid_points[k]);
            new_polygon[k].add(barycenter);
            new_polygon[k].add(mid_points[k-1>=0?k-1:N_mid_points-1]);
        }

        return new_polygon;

    }


    std::pair<std::vector <MC_polygon>,std::pair<MC_mesh_index_vector,std::pair<MC_int_vector,MC_int_vector> > > MC_polygon::subdivide_mid_edge(const MC_int_vector& constraint_edges) const
    {
        // in case there is no constraint
        if(constraint_edges.size()==0)
        {
            MC_mesh_index_vector m;
            MC_int_vector extra_index=MC_int_vector::linspace(size(),2*size()-1);
            MC_int_vector edge=MC_int_vector::linspace(0,size()-1);

            m.point_set()=v;
            int N=size();
            for(int k=0;k<N;++k)
            {
                MC_v3d x_middle=0.5*(v[k]+v[(k+1)%N]);
                m.point_set().add(x_middle);
            }
            // all triangles excepted the last on linking the mid_points
            for(int k=0;k<N;++k)
                m.connectivity().add(MC_int_vector( (k+1)%N,N+(k+N+1)%N,N+k ));
            //last one linking every mid_points
            m.connectivity().add(MC_int_vector::linspace(0,N-1)+N);


            return std::make_pair(m.get_polygon(),std::make_pair(m,std::make_pair(extra_index,edge) ) );
        }
        if(static_cast<int>(constraint_edges.to_set().size())==size()-1) //only one non subdivided edge
        {
            //find the non subdivided edge
            MC_int_pair non_subdivided_edge;
            int N=size();
            for(int k=0;k<N;++k)
                if(constraint_edges.find(k)==-1)
                    non_subdivided_edge=MC_int_pair(k,(k+1)%N);
            MC_v3d x_middle=0.5*(v[non_subdivided_edge[0]]+v[non_subdivided_edge[1]]);

            MC_int_vector extra_index(size());
            MC_int_vector edge=non_subdivided_edge[0];

            MC_mesh_index_vector m;
            m.point_set()=v;
            m.point_set().add(x_middle);

            std::vector <MC_polygon> v_polygon;
            for(int k=0;k<N;++k)
            {
                if(k!=non_subdivided_edge[0])
                {
                    MC_int_pair edge(k,(k+1)%N);
                    v_polygon.push_back(MC_polygon(v[edge[0]],v[edge[1]],x_middle));

                    m.connectivity().add(MC_int_vector(edge[0],edge[1],size()));
                }
            }
            return std::make_pair(v_polygon,std::make_pair(m,std::make_pair(extra_index,edge) ) );
        }

        //other case, only link the subdivided part
        MC_mesh_index_vector m;
        m.point_set()=v;

        MC_int_vector new_index_contour;
        int N=size();
        int k_new=0;
        MC_int_vector extra_index;
        MC_int_vector edge;
        for(int k=0;k<N;++k)
        {
            new_index_contour.add(k);
            if(constraint_edges.find(k)==-1)
            {
                extra_index.add(N+k_new);
                edge.add(k);
                new_index_contour.add(N+k_new);

                m.point_set().add(0.5*(v[k]+v[(k+1)%N]));
                k_new++;
            }
        }
        if(extra_index.size()==0)
        {
            m.connectivity()=MC_int_vector::linspace(0,size()-1);

            return std::make_pair(m.get_polygon(),std::make_pair(m,std::make_pair(extra_index,edge) ) );
        }


        //finish by the connectivity
        int N_new_points=m.point_set().size();
        int k_current=new_index_contour.find(extra_index[0]);
        int k_current0=k_current;
        MC_int_vector temp_poly=MC_int_vector(new_index_contour[k_current]);

        k_current=(k_current+1)%N_new_points;
        while(k_current!=k_current0)
        {
            temp_poly.add(new_index_contour[k_current]);
            if(extra_index.find(new_index_contour[k_current])!=-1)
            {
                m.connectivity().add(temp_poly);
                temp_poly.resize(0);
                temp_poly.add(new_index_contour[k_current]);
            }
            k_current=(k_current+1)%N_new_points;
        }
        temp_poly.add(new_index_contour[k_current0]);
        m.connectivity().add(temp_poly);

        if(extra_index.size()>=3)//fill the middle
            m.connectivity().add(extra_index);

        return std::make_pair(m.get_polygon(),std::make_pair(m,std::make_pair(extra_index,edge) ) );
    }
    std::pair<std::vector <MC_polygon>,std::pair<MC_mesh_index_vector,std::pair<MC_int_vector,MC_int_vector> > > MC_polygon::subdivide_barycenter_mid_edge(const MC_int_vector& constraint_edges) const
    {

        MC_mesh_index_vector m;
        MC_int_vector extra_index;
        MC_int_vector edge;

        // in case there is no constraint
        if(constraint_edges.size()==0)
        {
            extra_index=MC_int_vector::linspace(size()+1,2*size());
            edge=MC_int_vector::linspace(0,size()-1);

            m.point_set()=MC_v3d_vector(v)<<barycenter();
            int N=size();
            for(int k=0;k<N;++k)
            {
                MC_v3d x_middle=0.5*(v[k]+v[(k+1)%N]);
                m.point_set().add(x_middle);
                m.connectivity().add(MC_int_vector(N+1+(k+N+1)%N,N,N+1+k,(k+1)%N));
            }

            return std::make_pair(m.get_polygon(),std::make_pair(m,std::make_pair(extra_index,edge) ) );
        }



        MC_v3d bar=barycenter();
        m.point_set().add(v).add(bar);
        int N=size();

        MC_v3d_vector mid_point;
        MC_int_vector map_index_middle=MC_int_vector::zeros(N);int counter=0;
        for(int k=0;k<N;k++)
        {
            if(constraint_edges.find(k)==-1)
            {
                MC_v3d x=0.5*(v[k]+v[(k+1)%N]);
                mid_point.add(x);
                m.point_set().add(x);
                extra_index.add(N+1+counter);
                edge.add(k);

                map_index_middle[k]=counter++;

            }
            else
            {
                mid_point.add(MC_v3d(-1,-1,-1));
                map_index_middle[k]=-1;
            }
        }

        MC_polygon temp;
        for(int k=0;k<N;k++)
        {
            //classical quads
            if(constraint_edges.find(k)==-1 && constraint_edges.find((k-1)<0?N-1:k-1)==-1)
                m.connectivity().add(MC_int_vector(k,N+1+map_index_middle[k],N,N+1+map_index_middle[k-1>=0?k-1:map_index_middle.size()-1]));
            // previous ok but new one is a boundary
            else if(constraint_edges.find(k)!=-1 && constraint_edges.find((k-1)<0?N-1:k-1)==-1)
            {
                m.connectivity().add(MC_int_vector(k,N,N+1+map_index_middle[k-1>=0?k-1:map_index_middle.size()-1]));
                m.connectivity().add(MC_int_vector(k,(k+1)%N,N));
            }
            // previous is boundary, but next is ok
            else if(constraint_edges.find(k)==-1 && constraint_edges.find((k-1)<0?N-1:k-1)!=-1)
                m.connectivity().add(MC_int_vector(k,N+1+map_index_middle[k],N));
            // not ok, just triangle
            else if(constraint_edges.find(k)!=-1 && constraint_edges.find((k-1)<0?N-1:k-1)!=-1)
                m.connectivity().add(MC_int_vector(k,(k+1)%N,N));
        }
        return std::make_pair(m.get_polygon(),std::make_pair(m,std::make_pair(extra_index,edge) ) );

    }
    std::pair<std::vector <MC_polygon>,std::pair<MC_mesh_index_vector,std::pair<MC_int_vector,MC_int_vector> > > MC_polygon::subdivide_mixed_mid_edge(const MC_int_vector& constraint_edges) const
    {
        if(size()<3)
        {std::cout<<"Error in MC_polygon::subdivide_mixed_mid_edge(), polygon has size "<<size()<<std::endl;exit(-1);}
        if(size()==3)
            return subdivide_mid_edge(constraint_edges);
        else
            return subdivide_barycenter_mid_edge(constraint_edges);
    }

    bool MC_polygon::is_degenerated() const
    {
        std::set <MC_v3d,MC_v3d_less> set_poly;
        int N=size();
        for(int k=0;k<N;k++)
            if(set_poly.insert(v[k]).second==false)
                return true;
        return false;
    }

    std::pair<MC_polygon,std::pair<bool,bool> > MC_polygon::undegenerated() const
    {
        std::set <MC_v3d,MC_v3d_less> set_poly;
        int N=size();
        std::pair <MC_polygon,std::pair <bool,bool> > new_pol(MC_polygon(),std::pair<bool,bool>(true,true));
        for(int k=0;k<N;k++)
            if(set_poly.insert(v[k]).second==true)
                new_pol.first.add(v[k]);
        else
            new_pol.second.first=false;

        if(new_pol.first.size()>=3)
            new_pol.second.second=true;
        else
            new_pol.second.second=false;


        return new_pol;
    }








}