MC_v3d_vector.cpp

Go to the documentation of this file.

#include <MC_v3d_vector.hpp>
#include <MC_matrix.hpp>

#include <MC_double_vector.hpp>
#include <MC_int_vector.hpp>
#include <MC_connectivity_index.hpp>
#include <MC_mesh_index_vector.hpp>
#include <MC_polygon.hpp>
#include <MC_segment.hpp>


namespace mesh_conv
{

    MC_v3d_vector::MC_v3d_vector(){}
    MC_v3d_vector::MC_v3d_vector(const MC_v3d& v0){v.push_back(v0);}
    MC_v3d_vector::MC_v3d_vector(const MC_v3d& v0,const MC_v3d& v1){v.push_back(v0);v.push_back(v1);}
    MC_v3d_vector::MC_v3d_vector(const MC_v3d& v0,const MC_v3d& v1,const MC_v3d& v2){v.push_back(v0);v.push_back(v1);v.push_back(v2);}
    MC_v3d_vector::MC_v3d_vector(const MC_v3d& v0,const MC_v3d& v1,const MC_v3d& v2,const MC_v3d& v3){v.push_back(v0);v.push_back(v1);v.push_back(v2);v.push_back(v3);}
    MC_v3d_vector::MC_v3d_vector(const MC_v3d& v0,const MC_v3d& v1,const MC_v3d& v2,const MC_v3d& v3,const MC_v3d& v4){v.push_back(v0);v.push_back(v1);v.push_back(v2);v.push_back(v3);v.push_back(v4);}
    MC_v3d_vector::MC_v3d_vector(const int& new_size){v.resize(new_size);}


    int mesh_conv::MC_v3d_vector::size() const {return v.size();}
    void MC_v3d_vector::assert_bounds(const int& u) const
    {
        if(u<0 || u>=int(v.size()))
        {std::cout<<"Assert mesh_conv::MC_v3d_vector::assert_bounds("<<u<<") failed, because size="<<v.size()<<std::endl;exit(-1);}
    }

    const MC_v3d& MC_v3d_vector::operator()(const int& k_index) const {assert_bounds(k_index); return v[k_index];}
    MC_v3d& MC_v3d_vector::operator()(const int& k_index) {assert_bounds(k_index); return v[k_index];}
    const MC_v3d& MC_v3d_vector::operator[](const int& k_index) const {assert_bounds(k_index); return v[k_index];}
    MC_v3d& MC_v3d_vector::operator[](const int& k_index) {assert_bounds(k_index); return v[k_index];}

    MC_v3d_vector::MC_v3d_vector(const MC_v3d_vector& vec)
    {
        int N=vec.v.size();
        for(int k=0;k<N;k++)
            v.push_back(vec.v[k]);
    }
    MC_v3d_vector::MC_v3d_vector(const std::vector <MC_v3d>& vec)
    {
        int N=vec.size();
        for(int k=0;k<N;k++)
            v.push_back(vec[k]);
    }
    MC_v3d_vector::MC_v3d_vector(const std::set <MC_v3d,MC_v3d_less> set)
    {
        std::set <MC_v3d,mesh_conv::MC_v3d_less> :: const_iterator it;
        std::set <MC_v3d,mesh_conv::MC_v3d_less> :: const_iterator it_end=set.end();

        for(it=set.begin();it!=it_end;++it)
            v.push_back(*it);
    }
    MC_v3d_vector::MC_v3d_vector(const std::map <MC_v3d,int,MC_v3d_less> map)
    {
        std::map <MC_v3d,int,mesh_conv::MC_v3d_less> :: const_iterator it;
        std::map <MC_v3d,int,mesh_conv::MC_v3d_less> :: const_iterator it_end=map.end();

        for(it=map.begin();it!=map.end();++it)
            set(it->second,it->first);
    }
    MC_v3d_vector& MC_v3d_vector::add(const MC_v3d& x)
    {v.push_back(x);return *this;}
    MC_v3d_vector& MC_v3d_vector::add(const MC_v3d_vector& v_x)
    {
        int N=v_x.size();
        for(int k=0;k<N;k++)
            v.push_back(v_x[k]);
        return *this;
    }
    MC_v3d_vector& MC_v3d_vector::set(const int& k_index,const MC_v3d& x)
    {
        if(k_index<0)
        {std::cout<<"Error in MC_v3d_vector::set("<<k_index<<","<<x<<"), index must be >0"<<std::endl;exit(-1);}
        if(k_index>=size())
            v.resize(k_index+1);
        v[k_index]=x;
        return *this;
    }
    MC_v3d_vector& MC_v3d_vector::set(const MC_int_vector& k_index,const MC_v3d_vector& x)
    {
        MC_v3d_vector v_x=x;
        if(x.size()!=1 && x.size()!=k_index.size())
        {std::cout<<"Error in MC_v3d_vector::set(MC_int_vector,MC_v3d_vector), size are not correct: "<<k_index.size()<<","<<x.size()<<std::endl;exit(-1);}
        if(x.size()==1)
            v_x=MC_v3d_vector::zeros(k_index.size())+x;

        int N=v_x.size();
        for(int k=0;k<N;k++)
            set(k_index[k],v_x[k]);
        return *this;
    }
    MC_v3d_vector MC_v3d_vector::zeros(const int& new_size)
    {
        MC_v3d_vector new_vec;new_vec.resize(new_size);
        return new_vec;
    }
    void MC_v3d_vector::clear()
    {v.clear();}
    MC_v3d_vector& MC_v3d_vector::resize(const int& new_size)
    {
        int N=v.size();
        if(N==new_size) //same size
            return *this;

        if(new_size<0)
        {std::cout<<"Error in MC_v3d_vector::resize("<<new_size<<") new size must be >0"<<std::endl;exit(-1);}

        v.resize(new_size);
        return *this;
    }

