00001
00002 #include <MC_mesh_index_vector.hpp>
00003
00004 #include <MC_segment.hpp>
00005 #include <MC_polygon.hpp>
00006 #include <MC_triangle.hpp>
00007 #include <MC_matrix.hpp>
00008 #include <MC_int_vector.hpp>
00009 #include <MC_int_vector_vector.hpp>
00010 #include <MC_int_pair.hpp>
00011 #include <MC_double_vector_vector.hpp>
00012 #include <MC_io_off.hpp>
00013 #include <MC_helper_stl.hpp>
00014 #include <MC_string_helper.hpp>
00015
00016 #include <MC_io_off.hpp>
00017 #include <MC_io_obj.hpp>
00018 #include <MC_int_pair.hpp>
00019
00020 #include <sstream>
00021
00022
00023 #include <map>
00024 #include <set>
00025 #include <vector>
00026
00027
00028
00029
00030
00031 namespace mesh_conv
00032 {
00033
00034
00035 MC_mesh_index_vector::MC_mesh_index_vector(){}
00036
00037 int MC_mesh_index_vector::polygon_number() const {return connectivity_mesh.size();}
00038 int MC_mesh_index_vector::vertex_number() const {return point_set_mesh.size();}
00039
00040
00041 MC_mesh_index_vector& MC_mesh_index_vector::add_vertex(const MC_v3d& v) {point_set_mesh.add(v); return *this;}
00042 MC_mesh_index_vector& MC_mesh_index_vector::add_connectivity_index(const MC_int_vector& c) {connectivity_mesh.add(c); return *this;}
00043
00044
00045 MC_connectivity_index& MC_mesh_index_vector::connectivity() {return connectivity_mesh;}
00046 const MC_connectivity_index& MC_mesh_index_vector::connectivity() const {return connectivity_mesh;}
00047 MC_v3d_vector& MC_mesh_index_vector::point_set() {return point_set_mesh;}
00048 const MC_v3d_vector& MC_mesh_index_vector::point_set() const {return point_set_mesh;}
00049
00050 int MC_mesh_index_vector::polygon_size(const int& k_polygon) const
00051 {
00052
00053 if( k_polygon<0 || k_polygon>=connectivity_mesh.size() )
00054 {std::cout<<"Error in MC_mesh_index_vector::polygon_size("<<k_polygon<<"), with connectivity size="<<connectivity_mesh.size()<<std::endl;exit(-1);}
00055
00056 return connectivity_mesh[k_polygon].size();
00057 }
00058 MC_polygon MC_mesh_index_vector::get_polygon(const int& k_index) const
00059 {
00060
00061 if( k_index<0 || k_index>=connectivity_mesh.size() )
00062 {std::cout<<"Error in MC_mesh_index_vector::polygon_size("<<k_index<<"), with connectivity size="<<connectivity_mesh.size()<<std::endl;exit(-1);}
00063 return point_set_mesh(connectivity_mesh[k_index]);
00064 }
00065
00066 MC_mesh_index_vector MC_mesh_index_vector::build_cube()
00067 {
00068 MC_v3d_vector points = MC_v3d_vector(MC_v3d(0,0,0))<<
00069 MC_v3d(1,0,0)<<
00070 MC_v3d(1,1,0)<<
00071 MC_v3d(0,1,0)<<
00072 MC_v3d(1,0,1)<<
00073 MC_v3d(1,1,1)<<
00074 MC_v3d(0,1,1)<<
00075 MC_v3d(0,0,1);
00076
00077 MC_connectivity_index conn = MC_connectivity_index(MC_int_vector(0,1,2,3))<<
00078 MC_int_vector(1,4,5,2)<<
00079 MC_int_vector(2,5,6,3)<<
00080 MC_int_vector(3,6,7,0)<<
00081 MC_int_vector(0,7,4,1)<<
00082 MC_int_vector(6,5,4,7);
00083
00084
00085
00086 return MC_mesh_index_vector(points,conn);
00087 }
00088
00089 MC_mesh_index_vector::MC_mesh_index_vector(const MC_v3d_vector& _point_set,const MC_connectivity_index& _connectivity)
00090 {point_set_mesh=_point_set;connectivity_mesh=_connectivity;}
00091 MC_mesh_index_vector::MC_mesh_index_vector(const std::pair<MC_v3d_vector,MC_connectivity_index>& _mesh)
00092 {*this=MC_mesh_index_vector(_mesh.first,_mesh.second);}
00093
00094 MC_mesh_index_vector operator+(const MC_mesh_index_vector& vec,const MC_v3d& to_add)
00095 {return MC_mesh_index_vector(vec.point_set()+to_add,vec.connectivity());}
00096 MC_mesh_index_vector operator-(const MC_mesh_index_vector& vec,const MC_v3d& to_sub)
00097 {return MC_mesh_index_vector(vec.point_set()-to_sub,vec.connectivity());}
00098 MC_mesh_index_vector operator+(const MC_mesh_index_vector& vec,const MC_v3d_vector& to_add)
00099 {return MC_mesh_index_vector(vec.point_set()+to_add,vec.connectivity());}
00100 MC_mesh_index_vector operator-(const MC_mesh_index_vector& vec,const MC_v3d_vector& to_sub)
00101 {return MC_mesh_index_vector(vec.point_set()-to_sub,vec.connectivity());}
00102 MC_mesh_index_vector operator*(const MC_mesh_index_vector& vec,const double& to_mult)
00103 {return MC_mesh_index_vector(vec.point_set()*to_mult,vec.connectivity());}
00104 MC_mesh_index_vector operator*(const double& to_mult,const MC_mesh_index_vector& vec)
00105 {return vec*to_mult;}
00106 MC_mesh_index_vector operator/(const MC_mesh_index_vector& vec,const double& to_subdiv)
00107 {return MC_mesh_index_vector(vec.point_set()/to_subdiv,vec.connectivity());}
00108
00109 MC_mesh_index_vector& MC_mesh_index_vector::operator+=(const MC_v3d& to_add)
00110 {point_set_mesh+=to_add;return *this;}
00111 MC_mesh_index_vector& MC_mesh_index_vector::operator-=(const MC_v3d& to_sub)
00112 {point_set_mesh-=to_sub;return *this;}
00113 MC_mesh_index_vector& MC_mesh_index_vector::operator+=(const MC_v3d_vector& to_add)
00114 {point_set_mesh+=to_add;return *this;}
00115 MC_mesh_index_vector& MC_mesh_index_vector::operator-=(const MC_v3d_vector& to_sub)
00116 {point_set_mesh-=to_sub;return *this;}
00117 MC_mesh_index_vector& MC_mesh_index_vector::operator*=(const double& to_mult)
00118 {point_set_mesh*=to_mult;return *this;}
00119 MC_mesh_index_vector& MC_mesh_index_vector::operator/=(const double& to_subdiv)
00120 {point_set_mesh/=to_subdiv;return *this;}
00121
00122 MC_mesh_index_vector operator*(const MC_matrix& M,const MC_mesh_index_vector& mesh)
00123 {return MC_mesh_index_vector(M*mesh.point_set(),mesh.connectivity());}
00124 MC_mesh_index_vector& MC_mesh_index_vector::operator*=(const MC_matrix& M)
00125 {point_set_mesh*=M;return *this;}
00126
00127 MC_v3d_vector MC_mesh_index_vector::normal_vertex() const
00128 {
00129
00130
00131 MC_v3d_vector v_n(point_set().size());
00132 if(static_cast<int>(connectivity().star().size())!=point_set().size())
00133 {std::cout<<"Error in MC_v3d_vector MC_mesh_index_vector::normal_vertex(), one star has size "<<connectivity_mesh.star().size()<<", with "<<point_set_mesh.size()<<" vertices, probably need to compute one star"<<std::endl;exit(-1);}
00134 std::map<int,std::set<int> > star = connectivity().star();
00135
00136 std::map<int,std::set<int> > :: const_iterator it=star.begin();
00137 std::map<int,std::set<int> > :: const_iterator it_end=star.end();
00138 for(;it!=it_end;++it)
00139 {
00140
00141 int current_vertex=it->first;
00142 std::set <int> :: const_iterator it_poly=it->second.begin();
00143 std::set <int> :: const_iterator it_poly_end=it->second.end();
00144 for(;it_poly!=it_poly_end;++it_poly)
00145 {
00146 int current_poly=*it_poly;
00147 MC_v3d n=MC_polygon(point_set()(connectivity()(current_poly))).normal();
00148 v_n[current_vertex]+=n;
00149 }
00150 }
00151
00152 v_n=v_n.normalized();
00153 return v_n;
00154 }
00155
00156
00157 MC_v3d_vector MC_mesh_index_vector::normal_vertex(const MC_v3d_vector& point_set_mesh,const MC_connectivity_index& connectivity_mesh)
00158 {
00159
00160 MC_v3d_vector v_n(point_set_mesh.size());
00161 if(static_cast<int>(connectivity_mesh.star().size())!=point_set_mesh.size())
00162 {std::cout<<"Error in MC_v3d_vector MC_mesh_index_vector::normal_vertex(MC_v3d_vector,MC_connectivity_index), one star has size "<<connectivity_mesh.star().size()<<", with "<<point_set_mesh.size()<<" vertices, probably need to compute one star"<<std::endl;exit(-1);}
00163 std::map<int,std::set<int> > star = connectivity_mesh.star();
00164
00165 std::map<int,std::set<int> > :: const_iterator it=star.begin();
00166 std::map<int,std::set<int> > :: const_iterator it_end=star.end();
00167 for(;it!=it_end;++it)
00168 {
00169
00170 int current_vertex=it->first;
00171 std::set <int> :: const_iterator it_poly=it->second.begin();
00172 std::set <int> :: const_iterator it_poly_end=it->second.end();
00173 for(;it_poly!=it_poly_end;++it_poly)
00174 {
00175 int current_poly=*it_poly;
00176 MC_v3d n=MC_polygon(point_set_mesh(connectivity_mesh(current_poly))).normal();
00177 v_n[current_vertex]+=n;
00178 }
00179 }
00180
00181
00182
00183
00184
00185 v_n=v_n.normalized();
00186
00187 return v_n;
00188 }
00189 MC_double_vector MC_mesh_index_vector::normal_vertex(const MC_double_vector& point_set_mesh,const MC_connectivity_index& connectivity_mesh)
00190 {
00191 MC_double_vector v_n=MC_double_vector::zeros(point_set_mesh.size());
00192 if(static_cast<int>(connectivity_mesh.star().size())!=point_set_mesh.size()/3)
00193 {std::cout<<"Error in MC_v3d_vector MC_mesh_index_vector::normal_vertex(MC_v3d_vector,MC_connectivity_index), one star has size "<<connectivity_mesh.star().size()<<", with "<<point_set_mesh.size()<<" vertices, probably need to compute one star"<<std::endl;exit(-1);}
00194 std::map<int,std::set<int> > star = connectivity_mesh.star();
00195
00196 std::map<int,std::set<int> > :: const_iterator it=star.begin();
00197 std::map<int,std::set<int> > :: const_iterator it_end=star.end();
00198 for(;it!=it_end;++it)
00199 {
00200
00201
00202 int current_vertex=it->first;
00203 std::set <int> :: const_iterator it_poly=it->second.begin();
00204 std::set <int> :: const_iterator it_poly_end=it->second.end();
00205 for(;it_poly!=it_poly_end;++it_poly)
00206 {
00207 int current_poly=*it_poly;
00208 MC_int_vector index=connectivity_mesh(current_poly);int N_index=index.size();
00209 MC_polygon p;
00210 for(int k=0;k<N_index;++k)
00211 p.add(MC_v3d(point_set_mesh(3*index(k)+0),point_set_mesh(3*index(k)+1),point_set_mesh(3*index(k)+2)));
00212 MC_v3d n=p.normal();
00213
00214 for(int k_dim=0;k_dim<3;++k_dim)
00215 v_n[3*current_vertex+k_dim]+=n[k_dim];
00216 }
00217 }
00218
00219
00220 double epsilon=0.00000001;
00221 for(int k=0;k<point_set_mesh.size()/3;++k)
00222 {
00223 MC_v3d temp(v_n(3*k+0),v_n(3*k+1),v_n(3*k+2));
00224 double nn=temp.norm();
00225 if(nn>epsilon)
00226 {v_n(3*k+0)/=nn;v_n(3*k+1)/=nn;v_n(3*k+2)/=nn;}
00227 else
00228 {v_n(3*k+0)=1;v_n(3*k+1)=0;v_n(3*k+2)=0;}
00229 }
00230 return v_n;
00231 }
00232
00233 MC_v3d_vector MC_mesh_index_vector::normal_polygon() const
00234 {
00235 int N_poly=polygon_number();
00236 MC_v3d_vector n=MC_v3d_vector::zeros(N_poly);
00237 for(int k=0;k<N_poly;++k)
00238 n[k]=get_polygon(k).normal();
00239 return n;
00240 }
00241
00242
00243 MC_mesh_index_vector MC_mesh_index_vector::build_icosahedron()
00244 {
00245 double p = 0.5*(1+powf(5,0.5));
00246 double f = 1.0/powf(float(1.0+p*p),0.5f);
00247
00248
00249 MC_v3d_vector icosahedron_vertex=MC_v3d_vector(MC_v3d( p, 1, 0))
00250 <<MC_v3d(-p, 1, 0)
00251 <<MC_v3d( p,-1, 0)
00252 <<MC_v3d(-p,-1, 0)
00253 <<MC_v3d( 1, 0, p)
00254 <<MC_v3d( 1, 0,-p)
00255 <<MC_v3d(-1, 0, p)
00256 <<MC_v3d(-1, 0,-p)
00257 <<MC_v3d( 0, p, 1)
00258 <<MC_v3d( 0,-p, 1)
00259 <<MC_v3d( 0, p,-1)
00260 <<MC_v3d( 0,-p,-1) ;
00261
00262
00263 icosahedron_vertex*=f;
00264
00265
00266
00267 MC_connectivity_index icosahedron_tri=MC_connectivity_index(MC_int_vector(0, 8, 4))
00268 <<MC_int_vector( 0, 5,10)
00269 <<MC_int_vector( 2, 4, 9)
00270 <<MC_int_vector( 2,11, 5)
00271 <<MC_int_vector( 1, 6, 8)
00272 <<MC_int_vector( 1,10, 7)
00273 <<MC_int_vector( 3, 9, 6)
00274 <<MC_int_vector( 3, 7,11)
00275 <<MC_int_vector( 0,10, 8)
00276 <<MC_int_vector( 1, 8,10)
00277 <<MC_int_vector( 2, 9,11)
00278 <<MC_int_vector( 3,11, 9)
00279 <<MC_int_vector( 4, 2, 0)
00280 <<MC_int_vector( 5, 0, 2)
00281 <<MC_int_vector( 6, 1, 3)
00282 <<MC_int_vector( 7, 3, 1)
00283 <<MC_int_vector( 8, 6, 4)
00284 <<MC_int_vector( 9, 4, 6)
00285 <<MC_int_vector(10, 5, 7)
00286 <<MC_int_vector(11, 7, 5) ;
00287
00288
00289 return MC_mesh_index_vector(icosahedron_vertex,icosahedron_tri);
00290 }
00291
00292 MC_mesh_index_vector MC_mesh_index_vector::subdivide_mid_edge_unchanged_boundary(const std::set<MC_int_pair,MC_int_pair_less>& allowed_edge) const
00293 {
00294 std::set<MC_int_pair,MC_int_pair_less>::const_iterator allowed_edge_end=allowed_edge.end();
00295
00296 MC_mesh_index_vector new_mesh;
00297 new_mesh.point_set()=point_set();
00298
00299
00300 std::set <int> s=connectivity_mesh.boundary_polygon();
00301 std::set <MC_int_pair,MC_int_pair_less> bnd=connectivity_mesh.boundary();
00302 std::set <int> bnd_poly=connectivity_mesh.boundary_polygon();
00303
00304 std::map <MC_int_pair,int,MC_int_pair_less> newly_added_vertices;
00305
00306 MC_v3d_vector vertex=point_set_mesh;
00307 int N_polygon=polygon_number();
00308 for(int k=0;k<N_polygon;++k)
00309 {
00310 MC_int_vector poly=connectivity_mesh(k);
00311 std::pair<std::vector <MC_polygon>,std::pair<MC_mesh_index_vector,std::pair<MC_int_vector,MC_int_vector> > > current_poly;
00312 if(bnd_poly.find(k)==bnd_poly.end())
00313 current_poly=get_polygon(k).subdivide_mid_edge();
00314 else
00315 {
00316 MC_int_vector constraint_edge;
00317 int poly_size=poly.size();
00318 for(int k_poly=0;k_poly<poly_size;++k_poly)
00319 {
00320 MC_int_pair current_edge(poly[k_poly],poly[(k_poly+1)%poly_size]);
00321 if(bnd.find(current_edge)!=bnd.end())
00322 {
00323 std::set<MC_int_pair,MC_int_pair_less>::const_iterator temp_it=allowed_edge.find(current_edge);
00324 if(temp_it==allowed_edge_end)
00325 constraint_edge.add(k_poly);
00326 }
00327 }
00328 current_poly=get_polygon(k).subdivide_mid_edge(constraint_edge);
00329 }
00330
00331
00332
00333
00334 MC_mesh_index_vector local_mesh=current_poly.second.first;
00335 MC_int_vector extra_vertices=current_poly.second.second.first;
00336 MC_int_vector edge_of_extra_vertices=current_poly.second.second.second;
00337
00338
00339
00340
00341 std::map <int,int> extra_local_to_global;
00342 for(int k_edge=0;k_edge<edge_of_extra_vertices.size();++k_edge)
00343 {
00344 MC_int_pair index_edge_local(edge_of_extra_vertices[k_edge],(edge_of_extra_vertices[k_edge]+1)%poly.size());
00345 MC_int_pair index_edge_global(poly(index_edge_local[0]),poly(index_edge_local[1]));
00346
00347 std::pair <std::map<MC_int_pair,int,MC_int_pair_less>::iterator,bool> insert_it=newly_added_vertices.insert(std::make_pair(index_edge_global,new_mesh.point_set().size()));
00348 if(insert_it.second==true)
00349 {
00350 extra_local_to_global[extra_vertices[k_edge]]=new_mesh.point_set().size();
00351 new_mesh.point_set().add(local_mesh.point_set()[extra_vertices[k_edge]]);
00352 }
00353 else
00354 extra_local_to_global[extra_vertices[k_edge]]=insert_it.first->second;
00355 }
00356
00357 for(int k_new_poly=0;k_new_poly<local_mesh.polygon_number();++k_new_poly)
00358 {
00359 MC_int_vector local_poly=local_mesh.connectivity()(k_new_poly);
00360 MC_int_vector global_poly;
00361 for(int k_new_poly_vertex=0;k_new_poly_vertex<local_poly.size();++k_new_poly_vertex)
00362 {
00363 int u=local_poly[k_new_poly_vertex];
00364 if(u<poly.size())
00365 global_poly.add(poly[u]);
00366 else
00367 global_poly.add(extra_local_to_global[u]);
00368
00369 }
00370 new_mesh.connectivity().add(global_poly);
00371 }
00372 }
00373 return new_mesh;
00374 }
00375
00376
00377
00378 MC_mesh_index_vector MC_mesh_index_vector::subdivide_barycenter_mid_edge_unchanged_boundary(const std::set<MC_int_pair,MC_int_pair_less>& allowed_edge) const
00379 {
00380 MC_mesh_index_vector new_mesh;
00381 new_mesh.point_set()=point_set();
00382
00383
00384 std::set <int> s=connectivity_mesh.boundary_polygon();
00385 std::set <MC_int_pair,MC_int_pair_less> bnd=connectivity_mesh.boundary();
00386 std::set <int> bnd_poly=connectivity_mesh.boundary_polygon();
00387
00388 std::map <MC_int_pair,int,MC_int_pair_less> newly_added_vertices;
00389
00390 MC_v3d_vector vertex=point_set_mesh;
00391 int N_polygon=polygon_number();
00392 for(int k=0;k<N_polygon;++k)
00393 {
00394 MC_int_vector poly=connectivity_mesh(k);
00395 std::pair<std::vector <MC_polygon>,std::pair<MC_mesh_index_vector,std::pair<MC_int_vector,MC_int_vector> > > current_poly;
00396 if(bnd_poly.find(k)==bnd_poly.end())
00397 current_poly=get_polygon(k).subdivide_barycenter_mid_edge();
00398 else
00399 {
00400 MC_int_vector constraint_edge;
00401 int poly_size=poly.size();
00402 for(int k_poly=0;k_poly<poly_size;++k_poly)
00403 {
00404 MC_int_pair current_edge(poly[k_poly],poly[(k_poly+1)%poly_size]);
00405 if(bnd.find(current_edge)!=bnd.end())
00406 if(allowed_edge.find(current_edge)==allowed_edge.end())
00407 constraint_edge.add(k_poly);
00408 }
00409 current_poly=get_polygon(k).subdivide_barycenter_mid_edge(constraint_edge);
00410 }
00411
00412
00413
00414
00415 MC_mesh_index_vector local_mesh=current_poly.second.first;
00416 MC_int_vector extra_vertices=current_poly.second.second.first;
00417 MC_int_vector edge_of_extra_vertices=current_poly.second.second.second;
00418
00419
00420
00421
00422
00423 std::map <int,int> extra_local_to_global;
00424
00425
00426 extra_local_to_global[poly.size()]=new_mesh.point_set().size();
00427 new_mesh.point_set().add(local_mesh.point_set()[poly.size()]);
00428
00429 for(int k_edge=0;k_edge<edge_of_extra_vertices.size();++k_edge)
00430 {
00431 MC_int_pair index_edge_local(edge_of_extra_vertices[k_edge],(edge_of_extra_vertices[k_edge]+1)%poly.size());
00432 MC_int_pair index_edge_global(poly(index_edge_local[0]),poly(index_edge_local[1]));
00433
00434
00435 std::pair <std::map<MC_int_pair,int,MC_int_pair_less>::iterator,bool> insert_it=newly_added_vertices.insert(std::make_pair(index_edge_global,new_mesh.point_set().size()));
00436 if(insert_it.second==true)
00437 {
00438 extra_local_to_global[extra_vertices[k_edge]]=new_mesh.point_set().size();
00439 new_mesh.point_set().add(local_mesh.point_set()[extra_vertices[k_edge]]);
00440 }
00441 else
00442 extra_local_to_global[extra_vertices[k_edge]]=insert_it.first->second;
00443 }
00444
00445
00446
00447 for(int k_new_poly=0;k_new_poly<local_mesh.polygon_number();++k_new_poly)
00448 {
00449 MC_int_vector local_poly=local_mesh.connectivity()(k_new_poly);
00450 MC_int_vector global_poly;
00451 for(int k_new_poly_vertex=0;k_new_poly_vertex<local_poly.size();++k_new_poly_vertex)
00452 {
00453 int u=local_poly[k_new_poly_vertex];
00454 if(u<poly.size())
00455 global_poly.add(poly[u]);
00456 else
00457 global_poly.add(extra_local_to_global[u]);
00458
00459 }
00460 new_mesh.connectivity().add(global_poly);
00461 }
00462 }
00463 return new_mesh;
00464 }
00465
00466
00467 MC_mesh_index_vector MC_mesh_index_vector::subdivide_mixed_mid_edge_unchanged_boundary(const std::set<MC_int_pair,MC_int_pair_less>& allowed_edge) const
00468 {
00469 MC_mesh_index_vector new_mesh;
00470 new_mesh.point_set()=point_set();
00471
00472
00473 std::set <int> s=connectivity_mesh.boundary_polygon();
00474 std::set <MC_int_pair,MC_int_pair_less> bnd=connectivity_mesh.boundary();
00475 std::set <int> bnd_poly=connectivity_mesh.boundary_polygon();
00476
00477 std::map <MC_int_pair,int,MC_int_pair_less> newly_added_vertices;
00478
00479 MC_v3d_vector vertex=point_set_mesh;
00480 int N_polygon=polygon_number();
00481 for(int k=0;k<N_polygon;++k)
00482 {
00483 MC_int_vector poly=connectivity_mesh(k);
00484 std::pair<std::vector <MC_polygon>,std::pair<MC_mesh_index_vector,std::pair<MC_int_vector,MC_int_vector> > > current_poly;
00485 if(bnd_poly.find(k)==bnd_poly.end())
00486 current_poly=get_polygon(k).subdivide_mixed_mid_edge();
00487 else
00488 {
00489 MC_int_vector constraint_edge;
00490 int poly_size=poly.size();
00491 for(int k_poly=0;k_poly<poly_size;++k_poly)
00492 {
00493 MC_int_pair current_edge(poly[k_poly],poly[(k_poly+1)%poly_size]);
00494 if(bnd.find(current_edge)!=bnd.end())
00495 if(allowed_edge.find(current_edge)==allowed_edge.end())
00496 constraint_edge.add(k_poly);
00497 }
00498 current_poly=get_polygon(k).subdivide_mixed_mid_edge(constraint_edge);
00499 }
00500
00501
00502
00503
00504 MC_mesh_index_vector local_mesh=current_poly.second.first;
00505 MC_int_vector extra_vertices=current_poly.second.second.first;
00506 MC_int_vector edge_of_extra_vertices=current_poly.second.second.second;
00507
00508
00509
00510
00511
00512 std::map <int,int> extra_local_to_global;
00513
00514
00515 if(poly.size()!=3)
00516 {
00517 extra_local_to_global[poly.size()]=new_mesh.point_set().size();
00518 new_mesh.point_set().add(local_mesh.point_set()[poly.size()]);
00519 }
00520
00521 for(int k_edge=0;k_edge<edge_of_extra_vertices.size();++k_edge)
00522 {
00523 MC_int_pair index_edge_local(edge_of_extra_vertices[k_edge],(edge_of_extra_vertices[k_edge]+1)%poly.size());
00524 MC_int_pair index_edge_global(poly(index_edge_local[0]),poly(index_edge_local[1]));
00525
00526
00527 std::pair <std::map<MC_int_pair,int,MC_int_pair_less>::iterator,bool> insert_it=newly_added_vertices.insert(std::make_pair(index_edge_global,new_mesh.point_set().size()));
00528 if(insert_it.second==true)
00529 {
00530 extra_local_to_global[extra_vertices[k_edge]]=new_mesh.point_set().size();
00531 new_mesh.point_set().add(local_mesh.point_set()[extra_vertices[k_edge]]);
00532 }
00533 else
00534 extra_local_to_global[extra_vertices[k_edge]]=insert_it.first->second;
00535 }
00536
00537
00538
00539 for(int k_new_poly=0;k_new_poly<local_mesh.polygon_number();++k_new_poly)
00540 {
00541 MC_int_vector local_poly=local_mesh.connectivity()(k_new_poly);
00542 MC_int_vector global_poly;
00543 for(int k_new_poly_vertex=0;k_new_poly_vertex<local_poly.size();++k_new_poly_vertex)
00544 {
00545 int u=local_poly[k_new_poly_vertex];
00546 if(u<poly.size())
00547 global_poly.add(poly[u]);
00548 else
00549 global_poly.add(extra_local_to_global[u]);
00550
00551 }
00552 new_mesh.connectivity().add(global_poly);
00553 }
00554 }
00555 return new_mesh;
00556 }
00557
00558
00559
00560
00561 MC_mesh_index_vector MC_mesh_index_vector::subdivide_mid_edge() const
00562 {
00563
00564
00565
00566 MC_mesh_index_vector mesh2;
00567
00568
00569
00570
00571 std::vector <std::map <int,int> > index_mid_point(vertex_number());
00572
00573
00574
00575 mesh2.point_set()=point_set_mesh;
00576
00577
00578 int current_vertex=vertex_number();
00579 int N_polygon=polygon_number();
00580 for(int k_polygon=0;k_polygon<N_polygon;k_polygon++)
00581 {
00582 MC_int_vector index_poly=connectivity_mesh[k_polygon];
00583 int size_polygon=index_poly.size();
00584 for(int k_edge=0;k_edge<size_polygon;++k_edge)
00585 {
00586 int index_v0=index_poly[k_edge];
00587 int index_v1=index_poly[(k_edge+1)%size_polygon];
00588
00589 if( index_mid_point[index_v0].find(index_v1)==index_mid_point[index_v0].end() )
00590 {
00591
00592 index_mid_point[index_v0].insert(std::pair <int,int> (index_v1,current_vertex) );
00593 index_mid_point[index_v1].insert(std::pair <int,int> (index_v0,current_vertex) );
00594
00595
00596 mesh2.point_set().add(0.5*(point_set_mesh(index_v0)+point_set_mesh(index_v1)));
00597
00598 current_vertex++;
00599 }
00600 }
00601 }
00602
00603
00604 for(int k_polygon=0;k_polygon<N_polygon;++k_polygon)
00605 {
00606
00607 MC_int_vector index_poly=connectivity_mesh[k_polygon];
00608 int size_polygon=index_poly.size();
00609
00610 for(int k_edge=0;k_edge<size_polygon;++k_edge)
00611 {
00612 int index_v0 = index_poly(k_edge);
00613 int index_v1 = index_poly((k_edge+1)%size_polygon);
00614 int index_vm1 = index_poly((k_edge+size_polygon-1)%size_polygon);
00615
00616 int index_mid = index_mid_point[index_v0].find(index_v1)->second;
00617 int index_mid2 = index_mid_point[index_v0].find(index_vm1)->second;
00618
00619 mesh2.connectivity().add(MC_int_vector(index_v0,index_mid,index_mid2));
00620 }
00621
00622
00623 MC_int_vector temp_poly;
00624 for(int k_edge=0;k_edge<size_polygon;++k_edge)
00625 {
00626 int index_v0 = index_poly(k_edge);
00627 int index_v1 = index_poly((k_edge+1)%size_polygon);
00628
00629 int index_mid = index_mid_point[index_v0].find(index_v1)->second;
00630 temp_poly.add(index_mid);
00631 }
00632 mesh2.connectivity().add(temp_poly);
00633 }
00634
00635 return mesh2;
00636 }
00637
00638
00639 MC_mesh_index_vector MC_mesh_index_vector::subdivide_barycenter_mid_edge() const
00640 {
00641
00642
00643
00644 MC_mesh_index_vector mesh2;
00645
00646
00647
00648
00649 std::vector <std::map <int,int> > index_mid_point(vertex_number());
00650 std::map <int,int> barycenter_map;
00651
00652
00653 mesh2.point_set()=point_set_mesh;
00654
00655
00656 int current_vertex=vertex_number();
00657 int N_polygon=polygon_number();
00658 for(int k_polygon=0;k_polygon<N_polygon;k_polygon++)
00659 {
00660 MC_int_vector index_poly=connectivity_mesh[k_polygon];
00661 int size_polygon=index_poly.size();
00662 for(int k_edge=0;k_edge<size_polygon;++k_edge)
00663 {
00664 int index_v0=index_poly[k_edge];
00665 int index_v1=index_poly[(k_edge+1)%size_polygon];
00666
00667 if( index_mid_point[index_v0].find(index_v1)==index_mid_point[index_v0].end() )
00668 {
00669
00670 index_mid_point[index_v0].insert(std::pair <int,int> (index_v1,current_vertex) );
00671 index_mid_point[index_v1].insert(std::pair <int,int> (index_v0,current_vertex) );
00672
00673
00674 mesh2.point_set().add(0.5*(point_set_mesh(index_v0)+point_set_mesh(index_v1)));
00675
00676 current_vertex++;
00677 }
00678 }
00679
00680
00681 barycenter_map.insert(std::pair<int,int>(k_polygon,current_vertex++));
00682 mesh2.point_set().add(MC_v3d_vector::sum(mesh2.point_set()(index_poly))/size_polygon);
00683 }
00684
00685
00686 for(int k_polygon=0;k_polygon<N_polygon;++k_polygon)
00687 {
00688
00689 MC_int_vector index_poly=connectivity_mesh[k_polygon];
00690 int size_polygon=index_poly.size();
00691
00692 int index_barycenter=barycenter_map.find(k_polygon)->second;
00693 for(int k_edge=0;k_edge<size_polygon;++k_edge)
00694 {
00695 int index_v0 = index_poly(k_edge);
00696 int index_v1 = index_poly((k_edge+1)%size_polygon);
00697 int index_vm1 = index_poly((k_edge+size_polygon-1)%size_polygon);
00698
00699 int index_mid = index_mid_point[index_v0].find(index_v1)->second;
00700 int index_mid2 = index_mid_point[index_v0].find(index_vm1)->second;
00701
00702
00703 mesh2.connectivity().add(MC_int_vector(index_v0,index_mid,index_barycenter,index_mid2));
00704 }
00705 }
00706
00707 return mesh2;
00708 }
00709
00710 MC_mesh_index_vector MC_mesh_index_vector::build_sphere(const int& N_subdiv)
00711 {
00712 MC_mesh_index_vector new_sphere;
00713 new_sphere = MC_mesh_index_vector::build_icosahedron();
00714
00715 int k=0;
00716 int N_vertices_1=0,N_vertices_2=0;
00717
00718
00719 for(int k_subdiv=0;k_subdiv<N_subdiv;k_subdiv++)
00720 {
00721 N_vertices_1 = new_sphere.vertex_number();
00722 new_sphere=new_sphere.subdivide_mid_edge();
00723 N_vertices_2=new_sphere.vertex_number();
00724 for(k=N_vertices_1;k<N_vertices_2;k++)
00725 new_sphere.point_set()[k]=new_sphere.point_set()[k].normalized();
00726 }
00727 new_sphere.connectivity().update_neighbors();
00728 return new_sphere;
00729 }
00730
00731 MC_mesh_index_vector MC_mesh_index_vector::build_quad_sphere(const int& N_subdiv)
00732 {
00733 MC_mesh_index_vector new_sphere;
00734 new_sphere = (MC_mesh_index_vector::build_cube()-0.5*MC_v3d(1,1,1))*2/sqrt(3);
00735
00736 int k=0;
00737 int N_vertices_1=0,N_vertices_2=0;
00738
00739
00740 for(int k_subdiv=0;k_subdiv<N_subdiv;k_subdiv++)
00741 {
00742 N_vertices_1 = new_sphere.vertex_number();
00743 new_sphere=new_sphere.subdivide_barycenter_mid_edge();
00744 N_vertices_2=new_sphere.vertex_number();
00745 for(k=N_vertices_1;k<N_vertices_2;k++)
00746 new_sphere.point_set()[k]=new_sphere.point_set()[k].normalized();
00747 }
00748 return new_sphere;
00749
00750 }
00751
00752 MC_mesh_index_vector MC_mesh_index_vector::build_parametric_sphere(const int& N_1,const int& N_2)
00753 {
00754 MC_mesh_index_vector new_sphere;
00755 #define PI 3.14159
00756
00757 int k1=0,k2=0;
00758 double u1=0.0,u2=0.0;
00759 for(k1=0;k1<N_1;k1++)
00760 {
00761 u1=double(k1+1)/double(N_1+1);
00762 for(k2=0;k2<N_2;k2++)
00763 {
00764 u2=double(k2)/double(N_2);
00765
00766 new_sphere.add_vertex(MC_v3d(sin(PI*u1)*cos(2*PI*u2),sin(PI*u1)*sin(2*PI*u2),cos(PI*u1)));
00767 if(k1<N_1-1)
00768 new_sphere.connectivity().add(MC_int_vector(k2 +N_2*k1 ,
00769 (k2+1)%N_2 +N_2*k1 ,
00770 (k2+1)%N_2 +N_2*( (k1+1)%N_1 ),
00771 k2 +N_2*( (k1+1)%N_1 )));
00772 }
00773 }
00774
00775
00776 new_sphere.point_set().add(MC_v3d(0,0,1));
00777 new_sphere.point_set().add(MC_v3d(0,0,-1));
00778
00779
00780
00781 for(k2=0;k2<N_2;k2++)
00782 {
00783 new_sphere.connectivity().add(MC_int_vector(
00784 N_1*N_2,
00785 k2,
00786 (k2+1)%N_2
00787 ));
00788 new_sphere.connectivity().add(MC_int_vector(
00789 N_1*N_2+1,
00790 k2+N_2*(N_1-1),
00791 (k2+1)%N_2+N_2*(N_1-1)
00792 ));
00793 }
00794
00795 return new_sphere;
00796 }
00797
00798 MC_mesh_index_vector MC_mesh_index_vector::build_torus(const double& R0,const double& R1,const int& N_1,const int& N_2)
00799 {
00800
00801 MC_mesh_index_vector new_torus;
00802 #define PI 3.14159
00803
00804 int k1=0,k2=0;
00805 double u1=0.0,u2=0.0;
00806 for(k1=0;k1<N_1;k1++)
00807 {
00808 u1=double(k1)/double(N_1);
00809 for(k2=0;k2<N_2;k2++)
00810 {
00811 u2=double(k2)/double(N_2);
00812
00813 new_torus.point_set().add(MC_v3d(
00814 (R0+R1*cos(u1*2*PI))*cos(u2*2*PI),
00815 (R0+R1*cos(u1*2*PI))*sin(u2*2*PI),
00816 R1*sin(u1*2*PI)
00817 ));
00818
00819 new_torus.connectivity().add(MC_int_vector(
00820 k2 +N_2*k1 ,
00821 (k2+1)%N_2 +N_2*k1 ,
00822 (k2+1)%N_2 +N_2*( (k1+1)%N_1 ),
00823 k2 +N_2*( (k1+1)%N_1 )));
00824 }
00825 }
00826
00827 return new_torus;
00828
00829 }
00830
00831
00832
00833 MC_mesh_index_vector MC_mesh_index_vector::build_closed_cylinder(const int& N1,const int& N2,const bool& is_closed)
00834 {
00835 MC_mesh_index_vector cylinder;
00836
00837
00838 double u1=0.0,u2=0.0;
00839 int k1=0,k2=0;
00840 for(k1=0;k1<N1;k1++)
00841 {
00842 u1 = double(k1)/double(N1-1);
00843 for(k2=0;k2<N2;k2++)
00844 {
00845 u2 = double(k2)/double(N2);
00846 cylinder.point_set().add(MC_v3d(u1,cos(2*PI*u2),sin(2*PI*u2)));
00847 }
00848 }
00849
00850 for(k1=0;k1<N1-1;k1++)
00851 for(k2=0;k2<N2;k2++)
00852 cylinder.connectivity().add(MC_int_vector( k2+0 + N2*k1+0,
00853 (k2+1)%N2 + N2*k1+0,
00854 (k2+1)%N2 + N2*(k1+1),
00855 k2+0 + N2*(k1+1)));
00856
00857
00858 if(is_closed==true)
00859 {
00860
00861 cylinder.add_vertex(MC_v3d(0,0,0));
00862 cylinder.add_vertex(MC_v3d(1,0,0));
00863
00864 MC_int_vector closing;
00865 for(k2=0;k2<N2;++k2)
00866 cylinder.add_connectivity_index(MC_int_vector(k2,N2*N1,(k2+1)%N2));
00867 for(k2=0;k2<N2;++k2)
00868 cylinder.add_connectivity_index(MC_int_vector(N2*N1+1,N2*(N1-1)+k2,N2*(N1-1)+(k2+1)%N2));
00869 }
00870
00871 return cylinder;
00872
00873 }
00874
00875 MC_mesh_index_vector MC_mesh_index_vector::build_cone(const int& N_radius,const int& N_face)
00876 {
00877
00878 #define PI 3.14159
00879
00880
00881 if(N_radius==4 && N_face==2)
00882 {
00883 MC_mesh_index_vector cone;
00884 cone.point_set()=MC_v3d_vector(MC_v3d(1,1,0))<<MC_v3d(1,-1,0)<<MC_v3d(-1,-1,0)<<MC_v3d(-1,1,0)<<MC_v3d(0,0,1);
00885 cone.connectivity().add(MC_int_vector(0,3,2,1));
00886 cone.connectivity().add(MC_int_vector(0,1,4));
00887 cone.connectivity().add(MC_int_vector(1,2,4));
00888 cone.connectivity().add(MC_int_vector(2,3,4));
00889 cone.connectivity().add(MC_int_vector(3,0,4));
00890
00891 return cone;
00892 }
00893
00894 MC_mesh_index_vector cone;
00895
00896
00897 cone.add_vertex(MC_v3d(0,0,0));
00898 cone.add_vertex(MC_v3d(0,0,1));
00899
00900
00901 double r=1.0/(N_face-1);
00902 for(int k=0;k<N_radius;k++)
00903 {
00904 double alpha = 2*PI*double(k)/double(N_radius);
00905 MC_v3d x=MC_v3d(r*cos(alpha),r*sin(alpha),0.0);
00906 cone.add_vertex(x);
00907 }
00908
00909 for(int k2=1;k2<N_face;k2++)
00910 {
00911 double r=double(k2)/double(N_face-1);
00912 for(int k=0;k<N_radius;k++)
00913 {
00914 double alpha = 2*PI*double(k)/double(N_radius);
00915 MC_v3d x = MC_v3d(r*cos(alpha),r*sin(alpha),0.0);
00916 cone.add_vertex(x);
00917 }
00918 }
00919
00920 for(int k2=1;k2<N_face-1;k2++)
00921 {
00922 double r=double(N_face-1-k2)/double(N_face-1);
00923 double h=double(k2)/double(N_face-1);
00924 for(int k=0;k<N_radius;k++)
00925 {
00926 double alpha = 2*PI*double(k)/double(N_radius);
00927 MC_v3d x=MC_v3d(r*cos(alpha),r*sin(alpha),h);
00928 cone.add_vertex(x);
00929 }
00930 }
00931
00932
00933
00934 for(int k=0;k<N_radius;k++)
00935 cone.add_connectivity_index(MC_int_vector(0,(k+1)%N_radius+2,k+2));
00936 for(int k2=1;k2<N_face-1;k2++)
00937 {
00938 for(int k=0;k<N_radius;k++)
00939 {
00940 cone.add_connectivity_index(MC_int_vector((k+0) +(k2+0)*N_radius+2,
00941 (k+1)%N_radius+(k2+0)*N_radius+2,
00942 (k+1)%N_radius+(k2+1)*N_radius+2,
00943 (k+0) +(k2+1)*N_radius+2));
00944 }
00945 }
00946
00948 for(int k2=0;k2<N_face-2;k2++)
00949 for(int k=0;k<N_radius;k++)
00950 {
00951 cone.add_connectivity_index(MC_int_vector((k+0) +(k2+0)*N_radius+N_radius*(N_face-1)+2,
00952 (k+1)%N_radius+(k2+0)*N_radius+N_radius*(N_face-1)+2,
00953 (k+1)%N_radius+(k2+1)*N_radius+N_radius*(N_face-1)+2,
00954 (k+0) +(k2+1)*N_radius+N_radius*(N_face-1)+2));
00955 }
00956
00958 for(int k=0;k<N_radius;k++)
00959 {
00960 cone.add_connectivity_index(MC_int_vector( k+N_radius*(N_face-2)+N_radius*(N_face-1)+2,
00961 (k+1)%N_radius+N_radius*(N_face-2)+N_radius*(N_face-1)+2,
00962 1));
00963 }
00964
00965
00966
00967 return cone;
00968 }
00969
00970 MC_mesh_index_vector MC_mesh_index_vector::build_disc(const MC_v3d& center,const MC_v3d& normal,const double& radius,const int& N_radius,const int& N_interior)
00971 {
00972 double epsilon=0.0001;
00973 if(normal.norm()<epsilon)
00974 {std::cout<<"Error in MC_mesh_index_vector::build_disc(...), normal is null"<<std::endl;exit(-1);}
00975
00976 MC_v3d axis_1 = normal.cross(MC_v3d(1,0,0));
00977 if(axis_1.norm()<epsilon)
00978 {
00979 axis_1 = normal.cross(MC_v3d(0,1,0));
00980 if(axis_1.norm()<epsilon)
00981 {std::cout<<"Error in MC_mesh_index_vector::build_disc, something weird"<<std::endl; exit(-1);}
00982 }
00983
00984 axis_1 = axis_1.normalized();
00985 MC_v3d axis_2 = (axis_1.cross(normal)).normalized();
00986
00987 MC_mesh_index_vector disc;
00988
00989
00990 disc.add_vertex(center);
00991 for(int k_radius=0;k_radius<N_radius;k_radius++)
00992 {
00993 double theta = 2*PI*double(k_radius)/double(N_radius);
00994 for(int k_interior=1;k_interior<N_interior;k_interior++)
00995 {
00996 double r=radius*double(k_interior)/double(N_interior-1);
00997
00998 MC_v3d current=r*cos(theta)*axis_1 + r*sin(theta)*axis_2;
00999 disc.add_vertex(center+current);
01000 }
01001 }
01002
01003
01004
01005 for(int k_radius=0;k_radius<N_radius;k_radius++)
01006 {
01007 for(int k_interior=0;k_interior<N_interior-1;k_interior++)
01008 {
01009 if(k_interior==0)
01010 disc.add_connectivity_index(MC_int_vector( 0,k_radius*(N_interior-1)+1,((k_radius+1)%N_radius)*(N_interior-1)+1 ));
01011 else
01012 disc.add_connectivity_index(MC_int_vector((k_interior-1+0)+(k_radius+0)*(N_interior-1)+1,
01013 (k_interior-1+1)+(k_radius+0)*(N_interior-1)+1,
01014 (k_interior-1+1)+((k_radius+1)%N_radius)*(N_interior-1)+1,
01015 (k_interior-1+0)+((k_radius+1)%N_radius)*(N_interior-1)+1 ));
01016 }
01017 }
01018
01019
01020 return disc;
01021 }
01022
01023 MC_mesh_index_vector MC_mesh_index_vector::build_arrow(const MC_segment& dir,const double& radius_cylinder,const int& N_cylinder,const int& N_shape_cylinder,const double& radius_cone,const double& length_cone,const int& N_cone,const int& N_shape_cone)
01024 {
01025
01026
01027 double epsilon=0.0001;
01028
01029 MC_mesh_index_vector arrow;
01030 MC_v3d V_intermediate = dir.unit_vector()*(dir.length()-length_cone);
01031
01032 MC_mesh_index_vector cylinder = MC_mesh_index_vector::build_closed_cylinder(N_shape_cylinder,N_cylinder);
01033 cylinder.point_set().scale(MC_v3d(V_intermediate.norm(),radius_cylinder,radius_cylinder));
01034
01035 MC_matrix R=MC_matrix::identity(3);
01036 if(dir.unit_vector().cross(MC_v3d(1,0,0)).norm()>epsilon)
01037 R=MC_matrix::rotation_axis_to_axis(MC_v3d(1,0,0),dir.unit_vector());
01038 cylinder = R*cylinder;
01039
01040 MC_mesh_index_vector cone = MC_mesh_index_vector::build_cone(N_cone,N_shape_cone);
01041 cone.point_set().scale(MC_v3d(radius_cone,radius_cone,length_cone));
01042
01043 R=MC_matrix::identity(3);
01044 if(dir.unit_vector().cross(MC_v3d(0,0,1)).norm()>epsilon)
01045 R=MC_matrix::rotation_axis_to_axis(MC_v3d(0,0,1),dir.unit_vector()+MC_v3d(epsilon,0,0));
01046 cone=R*cone+V_intermediate;
01047
01048
01049 arrow = cone<<cylinder;
01050 arrow+=dir[0];
01051
01052 return arrow;
01053
01054
01055 }
01056
01057 MC_mesh_index_vector operator<<(const MC_mesh_index_vector& vec0,const MC_mesh_index_vector& vec1)
01058 {
01059 MC_mesh_index_vector new_vec;
01060
01061 new_vec.point_set()=vec0.point_set()<<vec1.point_set();
01062 new_vec.connectivity()=vec0.connectivity()<<(vec1.connectivity()+vec0.vertex_number());
01063
01064 return new_vec;
01065 }
01066
01067
01068 MC_mesh_index_vector& MC_mesh_index_vector::concatenation(const MC_mesh_index_vector& vec1)
01069 {
01070 int N_vertex=point_set().size();
01071 point_set().add(vec1.point_set());
01072 connectivity().add(vec1.connectivity()+N_vertex);
01073 return *this;
01074 }
01075
01076 std::pair <std::vector <MC_curve>,MC_int_vector_vector> MC_mesh_index_vector::boundary_curve() const
01077 {
01078
01079 MC_int_vector_vector full_index(MC_int_vector::connected_component(connectivity_mesh.boundary()));
01080
01081
01082 int N=full_index.size();
01083 std::vector <MC_curve> full_curve(N);
01084 for(int k=0;k<N;++k)
01085 full_curve[k]=point_set_mesh(full_index(k));
01086
01087 return std::pair <std::vector <MC_curve>,MC_int_vector_vector> (full_curve,full_index);
01088
01089 }
01090
01091
01092 std::vector <MC_curve> MC_mesh_index_vector::plane_intersection(const MC_v3d& n,const MC_v3d& x0) const
01093 {
01094 std::vector <MC_curve> curves;
01095 int k_polygon=0,N_polygon=polygon_number();
01096
01097 std::vector <MC_segment> s;int type=-1;
01098 MC_v3d_vector temp;
01099
01100
01101 for(k_polygon=0;k_polygon<N_polygon;k_polygon++)
01102 {
01103
01104 temp = get_polygon(k_polygon).plane_intersection(n,x0,&type);
01105
01106 if(type==1 || type==3)
01107 {
01108
01109 s.push_back(MC_segment(temp[0],temp[1]));
01110 }
01111 }
01112
01113
01114
01115 std::vector<MC_segment> s2=s;
01116 s.clear();
01117 for(int k=0,N=s2.size();k<N;++k)
01118 {
01119 bool is_unique=true;
01120 for(int k2=0;is_unique==true && k2<s.size();++k2)
01121 {
01122 if( (s2[k][0]==s[k2][0] && s2[k][1]==s[k2][1])
01123 ||
01124 (s2[k][1]==s[k2][0] && s2[k][0]==s[k2][1])
01125 )
01126 is_unique=false;
01127 }
01128 if(is_unique==true)
01129 s.push_back(s2[k]);
01130 }
01131
01132
01133
01134
01135 int segment_number=s.size();
01136 if(segment_number==0)
01137 return curves;
01138
01139
01140
01141 MC_int_vector is_segment_added;
01142 is_segment_added.resize(segment_number);
01143
01144 std::list <MC_v3d> current_curve;
01145
01146 current_curve.push_back(s[0][0]);
01147 current_curve.push_back(s[0][1]);
01148 is_segment_added[0]=1;
01149
01150
01151
01152
01153
01154 int k_segment=-1,is_addition=1;
01155 while(is_addition==1)
01156 {
01157
01158 is_addition=0;
01159
01160
01161 for(k_segment=0;k_segment<segment_number;k_segment++){
01162 if(is_segment_added[k_segment]==0)
01163 {
01164 if(s[k_segment][0]==current_curve.back())
01165 {current_curve.push_back(s[k_segment][1]);is_segment_added[k_segment]=1;is_addition=1;}
01166 else if(s[k_segment][1]==current_curve.back())
01167 {current_curve.push_back(s[k_segment][0]);is_segment_added[k_segment]=1;is_addition=1;}
01168 }
01169 }
01170
01171 for(k_segment=0;k_segment<segment_number;k_segment++){
01172 if(is_segment_added[k_segment]==0)
01173 {
01174 if(s[k_segment][0]==current_curve.front())
01175 {current_curve.push_front(s[k_segment][1]);is_segment_added[k_segment]=1;is_addition=1;}
01176 else if(s[k_segment][1]==current_curve.front())
01177 {current_curve.push_front(s[k_segment][0]);is_segment_added[k_segment]=1;is_addition=1;}
01178 }
01179 }
01180
01181
01182 if(is_addition==0)
01183 {
01184
01185 curves.push_back(current_curve);
01186 current_curve.clear();
01187
01188
01189 for(k_segment=0;is_addition==0 && k_segment<segment_number;k_segment++){
01190 if(is_segment_added[k_segment]==0)
01191 {
01192 is_addition=1;
01193 current_curve.push_back(s[k_segment][0]);current_curve.push_back(s[k_segment][1]);
01194 is_segment_added[k_segment]=1;
01195 }
01196 }
01197 }
01198
01199
01200 }
01201 return curves;
01202 }
01203
01204 MC_mesh_index_vector MC_mesh_index_vector::half_space_intersection(const MC_v3d& n,const MC_v3d& x0,int *type) const
01205 {
01206
01207 int k_polygon=0;
01208 int N_polygon=polygon_number();
01209
01210 std::vector <MC_polygon> polygon_soup;
01211
01212 if(type!=0)
01213 *type = 0;
01214 MC_polygon p; int t=0;
01215 MC_polygon q;
01216
01217
01218 for(k_polygon=0;k_polygon<N_polygon;++k_polygon)
01219 {
01220 p = get_polygon(k_polygon);
01221 q = p.half_space_intersection(n,x0,&t);
01222
01223 if(t!=2)
01224 polygon_soup.push_back(q);
01225 if(t!=0)
01226 if(type!=0)
01227 *type=1;
01228
01229
01230 }
01231
01232 MC_mesh_index_vector new_mesh=polygon_soup;
01233
01234 if(new_mesh.vertex_number()==0)
01235 if(type!=0)
01236 *type=2;
01237
01238
01239 return new_mesh;
01240 }
01241
01242 std::pair<MC_mesh_index_vector,std::pair<MC_int_vector_vector,MC_int_vector_vector> > MC_mesh_index_vector::added_polygon_soup(const std::vector<MC_polygon> poly_soup) const
01243 {
01244 MC_mesh_index_vector current_mesh=*this;
01245 std::map <MC_v3d,int,MC_v3d_less> v_3d_map = point_set().to_map();
01246
01247 std::pair <MC_int_vector_vector,MC_int_vector_vector> res;
01248 std::pair <MC_int_vector,MC_int_vector> temp;
01249
01250 int N_poly=poly_soup.size();
01251 for(int k=0;k<N_poly;k++)
01252 {
01253 bool is_added=true;
01254 MC_polygon p=poly_soup[k];
01255 if(p.is_degenerated()==true)
01256 {
01257 std::pair <MC_polygon,std::pair<bool,bool> > up=p.undegenerated();
01258 if(up.second.second==false)
01259 is_added=false;
01260 else
01261 p=up.first;
01262 }
01263 if(is_added==true)
01264 {
01265 temp=current_mesh.add_unique_polygon(p,&v_3d_map);
01266 res.first.add(temp.first); res.second.add(temp.second);
01267 }
01268 }
01269
01270 return std::pair<MC_mesh_index_vector,std::pair<MC_int_vector_vector,MC_int_vector_vector> > (current_mesh,res);
01271 }
01272
01273 std::pair <MC_int_vector,MC_int_vector> MC_mesh_index_vector::add_unique_polygon(const MC_polygon& polygon,std::map <MC_v3d,int,MC_v3d_less>* map_vertices)
01274 {
01275 MC_int_vector index_polygon;
01276 MC_int_vector is_new;
01277
01278
01279 std::map <MC_v3d,int,MC_v3d_less> :: iterator it;
01280
01281 int N_polygon=polygon.size();
01282 for(int k=0;k<N_polygon;++k)
01283 {
01284 std::pair<std::map <MC_v3d,int,MC_v3d_less>::iterator ,bool> it=map_vertices->insert(std::pair<MC_v3d,int>(polygon[k],vertex_number()));
01285
01286 if(it.second==true)
01287 {
01288 point_set().add(polygon[k]);
01289 index_polygon.add(vertex_number()-1);
01290 is_new.add(1);
01291 }
01292 else
01293 {
01294 index_polygon.add(it.first->second);
01295 is_new.add(0);
01296 }
01297 }
01298
01299 connectivity().add(index_polygon);
01300 return std::pair < MC_int_vector , MC_int_vector > (index_polygon,is_new);
01301 }
01302
01303 MC_mesh_index_vector::MC_mesh_index_vector(const std::vector<MC_polygon> poly_soup)
01304 {
01305 *this=this->added_polygon_soup(poly_soup).first;
01306 }
01307
01308 std::pair <MC_v3d_vector,std::pair <MC_int_vector,MC_double_vector> > MC_mesh_index_vector::segment_intersection(const MC_segment& s) const
01309 {
01310 MC_int_vector intersected_polygon;
01311 MC_double_vector segment_coord;
01312 MC_v3d_vector inter;
01313 int N=polygon_number();
01314 MC_polygon p;
01315 int type=0;
01316 for(int k=0;k<N;k++)
01317 {
01318 p=get_polygon(k);
01319 MC_v3d i = p.segment_intersection(s,&type);
01320
01321 if(type==1 || type==2 || type==3)
01322 {
01323 inter.add(i);
01324 intersected_polygon.add(k);
01325
01326 bool is_alignated=false;
01327 segment_coord.add(s.relative_position(i,&is_alignated));
01328 if(is_alignated==false)
01329 {std::cout<<"Something wrong in MC_mesh_index_vector::segment_intersection()"<<std::endl;exit(-1);}
01330 }
01331 }
01332
01333
01334 std::map <double,int> map_order=segment_coord.to_map();
01335 std::map <double,int>::const_iterator it=map_order.begin(),it_end=map_order.end();
01336
01337 std::pair <MC_v3d_vector,std::pair <MC_int_vector,MC_double_vector> > res;
01338 for(;it!=it_end;++it)
01339 {
01340 int index=it->second;
01341 res.first.add(inter[index]);
01342 res.second.first.add(intersected_polygon[index]);
01343 res.second.second.add(segment_coord[index]);
01344 }
01345
01346 return res;
01347 }
01348
01349 std::pair <MC_int_vector,std::pair <MC_v3d_vector,std::pair <MC_int_vector,MC_double_vector> > >
01350 MC_mesh_index_vector::segment_intersection(const std::vector <MC_mesh_index_vector>& v_mesh,const MC_segment& s)
01351 {
01352 std::vector <std::pair <MC_v3d_vector,std::pair <MC_int_vector,MC_double_vector> > > v_inter;
01353 int N=v_mesh.size();
01354
01355
01356 for(int k=0;k<N;++k)
01357 v_inter.push_back(v_mesh[k].segment_intersection(s));
01358
01359
01360 MC_int_vector temp_index;
01361 MC_v3d_vector temp_inter;
01362 MC_int_vector temp_poly;
01363 MC_double_vector temp_relative;
01364 for(int k=0;k<N;++k)
01365 {
01366 if(v_inter[k].first.size()>0)
01367 {
01368 temp_index.add(MC_int_vector::zeros(v_inter[k].first.size())+k);
01369 temp_inter.add(v_inter[k].first);
01370 temp_poly.add(v_inter[k].second.first);
01371 temp_relative.add(v_inter[k].second.second);
01372 }
01373 }
01374
01375
01376
01377 std::pair <MC_double_vector,MC_int_vector> ind=temp_relative.sort();
01378
01379 std::pair <MC_int_vector,std::pair <MC_v3d_vector,std::pair <MC_int_vector,MC_double_vector> > > res;
01380 res.first=temp_index(ind.second);
01381 res.second.first=temp_inter(ind.second);
01382 res.second.second.first=temp_poly(ind.second);
01383 res.second.second.second=ind.first;
01384
01385 return res;
01386 }
01387 MC_mesh_index_vector MC_mesh_index_vector::sweep_surface(const MC_curve& _c,const MC_curve& _pattern,const int& N_subdiv,const bool& is_closed,const MC_v3d_vector& e1)
01388 {
01389
01390 MC_curve c=_c.value_t(MC_double_vector::linspace(0,_c.size()-1, 1.0/static_cast<double>(N_subdiv-1)) ).first;
01391
01392 bool is_given_tangent=(c.size()==e1.size());
01393
01394 MC_curve pattern=_pattern;
01395 if(pattern.size()==0)
01396 pattern=MC_curve::build_circle(0.1,8,MC_v3d(0,0,1));
01397
01398 int N_curve=c.size();
01399 int k_curve=0;
01400
01401 MC_mesh_index_vector sweep;
01402 int N_pattern=pattern.size();
01403
01404
01405 MC_v3d_vector normal=c.diff_forward();
01406
01407
01408
01409 MC_matrix R0 =MC_matrix::identity(3);
01410 MC_matrix R0_old=MC_matrix::identity(3);
01411
01412 MC_matrix R1 =MC_matrix::identity(3);
01413 MC_matrix R1_old=MC_matrix::identity(3);
01414
01415
01416 MC_curve temp_pattern;
01417 int k_pattern=0;
01418
01419
01420
01421 MC_v3d ez(0,0,1),ex(1,0,0),ey(0,1,0);
01422 MC_v3d e0x,e1x,e1y,e0p;
01423 for(k_curve=0;k_curve<N_curve;++k_curve)
01424 {
01425
01426
01427
01428
01429
01430
01431 R0_old=R0;
01432 R0=MC_matrix::rotation_axis_to_axis(ez,normal[k_curve]);
01433 temp_pattern = R0*pattern;
01434
01435
01436
01437 if(is_given_tangent==false)
01438 {
01439 if(k_curve>=1)
01440 {
01441
01442 e0x = R1_old*R0_old*ex;
01443
01444 e1x = R0*ex;
01445 e1y = R0*ey;
01446
01447 e0p = (e0x.dot(e1x))*e1x+(e0x.dot(e1y))*e1y;
01448
01449
01450
01451
01452 R1=MC_matrix::rotation_axis_to_axis(e1x,e0p);
01453 temp_pattern = R1*temp_pattern;
01454
01455 }
01456 }
01457 else
01458 {
01459 e1x = R0*ex;
01460 R1=MC_matrix::rotation_axis_to_axis(e1x,e1[k_curve]);
01461 temp_pattern = R1*temp_pattern;
01462 }
01463
01464
01465
01466 temp_pattern=temp_pattern+c[k_curve];
01467
01468
01469
01470 R0_old=R0;R1_old=R1;
01471
01472
01473
01474
01475
01476
01477
01478
01479
01481 sweep.point_set().add(temp_pattern);
01482
01484 if(k_curve<N_curve-1)
01485 {
01486 for(k_pattern=0;k_pattern<N_pattern-1;k_pattern++)
01487 {
01488 sweep.connectivity().add(MC_int_vector(k_pattern+N_pattern*(k_curve),
01489 k_pattern+N_pattern*(k_curve+1),
01490 k_pattern+1+N_pattern*(k_curve+1),
01491 k_pattern+1+N_pattern*(k_curve)
01492 ));
01493 }
01494
01495 sweep.connectivity().add(MC_int_vector(k_pattern+N_pattern*(k_curve),
01496 k_pattern+N_pattern*(k_curve+1),
01497 0+N_pattern*(k_curve+1),
01498 0+N_pattern*(k_curve)
01499 ));
01500 }
01501
01502
01503
01504 }
01505
01506
01507
01508
01509
01510
01511
01512 if(is_closed==true)
01513 {
01514
01515 int k_subdiv=0;
01516 double s=0.0;
01517 int offset=0;
01518
01519 R0 =MC_matrix::identity(3);
01520 R0_old=MC_matrix::identity(3);
01521 R1 =MC_matrix::identity(3);
01522 R1_old=MC_matrix::identity(3);
01523
01524
01525 MC_v3d barycenter;
01526 for(k_curve=0;k_curve<N_curve;k_curve++)
01527 {
01528
01529
01530 R0_old=R0;
01531 R0=MC_matrix::rotation_axis_to_axis(ez,normal[k_curve]);
01532 temp_pattern = R0*pattern;
01533
01534
01535 if(is_given_tangent==false)
01536 {
01537 if(k_curve>=1)
01538 {
01539 e0x = R1_old*R0_old*ex;
01540
01541 e1x = R0*ex;
01542 e1y = R0*ey;
01543
01544 e0p = (e0x.dot(e1x))*e1x+(e0x.dot(e1y))*e1y;
01545
01546
01547 R1=MC_matrix::rotation_axis_to_axis(e1x,e0p);
01548 temp_pattern = R1*temp_pattern;
01549 }
01550 }
01551 else
01552 {
01553 e1x = R0*ex;
01554 R1=MC_matrix::rotation_axis_to_axis(e1x,e1[k_curve]);
01555 temp_pattern = R1*temp_pattern;
01556 }
01557
01558
01559 R0_old=R0;R1_old=R1;
01560
01561
01562
01563
01564
01565
01566
01567
01569 if(k_curve==0 || k_curve==N_curve-1)
01570 {
01571 if(k_curve==0)
01572 offset=N_curve*N_pattern;
01573 if(k_curve==N_curve-1)
01574 offset=N_curve*N_pattern+(N_subdiv-1)*N_pattern+1;
01575 for(k_subdiv=0;k_subdiv<N_subdiv-1;k_subdiv++)
01576 {
01577 s = double(k_subdiv+1)/double(N_subdiv+1);
01578 sweep.point_set().add((1-s)*temp_pattern+c[k_curve]);
01579
01580
01581 if(k_subdiv<N_subdiv-2)
01582 for(k_pattern=0;k_pattern<N_pattern;k_pattern++)
01583 {
01584 if(k_curve==0)
01585 sweep.connectivity().add(MC_int_vector(offset + k_pattern +N_pattern*(k_subdiv ),
01586 offset + (k_pattern+1)%N_pattern+N_pattern*(k_subdiv ),
01587 offset + (k_pattern+1)%N_pattern+N_pattern*(k_subdiv+1),
01588 offset + k_pattern +N_pattern*(k_subdiv+1)
01589 ));
01590 else
01591 sweep.connectivity().add(MC_int_vector(offset + k_pattern +N_pattern*(k_subdiv ),
01592 offset + k_pattern +N_pattern*(k_subdiv+1),
01593 offset + (k_pattern+1)%N_pattern+N_pattern*(k_subdiv+1),
01594 offset + (k_pattern+1)%N_pattern+N_pattern*(k_subdiv )
01595 ));
01596 }
01597 }
01598
01599 for(k_pattern=0;k_pattern<N_pattern;k_pattern++)
01600 {
01601 if(k_curve==0)
01602 sweep.connectivity().add(MC_int_vector(k_pattern +N_pattern*k_curve,
01603 (k_pattern+1)%N_pattern +N_pattern*k_curve,
01604 offset + (k_pattern+1)%N_pattern ,
01605 offset + k_pattern
01606 ));
01607 else
01608 sweep.connectivity().add(MC_int_vector(k_pattern +N_pattern*k_curve,
01609 offset + k_pattern ,
01610 offset + (k_pattern+1)%N_pattern ,
01611 (k_pattern+1)%N_pattern +N_pattern*k_curve
01612 ));
01613 }
01614
01615
01616 barycenter = (temp_pattern+c[k_curve]).barycenter();
01617
01618
01619 sweep.add_vertex(barycenter);
01620
01621
01622 for(k_pattern=0;k_pattern<N_pattern;k_pattern++)
01623 {
01624 if(k_curve==0)
01625 sweep.connectivity().add(MC_int_vector(offset + k_pattern + N_pattern*(N_subdiv-2),
01626 offset + (k_pattern+1)%N_pattern + N_pattern*(N_subdiv-2),
01627 offset + N_pattern*(N_subdiv-1)
01628 ));
01629 else
01630 sweep.connectivity().add(MC_int_vector(offset + k_pattern + N_pattern*(N_subdiv-2),
01631 offset + N_pattern*(N_subdiv-1),
01632 offset + (k_pattern+1)%N_pattern + N_pattern*(N_subdiv-2)
01633 ));
01634 }
01635 }
01636 }
01637
01638 }
01639
01640
01641
01642 return sweep;
01643 }
01644 double MC_mesh_index_vector::volume() const
01645 {
01646 double vol=0.0;
01647 int N_polygon=polygon_number();
01648 for(int k_polygon=0;k_polygon<N_polygon;++k_polygon)
01649 {
01650 MC_int_vector index_polygon=connectivity()(k_polygon);
01651 int N_vertex=index_polygon.size();
01652 int N_tri=N_vertex-2;
01653 for(int k_tri=0;k_tri<N_tri;++k_tri)
01654 {
01655 MC_v3d x0=point_set()(index_polygon[0]);
01656 MC_v3d x1=point_set()(index_polygon[k_tri+1]);
01657 MC_v3d x2=point_set()(index_polygon[k_tri+2]);
01658
01659 vol+= (x0[2]+x1[2]+x2[2])*
01660 ((x1[0]-x0[0])*(x2[1]-x0[1])-(x1[1]-x0[1])*(x2[0]-x0[0]));
01661 }
01662 }
01663 vol/=6.0;
01664 return vol;
01665 }
01666 MC_v3d_vector MC_mesh_index_vector::volume_gradient() const
01667 {
01668 int N_vertex=vertex_number();
01669 MC_v3d_vector grad(N_vertex);
01670 int N_polygon=polygon_number();
01671
01672 for(int k_polygon=0;k_polygon<N_polygon;++k_polygon)
01673 {
01674 MC_int_vector index_polygon=connectivity()(k_polygon);
01675 int N_vertex=index_polygon.size();
01676 int N_tri=N_vertex-2;
01677 for(int k_tri=0;k_tri<N_tri;++k_tri)
01678 {
01679 int index_0=index_polygon[0];
01680 int index_1=index_polygon[k_tri+1];
01681 int index_2=index_polygon[k_tri+2];
01682
01683 MC_v3d x0=point_set()(index_0);
01684 MC_v3d x1=point_set()(index_1);
01685 MC_v3d x2=point_set()(index_2);
01686
01687 grad[index_0][0] += 0.5* (x0[2]+x1[2]+x2[2]) * (x1[1]-x2[1]) * 1/3.0;
01688 grad[index_1][0] += 0.5* (x0[2]+x1[2]+x2[2]) * (x2[1]-x0[1]) * 1/3.0;
01689 grad[index_2][0] += 0.5* (x0[2]+x1[2]+x2[2]) * (x0[1]-x1[1]) * 1/3.0;
01690
01691 grad[index_0][1] += 0.5* (x0[2]+x1[2]+x2[2]) * (x2[0]-x1[0]) * 1/3.0;
01692 grad[index_1][1] += 0.5* (x0[2]+x1[2]+x2[2]) * (x0[0]-x2[0]) * 1/3.0;
01693 grad[index_2][1] += 0.5* (x0[2]+x1[2]+x2[2]) * (x1[0]-x0[0]) * 1/3.0;
01694
01695 grad[index_0][2] += 0.5*(x1[0]*x2[1]-x1[0]*x0[1]-x0[0]*x2[1]-x1[1]*x2[0]+x1[1]*x0[0]+x0[1]*x2[0]) * 1/3.0;
01696 grad[index_1][2] += 0.5*(x1[0]*x2[1]-x1[0]*x0[1]-x0[0]*x2[1]-x1[1]*x2[0]+x1[1]*x0[0]+x0[1]*x2[0]) * 1/3.0;
01697 grad[index_2][2] += 0.5*(x1[0]*x2[1]-x1[0]*x0[1]-x0[0]*x2[1]-x1[1]*x2[0]+x1[1]*x0[0]+x0[1]*x2[0]) * 1/3.0;
01698 }
01699 }
01700
01701 return grad;
01702 }
01703
01704
01705
01706
01707 MC_mesh_index_vector MC_mesh_index_vector::build_square(const int& N_1,const int& N_2)
01708 {
01709
01710 MC_mesh_index_vector mesh;
01711
01712 for(int k_1=0;k_1<N_1+2;++k_1){
01713 double alpha_1=static_cast<double>(k_1)/static_cast<double>(N_1+1);
01714 for(int k_2=0;k_2<N_2+2;++k_2){
01715 double alpha_2=static_cast<double>(k_2)/static_cast<double>(N_2+1);
01716 mesh.point_set().add(MC_v3d(alpha_1,alpha_2,0.0));
01717 }
01718 }
01719
01720 for(int k_1=0;k_1<N_1+2-1;++k_1){
01721 for(int k_2=0;k_2<N_2+2-1;++k_2){
01722 mesh.connectivity().add(MC_int_vector((k_1+0)*(N_2+2)+(k_2+0),
01723 (k_1+0)*(N_2+2)+(k_2+1),
01724 (k_1+1)*(N_2+2)+(k_2+1),
01725 (k_1+1)*(N_2+2)+(k_2+0)));
01726 }
01727 }
01728
01729 return mesh;
01730 }
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743 MC_double_vector_vector MC_mesh_index_vector::barycentric_coordinates(const MC_mesh_index_vector& mesh,const MC_v3d_vector& vertices,const MC_int_vector& belonging_polygon)
01744 {
01745 if(belonging_polygon.size()!=vertices.size())
01746 {std::cout<<"Error in MC_mesh_index_vector::barycentric_coordinates(), size are not compatible"<<std::endl;exit(-1);}
01747
01748 int N=vertices.size();
01749 int N_poly=mesh.polygon_number();
01750 MC_double_vector_vector bar=MC_double_vector_vector::empty(N);
01751 for(int k=0;k<N;++k)
01752 {
01753 int index_poly=belonging_polygon(k);
01754 if(index_poly<0 || index_poly>N_poly)
01755 {std::cout<<"Error in MC_mesh_index_vector::barycentric_coordinates() at index k="<<k<<"/"<<N<<" not correct with N_polygon="<<N_poly<<std::endl;exit(-1);}
01756 bar[k]=mesh.get_polygon(index_poly).barycentric_coordinates(vertices[k]);
01757 }
01758 return bar;
01759 }
01760
01761
01762
01763 std::pair <std::pair<MC_mesh_index_vector,MC_mesh_index_vector>,std::pair<MC_int_vector,MC_int_vector> > MC_mesh_index_vector::delete_polygon(const std::set<int> index_to_delete) const
01764 {
01765 int N_poly=polygon_number();
01766
01767 std::set <int> :: const_iterator index_to_delete_end=index_to_delete.end();
01768
01769 std::map <int,int> map_interior;
01770 std::map <int,int> map_exterior;
01771
01772 MC_v3d_vector point_set_interior;
01773 MC_v3d_vector point_set_exterior;
01774 MC_connectivity_index connectivity_exterior;
01775 MC_connectivity_index connectivity_interior;
01776
01777 MC_int_vector correspondance_interior;
01778 MC_int_vector correspondance_exterior;
01779
01780 for(int k_poly=0;k_poly<N_poly;++k_poly)
01781 {
01782 MC_int_vector current_poly_index=connectivity()(k_poly);
01783 int N_current=current_poly_index.size();
01784 if(index_to_delete.find(k_poly)==index_to_delete_end)
01785 {
01786 MC_int_vector exterior_polygon_index=MC_int_vector::zeros(N_current);
01787 for(int k_current=0;k_current<N_current;++k_current)
01788 {
01789 std::pair<std::map<int,int>::iterator,bool> inserted_exterior=map_exterior.insert(std::make_pair(current_poly_index[k_current],map_exterior.size()));
01790 if(inserted_exterior.second==true)
01791 {
01792 exterior_polygon_index[k_current]=map_exterior.size()-1;
01793 point_set_exterior.add(point_set()(current_poly_index[k_current]));
01794 correspondance_exterior.add(current_poly_index[k_current]);
01795 }
01796 else
01797 exterior_polygon_index[k_current]=inserted_exterior.first->second;
01798 }
01799 connectivity_exterior.add(exterior_polygon_index);
01800 }
01801 else
01802 {
01803 MC_int_vector interior_polygon_index=MC_int_vector::zeros(N_current);
01804 for(int k_current=0;k_current<N_current;++k_current)
01805 {
01806 std::pair<std::map<int,int>::iterator,bool> inserted_interior=map_interior.insert(std::make_pair(current_poly_index[k_current],map_interior.size()));
01807 if(inserted_interior.second==true)
01808 {
01809 interior_polygon_index[k_current]=map_interior.size()-1;
01810 point_set_interior.add(point_set()(current_poly_index[k_current]));
01811 correspondance_interior.add(current_poly_index[k_current]);
01812 }
01813 else
01814 interior_polygon_index[k_current]=inserted_interior.first->second;
01815 }
01816 connectivity_interior.add(interior_polygon_index);
01817 }
01818 }
01819
01820 return std::make_pair(std::make_pair(MC_mesh_index_vector(point_set_exterior,connectivity_exterior),MC_mesh_index_vector(point_set_interior,connectivity_interior)),std::make_pair(correspondance_exterior,correspondance_interior));
01821 }
01822
01823
01824
01825
01826
01827 std::pair<std::pair<MC_mesh_index_vector,MC_mesh_index_vector>,std::pair<MC_int_vector,MC_int_vector> > MC_mesh_index_vector::delete_vertex(const std::set<int> index_to_delete) const
01828 {
01829
01830 MC_mesh_index_vector mesh_deleted;
01831 std::map<int,int> map_vertex_deleted;
01832 std::map<int,int> :: iterator map_vertex_deleted_end=map_vertex_deleted.end();
01833 MC_int_vector correspondance_deleted;
01834
01835 MC_mesh_index_vector mesh_undeleted;
01836 std::map<int,int> map_vertex_undeleted;
01837 std::map<int,int> :: iterator map_vertex_undeleted_end=map_vertex_undeleted.end();
01838 MC_int_vector correspondance_undeleted;
01839
01840 int N_polygon=polygon_number();
01841 std::set<int> ::const_iterator it_index_to_delete_end=index_to_delete.end();
01842 for(int k_polygon=0;k_polygon<N_polygon;++k_polygon)
01843 {
01844 MC_int_vector current_polygon=connectivity()(k_polygon);
01845 int N_vertex=current_polygon.size();
01846
01847 bool is_deleted_polygon=false;
01848 for(int k_vertex=0;k_vertex<N_vertex;++k_vertex)
01849 if(index_to_delete.find(current_polygon[k_vertex])!=it_index_to_delete_end)
01850 is_deleted_polygon=true;
01851
01852
01853 if(is_deleted_polygon==true)
01854 {
01855 MC_int_vector temp_index;
01856 for(int k_vertex=0;k_vertex<N_vertex;++k_vertex)
01857 {
01858 std::pair<std::map<int,int>::iterator,bool> it=map_vertex_deleted.insert(std::make_pair(current_polygon[k_vertex],map_vertex_deleted.size()));
01859 if(it.second==true)
01860 {
01861 mesh_deleted.point_set().add(point_set()(current_polygon[k_vertex]));
01862 correspondance_deleted.add(it.first->first);
01863 }
01864 temp_index.add(it.first->second);
01865 }
01866 mesh_deleted.connectivity().add(temp_index);
01867 }
01868
01869 else
01870 {
01871 MC_int_vector temp_index;
01872 for(int k_vertex=0;k_vertex<N_vertex;++k_vertex)
01873 {
01874
01875 std::pair<std::map<int,int>::iterator,bool> it=map_vertex_undeleted.insert(std::make_pair(current_polygon[k_vertex],map_vertex_undeleted.size()));
01876 if(it.second==true)
01877 {
01878 mesh_undeleted.point_set().add(point_set()(current_polygon[k_vertex]));
01879 correspondance_undeleted.add(it.first->first);
01880 }
01881 temp_index.add(it.first->second);
01882 }
01883 mesh_undeleted.connectivity().add(temp_index);
01884 }
01885 }
01886
01887
01888 return std::make_pair(std::make_pair(mesh_undeleted,mesh_deleted),std::make_pair(correspondance_undeleted,correspondance_deleted));
01889
01890 }
01891
01892 std::pair<std::pair<MC_mesh_index_vector,MC_mesh_index_vector>,std::pair<MC_int_vector,MC_int_vector> > MC_mesh_index_vector::delete_boundary_polygon() const
01893 {
01894 std::set<int> border_polygon=connectivity_mesh.boundary_polygon();
01895 return delete_polygon(border_polygon);
01896 }
01897
01898
01899
01900 std::vector<MC_polygon> MC_mesh_index_vector::get_polygon() const
01901 {
01902
01903 int N=polygon_number();
01904 std::vector<MC_polygon> v_poly(N);
01905 for(int k=0;k<N;++k)
01906 v_poly[k]=get_polygon(k);
01907 return v_poly;
01908 }
01909
01910 MC_mesh_index_vector MC_mesh_index_vector::subdivide_mid_edge_unchanged_boundary() const
01911 {
01912 std::set<MC_int_pair,MC_int_pair_less> s;
01913 return subdivide_mid_edge_unchanged_boundary(s);
01914 }
01915 MC_mesh_index_vector MC_mesh_index_vector::subdivide_barycenter_mid_edge_unchanged_boundary() const
01916 {
01917 std::set<MC_int_pair,MC_int_pair_less> s;
01918 return subdivide_barycenter_mid_edge_unchanged_boundary(s);
01919 }
01920 MC_mesh_index_vector MC_mesh_index_vector::subdivide_mixed_mid_edge_unchanged_boundary() const
01921 {
01922 std::set<MC_int_pair,MC_int_pair_less> s;
01923 return subdivide_mixed_mid_edge_unchanged_boundary(s);
01924 }
01925
01926 MC_mesh_index_vector MC_mesh_index_vector::subdivide_mixed_mid_edge(const std::set<int>& polygon_to_subdivide) const
01927 {
01928 std::set<MC_int_pair,MC_int_pair_less> s;
01929 return subdivide_mixed_mid_edge(polygon_to_subdivide,s);
01930 }
01931 MC_mesh_index_vector MC_mesh_index_vector::subdivide_mixed_mid_edge(const std::set<int>& polygon_to_subdivide,const std::set<MC_int_pair,MC_int_pair_less>& allowed_edge) const
01932 {
01933
01934 MC_v3d_vector new_mesh_vector=point_set();
01935 MC_connectivity_index new_mesh_connectivity=connectivity().delete_index(MC_int_vector(polygon_to_subdivide)).first;
01936 MC_mesh_index_vector new_mesh(new_mesh_vector,new_mesh_connectivity);
01937
01938 std::pair <std::pair<MC_mesh_index_vector,MC_mesh_index_vector>,std::pair<MC_int_vector,MC_int_vector> > splitted=delete_polygon(polygon_to_subdivide);
01939
01940 std::map <int,int> new_to_original=splitted.second.second.to_map();
01941 std::map <int,int> original_to_new=MC_helper_stl::reverse_map(new_to_original);
01942
01943 std::set<MC_int_pair,MC_int_pair_less> new_allowed_edge;
01944 for(std::set<MC_int_pair,MC_int_pair_less>::const_iterator
01945 it=allowed_edge.begin(),
01946 it_end=allowed_edge.end();
01947 it!=it_end;++it)
01948 new_allowed_edge.insert(MC_int_pair(new_to_original[(*it)[0]],new_to_original[(*it)[1]]));
01949
01950 MC_mesh_index_vector subdivided=splitted.first.second.subdivide_mixed_mid_edge_unchanged_boundary(new_allowed_edge);
01951 int N_vertex_old=splitted.first.second.point_set().size();
01952 int N_vertex_new=subdivided.point_set().size();
01953
01954
01955
01956 for(int k=N_vertex_old;k<N_vertex_new;++k)
01957 new_mesh.point_set().add(subdivided.point_set()(k));
01958 int N_subdivided_poly=subdivided.polygon_number();
01959 for(int k=0;k<N_subdivided_poly;++k)
01960 {
01961 MC_int_vector new_poly=subdivided.connectivity()(k);
01962 int N_new_poly=new_poly.size();
01963 MC_int_vector new_poly_in_full_coordinate=MC_int_vector::zeros(N_new_poly);
01964 for(int k_index=0;k_index<N_new_poly;++k_index)
01965 {
01966 int u=new_poly[k_index];
01967 if(u<N_vertex_old)
01968 new_poly_in_full_coordinate[k_index]=splitted.second.second(u);
01969 else
01970 new_poly_in_full_coordinate[k_index]=(u-N_vertex_old)+point_set().size();
01971 }
01972 new_mesh.connectivity().add(new_poly_in_full_coordinate);
01973 }
01974
01975 new_mesh.connectivity().build_ring().build_star();
01976
01977 return new_mesh;
01978 }
01979 std::pair <MC_int_vector,MC_int_vector> MC_mesh_index_vector::polygons_inside_sphere(const MC_v3d& center,const double& radius) const
01980 {
01981 MC_int_vector index_inside,index_outside;
01982 for(int k_polygon=0,N_polygon=polygon_number();k_polygon<N_polygon;++k_polygon)
01983 {
01984 bool is_inside=false;
01985 MC_int_vector poly=connectivity()(k_polygon);
01986 for(int k=0,N=poly.size();is_inside==false&&k<N;++k)
01987 {
01988 if( (point_set()(poly[k])-center).norm()<=radius)
01989 {
01990 is_inside=true;
01991 index_inside.add(k_polygon);
01992
01993 }
01994 }
01995 if(is_inside==false)
01996 index_outside.add(k_polygon);
01997 }
01998 return std::make_pair(index_inside,index_outside);
01999 }
02000
02001 MC_mesh_index_vector MC_mesh_index_vector::laplacian_smoothing(const double& lambda,const int& steps,const bool& is_boundary_preserving) const
02002 {
02003
02004 if(steps==0)
02005 return *this;
02006
02007 MC_mesh_index_vector new_mesh=*this;
02008
02009
02010 std::map<int,std::set<int> > r=connectivity().ring();
02011
02012 std::set<int> bnd=connectivity().boundary_vertex();
02013
02014
02015 for(int k_step=0;k_step<steps;++k_step)
02016 {
02017 MC_v3d_vector temp=MC_v3d_vector::zeros(vertex_number());
02018 for(int k_vertex=0,N_vertex=vertex_number();k_vertex<N_vertex;++k_vertex)
02019 {
02020 MC_v3d x=new_mesh.point_set()(k_vertex);
02021
02022 if(is_boundary_preserving==false ||
02023 (is_boundary_preserving==true && bnd.find(k_vertex)==bnd.end()) )
02024 {
02025
02026 MC_v3d bar;
02027 std::set<int> current_ring=r[k_vertex];
02028 for(std::set<int>::const_iterator it=current_ring.begin(),it_end=current_ring.end();it!=it_end;++it)
02029 bar += new_mesh.point_set()(*it);
02030
02031 if(current_ring.size()>0)
02032 bar/=current_ring.size();
02033
02034 temp[k_vertex]=(1-lambda)*x+lambda*bar;
02035 }
02036 else
02037 temp[k_vertex]=x;
02038 }
02039 new_mesh.point_set()=temp;
02040 }
02041
02042 return new_mesh;
02043 }
02044
02045 MC_mesh_index_vector MC_mesh_index_vector::load_mesh_file(const std::string& filename)
02046 {
02047 MC_mesh_index_vector mesh;
02048
02049
02050 std::vector<std::string> token=MC_string_tokenizer::tokenize(filename,".");
02051 std::string extension=token[token.size()-1];
02052
02053 if(extension.find("off")!=std::string::npos)
02054 mesh=MC_io_off::read_off_file(filename);
02055 else if(extension.find("obj")!=std::string::npos)
02056 mesh=MC_io_obj::read_obj_file(filename);
02057 else
02058 std::cout<<"Warning in MC_mesh_index_vector::load_mesh_file("<<filename<<"), extension "<<extension<<" not recognized"<<std::endl;
02059
02060 return mesh;
02061 }
02062
02063 double MC_mesh_index_vector::area() const
02064 {
02065 double a=0.0;
02066 for(int k=0,N=polygon_number();k<N;++k)
02067 {
02068 MC_int_vector poly=connectivity_mesh(k);
02069
02070 MC_v3d x0=point_set_mesh(poly[0]);
02071 for(int k_tri=0,N_pol=poly.size();k_tri<N_pol-2;++k_tri)
02072 {
02073 MC_v3d x1=point_set_mesh(poly[k_tri+1]);
02074 MC_v3d x2=point_set_mesh(poly[k_tri+2]);
02075
02076 double C0=(x2[1]-x0[1])*(x1[2]-x0[2])-(x2[2]-x0[2])*(x1[1]-x0[1]);
02077 double C1=(x2[2]-x0[2])*(x1[0]-x0[0])-(x2[0]-x0[0])*(x1[2]-x0[2]);
02078 double C2=(x2[0]-x0[0])*(x1[1]-x0[1])-(x2[1]-x0[1])*(x1[0]-x0[0]);
02079
02080 double u=sqrt(C0*C0+C1*C1+C2*C2);
02081
02082 a+=0.5*u;
02083 }
02084 }
02085 return a;
02086 }
02087
02088 MC_v3d_vector MC_mesh_index_vector::area_gradient() const
02089 {
02090 MC_v3d_vector grad=MC_v3d_vector::zeros(vertex_number());
02091
02092 for(int k=0,N=polygon_number();k<N;++k)
02093 {
02094 MC_int_vector poly=connectivity_mesh(k);
02095
02096 MC_v3d x0=point_set_mesh(poly[0]);
02097 int index_0=poly[0];
02098 for(int k_tri=0,N_pol=poly.size();k_tri<N_pol-2;++k_tri)
02099 {
02100 MC_v3d x1=point_set_mesh(poly[k_tri+1]);
02101 MC_v3d x2=point_set_mesh(poly[k_tri+2]);
02102
02103
02104 int index_1=poly[k_tri+1];
02105 int index_2=poly[k_tri+2];
02106
02107 double C0=(x2[1]-x0[1])*(x1[2]-x0[2])-(x2[2]-x0[2])*(x1[1]-x0[1]);
02108 double C1=(x2[2]-x0[2])*(x1[0]-x0[0])-(x2[0]-x0[0])*(x1[2]-x0[2]);
02109 double C2=(x2[0]-x0[0])*(x1[1]-x0[1])-(x2[1]-x0[1])*(x1[0]-x0[0]);
02110
02111 double u=sqrt(C0*C0+C1*C1+C2*C2);
02112
02113
02114 grad[index_0][0] += ( (x1[2]-x2[2])*C1 + (x2[1]-x1[1])*C2 )/(2*u);
02115 grad[index_0][1] += ( (x2[2]-x1[2])*C0 + (x1[0]-x2[0])*C2 )/(2*u);
02116 grad[index_0][2] += ( (x1[1]-x2[1])*C0 + (x2[0]-x1[0])*C1 )/(2*u);
02117
02118 grad[index_1][0] += ( (x2[2]-x0[2])*C1 + (x0[1]-x2[1])*C2 )/(2*u);
02119 grad[index_1][1] += ( (x0[2]-x2[2])*C0 + (x2[0]-x0[0])*C2 )/(2*u);
02120 grad[index_1][2] += ( (x2[1]-x0[1])*C0 + (x0[0]-x2[0])*C1 )/(2*u);
02121
02122 grad[index_2][0] += ( (x0[2]-x1[2])*C1 + (x1[1]-x0[1])*C2 )/(2*u);
02123 grad[index_2][1] += ( (x1[2]-x0[2])*C0 + (x0[0]-x1[0])*C2 )/(2*u);
02124 grad[index_2][2] += ( (x0[1]-x1[1])*C0 + (x1[0]-x0[0])*C1 )/(2*u);
02125 }
02126 }
02127
02128 return grad;
02129 }
02130
02131
02132 MC_mesh_index_vector MC_mesh_index_vector::build_ball_point_set(const MC_v3d_vector& v,const double& radius,const int& N_subdiv)
02133 {
02134 MC_mesh_index_vector sph=MC_mesh_index_vector::build_sphere(N_subdiv)*radius;
02135 MC_mesh_index_vector visu_point_set;
02136 for(int k=0,N=v.size();k<N;++k)
02137 visu_point_set.concatenation(sph+v[k]);
02138 return visu_point_set;
02139 }
02140
02141 double MC_mesh_index_vector::average_edge_length() const
02142 {
02143 int counter=0;
02144 double avg_length=0.0;
02145 for(int k=0,N=polygon_number();k<N;++k)
02146 {
02147 MC_int_vector index=connectivity()(k);
02148 for(int k_v=0,N_v=index.size();k_v<N_v;++k_v)
02149 {
02150 avg_length += ( point_set()(index(k_v))-point_set()(index( (k_v+1)%N_v )) ).norm();
02151 ++counter;
02152 }
02153 }
02154 return avg_length/static_cast<double>(counter);
02155 }
02156
02157 MC_mesh_index_vector MC_mesh_index_vector::get_polygon_mesh() const
02158 {
02159 MC_mesh_index_vector exploded;
02160 for(int k=0,N=polygon_number();k<N;++k)
02161 {
02162 MC_polygon p=get_polygon(k);
02163 exploded.connectivity().add(MC_int_vector::linspace(0,p.size()-1)+exploded.vertex_number());
02164 exploded.point_set().add(p);
02165 }
02166 return exploded;
02167 }
02168
02169
02170
02171
02172 MC_mesh_index_vector MC_mesh_index_vector::build_wireframe(const MC_mesh_index_vector& mesh,const double& radius_cylinder,const int& N_circular,const int& N_subdiv_edges)
02173 {
02174 std::set<MC_int_pair,MC_int_pair_less> already_created_edge;
02175
02176 MC_curve pattern=MC_curve::build_circle(radius_cylinder,N_circular);
02177 MC_mesh_index_vector wireframe;
02178
02179 for(int k=0,N=mesh.polygon_number();k<N;++k)
02180 {
02181 MC_int_vector poly=mesh.connectivity()(k);
02182 for(int k2=0,N_poly=poly.size();k2<N_poly;++k2)
02183 {
02184 int index_0=poly[k2];
02185 int index_1=poly[(k2+1)%N_poly];
02186
02187 if(already_created_edge.insert(MC_int_pair(index_0,index_1)).second==true)
02188 {
02189 MC_v3d x0=mesh.point_set()(poly[k2]);
02190 MC_v3d x1=mesh.point_set()(poly[(k2+1)%N_poly]);
02191
02192 MC_curve c=MC_curve::build_line(x0,x1,N_subdiv_edges);
02193
02194 MC_mesh_index_vector sweep=MC_mesh_index_vector::sweep_surface(c,pattern,2,false);
02195 wireframe.concatenation(sweep);
02196 }
02197 }
02198 }
02199
02200 return wireframe;
02201 }
02202
02203 MC_matrix MC_mesh_index_vector::inertia() const
02204 {
02205 MC_matrix J(3);
02206 for(int k=0,N=polygon_number();k<N;++k)
02207 {
02208 std::vector<MC_triangle> tri=MC_triangle::triangulate(get_polygon(k));
02209 for(int k_tri=0,N_tri=tri.size();k_tri<N_tri;++k_tri)
02210 J+=tri[k_tri].inertia();
02211 }
02212
02213 return J;
02214 }
02215
02216 MC_mesh_index_vector MC_mesh_index_vector::build_sphere_point_set(const MC_v3d_vector& position,const double& radius,const unsigned int& N_subdiv)
02217 {
02218 MC_mesh_index_vector sph=radius*MC_mesh_index_vector::build_sphere(N_subdiv);
02219 MC_mesh_index_vector v_sph;
02220 for(int k=0,N=position.size();k<N;++k)
02221 v_sph.concatenation(sph+position[k]);
02222
02223 return v_sph;
02224 }
02225
02226
02227
02228 MC_curve MC_mesh_index_vector::intersection_curve(const MC_mesh_index_vector& mesh,const MC_segment& s)
02229 {
02230
02231 MC_v3d_vector normal_mesh=mesh.get_polygon(0).normal();
02232 MC_v3d_vector normal_plane=(s[1]-s[0]).normalized().cross(normal_mesh).normalized();
02233
02234
02235
02236 std::vector <MC_curve> curve_inter_plane=mesh.plane_intersection(normal_plane,s[0]);
02237
02238
02239
02240 MC_io_off::write_off_file("mesh_to_cut.off",mesh);
02241 std::vector<MC_segment> v_s;v_s.push_back(s);
02242 MC_io_off::write_off_file("line_in_mesh.off",MC_mesh_index_vector::build_segment(v_s,0.005,10));
02243
02244 std::cout<<" ! "<<curve_inter_plane.size()<<std::endl;
02245
02246 std::cout<<"-- "<<curve_inter_plane[0]<<std::endl;;
02247
02248 if(curve_inter_plane.size()<1)
02249 {std::cout<<"Something strange in MC_mesh_index_vector::intersection_curve()"<<std::endl;exit(-1);}
02250
02251 return curve_inter_plane[0];
02252
02253 }
02254
02255 MC_mesh_index_vector MC_mesh_index_vector::build_segment(const std::vector<MC_segment>& v_seg,const double& radius,const double& N_circular)
02256 {
02257 MC_curve pattern=MC_curve::build_circle(radius,N_circular);
02258 MC_mesh_index_vector wireframe;
02259 unsigned int N_subdiv_edges=2;
02260 for(int k=0,N=v_seg.size();k<N;++k)
02261 {
02262 const MC_v3d& x0=v_seg[k][0];
02263 const MC_v3d& x1=v_seg[k][1];
02264
02265 MC_curve c=MC_curve::build_line(x0,x1,N_subdiv_edges);
02266 MC_mesh_index_vector sweep=MC_mesh_index_vector::sweep_surface(c,pattern,8,false);
02267 wireframe.concatenation(sweep);
02268 }
02269 return wireframe;
02270 }
02271
02272 MC_mesh_index_vector MC_mesh_index_vector::texture_convert_planar_xy(const MC_mesh_index_vector& mesh)
02273 {
02274
02275 MC_v3d_vector texture_coord=mesh.point_set();
02276 texture_coord=texture_coord.scaled_to_unit_and_center().first;
02277 return MC_mesh_index_vector(texture_coord,mesh.connectivity());
02278 }
02279 MC_mesh_index_vector MC_mesh_index_vector::texture_convert_planar_xy_same_connectivity(const MC_mesh_index_vector& mesh_tex,const MC_mesh_index_vector& mesh_ref)
02280 {
02281 MC_v3d_vector new_texture;
02282 MC_mesh_index_vector tex=texture_convert_planar_xy(mesh_tex);
02283 for(unsigned int k_poly=0,N_poly=mesh_ref.polygon_number();k_poly<N_poly;++k_poly)
02284 for(unsigned int k=0,N=mesh_ref.connectivity()[k_poly].size();k<N;++k)
02285 new_texture.set( mesh_ref.connectivity()[k_poly][k],tex.point_set()(tex.connectivity()[k_poly][k]) );
02286 tex=MC_mesh_index_vector(new_texture,mesh_ref.connectivity());
02287 return tex;
02288 }
02289
02290 std::pair<MC_mesh_index_vector,std::pair<MC_int_vector_vector,MC_int_vector_vector> > MC_mesh_index_vector::build_from_polygon_soup(const std::vector<MC_polygon> poly_soup)
02291 {
02292 MC_mesh_index_vector m;
02293 return m.added_polygon_soup(poly_soup);
02294 }
02295
02296 MC_mesh_index_vector MC_mesh_index_vector::build_square(const int& N1,const int& N2,const MC_v3d& tangent_1,const MC_v3d& tangent_2,const MC_v3d& center,const double& L1,const double& L2)
02297 {
02298 MC_mesh_index_vector m=build_square(N1,N2);
02299
02300
02301 m.point_set() += MC_v3d(-0.5,-0.5,0.0);
02302
02303
02304 m.point_set()=m.point_set().scale(MC_v3d(L1,L2,1.0));
02305
02306
02307
02308
02309 MC_matrix R1=MC_matrix::rotation_axis_to_axis(MC_v3d(1,0,0),tangent_1);
02310 MC_matrix R2=MC_matrix::rotation_axis_to_axis(R1*MC_v3d(0,1,0),tangent_2);
02311 m=R2*R1*m;
02312
02313
02314 m.point_set() += center;
02315
02316 return m;
02317 }
02318
02319 }
