Merge pull request #186 from dyollb/fix_poisson_disk_sampler

fix bug in PoissonDiskSampler, add test

docs/docs/surface/algorithms/surface_sampling.md

@@ -8,9 +8,7 @@ Currently the algorithm only works on manifold meshes.
 
 The algorithm has a few parameters that roughly correspond to the algorithm of Bridson's 2007 [Fast Poisson Disk Sampling in Arbitrary Dimensions](https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf).
 
-Additionally, you can specify points around samples should be avoided (shown in red below.) By default, samples will avoid these points with the same radius `r` used in the rest of the algorithm. You can optionally specify a "radius of avoidance" for these points, where the radius of avoidance is given in multiples of `r`. 
-
-The radius of avoidance can be further be specified to be a radius in 3D space, or a radius in terms of distance along the surface. The former will produce a radius of avoidance that will appear perfectly round and is likely more visually pleasing, but for very large radii may occlude samples from opposite sides of the mesh. The latter will restrict the radius of avoidance to only be along the surface, but such a metric ball will not appear perfectly round, especially in areas with very large and sudden changes in curvature.
+Additionally, you can specify points around samples should be avoided (shown in red below.)
 
 ![poisson disk sample with point of avoidance](/media/poisson_disk_sample.png)
 
@@ -23,15 +21,9 @@ The radius of avoidance can be further be specified to be a radius in 3D space,
 
     The mesh and geometry cannot be changed after construction.
 
-??? func "`#!cpp std::vector<SurfacePoint> PoissonDiskSampler::sample(double rCoef = 1.0, int kCandidates = 30, std::vector<SurfacePoint> pointsToAvoid = std::vector<SurfacePoint>(), int rAvoidance = 1, bool use3DAvoidanceRadius = true);`"
+??? func "`#!cpp std::vector<SurfacePoint> sample(const PoissonDiskOptions& options = PoissonDiskOptions())`"
 
     Poisson disk-sample the surface mesh.
-
-    - `rCoef`: corresponds to the minimum distance between samples, expressed as a multiple of the mean edge length. The actual minimum distance is computed as `r = rCoef * meanEdgeLength`
-    - `kCandidates`: the number of candidate points chosen from the (r,2r)-annulus around each sample.
-    - `pointsToAvoid`: SurfacePoints which samples should avoid.
-    - `rAvoidance`: the radius of avoidance around each point to avoid, expressed as a multiple of `r`.
-    - `use3DAvoidanceRadius`: If true, the radius of avoidance will specify a solid ball in 3D space around which samples are avoided. Otherwise, samples are avoided within a ball _on the surface_.
 
 ### Example
 
@@ -51,4 +43,23 @@ std::tie(mesh, geometry) = readManifoldSurfaceMesh(filename);
 // construct a solver
 PoissonDiskSampler poissonSampler(*mesh, *geometry);
 std::vector<SurfacePoint> samples = poissonSampler.sample(); // sample using default parameters
-```
+
+// Sample with some different parameters.
+PoissonDiskOptions sampleOptions;
+sampleOptions.minDist = 2.;
+std::vector<SurfacePoint> newSamples = poissonSampler.sample(sampleOptions);
+```
+## Helper Types
+### Options
+
+Options are passed in to `options` via a `PoissonDiskOptions` object.
+
+| Field | Default value |Meaning|
+|---|---|---|
+| `#!cpp double minDist`| `1` | The minimum distance `r` between samples, expressed in world-space units. |
+| `#!cpp int kCandidates`| `30` | The number of candidate points chosen from the (`r`,2`r`)-annulus around each sample. |
+| `#!cpp std::vector<SurfacePoint> pointsToAvoid`| `std::vector<SurfacePoint>()` | Points which samples should avoid. |
+| `#!cpp double minDistAvoidance`| `1` | The radius of avoidance around each point to avoid, expressed in world-space units. |
+| `#!cpp bool use3DAvoidance`| `true` | If true, the radius of avoidance will specify a solid ball in 3D space around which samples are avoided. Otherwise, samples are avoided within a geodesic ball on the surface. |
+
+Using the `use3DAvoidance` option, the radius of avoidance `minDistAvoidance` can specify either the radius in 3D space, or in terms of distance along the surface. The former will produce a radius of avoidance that will appear perfectly round, and is likely more visually pleasing, but for very large radii may occlude samples from opposite sides of the mesh. The latter will restrict the radius of avoidance to only be along the surface, but such a metric ball may not appear perfectly round, especially in areas with very large changes in curvature.