    MC_v3d_vector MC_v3d_vector::operator-() const
    {
        int N=size();
        MC_v3d_vector new_vec;new_vec.resize(size());
        for(int k=0;k<N;++k)
            new_vec[k]=-v[k];
        return new_vec;
    }

    MC_v3d_vector operator+(const MC_v3d_vector& vec,const MC_v3d& to_add)
    {
        int N=vec.size();
        MC_v3d_vector new_vec;new_vec.resize(N);
        for(int k=0;k<N;++k)
            new_vec[k]=vec[k]+to_add;
        return new_vec;
    }
    MC_v3d_vector operator+(const MC_v3d& to_add,const MC_v3d_vector& vec){return vec+to_add;}
    MC_v3d_vector operator-(const MC_v3d_vector& vec,const MC_v3d& to_sub)
    {
        int N=vec.size();
        MC_v3d_vector new_vec;new_vec.resize(N);
        for(int k=0;k<N;++k)
            new_vec[k]=vec[k]-to_sub;
        return new_vec;
    }
    MC_v3d_vector operator*(const MC_v3d_vector& vec,const double& to_mult)
    {
        int N=vec.size();
        MC_v3d_vector new_vec;new_vec.resize(N);
        for(int k=0;k<N;++k)
            new_vec[k]=vec[k]*to_mult;
        return new_vec;
    }
    MC_v3d_vector operator/(const MC_v3d_vector& vec,const double& to_subdiv)
    {
        double epsilon=0.000001;
        if(fabs(to_subdiv)<epsilon)
        {std::cout<<"Error in MC_v3d_vector("<<vec<<"/"<<to_subdiv<<", divide by zero"<<std::endl;exit(-1);}

        int N=vec.size();
        MC_v3d_vector new_vec;new_vec.resize(N);
        for(int k=0;k<N;++k)
            new_vec[k]=vec[k]/to_subdiv;
        return new_vec;
    }
    MC_v3d_vector operator+(const MC_v3d_vector& vec,const MC_v3d_vector& to_add)
    {
        if(to_add.size()==1)
            return vec+to_add.first();
        else
        {
            if(vec.size()!=to_add.size())
            {std::cout<<"Error in MC_v3d_vector::operator+(MC_v3d_vector,MC_v3d_vector), size are not compatible "<<vec.size()<<";"<<to_add.size()<<std::endl;exit(-1);}


            int N=vec.size();
            MC_v3d_vector new_vec=MC_v3d_vector::zeros(N);
            for(int k=0;k<N;++k)
                new_vec[k]=vec[k]+to_add[k];
            return new_vec;
        }
    }
    MC_v3d_vector operator-(const MC_v3d_vector& vec,const MC_v3d_vector& to_sub)
    {
        if(to_sub.size()==1)
            return vec-to_sub.first();
        else
        {
            if(vec.size()!=to_sub.size())
            {std::cout<<"Error in MC_v3d_vector::operator-(MC_v3d_vector,MC_v3d_vector), size are not compatible "<<vec.size()<<";"<<to_sub.size()<<std::endl;exit(-1);}


            int N=vec.size();
            MC_v3d_vector new_vec=MC_v3d_vector::zeros(N);
            for(int k=0;k<N;++k)
                new_vec[k]=vec[k]-to_sub[k];
            return new_vec;
        }
    }
    MC_v3d_vector& MC_v3d_vector::operator+=(const MC_v3d& to_add)
    {
        int N=size();
        for(int k=0;k<N;++k)
            v[k]+=to_add;
        return *this;
    }
    MC_v3d_vector& MC_v3d_vector::operator-=(const MC_v3d& to_sub)
    {
        int N=size();
        for(int k=0;k<N;++k)
            v[k]-=to_sub;
        return *this;
    }
    MC_v3d_vector& MC_v3d_vector::operator*=(const double& to_mult)
    {
        int N=size();
        for(int k=0;k<N;++k)
            v[k]*=to_mult;
        return *this;
    }
    MC_v3d_vector& MC_v3d_vector::operator/=(const double& to_subdiv)
    {
        double epsilon=0.000001;
        if(fabs(to_subdiv)<epsilon)
        {std::cout<<"Error in MC_v3d_vector/="<<to_subdiv<<", divide by zero"<<std::endl;exit(-1);}

        int N=size();
        for(int k=0;k<N;++k)
            v[k]/=to_subdiv;
        return *this;
    }
    MC_v3d_vector& MC_v3d_vector::operator+=(const MC_v3d_vector& to_add)
    {
        if(to_add.size()==1)
            return *this+=to_add.first();
        else
        {
            if(size()!=to_add.size())
            {std::cout<<"Error in MC_v3d_vector::operator+=(MC_v3d_vector), size are not compatible "<<size()<<";"<<to_add.size()<<std::endl;exit(-1);}

            int N=size();
            for(int k=0;k<N;++k)
                v[k]+=to_add[k];
            return *this;
        }
    }
    MC_v3d_vector& MC_v3d_vector::operator-=(const MC_v3d_vector& to_sub)
    {
        if(to_sub.size()==1)
            return *this-=to_sub.first();
        else
        {
            if(size()!=to_sub.size())
            {std::cout<<"Error in MC_v3d_vector::operator-=(MC_v3d_vector), size are not compatible "<<size()<<";"<<to_sub.size()<<std::endl;exit(-1);}

            int N=size();
            for(int k=0;k<N;++k)
                v[k]-=to_sub[k];
            return *this;
        }
    }
    std::ostream& operator<<(std::ostream& output,const MC_v3d_vector& in)
    {
        int N=in.size();
        for(int k=0;k<N;++k)
            output<<in[k];
        return output;
    }

