MC_curve.cpp

Go to the documentation of this file.

#include <MC_v4d.hpp>
#include <MC_int_vector.hpp>
#include <MC_matrix.hpp>
#include <MC_int_pair_unique.hpp>
#include <MC_segment.hpp>
#include <MC_double_vector_vector.hpp>
#include <MC_string_helper.hpp>

#include <fstream>

#include <MC_curve.hpp>
#include <map>



namespace mesh_conv
{


    MC_curve::MC_curve():MC_v3d_vector(){}
    MC_curve::MC_curve(const MC_v3d& v0):MC_v3d_vector(v0){}
    MC_curve::MC_curve(const MC_v3d& v0,const MC_v3d& v1):MC_v3d_vector(v0,v1){}
    MC_curve::MC_curve(const MC_v3d& v0,const MC_v3d& v1,const MC_v3d& v2):MC_v3d_vector(v0,v1,v2){}
    MC_curve::MC_curve(const MC_v3d& v0,const MC_v3d& v1,const MC_v3d& v2,const MC_v3d& v3):MC_v3d_vector(v0,v1,v2,v3){}
    MC_curve::MC_curve(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_curve::MC_curve(const MC_v3d_vector& vec)
            :MC_v3d_vector(vec){}
    MC_curve::MC_curve(const MC_curve& vec)
            :MC_v3d_vector(static_cast<MC_v3d_vector>(vec)){}
    MC_curve::MC_curve(const MC_double_vector& x_vector,const MC_double_vector& y_vector,const MC_double_vector& z_vector)
            :MC_v3d_vector(x_vector,y_vector,z_vector){}


    MC_curve MC_curve::read_lin(std::istream& stream)
    {
        MC_curve curve;

        MC_v3d_vector coords;
        std::string buffer;

        while(stream.good()==true)
        {
            stream>>buffer;
            if(stream.good()==true)
            {
                if(buffer.compare("v")==0)
                {
                    MC_v3d x; stream>>x[0]; stream>>x[1]; stream>>x[2];
                    coords.add(x);
                }
                else if(buffer.compare("s")==0)
                {
                    int u0=-1,u1=-1;
                    stream>>u0; stream>>u1;
                    if(u0-1<0 || u1-1<0 || u0-1>=coords.size() || u1-1>=coords.size())
                    {std::cout<<"Error in MC_curve::load_lin(stream), size is not correct ("<<u0<<","<<u1<<") ; and coords size is "<<coords.size()<<std::endl;exit(-1);}

                    if(curve.size()>0 && (coords[u0-1]-curve[curve.size()-1]).norm()>0.00001)
                    {std::cout<<"Something strange in Curve_3D::load_lin(stream), doesn't seems to be one line"<<std::endl;}

                    if(curve.size()==0)
                        curve.add(coords[u0-1]);

                    curve.add(coords[u1-1]);

                }
            }
        }

        if(curve.size()<=0)
        {std::cout<<"Something strange in Curve_3D::load_lin(stream), line is empty"<<std::endl;}

        return curve;
    }

//    MC_curve MC_curve::read_vect_curve(std::istream& stream)
//    {
//
//        MC_curve curve;
//
//        MC_v3d_vector coords;
//        std::string buffer;
//        const int SIZE_BUFFER=2048;
//        char c_buffer[SIZE_BUFFER]={'\0'};
//
//        while(stream.good()==true)
//        {
//            stream>>buffer;
//            if(buffer.compare("VECT")==0)
//            {
//
//                while(stream.good()==true)
//                {
//                    stream>>buffer;
//                    if(stream.good()==true)
//                    {
//
//                        // jump comments
//                        if(buffer[0]=='#')
//                            stream.getline(c_buffer,SIZE_BUFFER);
//
//                        // check if only one line
//                        else if(buffer[0]!='1')
//                        {std::cout<<"Error in MC_curve::load_vect(stream), more than one line "<<buffer<<std::endl;exit(-1);}
//                        else
//                        {
//                            // get number of vertices color ...
//                            int N_total_vertices=0,N_total_color=0,N_vertices,N_color=0;
//                            stream>>N_total_vertices;
//                            stream>>N_total_color;
//                            stream>>N_vertices;
//                            stream>>N_color;
//
//
//                            if(N_total_vertices!=N_vertices)
//                            {std::cout<<"Error in Curve_3D::load_vect(stream), N_total_vertices!=N_vertices ("<<N_total_vertices<<","<<N_vertices<<") for only one line"<<std::endl;exit(-1);}
//
//                            for(int k=0;k<N_vertices;k++)
//                            {
//                                MC_v3d x;
//                                for(int k_dim=0;k_dim<3;k_dim++)
//                                    stream>>x[k_dim];
//                                curve.add(x);
//                            }
//                            for(int k=0;k<N_color;k++)
//                            {
//                                MC_v4d useless_color;
//                                for(int k_dim=0;k_dim<4;k_dim++)
//                                    stream>>useless_color[k_dim];
//                            }
//                        }
//                    }
//                }
//
//            }
//        }
//
//
//        return curve;
//    }

    double MC_curve::length() const
    {
        double L=0.0;
        int N=size();
        for(int k=0;k<N-1;++k)
            L += (v[k+1]-v[k]).norm();
        return L;
    }

    std::pair <MC_curve,bool> MC_curve::unclosed() const
    {
        if(size()<2)
        {std::cout<<"Error in MC_curve::unclosed(), size()="<<size()<<std::endl;exit(-1);}

        MC_curve res;
        if(last()==first())
            return std::pair <MC_curve,bool> ((*this)(MC_int_vector::linspace(0,size()-2)),true);
        else
            return std::pair <MC_curve,bool> (*this,false);
    }
    std::pair <MC_curve,bool> MC_curve::closed() const
    {
        if(size()<2)
        {std::cout<<"Error in MC_curve::closed(), size()="<<size()<<std::endl;exit(-1);}

        MC_curve res;
        if(last()!=first())
            return std::pair <MC_curve,bool> ((*this)<<first(),true);
        else
            return std::pair <MC_curve,bool> (*this,false);
    }

    MC_curve::MC_curve(const std::list<MC_v3d>& input)
    {
        resize(input.size());
        std::list <MC_v3d> :: const_iterator it,it_end=input.end();
        int k=0;
        for(it=input.begin();it!=it_end;++it,++k)
            v[k]=*it;
    }


    MC_curve MC_curve::build_circle(const double& R,const int& N,const MC_v3d& normal)
    {
#define PI 3.14159

        MC_curve circle(N);

        double epsilon=0.0001;
        if(normal.norm()<epsilon)
        {std::cout<<"Error in MC_curve::build_circle(...), normal is null"<<std::endl;exit(-1);}

        MC_v3d axis_1 = normal.cross(MC_v3d(1,0,0));
        if(axis_1.norm()<epsilon)
        {
            axis_1 = normal.cross(MC_v3d(0,1,0));
            if(axis_1.norm()<epsilon)
            {std::cout<<"Error in MC_mesh_index_vector::build_disc, something weird"<<std::endl; exit(-1);}
        }
        axis_1 = axis_1.normalized();
        MC_v3d axis_2 = (axis_1.cross(normal)).normalized();

        // vertices
        for(int k_radius=0;k_radius<N;k_radius++)
        {
            double theta = 2*PI*double(k_radius)/double(N);
            circle[k_radius]=R*sin(theta)*axis_1 + R*cos(theta)*axis_2;
        }
        return circle;
    }
    MC_curve MC_curve::build_ellipsoid(const double& R1,const double& R2,const int& N,const MC_v3d& e0,const MC_v3d& e1,const double& theta_1,const double& theta_2)
    {
        MC_curve c(N);
        for(int k=0;k<N;k++)
        {
            double theta=theta_1+k/double(N)*(theta_2-theta_1);
            c[k]=MC_v3d(R1*cos(theta),R2*sin(theta),0.0);
        }

        if(e0.norm()>0.0001 && e1.norm()>0.0001)
        {
            MC_v3d u0=e0.normalized();
            MC_v3d u1=e1.normalized();
            if( abs(u0.dot(u1))>0.0001 )
            {std::cout<<"Error in MC_curve::get_ellipsoid() e0 not perp to e1 :"<<u0<<","<<u1<<std::endl;exit(-1);}

            MC_v3d u2=u0.cross(u1);
            MC_matrix R(u0,u1,u2);
            for(int k=0;k<N;k++)
                c(k)=R*c(k);
        }

        return c;
    }
    MC_curve MC_curve::build_heart(const double& R,const int& N)
    {
        double Rx=1,Ry=1;double Lx=2*Rx,Ly=3;

        MC_curve heart;
        MC_curve circle_0,circle_1;
        MC_curve line_0,line_1;


        line_0=MC_curve::build_line(MC_v3d(0,0,0),MC_v3d(Lx,Ly,0),N);
        circle_0 = build_ellipsoid(Rx,Ry,N,MC_v3d(),MC_v3d(),0.0,PI);circle_0=circle_0+MC_v3d(Lx-Rx,Ly,0);
        circle_1 = circle_0+MC_v3d(-2*Rx,0,0);
        line_1=MC_curve::build_line(MC_v3d(-Lx,Ly,0),MC_v3d(0,0,0),N);

        heart.add(line_0);heart.add(circle_0);heart.add(circle_1);heart.add(line_1);
        heart = heart*R;

        heart=heart-heart.barycenter();

        return heart;
    }
    MC_curve MC_curve::build_line(const MC_v3d& X0,const MC_v3d& X1,const int& N)
    {
        MC_curve line=MC_curve::zeros(N);
        for(int k=0;k<N;++k)
        {
            double alpha = static_cast<double>(k)/static_cast<double>(N-1);
            line[k]=(1-alpha)*X0+alpha*X1;
        }
        return line;
    }
    MC_curve MC_curve::build_B_spline(const MC_curve& control_polygon,const int& N_subdiv)
    {
        MC_curve spline;
        MC_matrix M(4,4);
        M(0,0)=-1;M(1,0)= 3;M(2,0)=-3;M(3,0)= 1;
        M(0,1)= 3;M(1,1)=-6;M(2,1)= 0;M(3,1)= 4;
        M(0,2)=-3;M(1,2)= 3;M(2,2)= 3;M(3,2)= 1;
        M(0,3)= 1;M(1,3)= 0;M(2,3)= 0;M(3,3)= 0;
        M/=6.0;



        int k_t=0;
        int k=0;
        if(control_polygon.size()<4)
        {std::cout<<"Warning MC_curve::build_B_spline(), must be at least 4 control points"<<std::endl;return control_polygon;}

        int N=control_polygon.size();
        MC_double_vector Px,Py,Pz;
        MC_double_vector T=MC_double_vector::zeros(4);
        double t=0.0;
        int k_section=0;

        double Sx=0.0,Sy=0.0,Sz=0.0;

        MC_double_vector temp_x,temp_y,temp_z;
        for(k_section=0;k_section<N-3;k_section++)
        {
            // polygon
            Px.resize(0);Py.resize(0);Pz.resize(0);

            Px.add(control_polygon[k_section+0][0]);
            Px.add(control_polygon[k_section+1][0]);
            Px.add(control_polygon[k_section+2][0]);
            Px.add(control_polygon[k_section+3][0]);

            Py.add(control_polygon[k_section+0][1]);
            Py.add(control_polygon[k_section+1][1]);
            Py.add(control_polygon[k_section+2][1]);
            Py.add(control_polygon[k_section+3][1]);

            Pz.add(control_polygon[k_section+0][2]);
            Pz.add(control_polygon[k_section+1][2]);
            Pz.add(control_polygon[k_section+2][2]);
            Pz.add(control_polygon[k_section+3][2]);


            //kernel
            temp_x = M*Px;
            temp_y = M*Py;
            temp_z = M*Pz;

            for(k_t=0;k_t<N_subdiv;k_t++)
            {
                t = double(k_t)/double(N_subdiv);
                // time
                T[0]=t*t*t;T[1]=t*t;T[2]=t;T[3]=1;

                // spline
                Sx=0;Sy=0;Sz=0;
                for(k=0;k<4;k++)
                {
                    Sx += T[k]*temp_x[k];
                    Sy += T[k]*temp_y[k];
                    Sz += T[k]*temp_z[k];
                }

                spline.add(MC_v3d(Sx,Sy,Sz));
            }
        }

        return spline;
    }
    MC_v3d_vector MC_curve::diff_forward() const
    {

        int N=size();

        //special case
        if(N==1)
            return MC_v3d_vector(MC_v3d(0,0,0));
        if(N==2)
            return MC_v3d_vector(v[1]-v[0],v[1]-v[0]);

        MC_v3d_vector diff(N);;
        for(int k=0;k<N-1;++k)
            diff[k]=v[k+1]-v[k];

        //last one
        diff[N-1]=v[N-1]-v[N-2];

        return diff;
    }
    MC_v3d_vector MC_curve::diff_forward_close() const
    {
        int N=size();

        //special case
        if(N==1)
            return MC_v3d_vector(MC_v3d(0,0,0));
        if(N==2)
            return MC_v3d_vector(v[1]-v[0],v[0]-v[1]);

        MC_v3d_vector diff(N);;
        for(int k=0;k<N-1;++k)
            diff[k]=v[k+1]-v[k];

        //last one
        diff[N-1]=v[0]-v[N-1];

        return diff;
    }


    std::pair<MC_curve,std::pair<std::vector<std::pair <MC_int_pair_unique,MC_double_vector> >,MC_double_vector> > MC_curve::sample_linear(const double& d_L) const
    {

        //tune the sampling to fine a correct number of iteration
        double total_length=length();
        double N_theory = total_length/d_L;
        int N_real=static_cast<int>(N_theory+1);
        MC_double_vector t_sampling=MC_double_vector::sample(0,1,N_real);

        //std::cout<<"t_sampling "<<N_real<<" | "<<t_sampling<<std::endl;
        //samples
        std::pair <MC_v3d_vector,std::vector<std::pair<MC_int_pair_unique,MC_double_vector> > > sampling=value(t_sampling);

        return std::make_pair(sampling.first,std::make_pair(sampling.second,t_sampling));

//
//        std::pair <MC_curve,std::vector <std::pair<MC_int_pair_unique,MC_double_vector> > >res;
//        //Curve_3D new_curve;
//
//        int N=size();
//        // resample each segment
//        for(int k=0;k<N-1;k++)
//        {
//            int k2=k+1;
//            MC_segment s((*this)(k),(*this)(k2));
//            while(k2<N-1 && s.length()<d_L)
//            {
//                k2++;
//                if(k2<N-1)
//                    s=MC_segment((*this)(k),(*this)(k2));
//            }
//
//
//            std::pair <MC_curve,MC_double_vector> temp=MC_segment((*this)(k),(*this)(k2)).linear_sampling_intervals(d_L);
//
//            //std::cout<<"kk : "<<k<<" ! "<<temp.first<<std::endl;
//            if(k==0)//add everything
//                res.first.add(temp.first);
//            else
//                res.first.add(temp.first(MC_int_vector::linspace(1,temp.first.size()-1)));
//
//
//            // now check the t values
//            std::pair <MC_int_pair_unique,MC_double_vector> t_vals;
//            t_vals.first=MC_int_pair_unique(k,k2);
//            for(int k3=0;k3<temp.first.size();k3++)
//            {
//                if((k3==0 && k==0) || (k3!=0 && k!=0))
//                    t_vals.second.add(temp.second[k3]);
//            }
//            res.second.push_back(t_vals);
//
//            k=k2-1;
//        }
//        return res;
    }

    std::vector <MC_segment> MC_curve::segment() const
    {
        int N=size();
        std::vector <MC_segment> seg(N-1);
        for(int k=0;k<N-1;++k)
            seg[k]=MC_segment((*this)(k),(*this)(k+1));
        return seg;
    }
    MC_double_vector MC_curve::relativ_position() const
    {
        if(size()==0)
        {std::cout<<"Error in MC_curve::relativ_position(), size=0"<<std::endl;exit(-1);}
        if(size()==1)
            return MC_double_vector(0.0);
        if(size()==2)
            return MC_double_vector(0.0,1.0);

        int N=size();
        std::vector <MC_segment> seg=segment();
        double length_total=length();
        MC_double_vector t_relativ=MC_double_vector::zeros(N);t_relativ[0]=0.0;
        double length_cumulativ=0.0;
        for(int k=0;k<N-1;++k)
        {
            length_cumulativ+=seg[k].length();
            t_relativ[k+1]=length_cumulativ/length_total;
        }
        return t_relativ;
    }
    std::pair <MC_v3d,std::pair<MC_int_pair_unique,double> > MC_curve::value(const double& t) const
    {
        double epsilon=0.00001;
        int N=size();

        if(t<0 || t>1)
        {std::cout<<"Error in MC_curve::value("<<t<<"), value must be in [0,1]"<<std::endl;exit(-1);}

        if(size()==0)
        {std::cout<<"Error in MC_curve::value("<<t<<"), size=0"<<std::endl;exit(-1);}
        if(size()==1)
            return std::make_pair((*this)[0],std::make_pair(MC_int_pair_unique(0,0),0));

        //special cases
        if(fabs(t)<epsilon)
            return std::make_pair(first(),std::make_pair(MC_int_pair_unique(0,1),0));
        if(fabs(t-1.0)<epsilon)
            return std::make_pair(last(),std::make_pair(MC_int_pair_unique(N-2,N-1),1));

        //get the relative position of the vertices
        MC_double_vector vt=relativ_position();

        //get the surrounding vertices
        int closest_inferior=-1;
        for(int k=1;closest_inferior==-1 && k<N;++k)
        {
            if(fabs(t-vt[k])<epsilon)
                return std::make_pair((*this)[k],std::make_pair(MC_int_pair_unique(k-1,k),1.0));
            if(t>vt[k-1] && t<vt[k])
                closest_inferior=k-1;
        }

        // now linear interpolation
        double t0=vt[closest_inferior];
        double t1=vt[closest_inferior+1];
        MC_v3d y0=(*this)[closest_inferior];
        MC_v3d y1=(*this)[closest_inferior+1];

        double alpha=(t-t0)/(t1-t0);
        MC_v3d y= (1-alpha)*y0+alpha*y1;

        return std::make_pair(y,std::make_pair(MC_int_pair_unique(closest_inferior,closest_inferior+1),alpha));

    }
    std::pair <MC_v3d_vector,std::vector<std::pair<MC_int_pair_unique,MC_double_vector> > > MC_curve::value(const MC_double_vector& t) const
    {
        std::vector <std::pair <MC_int_pair_unique,MC_double_vector> > index_segment(t.size());
        MC_v3d_vector y=MC_v3d_vector::zeros(t.size());

        int N=t.size();
        for(int k=0;k<N;++k)
        {
            std::pair <MC_v3d,std::pair <MC_int_pair_unique,double> > temp=value(t[k]);
            y[k]=temp.first;
            index_segment[k].first=temp.second.first;
            index_segment[k].second=temp.second.second;
        }
        return make_pair(y,index_segment);
    }

    std::pair <MC_v3d,double> MC_curve::value(const int& k_segment,const double& _relative_position) const
    {
        int N=size();
        if(k_segment>=N-1 || k_segment<0)
        {std::cout<<"Error in MC_curve::value("<<k_segment<<","<<_relative_position<<"). k_segment is not correct for a curve of size "<<N<<std::endl;exit(-1);}
        MC_segment s((*this)[k_segment],(*this)[k_segment+1]);
        MC_v3d x=s.value(_relative_position);

        MC_double_vector relativ=relativ_position();
        double r=relativ[k_segment]*(1-_relative_position)+relativ[k_segment+1]*_relative_position;

        return std::make_pair(x,r);
    }
    std::pair <MC_v3d_vector,MC_double_vector> MC_curve::value(const MC_int_vector& k_segment,const MC_double_vector& _relative_position) const
    {
        std::vector <MC_segment> s=segment();
        MC_double_vector relativ=relativ_position();

        int N_in=k_segment.size();
        int N_curve=size();
        if(k_segment.size()!=_relative_position.size())
        {std::cout<<"Error in MC_curve::value("<<k_segment<<","<<_relative_position<<"), size are not compatible ("<<k_segment.size()<<"!="<<_relative_position.size()<<")"<<std::endl;exit(-1);}

        MC_v3d_vector x=MC_v3d_vector::zeros(N_in);
        MC_double_vector t=MC_double_vector::zeros(N_in);
        for(int k=0;k<N_in;++k)
        {
            int current_seg=k_segment[k];
            if(current_seg==N_curve-1 && std::fabs(_relative_position[k])<1e-4)
            {x[k]=last(); t[k]=1.0;}
            else if(current_seg<0 || current_seg>=N_curve-1)
            {std::cout<<"Error in MC_curve::value(), at entry "<<k<<"/"<<N_in<<", value is not correct: "<<current_seg<<", for curve size="<<N_curve<<std::endl;exit(-1);}
            else
            {
                x[k]=s[current_seg].value(_relative_position[k]);
                t[k]=relativ[current_seg]*(1-_relative_position[k])+relativ[current_seg+1]*_relative_position[k];
            }
        }

        return std::make_pair(x,t);
    }

    MC_curve MC_curve::build_square(const int& N)
    {
        MC_curve square;
        for(int k=0;k<N-1;++k)
        {
            double alpha=static_cast<double>(k)/static_cast<double>(N-1);
            MC_v3d x=MC_v3d(alpha,0,0);
            square.add(x);
        }
        for(int k=0;k<N-1;++k)
        {
            double alpha=static_cast<double>(k)/static_cast<double>(N-1);
            MC_v3d x=MC_v3d(1,alpha,0);
            square.add(x);
        }
        for(int k=0;k<N-1;++k)
        {
            double alpha=static_cast<double>(k)/static_cast<double>(N-1);
            MC_v3d x=MC_v3d(1-alpha,1,0);
            square.add(x);
        }
        for(int k=0;k<N-1;++k)
        {
            double alpha=static_cast<double>(k)/static_cast<double>(N-1);
            MC_v3d x=MC_v3d(0,1-alpha,0);
            square.add(x);
        }
        return square;
    }


    std::pair <MC_v3d,std::pair<int,double> > MC_curve::closest_point(const MC_v3d& x) const
    {

        std::vector <MC_segment> s=segment();
        if(s.size()<=1)
        {
            std::cout<<"Warning, MC_curve::closest_point, curve is empty"<<std::endl;
            return std::make_pair(MC_v3d(-1,-1,-1),std::make_pair(-1,-1));
        }
        MC_v3d saved_closest=s[0].closest_point(x);
        double min_dist=(saved_closest-x).norm();
        std::pair<int,double> param_to_save;
        param_to_save.first=0;
        param_to_save.second=s[0].relative_position(saved_closest);

        for(unsigned int k=1;k<s.size();++k)
        {
            MC_v3d current_closest=s[k].closest_point(x);
            //std::cout<< k <<" "<<s[k]<<std::endl;
            double current_distance=(current_closest-x).norm();
            if(current_distance<min_dist)
            {
                min_dist=current_distance;

                saved_closest=current_closest;
                param_to_save.first=k;
                param_to_save.second=s[k].relative_position(current_closest);

            }
        }

        return std::make_pair(saved_closest,param_to_save);
    }

    MC_curve::MC_curve(const MC_segment& seg):MC_v3d_vector(seg[0],seg[1])
    {}

    MC_v3d MC_curve::barycenter() const
    {
        MC_v3d barycenter;
        double total_length=length();

        for(int k=0,N=size();k<N-1;++k)
        {
            MC_segment seg=MC_segment((*this)[k],(*this)[k+1]);
            double L=seg.length();

            MC_v3d p=0.5*(seg[1]+seg[0]);

            barycenter += L*p;
        }
        barycenter/=total_length;
        return barycenter;
    }

    MC_curve MC_curve::skinning(const MC_curve& skeleton_old,const MC_curve& skeleton_new) const
    {
        if(skeleton_old.size()!=skeleton_new.size())
        {std::cout<<"Error in MC_curve::skinning, skeleton does not have the same size"<<std::endl;exit(-1);}


        //weights
        int N=skeleton_new.size();
        int N_current=size();

        std::vector<MC_segment> segment_old=skeleton_old.segment();
        std::vector<MC_segment> segment_new=skeleton_new.segment();

        MC_double_vector_vector w;w.resize(N_current);
        MC_v3d_vector detail=MC_v3d_vector::zeros(N_current);

        for(int k=0;k<N_current;++k)
        {
            w[k].resize(N-1);
            for(int k2=0;k2<N-1;++k2)
            {
                MC_v3d closest=segment_old[k2].closest_point(v[k]);
                w[k][k2]=(closest-v[k]).norm();
            }
            //normalize
            w[k]=w[k]/MC_double_vector::sum(w[k]);
        }


        std::pair<MC_int_vector,MC_double_vector> barycentric;
        for(int k=0;k<N_current;++k)
        {
            std::pair<MC_v3d,std::pair<int,double> > bar = skeleton_old.closest_point(v[k]);
            barycentric.first.add(bar.second.first);
            barycentric.second.add(bar.second.second);
            detail[k]=-(bar.first-v[k]);
        }


        std::vector<MC_matrix> R(N);
        for(int k=0;k<N-1;++k)
        {
            MC_v3d x0=segment_old[k][0];
            MC_v3d x1=segment_old[k][1];

            MC_v3d y0=segment_new[k][0];
            MC_v3d y1=segment_new[k][1];

            MC_v3d u0=x1-x0;
            MC_v3d u1=y1-y0;

            R[k]=MC_matrix::rotation_axis_to_axis(u0,u1);
        }
        std::vector<MC_matrix> T(N_current);
        for(int k=0;k<N_current;++k)
        {
            T[k]=MC_matrix(3);
            for(int k2=0;k2<N-1;++k2)
                T[k]+=w[k][k2]*R[k2].inverted();
        }

        MC_curve curve_new=skeleton_new.value(barycentric.first,barycentric.second).first;
        for(int k=0;k<N_current;++k)
            curve_new[k] += T[k]*detail[k];

        return curve_new;
    }

    void MC_curve::write_vect(std::ostream& stream,const std::vector<MC_curve>& v_curve,const MC_v3d_vector& _v_color)
    {
        if(stream.good()!=true)
        {std::cout<<"Error in MC_curve::write_vect(), ostream is not correct for writing"<<std::endl;exit(-1);}

        if(static_cast<unsigned int>(v_curve.size())!=static_cast<unsigned int>(_v_color.size()) && _v_color.size()!=0)
        {std::cout<<"Error in MC_curve::write_vect(), v_curve and v_color have non coherent size "<<v_curve.size()<<"!="<<_v_color.size()<<std::endl;exit(-1);}

        MC_v3d_vector v_color=_v_color;
        if(_v_color.size()==0)
            v_color=MC_v3d_vector::zeros(v_curve.size());

        unsigned int N=v_curve.size();
        unsigned int counter_vertices=0;
        for(unsigned int k=0;k<N;++k)
            counter_vertices+=v_curve[k].size();

        //header
        stream<<"VECT"<<std::endl;
        stream<<v_curve.size()<<" "<<counter_vertices<<" "<<v_color.size()<<std::endl;
        for(unsigned int k=0;k<N;++k)
            stream<<v_curve[k].size()<<" ";
        stream<<std::endl;
        for(unsigned int k=0;k<N;++k)
            stream<<"1"<<" ";
        stream<<std::endl<<std::endl;

        //vertices
        for(unsigned int k=0;k<N;++k)
        {
            const MC_curve& curve_current=v_curve[k];
            int N_current=curve_current.size();
            for(int kc=0;kc<N_current;++kc)
                stream<<curve_current[kc][0]<<" "<<curve_current[kc][1]<<" "<<curve_current[kc][2]<<"  ";
            stream<<std::endl;
        }
        stream<<std::endl;

        //color
        for(unsigned int k=0;k<N;++k)
        {
            const MC_v3d& current_color=v_color[k];
            stream<<current_color[0]<<" "<<current_color[1]<<" "<<current_color[2]<<" 1"<<std::endl;
        }

    }
    void MC_curve::write_vect(const std::string& filename,const std::vector<MC_curve>& v_curve,const MC_v3d_vector& v_color)
    {
        std::ofstream fid(filename.c_str(),std::ios::out);
        if(fid.good()!=true)
        {std::cout<<"Error in MC_curve::write_vect(), cannot open file "<<filename<<std::endl;exit(-1);}
        write_vect(fid,v_curve,v_color);
        fid.close();
    }

    void MC_curve::write_vect(const std::string& filename) const
    {
        std::vector<MC_curve> v_c;v_c.push_back(*this);
        write_vect(filename,v_c);
    }

     std::pair<std::vector<MC_curve>,MC_v3d_vector> MC_curve::read_vect(std::istream& stream)
     {
         if(stream.good()!=true)
         {std::cout<<"Error in MC_curve::read_vect(), stream is not open"<<std::endl;exit(-1);}

         std::vector<MC_curve> v_curve;
         MC_v3d_vector v_color;

         std::string buffer;
         //read header
         bool is_loop=true;
         while(stream.good() && is_loop==true)
         {
             std::getline(stream,buffer);
             if(buffer.length()>0
                && buffer[0]!='#'
                && buffer.find("VECT")!=std::string::npos)
                 is_loop=false;
         }
         if(stream.good()==false && is_loop==true)
         {std::cout<<"ERROR in MC_io_off::read_off(istream), cannot find OFF header"<<std::endl;exit(-1);}


         unsigned int N_polyline=0;
         unsigned int N_vertex=0;
         unsigned int N_color=0;

         is_loop=true;
         //read N_polyline, N_vertex, N_color
         while(stream.good() && is_loop==true)
         {
             std::getline(stream,buffer);
             if(buffer.length()>0 && buffer[0]!='#')
             {
                 std::vector <std::string> token=MC_string_converter::delete_empty(MC_string_tokenizer::tokenize(buffer));
                 if(token.size()>=3)
                 {
                     bool polyline_converted=false;
                     bool vertex_converted=false;
                     bool color_converted=false;
                     N_polyline  = MC_string_converter::value_of<unsigned int>(token[0],&polyline_converted);
                     N_vertex = MC_string_converter::value_of<unsigned int>(token[1],&vertex_converted);
                     N_color    = MC_string_converter::value_of<unsigned int>(token[2],&color_converted);

                     if(polyline_converted==true && vertex_converted==true && color_converted==true)
                         is_loop=false;
                 }
             }
         }

//         if(N_polyline<=0)
//         {
//             std::vector<MC_curve> v_empty;
//             return std::make_pair(v_empty,MC_v3d_vector());
//         }

         //read size
         MC_int_vector polyline_size=MC_int_vector::zeros(N_polyline);
         is_loop=true;
         while(stream.good() && is_loop==true)
         {
             std::getline(stream,buffer);
             if(buffer.length()>0 && buffer[0]!='#')
             {
                 std::vector <std::string> token=MC_string_converter::delete_empty(MC_string_tokenizer::tokenize(buffer));
                 if(token.size()>=N_polyline)
                 {
                     for(unsigned int k=0;k<N_polyline;++k)
                     {
                         bool polysize_converted=false;
                         polyline_size[k]=MC_string_converter::value_of<unsigned int>(token[k],&polysize_converted);
                         if(polysize_converted==true)
                             is_loop=false;
                     }
                 }
                 else
                 {std::cout<<"Error in MC_curve::read_vect(), inconsistent vect size: expected "<<N_polyline<<" entries while reading "<<buffer<<std::endl; exit(-1);}
             }
         }

         //read color
         MC_int_vector color_index=MC_int_vector::zeros(N_polyline);
         is_loop=true;
         while(stream.good() && is_loop==true)
         {
             std::getline(stream,buffer);
             if(buffer.length()>0 && buffer[0]!='#')
             {
                 std::vector <std::string> token=MC_string_converter::delete_empty(MC_string_tokenizer::tokenize(buffer));
                 if(token.size()>=N_polyline)
                 {
                     for(unsigned int k=0;k<N_polyline;++k)
                     {
                         bool color_converted=false;
                         color_index[k]=MC_string_converter::value_of<unsigned int>(token[k],&color_converted);
                         if(color_converted==true)
                             is_loop=false;
                     }
                 }
                 else
                 {std::cout<<"Error in MC_curve::read_vect(), inconsistent color size: expected "<<N_polyline<<" entries while reading "<<buffer<<std::endl; exit(-1);}
             }
         }


         //read vertices
         v_curve.resize(N_polyline);
         for(unsigned int k_poly=0;k_poly<N_polyline;++k_poly)
         {
             v_curve[k_poly].resize(polyline_size[k_poly]);

             if(stream.good()==false)
             {std::cout<<"Error in MC_curve::read_vect() at polyline"<<k_poly<<" EOF found"<<std::endl;exit(-1);}

             std::getline(stream,buffer);
             if(buffer.size()>0 && buffer[0]!='#')
             {
                 std::vector<std::string> v_token=MC_string_converter::delete_empty(MC_string_tokenizer::tokenize(buffer));
                 if(v_token.size()<3*static_cast<unsigned int>(polyline_size[k_poly]))
                 {std::cout<<"Error in MC_curve::read_vect(), something wrong at polyline "<<k_poly<<", size of token="<<v_token.size()<<", "<<buffer<<std::endl;exit(-1);}


                 for(int k_vertex=0,N_vertex=polyline_size[k_poly];k_vertex<N_vertex;++k_vertex)
                 {

                     bool converted_x0=false,converted_x1=false,converted_x2=false;
                     double x0=MC_string_converter::value_of<double>(v_token[3*k_vertex+0],&converted_x0);
                     double x1=MC_string_converter::value_of<double>(v_token[3*k_vertex+1],&converted_x1);
                     double x2=MC_string_converter::value_of<double>(v_token[3*k_vertex+2],&converted_x2);

                     if(converted_x0==false || converted_x1==false || converted_x2==false)
                     {std::cout<<"Something wrong in MC_curve::read_vect(), cannot convert at vertex "<<k_vertex<<", of polyline "<<k_poly<<" : "<<buffer<<std::endl;exit(-1);}

                     v_curve[k_poly][k_vertex]=MC_v3d(x0,x1,x2);
                 }
             }
             else
                 --k_poly;
         }

         //read color
         v_color.resize(N_polyline);
         for(unsigned int k_poly=0;k_poly<N_polyline;++k_poly)
         {
             if(color_index[k_poly]==0)
             {
                 if(k_poly>0)
                     v_color[k_poly]=v_color[k_poly-1];
                 else
                 {std::cout<<"Error in MC_curve::read_vect(), something wrong in reading the first color: no color"<<std::endl;exit(-1);}
             }
             else
             {
                 is_loop=true;
                 //read the actual entry
                 while(stream.good() && is_loop==true)
                 {
                     std::getline(stream,buffer);
                     if(buffer.length()>0 && buffer[0]!='#')
                     {
                         std::vector <std::string> token=MC_string_converter::delete_empty(MC_string_tokenizer::tokenize(buffer));
                         if(token.size()>=3)
                         {
                             bool r_converted=false,g_converted=false,b_converted=false;
                             double r=MC_string_converter::value_of<double>(token[0],&r_converted);
                             double g=MC_string_converter::value_of<double>(token[1],&g_converted);
                             double b=MC_string_converter::value_of<double>(token[2],&b_converted);


                             if(r_converted==false || g_converted==false || b_converted==false)
                             {std::cout<<"Something wrong in MC_curve::read_vect(), cannot convert at color at polyline "<<k_poly<<" : "<<buffer<<std::endl;exit(-1);}

                             v_color[k_poly]=MC_v3d(r,g,b);
                             is_loop=false;

                         }
                         else
                         {std::cout<<"Error in MC_curve::read_vect(), inconsistent color size: expected "<<3<<" entries while reading "<<buffer<<std::endl; exit(-1);}
                     }
                 }
             }
         }



         return std::make_pair(v_curve,v_color);

     }

     std::pair<std::vector<MC_curve>,MC_v3d_vector> MC_curve::read_vect(const std::string& filename)
     {
         std::ifstream ifile(filename.c_str(),std::ios::in);
         if(ifile.good()!=true)
         {std::cout<<"Error in MC_curve::read_vect("<<filename<<"), cannot open the file"<<std::endl;exit(-1);}

         std::pair<std::vector<MC_curve>,MC_v3d_vector> input=read_vect(ifile);

         ifile.close();
         return input;
     }
     MC_curve MC_curve::read_vect_curve(std::istream& stream)
     {
         std::pair<std::vector<MC_curve>,MC_v3d_vector> input=read_vect(stream);
         if(input.first.size()<1)
         {std::cout<<"Error in MC_curve::read_vect_curve(), no curve found"<<std::endl;exit(-1);}
         return input.first[0];
     }
     MC_curve MC_curve::read_vect_curve(const std::string& filename)
     {
         std::ifstream ifile(filename.c_str(),std::ios::in);
         if(ifile.good()!=true)
         {std::cout<<"Error in MC_curve::read_vect("<<filename<<"), cannot open the file"<<std::endl;exit(-1);}

         MC_curve input=read_vect_curve(ifile);

         ifile.close();
         return input;
     }
     std::pair<MC_v3d_vector,MC_double_vector> MC_curve::value_t(const MC_double_vector t) const
     {
         MC_int_vector k_seg=MC_int_vector::zeros(t.size());
         MC_double_vector t_seg=MC_double_vector::zeros(t.size());

         for(int k=0,N=t.size();k<N;++k)
         {
             k_seg[k]=static_cast<int>(t[k]);
             t_seg[k]=t[k]-k_seg[k];
         }
         return value(k_seg,t_seg);
     }
     MC_curve MC_curve::isometric_deformation(const MC_int_vector& constraint_index,const MC_v3d_vector& constraint_position,const MC_curve& original_curve,const double& blending_factor) const
     {
         MC_curve c=*this;
         MC_double_vector L=MC_v3d_vector::norm(original_curve.diff_forward_close());

         int N=size();

         MC_curve c1=*this;
         MC_curve c2=*this;

         std::map<int,MC_v3d> constraint_map;
         for(int k=0,N_constraint=constraint_index.size();k<N_constraint;++k)
             constraint_map.insert(std::make_pair(constraint_index[k],constraint_position[k]));
         std::map<int,MC_v3d>::const_iterator constraint_end=constraint_map.end();

         for(std::map<int,MC_v3d>::const_iterator it=constraint_map.begin();it!=constraint_end;++it)
         {
             c1[it->first]=it->second;
             c2[it->first]=it->second;
         }

         //iteration loop

         for(int k_loop=0;k_loop<3*N;++k_loop)
         {
             //forward
             for(int k=0;k<N;++k)
             {
                 std::map<int,MC_v3d>::const_iterator it_current=constraint_map.find(k);
                 if(it_current==constraint_end)
                 {
                     MC_v3d dir0=(original_curve[k]-original_curve[(k-1+N)%N]).normalized();
                     MC_v3d dir=(c1[k]-c1[(k-1+N)%N]).normalized();
                     dir=(dir+blending_factor*dir0).normalized();
                     c1[k] = c1[(k-1+N)%N] + L[k]*dir;
                 }
             }
             //backward
             for(int k=N-1;k>=0;--k)
             {
                 std::map<int,MC_v3d>::const_iterator it_current=constraint_map.find(k);
                 if(it_current==constraint_end)
                 {
                     MC_v3d dir0=(original_curve[k]-original_curve[(k+1)%N]).normalized();
                     MC_v3d dir=(c2[k]-c2[(k+1)%N]).normalized();
                     dir=(dir+blending_factor*dir0).normalized();
                     c2[k] = c2[(k+1)%N] + L[k]*dir;
                 }
             }
             //average
             c1=0.5*(c1+c2);
             c2=c1;
         }

         return c1;

     }

 }