From e94cfb5a43278ca206c690eeee98869b01550962 Mon Sep 17 00:00:00 2001 From: cignoni Date: Fri, 22 Mar 2013 17:06:41 +0000 Subject: [PATCH] Refactored a bit the extended marching cube core. Cleaned up a bit the trivial walker to be used by both of them. Updated the sample for marching cube. --- .../trimesh_isosurface/trimesh_isosurface.cpp | 16 +- .../algorithms/create/emc_lookup_table.h | 12 +- .../create/extended_marching_cubes.h | 103 ++++++----- .../algorithms/create/mc_trivial_walker.h | 175 +++++++++--------- 4 files changed, 156 insertions(+), 150 deletions(-) diff --git a/apps/sample/trimesh_isosurface/trimesh_isosurface.cpp b/apps/sample/trimesh_isosurface/trimesh_isosurface.cpp index f056ba13..d3b30465 100644 --- a/apps/sample/trimesh_isosurface/trimesh_isosurface.cpp +++ b/apps/sample/trimesh_isosurface/trimesh_isosurface.cpp @@ -36,9 +36,9 @@ class MyFace; class MyVertex; struct MyUsedTypes : public UsedTypes< Use ::AsVertexType, - Use ::AsFaceType>{}; + Use ::AsFaceType>{}; -class MyVertex : public Vertex< MyUsedTypes, vertex::Coord3f>{}; +class MyVertex : public Vertex< MyUsedTypes, vertex::Coord3f, vertex::Normal3f, vertex::BitFlags>{}; class MyFace : public Face< MyUsedTypes, face::VertexRef, face::BitFlags> {}; class MyMesh : public vcg::tri::TriMesh< std::vector< MyVertex>, std::vector< MyFace > > {}; @@ -49,15 +49,15 @@ typedef SimpleVolume MyVolume; int main(int /*argc*/ , char **/*argv*/) { - MyVolume volume; - + MyVolume volume; + typedef vcg::tri::TrivialWalker MyWalker; - typedef vcg::tri::MarchingCubes MyMarchingCubes; - MyWalker walker; - + typedef vcg::tri::MarchingCubes MyMarchingCubes; + MyWalker walker; + // Simple initialization of the volume with some cool perlin noise - volume.Init(Point3i(64,64,64)); + volume.Init(Point3i(64,64,64)); for(int i=0;i<64;i++) for(int j=0;j<64;j++) for(int k=0;k<64;k++) diff --git a/vcg/complex/algorithms/create/emc_lookup_table.h b/vcg/complex/algorithms/create/emc_lookup_table.h index fde806dd..962ad84f 100644 --- a/vcg/complex/algorithms/create/emc_lookup_table.h +++ b/vcg/complex/algorithms/create/emc_lookup_table.h @@ -4,7 +4,7 @@ * Copyright (C) 2002 by Computer Graphics Group, RWTH Aachen * * www.rwth-graphics.de * * * -*---------------------------------------------------------------------------* +*---------------------------------------------------------------------------* * * * License * * * @@ -23,7 +23,7 @@ * * \*===========================================================================*/ -//== TABLES ================================================================== +//== TABLES ================================================================== #ifndef __VCG_EMC_LOOK_UP_TABLE #define __VCG_EMC_LOOK_UP_TABLE @@ -35,7 +35,7 @@ namespace vcg class EMCLookUpTable { public: - static const int EdgeTable(unsigned char cubetype) + static int EdgeTable(unsigned char cubetype) { static const int edgeTable[256]= { @@ -70,11 +70,11 @@ namespace vcg 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190, 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, - 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 + 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 }; return edgeTable[cubetype]; }; // end of EdgeTable - + //----------------------------------------------------------------------------- static int* TriTable(unsigned char cubetype, int u) @@ -904,7 +904,7 @@ namespace vcg //----------------------------------------------------------------------------- - static const int PolyTable(unsigned int cubetype, int u) + static int PolyTable(unsigned int cubetype, int u) { static const int polyTable[8][16] = { diff --git a/vcg/complex/algorithms/create/extended_marching_cubes.h b/vcg/complex/algorithms/create/extended_marching_cubes.h index 755c1d1e..3270a1bc 100644 --- a/vcg/complex/algorithms/create/extended_marching_cubes.h +++ b/vcg/complex/algorithms/create/extended_marching_cubes.h @@ -348,13 +348,13 @@ namespace vcg //--A[i][2] = normals[i][2]; //--b[i] = (points[i] * normals[i]); A(i,0) = normals[i][0]; - A(i,0) = normals[i][1]; - A(i,0) = normals[i][2]; + A(i,1) = normals[i][1]; + A(i,2) = normals[i][2]; b(i) = (points[i] * normals[i]); } // SVD of matrix A - Eigen::JacobiSVD svd(A); + Eigen::JacobiSVD svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV); Eigen::MatrixXd sol(3,1); sol=svd.solve(b); @@ -386,7 +386,7 @@ namespace vcg // transform x to world coords //--CoordType point((ScalarType) x[0], (ScalarType) x[1], (ScalarType) x[2]); - CoordType point((ScalarType) sol[0], (ScalarType) sol[1], (ScalarType) sol[2]); + CoordType point((ScalarType) sol(0), (ScalarType) sol(1), (ScalarType) sol(2)); point += center; // Safety check if the feature point found by svd is @@ -411,58 +411,59 @@ namespace vcg */ void FlipEdges() { - size_t i; - std::vector< LightEdge > edges; - FaceIterator f_iter = _mesh->face.begin(); - FaceIterator f_end = _mesh->face.end(); - for (i=0; f_iter!=f_end; f_iter++, i++) + std::vector< LightEdge > edges; + for (FaceIterator fi = _mesh->face.begin(); fi!=_mesh->face.end(); fi++) + { + size_t i = tri::Index(*_mesh,*fi); + if (fi->V(1) > fi->V(0)) edges.push_back( LightEdge(i,0) ); + if (fi->V(2) > fi->V(1)) edges.push_back( LightEdge(i,1) ); + if (fi->V(0) > fi->V(2)) edges.push_back( LightEdge(i,2) ); + } + vcg::tri::UpdateTopology< TRIMESH_TYPE >::FaceFace( *_mesh ); + + // Select all the triangles that has a vertex shared with a non manifold edge. + int nonManifEdge = tri::Clean< TRIMESH_TYPE >::CountNonManifoldEdgeFF(*_mesh,true); + if(nonManifEdge >0) + tri::UpdateSelection< TRIMESH_TYPE >::FaceFromVertexLoose(*_mesh); + //qDebug("Got %i non manif edges",nonManifEdge); + + typename std::vector< LightEdge >::iterator e_it = edges.begin(); + typename std::vector< LightEdge >::iterator e_end = edges.end(); + + FacePointer g, f; + int w, z; + for( ; e_it!=e_end; e_it++) + { + f = &_mesh->face[e_it->face]; + z = (int) e_it->edge; + + // v2------v1 swap the diagonal only if v2 and v3 are feature and v0 and v1 are not. + // | / | + // | / | + // v0------v3 + if (!(f->IsS()) && vcg::face::CheckFlipEdge< FaceType >(*f, z)) { - if (f_iter->V(1) > f_iter->V(0)) edges.push_back( LightEdge(i,0) ); - if (f_iter->V(2) > f_iter->V(1)) edges.push_back( LightEdge(i,1) ); - if (f_iter->V(0) > f_iter->V(2)) edges.push_back( LightEdge(i,2) ); - } - vcg::tri::UpdateTopology< TRIMESH_TYPE >::FaceFace( *_mesh ); + VertexPointer v0, v1, v2, v3; + v0 = f->V(z); + v1 = f->V1(z); + v2 = f->V2(z); + g = f->FFp(z); + w = f->FFi(z); + v3 = g->V2(w); + bool b0, b1, b2, b3; + b0 = !v0->IsUserBit(_featureFlag) ; + b1 = !v1->IsUserBit(_featureFlag) ; + b2 = v2->IsUserBit(_featureFlag) ; + b3 = v3->IsUserBit(_featureFlag) ; + if( b0 && b1 && b2 && b3) + vcg::face::FlipEdge< FaceType >(*f, z); - // Select all the triangles that has a vertex shared with a non manifold edge. - int nonManifEdge = tri::Clean< TRIMESH_TYPE >::CountNonManifoldEdgeFF(*_mesh,true); - if(nonManifEdge >0) - tri::UpdateSelection< TRIMESH_TYPE >::FaceFromVertexLoose(*_mesh); - //qDebug("Got %i non manif edges",nonManifEdge); + } // end if (vcg::face::CheckFlipEdge< _Face >(*f, z)) + } // end for( ; e_it!=e_end; e_it++) - typename std::vector< LightEdge >::iterator e_it = edges.begin(); - typename std::vector< LightEdge >::iterator e_end = edges.end(); + } //end of FlipEdges - FacePointer g, f; - int w, z; - for( ; e_it!=e_end; e_it++) - { - f = &_mesh->face[e_it->face]; - z = (int) e_it->edge; -// v2------v1 swap the diagonal only if v2 and v3 are feature and v0 and v1 are not. -// | / | -// | / | -// v0------v3 - if (!(f->IsS()) && vcg::face::CheckFlipEdge< FaceType >(*f, z)) - { - VertexPointer v0, v1, v2, v3; - v0 = f->V(z); - v1 = f->V1(z); - v2 = f->V2(z); - g = f->FFp(z); - w = f->FFi(z); - v3 = g->V2(w); - bool b0, b1, b2, b3; - b0 = !v0->IsUserBit(_featureFlag) ; - b1 = !v1->IsUserBit(_featureFlag) ; - b2 = v2->IsUserBit(_featureFlag) ; - b3 = v3->IsUserBit(_featureFlag) ; - if( b0 && b1 && b2 && b3) - vcg::face::FlipEdge< FaceType >(*f, z); - - } // end if (vcg::face::CheckFlipEdge< _Face >(*f, z)) - } // end for( ; e_it!=e_end; e_it++) - }; //end of FlipEdges }; // end of class ExtendedMarchingCubes // /*! @} */ // end of Doxygen documentation diff --git a/vcg/complex/algorithms/create/mc_trivial_walker.h b/vcg/complex/algorithms/create/mc_trivial_walker.h index d721f6d5..775ff0a9 100644 --- a/vcg/complex/algorithms/create/mc_trivial_walker.h +++ b/vcg/complex/algorithms/create/mc_trivial_walker.h @@ -8,7 +8,7 @@ * \ * * All rights reserved. * * * -* This program is free software; you can redistribute it and/or modify * +* This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * @@ -26,7 +26,7 @@ namespace vcg { -// Very simple volume class. +// Very simple volume class. // just an example of the interface that the trivial walker expects template @@ -35,11 +35,11 @@ class SimpleVolume public: typedef VOX_TYPE VoxelType; std::vector Vol; - + Point3i sz; /// Dimensioni griglia come numero di celle per lato - - const Point3i &ISize() {return sz;}; /// Dimensioni griglia come numero di celle per lato - + + const Point3i &ISize() {return sz;} /// Dimensioni griglia come numero di celle per lato + void Init(Point3i _sz) { sz=_sz; @@ -47,39 +47,45 @@ public: } float Val(const int &x,const int &y,const int &z) const { - return cV(x,y,z).V(); - //else return numeric_limits::quiet_NaN( ); + return cV(x,y,z).V(); + //else return numeric_limits::quiet_NaN( ); } float &Val(const int &x,const int &y,const int &z) { - return V(x,y,z).V(); - //else return numeric_limits::quiet_NaN( ); + return V(x,y,z).V(); + //else return numeric_limits::quiet_NaN( ); } VOX_TYPE &V(const int &x,const int &y,const int &z) { - return Vol[x+y*sz[0]+z*sz[0]*sz[1]]; + return Vol[x+y*sz[0]+z*sz[0]*sz[1]]; + } + + VOX_TYPE &V(const Point3i &pi) { + return Vol[ pi[0] + pi[1]*sz[0] + pi[2]*sz[0]*sz[1] ]; } const VOX_TYPE &cV(const int &x,const int &y,const int &z) const { - return Vol[x+y*sz[0]+z*sz[0]*sz[1]]; + return Vol[x+y*sz[0]+z*sz[0]*sz[1]]; } -typedef enum { XAxis=0,YAxis=1,ZAxis=2} VolumeAxis; + typedef enum { XAxis=0,YAxis=1,ZAxis=2} VolumeAxis; -template < class VertexPointerType, VolumeAxis AxisVal > - void GetIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr) -{ - float f1 = Val(p1.X(), p1.Y(), p1.Z())-thr; - float f2 = Val(p2.X(), p2.Y(), p2.Z())-thr; - float u = (float) f1/(f1-f2); - if(AxisVal==XAxis) v->P().X() = (float) p1.X()*(1-u) + u*p2.X(); - else v->P().X() = (float) p1.X(); - if(AxisVal==YAxis) v->P().Y() = (float) p1.Y()*(1-u) + u*p2.Y(); - else v->P().Y() = (float) p1.Y(); - if(AxisVal==ZAxis) v->P().Z() = (float) p1.Z()*(1-u) + u*p2.Z(); - else v->P().Z() = (float) p1.Z(); -} + template < class VertexPointerType, VolumeAxis AxisVal > + void GetIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr) + { + float f1 = V(p1).V()-thr; + float f2 = V(p2).V()-thr; + float u = (float) f1/(f1-f2); + if(AxisVal==XAxis) v->P().X() = (float) p1.X()*(1-u) + u*p2.X(); + else v->P().X() = (float) p1.X(); + if(AxisVal==YAxis) v->P().Y() = (float) p1.Y()*(1-u) + u*p2.Y(); + else v->P().Y() = (float) p1.Y(); + if(AxisVal==ZAxis) v->P().Z() = (float) p1.Z()*(1-u) + u*p2.Z(); + else v->P().Z() = (float) p1.Z(); + + if(VoxelType::HasNormal()) v->N() = V(p1).N()*(1-u) + V(p2).N()*u; + } template < class VertexPointerType > void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr) @@ -93,32 +99,31 @@ template < class VertexPointerType > void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr) { GetIntercept(p1,p2,v,thr); } }; -template -class RawVolumeImporter -{ -public: - enum DataType -{ - // Funzioni superiori - UNDEF=0, - BYTE=1, - SHORT=2, - FLOAT=3 -}; - -static bool Open(const char *filename, VolumeType &V, Point3i sz, DataType d) -{ -return true; -} -}; class SimpleVoxel { private: float _v; public: - float &V() {return _v;}; - float V() const {return _v;}; + float &V() {return _v;} + float V() const {return _v;} + static bool HasNormal() {return false;} + vcg::Point3f N() const {return Point3f(0,0,0);} + vcg::Point3f &N() { static Point3f _p(0,0,0); return _p;} +}; + +class SimpleVoxelWithNormal +{ +private: + float _v; + vcg::Point3f _n; +public: + float &V() {return _v;} + float V() const {return _v;} + vcg::Point3f &N() {return _n;} + vcg::Point3f N() const {return _n;} + static bool HasNormal() {return true;} + }; @@ -126,21 +131,21 @@ namespace tri { // La classe Walker implementa la politica di visita del volume; conoscendo l'ordine di visita del volume -// Ë conveniente che il Walker stesso si faccia carico del caching dei dati utilizzati durante l'esecuzione +// Ë conveniente che il Walker stesso si faccia carico del caching dei dati utilizzati durante l'esecuzione // degli algoritmi MarchingCubes ed ExtendedMarchingCubes, in particolare il calcolo del volume ai vertici // delle celle e delle intersezioni della superficie con le celle. In questo esempio il volume da processare // viene suddiviso in fette; in questo modo se il volume ha dimensione h*l*w (rispettivamente altezza, // larghezza e profondit‡), lo spazio richiesto per il caching dei vertici gi‡ allocati passa da O(h*l*w) -// a O(h*l). +// a O(h*l). template class TrivialWalker { private: - typedef int VertexIndex; + typedef int VertexIndex; typedef typename MeshType::ScalarType ScalarType; typedef typename MeshType::VertexPointer VertexPointer; - public: + public: // bbox is the portion of the volume to be computed // resolution determine the sampling step: @@ -148,37 +153,37 @@ private: void Init(VolumeType &volume) - { - _bbox = Box3i(Point3i(0,0,0),volume.ISize()); - _slice_dimension = _bbox.DimX()*_bbox.DimZ(); + { + _bbox = Box3i(Point3i(0,0,0),volume.ISize()); + _slice_dimension = _bbox.DimX()*_bbox.DimZ(); _x_cs = new VertexIndex[ _slice_dimension ]; _y_cs = new VertexIndex[ _slice_dimension ]; _z_cs = new VertexIndex[ _slice_dimension ]; _x_ns = new VertexIndex[ _slice_dimension ]; _z_ns = new VertexIndex[ _slice_dimension ]; - + }; ~TrivialWalker() {_thr=0;} - - template + + template void BuildMesh(MeshType &mesh, VolumeType &volume, EXTRACTOR_TYPE &extractor, const float threshold, vcg::CallBackPos * cb=0) - { + { Init(volume); - _volume = &volume; - _mesh = &mesh; - _mesh->Clear(); + _volume = &volume; + _mesh = &mesh; + _mesh->Clear(); _thr=threshold; - vcg::Point3i p1, p2; + vcg::Point3i p1, p2; Begin(); extractor.Initialize(); for (int j=_bbox.min.Y(); j<(_bbox.max.Y()-1)-1; j+=1) - { + { - if(cb && ((j%10)==0) ) cb(j*_bbox.DimY()/100.0,"Marching volume"); + if(cb && ((j%10)==0) ) cb(j*_bbox.DimY()/100.0,"Marching volume"); for (int i=_bbox.min.X(); i<(_bbox.max.X()-1)-1; i+=1) { @@ -186,7 +191,7 @@ private: { p1.X()=i; p1.Y()=j; p1.Z()=k; p2.X()=i+1; p2.Y()=j+1; p2.Z()=k+1; - extractor.ProcessCell(p1, p2); + extractor.ProcessCell(p1, p2); } } NextSlice(); @@ -198,11 +203,11 @@ private: float V(int pi, int pj, int pk) { - return _volume->Val(pi, pj, pk)-_thr; + return _volume->Val(pi, pj, pk)-_thr; } bool Exist(const vcg::Point3i &p0, const vcg::Point3i &p1, VertexPointer &v) - { + { int pos = p0.X()+p0.Z()*_bbox.max.X(); int vidx; @@ -219,8 +224,8 @@ private: return v!=NULL; } - void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) - { + void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) + { int i = p1.X() - _bbox.min.X(); int z = p1.Z() - _bbox.min.Z(); VertexIndex index = i+z*_bbox.max.X(); @@ -231,7 +236,7 @@ private: { _x_cs[index] = (VertexIndex) _mesh->vert.size(); pos = _x_cs[index]; - Allocator::AddVertices( *_mesh, 1 ); + Allocator::AddVertices( *_mesh, 1 ); v = &_mesh->vert[pos]; _volume->GetXIntercept(p1, p2, v, _thr); return; @@ -249,10 +254,10 @@ private: return; } } - assert(pos >=0 && size_t(pos)< _mesh->vert.size()); + assert(pos >=0 && size_t(pos)< _mesh->vert.size()); v = &_mesh->vert[pos]; } - void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) + void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) { int i = p1.X() - _bbox.min.X(); int z = p1.Z() - _bbox.min.Z(); @@ -268,7 +273,7 @@ private: } v = &_mesh->vert[pos]; } - void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) + void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) { int i = p1.X() - _bbox.min.X(); int z = p1.Z() - _bbox.min.Z(); @@ -306,26 +311,26 @@ protected: int _slice_dimension; int _current_slice; - + VertexIndex *_x_cs; // indici dell'intersezioni della superficie lungo gli Xedge della fetta corrente VertexIndex *_y_cs; // indici dell'intersezioni della superficie lungo gli Yedge della fetta corrente VertexIndex *_z_cs; // indici dell'intersezioni della superficie lungo gli Zedge della fetta corrente - VertexIndex *_x_ns; // indici dell'intersezioni della superficie lungo gli Xedge della prossima fetta - VertexIndex *_z_ns; // indici dell'intersezioni della superficie lungo gli Zedge della prossima fetta + VertexIndex *_x_ns; // indici dell'intersezioni della superficie lungo gli Xedge della prossima fetta + VertexIndex *_z_ns; // indici dell'intersezioni della superficie lungo gli Zedge della prossima fetta MeshType *_mesh; VolumeType *_volume; float _thr; - void NextSlice() - { - memset(_x_cs, -1, _slice_dimension*sizeof(VertexIndex)); - memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex)); - memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex)); + void NextSlice() + { + memset(_x_cs, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex)); std::swap(_x_cs, _x_ns); - std::swap(_z_cs, _z_ns); - + std::swap(_z_cs, _z_ns); + _current_slice += 1; } @@ -338,9 +343,9 @@ protected: memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex)); memset(_x_ns, -1, _slice_dimension*sizeof(VertexIndex)); memset(_z_ns, -1, _slice_dimension*sizeof(VertexIndex)); - + } }; -} // end namespace -} // end namespace +} // end namespace +} // end namespace #endif // __VCGTEST_WALKER