    const MC_v3d& MC_v3d_vector::first() const
    {
        if(size()>0)
            return v[0];
        else
        {
            std::cout<<"Error, call MC_v3d_vector::first() with size=0"<<std::endl;exit(-1);
            exit(-1);
        }
    }
    MC_v3d& MC_v3d_vector::first()
    {
        if(size()>0)
            return v[0];
        else
        {
            std::cout<<"Error, call MC_v3d_vector::first() with size=0"<<std::endl;exit(-1);
            exit(-1);
        }
    }

    const MC_v3d& MC_v3d_vector::last() const
    {
        if(size()>0)
            return v[size()-1];
        else
        {
            std::cout<<"Error, call MC_v3d_vector::last() with size=0"<<std::endl;exit(-1);
        }
    }
    MC_v3d& MC_v3d_vector::last()
    {
        if(size()>0)
            return v[size()-1];
        else
        {
            std::cout<<"Error, call MC_v3d_vector::last() with size=0"<<std::endl;exit(-1);
        }
    }
    MC_v3d MC_v3d_vector::sum(const MC_v3d_vector& vec)
    {
        MC_v3d x;
        int N=vec.size();
        for(int k=0;k<N;++k)
            x+=vec[k];
        return x;
    }

    MC_v3d_vector MC_v3d_vector::operator()(const MC_int_vector& index) const
    {

        int N=size();
        int N_index=index.size();

        MC_v3d_vector vec=MC_v3d_vector::zeros(N_index);
        for(int k_index=0;k_index<N_index;k_index++)
        {
            int u_index=index[k_index];
            if(u_index<0 || u_index>=N)
            {std::cout<<"Error in MC_v3d_vector(MC_int_vector index), index["<<k_index<<"]="<<u_index<<", while size of vector is "<<N<<std::endl;exit(-1);}

            vec[k_index]=v[u_index];
        }
        return vec;
    }

