00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054 #include<vector>
00055 #include <algorithm>
00056 #include<vcg/space/point3.h>
00057 #include<vcg/space/plane3.h>
00058 #include<vcg/space/segment3.h>
00059 #include<vcg/space/intersection3.h>
00060 #include<vcg/complex/edgemesh/allocate.h>
00061 #include<vcg/complex/trimesh/allocate.h>
00062 #include<vcg/complex/trimesh/subset.h>
00063 #include<vcg/complex/trimesh/closest.h>
00064 #include<vcg/complex/trimesh/base.h>
00065
00066 #ifndef __VCGLIB_INTERSECTION_TRI_MESH
00067 #define __VCGLIB_INTERSECTION_TRI_MESH
00068
00069 namespace vcg{
00070
00076 template < typename GridType,typename ScalarType>
00077 bool Intersect( GridType & grid,Plane3<ScalarType> plane, std::vector<typename GridType::Cell *> &cells){
00078 Point3d p,_d;
00079 Plane3d pl;
00080 _d.Import(plane.Direction());
00081 pl.SetDirection(_d);
00082 pl.SetOffset(plane.Offset());
00083 for( int ax = 0; ax <3; ++ax)
00084 { int axis = ax;
00085 int axis0 = (axis+1)%3;
00086 int axis1 = (axis+2)%3;
00087 int i,j;
00088 Point3i pi;
00089
00090 Segment3<double> seg;
00091 seg.P0().Import(grid.bbox.min);
00092 seg.P1().Import(grid.bbox.min);
00093 seg.P1()[axis] = grid.bbox.max[axis];
00094
00095 for(i = 0 ; i <= grid.siz[axis0]; ++i){
00096 for(j = 0 ; j <= grid.siz[axis1]; ++j)
00097 {
00098 seg.P0()[axis0] = grid.bbox.min[axis0]+ (i+0.01) * grid.voxel[axis0] ;
00099 seg.P1()[axis0] = grid.bbox.min[axis0]+ (i+0.01) * grid.voxel[axis0];
00100 seg.P0()[axis1] = grid.bbox.min[axis1]+ (j+0.01) * grid.voxel[axis1];
00101 seg.P1()[axis1] = grid.bbox.min[axis1]+ (j+0.01) * grid.voxel[axis1];
00102 if ( IntersectionPlaneSegmentEpsilon(pl,seg,p))
00103 {
00104 pi[axis] = std::min(std::max(0,(int)floor((p[axis ]-grid.bbox.min[axis])/grid.voxel[axis])),grid.siz[axis]);
00105 pi[axis0] = i;
00106 pi[axis1] = j;
00107 grid.Grid(pi,axis,cells);
00108 }
00109 }
00110 }
00111 }
00112 sort(cells.begin(),cells.end());
00113 cells.erase(unique(cells.begin(),cells.end()),cells.end());
00114
00115 return false;
00116 }
00117
00123 template < typename TriMeshType, typename EdgeMeshType, class ScalarType, class IndexingType >
00124 bool Intersection(
00125 Plane3<ScalarType> pl,
00126 EdgeMeshType & em,
00127 double& ave_length,
00128 IndexingType *grid,
00129 typename std::vector< typename IndexingType::Cell* >& cells)
00130 {
00131 typedef typename TriMeshType::FaceContainer FaceContainer;
00132 typedef IndexingType GridType;
00133 typename EdgeMeshType::VertexIterator vi;
00134 typename TriMeshType::FaceIterator fi;
00135 std::vector<typename TriMeshType::FaceType*> v;
00136 v.clear();
00137 Intersect(*grid,pl,cells);
00138 Segment3<ScalarType> seg;
00139 ave_length = 0.0;
00140 typename std::vector<typename GridType::Cell*>::iterator ic;
00141 typename GridType::Cell fs,ls;
00142 for(ic = cells.begin(); ic != cells.end();++ic)
00143 {
00144 grid->Grid(*ic,fs,ls);
00145 typename GridType::Link * lk = fs;
00146 while(lk != ls){
00147 typename TriMeshType::FaceType & face = *(lk->Elem());
00148 if(!face.IsS())
00149 {
00150 face.SetS();
00151 v.push_back(&face);
00152 if(vcg::IntersectionPlaneTriangle(pl,face,seg))
00153 {
00154 face.SetS();
00155
00156 ave_length+=seg.Length();
00157 vcg::edg::Allocator<EdgeMeshType>::AddEdges(em,1);
00158 vi = vcg::edg::Allocator<EdgeMeshType>::AddVertices(em,2);
00159 (*vi).P() = seg.P0();
00160 em.edges.back().V(0) = &(*vi);
00161 vi++;
00162 (*vi).P() = seg.P1();
00163 em.edges.back().V(1) = &(*vi);
00164 }
00165 }
00166 lk++;
00167 }
00168 }
00169 ave_length/=em.en;
00170 typename std::vector<typename TriMeshType::FaceType*>::iterator v_i;
00171 for(v_i=v.begin(); v_i!=v.end(); ++v_i) (*v_i)->ClearS();
00172
00173 return true;
00174 }
00175
00186 template < typename TriMeshType, typename EdgeMeshType, class ScalarType >
00187 bool Intersection(TriMeshType & m,
00188 Plane3<ScalarType> pl,
00189 EdgeMeshType & em)
00190 {
00191 typename EdgeMeshType::VertexIterator vi;
00192 typename TriMeshType::FaceIterator fi;
00193 em.Clear();
00194 Segment3<ScalarType> seg;
00195 for(fi=m.face.begin();fi!=m.face.end();++fi)
00196 if(!(*fi).IsD())
00197 {
00198 if(vcg::IntersectionPlaneTriangle(pl,*fi,seg))
00199 {
00200 vcg::edg::Allocator<EdgeMeshType>::AddEdges(em,1);
00201 vi = vcg::edg::Allocator<EdgeMeshType>::AddVertices(em,2);
00202 (*vi).P() = seg.P0();
00203 em.edges.back().V(0) = &(*vi);
00204 vi++;
00205 (*vi).P() = seg.P1();
00206 em.edges.back().V(1) = &(*vi);
00207 }
00208 }
00209
00210 return true;
00211 }
00212
00220 template < typename TriMeshType, class ScalarType, class IndexingType >
00221 bool Intersection(Plane3<ScalarType> pl,
00222 IndexingType *grid,
00223 typename std::vector<typename TriMeshType::FaceType*> &v)
00224 {
00225 typedef typename TriMeshType::FaceContainer FaceContainer;
00226 typedef IndexingType GridType;
00227 typename TriMeshType::FaceIterator fi;
00228 v.clear();
00229 typename std::vector< typename GridType::Cell* > cells;
00230 Intersect(*grid,pl,cells);
00231 typename std::vector<typename GridType::Cell*>::iterator ic;
00232 typename GridType::Cell fs,ls;
00233
00234 for(ic = cells.begin(); ic != cells.end();++ic)
00235 {
00236 grid->Grid(*ic,fs,ls);
00237 typename GridType::Link * lk = fs;
00238 while(lk != ls){
00239 typename TriMeshType::FaceType & face = *(lk->Elem());
00240 v.push_back(&face);
00241 lk++;
00242 }
00243 }
00244 return true;
00245 }
00246
00250 template < typename TriMeshType, class ScalarType>
00251 bool IntersectionRayMesh(
00252 TriMeshType * m,
00253 const Line3<ScalarType> & ray,
00254 Point3<ScalarType> & hitPoint)
00255 {
00256
00257 typename TriMeshType::FaceIterator fi;
00258 bool hit=false;
00259
00260 if(m==0) return false;
00261
00262
00263
00264
00265 ScalarType bar1,bar2,dist;
00266 Point3<ScalarType> p1;
00267 Point3<ScalarType> p2;
00268 Point3<ScalarType> p3;
00269 for(fi = m->face.begin(); fi != m->face.end(); ++fi)
00270 {
00271 p1=vcg::Point3<ScalarType>( (*fi).P(0).X() ,(*fi).P(0).Y(),(*fi).P(0).Z() );
00272 p2=vcg::Point3<ScalarType>( (*fi).P(1).X() ,(*fi).P(1).Y(),(*fi).P(1).Z() );
00273 p3=vcg::Point3<ScalarType>( (*fi).P(2).X() ,(*fi).P(2).Y(),(*fi).P(2).Z() );
00274 if(IntersectionLineTriangle<ScalarType>(ray,p1,p2,p3,bar1,bar2,dist))
00275 {
00276 hitPoint= p1*(1-bar1-bar2) + p2*bar1 + p3*bar2;
00277 hit=true;
00278 }
00279 }
00280
00281 return hit;
00282 }
00283
00288 template < typename TriMeshType, class ScalarType>
00289 bool IntersectionRayMesh(
00290 TriMeshType * m,
00291 const Line3<ScalarType> & ray,
00292 Point3<ScalarType> & hitPoint,
00293 ScalarType &bar1,
00294 ScalarType &bar2,
00295 ScalarType &bar3,
00296 typename TriMeshType::FacePointer fp
00297 )
00298 {
00299
00300 typename TriMeshType::FaceIterator fi;
00301 bool hit=false;
00302
00303 if(m==0) return false;
00304
00305
00306
00307
00308 ScalarType dist;
00309 Point3<ScalarType> p1;
00310 Point3<ScalarType> p2;
00311 Point3<ScalarType> p3;
00312 for(fi = m->face.begin(); fi != m->face.end(); ++fi)
00313 {
00314 p1=vcg::Point3<ScalarType>( (*fi).P(0).X() ,(*fi).P(0).Y(),(*fi).P(0).Z() );
00315 p2=vcg::Point3<ScalarType>( (*fi).P(1).X() ,(*fi).P(1).Y(),(*fi).P(1).Z() );
00316 p3=vcg::Point3<ScalarType>( (*fi).P(2).X() ,(*fi).P(2).Y(),(*fi).P(2).Z() );
00317 if(IntersectionLineTriangle<ScalarType>(ray,p1,p2,p3,bar1,bar2,dist))
00318 {
00319 bar3 = (1-bar1-bar2);
00320 hitPoint= p1*bar3 + p2*bar1 + p3*bar2;
00321 fp = &(*fi);
00322 hit=true;
00323 }
00324 }
00325
00326 return hit;
00327 }
00328
00337 template < typename TriMeshType, class ScalarType>
00338 void IntersectionBallMesh( TriMeshType & m, const vcg::Sphere3<ScalarType> &ball, TriMeshType & res,
00339 float tol = 0){
00340
00341 typename TriMeshType::VertexIterator v0,v1,v2;
00342 typename TriMeshType::FaceIterator fi;
00343 std::vector<typename TriMeshType:: FaceType*> closests;
00344 vcg::Point3<ScalarType> witness;
00345 std::pair<ScalarType, ScalarType> info;
00346
00347 if(tol == 0) tol = M_PI * ball.Radius() * ball.Radius() / 100000;
00348
00349 for(fi = m.face.begin(); fi != m.face.end(); ++fi)
00350 if(!(*fi).IsD() && IntersectionSphereTriangle<ScalarType>(ball ,(*fi), witness , &info))
00351 closests.push_back(&(*fi));
00352
00353 res.Clear();
00354 SubSet(res,closests);
00355 int i =0;
00356 while(i<res.fn){
00357 bool allIn = ( ball.IsIn(res.face[i].P(0)) && ball.IsIn(res.face[i].P(1))&&ball.IsIn(res.face[i].P(2)));
00358 if( IntersectionSphereTriangle<ScalarType>(ball ,res.face[i], witness , &info) && !allIn){
00359 if(vcg::DoubleArea(res.face[i]) > tol)
00360 {
00361
00362 v0 = vcg::tri::Allocator<TriMeshType>::AddVertices(res,3);
00363 fi = vcg::tri::Allocator<TriMeshType>::AddFaces(res,4);
00364
00365 v1 = v0; ++v1;
00366 v2 = v1; ++v2;
00367 (*v0).P() = (res.face[i].P(0) + res.face[i].P(1))*0.5;
00368 (*v1).P() = (res.face[i].P(1) + res.face[i].P(2))*0.5;
00369 (*v2).P() = (res.face[i].P(2) + res.face[i].P(0))*0.5;
00370
00371 (*fi).V(0) = res.face[i].V(0);
00372 (*fi).V(1) = &(*v0);
00373 (*fi).V(2) = &(*v2);
00374 ++fi;
00375
00376 (*fi).V(0) = res.face[i].V(1);
00377 (*fi).V(1) = &(*v1);
00378 (*fi).V(2) = &(*v0);
00379 ++fi;
00380
00381 (*fi).V(0) = &(*v0);
00382 (*fi).V(1) = &(*v1);
00383 (*fi).V(2) = &(*v2);
00384 ++fi;
00385
00386 (*fi).V(0) = &(*v2);
00387 (*fi).V(1) = &(*v1);
00388 (*fi).V(2) = res.face[i].V(2) ;
00389
00390 vcg::tri::Allocator<TriMeshType>::DeleteFace(res,res.face[i]);
00391 }
00392 }
00393
00394 if(info.first > 0.0)
00395 vcg::tri::Allocator<TriMeshType>::DeleteFace(res,res.face[i]);
00396 ++i;
00397 }
00398 }
00399
00400
00401 template < typename TriMeshType, class ScalarType, class IndexingType>
00402 void IntersectionBallMesh( IndexingType * grid, TriMeshType & m, const vcg::Sphere3<ScalarType> &ball, TriMeshType & res,
00403 float tol = 0){
00404
00405 typename TriMeshType::VertexIterator v0,v1,v2;
00406 typename std::vector<typename TriMeshType::FacePointer >::iterator cfi;
00407 typename TriMeshType::FaceIterator fi;
00408 std::vector<typename TriMeshType:: FaceType*> closestsF,closests;
00409 vcg::Point3<ScalarType> witness;
00410 std::vector<vcg::Point3<ScalarType> > witnesses;
00411 std::vector<ScalarType> distances;
00412 std::pair<ScalarType, ScalarType> info;
00413
00414 if(tol == 0) tol = M_PI * ball.Radius() * ball.Radius() / 100000;
00415
00416 vcg::tri::GetInSphereFace(m,*grid, ball.Center(), ball.Radius(),closestsF,distances,witnesses);
00417 for(cfi =closestsF.begin(); cfi != closestsF.end(); ++cfi)
00418 if(!(**cfi).IsD() && IntersectionSphereTriangle<ScalarType>(ball ,(**cfi), witness , &info))
00419 closests.push_back(&(**cfi));
00420
00421 res.Clear();
00422 SubSet(res,closests);
00423 int i =0;
00424 while(i<res.fn){
00425 bool allIn = ( ball.IsIn(res.face[i].P(0)) && ball.IsIn(res.face[i].P(1))&&ball.IsIn(res.face[i].P(2)));
00426 if( IntersectionSphereTriangle<ScalarType>(ball ,res.face[i], witness , &info) && !allIn){
00427 if(vcg::DoubleArea(res.face[i]) > tol)
00428 {
00429
00430 v0 = vcg::tri::Allocator<TriMeshType>::AddVertices(res,3);
00431 fi = vcg::tri::Allocator<TriMeshType>::AddFaces(res,4);
00432
00433 v1 = v0; ++v1;
00434 v2 = v1; ++v2;
00435 (*v0).P() = (res.face[i].P(0) + res.face[i].P(1))*0.5;
00436 (*v1).P() = (res.face[i].P(1) + res.face[i].P(2))*0.5;
00437 (*v2).P() = (res.face[i].P(2) + res.face[i].P(0))*0.5;
00438
00439 (*fi).V(0) = res.face[i].V(0);
00440 (*fi).V(1) = &(*v0);
00441 (*fi).V(2) = &(*v2);
00442 ++fi;
00443
00444 (*fi).V(0) = res.face[i].V(1);
00445 (*fi).V(1) = &(*v1);
00446 (*fi).V(2) = &(*v0);
00447 ++fi;
00448
00449 (*fi).V(0) = &(*v0);
00450 (*fi).V(1) = &(*v1);
00451 (*fi).V(2) = &(*v2);
00452 ++fi;
00453
00454 (*fi).V(0) = &(*v2);
00455 (*fi).V(1) = &(*v1);
00456 (*fi).V(2) = res.face[i].V(2) ;
00457
00458 vcg::tri::Allocator<TriMeshType>::DeleteFace(res,res.face[i]);
00459 }
00460 }
00461
00462 if(info.first > 0.0)
00463 vcg::tri::Allocator<TriMeshType>::DeleteFace(res,res.face[i]);
00464 ++i;
00465 }
00466 }
00467
00469 }
00470 #endif
