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

Many Changes to point sampling: 

- Heavy refactoring (typename changes, declaration etc) for clarity and shortness
- add a reset() method to the sampler interface
- changed the way used to get a weight in the sampling process now the variance is used to map the 'metric' along 1 and <variance>; Made it uniform in poisson and montecarlo sampling.
- changed the way in which the weight is used/passed: now with attribute!
- added exact number poisson disk pruning.
- stats are always computed (no performance impact clearer code)
This commit is contained in:
Paolo Cignoni 2014-05-15 10:35:08 +00:00
parent 4fc8ed68ba
commit 3b7753ef20
1 changed files with 315 additions and 195 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

462
vcg/complex/algorithms/point_sampling.h
View File

@ -39,6 +39,7 @@ sampling strategies (montecarlo, stratified etc).
#include <vcg/math/random_generator.h> #include <vcg/math/random_generator.h>
#include <vcg/complex/algorithms/closest.h> #include <vcg/complex/algorithms/closest.h>
#include <vcg/space/index/spatial_hashing.h> #include <vcg/space/index/spatial_hashing.h>
#include <vcg/complex/algorithms/hole.h>
#include <vcg/complex/algorithms/stat.h> #include <vcg/complex/algorithms/stat.h>
#include <vcg/complex/algorithms/create/platonic.h> #include <vcg/complex/algorithms/create/platonic.h>
#include <vcg/complex/algorithms/update/topology.h> #include <vcg/complex/algorithms/update/topology.h>
@ -75,6 +76,11 @@ public:
typedef typename MeshType::VertexType VertexType; typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::FaceType FaceType; typedef typename MeshType::FaceType FaceType;
void reset()
{
sampleVec->clear();
}
TrivialSampler() TrivialSampler()
{ {
sampleVec = new std::vector<CoordType>(); sampleVec = new std::vector<CoordType>();
@ -84,8 +90,8 @@ public:
TrivialSampler(std::vector<CoordType> &Vec) TrivialSampler(std::vector<CoordType> &Vec)
{ {
sampleVec = &Vec; sampleVec = &Vec;
sampleVec->clear();
vectorOwner=false; vectorOwner=false;
reset();
} }
~TrivialSampler() ~TrivialSampler()
@ -126,6 +132,11 @@ public:
TrivialPointerSampler() {} TrivialPointerSampler() {}
~TrivialPointerSampler() {} ~TrivialPointerSampler() {}
void reset()
{
sampleVec.clear();
}
public: public:
std::vector<VertexType *> sampleVec; std::vector<VertexType *> sampleVec;
@ -149,6 +160,10 @@ public:
MeshSampler(MeshType &_m):m(_m){} MeshSampler(MeshType &_m):m(_m){}
MeshType &m; MeshType &m;
void reset()
{
m.Clear();
}
void AddVert(const VertexType &p) void AddVert(const VertexType &p)
{ {
@ -180,18 +195,18 @@ The class is templated over the PointSampler object that allows to customize the
**/ **/
template <class MetroMesh, class VertexSampler = TrivialSampler< MetroMesh> > template <class MeshType, class VertexSampler = TrivialSampler< MeshType> >
class SurfaceSampling class SurfaceSampling
{ {
typedef typename MetroMesh::CoordType CoordType; typedef typename MeshType::CoordType CoordType;
typedef typename MetroMesh::ScalarType ScalarType; typedef typename MeshType::ScalarType ScalarType;
typedef typename MetroMesh::VertexType VertexType; typedef typename MeshType::VertexType VertexType;
typedef typename MetroMesh::VertexPointer VertexPointer; typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MetroMesh::VertexIterator VertexIterator; typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MetroMesh::FacePointer FacePointer; typedef typename MeshType::FacePointer FacePointer;
typedef typename MetroMesh::FaceIterator FaceIterator; typedef typename MeshType::FaceIterator FaceIterator;
typedef typename MetroMesh::FaceType FaceType; typedef typename MeshType::FaceType FaceType;
typedef typename MetroMesh::FaceContainer FaceContainer; typedef typename MeshType::FaceContainer FaceContainer;
typedef typename vcg::Box3<ScalarType> BoxType; typedef typename vcg::Box3<ScalarType> BoxType;
typedef typename vcg::SpatialHashTable<FaceType, ScalarType> MeshSHT; typedef typename vcg::SpatialHashTable<FaceType, ScalarType> MeshSHT;
@ -201,6 +216,8 @@ class SurfaceSampling
typedef typename vcg::SpatialHashTable<VertexType, ScalarType> SampleSHT; typedef typename vcg::SpatialHashTable<VertexType, ScalarType> SampleSHT;
typedef typename vcg::SpatialHashTable<VertexType, ScalarType>::CellIterator SampleSHTIterator; typedef typename vcg::SpatialHashTable<VertexType, ScalarType>::CellIterator SampleSHTIterator;
typedef typename MeshType::template PerVertexAttributeHandle<float> PerVertexFloatAttribute;
public: public:
static math::MarsenneTwisterRNG &SamplingRandomGenerator() static math::MarsenneTwisterRNG &SamplingRandomGenerator()
@ -338,7 +355,7 @@ static int Poisson(double lambda)
} }
static void AllVertex(MetroMesh & m, VertexSampler &ps) static void AllVertex(MeshType & m, VertexSampler &ps)
{ {
VertexIterator vi; VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi) for(vi=m.vert.begin();vi!=m.vert.end();++vi)
@ -358,7 +375,7 @@ static void AllVertex(MetroMesh & m, VertexSampler &ps)
/// 2) shuffle vertices. /// 2) shuffle vertices.
/// 3) for each vertices choose it if rand > thr; /// 3) for each vertices choose it if rand > thr;
static void VertexWeighted(MetroMesh & m, VertexSampler &ps, int sampleNum) static void VertexWeighted(MeshType & m, VertexSampler &ps, int sampleNum)
{ {
ScalarType qSum = 0; ScalarType qSum = 0;
VertexIterator vi; VertexIterator vi;
@ -396,7 +413,7 @@ static void VertexWeighted(MetroMesh & m, VertexSampler &ps, int sampleNum)
/// Sample the vertices in a uniform way. Each vertex has a probability of being chosen /// Sample the vertices in a uniform way. Each vertex has a probability of being chosen
/// that is proportional to the area it represent. /// that is proportional to the area it represent.
static void VertexAreaUniform(MetroMesh & m, VertexSampler &ps, int sampleNum) static void VertexAreaUniform(MeshType & m, VertexSampler &ps, int sampleNum)
{ {
VertexIterator vi; VertexIterator vi;
for(vi = m.vert.begin(); vi != m.vert.end(); ++vi) for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
@ -416,10 +433,9 @@ static void VertexAreaUniform(MetroMesh & m, VertexSampler &ps, int sampleNum)
VertexWeighted(m,ps,sampleNum); VertexWeighted(m,ps,sampleNum);
} }
static void FillAndShuffleFacePointerVector(MetroMesh & m, std::vector<FacePointer> &faceVec) static void FillAndShuffleFacePointerVector(MeshType & m, std::vector<FacePointer> &faceVec)
{ {
FaceIterator fi; for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)
for(fi=m.face.begin();fi!=m.face.end();++fi)
if(!(*fi).IsD()) faceVec.push_back(&*fi); if(!(*fi).IsD()) faceVec.push_back(&*fi);
assert((int)faceVec.size()==m.fn); assert((int)faceVec.size()==m.fn);
@ -427,10 +443,9 @@ static void FillAndShuffleFacePointerVector(MetroMesh & m, std::vector<FacePoint
unsigned int (*p_myrandom)(unsigned int) = RandomInt; unsigned int (*p_myrandom)(unsigned int) = RandomInt;
std::random_shuffle(faceVec.begin(),faceVec.end(), p_myrandom); std::random_shuffle(faceVec.begin(),faceVec.end(), p_myrandom);
} }
static void FillAndShuffleVertexPointerVector(MetroMesh & m, std::vector<VertexPointer> &vertVec) static void FillAndShuffleVertexPointerVector(MeshType & m, std::vector<VertexPointer> &vertVec)
{ {
VertexIterator vi; for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
if(!(*vi).IsD()) vertVec.push_back(&*vi); if(!(*vi).IsD()) vertVec.push_back(&*vi);
assert((int)vertVec.size()==m.vn); assert((int)vertVec.size()==m.vn);
@ -440,7 +455,7 @@ static void FillAndShuffleVertexPointerVector(MetroMesh & m, std::vector<VertexP
} }
/// Sample the vertices in a uniform way. Each vertex has the same probabiltiy of being chosen. /// Sample the vertices in a uniform way. Each vertex has the same probabiltiy of being chosen.
static void VertexUniform(MetroMesh & m, VertexSampler &ps, int sampleNum) static void VertexUniform(MeshType & m, VertexSampler &ps, int sampleNum)
{ {
if(sampleNum>=m.vn) { if(sampleNum>=m.vn) {
AllVertex(m,ps); AllVertex(m,ps);
@ -459,9 +474,9 @@ static void VertexUniform(MetroMesh & m, VertexSampler &ps, int sampleNum)
/// It assumes that the border flag have been set over the mesh both for vertex and for faces. /// It assumes that the border flag have been set over the mesh both for vertex and for faces.
/// All the vertices on the border where the surface forms an angle smaller than the given threshold are sampled. /// All the vertices on the border where the surface forms an angle smaller than the given threshold are sampled.
/// ///
static void VertexBorderCorner(MetroMesh & m, VertexSampler &ps, float angleRad) static void VertexBorderCorner(MeshType & m, VertexSampler &ps, float angleRad)
{ {
typename MetroMesh::template PerVertexAttributeHandle <float> angleSumH = tri::Allocator<MetroMesh>:: template GetPerVertexAttribute<float> (m); typename MeshType::template PerVertexAttributeHandle <float> angleSumH = tri::Allocator<MeshType>:: template GetPerVertexAttribute<float> (m);
for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
angleSumH[vi]=0; angleSumH[vi]=0;
@ -480,7 +495,7 @@ static void VertexBorderCorner(MetroMesh & m, VertexSampler &ps, float angleRad)
ps.AddVert(*vi); ps.AddVert(*vi);
} }
tri::Allocator<MetroMesh>:: template DeletePerVertexAttribute<float> (m,angleSumH); tri::Allocator<MeshType>:: template DeletePerVertexAttribute<float> (m,angleSumH);
} }
/// \brief Sample all the border vertices /// \brief Sample all the border vertices
@ -488,7 +503,7 @@ static void VertexBorderCorner(MetroMesh & m, VertexSampler &ps, float angleRad)
/// It assumes that the border flag have been set over the mesh. /// It assumes that the border flag have been set over the mesh.
/// All the vertices on the border are sampled. /// All the vertices on the border are sampled.
/// ///
static void VertexBorder(MetroMesh & m, VertexSampler &ps) static void VertexBorder(MeshType & m, VertexSampler &ps)
{ {
VertexBorderCorner(m,ps,std::numeric_limits<ScalarType>::max()); VertexBorderCorner(m,ps,std::numeric_limits<ScalarType>::max());
} }
@ -499,14 +514,14 @@ static void VertexBorder(MetroMesh & m, VertexSampler &ps)
/// tri::UpdateFlags<MeshType>::FaceFauxCrease(mesh,creaseAngleRad); /// tri::UpdateFlags<MeshType>::FaceFauxCrease(mesh,creaseAngleRad);
/// Then it chooses all the vertices where there are at least three non faux edges. /// Then it chooses all the vertices where there are at least three non faux edges.
/// ///
static void VertexCrease(MetroMesh & m, VertexSampler &ps) static void VertexCrease(MeshType & m, VertexSampler &ps)
{ {
typedef typename UpdateTopology<MetroMesh>::PEdge SimpleEdge; typedef typename UpdateTopology<MeshType>::PEdge SimpleEdge;
std::vector< SimpleEdge > Edges; std::vector< SimpleEdge > Edges;
typename std::vector< SimpleEdge >::iterator ei; typename std::vector< SimpleEdge >::iterator ei;
UpdateTopology<MetroMesh>::FillUniqueEdgeVector(m,Edges,false); UpdateTopology<MeshType>::FillUniqueEdgeVector(m,Edges,false);
typename MetroMesh::template PerVertexAttributeHandle <int> hv = tri::Allocator<MetroMesh>:: template GetPerVertexAttribute<int> (m); typename MeshType::template PerVertexAttributeHandle <int> hv = tri::Allocator<MeshType>:: template GetPerVertexAttribute<int> (m);
for(ei=Edges.begin(); ei!=Edges.end(); ++ei) for(ei=Edges.begin(); ei!=Edges.end(); ++ei)
{ {
@ -522,7 +537,7 @@ static void VertexCrease(MetroMesh & m, VertexSampler &ps)
} }
static void FaceUniform(MetroMesh & m, VertexSampler &ps, int sampleNum) static void FaceUniform(MeshType & m, VertexSampler &ps, int sampleNum)
{ {
if(sampleNum>=m.fn) { if(sampleNum>=m.fn) {
AllFace(m,ps); AllFace(m,ps);
@ -536,7 +551,7 @@ static void FaceUniform(MetroMesh & m, VertexSampler &ps, int sampleNum)
ps.AddFace(*faceVec[i],Barycenter(*faceVec[i])); ps.AddFace(*faceVec[i],Barycenter(*faceVec[i]));
} }
static void AllFace(MetroMesh & m, VertexSampler &ps) static void AllFace(MeshType & m, VertexSampler &ps)
{ {
FaceIterator fi; FaceIterator fi;
for(fi=m.face.begin();fi!=m.face.end();++fi) for(fi=m.face.begin();fi!=m.face.end();++fi)
@ -547,13 +562,13 @@ static void AllFace(MetroMesh & m, VertexSampler &ps)
} }
static void AllEdge(MetroMesh & m, VertexSampler &ps) static void AllEdge(MeshType & m, VertexSampler &ps)
{ {
// Edge sampling. // Edge sampling.
typedef typename UpdateTopology<MetroMesh>::PEdge SimpleEdge; typedef typename UpdateTopology<MeshType>::PEdge SimpleEdge;
std::vector< SimpleEdge > Edges; std::vector< SimpleEdge > Edges;
typename std::vector< SimpleEdge >::iterator ei; typename std::vector< SimpleEdge >::iterator ei;
UpdateTopology<MetroMesh>::FillUniqueEdgeVector(m,Edges); UpdateTopology<MeshType>::FillUniqueEdgeVector(m,Edges);
for(ei=Edges.begin(); ei!=Edges.end(); ++ei) for(ei=Edges.begin(); ei!=Edges.end(); ++ei)
ps.AddFace(*(*ei).f,ei->EdgeBarycentricToFaceBarycentric(0.5)); ps.AddFace(*(*ei).f,ei->EdgeBarycentricToFaceBarycentric(0.5));
@ -563,11 +578,11 @@ static void AllEdge(MetroMesh & m, VertexSampler &ps)
// Each edge is subdivided in a number of pieces proprtional to its length // Each edge is subdivided in a number of pieces proprtional to its length
// Sample are choosen without touching the vertices. // Sample are choosen without touching the vertices.
static void EdgeUniform(MetroMesh & m, VertexSampler &ps,int sampleNum, bool sampleFauxEdge=true) static void EdgeUniform(MeshType & m, VertexSampler &ps,int sampleNum, bool sampleFauxEdge=true)
{ {
typedef typename UpdateTopology<MetroMesh>::PEdge SimpleEdge; typedef typename UpdateTopology<MeshType>::PEdge SimpleEdge;
std::vector< SimpleEdge > Edges; std::vector< SimpleEdge > Edges;
UpdateTopology<MetroMesh>::FillUniqueEdgeVector(m,Edges,sampleFauxEdge); UpdateTopology<MeshType>::FillUniqueEdgeVector(m,Edges,sampleFauxEdge);
// First loop compute total edge length; // First loop compute total edge length;
float edgeSum=0; float edgeSum=0;
typename std::vector< SimpleEdge >::iterator ei; typename std::vector< SimpleEdge >::iterator ei;
@ -607,9 +622,9 @@ static CoordType RandomPointInTriangle(const FaceType &f)
return f.cP(0)*u[0] + f.cP(1)*u[1] + f.cP(2)*u[2]; return f.cP(0)*u[0] + f.cP(1)*u[1] + f.cP(2)*u[2];
} }
static void StratifiedMontecarlo(MetroMesh & m, VertexSampler &ps,int sampleNum) static void StratifiedMontecarlo(MeshType & m, VertexSampler &ps,int sampleNum)
{ {
ScalarType area = Stat<MetroMesh>::ComputeMeshArea(m); ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
ScalarType samplePerAreaUnit = sampleNum/area; ScalarType samplePerAreaUnit = sampleNum/area;
// Montecarlo sampling. // Montecarlo sampling.
double floatSampleNum = 0.0; double floatSampleNum = 0.0;
@ -643,9 +658,9 @@ static void StratifiedMontecarlo(MetroMesh & m, VertexSampler &ps,int sampleNum)
*/ */
static void MontecarloPoisson(MetroMesh & m, VertexSampler &ps,int sampleNum) static void MontecarloPoisson(MeshType & m, VertexSampler &ps,int sampleNum)
{ {
ScalarType area = Stat<MetroMesh>::ComputeMeshArea(m); ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
ScalarType samplePerAreaUnit = sampleNum/area; ScalarType samplePerAreaUnit = sampleNum/area;
FaceIterator fi; FaceIterator fi;
@ -669,11 +684,11 @@ static void MontecarloPoisson(MetroMesh & m, VertexSampler &ps,int sampleNum)
and actually shooting sample over this line and actually shooting sample over this line
*/ */
static void EdgeMontecarlo(MetroMesh & m, VertexSampler &ps, int sampleNum, bool sampleAllEdges) static void EdgeMontecarlo(MeshType & m, VertexSampler &ps, int sampleNum, bool sampleAllEdges)
{ {
typedef typename UpdateTopology<MetroMesh>::PEdge SimpleEdge; typedef typename UpdateTopology<MeshType>::PEdge SimpleEdge;
std::vector< SimpleEdge > Edges; std::vector< SimpleEdge > Edges;
UpdateTopology<MetroMesh>::FillUniqueEdgeVector(m,Edges,sampleAllEdges); UpdateTopology<MeshType>::FillUniqueEdgeVector(m,Edges,sampleAllEdges);
assert(!Edges.empty()); assert(!Edges.empty());
@ -710,7 +725,7 @@ static void EdgeMontecarlo(MetroMesh & m, VertexSampler &ps, int sampleNum, bool
and actually shooting sample over this line and actually shooting sample over this line
*/ */
static void Montecarlo(MetroMesh & m, VertexSampler &ps,int sampleNum) static void Montecarlo(MeshType & m, VertexSampler &ps,int sampleNum)
{ {
typedef std::pair<ScalarType, FacePointer> IntervalType; typedef std::pair<ScalarType, FacePointer> IntervalType;
std::vector< IntervalType > intervals (m.fn+1); std::vector< IntervalType > intervals (m.fn+1);
@ -739,33 +754,38 @@ static void Montecarlo(MetroMesh & m, VertexSampler &ps,int sampleNum)
} }
} }
static ScalarType WeightedArea(FaceType f) static ScalarType WeightedArea(FaceType &f, PerVertexFloatAttribute &wH)
{ {
ScalarType averageQ = ( f.V(0)->Q() + f.V(1)->Q() + f.V(2)->Q() ) /3.0; ScalarType averageQ = ( wH[f.V(0)] + wH[f.V(1)] + wH[f.V(2)] )/3.0;
return DoubleArea(f)*averageQ/2.0; return averageQ*averageQ*DoubleArea(f)/2.0;
} }
/// Compute a sampling of the surface that is weighted by the quality /// Compute a sampling of the surface that is weighted by the quality and a variance
/// the area of each face is multiplied by the average of the quality of the vertices. ///
/// So the a face with a zero quality on all its vertices is never sampled and a face with average quality 2 get twice the samples of a face with the same area but with an average quality of 1; /// We use the quality as linear distortion of density.
static void WeightedMontecarlo(MetroMesh & m, VertexSampler &ps, int sampleNum) /// We consider each triangle as scaled between 1 and 1/variance linearly according quality.
///
/// In practice with variance 2 the average distance between sample will double where the quality is maxima.
/// If you have two same area region A with q==-1 and B with q==1, if variance==2 the A will have 4 times more samples than B
///
static void WeightedMontecarlo(MeshType & m, VertexSampler &ps,int sampleNum, float variance)
{ {
assert(tri::HasPerVertexQuality(m)); tri::RequirePerVertexQuality(m);
tri::RequireCompactness(m);
PerVertexFloatAttribute rH = tri::Allocator<MeshType>:: template GetPerVertexAttribute<float> (m,"radius");
InitRadiusHandleFromQuality(m, rH, 1.0, variance, true);
ScalarType weightedArea = 0; ScalarType weightedArea = 0;
FaceIterator fi; for(FaceIterator fi = m.face.begin(); fi != m.face.end(); ++fi)
for(fi = m.face.begin(); fi != m.face.end(); ++fi) weightedArea += WeightedArea(*fi,rH);
if(!(*fi).IsD())
weightedArea += WeightedArea(*fi);
ScalarType samplePerAreaUnit = sampleNum/weightedArea; ScalarType samplePerAreaUnit = sampleNum/weightedArea;
// Montecarlo sampling. // Montecarlo sampling.
double floatSampleNum = 0.0; double floatSampleNum = 0.0;
for(fi=m.face.begin(); fi != m.face.end(); fi++) for(FaceIterator fi=m.face.begin(); fi != m.face.end(); fi++)
if(!(*fi).IsD())
{ {
// compute # samples in the current face (taking into account of the remainders) // compute # samples in the current face (taking into account of the remainders)
floatSampleNum += WeightedArea(*fi) * samplePerAreaUnit; floatSampleNum += WeightedArea(*fi,rH) * samplePerAreaUnit;
int faceSampleNum = (int) floatSampleNum; int faceSampleNum = (int) floatSampleNum;
// for every sample p_i in T... // for every sample p_i in T...
@ -840,14 +860,14 @@ static int SingleFaceSubdivision(int sampleNum, const CoordType & v0, const Coor
/// Compute a sampling of the surface where the points are regularly scattered over the face surface using a recursive longest-edge subdivision rule. /// Compute a sampling of the surface where the points are regularly scattered over the face surface using a recursive longest-edge subdivision rule.
static void FaceSubdivision(MetroMesh & m, VertexSampler &ps,int sampleNum, bool randSample) static void FaceSubdivision(MeshType & m, VertexSampler &ps,int sampleNum, bool randSample)
{ {
ScalarType area = Stat<MetroMesh>::ComputeMeshArea(m); ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
ScalarType samplePerAreaUnit = sampleNum/area; ScalarType samplePerAreaUnit = sampleNum/area;
std::vector<FacePointer> faceVec; std::vector<FacePointer> faceVec;
FillAndShuffleFacePointerVector(m,faceVec); FillAndShuffleFacePointerVector(m,faceVec);
vcg::tri::UpdateNormal<MetroMesh>::PerFaceNormalized(m); vcg::tri::UpdateNormal<MeshType>::PerFaceNormalized(m);
double floatSampleNum = 0.0; double floatSampleNum = 0.0;
int faceSampleNum; int faceSampleNum;
// Subdivision sampling. // Subdivision sampling.
@ -931,14 +951,14 @@ static int SingleFaceSubdivisionOld(int sampleNum, const CoordType & v0, const C
/// Compute a sampling of the surface where the points are regularly scattered over the face surface using a recursive longest-edge subdivision rule. /// Compute a sampling of the surface where the points are regularly scattered over the face surface using a recursive longest-edge subdivision rule.
static void FaceSubdivisionOld(MetroMesh & m, VertexSampler &ps,int sampleNum, bool randSample) static void FaceSubdivisionOld(MeshType & m, VertexSampler &ps,int sampleNum, bool randSample)
{ {
ScalarType area = Stat<MetroMesh>::ComputeMeshArea(m); ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
ScalarType samplePerAreaUnit = sampleNum/area; ScalarType samplePerAreaUnit = sampleNum/area;
std::vector<FacePointer> faceVec; std::vector<FacePointer> faceVec;
FillAndShuffleFacePointerVector(m,faceVec); FillAndShuffleFacePointerVector(m,faceVec);
tri::UpdateNormal<MetroMesh>::PerFaceNormalized(m); tri::UpdateNormal<MeshType>::PerFaceNormalized(m);
double floatSampleNum = 0.0; double floatSampleNum = 0.0;
int faceSampleNum; int faceSampleNum;
// Subdivision sampling. // Subdivision sampling.
@ -1038,11 +1058,11 @@ static int SingleFaceSimilarDual(FacePointer fp, VertexSampler &ps, int n_sample
//template <class MetroMesh> //template <class MeshType>
//void Sampling<MetroMesh>::SimilarFaceSampling() //void Sampling<MeshType>::SimilarFaceSampling()
static void FaceSimilar(MetroMesh & m, VertexSampler &ps,int sampleNum, bool dualFlag, bool randomFlag) static void FaceSimilar(MeshType & m, VertexSampler &ps,int sampleNum, bool dualFlag, bool randomFlag)
{ {
ScalarType area = Stat<MetroMesh>::ComputeMeshArea(m); ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
ScalarType samplePerAreaUnit = sampleNum/area; ScalarType samplePerAreaUnit = sampleNum/area;
// Similar Triangles sampling. // Similar Triangles sampling.
@ -1081,13 +1101,13 @@ static void FaceSimilar(MetroMesh & m, VertexSampler &ps,int sampleNum, bool dua
// This function does rasterization with a safety buffer area, thus accounting some points actually outside triangle area // This function does rasterization with a safety buffer area, thus accounting some points actually outside triangle area
// The safety area samples are generated according to face flag BORDER which should be true for texture space border edges // The safety area samples are generated according to face flag BORDER which should be true for texture space border edges
// Use correctSafePointsBaryCoords = true to map safety texels to closest point barycentric coords (on edge). // Use correctSafePointsBaryCoords = true to map safety texels to closest point barycentric coords (on edge).
static void SingleFaceRaster(typename MetroMesh::FaceType &f, VertexSampler &ps, static void SingleFaceRaster(typename MeshType::FaceType &f, VertexSampler &ps,
const Point2<typename MetroMesh::ScalarType> & v0, const Point2<typename MeshType::ScalarType> & v0,
const Point2<typename MetroMesh::ScalarType> & v1, const Point2<typename MeshType::ScalarType> & v1,
const Point2<typename MetroMesh::ScalarType> & v2, const Point2<typename MeshType::ScalarType> & v2,
bool correctSafePointsBaryCoords=true) bool correctSafePointsBaryCoords=true)
{ {
typedef typename MetroMesh::ScalarType S; typedef typename MeshType::ScalarType S;
// Calcolo bounding box // Calcolo bounding box
Box2i bbox; Box2i bbox;
Box2<S> bboxf; Box2<S> bboxf;
@ -1153,7 +1173,7 @@ static void FaceSimilar(MetroMesh & m, VertexSampler &ps,int sampleNum, bool dua
{ {
if( ((n[0]>=0 && n[1]>=0 && n[2]>=0) || (n[0]<=0 && n[1]<=0 && n[2]<=0)) && (de != 0)) if( ((n[0]>=0 && n[1]>=0 && n[2]>=0) || (n[0]<=0 && n[1]<=0 && n[2]<=0)) && (de != 0))
{ {
typename MetroMesh::CoordType baryCoord; typename MeshType::CoordType baryCoord;
baryCoord[0] = double(-y*v1[0]+v2[0]*y+v1[1]*x-v2[0]*v1[1]+v1[0]*v2[1]-x*v2[1])/de; baryCoord[0] = double(-y*v1[0]+v2[0]*y+v1[1]*x-v2[0]*v1[1]+v1[0]*v2[1]-x*v2[1])/de;
baryCoord[1] = -double( x*v0[1]-x*v2[1]-v0[0]*y+v0[0]*v2[1]-v2[0]*v0[1]+v2[0]*y)/de; baryCoord[1] = -double( x*v0[1]-x*v2[1]-v0[0]*y+v0[0]*v2[1]-v2[0]*v0[1]+v2[0]*y)/de;
baryCoord[2] = 1-baryCoord[0]-baryCoord[1]; baryCoord[2] = 1-baryCoord[0]-baryCoord[1];
@ -1192,7 +1212,7 @@ static void FaceSimilar(MetroMesh & m, VertexSampler &ps,int sampleNum, bool dua
if (closeEdge >= 0) if (closeEdge >= 0)
{ {
typename MetroMesh::CoordType baryCoord; typename MeshType::CoordType baryCoord;
if (correctSafePointsBaryCoords) if (correctSafePointsBaryCoords)
{ {
// Add x,y sample with closePoint barycentric coords (on edge) // Add x,y sample with closePoint barycentric coords (on edge)
@ -1224,7 +1244,7 @@ static bool checkPoissonDisk(SampleSHT & sht, const Point3<ScalarType> & p, Scal
{ {
// get the samples closest to the given one // get the samples closest to the given one
std::vector<VertexType*> closests; std::vector<VertexType*> closests;
typedef VertTmark<MetroMesh> MarkerVert; typedef VertTmark<MeshType> MarkerVert;
static MarkerVert mv; static MarkerVert mv;
Box3f bb(p-Point3f(radius,radius,radius),p+Point3f(radius,radius,radius)); Box3f bb(p-Point3f(radius,radius,radius),p+Point3f(radius,radius,radius));
@ -1251,7 +1271,6 @@ struct PoissonDiskParam
preGenFlag = false; preGenFlag = false;
preGenMesh = NULL; preGenMesh = NULL;
geodesicDistanceFlag = false; geodesicDistanceFlag = false;
pds=NULL;
} }
struct Stat struct Stat
@ -1272,11 +1291,13 @@ struct PoissonDiskParam
bool adaptiveRadiusFlag; bool adaptiveRadiusFlag;
float radiusVariance; float radiusVariance;
bool invertQuality; bool invertQuality;
bool preGenFlag; // when generating a poisson distribution, you can initialize the set of computed points with ALL the vertices of another mesh. Useful for building progressive refinements. bool preGenFlag; // when generating a poisson distribution, you can initialize the set of computed points with
MetroMesh *preGenMesh; // ALL the vertices of another mesh. Useful for building progressive//prioritize refinements.
MeshType *preGenMesh; // There are two ways of passing the pregen vertexes to the pruning, 1) is with a mesh pointer
// 2) with a per vertex attribute.
int MAXLEVELS; int MAXLEVELS;
Stat *pds; Stat pds;
}; };
@ -1315,9 +1336,9 @@ static VertexPointer getBestPrecomputedMontecarloSample(Point3i &cell, Montecarl
} }
static ScalarType ComputePoissonDiskRadius(MetroMesh &origMesh, int sampleNum) static ScalarType ComputePoissonDiskRadius(MeshType &origMesh, int sampleNum)
{ {
ScalarType meshArea = Stat<MetroMesh>::ComputeMeshArea(origMesh); ScalarType meshArea = Stat<MeshType>::ComputeMeshArea(origMesh);
// Manage approximately the PointCloud Case, use the half a area of the bbox. // Manage approximately the PointCloud Case, use the half a area of the bbox.
// TODO: If you had the radius a much better approximation could be done. // TODO: If you had the radius a much better approximation could be done.
if(meshArea ==0) if(meshArea ==0)
@ -1330,28 +1351,30 @@ static ScalarType ComputePoissonDiskRadius(MetroMesh &origMesh, int sampleNum)
return diskRadius; return diskRadius;
} }
static int ComputePoissonSampleNum(MetroMesh &origMesh, ScalarType diskRadius) static int ComputePoissonSampleNum(MeshType &origMesh, ScalarType diskRadius)
{ {
ScalarType meshArea = Stat<MetroMesh>::ComputeMeshArea(origMesh); ScalarType meshArea = Stat<MeshType>::ComputeMeshArea(origMesh);
int sampleNum = meshArea / (diskRadius*diskRadius *M_PI *0.7) ; // 0.7 is a density factor int sampleNum = meshArea / (diskRadius*diskRadius *M_PI *0.7) ; // 0.7 is a density factor
return sampleNum; return sampleNum;
} }
// Note that this function actually CHANGE the quality of the montecarlo samples so that it represents the expected radius of each sample /// When performing an adptive pruning for each sample we expect a varying radius to be removed.
// the expected radius of the sample is computed so that it linearly maps the quality between diskradius / variance and diskradius*variance /// The radius is a PerVertex attribute that we compute from the current quality
///
/// the expected radius of the sample is computed so that
/// it linearly maps the quality between diskradius and diskradius*variance
/// in other words the radius
static void ComputePoissonSampleRadii(MetroMesh &sampleMesh, ScalarType diskRadius, ScalarType radiusVariance, bool invert) static void InitRadiusHandleFromQuality(MeshType &sampleMesh, PerVertexFloatAttribute &rH, ScalarType diskRadius, ScalarType radiusVariance, bool invert)
{ {
VertexIterator vi; std::pair<float,float> minmax = tri::Stat<MeshType>::ComputePerVertexQualityMinMax( sampleMesh);
std::pair<float,float> minmax = tri::Stat<MetroMesh>::ComputePerVertexQualityMinMax( sampleMesh); float minRad = diskRadius ;
float minRad = diskRadius / radiusVariance;
float maxRad = diskRadius * radiusVariance; float maxRad = diskRadius * radiusVariance;
float deltaQ = minmax.second-minmax.first; float deltaQ = minmax.second-minmax.first;
float deltaRad = maxRad-minRad; float deltaRad = maxRad-minRad;
for (vi = sampleMesh.vert.begin(); vi != sampleMesh.vert.end(); vi++) for (VertexIterator vi = sampleMesh.vert.begin(); vi != sampleMesh.vert.end(); vi++)
if(!(*vi).IsD())
{ {
(*vi).Q() = minRad + deltaRad*((invert ? minmax.second - (*vi).Q() : (*vi).Q() - minmax.first )/deltaQ); rH[*vi] = minRad + deltaRad*((invert ? minmax.second - (*vi).Q() : (*vi).Q() - minmax.first )/deltaQ);
} }
} }
@ -1360,8 +1383,8 @@ static void ComputePoissonSampleRadii(MetroMesh &sampleMesh, ScalarType diskRadi
// This radius implies that when we pick a sample in a cell all that cell probably will not be touched again. // This radius implies that when we pick a sample in a cell all that cell probably will not be touched again.
// Howvever we must ensure that we do not put too many vertices inside each hash cell // Howvever we must ensure that we do not put too many vertices inside each hash cell
static void InitSpatialHashTable(MetroMesh &montecarloMesh, MontecarloSHT &montecarloSHT, ScalarType diskRadius, static void InitSpatialHashTable(MeshType &montecarloMesh, MontecarloSHT &montecarloSHT, ScalarType diskRadius,
const struct PoissonDiskParam pp=PoissonDiskParam()) struct PoissonDiskParam pp=PoissonDiskParam())
{ {
ScalarType cellsize = 2.0f* diskRadius / sqrt(3.0); ScalarType cellsize = 2.0f* diskRadius / sqrt(3.0);
float occupancyRatio=0; float occupancyRatio=0;
@ -1386,11 +1409,8 @@ static void InitSpatialHashTable(MetroMesh &montecarloMesh, MontecarloSHT &monte
} }
montecarloSHT.UpdateAllocatedCells(); montecarloSHT.UpdateAllocatedCells();
if(pp.pds) pp.pds.gridSize = gridsize;
{ pp.pds.gridCellNum = (int)montecarloSHT.AllocatedCells.size();
pp.pds->gridSize = gridsize;
pp.pds->gridCellNum = (int)montecarloSHT.AllocatedCells.size();
}
cellsize/=2.0f; cellsize/=2.0f;
occupancyRatio = float(montecarloMesh.vn) / float(montecarloSHT.AllocatedCells.size()); occupancyRatio = float(montecarloMesh.vn) / float(montecarloSHT.AllocatedCells.size());
// qDebug(" %i / %i = %6.3f", montecarloMesh.vn , montecarloSHT.AllocatedCells.size(),occupancyRatio); // qDebug(" %i / %i = %6.3f", montecarloMesh.vn , montecarloSHT.AllocatedCells.size(),occupancyRatio);
@ -1398,10 +1418,69 @@ static void InitSpatialHashTable(MetroMesh &montecarloMesh, MontecarloSHT &monte
while( occupancyRatio> 100); while( occupancyRatio> 100);
} }
// Trivial approach that puts all the samples in a UG and removes all the ones that surely do not fit the static void PoissonDiskPruningByNumber(VertexSampler &ps, MeshType &m,
static void PoissonDiskPruning(VertexSampler &ps, MetroMesh &montecarloMesh, int sampleNum, ScalarType &diskRadius,
ScalarType diskRadius, const struct PoissonDiskParam pp=PoissonDiskParam()) PoissonDiskParam &pp,
float tolerance=0.04,
int maxIter=20)
{ {
size_t sampleNumMin = int(float(sampleNum)*(1.0f-tolerance));
size_t sampleNumMax = int(float(sampleNum)*(1.0f+tolerance));
float RangeMinRad = m.bbox.Diag()/50.0;
float RangeMaxRad = m.bbox.Diag()/50.0;
size_t RangeMinRadNum;
size_t RangeMaxRadNum;
// Note RangeMinRad < RangeMaxRad
// but RangeMinRadNum > sampleNum > RangeMaxRadNum
do {
ps.reset();
RangeMinRad/=2.0f;
PoissonDiskPruning(ps, m ,RangeMinRad,pp);
RangeMinRadNum = pp.pds.sampleNum;
// qDebug("PoissonDiskPruning Iteratin Min (%6.3f:%5i) instead of %i",RangeMinRad,RangeMinRadNum,sampleNum);
} while(RangeMinRadNum < sampleNum); // if the number of sample is still smaller you have to make radius larger.
do {
ps.reset();
RangeMaxRad*=2.0f;
PoissonDiskPruning(ps, m ,RangeMaxRad,pp);
RangeMaxRadNum = pp.pds.sampleNum;
// qDebug("PoissonDiskPruning Iteratin Max (%6.3f:%5i) instead of %i",RangeMaxRad,RangeMaxRadNum,sampleNum);
} while(RangeMaxRadNum > sampleNum);
float curRadius;
int iterCnt=0;
while(iterCnt<maxIter &&
(pp.pds.sampleNum < sampleNumMin || pp.pds.sampleNum > sampleNumMax) )
{
iterCnt++;
ps.reset();
curRadius=(RangeMaxRad+RangeMinRad)/2.0f;
PoissonDiskPruning(ps, m ,curRadius,pp);
qDebug("PoissonDiskPruning Iteratin (%6.3f:%5i %6.3f:%5i) Cur Radius %f -> %i sample instead of %i",RangeMinRad,RangeMinRadNum,RangeMaxRad,RangeMaxRadNum,curRadius,pp.pds.sampleNum,sampleNum);
if(pp.pds.sampleNum > sampleNum){
RangeMinRad = curRadius;
RangeMinRadNum = pp.pds.sampleNum;
}
if(pp.pds.sampleNum < sampleNum){
RangeMaxRad = curRadius;
RangeMaxRadNum = pp.pds.sampleNum;
}
}
diskRadius = curRadius;
}
/// This is the main function that is used to build a poisson distribuition
/// starting from a dense sample cloud.
/// Trivial approach that puts all the samples in a hashed UG and randomly choose a sample
/// and remove all the points in the sphere centered on the chosen sample
static void PoissonDiskPruning(VertexSampler &ps, MeshType &montecarloMesh,
ScalarType diskRadius, PoissonDiskParam &pp)
{
tri::RequireCompactness(montecarloMesh);
if(pp.adaptiveRadiusFlag) if(pp.adaptiveRadiusFlag)
tri::RequirePerVertexQuality(montecarloMesh); tri::RequirePerVertexQuality(montecarloMesh);
int t0 = clock(); int t0 = clock();
@ -1409,27 +1488,41 @@ static void PoissonDiskPruning(VertexSampler &ps, MetroMesh &montecarloMesh,
MontecarloSHT montecarloSHT; MontecarloSHT montecarloSHT;
InitSpatialHashTable(montecarloMesh,montecarloSHT,diskRadius,pp); InitSpatialHashTable(montecarloMesh,montecarloSHT,diskRadius,pp);
// if we are doing variable density sampling we have to prepare the random samples quality with the correct expected radii. // if we are doing variable density sampling we have to prepare the handle that keeps the the random samples expected radii.
// At this point we just assume that there is the quality values as sampled from the base mesh // At this point we just assume that there is the quality values as sampled from the base mesh
PerVertexFloatAttribute rH = tri::Allocator<MeshType>:: template GetPerVertexAttribute<float> (montecarloMesh,"radius");
if(pp.adaptiveRadiusFlag) if(pp.adaptiveRadiusFlag)
ComputePoissonSampleRadii(montecarloMesh, diskRadius, pp.radiusVariance, pp.invertQuality); InitRadiusHandleFromQuality(montecarloMesh, rH, diskRadius, pp.radiusVariance, pp.invertQuality);
unsigned int (*p_myrandom)(unsigned int) = RandomInt; unsigned int (*p_myrandom)(unsigned int) = RandomInt;
std::random_shuffle(montecarloSHT.AllocatedCells.begin(),montecarloSHT.AllocatedCells.end(), p_myrandom); std::random_shuffle(montecarloSHT.AllocatedCells.begin(),montecarloSHT.AllocatedCells.end(), p_myrandom);
int t1 = clock(); int t1 = clock();
if(pp.pds) { pp.pds.montecarloSampleNum = montecarloMesh.vn;
pp.pds->montecarloSampleNum = montecarloMesh.vn; pp.pds.sampleNum =0;
}
int removedCnt=0; int removedCnt=0;
if(pp.preGenFlag && pp.preGenMesh !=0)
{
// Initial pass for pruning the Hashed grid with the an eventual pre initialized set of samples // Initial pass for pruning the Hashed grid with the an eventual pre initialized set of samples
for(VertexIterator vi =pp.preGenMesh->vert.begin(); vi!=pp.preGenMesh->vert.end();++vi) if(pp.preGenFlag)
if(!(*vi).IsD())
{ {
if(pp.preGenMesh==0)
{
typename MeshType::template PerVertexAttributeHandle<bool> fixed;
fixed = tri::Allocator<MeshType>:: template GetPerVertexAttribute<bool> (montecarloMesh,"fixed");
for(VertexIterator vi=montecarloMesh.vert.begin();vi!=montecarloMesh.vert.end();++vi)
if(fixed[*vi]) {
pp.pds.sampleNum++;
ps.AddVert(*vi); ps.AddVert(*vi);
removedCnt += montecarloSHT.RemoveInSphere(vi->cP(),diskRadius); removedCnt += montecarloSHT.RemoveInSphere(vi->cP(),diskRadius);
} }
}
else
{
for(VertexIterator vi =pp.preGenMesh->vert.begin(); vi!=pp.preGenMesh->vert.end();++vi)
{
ps.AddVert(*vi);
pp.pds.sampleNum++;
removedCnt += montecarloSHT.RemoveInSphere(vi->cP(),diskRadius);
}
}
montecarloSHT.UpdateAllocatedCells(); montecarloSHT.UpdateAllocatedCells();
} }
vertex::ApproximateGeodesicDistanceFunctor<VertexType> GDF; vertex::ApproximateGeodesicDistanceFunctor<VertexType> GDF;
@ -1447,20 +1540,18 @@ static void PoissonDiskPruning(VertexSampler &ps, MetroMesh &montecarloMesh,
sp = getSampleFromCell(montecarloSHT.AllocatedCells[i], montecarloSHT); sp = getSampleFromCell(montecarloSHT.AllocatedCells[i], montecarloSHT);
if(pp.adaptiveRadiusFlag) if(pp.adaptiveRadiusFlag)
currentRadius = sp->Q(); currentRadius = rH[sp];
ps.AddVert(*sp); ps.AddVert(*sp);
pp.pds.sampleNum++;
if(pp.geodesicDistanceFlag) removedCnt += montecarloSHT.RemoveInSphereNormal(sp->cP(),sp->cN(),GDF,currentRadius); if(pp.geodesicDistanceFlag) removedCnt += montecarloSHT.RemoveInSphereNormal(sp->cP(),sp->cN(),GDF,currentRadius);
else removedCnt += montecarloSHT.RemoveInSphere(sp->cP(),currentRadius); else removedCnt += montecarloSHT.RemoveInSphere(sp->cP(),currentRadius);
} }
montecarloSHT.UpdateAllocatedCells(); montecarloSHT.UpdateAllocatedCells();
} }
int t2 = clock(); int t2 = clock();
if(pp.pds) pp.pds.gridTime = t1-t0;
{ pp.pds.pruneTime = t2-t1;
pp.pds->gridTime = t1-t0;
pp.pds->pruneTime = t2-t1;
}
} }
/** Compute a Poisson-disk sampling of the surface. /** Compute a Poisson-disk sampling of the surface.
@ -1473,7 +1564,7 @@ static void PoissonDiskPruning(VertexSampler &ps, MetroMesh &montecarloMesh,
* IEEE Symposium on Interactive Ray Tracing, 2007, * IEEE Symposium on Interactive Ray Tracing, 2007,
* 10-12 Sept. 2007, pp. 129-132. * 10-12 Sept. 2007, pp. 129-132.
*/ */
static void HierarchicalPoissonDisk(MetroMesh &origMesh, VertexSampler &ps, MetroMesh &montecarloMesh, ScalarType diskRadius, const struct PoissonDiskParam pp=PoissonDiskParam()) static void HierarchicalPoissonDisk(MeshType &origMesh, VertexSampler &ps, MeshType &montecarloMesh, ScalarType diskRadius, const struct PoissonDiskParam pp=PoissonDiskParam())
{ {
// int t0=clock(); // int t0=clock();
// spatial index of montecarlo samples - used to choose a new sample to insert // spatial index of montecarlo samples - used to choose a new sample to insert
@ -1520,8 +1611,9 @@ static void HierarchicalPoissonDisk(MetroMesh &origMesh, VertexSampler &ps, Metr
montecarloSHTVec[0].UpdateAllocatedCells(); montecarloSHTVec[0].UpdateAllocatedCells();
// if we are doing variable density sampling we have to prepare the random samples quality with the correct expected radii. // if we are doing variable density sampling we have to prepare the random samples quality with the correct expected radii.
PerVertexFloatAttribute rH = tri::Allocator<MeshType>:: template GetPerVertexAttribute<float> (montecarloMesh,"radius");
if(pp.adaptiveRadiusFlag) if(pp.adaptiveRadiusFlag)
ComputePoissonSampleRadii(montecarloMesh, diskRadius, pp.radiusVariance, pp.invertQuality); InitRadiusHandleFromQuality(montecarloMesh, rH, diskRadius, pp.radiusVariance, pp.invertQuality);
do do
{ {
@ -1556,7 +1648,7 @@ static void HierarchicalPoissonDisk(MetroMesh &origMesh, VertexSampler &ps, Metr
VertexPointer sp = (*hi).second; VertexPointer sp = (*hi).second;
// vr spans between 3.0*r and r / 4.0 according to vertex quality // vr spans between 3.0*r and r / 4.0 according to vertex quality
ScalarType sampleRadius = diskRadius; ScalarType sampleRadius = diskRadius;
if(pp.adaptiveRadiusFlag) sampleRadius = sp->Q(); if(pp.adaptiveRadiusFlag) sampleRadius = rH[sp];
if (checkPoissonDisk(checkSHT, sp->cP(), sampleRadius)) if (checkPoissonDisk(checkSHT, sp->cP(), sampleRadius))
{ {
ps.AddVert(*sp); ps.AddVert(*sp);
@ -1580,8 +1672,8 @@ static void HierarchicalPoissonDisk(MetroMesh &origMesh, VertexSampler &ps, Metr
} while(level < 5); } while(level < 5);
} }
//template <class MetroMesh> //template <class MeshType>
//void Sampling<MetroMesh>::SimilarFaceSampling() //void Sampling<MeshType>::SimilarFaceSampling()
// This function also generates samples outside faces if those affects faces in texture space. // This function also generates samples outside faces if those affects faces in texture space.
// Use correctSafePointsBaryCoords = true to map safety texels to closest point barycentric coords (on edge) // Use correctSafePointsBaryCoords = true to map safety texels to closest point barycentric coords (on edge)
@ -1593,7 +1685,7 @@ static void HierarchicalPoissonDisk(MetroMesh &origMesh, VertexSampler &ps, Metr
// Else make sure to update border flags from texture space FFadj // Else make sure to update border flags from texture space FFadj
// vcg::tri::UpdateTopology<Mesh>::FaceFaceFromTexCoord(m); // vcg::tri::UpdateTopology<Mesh>::FaceFaceFromTexCoord(m);
// vcg::tri::UpdateFlags<Mesh>::FaceBorderFromFF(m); // vcg::tri::UpdateFlags<Mesh>::FaceBorderFromFF(m);
static void Texture(MetroMesh & m, VertexSampler &ps, int textureWidth, int textureHeight, bool correctSafePointsBaryCoords=true) static void Texture(MeshType & m, VertexSampler &ps, int textureWidth, int textureHeight, bool correctSafePointsBaryCoords=true)
{ {
FaceIterator fi; FaceIterator fi;
@ -1617,11 +1709,11 @@ class RRParam
public: public:
float offset; float offset;
float minDiag; float minDiag;
tri::FaceTmark<MetroMesh> markerFunctor; tri::FaceTmark<MeshType> markerFunctor;
TriMeshGrid gM; TriMeshGrid gM;
}; };
static void RegularRecursiveOffset(MetroMesh & m, std::vector<Point3f> &pvec, ScalarType offset, float minDiag) static void RegularRecursiveOffset(MeshType & m, std::vector<Point3f> &pvec, ScalarType offset, float minDiag)
{ {
Box3<ScalarType> bb=m.bbox; Box3<ScalarType> bb=m.bbox;
bb.Offset(offset*2.0); bb.Offset(offset*2.0);
@ -1638,7 +1730,7 @@ static void RegularRecursiveOffset(MetroMesh & m, std::vector<Point3f> &pvec, Sc
SubdivideAndSample(m, pvec, bb, rrp, bb.Diag()); SubdivideAndSample(m, pvec, bb, rrp, bb.Diag());
} }
static void SubdivideAndSample(MetroMesh & m, std::vector<Point3f> &pvec, const Box3<ScalarType> bb, RRParam &rrp, float curDiag) static void SubdivideAndSample(MeshType & m, std::vector<Point3f> &pvec, const Box3<ScalarType> bb, RRParam &rrp, float curDiag)
{ {
Point3f startPt = bb.Center(); Point3f startPt = bb.Center();
@ -1679,6 +1771,22 @@ static void SubdivideAndSample(MetroMesh & m, std::vector<Point3f> &pvec, const
} }
}; // end sampling class }; // end sampling class
template <class MeshType>
typename MeshType::ScalarType ComputePoissonDiskRadius(MeshType &origMesh, int sampleNum)
{
typedef typename MeshType::ScalarType ScalarType;
ScalarType meshArea = Stat<MeshType>::ComputeMeshArea(origMesh);
// Manage approximately the PointCloud Case, use the half a area of the bbox.
// TODO: If you had the radius a much better approximation could be done.
if(meshArea ==0)
{
meshArea = (origMesh.bbox.DimX()*origMesh.bbox.DimY() +
origMesh.bbox.DimX()*origMesh.bbox.DimZ() +
origMesh.bbox.DimY()*origMesh.bbox.DimZ());
}
ScalarType diskRadius = sqrt(meshArea / (0.7 * M_PI * sampleNum)); // 0.7 is a density factor
return diskRadius;
}
@ -1715,17 +1823,14 @@ void PoissonSampling(MeshType &m, // the mesh that has to be sampled
{ {
typedef tri::TrivialSampler<MeshType> BaseSampler; typedef tri::TrivialSampler<MeshType> BaseSampler;
typedef tri::MeshSampler<MeshType> MontecarloSampler; typedef tri::MeshSampler<MeshType> MontecarloSampler;
tri::RequireVFAdjacency(m);
typename tri::SurfaceSampling<MeshType, BaseSampler>::PoissonDiskParam pp; typename tri::SurfaceSampling<MeshType, BaseSampler>::PoissonDiskParam pp;
typename tri::SurfaceSampling<MeshType, BaseSampler>::PoissonDiskParam::Stat stat;
pp.pds = &stat;
int t0=clock(); int t0=clock();
if(sampleNum>0) radius = tri::SurfaceSampling<MeshType,BaseSampler>::ComputePoissonDiskRadius(m,sampleNum); if(sampleNum>0) radius = tri::SurfaceSampling<MeshType,BaseSampler>::ComputePoissonDiskRadius(m,sampleNum);
if(radius>0 && sampleNum==0) sampleNum = tri::SurfaceSampling<MeshType,BaseSampler>::ComputePoissonSampleNum(m,radius); if(radius>0 && sampleNum==0) sampleNum = tri::SurfaceSampling<MeshType,BaseSampler>::ComputePoissonSampleNum(m,radius);
pp.pds->sampleNum = sampleNum; pp.pds.sampleNum = sampleNum;
poissonSamples.clear(); poissonSamples.clear();
// std::vector<Point3f> MontecarloSamples; // std::vector<Point3f> MontecarloSamples;
MeshType MontecarloMesh; MeshType MontecarloMesh;
@ -1738,7 +1843,7 @@ void PoissonSampling(MeshType &m, // the mesh that has to be sampled
tri::UpdateBounding<MeshType>::Box(MontecarloMesh); tri::UpdateBounding<MeshType>::Box(MontecarloMesh);
// tri::Build(MontecarloMesh, MontecarloSamples); // tri::Build(MontecarloMesh, MontecarloSamples);
int t1=clock(); int t1=clock();
pp.pds->montecarloTime = t1-t0; pp.pds.montecarloTime = t1-t0;
if(radiusVariance !=1) if(radiusVariance !=1)
{ {
@ -1747,42 +1852,22 @@ void PoissonSampling(MeshType &m, // the mesh that has to be sampled
} }
tri::SurfaceSampling<MeshType,BaseSampler>::PoissonDiskPruning(pdSampler, MontecarloMesh, radius,pp); tri::SurfaceSampling<MeshType,BaseSampler>::PoissonDiskPruning(pdSampler, MontecarloMesh, radius,pp);
int t2=clock(); int t2=clock();
pp.pds->totalTime = t2-t0; pp.pds.totalTime = t2-t0;
}
/// \brief Very simple wrapping for the Poisson Disk Pruning
///
/// This function simply takes a mesh and a radius and returns a vector of points
//
template <class MeshType>
void PoissonPruning(MeshType &m, // the mesh that has to be pruned
std::vector<Point3f> &poissonSamples, // the vector that will contain the chosen set of points
float & radius, int sampleNum=0)
{
typedef tri::TrivialSampler<MeshType> BaseSampler;
typename tri::SurfaceSampling<MeshType, BaseSampler>::PoissonDiskParam pp;
if(sampleNum>0 && radius==0)
radius = tri::SurfaceSampling<MeshType,BaseSampler>::ComputePoissonDiskRadius(m,sampleNum);
tri::UpdateBounding<MeshType>::Box(m);
BaseSampler pdSampler(poissonSamples);
tri::SurfaceSampling<MeshType,BaseSampler>::PoissonDiskPruning(pdSampler, m, radius,pp);
} }
/// \brief Low level wrapper for Poisson Disk Pruning
/// \brief Very simple wrapping for the Poisson Disk Pruning
/// ///
/// This function simply takes a mesh and a radius and returns a vector of points /// This function simply takes a mesh and a radius and returns a vector of vertex pointers listing the "surviving" points.
// //
template <class MeshType> template <class MeshType>
void PoissonPruning(MeshType &m, // the mesh that has to be pruned void PoissonPruning(MeshType &m, // the mesh that has to be pruned
std::vector<typename MeshType::VertexPointer> &poissonSamples, // the vector that will contain the chosen set of points std::vector<typename MeshType::VertexPointer> &poissonSamples, // the vector that will contain the chosen set of points
float & radius, int sampleNum=0, unsigned int randSeed=0) float radius, unsigned int randSeed=0)
{ {
typedef tri::TrivialPointerSampler<MeshType> BaseSampler; typedef tri::TrivialPointerSampler<MeshType> BaseSampler;
typename tri::SurfaceSampling<MeshType, BaseSampler>::PoissonDiskParam pp; typename tri::SurfaceSampling<MeshType, BaseSampler>::PoissonDiskParam pp;
if(randSeed !=0) if(randSeed !=0)
tri::SurfaceSampling<MeshType,BaseSampler>::SamplingRandomGenerator().initialize(randSeed); tri::SurfaceSampling<MeshType,BaseSampler>::SamplingRandomGenerator().initialize(randSeed);
if(sampleNum>0 && radius==0)
radius = tri::SurfaceSampling<MeshType,BaseSampler>::ComputePoissonDiskRadius(m,sampleNum);
tri::UpdateBounding<MeshType>::Box(m); tri::UpdateBounding<MeshType>::Box(m);
BaseSampler pdSampler; BaseSampler pdSampler;
@ -1790,38 +1875,73 @@ void PoissonPruning(MeshType &m, // the mesh that has to be pruned
std::swap(pdSampler.sampleVec,poissonSamples); std::swap(pdSampler.sampleVec,poissonSamples);
} }
/// \brief Low level wrapper for Poisson Disk Pruning
///
/// This function simply takes a mesh and a radius and returns a vector
/// of vertex pointers listing the "surviving" points.
//
template <class MeshType>
void PoissonPruning(MeshType &m, // the mesh that has to be pruned
std::vector<Point3f> &poissonSamples, // the vector that will contain the chosen set of points
float radius, unsigned int randSeed=0)
{
std::vector<typename MeshType::VertexPointer> poissonSamplesVP;
PoissonPruning(m,poissonSamplesVP,radius,randSeed);
poissonSamples.resize(poissonSamplesVP.size());
for(size_t i=0;i<poissonSamplesVP.size();++i)
poissonSamples[i]=poissonSamplesVP[i]->P();
}
/// \brief Very simple wrapping for the Exact Poisson Disk Pruning /// \brief Very simple wrapping for the Exact Poisson Disk Pruning
/// ///
/// This function simply takes a mesh and an expected number of points and returns /// This function simply takes a mesh and an expected number of points and returns
/// vector of points. It performs multiple attempts with varius radii to correctly get the expected number of samples. /// vector of points. It performs multiple attempts with varius radii to correctly get the expected number of samples.
/// It is obviously much slower than the other versions... /// It is obviously much slower than the other versions...
template <class MeshType> template <class MeshType>
void PoissonPruningExact(MeshType &m, // the mesh that has to be pruned void PoissonPruningExact(MeshType &m, /// the mesh that has to be pruned
std::vector<Point3f> &poissonSamples, // the vector that will contain the chosen set of points std::vector<typename MeshType::VertexPointer> &poissonSamples, /// the vector that will contain the chosen set of points
float & radius, float & radius,
int sampleNum, int sampleNum,
float tolerance=0.04, float tolerance=0.04,
int maxIter=20) int maxIter=20)
{ {
int sampleNumMin = int(float(sampleNum)*(1.0f-tolerance)); size_t sampleNumMin = int(float(sampleNum)*(1.0f-tolerance));
int sampleNumMax = int(float(sampleNum)*(1.0f+tolerance)); size_t sampleNumMax = int(float(sampleNum)*(1.0f+tolerance));
float RadiusRangeMin = m.bbox.Diag()/100.0; float RangeMinRad = m.bbox.Diag()/10.0f;
float RadiusRangeMax = m.bbox.Diag(); float RangeMaxRad = m.bbox.Diag()/10.0f;
std::vector<Point3f> poissonSamplesTmp; size_t RangeMinNum;
size_t RangeMaxNum;
std::vector<typename MeshType::VertexPointer> poissonSamplesTmp;
do
{
RangeMinRad/=2.0f;
PoissonPruning(m,poissonSamplesTmp,RangeMinRad);
RangeMinNum = poissonSamplesTmp.size();
} while(RangeMinNum > sampleNumMin);
do
{
RangeMaxRad*=2.0f;
PoissonPruning(m,poissonSamplesTmp,RangeMaxRad);
RangeMaxNum = poissonSamplesTmp.size();
} while(RangeMaxNum < sampleNumMax);
float curRadius; float curRadius;
int iterCnt=0; int iterCnt=0;
while(iterCnt<maxIter && while(iterCnt<maxIter &&
(poissonSamplesTmp.size() < sampleNumMin || poissonSamplesTmp.size() > sampleNumMax) ) (poissonSamplesTmp.size() < sampleNumMin || poissonSamplesTmp.size() > sampleNumMax) )
{ {
curRadius=(RadiusRangeMax+RadiusRangeMin)/2.0f; curRadius=(RangeMaxRad+RangeMinRad)/2.0f;
PoissonPruning(m,poissonSamplesTmp,curRadius); PoissonPruning(m,poissonSamplesTmp,curRadius);
//qDebug("(%f %f) Cur Radius %f -> %i sample instead of %i",RadiusRangeMin,RadiusRangeMax,curRadius,poissonSamplesTmp.size(),sampleNum); qDebug("(%6.3f:%5i %6.3f:%5i) Cur Radius %f -> %i sample instead of %i",RangeMinRad,RangeMinNum,RangeMaxRad,RangeMaxNum,curRadius,poissonSamplesTmp.size(),sampleNum);
if(poissonSamplesTmp.size() > sampleNum) if(poissonSamplesTmp.size() > sampleNum)
RadiusRangeMin = curRadius; RangeMinRad = curRadius;
if(poissonSamplesTmp.size() < sampleNum) if(poissonSamplesTmp.size() < sampleNum)
RadiusRangeMax = curRadius; RangeMaxRad = curRadius;
} }
swap(poissonSamples,poissonSamplesTmp); swap(poissonSamples,poissonSamplesTmp);