    MC_v3d_vector MC_v3d_vector::operator[](const MC_int_vector& index) const
    {return (*this)(index);}

    MC_v3d_vector operator<<(const MC_v3d_vector& vec,const MC_v3d& value)
    {return MC_v3d_vector(vec).add(value);}

    MC_v3d_vector operator<<(const MC_v3d_vector& vec,const MC_v3d_vector& value)
    {return MC_v3d_vector(vec).add(value);}


    MC_double_vector MC_v3d_vector::component(const int& k_dim) const
    {
        if(k_dim<0 || k_dim>3)
        {std::cout<<"Error in MC_v3d_vector::component("<<k_dim<<"), k_dim must be in [0,2]"<<std::endl;exit(-1);}
        int N=size();
        MC_double_vector vec=MC_double_vector::zeros(N);
        for(int k=0;k<N;k++)
            vec[k]=v[k][k_dim];
        return vec;
    }

    MC_v3d_vector::MC_v3d_vector(const MC_double_vector& x_vector,const MC_double_vector& y_vector,const MC_double_vector& z_vector)
    {
        // check size
        if(x_vector.size()!=y_vector.size() || x_vector.size()!=z_vector.size())
        {std::cout<<"Error in MC_v3d_vector::MC_v3d_vector(MC_double_vector("<<x_vector.size()<<"),MC_double_vector("<<y_vector.size()<<"),MC_double_vector("<<z_vector.size()<<")), size are not consistent"<<std::endl;exit(-1);}

        int N=x_vector.size();
        v.resize(N);
        for(int k=0;k<N;k++)
        {
            v[k][0]=x_vector[k];
            v[k][1]=y_vector[k];
            v[k][2]=z_vector[k];
        }
    }

    MC_v3d_vector operator*(const MC_matrix& M,const MC_v3d_vector& vec)
    {
        int N=vec.size();
        MC_v3d_vector new_vector(N);
        for(int k=0;k<N;++k)
        {new_vector[k]=M*vec[k];}

        return new_vector;
    }
    MC_v3d_vector& MC_v3d_vector::operator*=(const MC_matrix& M)
    {
        int N=size();
        for(int k=0;k<N;k++)
            v[k]*=M;
        return *this;
    }

    MC_v3d_vector operator*(const MC_double_vector w,const MC_v3d_vector& vec)
    {
        int N=vec.size();
        MC_v3d_vector res(N);

        if(w.size()==1 && N!=1)
        {
            double dw=w.first();
            for(int k=0;k<N;++k)
                res[k]=dw*vec[k];
            return res;
        }
        else if(w.size()==N)
        {
            for(int k=0;k<N;++k)
                res[k]=w[k]*vec[k];
            return res;
        }
        else
        {std::cout<<"Error in MC_v3d_vector::operator*(MC_double_vector("<<w.size()<<"),MC_v3d_vector("<<vec.size()<<")), size are not compatible"<<std::endl;exit(-1);}
    }
    MC_v3d_vector operator*(const MC_v3d_vector& vec,const MC_double_vector& w)
    {return w*vec;}
    MC_v3d_vector& MC_v3d_vector::operator*=(const MC_double_vector& w)
    {
        int N=size();
        if(w.size()==1 && N!=1)
        {
            (*this)*=w.first();
            return *this;
        }
        else if(w.size()==N)
        {
            for(int k=0;k<N;++k)
                v[k]*=w[k];
            return *this;
        }
        else
        {std::cout<<"Error in MC_v3d_vector*=(MC_double_vector), size are not compatible: "<<size()<<"!="<<w.size()<<std::endl;exit(-1);}
    }

    MC_v3d_vector MC_v3d_vector::normalized() const
    {
        int N=size();MC_v3d_vector res(N);
        for(int k=0;k<N;++k)
            res[k]=v[k].normalized();
        return res;
    }

    MC_v3d_vector& MC_v3d_vector::scale(const MC_v3d_vector& s)
    {
        int N=size();
        for(int k=0;k<N;++k)
            v[k].scale(s);
        return *this;
    }

