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 #ifndef __VCG_EXTENDED_MARCHING_CUBES
00027 #define __VCG_EXTENDED_MARCHING_CUBES
00028
00029 #include <float.h>
00030 #include <assert.h>
00031 #include <vector>
00032 #include <vcg/math/base.h>
00033 #include <vcg/math/matrix.h>
00034 #include <vcg/math/lin_algebra.h>
00035 #include <vcg/simplex/face/topology.h>
00036 #include <vcg/complex/trimesh/update/edges.h>
00037 #include <vcg/complex/trimesh/update/normal.h>
00038 #include <vcg/complex/trimesh/update/topology.h>
00039 #include <vcg/complex/trimesh/allocate.h>
00040 #include <vcg/space/point3.h>
00041 #include "emc_lookup_table.h"
00042
00043 namespace vcg
00044 {
00045 namespace tri
00046 {
00047
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00066
00073 template<class TRIMESH_TYPE, class WALKER_TYPE>
00074 class ExtendedMarchingCubes
00075 {
00076 public:
00077 #if defined(__GNUC__)
00078 typedef unsigned int size_t;
00079 #else
00080 #ifdef _WIN64
00081 typedef unsigned __int64 size_t;
00082 #else
00083 typedef _W64 unsigned int size_t;
00084 #endif
00085 #endif
00086 typedef typename vcg::tri::Allocator< TRIMESH_TYPE > AllocatorType;
00087 typedef typename TRIMESH_TYPE::ScalarType ScalarType;
00088 typedef typename TRIMESH_TYPE::VertexType VertexType;
00089 typedef typename TRIMESH_TYPE::VertexPointer VertexPointer;
00090 typedef typename TRIMESH_TYPE::VertexIterator VertexIterator;
00091 typedef typename TRIMESH_TYPE::FaceType FaceType;
00092 typedef typename TRIMESH_TYPE::FacePointer FacePointer;
00093 typedef typename TRIMESH_TYPE::FaceIterator FaceIterator;
00094 typedef typename TRIMESH_TYPE::CoordType CoordType;
00095 typedef typename TRIMESH_TYPE::CoordType* CoordPointer;
00096
00097 typedef struct
00098 {
00099 size_t face, edge;
00100 } LightEdge;
00101
00109 ExtendedMarchingCubes(TRIMESH_TYPE &mesh, WALKER_TYPE &walker, ScalarType angle=30)
00110 {
00111 _mesh = &mesh;
00112 _walker = &walker;
00113 _featureAngle = vcg::math::ToRad(angle);
00114 _initialized = _finalized = false;
00115 };
00116
00121 void Initialize()
00122 {
00123 assert(!_initialized && !_finalized);
00124 _featureFlag = VertexType::NewBitFlag();
00125 _initialized = true;
00126 };
00127
00132 void Finalize()
00133 {
00134 assert(_initialized && !_finalized);
00135 FlipEdges();
00136
00137 VertexIterator v_iter = _mesh->vert.begin();
00138 VertexIterator v_end = _mesh->vert.end();
00139 for ( ; v_iter!=v_end; v_iter++)
00140 v_iter->ClearUserBit( _featureFlag );
00141 VertexType::DeleteBitFlag( _featureFlag );
00142 _featureFlag = 0;
00143 _mesh = NULL;
00144 _walker = NULL;
00145 _finalized = true;
00146 };
00147
00154 void ProcessCell(const vcg::Point3i &min, const vcg::Point3i &max)
00155 {
00156 assert(_initialized && !_finalized);
00157 assert(min[0]<max[0] && min[1]<max[1] && min[2]<max[2]);
00158 _corners[0].X()=min.X(); _corners[0].Y()=min.Y(); _corners[0].Z()=min.Z();
00159 _corners[1].X()=max.X(); _corners[1].Y()=min.Y(); _corners[1].Z()=min.Z();
00160 _corners[2].X()=max.X(); _corners[2].Y()=max.Y(); _corners[2].Z()=min.Z();
00161 _corners[3].X()=min.X(); _corners[3].Y()=max.Y(); _corners[3].Z()=min.Z();
00162 _corners[4].X()=min.X(); _corners[4].Y()=min.Y(); _corners[4].Z()=max.Z();
00163 _corners[5].X()=max.X(); _corners[5].Y()=min.Y(); _corners[5].Z()=max.Z();
00164 _corners[6].X()=max.X(); _corners[6].Y()=max.Y(); _corners[6].Z()=max.Z();
00165 _corners[7].X()=min.X(); _corners[7].Y()=max.Y(); _corners[7].Z()=max.Z();
00166
00167 unsigned char cubetype = 0;
00168 if ((_field[0] = _walker->V(_corners[0].X(), _corners[0].Y(), _corners[0].Z())) >= 0) cubetype+= 1;
00169 if ((_field[1] = _walker->V(_corners[1].X(), _corners[1].Y(), _corners[1].Z())) >= 0) cubetype+= 2;
00170 if ((_field[2] = _walker->V(_corners[2].X(), _corners[2].Y(), _corners[2].Z())) >= 0) cubetype+= 4;
00171 if ((_field[3] = _walker->V(_corners[3].X(), _corners[3].Y(), _corners[3].Z())) >= 0) cubetype+= 8;
00172 if ((_field[4] = _walker->V(_corners[4].X(), _corners[4].Y(), _corners[4].Z())) >= 0) cubetype+= 16;
00173 if ((_field[5] = _walker->V(_corners[5].X(), _corners[5].Y(), _corners[5].Z())) >= 0) cubetype+= 32;
00174 if ((_field[6] = _walker->V(_corners[6].X(), _corners[6].Y(), _corners[6].Z())) >= 0) cubetype+= 64;
00175 if ((_field[7] = _walker->V(_corners[7].X(), _corners[7].Y(), _corners[7].Z())) >= 0) cubetype+=128;
00176
00177 if (cubetype==0 || cubetype==255)
00178 return;
00179
00180 size_t vertices_idx[12];
00181 memset(vertices_idx, -1, 12*sizeof(size_t));
00182 int code = EMCLookUpTable::EdgeTable(cubetype);
00183 VertexPointer vp = NULL;
00184 if ( 1&code ) { _walker->GetXIntercept(_corners[0], _corners[1], vp); vertices_idx[ 0] = vp - &_mesh->vert[0]; }
00185 if ( 2&code ) { _walker->GetYIntercept(_corners[1], _corners[2], vp); vertices_idx[ 1] = vp - &_mesh->vert[0]; }
00186 if ( 4&code ) { _walker->GetXIntercept(_corners[3], _corners[2], vp); vertices_idx[ 2] = vp - &_mesh->vert[0]; }
00187 if ( 8&code ) { _walker->GetYIntercept(_corners[0], _corners[3], vp); vertices_idx[ 3] = vp - &_mesh->vert[0]; }
00188 if ( 16&code ) { _walker->GetXIntercept(_corners[4], _corners[5], vp); vertices_idx[ 4] = vp - &_mesh->vert[0]; }
00189 if ( 32&code ) { _walker->GetYIntercept(_corners[5], _corners[6], vp); vertices_idx[ 5] = vp - &_mesh->vert[0]; }
00190 if ( 64&code ) { _walker->GetXIntercept(_corners[7], _corners[6], vp); vertices_idx[ 6] = vp - &_mesh->vert[0]; }
00191 if ( 128&code ) { _walker->GetYIntercept(_corners[4], _corners[7], vp); vertices_idx[ 7] = vp - &_mesh->vert[0]; }
00192 if ( 256&code ) { _walker->GetZIntercept(_corners[0], _corners[4], vp); vertices_idx[ 8] = vp - &_mesh->vert[0]; }
00193 if ( 512&code ) { _walker->GetZIntercept(_corners[1], _corners[5], vp); vertices_idx[ 9] = vp - &_mesh->vert[0]; }
00194 if (1024&code ) { _walker->GetZIntercept(_corners[2], _corners[6], vp); vertices_idx[10] = vp - &_mesh->vert[0]; }
00195 if (2048&code ) { _walker->GetZIntercept(_corners[3], _corners[7], vp); vertices_idx[11] = vp - &_mesh->vert[0]; }
00196
00197 int m, n, vertices_num;
00198 int components = EMCLookUpTable::TriTable(cubetype, 1)[0];
00199 int *indices = &EMCLookUpTable::TriTable(cubetype, 1)[components+1];
00200
00201 std::vector< size_t > vertices_list;
00202 for (m=1; m<=components; m++)
00203 {
00204
00205 vertices_num = EMCLookUpTable::TriTable(cubetype, 1)[m];
00206
00207
00208 vertices_list.clear();
00209 for (n=0; n<vertices_num; ++n)
00210 vertices_list.push_back( vertices_idx[ indices[n] ] );
00211
00212 VertexPointer feature = FindFeature( vertices_list );
00213 if (feature != NULL)
00214 {
00215
00216 size_t feature_idx = feature - &_mesh->vert[0];
00217 size_t face_idx = _mesh->face.size();
00218 vertices_list.push_back( vertices_list[0] );
00219 AllocatorType::AddFaces(*_mesh, (int) vertices_num);
00220 for (int j=0; j<vertices_num; ++j, face_idx++)
00221 {
00222 _mesh->face[face_idx].V(0) = &_mesh->vert[ vertices_list[j ] ];
00223 _mesh->face[face_idx].V(1) = &_mesh->vert[ vertices_list[j+1] ];
00224 _mesh->face[face_idx].V(2) = &_mesh->vert[ feature_idx ];
00225 }
00226 }
00227 else
00228 {
00229
00230 for (int j=0; EMCLookUpTable::PolyTable(vertices_num, j) != -1; j+=3)
00231 {
00232 size_t face_idx = _mesh->face.size();
00233 AllocatorType::AddFaces(*_mesh, 1);
00234
00235
00236
00237 _mesh->face[ face_idx].V(0) = &_mesh->vert[ vertices_idx[ indices[ EMCLookUpTable::PolyTable(vertices_num, j ) ] ] ];
00238 _mesh->face[ face_idx].V(1) = &_mesh->vert[ vertices_idx[ indices[ EMCLookUpTable::PolyTable(vertices_num, j+1) ] ] ];
00239 _mesh->face[ face_idx].V(2) = &_mesh->vert[ vertices_idx[ indices[ EMCLookUpTable::PolyTable(vertices_num, j+2) ] ] ];
00240 }
00241 }
00242 indices += vertices_num;
00243
00244 }
00245 };
00246
00247 private:
00250 WALKER_TYPE *_walker;
00253 TRIMESH_TYPE *_mesh;
00256 bool _initialized;;
00259 bool _finalized;
00263 ScalarType _featureAngle;
00267 int _featureFlag;
00271 vcg::Point3i _corners[8];
00275 ScalarType _field[8];
00276
00277
00283 VertexPointer FindFeature(const std::vector<size_t> &vertices_idx)
00284 {
00285 unsigned int i, j, rank;
00286 size_t vertices_num = (size_t) vertices_idx.size();
00287
00288 CoordType *points = new CoordType[ vertices_num ];
00289 CoordType *normals = new CoordType[ vertices_num ];
00290
00291 for (i=0; i<vertices_num; i++)
00292 {
00293 points[i] = _mesh->vert[ vertices_idx[i] ].P();
00294 normals[i] = _mesh->vert[ vertices_idx[i] ].N();
00295 }
00296
00297
00298 CoordType center((ScalarType) 0.0, (ScalarType) 0.0, (ScalarType) 0.0);
00299 for (i=0; i<vertices_num; ++i)
00300 center += points[i];
00301 center /= (ScalarType) vertices_num;
00302 for (i=0; i<vertices_num; ++i)
00303 points[i] -= center;
00304
00305
00306 double c, minC, maxC;
00307 CoordType axis;
00308 for (minC=1.0, i=0; i<vertices_num-1; ++i)
00309 {
00310 for (j=i+1; j<vertices_num; ++j)
00311 {
00312 c = normals[i]*normals[j];
00313 if (c < minC)
00314 {
00315 minC = c;
00316 axis = normals[i] ^ normals[j];
00317 }
00318 }
00319 }
00320
00321 if (minC > cos(_featureAngle))
00322 return NULL;
00323
00324
00325 axis.Normalize();
00326 for (minC=1.0, maxC=-1.0, i=0; i<vertices_num; ++i)
00327 {
00328 c = axis * normals[i];
00329 if (c < minC) minC = c;
00330 if (c > maxC) maxC = c;
00331 }
00332 c = vcg::math::Max< double >(fabs(minC), fabs(maxC));
00333 c = sqrt(1.0-c*c);
00334 rank = (c > cos(_featureAngle) ? 2 : 3);
00335
00336
00337 vcg::ndim::Matrix<double> A((unsigned int) vertices_num, 3);
00338 double *b = new double[ vertices_num ];
00339 for (i=0; i<vertices_num; ++i)
00340 {
00341 A[i][0] = normals[i][0];
00342 A[i][1] = normals[i][1];
00343 A[i][2] = normals[i][2];
00344 b[i] = (points[i] * normals[i]);
00345 }
00346
00347
00348 vcg::ndim::Matrix<double> V(3, 3);
00349 double *w = new double[vertices_num];
00350 vcg::SingularValueDecomposition< typename vcg::ndim::Matrix<double> > (A, w, V, LeaveUnsorted, 100);
00351
00352
00353 if (rank == 2)
00354 {
00355 double smin = DBL_MAX;
00356 unsigned int sminid = 0;
00357 unsigned int srank = vcg::math::Min< unsigned int >(vertices_num, 3u);
00358
00359 for (i=0; i<srank; ++i)
00360 {
00361 if (w[i] < smin)
00362 {
00363 smin = w[i];
00364 sminid = i;
00365 }
00366 }
00367 w[sminid] = 0.0;
00368 }
00369
00370
00371 double *x = new double[3];
00372 vcg::SingularValueBacksubstitution< vcg::ndim::Matrix<double> >(A, w, V, x, b);
00373
00374
00375 CoordType point((ScalarType) x[0], (ScalarType) x[1], (ScalarType) x[2]);
00376 point += center;
00377
00378
00379 VertexPointer mean_point = &*AllocatorType::AddVertices( *_mesh, 1);
00380 mean_point->SetUserBit(_featureFlag);
00381 mean_point->P() = point;
00382 mean_point->N().SetZero();
00383 delete []x;
00384 delete []points;
00385 delete []normals;
00386 return mean_point;
00387 }
00388
00394 void FlipEdges()
00395 {
00396 size_t i;
00397 std::vector< LightEdge > edges;
00398 FaceIterator f_iter = _mesh->face.begin();
00399 FaceIterator f_end = _mesh->face.end();
00400 for (i=0; f_iter!=f_end; f_iter++, i++)
00401 {
00402 if (f_iter->V(1) > f_iter->V(0))
00403 {
00404 LightEdge le;
00405 le.face = i;
00406 le.edge = 0;
00407 edges.push_back( le );
00408 }
00409 if (f_iter->V(2) > f_iter->V(1))
00410 {
00411 LightEdge le;
00412 le.face = i;
00413 le.edge = 1;
00414 edges.push_back( LightEdge(le));
00415 }
00416 if (f_iter->V(0) > f_iter->V(2))
00417 {
00418 LightEdge le;
00419 le.face = i;
00420 le.edge = 2;
00421 edges.push_back( le );
00422 }
00423 }
00424 vcg::tri::UpdateTopology< TRIMESH_TYPE >::VertexFace( *_mesh );
00425 vcg::tri::UpdateTopology< TRIMESH_TYPE >::FaceFace( *_mesh );
00426
00427 typename std::vector< LightEdge >::iterator e_it = edges.begin();
00428 typename std::vector< LightEdge >::iterator e_end = edges.end();
00429
00430 FacePointer g, f;
00431 int w, z;
00432 for( ; e_it!=e_end; e_it++)
00433 {
00434 f = &_mesh->face[e_it->face];
00435 z = (int) e_it->edge;
00436
00437 if (vcg::face::CheckFlipEdge< FaceType >(*f, z))
00438 {
00439 VertexPointer v0, v1, v2, v3;
00440 v0 = f->V(z);
00441 v1 = f->V1(z);
00442 v2 = f->V2(z);
00443 g = f->FFp(z);
00444 w = f->FFi(z);
00445 v3 = g->V2(w);
00446 bool b0, b1, b2, b3;
00447 b0 = !v0->IsUserBit(_featureFlag);
00448 b1 = !v1->IsUserBit(_featureFlag);
00449 b2 = v2->IsUserBit(_featureFlag);
00450 b3 = v3->IsUserBit(_featureFlag);
00451 if( b0 && b1 && b2 && b3)
00452 vcg::face::FlipEdge< FaceType >(*f, z);
00453
00454 }
00455 }
00456 };
00457 };
00458
00459
00460
00461 }
00462 };
00463
00464 #endif // __VCG_EXTENDED_MARCHING_CUBES
