diff --git a/vcg/complex/algorithms/point_sampling.h b/vcg/complex/algorithms/point_sampling.h index addbed1e..080bef46 100644 --- a/vcg/complex/algorithms/point_sampling.h +++ b/vcg/complex/algorithms/point_sampling.h @@ -39,6 +39,7 @@ sampling strategies (montecarlo, stratified etc). #include #include #include +#include #include #include #include @@ -75,6 +76,11 @@ public: typedef typename MeshType::VertexType VertexType; typedef typename MeshType::FaceType FaceType; + void reset() + { + sampleVec->clear(); + } + TrivialSampler() { sampleVec = new std::vector(); @@ -84,8 +90,8 @@ public: TrivialSampler(std::vector &Vec) { sampleVec = &Vec; - sampleVec->clear(); vectorOwner=false; + reset(); } ~TrivialSampler() @@ -126,6 +132,11 @@ public: TrivialPointerSampler() {} ~TrivialPointerSampler() {} + void reset() + { + sampleVec.clear(); + } + public: std::vector sampleVec; @@ -149,6 +160,10 @@ public: MeshSampler(MeshType &_m):m(_m){} MeshType &m; + void reset() + { + m.Clear(); + } void AddVert(const VertexType &p) { @@ -180,18 +195,18 @@ The class is templated over the PointSampler object that allows to customize the **/ -template > +template > class SurfaceSampling { - typedef typename MetroMesh::CoordType CoordType; - typedef typename MetroMesh::ScalarType ScalarType; - typedef typename MetroMesh::VertexType VertexType; - typedef typename MetroMesh::VertexPointer VertexPointer; - typedef typename MetroMesh::VertexIterator VertexIterator; - typedef typename MetroMesh::FacePointer FacePointer; - typedef typename MetroMesh::FaceIterator FaceIterator; - typedef typename MetroMesh::FaceType FaceType; - typedef typename MetroMesh::FaceContainer FaceContainer; + typedef typename MeshType::CoordType CoordType; + typedef typename MeshType::ScalarType ScalarType; + typedef typename MeshType::VertexType VertexType; + typedef typename MeshType::VertexPointer VertexPointer; + typedef typename MeshType::VertexIterator VertexIterator; + typedef typename MeshType::FacePointer FacePointer; + typedef typename MeshType::FaceIterator FaceIterator; + typedef typename MeshType::FaceType FaceType; + typedef typename MeshType::FaceContainer FaceContainer; typedef typename vcg::Box3 BoxType; typedef typename vcg::SpatialHashTable MeshSHT; @@ -201,6 +216,8 @@ class SurfaceSampling typedef typename vcg::SpatialHashTable SampleSHT; typedef typename vcg::SpatialHashTable::CellIterator SampleSHTIterator; + typedef typename MeshType::template PerVertexAttributeHandle PerVertexFloatAttribute; + public: 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; for(vi=m.vert.begin();vi!=m.vert.end();++vi) @@ -358,7 +375,7 @@ static void AllVertex(MetroMesh & m, VertexSampler &ps) /// 2) shuffle vertices. /// 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; 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 /// 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; 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); } -static void FillAndShuffleFacePointerVector(MetroMesh & m, std::vector &faceVec) +static void FillAndShuffleFacePointerVector(MeshType & m, std::vector &faceVec) { - FaceIterator fi; - for(fi=m.face.begin();fi!=m.face.end();++fi) + for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) faceVec.push_back(&*fi); assert((int)faceVec.size()==m.fn); @@ -427,10 +443,9 @@ static void FillAndShuffleFacePointerVector(MetroMesh & m, std::vector &vertVec) +static void FillAndShuffleVertexPointerVector(MeshType & m, std::vector &vertVec) { - VertexIterator vi; - for(vi=m.vert.begin();vi!=m.vert.end();++vi) + for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD()) vertVec.push_back(&*vi); assert((int)vertVec.size()==m.vn); @@ -440,7 +455,7 @@ static void FillAndShuffleVertexPointerVector(MetroMesh & m, std::vector=m.vn) { 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. /// 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 angleSumH = tri::Allocator:: template GetPerVertexAttribute (m); + typename MeshType::template PerVertexAttributeHandle angleSumH = tri::Allocator:: template GetPerVertexAttribute (m); for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) angleSumH[vi]=0; @@ -480,7 +495,7 @@ static void VertexBorderCorner(MetroMesh & m, VertexSampler &ps, float angleRad) ps.AddVert(*vi); } - tri::Allocator:: template DeletePerVertexAttribute (m,angleSumH); + tri::Allocator:: template DeletePerVertexAttribute (m,angleSumH); } /// \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. /// 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::max()); } @@ -499,14 +514,14 @@ static void VertexBorder(MetroMesh & m, VertexSampler &ps) /// tri::UpdateFlags::FaceFauxCrease(mesh,creaseAngleRad); /// 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::PEdge SimpleEdge; + typedef typename UpdateTopology::PEdge SimpleEdge; std::vector< SimpleEdge > Edges; typename std::vector< SimpleEdge >::iterator ei; - UpdateTopology::FillUniqueEdgeVector(m,Edges,false); + UpdateTopology::FillUniqueEdgeVector(m,Edges,false); - typename MetroMesh::template PerVertexAttributeHandle hv = tri::Allocator:: template GetPerVertexAttribute (m); + typename MeshType::template PerVertexAttributeHandle hv = tri::Allocator:: template GetPerVertexAttribute (m); 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) { AllFace(m,ps); @@ -536,7 +551,7 @@ static void FaceUniform(MetroMesh & m, VertexSampler &ps, int sampleNum) ps.AddFace(*faceVec[i],Barycenter(*faceVec[i])); } -static void AllFace(MetroMesh & m, VertexSampler &ps) +static void AllFace(MeshType & m, VertexSampler &ps) { FaceIterator 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. - typedef typename UpdateTopology::PEdge SimpleEdge; + typedef typename UpdateTopology::PEdge SimpleEdge; std::vector< SimpleEdge > Edges; typename std::vector< SimpleEdge >::iterator ei; - UpdateTopology::FillUniqueEdgeVector(m,Edges); + UpdateTopology::FillUniqueEdgeVector(m,Edges); for(ei=Edges.begin(); ei!=Edges.end(); ++ei) 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 // 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::PEdge SimpleEdge; + typedef typename UpdateTopology::PEdge SimpleEdge; std::vector< SimpleEdge > Edges; - UpdateTopology::FillUniqueEdgeVector(m,Edges,sampleFauxEdge); + UpdateTopology::FillUniqueEdgeVector(m,Edges,sampleFauxEdge); // First loop compute total edge length; float edgeSum=0; 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]; } -static void StratifiedMontecarlo(MetroMesh & m, VertexSampler &ps,int sampleNum) +static void StratifiedMontecarlo(MeshType & m, VertexSampler &ps,int sampleNum) { - ScalarType area = Stat::ComputeMeshArea(m); + ScalarType area = Stat::ComputeMeshArea(m); ScalarType samplePerAreaUnit = sampleNum/area; // Montecarlo sampling. 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::ComputeMeshArea(m); + ScalarType area = Stat::ComputeMeshArea(m); ScalarType samplePerAreaUnit = sampleNum/area; FaceIterator fi; @@ -669,11 +684,11 @@ static void MontecarloPoisson(MetroMesh & m, VertexSampler &ps,int sampleNum) 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::PEdge SimpleEdge; + typedef typename UpdateTopology::PEdge SimpleEdge; std::vector< SimpleEdge > Edges; - UpdateTopology::FillUniqueEdgeVector(m,Edges,sampleAllEdges); + UpdateTopology::FillUniqueEdgeVector(m,Edges,sampleAllEdges); assert(!Edges.empty()); @@ -710,7 +725,7 @@ static void EdgeMontecarlo(MetroMesh & m, VertexSampler &ps, int sampleNum, bool 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 IntervalType; std::vector< IntervalType > intervals (m.fn+1); @@ -739,40 +754,45 @@ 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; - return DoubleArea(f)*averageQ/2.0; + ScalarType averageQ = ( wH[f.V(0)] + wH[f.V(1)] + wH[f.V(2)] )/3.0; + return averageQ*averageQ*DoubleArea(f)/2.0; } -/// Compute a sampling of the surface that is weighted by the quality -/// 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; -static void WeightedMontecarlo(MetroMesh & m, VertexSampler &ps, int sampleNum) +/// Compute a sampling of the surface that is weighted by the quality and a variance +/// +/// We use the quality as linear distortion of density. +/// 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:: template GetPerVertexAttribute (m,"radius"); + InitRadiusHandleFromQuality(m, rH, 1.0, variance, true); - ScalarType weightedArea = 0; - FaceIterator fi; - for(fi = m.face.begin(); fi != m.face.end(); ++fi) - if(!(*fi).IsD()) - weightedArea += WeightedArea(*fi); + ScalarType weightedArea = 0; + for(FaceIterator fi = m.face.begin(); fi != m.face.end(); ++fi) + weightedArea += WeightedArea(*fi,rH); - ScalarType samplePerAreaUnit = sampleNum/weightedArea; - // Montecarlo sampling. - double floatSampleNum = 0.0; - for(fi=m.face.begin(); fi != m.face.end(); fi++) - if(!(*fi).IsD()) + ScalarType samplePerAreaUnit = sampleNum/weightedArea; + // Montecarlo sampling. + double floatSampleNum = 0.0; + for(FaceIterator fi=m.face.begin(); fi != m.face.end(); fi++) { - // compute # samples in the current face (taking into account of the remainders) - floatSampleNum += WeightedArea(*fi) * samplePerAreaUnit; - int faceSampleNum = (int) floatSampleNum; + // compute # samples in the current face (taking into account of the remainders) + floatSampleNum += WeightedArea(*fi,rH) * samplePerAreaUnit; + int faceSampleNum = (int) floatSampleNum; - // for every sample p_i in T... - for(int i=0; i < faceSampleNum; i++) - ps.AddFace(*fi,RandomBarycentric()); + // for every sample p_i in T... + for(int i=0; i < faceSampleNum; i++) + ps.AddFace(*fi,RandomBarycentric()); - floatSampleNum -= (double) faceSampleNum; + floatSampleNum -= (double) faceSampleNum; } } @@ -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. -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::ComputeMeshArea(m); + ScalarType area = Stat::ComputeMeshArea(m); ScalarType samplePerAreaUnit = sampleNum/area; std::vector faceVec; FillAndShuffleFacePointerVector(m,faceVec); - vcg::tri::UpdateNormal::PerFaceNormalized(m); + vcg::tri::UpdateNormal::PerFaceNormalized(m); double floatSampleNum = 0.0; int faceSampleNum; // 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. -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::ComputeMeshArea(m); + ScalarType area = Stat::ComputeMeshArea(m); ScalarType samplePerAreaUnit = sampleNum/area; std::vector faceVec; FillAndShuffleFacePointerVector(m,faceVec); - tri::UpdateNormal::PerFaceNormalized(m); + tri::UpdateNormal::PerFaceNormalized(m); double floatSampleNum = 0.0; int faceSampleNum; // Subdivision sampling. @@ -1038,11 +1058,11 @@ static int SingleFaceSimilarDual(FacePointer fp, VertexSampler &ps, int n_sample -//template -//void Sampling::SimilarFaceSampling() -static void FaceSimilar(MetroMesh & m, VertexSampler &ps,int sampleNum, bool dualFlag, bool randomFlag) +//template +//void Sampling::SimilarFaceSampling() +static void FaceSimilar(MeshType & m, VertexSampler &ps,int sampleNum, bool dualFlag, bool randomFlag) { - ScalarType area = Stat::ComputeMeshArea(m); + ScalarType area = Stat::ComputeMeshArea(m); ScalarType samplePerAreaUnit = sampleNum/area; // 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 // 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). - static void SingleFaceRaster(typename MetroMesh::FaceType &f, VertexSampler &ps, - const Point2 & v0, - const Point2 & v1, - const Point2 & v2, + static void SingleFaceRaster(typename MeshType::FaceType &f, VertexSampler &ps, + const Point2 & v0, + const Point2 & v1, + const Point2 & v2, bool correctSafePointsBaryCoords=true) { - typedef typename MetroMesh::ScalarType S; + typedef typename MeshType::ScalarType S; // Calcolo bounding box Box2i bbox; Box2 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)) { - 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[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]; @@ -1192,7 +1212,7 @@ static void FaceSimilar(MetroMesh & m, VertexSampler &ps,int sampleNum, bool dua if (closeEdge >= 0) { - typename MetroMesh::CoordType baryCoord; + typename MeshType::CoordType baryCoord; if (correctSafePointsBaryCoords) { // Add x,y sample with closePoint barycentric coords (on edge) @@ -1224,7 +1244,7 @@ static bool checkPoissonDisk(SampleSHT & sht, const Point3 & p, Scal { // get the samples closest to the given one std::vector closests; - typedef VertTmark MarkerVert; + typedef VertTmark MarkerVert; static MarkerVert mv; Box3f bb(p-Point3f(radius,radius,radius),p+Point3f(radius,radius,radius)); @@ -1251,7 +1271,6 @@ struct PoissonDiskParam preGenFlag = false; preGenMesh = NULL; geodesicDistanceFlag = false; - pds=NULL; } struct Stat @@ -1272,11 +1291,13 @@ struct PoissonDiskParam bool adaptiveRadiusFlag; float radiusVariance; 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. - MetroMesh *preGenMesh; + 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//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; - 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::ComputeMeshArea(origMesh); + ScalarType meshArea = Stat::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) @@ -1330,29 +1351,31 @@ static ScalarType ComputePoissonDiskRadius(MetroMesh &origMesh, int sampleNum) return diskRadius; } -static int ComputePoissonSampleNum(MetroMesh &origMesh, ScalarType diskRadius) +static int ComputePoissonSampleNum(MeshType &origMesh, ScalarType diskRadius) { - ScalarType meshArea = Stat::ComputeMeshArea(origMesh); + ScalarType meshArea = Stat::ComputeMeshArea(origMesh); int sampleNum = meshArea / (diskRadius*diskRadius *M_PI *0.7) ; // 0.7 is a density factor return sampleNum; } -// Note that this function actually CHANGE the quality of the montecarlo samples so that it represents the expected radius of each sample -// the expected radius of the sample is computed so that it linearly maps the quality between diskradius / variance and diskradius*variance +/// When performing an adptive pruning for each sample we expect a varying radius to be removed. +/// 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 minmax = tri::Stat::ComputePerVertexQualityMinMax( sampleMesh); - float minRad = diskRadius / radiusVariance; - float maxRad = diskRadius * radiusVariance; - float deltaQ = minmax.second-minmax.first; - float deltaRad = maxRad-minRad; - for (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); - } + std::pair minmax = tri::Stat::ComputePerVertexQualityMinMax( sampleMesh); + float minRad = diskRadius ; + float maxRad = diskRadius * radiusVariance; + float deltaQ = minmax.second-minmax.first; + float deltaRad = maxRad-minRad; + for (VertexIterator vi = sampleMesh.vert.begin(); vi != sampleMesh.vert.end(); vi++) + { + rH[*vi] = minRad + deltaRad*((invert ? minmax.second - (*vi).Q() : (*vi).Q() - minmax.first )/deltaQ); + } } // initialize spatial hash table for searching @@ -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. // 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, - const struct PoissonDiskParam pp=PoissonDiskParam()) +static void InitSpatialHashTable(MeshType &montecarloMesh, MontecarloSHT &montecarloSHT, ScalarType diskRadius, + struct PoissonDiskParam pp=PoissonDiskParam()) { ScalarType cellsize = 2.0f* diskRadius / sqrt(3.0); float occupancyRatio=0; @@ -1380,56 +1403,126 @@ static void InitSpatialHashTable(MetroMesh &montecarloMesh, MontecarloSHT &monte montecarloSHT.InitEmpty(bb, gridsize); for (VertexIterator vi = montecarloMesh.vert.begin(); vi != montecarloMesh.vert.end(); vi++) - if(!(*vi).IsD()) - { - montecarloSHT.Add(&(*vi)); - } + if(!(*vi).IsD()) + { + montecarloSHT.Add(&(*vi)); + } 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; 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); } 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 PoissonDiskPruning(VertexSampler &ps, MetroMesh &montecarloMesh, - ScalarType diskRadius, const struct PoissonDiskParam pp=PoissonDiskParam()) +static void PoissonDiskPruningByNumber(VertexSampler &ps, MeshType &m, + int sampleNum, ScalarType &diskRadius, + PoissonDiskParam &pp, + float tolerance=0.04, + int maxIter=20) + { - if(pp.adaptiveRadiusFlag) - tri::RequirePerVertexQuality(montecarloMesh); - int t0 = clock(); + 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 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) + tri::RequirePerVertexQuality(montecarloMesh); + int t0 = clock(); // spatial index of montecarlo samples - used to choose a new sample to insert MontecarloSHT montecarloSHT; 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 + PerVertexFloatAttribute rH = tri::Allocator:: template GetPerVertexAttribute (montecarloMesh,"radius"); 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; std::random_shuffle(montecarloSHT.AllocatedCells.begin(),montecarloSHT.AllocatedCells.end(), p_myrandom); int t1 = clock(); - if(pp.pds) { - pp.pds->montecarloSampleNum = montecarloMesh.vn; - } + pp.pds.montecarloSampleNum = montecarloMesh.vn; + pp.pds.sampleNum =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 + if(pp.preGenFlag) { - // 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(!(*vi).IsD()) - { + if(pp.preGenMesh==0) + { + typename MeshType::template PerVertexAttributeHandle fixed; + fixed = tri::Allocator:: template GetPerVertexAttribute (montecarloMesh,"fixed"); + for(VertexIterator vi=montecarloMesh.vert.begin();vi!=montecarloMesh.vert.end();++vi) + if(fixed[*vi]) { + pp.pds.sampleNum++; ps.AddVert(*vi); 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(); } vertex::ApproximateGeodesicDistanceFunctor GDF; @@ -1447,20 +1540,18 @@ static void PoissonDiskPruning(VertexSampler &ps, MetroMesh &montecarloMesh, sp = getSampleFromCell(montecarloSHT.AllocatedCells[i], montecarloSHT); if(pp.adaptiveRadiusFlag) - currentRadius = sp->Q(); + currentRadius = rH[sp]; ps.AddVert(*sp); + pp.pds.sampleNum++; if(pp.geodesicDistanceFlag) removedCnt += montecarloSHT.RemoveInSphereNormal(sp->cP(),sp->cN(),GDF,currentRadius); else removedCnt += montecarloSHT.RemoveInSphere(sp->cP(),currentRadius); } montecarloSHT.UpdateAllocatedCells(); } 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. @@ -1473,7 +1564,7 @@ static void PoissonDiskPruning(VertexSampler &ps, MetroMesh &montecarloMesh, * IEEE Symposium on Interactive Ray Tracing, 2007, * 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(); // 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(); // if we are doing variable density sampling we have to prepare the random samples quality with the correct expected radii. + PerVertexFloatAttribute rH = tri::Allocator:: template GetPerVertexAttribute (montecarloMesh,"radius"); if(pp.adaptiveRadiusFlag) - ComputePoissonSampleRadii(montecarloMesh, diskRadius, pp.radiusVariance, pp.invertQuality); + InitRadiusHandleFromQuality(montecarloMesh, rH, diskRadius, pp.radiusVariance, pp.invertQuality); do { @@ -1556,7 +1648,7 @@ static void HierarchicalPoissonDisk(MetroMesh &origMesh, VertexSampler &ps, Metr VertexPointer sp = (*hi).second; // vr spans between 3.0*r and r / 4.0 according to vertex quality ScalarType sampleRadius = diskRadius; - if(pp.adaptiveRadiusFlag) sampleRadius = sp->Q(); + if(pp.adaptiveRadiusFlag) sampleRadius = rH[sp]; if (checkPoissonDisk(checkSHT, sp->cP(), sampleRadius)) { ps.AddVert(*sp); @@ -1580,8 +1672,8 @@ static void HierarchicalPoissonDisk(MetroMesh &origMesh, VertexSampler &ps, Metr } while(level < 5); } -//template -//void Sampling::SimilarFaceSampling() +//template +//void Sampling::SimilarFaceSampling() // 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) @@ -1593,7 +1685,7 @@ static void HierarchicalPoissonDisk(MetroMesh &origMesh, VertexSampler &ps, Metr // Else make sure to update border flags from texture space FFadj // vcg::tri::UpdateTopology::FaceFaceFromTexCoord(m); // vcg::tri::UpdateFlags::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; @@ -1617,11 +1709,11 @@ class RRParam public: float offset; float minDiag; -tri::FaceTmark markerFunctor; +tri::FaceTmark markerFunctor; TriMeshGrid gM; }; -static void RegularRecursiveOffset(MetroMesh & m, std::vector &pvec, ScalarType offset, float minDiag) +static void RegularRecursiveOffset(MeshType & m, std::vector &pvec, ScalarType offset, float minDiag) { Box3 bb=m.bbox; bb.Offset(offset*2.0); @@ -1638,7 +1730,7 @@ static void RegularRecursiveOffset(MetroMesh & m, std::vector &pvec, Sc SubdivideAndSample(m, pvec, bb, rrp, bb.Diag()); } -static void SubdivideAndSample(MetroMesh & m, std::vector &pvec, const Box3 bb, RRParam &rrp, float curDiag) +static void SubdivideAndSample(MeshType & m, std::vector &pvec, const Box3 bb, RRParam &rrp, float curDiag) { Point3f startPt = bb.Center(); @@ -1679,6 +1771,22 @@ static void SubdivideAndSample(MetroMesh & m, std::vector &pvec, const } }; // end sampling class +template +typename MeshType::ScalarType ComputePoissonDiskRadius(MeshType &origMesh, int sampleNum) +{ + typedef typename MeshType::ScalarType ScalarType; + ScalarType meshArea = Stat::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 BaseSampler; typedef tri::MeshSampler MontecarloSampler; - tri::RequireVFAdjacency(m); typename tri::SurfaceSampling::PoissonDiskParam pp; - typename tri::SurfaceSampling::PoissonDiskParam::Stat stat; - pp.pds = &stat; int t0=clock(); if(sampleNum>0) radius = tri::SurfaceSampling::ComputePoissonDiskRadius(m,sampleNum); if(radius>0 && sampleNum==0) sampleNum = tri::SurfaceSampling::ComputePoissonSampleNum(m,radius); - pp.pds->sampleNum = sampleNum; + pp.pds.sampleNum = sampleNum; poissonSamples.clear(); // std::vector MontecarloSamples; MeshType MontecarloMesh; @@ -1738,7 +1843,7 @@ void PoissonSampling(MeshType &m, // the mesh that has to be sampled tri::UpdateBounding::Box(MontecarloMesh); // tri::Build(MontecarloMesh, MontecarloSamples); int t1=clock(); - pp.pds->montecarloTime = t1-t0; + pp.pds.montecarloTime = t1-t0; if(radiusVariance !=1) { @@ -1747,42 +1852,22 @@ void PoissonSampling(MeshType &m, // the mesh that has to be sampled } tri::SurfaceSampling::PoissonDiskPruning(pdSampler, MontecarloMesh, radius,pp); int t2=clock(); - 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 -void PoissonPruning(MeshType &m, // the mesh that has to be pruned - std::vector &poissonSamples, // the vector that will contain the chosen set of points - float & radius, int sampleNum=0) -{ - typedef tri::TrivialSampler BaseSampler; - typename tri::SurfaceSampling::PoissonDiskParam pp; - if(sampleNum>0 && radius==0) - radius = tri::SurfaceSampling::ComputePoissonDiskRadius(m,sampleNum); - tri::UpdateBounding::Box(m); - BaseSampler pdSampler(poissonSamples); - tri::SurfaceSampling::PoissonDiskPruning(pdSampler, m, radius,pp); + pp.pds.totalTime = t2-t0; } - -/// \brief Very simple wrapping for the Poisson Disk Pruning +/// \brief Low level wrapper for 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 void PoissonPruning(MeshType &m, // the mesh that has to be pruned std::vector &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 BaseSampler; typename tri::SurfaceSampling::PoissonDiskParam pp; if(randSeed !=0) tri::SurfaceSampling::SamplingRandomGenerator().initialize(randSeed); - if(sampleNum>0 && radius==0) - radius = tri::SurfaceSampling::ComputePoissonDiskRadius(m,sampleNum); tri::UpdateBounding::Box(m); BaseSampler pdSampler; @@ -1790,38 +1875,73 @@ void PoissonPruning(MeshType &m, // the mesh that has to be pruned 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 +void PoissonPruning(MeshType &m, // the mesh that has to be pruned + std::vector &poissonSamples, // the vector that will contain the chosen set of points + float radius, unsigned int randSeed=0) +{ + std::vector poissonSamplesVP; + PoissonPruning(m,poissonSamplesVP,radius,randSeed); + poissonSamples.resize(poissonSamplesVP.size()); + for(size_t i=0;iP(); +} + + + /// \brief Very simple wrapping for the Exact Poisson Disk Pruning /// /// 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. /// It is obviously much slower than the other versions... - - template -void PoissonPruningExact(MeshType &m, // the mesh that has to be pruned - std::vector &poissonSamples, // the vector that will contain the chosen set of points +void PoissonPruningExact(MeshType &m, /// the mesh that has to be pruned + std::vector &poissonSamples, /// the vector that will contain the chosen set of points float & radius, int sampleNum, float tolerance=0.04, int maxIter=20) { - int sampleNumMin = int(float(sampleNum)*(1.0f-tolerance)); - int sampleNumMax = int(float(sampleNum)*(1.0f+tolerance)); - float RadiusRangeMin = m.bbox.Diag()/100.0; - float RadiusRangeMax = m.bbox.Diag(); - std::vector poissonSamplesTmp; + size_t sampleNumMin = int(float(sampleNum)*(1.0f-tolerance)); + size_t sampleNumMax = int(float(sampleNum)*(1.0f+tolerance)); + float RangeMinRad = m.bbox.Diag()/10.0f; + float RangeMaxRad = m.bbox.Diag()/10.0f; + size_t RangeMinNum; + size_t RangeMaxNum; + std::vector 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; int iterCnt=0; while(iterCnt sampleNumMax) ) { - curRadius=(RadiusRangeMax+RadiusRangeMin)/2.0f; + curRadius=(RangeMaxRad+RangeMinRad)/2.0f; 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) - RadiusRangeMin = curRadius; + RangeMinRad = curRadius; if(poissonSamplesTmp.size() < sampleNum) - RadiusRangeMax = curRadius; + RangeMaxRad = curRadius; } swap(poissonSamples,poissonSamplesTmp);