    std::map <MC_v3d,int,MC_v3d_less> MC_v3d_vector::to_map() const
    {
        std::map <MC_v3d,int,MC_v3d_less> map_point;
        int N=size();
        for(int k=0;k<N;k++)
            map_point.insert(std::pair<MC_v3d,int>(v[k],k));
        return map_point;
    }
    std::set <MC_v3d,MC_v3d_less> MC_v3d_vector::to_set() const
    {
        std::set <MC_v3d,MC_v3d_less> set_point;
        int N=size();
        for(int k=0;k<N;k++)
            set_point.insert(v[k]);
        return set_point;
    }
    std::list <MC_v3d> MC_v3d_vector::to_list() const
    {
        std::list <MC_v3d> list_point;
        int N=size();
        for(int k=0;k<N;k++)
            list_point.push_back(v[k]);

        return list_point;
    }

    MC_v3d MC_v3d_vector::barycenter() const
    {
        return MC_v3d_vector::sum(*this)/size();
    }
    MC_v3d_vector operator*(const double& to_mult,const MC_v3d_vector& vec)
    {return vec*to_mult;}

    const MC_v3d* MC_v3d_vector::pointer() const{return &v[0];}
    MC_v3d* MC_v3d_vector::pointer_unprotected() {return &v[0];}

    std::pair <MC_v3d_vector,double> MC_v3d_vector::scaled_to_unit() const
    {
        if(size()==0)
        {std::cout<<"Warning in MC_v3d_vector::scaled_to_unit(), size=0"<<std::endl;return std::pair<MC_v3d_vector,double>(*this,1.0);}
        std::pair<MC_v3d,MC_v3d> bb=bounding_box_elements();
        MC_v3d L=bounding_box_length();
        double Lmax=std::max(std::max(L[0],L[1]),L[2]);
        if(Lmax<0.000001)
        {std::cout<<"Error in MC_v3d_vector::scaled_to_unit(), vector has zero length"<<std::endl;exit(-1);}
        double factor=1/Lmax;
        return std::pair <MC_v3d_vector,double> ((*this)*factor,factor);
    }
    std::pair <MC_v3d,MC_v3d> MC_v3d_vector::bounding_box_elements() const
    {
        MC_v3d x_max=first(),x_min=first();
        int N=size();
        for(int k=0;k<N;++k)
        {
            MC_v3d x=v[k];

            for(int k_dim=0;k_dim<3;++k_dim)
            {
                if(x_max[k_dim]<x[k_dim])
                    x_max[k_dim]=x[k_dim];
                if(x_min[k_dim]>x[k_dim])
                    x_min[k_dim]=x[k_dim];
            }
        }
        return std::pair <MC_v3d,MC_v3d> (x_min,x_max);
    }

    MC_v3d MC_v3d_vector::bounding_box_length() const
    {
        std::pair <MC_v3d,MC_v3d> bb=bounding_box_elements();
        return bb.second-bb.first;
    }




    MC_double_vector MC_v3d_vector::norm(const MC_v3d_vector& vec)
    {
        int N=vec.size();
        MC_double_vector n=MC_double_vector::zeros(N);
        for(int k=0;k<N;++k)
            n[k]=vec[k].norm();
        return n;
    }

    std::ostream& MC_v3d_vector::export_stream(std::ostream& output) const
    {
        int N=size();
        for(int k=0;k<N;++k)
            output<<v[k][0]<<" "<<v[k][1]<<" "<<v[k][2]<<" ";
        return output;
    }
    MC_v3d_vector& MC_v3d_vector::read_stream(std::istream& input)
    {
        while(input.good())
        {
            MC_v3d temp(0,0,0);
            for(int k_dim=0;k_dim<3;++k_dim)
                input>>temp[k_dim];
            if(input.good())
                add(temp);
        }
        return *this;
    }

