manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

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
This commit is contained in:
Paolo Cignoni 2014-03-18 11:27:46 +00:00
parent 185c0f7152
commit c3f7b86500
3 changed files with 101 additions and 51 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
apps/sample/trimesh_voronoi/trimesh_voronoi.cpp
View File

@ -78,7 +78,7 @@ int main( int argc, char **argv )
float radius;
tri::PoissonSampling<MyMesh>(baseMesh,pointVec,sampleNum,radius);
vector<MyVertex *> seedVec;
tri::VoronoiProcessing<MyMesh>::PreprocessForVoronoi(baseMesh,radius);
tri::VoronoiProcessing<MyMesh>::PreprocessForVoronoi(baseMesh,radius,vpp);
tri::VoronoiProcessing<MyMesh>::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<MyVertex *> sources;
sources= tri::Allocator<MyMesh>::GetPerVertexAttribute<MyVertex *> (baseMesh,"sources");
tri::Geodesic<MyMesh>::Compute(baseMesh, seedVec, df,std::numeric_limits<float>::max(),0,&sources);
tri::io::ExporterPLY<MyMesh>::Save(baseMesh,"baseMesh.ply",tri::io::Mask::IOM_VERTCOLOR | tri::io::Mask::IOM_VERTQUALITY );
if(tri::VoronoiProcessing<MyMesh>::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<MyMesh>::ConvertDelaunayTriangulationToMesh(baseMesh,delaunayMesh,seedVec,true);
tri::io::ExporterPLY<MyMesh>::Save(delaunayMesh,"delaunayBaseMesh.ply",tri::io::Mask::IOM_VERTCOLOR | tri::io::Mask::IOM_VERTFLAGS,false );
tri::VoronoiProcessing<MyMesh>::RelaxRefineTriangulationSpring(baseMesh,delaunayMesh,refineStep,relaxStep);

29
apps/sample/trimesh_voronoisampling/trimesh_voronoisampling.cpp
View File

@ -42,7 +42,7 @@ struct MyUsedTypes : public UsedTypes< Use<MyVertex> ::AsVertexType,
Use<MyFace> ::AsFaceType>{};
class MyVertex : public Vertex<MyUsedTypes, vertex::Coord3f, vertex::Normal3f, vertex::VFAdj, vertex::Qualityf, vertex::Color4b, vertex::BitFlags >{};
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<MyVertex>, vector<MyEdge>, vector<MyFace> > {};
@ -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<MyMesh>::Open(baseMesh,argv[1]);
if(ret!=0)
{
@ -85,11 +85,13 @@ int main( int argc, char **argv )
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);
tri::SurfaceSampling<MyMesh,tri::TrivialSampler<MyMesh> >::PoissonDiskParam pp;
float radius = tri::SurfaceSampling<MyMesh,tri::TrivialSampler<MyMesh> >::ComputePoissonDiskRadius(baseMesh,sampleNum);
tri::VoronoiProcessing<MyMesh>::PreprocessForVoronoi(baseMesh,radius,vpp);
tri::UpdateFlags<MyMesh>::FaceBorderFromVF(baseMesh);
tri::UpdateFlags<MyMesh>::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<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 actualIter = tri::VoronoiProcessing<MyMesh, tri::EuclideanDistance<MyMesh> >::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<MyMesh, tri::EuclideanDistance<MyMesh> >::ConvertVoronoiDiagramToMesh(baseMesh,voroMesh,voroPoly, seedVec, vpp);
if(tri::VoronoiProcessing<MyMesh>::CheckVoronoiTopology(baseMesh,seedVec))
tri::VoronoiProcessing<MyMesh>::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<MyMesh>::Save(baseMesh,"base.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY );
tri::io::ExporterPLY<MyMesh>::Save(voroMesh,"voroMesh.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_FLAGS );
tri::io::ExporterPLY<MyMesh>::Save(voroPoly,"voroPoly.ply",tri::io::Mask::IOM_VERTCOLOR| tri::io::Mask::IOM_EDGEINDEX ,false);
tri::VoronoiProcessing<MyMesh>::ConvertDelaunayTriangulationToMesh(baseMesh,delaMesh, seedVec);
tri::io::ExporterPLY<MyMesh>::Save(delaMesh,"delaMesh.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY );
tri::VoronoiProcessing<MyMesh>::RelaxRefineTriangulationSpring(baseMesh,delaMesh,2,10);
tri::io::ExporterPLY<MyMesh>::Save(delaMesh,"delaMeshRef.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY );
// tri::io::ImporterPLY<MyMesh>::Open(baseMesh,argv[1]);
@ -186,6 +195,6 @@ int main( int argc, char **argv )
// tri::io::ExporterPLY<MyMesh>::Save(outMesh,"outW.ply",tri::io::Mask::IOM_VERTCOLOR );
// tri::io::ExporterPLY<MyMesh>::Save(polyMesh,"polyW.ply",tri::io::Mask::IOM_VERTCOLOR | tri::io::Mask::IOM_EDGEINDEX,false);
// tri::io::ExporterDXF<MyMesh>::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;
}

116
vcg/complex/algorithms/voronoi_processing.h
View File

@ -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<VertexPointer> PerVertexPointerHandle;
typedef typename MeshType::template PerVertexAttributeHandle<bool> PerVertexBoolHandle;
typedef typename MeshType::template PerFaceAttributeHandle<VertexPointer> PerFacePointerHandle;
// Given a vector of point3f it finds the closest vertices on the mesh.
static void SeedToVertexConversion(MeshType &m,std::vector<CoordType> &seedPVec,std::vector<VertexType *> &seedVVec, bool compactFlag = true)
@ -152,8 +163,6 @@ static void SeedToVertexConversion(MeshType &m,std::vector<CoordType> &seedPVec,
}
}
typedef typename MeshType::template PerVertexAttributeHandle<VertexPointer> PerVertexPointerHandle;
typedef typename MeshType::template PerFaceAttributeHandle<VertexPointer> PerFacePointerHandle;
static void ComputePerVertexSources(MeshType &m, std::vector<VertexType *> &seedVec, DistanceFunctor &df)
{
@ -418,8 +427,7 @@ static void ConvertVoronoiDiagramToMesh(MeshType &m,
VoronoiProcessingParameter &vpp )
{
tri::RequirePerVertexAttribute(m,"sources");
typename MeshType::template PerVertexAttributeHandle<VertexPointer> sources;
sources = tri::Allocator<MeshType>:: template GetPerVertexAttribute<VertexPointer> (m,"sources");
PerVertexPointerHandle sources = tri::Allocator<MeshType>:: template GetPerVertexAttribute<VertexPointer> (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))<seedVec.size())
{
// if(b0==b1)
if(!seedVec[tri::Index(outMesh,fi->V(i))]->IsS())
if(face::FFLinkCondition(*fi, i))
{
printf("collapse %i\n",tri::Index(outMesh,fi->V(i)));
// tri::io::ExporterPLY<MeshType>::Save(outMesh,"pre.ply");
face::FFEdgeCollapse(outMesh, *fi,i);
// tri::io::ExporterPLY<MeshType>::Save(outMesh,"post.ply");
break;
}
}
@ -889,10 +894,8 @@ static bool QuadricRelax(MeshType &m, std::vector<VertexType *> &seedVec, std::v
DistanceFunctor &df, VoronoiProcessingParameter &vpp)
{
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");
PerVertexPointerHandle sources = tri::Allocator<MeshType>:: template GetPerVertexAttribute<VertexPointer> (m,"sources");
PerVertexBoolHandle fixed = tri::Allocator<MeshType>:: template GetPerVertexAttribute<bool> (m,"fixed");
QuadricSumDistance dz;
std::vector<QuadricSumDistance> dVec(m.vert.size(),dz);
@ -900,7 +903,14 @@ static bool QuadricRelax(MeshType &m, std::vector<VertexType *> &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<VertexType *> &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<VertexType *> &vertToFixVec)
{
typename MeshType::template PerVertexAttributeHandle<bool> fixed;
@ -1059,8 +1073,8 @@ static int VoronoiRelaxing(MeshType &m, std::vector<VertexType *> &seedVec,
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)
int iter;
for(iter=0;iter<relaxIter;++iter)
{
if(cb) cb(iter*100/relaxIter,"Voronoi Lloyd Relaxation: First Partitioning");
@ -1103,9 +1117,14 @@ static int VoronoiRelaxing(MeshType &m, std::vector<VertexType *> &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<MeshType>::Compute(m, seedVec, df,std::numeric_limits<ScalarType>::max(),0,&sources);
return iter;
}
@ -1206,6 +1225,15 @@ static bool CheckVoronoiTopology(MeshType& m,std::vector<VertexType *> &seedVec)
std::map<VertexPointer, int> 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<m.vn;++i)
{
VertexPointer vp = sources[i];
int seedInd = seedMap[vp];
if(seedInd <0)
return false;
}
std::vector<MeshType *> 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<VertexType *> &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<MeshType>::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<VertexType *> &seedVec, std::
seedMap[&(m.vert[i])]=-1;
for(size_t i=0;i<seedVec.size();++i)
seedMap[seedVec[i]]=i;
for(size_t i=0;i<seedVec.size();++i)
assert(tri::Index(m,seedVec[i])>=0 && tri::Index(m,seedVec[i])<m.vn);
}
/// \brief Build a mesh of the Delaunay triangulation induced by the given seeds
@ -1301,8 +1331,7 @@ static void ConvertDelaunayTriangulationToMesh(MeshType &m,
tri::RequireCompactness(m);
tri::RequireVFAdjacency(m);
typename MeshType::template PerVertexAttributeHandle<VertexPointer> sources;
sources = tri::Allocator<MeshType>:: template GetPerVertexAttribute<VertexPointer> (m,"sources");
PerVertexPointerHandle sources = tri::Allocator<MeshType>:: template GetPerVertexAttribute<VertexPointer> (m,"sources");
outMesh.Clear();
tri::UpdateTopology<MeshType>::FaceFace(m);
@ -1364,7 +1393,8 @@ static void ConvertDelaunayTriangulationToMesh(MeshType &m,
template <class MidPointType >
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<maxSubDiv;++i)
{
bool ret = tri::Refine<MeshType, MidPointType >(m,mid,min(edgeLen*2.0f,radius/5.0f));
bool ret = tri::Refine<MeshType, MidPointType >(m,mid,min(edgeLen*2.0f,radius/vpp.refinementRatio));
if(!ret) break;
}
tri::Allocator<MeshType>::CompactEveryVector(m);
tri::UpdateTopology<MeshType>::VertexFace(m);
}
static void PreprocessForVoronoi(MeshType &m, float radius)
static void PreprocessForVoronoi(MeshType &m, float radius, VoronoiProcessingParameter &vpp)
{
tri::MidPoint<MeshType> mid(m);
PreprocessForVoronoi<tri::MidPoint<MeshType> >(m, radius,mid);
tri::MidPoint<MeshType> mid(&m);
PreprocessForVoronoi<tri::MidPoint<MeshType> >(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<MeshType>::PerVertexNormalizedPerFaceNormalized(m);
typedef GridStaticPtr<FaceType, ScalarType> TriMeshGrid;
TriMeshGrid ug;
ug.Set(m.face.begin(),m.face.end());
typedef typename vcg::SpatialHashTable<VertexType, ScalarType> HashVertexGrid;
HashVertexGrid HG;
HG.Set(m.vert.begin(),m.vert.end());
PerVertexBoolHandle fixed = tri::Allocator<MeshType>:: template GetPerVertexAttribute<bool> (m,"fixed");
const ScalarType maxDist = m.bbox.Diag()/4.f;
for(int kk=0;kk<refineStep;kk++)
{
@ -1414,7 +1451,9 @@ static void RelaxRefineTriangulationSpring(MeshType &m, MeshType &delaMesh, int
if(nonManif) return;
tri::Refine<MeshType, tri::MidPoint<MeshType> >(delaMesh,tri::MidPoint<MeshType>(&delaMesh));
}
tri::UpdateTopology<MeshType>::VertexFace(delaMesh);
tri::UpdateTopology<MeshType>::VertexFace(delaMesh);
const float dist_upper_bound=m.bbox.Diag()/10.0;
float dist;
for(int k=0;k<relaxStep;k++)
{
@ -1441,16 +1480,21 @@ static void RelaxRefineTriangulationSpring(MeshType &m, MeshType &delaMesh, int
bool changed=false;
for(int i=0;i<delaMesh.vn;++i)
{
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;
VertexPointer vp = tri::GetClosestVertex<MeshType,HashVertexGrid>(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
}
}