From c3f7b8650054b7cddf64d712cacd8d776cfd4941 Mon Sep 17 00:00:00 2001 From: cignoni Date: Tue, 18 Mar 2014 11:27:46 +0000 Subject: [PATCH] Some changes to the voronoi processing class. Now it performs Loyyd relaxation on constrained elements only keeping into account the constrained set. In other words sample on the boundary are relaxed only keeping into account of he other boundary vertexes This will result in much better distributions of samples on the boundaries. Improved also boundary management in the refinement/spring relaxing. Added a parameter for controlling the preprocessing refinment --- .../trimesh_voronoi/trimesh_voronoi.cpp | 7 +- .../trimesh_voronoisampling.cpp | 29 +++-- vcg/complex/algorithms/voronoi_processing.h | 116 ++++++++++++------ 3 files changed, 101 insertions(+), 51 deletions(-) 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 } }