    std::vector <MC_segment> MC_v3d_vector::bounding_box_segment(const double& length_ratio) const
    {
        std::pair <MC_v3d,MC_v3d> bb=bounding_box_elements();

        double min_length=std::min(std::min(
                std::fabs(bb.second[0]-bb.first[0]),
                std::fabs(bb.second[1]-bb.first[1])),
                                   std::fabs(bb.second[2]-bb.first[2]));

        double actual_length=0.5*min_length*length_ratio;

        std::vector <MC_segment> v_seg;

        double epsilon=0.00001;
        if( std::fabs(length_ratio-1.0)>epsilon )
        {
            v_seg.push_back(bb.first+MC_segment(MC_v3d(0,0,0),MC_v3d(actual_length,0,0)));
            v_seg.push_back(bb.first+MC_segment(MC_v3d(0,0,0),MC_v3d(0,actual_length,0)));
            v_seg.push_back(bb.first+MC_segment(MC_v3d(0,0,0),MC_v3d(0,0,actual_length)));

            v_seg.push_back(bb.first.mask(0,1,1)+bb.second.mask(1,0,0)+MC_segment(MC_v3d(0,0,0),MC_v3d(-actual_length,0,0)));
            v_seg.push_back(bb.first.mask(0,1,1)+bb.second.mask(1,0,0)+MC_segment(MC_v3d(0,0,0),MC_v3d(0,actual_length,0)));
            v_seg.push_back(bb.first.mask(0,1,1)+bb.second.mask(1,0,0)+MC_segment(MC_v3d(0,0,0),MC_v3d(0,0,actual_length)));

            v_seg.push_back(bb.first.mask(1,0,1)+bb.second.mask(0,1,0)+MC_segment(MC_v3d(0,0,0),MC_v3d(actual_length,0,0)));
            v_seg.push_back(bb.first.mask(1,0,1)+bb.second.mask(0,1,0)+MC_segment(MC_v3d(0,0,0),MC_v3d(0,-actual_length,0)));
            v_seg.push_back(bb.first.mask(1,0,1)+bb.second.mask(0,1,0)+MC_segment(MC_v3d(0,0,0),MC_v3d(0,0,actual_length)));

            v_seg.push_back(bb.first.mask(1,1,0)+bb.second.mask(0,0,1)+MC_segment(MC_v3d(0,0,0),MC_v3d(actual_length,0,0)));
            v_seg.push_back(bb.first.mask(1,1,0)+bb.second.mask(0,0,1)+MC_segment(MC_v3d(0,0,0),MC_v3d(0,actual_length,0)));
            v_seg.push_back(bb.first.mask(1,1,0)+bb.second.mask(0,0,1)+MC_segment(MC_v3d(0,0,0),MC_v3d(0,0,-actual_length)));


            v_seg.push_back(bb.first.mask(0,0,1)+bb.second.mask(1,1,0)+MC_segment(MC_v3d(0,0,0),MC_v3d(-actual_length,0,0)));
            v_seg.push_back(bb.first.mask(0,0,1)+bb.second.mask(1,1,0)+MC_segment(MC_v3d(0,0,0),MC_v3d(0,-actual_length,0)));
            v_seg.push_back(bb.first.mask(0,0,1)+bb.second.mask(1,1,0)+MC_segment(MC_v3d(0,0,0),MC_v3d(0,0,actual_length)));

            v_seg.push_back(bb.first.mask(0,1,0)+bb.second.mask(1,0,1)+MC_segment(MC_v3d(0,0,0),MC_v3d(-actual_length,0,0)));
            v_seg.push_back(bb.first.mask(0,1,0)+bb.second.mask(1,0,1)+MC_segment(MC_v3d(0,0,0),MC_v3d(0,actual_length,0)));
            v_seg.push_back(bb.first.mask(0,1,0)+bb.second.mask(1,0,1)+MC_segment(MC_v3d(0,0,0),MC_v3d(0,0,-actual_length)));

            v_seg.push_back(bb.first.mask(1,0,0)+bb.second.mask(0,1,1)+MC_segment(MC_v3d(0,0,0),MC_v3d(actual_length,0,0)));
            v_seg.push_back(bb.first.mask(1,0,0)+bb.second.mask(0,1,1)+MC_segment(MC_v3d(0,0,0),MC_v3d(0,-actual_length,0)));
            v_seg.push_back(bb.first.mask(1,0,0)+bb.second.mask(0,1,1)+MC_segment(MC_v3d(0,0,0),MC_v3d(0,0,-actual_length)));

            v_seg.push_back(bb.second+MC_segment(MC_v3d(0,0,0),MC_v3d(-actual_length,0,0)));
            v_seg.push_back(bb.second+MC_segment(MC_v3d(0,0,0),MC_v3d(0,-actual_length,0)));
            v_seg.push_back(bb.second+MC_segment(MC_v3d(0,0,0),MC_v3d(0,0,-actual_length)));
        }
        else
        {
            v_seg.push_back(MC_segment(bb.first,bb.first.mask(0,1,1)+bb.second.mask(1,0,0)));
            v_seg.push_back(MC_segment(bb.first,bb.first.mask(1,0,1)+bb.second.mask(0,1,0)));
            v_seg.push_back(MC_segment(bb.first,bb.first.mask(1,1,0)+bb.second.mask(0,0,1)));

            v_seg.push_back(MC_segment(bb.second,bb.first.mask(0,0,1)+bb.second.mask(1,1,0)));
            v_seg.push_back(MC_segment(bb.second,bb.first.mask(0,1,0)+bb.second.mask(1,0,1)));
            v_seg.push_back(MC_segment(bb.second,bb.first.mask(1,0,0)+bb.second.mask(0,1,1)));

            v_seg.push_back(MC_segment(bb.first.mask(0,1,1)+bb.second.mask(1,0,0),bb.first.mask(0,0,1)+bb.second.mask(1,1,0)));
            v_seg.push_back(MC_segment(bb.first.mask(0,1,1)+bb.second.mask(1,0,0),bb.first.mask(0,1,0)+bb.second.mask(1,0,1)));

            v_seg.push_back(MC_segment(bb.first.mask(1,0,1)+bb.second.mask(0,1,0),bb.first.mask(0,0,1)+bb.second.mask(1,1,0)));
            v_seg.push_back(MC_segment(bb.first.mask(1,0,1)+bb.second.mask(0,1,0),bb.first.mask(1,0,0)+bb.second.mask(0,1,1)));

            v_seg.push_back(MC_segment(bb.first.mask(1,1,0)+bb.second.mask(0,0,1),bb.first.mask(1,0,0)+bb.second.mask(0,1,1)));
            v_seg.push_back(MC_segment(bb.first.mask(1,1,0)+bb.second.mask(0,0,1),bb.first.mask(0,1,0)+bb.second.mask(1,0,1)));

        }

        return v_seg;
    }

