Added the possibility of constraining the movement of some seeds onto a specific domain. Now during relaxation you can for example fix some seeds and constrain some other seeds over linear features defined as subset of vertices.

apps/sample/trimesh_voronoisampling/trimesh_voronoisampling.cpp

@@ -81,15 +81,47 @@ int main( int argc, char **argv )
     return -1;
   }
 
+  tri::Clean<MyMesh>::RemoveUnreferencedVertex(baseMesh);
+  tri::Allocator<MyMesh>::CompactEveryVector(baseMesh);
   tri::UpdateTopology<MyMesh>::VertexFace(baseMesh);
   tri::UpdateFlags<MyMesh>::FaceBorderFromVF(baseMesh);
+  tri::UpdateFlags<MyMesh>::VertexBorderFromFace(baseMesh);
 
-  // -- Build the mesh with corners
+  tri::SurfaceSampling<MyMesh,tri::TrivialSampler<MyMesh> >::PoissonDiskParam pp;
-  MyMesh cornerMesh;
+  float radius = tri::SurfaceSampling<MyMesh,tri::TrivialSampler<MyMesh> >::ComputePoissonDiskRadius(baseMesh,sampleNum);
+
+
+  // -- Build a sampling with just corners (Poisson filtered)
+  MyMesh poissonCornerMesh;
   std::vector<Point3f> sampleVec;
   tri::TrivialSampler<MyMesh> mps(sampleVec);
   tri::SurfaceSampling<MyMesh,tri::TrivialSampler<MyMesh> >::VertexBorderCorner(baseMesh,mps,math::ToRad(150.f));
-  tri::Build(cornerMesh,sampleVec);
+  tri::Build(poissonCornerMesh,sampleVec);
+  tri::io::ExporterPLY<MyMesh>::Save(poissonCornerMesh,"cornerMesh.ply");
+  sampleVec.clear();
+  tri::SurfaceSampling<MyMesh,tri::TrivialSampler<MyMesh> >::PoissonDiskPruning(mps, poissonCornerMesh, radius, pp);
+  tri::Build(poissonCornerMesh,sampleVec);
+  tri::io::ExporterPLY<MyMesh>::Save(poissonCornerMesh,"poissonCornerMesh.ply");
+
+  // Now save the corner as Fixed Seeds for later...
+  std::vector<MyVertex *> fixedSeedVec;
+  tri::VoronoiProcessing<MyMesh>::SeedToVertexConversion(baseMesh,sampleVec,fixedSeedVec);
+  tri::VoronoiProcessing<MyMesh, tri::EuclideanDistance<MyMesh> >::FixVertexVector(baseMesh,fixedSeedVec);
+
+  // -- Build a sampling with points on the border
+  MyMesh borderMesh,poissonBorderMesh;
+  sampleVec.clear();
+  tri::SurfaceSampling<MyMesh,tri::TrivialSampler<MyMesh> >::VertexBorder(baseMesh,mps);
+  tri::Build(borderMesh,sampleVec);
+  tri::io::ExporterPLY<MyMesh>::Save(borderMesh,"borderMesh.ply");
+
+  // -- and then prune the border sampling with poisson strategy using the precomputed corner vertexes.
+  pp.preGenMesh = &poissonCornerMesh;
+  pp.preGenFlag=true;
+  sampleVec.clear();
+  tri::SurfaceSampling<MyMesh,tri::TrivialSampler<MyMesh> >::PoissonDiskPruning(mps, borderMesh, radius, pp);
+  tri::Build(poissonBorderMesh,sampleVec);
+  tri::io::ExporterPLY<MyMesh>::Save(poissonBorderMesh,"PoissonEdgeMesh.ply");
 
   // -- Build the montercarlo sampling of the surface
   MyMesh MontecarloSurfaceMesh;
@@ -98,12 +130,9 @@ int main( int argc, char **argv )
   tri::Build(MontecarloSurfaceMesh,sampleVec);
   tri::io::ExporterPLY<MyMesh>::Save(MontecarloSurfaceMesh,"MontecarloSurfaceMesh.ply");
 