include/geometrycentral/surface/poisson_disk_sampler.h

@@ -7,6 +7,14 @@
 namespace geometrycentral {
 namespace surface {
 
+struct PoissonDiskOptions {
+  double minDist = 1.0; // minimum distance r between samples
+  int kCandidates = 30; // number of candidate points chosen from the (r,2r)-annulus around each sample
+  std::vector<SurfacePoint> pointsToAvoid;
+  double minDistAvoidance = 1.0; // radius of avoidance
+  bool use3DAvoidance = true;
+};
+
 class PoissonDiskSampler {
 
 public:
@@ -18,22 +26,18 @@ class PoissonDiskSampler {
   PoissonDiskSampler(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geometry);
 
   // ===== The main function: Sample the surface using Poisson Disk Sampling.
-  std::vector<SurfacePoint> sample(double rCoef = 1.0, int kCandidates = 30,
-                                   std::vector<SurfacePoint> pointsToAvoid = std::vector<SurfacePoint>(),
-                                   int rAvoidance = 1, bool use3DAvoidanceRadius = true);
+  std::vector<SurfacePoint> sample(const PoissonDiskOptions& options = PoissonDiskOptions());
 
 private:
   // ===== Members
   ManifoldSurfaceMesh& mesh;
   VertexPositionGeometry& geometry;
 
-  double rCoef;    // minimum distance between samples, expressed as a multiple of the mean edge length
   int kCandidates; // number of candidate points chosen from the (r,2r)-annulus around each sample
   std::vector<SurfacePoint> pointsToAvoid;
 
-  double rMinDist;       // the actual minimum distance between samples
-  double meanEdgeLength; // the mean edge length
-  double sideLength;     // side length of each bucket
+  double rMinDist;   // the minimum distance between samples
+  double sideLength; // side length of each bucket
 
   std::vector<Face> faceFromEachComponent; // holds one face for each connected component in the mesh
   std::vector<SurfacePoint> activeList;    // holds candidate points
@@ -51,7 +55,7 @@ class PoissonDiskSampler {
   void addNewSample(const SurfacePoint& sample);
   SpatialKey positionKey(const Vector3& position) const;
   void addPointToSpatialLookup(const Vector3& newPos);
-  void addPointToSpatialLookupWithRadius(const SurfacePoint& newPoint, int R = 0, bool use3DAvoidanceRadius = true);
+  void addPointToSpatialLookupWithRadius(const SurfacePoint& newPoint, double radius = 0, bool use3DAvoidance = true);
   bool isCandidateValid(const SurfacePoint& candidate) const;
   void sampleOnConnectedComponent(const Face& f);
   void clearData();

src/surface/poisson_disk_sampler.cpp

@@ -6,15 +6,6 @@ namespace surface {
 PoissonDiskSampler::PoissonDiskSampler(ManifoldSurfaceMesh& mesh_, VertexPositionGeometry& geometry_)
     : mesh(mesh_), geometry(geometry_) {
 
-  // Compute mean edge length.
-  geometry.requireEdgeLengths();
-  meanEdgeLength = 0.;
-  for (Edge e : mesh.edges()) {
-    meanEdgeLength += geometry.edgeLengths[e];
-  }
-  meanEdgeLength /= mesh.nEdges();
-  geometry.unrequireEdgeLengths();
-
   // Prepare spatial lookup
   Vector3 bboxMin, bboxMax;
   std::tie(bboxMin, bboxMax) = boundingBox();
@@ -90,7 +81,7 @@ SurfacePoint PoissonDiskSampler::generateCandidate(const SurfacePoint& xi) const
   TraceGeodesicResult trace;
 
   Vector2 dir = Vector2::fromAngle(randomReal(0., 2. * PI));
-  double dist = std::sqrt(randomReal(rMinDist, 2. * rMinDist));
+  double dist = randomReal(rMinDist, 2. * rMinDist);
   trace = traceGeodesic(geometry, xi, dist * dir);
 
   pathEndpoint = trace.endPoint;
@@ -124,17 +115,15 @@ void PoissonDiskSampler::addPointToSpatialLookup(const Vector3& newPos) {
 }
 
 /*
- * Mark all buckets with a radius of <R> buckets as occupied, as well.
+ * Mark all buckets within a radius of <radius> as occupied.
  */
-void PoissonDiskSampler::addPointToSpatialLookupWithRadius(const SurfacePoint& newPoint, int R,
-                                                           bool use3DAvoidanceRadius) {
+void PoissonDiskSampler::addPointToSpatialLookupWithRadius(const SurfacePoint& newPoint, double radius,
+                                                           bool use3DAvoidance) {
 
   Vector3 newPos = newPoint.interpolate(geometry.vertexPositions);
   addPointToSpatialLookup(newPos);
 
-  if (R <= 0) return;
-
-  if (!use3DAvoidanceRadius) {
+  if (!use3DAvoidance) {
     // This places fictitious points in a metric ball approximately of radius R*r centered at <newPoint>.
     // The solid ball is built by constructing R layers, radially outward; each layer is spaced r apart, and
     // points in each layer are spaced approx. r apart around the circular layer. The idea is that no point can be added
@@ -144,21 +133,22 @@ void PoissonDiskSampler::addPointToSpatialLookupWithRadius(const SurfacePoint& n
     // "spiky" meshes.
     SurfacePoint pathEndpoint;
     TraceGeodesicResult trace;
-    for (int rIter = 0; rIter <= R; rIter++) {
-      double dist = rIter * rMinDist;
-      for (double theta = 0.; theta <= 2. * PI; theta += rMinDist / dist / 2.) {
-        trace = traceGeodesic(geometry, newPoint, dist * Vector2::fromAngle(theta));
+    double r = rMinDist;
+    while (r < radius) {
+      for (double theta = 0.; theta <= 2. * PI; theta += rMinDist / r / 2.) {
+        trace = traceGeodesic(geometry, newPoint, r * Vector2::fromAngle(theta));
         pathEndpoint = trace.endPoint;
         addPointToSpatialLookup(pathEndpoint.interpolate(geometry.vertexPositions));
       }
+      r += rMinDist;
     }
   } else {
     // This places fictitious points in a *solid 3D ball* approximately of radius R*r centered at <newPoint>.
     // The solid ball is built by constructing R layers, radially outward; each layer is spaced r apart, and
     // points in each layer are spaced approx. r apart in a "grid" (a mapping from the cylinder to the sphere that has
     // been scaled s.t. projected points end up being approx. r apart.)
-    for (int rIter = 0; rIter <= R; rIter++) {
-      double r = rIter * rMinDist;
+    double r = rMinDist;
+    while (r < radius) {
       for (double z = -0.99; z <= 0.99; z += rMinDist) {
         double coeff = std::sqrt(1. - z * z);
         for (double theta = 0.0; theta <= 2.0 * PI; theta += rMinDist / coeff) {
@@ -167,6 +157,7 @@ void PoissonDiskSampler::addPointToSpatialLookupWithRadius(const SurfacePoint& n
           addPointToSpatialLookup(pos);
         }
       }
+      r += rMinDist;
     }
   }
 }
@@ -256,22 +247,19 @@ void PoissonDiskSampler::clearData() {
 /*
  * The final function.
  */
-std::vector<SurfacePoint> PoissonDiskSampler::sample(double rCoef_, int kCandidates_,
-                                                     std::vector<SurfacePoint> pointsToAvoid_, int rAvoidance,
-                                                     bool use3DAvoidanceRadius) {
+std::vector<SurfacePoint> PoissonDiskSampler::sample(const PoissonDiskOptions& options) {
 
   clearData();
 
   // Set parameters.
-  rCoef = rCoef_;
-  rMinDist = rCoef * meanEdgeLength;
+  rMinDist = options.minDist;
   sideLength = rMinDist / std::sqrt(3.0);
-  kCandidates = kCandidates_;
-  pointsToAvoid = pointsToAvoid_;
+  kCandidates = options.kCandidates;
+  pointsToAvoid = options.pointsToAvoid;
 
   // Add points to avoid.
   for (const SurfacePoint& pt : pointsToAvoid) {
-    addPointToSpatialLookupWithRadius(pt, rAvoidance - 1, use3DAvoidanceRadius);
+    addPointToSpatialLookupWithRadius(pt, options.minDistAvoidance, options.use3DAvoidance);
   }
 
   // Carry out sampling process for each connected component.

test/CMakeLists.txt

@@ -97,16 +97,17 @@ add_definitions(-DGC_TEST_ASSETS_ABS_PATH=\"${ABS_ASSETS_DIR}\")
 # Build the tests
 set(TEST_SRCS
   src/main_test.cpp
-  src/load_test_meshes.cpp
   src/eigen_interop_helpers_test.cpp
+  src/halfedge_geometry_test.cpp
   src/halfedge_mesh_test.cpp
   src/halfedge_mutation_test.cpp
-  src/halfedge_geometry_test.cpp
-  src/surface_misc_test.cpp
-  src/point_cloud_test.cpp
+  src/intrinsic_triangulation_test.cpp
   src/linear_algebra_test.cpp
+  src/load_test_meshes.cpp
+  src/point_cloud_test.cpp
+  src/poisson_disk_sampler_test.cpp
   src/stl_reader_test.cpp
-  src/intrinsic_triangulation_test.cpp
+  src/surface_misc_test.cpp
 )
 
 add_executable(geometry-central-test "${TEST_SRCS}")

test/src/poisson_disk_sampler_test.cpp

@@ -0,0 +1,52 @@
+
+#include <gtest/gtest.h>
+
+#include <geometrycentral/surface/poisson_disk_sampler.h>
+#include <geometrycentral/surface/simple_polygon_mesh.h>
+#include <geometrycentral/surface/surface_mesh_factories.h>
+
+
+using namespace geometrycentral;
+using namespace geometrycentral::surface;
+
+
+size_t sampleSquareDisk(double width, double sampling_distance) {
+  double const PTS[4][3] = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {1, 1, 0}};
+  unsigned const TRIS[2][3] = {{0, 1, 2}, {0, 2, 3}};
+
+  SimplePolygonMesh simpleMesh;
+  for (const auto& t : TRIS) simpleMesh.polygons.push_back({t[0], t[1], t[2]});
+
+  for (const auto& p : PTS) simpleMesh.vertexCoordinates.push_back({width * p[0], width * p[1], width * p[2]});
+
+  std::unique_ptr<ManifoldSurfaceMesh> mesh;
+  std::unique_ptr<VertexPositionGeometry> geometry;
+  std::tie(mesh, geometry) = makeManifoldSurfaceMeshAndGeometry(simpleMesh.polygons, simpleMesh.vertexCoordinates);
+
+  PoissonDiskOptions options;
+  options.minDist = sampling_distance;
+
+  // make tests reproducible
+  geometrycentral::util_mersenne_twister.seed(101);
+
+  PoissonDiskSampler sampler(*mesh, *geometry);
+  auto samples = sampler.sample(options);
+  return samples.size();
+}
+
+
+class PoissonDiskSamplerSuite : public ::testing::Test {};
+
+TEST_F(PoissonDiskSamplerSuite, PoissonDiskSamplerConstructor) {
+  // PoissonDiskSampler doesn't allow to set a random seed.
+  // To prevent flaky test failures we 'average' over 10 iterations.
+  size_t n1 = 0, n2 = 0;
+  for (size_t iter = 0; iter < 10; iter++) {
+    n1 += sampleSquareDisk(1.0, 0.1);
+    n2 += sampleSquareDisk(100.0, 100.0 * 0.1);
+  }
+
+  EXPECT_GT(n1, 600);
+  EXPECT_LT(n1, 720);
+  EXPECT_NEAR(n1, n2, 100);
+}