    std::pair<MC_v3d,int> MC_v3d_vector::closest(const MC_v3d& x) const
    {
        double min_dist=0;
        int k_min=-1;
        MC_v3d record_min;
        for(int k=0,N=size();k<N;++k)
        {
            const MC_v3d& y=v[k];
            double current_dist=(y-x).norm();
            if(current_dist<min_dist || k==0)
            {
                min_dist=current_dist;
                k_min=k;
                record_min=y;
            }
        }
        return std::make_pair(record_min,k_min);
    }


    std::pair<MC_v3d_vector,MC_matrix> MC_v3d_vector::scaled_to_unit_and_center() const
    {
        std::pair<MC_v3d_vector,double> scaling=scaled_to_unit();
        std::pair<MC_v3d,MC_v3d> bb=scaling.first.bounding_box_elements();
        scaling.first=scaling.first-bb.first;
        MC_matrix M(4);
        for(int k=0;k<3;++k)
            M(k,k)=scaling.second;
        M.set_translation(-bb.first);
        return std::make_pair(scaling.first,M);
    }

    MC_v3d_vector MC_v3d_vector::inverted_position() const
    {
        MC_int_vector index=MC_int_vector::linspace(size()-1,0,-1);
        return (*this)(index);
    }

    MC_v3d_vector MC_v3d_vector::build_from_concatenated_vector(const MC_double_vector& X)
    {
        if(X.size()%3!=0)
        {std::cout<<"Error in MC_v3d_vector::build_from_concatenated_vector(), size is not correct"<<std::endl;exit(-1);}
        MC_v3d_vector y=MC_v3d_vector::zeros(X.size()/3);
        for(int k=0,N=X.size()/3;k<N;++k)
            for(int k_dim=0;k_dim<3;++k_dim)
                y[k][k_dim]=X[3*k+k_dim];
        return y;
    }

    MC_v3d_vector MC_v3d_vector::removed(const MC_int_vector& index_to_remove) const
    {
        MC_v3d_vector new_vec;
        std::set<int> index_set=index_to_remove.to_set();
        for(int k=0,N=size();k<N;++k)
        {
            if(index_set.find(k)==index_set.end())
                new_vec.add(v[k]);
        }
        return new_vec;
    }
    MC_v3d_vector MC_v3d_vector::build_duplicate(const unsigned int& N,const MC_v3d& value_to_duplicate)
    {
        MC_v3d_vector vec;
        vec.resize(N);
        for(int k=0;k<N;++k)
            vec[k]=value_to_duplicate;
        return vec;
    }

}