diff --git a/apps/sample/trimesh_voronoi/trimesh_voronoi.cpp b/apps/sample/trimesh_voronoi/trimesh_voronoi.cpp index eeb9037c..8b76f940 100644 --- a/apps/sample/trimesh_voronoi/trimesh_voronoi.cpp +++ b/apps/sample/trimesh_voronoi/trimesh_voronoi.cpp @@ -78,7 +78,7 @@ int main( int argc, char **argv ) float radius; tri::PoissonSampling(baseMesh,pointVec,sampleNum,radius); vector seedVec; - tri::VoronoiProcessing::PreprocessForVoronoi(baseMesh,radius); + tri::VoronoiProcessing::PreprocessForVoronoi(baseMesh,radius,vpp); tri::VoronoiProcessing::SeedToVertexConversion(baseMesh,pointVec,seedVec); int t2=clock(); @@ -93,10 +93,6 @@ int main( int argc, char **argv ) printf("relaxed %lu seeds for %i(up to %i) iterations in %f secs\n", seedVec.size(), actualIter, iterNum,float(t3-t2)/CLOCKS_PER_SEC); - MyMesh::PerVertexAttributeHandle sources; - sources= tri::Allocator::GetPerVertexAttribute (baseMesh,"sources"); - tri::Geodesic::Compute(baseMesh, seedVec, df,std::numeric_limits::max(),0,&sources); - tri::io::ExporterPLY::Save(baseMesh,"baseMesh.ply",tri::io::Mask::IOM_VERTCOLOR | tri::io::Mask::IOM_VERTQUALITY ); if(tri::VoronoiProcessing::CheckVoronoiTopology(baseMesh,seedVec)) { @@ -107,6 +103,7 @@ int main( int argc, char **argv ) printf("WARNING some voronoi region are not disk like; the resulting delaunay triangulation is not manifold.\n"); refineStep=1; } + tri::VoronoiProcessing::ConvertDelaunayTriangulationToMesh(baseMesh,delaunayMesh,seedVec,true); tri::io::ExporterPLY::Save(delaunayMesh,"delaunayBaseMesh.ply",tri::io::Mask::IOM_VERTCOLOR | tri::io::Mask::IOM_VERTFLAGS,false ); tri::VoronoiProcessing::RelaxRefineTriangulationSpring(baseMesh,delaunayMesh,refineStep,relaxStep); diff --git a/apps/sample/trimesh_voronoisampling/trimesh_voronoisampling.cpp b/apps/sample/trimesh_voronoisampling/trimesh_voronoisampling.cpp index ae1f688f..99334734 100644 --- a/apps/sample/trimesh_voronoisampling/trimesh_voronoisampling.cpp +++ b/apps/sample/trimesh_voronoisampling/trimesh_voronoisampling.cpp @@ -42,7 +42,7 @@ struct MyUsedTypes : public UsedTypes< Use ::AsVertexType, Use ::AsFaceType>{}; class MyVertex : public Vertex{}; -class MyFace : public Face< MyUsedTypes, face::VertexRef, face::Normal3f, face::BitFlags, face::VFAdj, face::FFAdj > {}; +class MyFace : public Face< MyUsedTypes, face::VertexRef, face::Normal3f, face::BitFlags, face::Mark, face::VFAdj, face::FFAdj > {}; class MyEdge : public Edge< MyUsedTypes, edge::VertexRef, edge::BitFlags>{}; class MyMesh : public tri::TriMesh< vector, vector, vector > {}; @@ -70,9 +70,9 @@ int main( int argc, char **argv ) int sampleNum = atoi(argv[2]); int iterNum = atoi(argv[3]); - bool fixCornerFlag=true; + bool fixCornerFlag=false; - printf("Reading %s and sampling %i points with %i iteration",argv[1],sampleNum,iterNum); + printf("Reading %s and sampling %i points with %i iteration\n",argv[1],sampleNum,iterNum); int ret= tri::io::ImporterPLY::Open(baseMesh,argv[1]); if(ret!=0) { @@ -85,11 +85,13 @@ int main( int argc, char **argv ) tri::Clean::RemoveUnreferencedVertex(baseMesh); tri::Allocator::CompactEveryVector(baseMesh); tri::UpdateTopology::VertexFace(baseMesh); - tri::UpdateFlags::FaceBorderFromVF(baseMesh); - tri::UpdateFlags::VertexBorderFromFace(baseMesh); tri::SurfaceSampling >::PoissonDiskParam pp; float radius = tri::SurfaceSampling >::ComputePoissonDiskRadius(baseMesh,sampleNum); + tri::VoronoiProcessing::PreprocessForVoronoi(baseMesh,radius,vpp); + + tri::UpdateFlags::FaceBorderFromVF(baseMesh); + tri::UpdateFlags::VertexBorderFromFace(baseMesh); // -- Build a sampling with just corners (Poisson filtered) @@ -153,26 +155,33 @@ int main( int argc, char **argv ) baseMesh.vert[i].SetS(); } - vpp.deleteUnreachedRegionFlag=true; - vpp.triangulateRegion = true; +// vpp.deleteUnreachedRegionFlag=true; + vpp.deleteUnreachedRegionFlag=false; + vpp.triangulateRegion = false; vpp.geodesicRelaxFlag=false; vpp.constrainSelectedSeed=true; tri::EuclideanDistance dd; int t0=clock(); // And now, at last, the relaxing procedure! - tri::VoronoiProcessing >::VoronoiRelaxing(baseMesh, seedVec, iterNum, dd, vpp); + int actualIter = tri::VoronoiProcessing >::VoronoiRelaxing(baseMesh, seedVec, iterNum, dd, vpp); int t1=clock(); MyMesh voroMesh, voroPoly, delaMesh; // Get the result in some pleasant form converting it to a real voronoi diagram. - tri::VoronoiProcessing >::ConvertVoronoiDiagramToMesh(baseMesh,voroMesh,voroPoly, seedVec, vpp); + if(tri::VoronoiProcessing::CheckVoronoiTopology(baseMesh,seedVec)) + tri::VoronoiProcessing::ConvertVoronoiDiagramToMesh(baseMesh,voroMesh,voroPoly,seedVec, vpp); + else + printf("WARNING some voronoi region are not disk like; the resulting delaunay triangulation is not manifold.\n"); + tri::io::ExporterPLY::Save(baseMesh,"base.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY ); tri::io::ExporterPLY::Save(voroMesh,"voroMesh.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_FLAGS ); tri::io::ExporterPLY::Save(voroPoly,"voroPoly.ply",tri::io::Mask::IOM_VERTCOLOR| tri::io::Mask::IOM_EDGEINDEX ,false); tri::VoronoiProcessing::ConvertDelaunayTriangulationToMesh(baseMesh,delaMesh, seedVec); tri::io::ExporterPLY::Save(delaMesh,"delaMesh.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY ); + tri::VoronoiProcessing::RelaxRefineTriangulationSpring(baseMesh,delaMesh,2,10); + tri::io::ExporterPLY::Save(delaMesh,"delaMeshRef.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY ); // tri::io::ImporterPLY::Open(baseMesh,argv[1]); @@ -186,6 +195,6 @@ int main( int argc, char **argv ) // tri::io::ExporterPLY::Save(outMesh,"outW.ply",tri::io::Mask::IOM_VERTCOLOR ); // tri::io::ExporterPLY::Save(polyMesh,"polyW.ply",tri::io::Mask::IOM_VERTCOLOR | tri::io::Mask::IOM_EDGEINDEX,false); // tri::io::ExporterDXF::Save(polyMesh,"outW.dxf"); - printf("Completed! %i iterations in %f sec for %lu seeds \n",iterNum,float(t1-t0)/CLOCKS_PER_SEC,seedVec.size()); + printf("Completed! %i (%i) iterations in %f sec for %lu seeds \n",actualIter, iterNum,float(t1-t0)/CLOCKS_PER_SEC,seedVec.size()); return 0; } diff --git a/vcg/complex/algorithms/voronoi_processing.h b/vcg/complex/algorithms/voronoi_processing.h index 0df1b6f6..0e418836 100644 --- a/vcg/complex/algorithms/voronoi_processing.h +++ b/vcg/complex/algorithms/voronoi_processing.h @@ -76,6 +76,7 @@ struct VoronoiProcessingParameter triangulateRegion=false; unbiasedSeedFlag = true; geodesicRelaxFlag = true; + refinementRatio = 5.0f; } int colorStrategy; @@ -93,6 +94,11 @@ struct VoronoiProcessingParameter bool preserveFixedSeed; /// If true the 'fixed' seeds are not moved during relaxation. /// \see FixVertexVector function to see how to fix a set of seeds. + float refinementRatio; /// It defines how much the input mesh has to be refined in order to have a supporting + /// triangulation that is dense enough to well approximate the voronoi diagram. + /// reasonable values are in the range 4..10. It is used by PreprocessForVoronoi and this value + /// says how many triangles you should expect in a voronoi region of a given radius. + // Convertion to Voronoi Diagram Parameters bool triangulateRegion; /// If true when building the voronoi diagram mesh each region is a @@ -121,6 +127,11 @@ class VoronoiProcessing public: + typedef typename MeshType::template PerVertexAttributeHandle PerVertexPointerHandle; + typedef typename MeshType::template PerVertexAttributeHandle PerVertexBoolHandle; + typedef typename MeshType::template PerFaceAttributeHandle PerFacePointerHandle; + + // Given a vector of point3f it finds the closest vertices on the mesh. static void SeedToVertexConversion(MeshType &m,std::vector &seedPVec,std::vector &seedVVec, bool compactFlag = true) @@ -152,8 +163,6 @@ static void SeedToVertexConversion(MeshType &m,std::vector &seedPVec, } } -typedef typename MeshType::template PerVertexAttributeHandle PerVertexPointerHandle; -typedef typename MeshType::template PerFaceAttributeHandle PerFacePointerHandle; static void ComputePerVertexSources(MeshType &m, std::vector &seedVec, DistanceFunctor &df) { @@ -418,8 +427,7 @@ static void ConvertVoronoiDiagramToMesh(MeshType &m, VoronoiProcessingParameter &vpp ) { tri::RequirePerVertexAttribute(m,"sources"); - typename MeshType::template PerVertexAttributeHandle sources; - sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); + PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); outMesh.Clear(); outPoly.Clear(); @@ -630,9 +638,10 @@ static void ConvertVoronoiDiagramToMesh(MeshType &m, // ******************* star to tri conversion ********* + // If requested the voronoi regions are converted from a star arragned polygon + // with vertex on the seed to a simple triangulated polygon by mean of a simple edge collapse if(vpp.triangulateRegion) { -// printf("Seedvec.size %i\n",seedVec.size()); for(FaceIterator fi=outMesh.face.begin();fi!=outMesh.face.end();++fi) if(!fi->IsD()) { for(int i=0;i<3;++i) @@ -642,14 +651,10 @@ static void ConvertVoronoiDiagramToMesh(MeshType &m, if( ((b0 && b1) || (fi->IsF(i) && !b0 && !b1) ) && tri::Index(outMesh,fi->V(i))V(i))]->IsS()) if(face::FFLinkCondition(*fi, i)) { - printf("collapse %i\n",tri::Index(outMesh,fi->V(i))); -// tri::io::ExporterPLY::Save(outMesh,"pre.ply"); face::FFEdgeCollapse(outMesh, *fi,i); -// tri::io::ExporterPLY::Save(outMesh,"post.ply"); break; } } @@ -889,10 +894,8 @@ static bool QuadricRelax(MeshType &m, std::vector &seedVec, std::v DistanceFunctor &df, VoronoiProcessingParameter &vpp) { newSeeds.clear(); - typename MeshType::template PerVertexAttributeHandle sources; - sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); - typename MeshType::template PerVertexAttributeHandle fixed; - fixed = tri::Allocator:: template GetPerVertexAttribute (m,"fixed"); + PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); + PerVertexBoolHandle fixed = tri::Allocator:: template GetPerVertexAttribute (m,"fixed"); QuadricSumDistance dz; std::vector dVec(m.vert.size(),dz); @@ -900,7 +903,14 @@ static bool QuadricRelax(MeshType &m, std::vector &seedVec, std::v { assert(sources[vi]!=0); int seedIndex = tri::Index(m,sources[vi]); - dVec[seedIndex].AddPoint(vi->P()); + // When constraining seeds movement we move selected seeds only onto other selected vertices + if(vpp.constrainSelectedSeed) + { // So we sum only the contribs of the selected vertices + if( (sources[vi]->IsS() && vi->IsS()) || (!sources[vi]->IsS())) + dVec[seedIndex].AddPoint(vi->P()); + } + else + dVec[seedIndex].AddPoint(vi->P()); } // Search the local maxima for each region and use them as new seeds @@ -1028,7 +1038,11 @@ static void PruneSeedByRegionArea(std::vector &seedVec, swap(seedVec,newSeedVec); } -/// Mark a vector of seeds to be fixed. +/// \brief Mark a vector of seeds to be fixed. +/// +/// Vertex pointers must belong to the mesh. +/// The framework use a boolean attribute called "fixed" to store this info. +/// static void FixVertexVector(MeshType &m, std::vector &vertToFixVec) { typename MeshType::template PerVertexAttributeHandle fixed; @@ -1059,8 +1073,8 @@ static int VoronoiRelaxing(MeshType &m, std::vector &seedVec, sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); typename MeshType::template PerVertexAttributeHandle fixed; fixed = tri::Allocator:: template GetPerVertexAttribute (m,"fixed"); - - for(int iter=0;iter &seedVec, newSeedVec[i]->C() = Color4b::White; swap(newSeedVec,seedVec); - if(!changed) relaxIter = iter; + if(!changed) break; } - return relaxIter; + + // Last run: Needed if we have changed the seed set to leave the sources handle correct. + if(iter==relaxIter) + tri::Geodesic::Compute(m, seedVec, df,std::numeric_limits::max(),0,&sources); + + return iter; } @@ -1206,6 +1225,15 @@ static bool CheckVoronoiTopology(MeshType& m,std::vector &seedVec) std::map seedMap; // It says if a given vertex of m is a seed (and its index in seedVec) BuildSeedMap(m,seedVec,seedMap); + // Very basic check: each vertex must have a source that is a seed. + for(int i=0;i regionVec(seedVec.size(),0); for(int i=0; i< seedVec.size();i++) regionVec[i] = new MeshType; @@ -1214,7 +1242,7 @@ static bool CheckVoronoiTopology(MeshType& m,std::vector &seedVec) int vi0 = seedMap[sources[m.face[i].V(0)]]; int vi1 = seedMap[sources[m.face[i].V(1)]]; int vi2 = seedMap[sources[m.face[i].V(2)]]; - + assert(vi0>=0 && vi1>=0 && vi2>=0); tri::Allocator::AddFace(*regionVec[vi0], m.face[i].cP(0),m.face[i].cP(1),m.face[i].cP(2)); if(vi1 != vi0) @@ -1281,6 +1309,8 @@ static void BuildSeedMap(MeshType &m, std::vector &seedVec, std:: seedMap[&(m.vert[i])]=-1; for(size_t i=0;i=0 && tri::Index(m,seedVec[i]) sources; - sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); + PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); outMesh.Clear(); tri::UpdateTopology::FaceFace(m); @@ -1364,7 +1393,8 @@ static void ConvertDelaunayTriangulationToMesh(MeshType &m, template static void PreprocessForVoronoi(MeshType &m, float radius, - MidPointType mid) + MidPointType mid, + VoronoiProcessingParameter &vpp) { const int maxSubDiv = 10; tri::RequireFFAdjacency(m); @@ -1374,17 +1404,17 @@ static void PreprocessForVoronoi(MeshType &m, float radius, for(int i=0;i(m,mid,min(edgeLen*2.0f,radius/5.0f)); + bool ret = tri::Refine(m,mid,min(edgeLen*2.0f,radius/vpp.refinementRatio)); if(!ret) break; } tri::Allocator::CompactEveryVector(m); tri::UpdateTopology::VertexFace(m); } -static void PreprocessForVoronoi(MeshType &m, float radius) +static void PreprocessForVoronoi(MeshType &m, float radius, VoronoiProcessingParameter &vpp) { - tri::MidPoint mid(m); - PreprocessForVoronoi >(m, radius,mid); + tri::MidPoint mid(&m); + PreprocessForVoronoi >(m, radius,mid,vpp); } static void RelaxRefineTriangulationSpring(MeshType &m, MeshType &delaMesh, int refineStep=3, int relaxStep=10 ) @@ -1399,10 +1429,17 @@ static void RelaxRefineTriangulationSpring(MeshType &m, MeshType &delaMesh, int const float eulerStep = 0.1f; tri::UpdateNormal::PerVertexNormalizedPerFaceNormalized(m); + typedef GridStaticPtr TriMeshGrid; TriMeshGrid ug; ug.Set(m.face.begin(),m.face.end()); + typedef typename vcg::SpatialHashTable HashVertexGrid; + HashVertexGrid HG; + HG.Set(m.vert.begin(),m.vert.end()); + + PerVertexBoolHandle fixed = tri::Allocator:: template GetPerVertexAttribute (m,"fixed"); + const ScalarType maxDist = m.bbox.Diag()/4.f; for(int kk=0;kk >(delaMesh,tri::MidPoint(&delaMesh)); } - tri::UpdateTopology::VertexFace(delaMesh); + tri::UpdateTopology::VertexFace(delaMesh); + const float dist_upper_bound=m.bbox.Diag()/10.0; + float dist; for(int k=0;k avgLenVec[i]*convergenceThr) changed = true; - delaMesh.vert[i].P()=closest; + VertexPointer vp = tri::GetClosestVertex(m, HG, delaMesh.vert[i].P(), dist_upper_bound, dist); + if(!fixed[vp] && !(vp->IsS())) // update only non fixed vertices + { + delaMesh.vert[i].P() += (avgForce[i]*eulerStep); + CoordType closest; + float dist; + tri::GetClosestFaceBase(m,ug,delaMesh.vert[i].cP(), maxDist,dist,closest); + assert(dist!=maxDist); + if(Distance(closest,delaMesh.vert[i].P()) > avgLenVec[i]*convergenceThr) changed = true; + delaMesh.vert[i].P()=closest; + } } + if(!changed) k=relaxStep; - } + } // end for k } }