From 7bd8b4f19b54daeacc2135b75d12d55d9f427b24 Mon Sep 17 00:00:00 2001 From: maxcorsini Date: Fri, 16 Jan 2009 11:30:19 +0000 Subject: [PATCH] poissondisk sampling completely restructure for performance --- vcg/complex/trimesh/point_sampling.h | 248 +++++++++------------------ 1 file changed, 83 insertions(+), 165 deletions(-) diff --git a/vcg/complex/trimesh/point_sampling.h b/vcg/complex/trimesh/point_sampling.h index 4ac739f9..a5e92fcf 100644 --- a/vcg/complex/trimesh/point_sampling.h +++ b/vcg/complex/trimesh/point_sampling.h @@ -88,6 +88,13 @@ class SurfaceSampling typedef typename MetroMesh::FaceType FaceType; typedef typename MetroMesh::FaceContainer FaceContainer; + typedef typename vcg::SpatialHashTable MeshSHT; + typedef typename vcg::SpatialHashTable::CellIterator MeshSHTIterator; + typedef typename vcg::SpatialHashTable MontecarloSHT; + typedef typename vcg::SpatialHashTable::CellIterator MontecarloSHTIterator; + typedef typename vcg::SpatialHashTable SampleSHT; + typedef typename vcg::SpatialHashTable::CellIterator SampleSHTIterator; + public: static math::MarsenneTwisterRNG &SamplingRandomGenerator() @@ -658,105 +665,33 @@ static CoordType RandomBox(vcg::Box3 box) return p; } -// naive projection generates a sample inside the given box, and project it in the -// faces intersected by this box. It always returns a valid sample. -static vcg::Point3 naiveProjection(vcg::Box3 box, - MetroMesh &mesh, - vcg::SpatialHashTable & sht) +static bool generatePoissonDiskSample(Point3i *cell, MontecarloSHT samplepool, vcg::Point3 & p) { - vcg::Point3 p; + MontecarloSHTIterator cellBegin; + MontecarloSHTIterator cellEnd; + MontecarloSHTIterator cellIt; + samplepool.Grid(*cell, cellBegin, cellEnd); - p = RandomBox(box); + std::vector samples; + for (cellIt = cellBegin; cellIt != cellEnd; cellIt++) + samples.push_back(*cellIt); - // projection on the faces - CoordType closestp = p; - ScalarType maxdist = box.DimX(); - ScalarType mindist; - - // THIS NOT WORK... - //vcg::tri::GetClosestFace >( - // mesh, sht, p, maxdist, mindist, closestp); - -/* THIS NOT WORK... - typedef vcg::tri::FaceTmark Marker; - Marker markerFunctor; - markerFunctor.SetMesh(&mesh); - FaceType *nearestF=0; - vcg::face::PointDistanceBaseFunctor PDistFunct; - nearestF = sht.GetClosest(PDistFunct,markerFunctor,p,maxdist,mindist,closestp); -*/ - return closestp; -} - -// Montecarlo "reduced" generates a sample on the faces intersected by the given box -static vcg::Point3 naiveMontecarloReduced(vcg::Box3 box, - std::vector faces) -{ - vcg::Point3 p; - //...TODO... - return p; -} - -// naive projection generates a sample inside the given box, and project it on the -// faces intersected by this box. If the projected point lies outside the box -// the sample is re-generated. -// false is returned if the maximum number of attempts has been reached. -static bool naiveProjectionInside(vcg::Box3 box, - MetroMesh &mesh, - vcg::SpatialHashTable & sht, - int maxattemps, vcg::Point3 &p) -{ - bool flaginside = false; - int k=1; - do + if (samples.size() >= 1) { - p = naiveProjection(box, mesh, sht); - - // check if p is inside the box - flaginside = box.IsIn(p); - k++; - - } while (!flaginside && k < maxattemps); - - return flaginside; + int index = RandomInt(samples.size()); + p = samples[index]->P(); + return true; + } + else + return false; } -// naive surface generation generates a sample on the faces intersected by the given box -// if the point lies outside the box the sample is re-generated. -// false is returned if the maximum number of attempts has been reached. -static bool naiveMontecarloReducedInside(vcg::Box3 box, std::vector faces, - int maxattemps, vcg::Point3 &p) -{ - //...TODO... - return true; -} - -// [erfect projection generates a sample inside the given box, and project it -// on the surface inside this box. It always returns a valid sample. -static vcg::Point3 perfectProjection(vcg::Box3 box, std::vector faces) -{ - vcg::Point3 p; - //...TODO... - return p; -} - -// perfect Montecarlo Reduced generates a sample on the faces contained in the given box. -// Intersecting faces are treated properly. It always returns a valid sample. -static vcg::Point3 perfectMontecarloReduced(vcg::Box3 box, std::vector faces) -{ - vcg::Point3 p; - //...TODO... - return p; -} - - - // check the radius constrain -static bool checkPoissonDisk(vcg::SpatialHashTable sht, Point3 p, ScalarType radius) +static bool checkPoissonDisk(SampleSHT sht, Point3 p, ScalarType radius) { - typename vcg::SpatialHashTable::CellIterator itBegin; - typename vcg::SpatialHashTable::CellIterator itEnd; - typename vcg::SpatialHashTable::CellIterator it; + SampleSHTIterator itBegin; + SampleSHTIterator itEnd; + SampleSHTIterator it; // get the samples closest to the given one sht.Grid(p, itBegin, itEnd); @@ -784,18 +719,18 @@ static bool checkPoissonDisk(vcg::SpatialHashTable sht, * IEEE Symposium on Interactive Ray Tracing, 2007, * 10-12 Sept. 2007, pp. 129-132. */ -static void Poissondisk(MetroMesh &m, VertexSampler &ps, int sampleNum, int version) +static void Poissondisk(MetroMesh &m, VertexSampler &ps, MetroMesh & montecarloMesh, int sampleNum) { const int MAXLEVELS = 10; // maximum level of subdivision // spatial index of mesh face - used to search where to place the samples - vcg::SpatialHashTable searchSHT; - typename vcg::SpatialHashTable::CellIterator cellBegin; - typename vcg::SpatialHashTable::CellIterator cellEnd; - typename vcg::SpatialHashTable::CellIterator cellIt; + //MeshSHT searchSHT; + + // spatial index of montecarlo samples - used to choose a new sample to insert + MontecarloSHT montecarloSHT; // spatial hash table of the generated samples - used to check the radius constrain - vcg::SpatialHashTable checkSHT; + SampleSHT checkSHT; // initialize spatial hash table for searching ScalarType meshArea = Stat::ComputeMeshArea(m); @@ -806,37 +741,37 @@ static void Poissondisk(MetroMesh &m, VertexSampler &ps, int sampleNum, int vers int sizeX = m.bbox.DimX() / r; int sizeY = m.bbox.DimY() / r; int sizeZ = m.bbox.DimZ() / r; - Point3i gridsize(sizeX, sizeY, sizeZ); +/* searchSHT.InitEmpty(m.bbox, gridsize); FaceIterator fi; for (fi = m.face.begin(); fi != m.face.end(); fi++) searchSHT.Add(&(*fi)); +*/ + // initialize spatial hash to index pre-generated samples + VertexIterator vi; + montecarloSHT.InitEmpty(m.bbox, gridsize); + for (vi = m.vert.begin(); vi != m.vert.end(); vi++) + montecarloSHT.Add(&(*vi)); - // initialize spatial hash table for checking + // initialize spatial hash table for check poisson-disk radius constrain checkSHT.InitEmpty(m.bbox, gridsize); - // sampling algorithm ("Projection-based") - // --------------------------------------------------- + // sampling algorithm + // ------------------ // + // - generate millions of samples using montecarlo algorithm // - extract a cell (C) from the active cell list (with probability proportional to cell's volume) // - with a probability proportional to the intersection between the surface and the cell - // generate a sample inside C and project it on the mesh + // generate a sample inside C by choosing one of the contained pre-generated samples // - if the sample violates the radius constrain discard it, and add the cell to the cells-to-subdivide list // - iterate until the active cell list is empty or a pre-defined number of subdivisions is reached // - // sampling algorithm ("Surface-based") - // ------------------------------------------------ - // - // - extract a cell (C) from the active cell list (with probability proportional to the cell's volume) - // - with a probability proportional to the intersection between the surface and the cell - // generate a sample on the triangles inside C - // - if the sample violated the radius constrain discard it, and add the cell to the cells-to-subdivide list - // - iterate until the active cell list is empty or a pre-defined number of subdivisions is reached - // - std::vector activeCells; + std::vector cellsToSubdivide; + typename std::vector::iterator cellIt; + std::set nextFaces; // faces to add to the next level of subdivision typename std::set::iterator nextFacesIt; typename std::vector::iterator it; @@ -854,7 +789,7 @@ static void Poissondisk(MetroMesh &m, VertexSampler &ps, int sampleNum, int vers /////////////////////////////////////////////////////////////////////////////////////////////////// // create active cell list - for (it = searchSHT.AllocatedCells.begin(); it != searchSHT.AllocatedCells.end(); it++) + for (it = montecarloSHT.AllocatedCells.begin(); it != montecarloSHT.AllocatedCells.end(); it++) { activeCells.push_back(&(*it)); } @@ -865,7 +800,7 @@ static void Poissondisk(MetroMesh &m, VertexSampler &ps, int sampleNum, int vers int ncell = static_cast(activeCells.size()); int index,index2; Point3i *temp; - for (int i = 0; i < 129248; i++) + for (int i = 0; i < ncell/2; i++) { index = RandomInt(ncell); index2 = RandomInt(ncell); @@ -884,47 +819,8 @@ static void Poissondisk(MetroMesh &m, VertexSampler &ps, int sampleNum, int vers { currentCell = activeCells[i]; - // calculate box and contained faces - searchSHT.Grid(*currentCell, cellBegin, cellEnd); - searchSHT.IPiToBox(*currentCell, currentBox); - - faces.clear(); - for (cellIt = cellBegin; cellIt != cellEnd; ++cellIt) - { - faces.push_back(*cellIt); - } - - validsample = true; - if (version == 0) - { - // naive projection method (unsafe) - s = naiveProjection(currentBox, m, searchSHT); - } - else if (version == 1) - { - // naive montecarlo reduced (unsafe) - s = naiveMontecarloReduced(currentBox, faces); - } - else if (version == 2) - { - // naive projection method (safer, but still not perfect) - validsample = naiveProjectionInside(currentBox, m, searchSHT, 10, s); - } - else if (version == 3) - { - // naive montecarlo reduced (safer, but still not perfect) - validsample = naiveMontecarloReducedInside(currentBox, faces, 10, s); - } - else if (version == 4) - { - // perfect projection method - s = perfectProjection(currentBox, faces); - } - else if (version == 5) - { - // perfect montecarlo reduced - s = perfectMontecarloReduced(currentBox, faces); - } + // generate a sample chosen from the pre-generated one + validsample = generatePoissonDiskSample(currentCell, montecarloSHT, s); if (validsample) { @@ -946,30 +842,52 @@ static void Poissondisk(MetroMesh &m, VertexSampler &ps, int sampleNum, int vers } else { - // add these faces to the faces for the next level of subdivision - for (facesIt = faces.begin(); facesIt != faces.end(); facesIt++) - nextFaces.insert(*facesIt); + // subdivide this cell + cellsToSubdivide.push_back(currentCell); } } - activeCells.clear(); // proceed to the next level of subdivision - searchSHT.Clear(); + /////////////////////////////////////////////////////////////////////////// + + // cleaning spatial index data structures + //searchSHT.Clear(); + montecarloSHT.Clear(); + + // increase grid resolution gridsize[0] *= 2; gridsize[1] *= 2; gridsize[2] *= 2; - searchSHT.InitEmpty(m.bbox, gridsize); - for (nextFacesIt = nextFaces.begin(); nextFacesIt != nextFaces.end(); nextFacesIt++) - searchSHT.Add(*nextFacesIt); + //searchSHT.InitEmpty(m.bbox, gridsize); + montecarloSHT.InitEmpty(m.bbox, gridsize); + for (cellIt = cellsToSubdivide.begin(); cellIt != cellsToSubdivide.end(); cellIt++) + { + MeshSHTIterator faceBegin,faceEnd,faceIt; + MontecarloSHTIterator ptBegin,ptEnd,ptIt; + + Point3i *cell = *cellIt; + + // pre-generated samples for the next level of subdivision + montecarloSHT.Grid(*cell, ptBegin, ptEnd); + for (ptIt = ptBegin; ptIt != ptEnd; ptIt++) + montecarloSHT.Add(*ptIt); + +/* + // faces of the mesh for the next level of subdivision + searchSHT.Grid(*cell, faceBegin, faceEnd); + for (faceIt = faceBegin; faceIt != faceEnd; faceIt++) + searchSHT.Add(*faceIt); +*/ + } + level++; - } while(level < 2); + } while(level < 1); - } //template