-  // -- Prune the montecarlo sampling with poisson strategy using the precomputed corner vertexes.
+  // -- Prune the montecarlo sampling with poisson strategy using the precomputed vertexes on the border.
-  tri::SurfaceSampling<MyMesh,tri::TrivialSampler<MyMesh> >::PoissonDiskParam pp;
+  pp.preGenMesh = &poissonBorderMesh;
-  pp.preGenMesh = &cornerMesh;
-  pp.preGenFlag=true;
   sampleVec.clear();
-  float radius = tri::SurfaceSampling<MyMesh,tri::TrivialSampler<MyMesh> >::ComputePoissonDiskRadius(baseMesh,sampleNum);
   tri::SurfaceSampling<MyMesh,tri::TrivialSampler<MyMesh> >::PoissonDiskPruning(mps, MontecarloSurfaceMesh, radius, pp);
   MyMesh PoissonMesh;
   tri::Build(PoissonMesh,sampleVec);
@@ -111,26 +140,32 @@ int main( int argc, char **argv )
 
   std::vector<MyVertex *> seedVec;
   tri::VoronoiProcessing<MyMesh>::SeedToVertexConversion(baseMesh,sampleVec,seedVec);
-  float eps = baseMesh.bbox.Diag()/10000.0f;
+
-  for(size_t i=0;i<cornerMesh.vert.size();++i)
+  // Select all the vertexes on the border to define a constrained domain.
-  {
+  // In our case we select the border vertexes to make sure that the seeds on the border
-    for(size_t j=0;j<seedVec.size();++j)
+  // relax themselves remaining on the border
-      if(Distance(cornerMesh.vert[i].P(),seedVec[j]->P()) < eps)
+  for(size_t i=0;i<baseMesh.vert.size();++i){
-        seedVec[j]->SetS();
+    if(baseMesh.vert[i].IsB())
+        baseMesh.vert[i].SetS();
   }
 
   tri::VoronoiProcessingParameter vpp;
   vpp.deleteUnreachedRegionFlag=true;
-  vpp.fixSelectedSeed=true;
+  vpp.preserveFixedSeed=true;
   vpp.collapseShortEdge=true;
   vpp.collapseShortEdgePerc=collapseShortEdgePerc;
   vpp.triangulateRegion = true;
   vpp.unbiasedSeedFlag =true;
+  vpp.geodesicRelaxFlag=false;
+  vpp.constrainSelectedSeed=true;
 
   tri::EuclideanDistance<MyMesh> dd;
   int t0=clock();
+  // And now, at last, the relaxing procedure!
   tri::VoronoiProcessing<MyMesh, tri::EuclideanDistance<MyMesh> >::VoronoiRelaxing(baseMesh, seedVec, iterNum, dd, vpp);
   int t1=clock();
+
+  // Get the result in some pleasant form converting it to a real voronoi diagram.
   tri::VoronoiProcessing<MyMesh, tri::EuclideanDistance<MyMesh> >::ConvertVoronoiDiagramToMesh(baseMesh,outMesh,polyMesh, seedVec, dd, vpp);
   tri::io::ExporterPLY<MyMesh>::Save(baseMesh,"base.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY );
   tri::io::ExporterPLY<MyMesh>::Save(outMesh,"out.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_FLAGS );

vcg/complex/algorithms/voronoi_clustering.h

@@ -65,7 +65,8 @@ struct VoronoiProcessingParameter
     colorStrategy = DistanceFromSeed;
     areaThresholdPerc=0;
     deleteUnreachedRegionFlag=false;
-    fixSelectedSeed=false;
+    constrainSelectedSeed=false;
+    preserveFixedSeed=false;
     collapseShortEdge=false;
     collapseShortEdgePerc = 0.01f;
     triangulateRegion=false;
@@ -76,9 +77,15 @@ struct VoronoiProcessingParameter
   float areaThresholdPerc;
   bool deleteUnreachedRegionFlag;
   bool unbiasedSeedFlag;
-  bool fixSelectedSeed;   /// the vertexes that are selected are used as fixed seeds:
+  bool constrainSelectedSeed;   /// If true the selected vertexes define a constraining domain:
-                          /// They will not move during relaxing
+                                /// During relaxation all selected seeds are constrained to move
-                          /// and they will be always included in the final triangulation (if on a border).
+                                /// only on other selected vertices.
+                                /// In this way you can constrain some seed to move only on certain
+                                /// domains, for example moving only along some linear features
+                                /// like border of creases.
+  bool preserveFixedSeed;       /// If true the 'fixed' seeds are not moved during relaxation.
+                                /// \see FixVertexVector function to see how to fix a set of seeds.
+
   bool triangulateRegion;
   bool collapseShortEdge;
   float collapseShortEdgePerc;
@@ -819,8 +826,11 @@ struct QuadricSumDistance
   }
 };
 
-/// Find the new position according to the geodesic rule.
+/// Find the new position
-/// For each region, given the frontiers, it chooses the point with the highest distance from the frontier
+/// For each region it search the vertex that minimize the sum of the squared distance
+/// from all the points of the region.
+/// It uses a vector of QuadricSumDistances
+/// (for simplicity of the size of the vertex but only the ones of the seed are used).
 ///
 static void QuadricRelax(MeshType &m, std::vector<VertexType *> &seedVec, std::vector<VertexPointer> &frontierVec,
                           std::vector<VertexType *> &newSeeds,
@@ -829,15 +839,18 @@ static void QuadricRelax(MeshType &m, std::vector<VertexType *> &seedVec, std::v
   newSeeds.clear();
   typename MeshType::template PerVertexAttributeHandle<VertexPointer> sources;
   sources = tri::Allocator<MeshType>:: template GetPerVertexAttribute<VertexPointer> (m,"sources");
+  typename MeshType::template PerVertexAttributeHandle<bool> fixed;
+  fixed = tri::Allocator<MeshType>:: template GetPerVertexAttribute<bool> (m,"fixed");
+
   QuadricSumDistance dz;
   std::vector<QuadricSumDistance> dVec(m.vert.size(),dz);
-  assert((int)m.vert.size()==m.vn);
   for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
   {
     assert(sources[vi]!=0);
     int seedIndex = tri::Index(m,sources[vi]);
     dVec[seedIndex].AddPoint(vi->P());
   }
+
   // Search the local maxima for each region and use them as new seeds
   std::pair<float,VertexPointer> zz(std::numeric_limits<ScalarType>::max(), static_cast<VertexPointer>(0));
   std::vector< std::pair<float,VertexPointer> > seedMaximaVec(m.vert.size(),zz);
@@ -847,22 +860,27 @@ static void QuadricRelax(MeshType &m, std::vector<VertexType *> &seedVec, std::v
     int seedIndex = tri::Index(m,sources[vi]);
     ScalarType val = dVec[seedIndex].Eval(vi->P());
     vi->Q()=val;
-    if(seedMaximaVec[seedIndex].first > val)
+    // if constrainSelectedSeed we search only among selected vertices
+    if(!vpp.constrainSelectedSeed || !sources[vi]->IsS() || vi->IsS())
     {
-      seedMaximaVec[seedIndex].first = val;
+      if(seedMaximaVec[seedIndex].first > val)
-      seedMaximaVec[seedIndex].second = &*vi;
+      {
+        seedMaximaVec[seedIndex].first = val;
+        seedMaximaVec[seedIndex].second = &*vi;
+      }
     }
   }
 
   tri::UpdateColor<MeshType>::PerVertexQualityRamp(m);
 //  tri::io::ExporterPLY<MeshType>::Save(m,"last.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY );
 
-  // update the seedvector with the new maxima (For the vertex not selected)
+  // update the seedvector with the new maxima (For the vertex not fixed)
   for(size_t i=0;i<m.vert.size();++i)
-    if(seedMaximaVec[i].second) // only  seeds entries have a non zero pointer
+    if(seedMaximaVec[i].second) // only updated entries have a non zero pointer
     {
-      if(vpp.fixSelectedSeed && sources[seedMaximaVec[i].second]->IsS())
+      VertexPointer curSrc = sources[seedMaximaVec[i].second];
-        newSeeds.push_back(sources[seedMaximaVec[i].second]);
+      if(vpp.preserveFixedSeed && fixed[curSrc])
+        newSeeds.push_back(curSrc);
       else
         newSeeds.push_back(seedMaximaVec[i].second);
     }
@@ -878,6 +896,8 @@ static void GeodesicRelax(MeshType &m, std::vector<VertexType *> &seedVec, std::
   newSeeds.clear();
   typename MeshType::template PerVertexAttributeHandle<VertexPointer> sources;
   sources = tri::Allocator<MeshType>:: template GetPerVertexAttribute<VertexPointer> (m,"sources");
+  typename MeshType::template PerVertexAttributeHandle<bool> fixed;
+  fixed = tri::Allocator<MeshType>:: template GetPerVertexAttribute<bool> (m,"fixed");
 
   std::vector<typename tri::Geodesic<MeshType>::VertDist> biasedFrontierVec;
   BuildBiasedSeedVec(m,df,seedVec,frontierVec,biasedFrontierVec,vpp);
@@ -895,21 +915,29 @@ static void GeodesicRelax(MeshType &m, std::vector<VertexType *> &seedVec, std::
   {
     assert(sources[vi]!=0);
     int seedIndex = tri::Index(m,sources[vi]);
-    if(seedMaximaVec[seedIndex].first < (*vi).Q())
+    if(!vpp.constrainSelectedSeed || !sources[vi]->IsS() || vi->IsS())
     {
-      seedMaximaVec[seedIndex].first=(*vi).Q();
+      if(seedMaximaVec[seedIndex].first < (*vi).Q())
-      seedMaximaVec[seedIndex].second=&*vi;
+      {
+        seedMaximaVec[seedIndex].first=(*vi).Q();
+        seedMaximaVec[seedIndex].second=&*vi;
+      }
     }
   }
 
   // update the seedvector with the new maxima (For the vertex not selected)
   for(size_t i=0;i<seedMaximaVec.size();++i)
-    if(seedMaximaVec[i].second)
+    if(seedMaximaVec[i].second)// only updated entries have a non zero pointer
     {
-      if(vpp.fixSelectedSeed && sources[seedMaximaVec[i].second]->IsS())
+      VertexPointer curSrc = sources[seedMaximaVec[i].second];
-        newSeeds.push_back(sources[seedMaximaVec[i].second]);
+      if(vpp.preserveFixedSeed && fixed[curSrc])
+        newSeeds.push_back(curSrc);
       else
         newSeeds.push_back(seedMaximaVec[i].second);
+//      if(vpp.fixSelectedSeed && sources[seedMaximaVec[i].second]->IsS())
+//        newSeeds.push_back(sources[seedMaximaVec[i].second]);
+//      else
+//        newSeeds.push_back(seedMaximaVec[i].second);
     }
 }
 
@@ -934,18 +962,37 @@ static void PruneSeedByRegionArea(std::vector<VertexType *> &seedVec,
   swap(seedVec,newSeedVec);
 }
 
+/// Mark a vector of seeds to be fixed.
+static void FixVertexVector(MeshType &m, std::vector<VertexType *> &vertToFixVec)
+{
+  typename MeshType::template PerVertexAttributeHandle<bool> fixed;
+  fixed = tri::Allocator<MeshType>:: template GetPerVertexAttribute<bool> (m,"fixed");
+  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
+    fixed[vi]=false;
+  for(size_t i=0;i<vertToFixVec.size();++i)
+    fixed[vertToFixVec[i]]=true;
+}
+
 /// \brief Perform a Lloyd relaxation cycle over a mesh
 ///
 ///
 
-static void VoronoiRelaxing(MeshType &m, std::vector<VertexType *> &seedVec, int relaxIter, DistanceFunctor &df,
+static void VoronoiRelaxing(MeshType &m, std::vector<VertexType *> &seedVec,
-                            VoronoiProcessingParameter &vpp, vcg::CallBackPos *cb=0)
+                            int relaxIter, DistanceFunctor &df,
+                            VoronoiProcessingParameter &vpp,
+                            vcg::CallBackPos *cb=0)
 {
   tri::RequireVFAdjacency(m);
+  tri::RequireCompactness(m);
+  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
+    assert(vi->VFp());
+
   tri::UpdateFlags<MeshType>::FaceBorderFromVF(m);
   tri::UpdateFlags<MeshType>::VertexBorderFromFace(m);
   typename MeshType::template PerVertexAttributeHandle<VertexPointer> sources;
   sources = tri::Allocator<MeshType>:: template GetPerVertexAttribute<VertexPointer> (m,"sources");
+  typename MeshType::template PerVertexAttributeHandle<bool> fixed;
+  fixed = tri::Allocator<MeshType>:: template GetPerVertexAttribute<bool> (m,"fixed");
 
   for(int iter=0;iter<relaxIter;++iter)
   {