poission: add more thoroughly documented poisson sampler

geom/sample.cpp

@@ -15,57 +15,3 @@
  */
 
 #include "sample.hpp"
-
-#include "../point.hpp"
-
-#include <random>
-
-
-///////////////////////////////////////////////////////////////////////////////
-std::vector<util::point2f>
-util::geom::poisson_sample (util::extent2i area, float distance, int samples)
-{
-    std::vector<util::point2f> selected;
-
-    std::vector<util::point2f> candidates (samples);
-    std::uniform_real_distribution<float> dist_w (0, area.w);
-    std::uniform_real_distribution<float> dist_h (0, area.h);
-    std::random_device rd;
-    std::default_random_engine gen;
-
-    selected.push_back ({
-        dist_w (gen),
-        dist_h (gen),
-    });
-
-    auto min_distance = [&selected] (auto q) {
-        float min_d = INFINITY;
-
-        for (auto s: selected)
-            if (auto s_d = util::distance2 (s, q); s_d < min_d)
-                min_d = s_d;
-
-        return min_d;
-    };
-
-    while (1) {
-        util::point2f p;
-        float d = -INFINITY;
-
-        for (int i = 0; i < samples; ++i) {
-            util::point2f candidate { dist_w (gen), dist_h (gen) };
-
-            if (auto p_d = min_distance (candidate); p_d > d) {
-                p = candidate;
-                d = p_d;
-            }
-        }
-
-        if (std::sqrt (d) < distance)
-            break;
-
-        selected.push_back (p);
-    }
-
-    return selected;
-}

geom/sample.hpp

@@ -77,6 +77,84 @@ namespace util::geom {
 
     std::vector<util::point2f>
     poisson_sample (util::extent2i, float distance, int samples);
+
+
+    namespace surface {
+        // a generator of samples that lie on the surface of a shape
+        template <typename ShapeT>
+        class sampler;
+
+        template <typename ShapeT>
+        sampler (ShapeT) -> sampler<std::decay_t<ShapeT>>;
+
+        /// approximate a poisson disc sampling through mitchell's best candidate.
+        ///
+        /// try to keep adding a new point to a set. each new point is the
+        /// best of a set of candidates. the 'best' is the point that is
+        /// furthest from all selected points.
+        ///
+        /// \return a vector of the computed points
+        template <typename SamplerT, typename GeneratorT>
+        auto
+        poisson (SamplerT const &target,
+                 GeneratorT &&gen,
+                 float minimum_distance,
+                 size_t candidates_count)
+
+        {
+            using point_type = decltype (target (gen));
+            std::vector<point_type> selected;
+            std::vector<point_type> candidates;
+
+            // prime the found elements list with an initial point we can
+            // perform distance calculations on
+            selected.push_back (target (gen));
+
+            // keep trying to add one more new point
+            while (1) {
+                // generate a group of candidate points
+                candidates.clear ();
+                std::generate_n (
+                    std::back_inserter (candidates),
+                    candidates_count,
+                    [&] (void) {
+                        return target (gen);
+                    }
+                );
+
+                // find the point whose minimum distance to the existing
+                // points is the greatest (while also being greater than the
+                // required minimum distance);
+                float best_distance = -INFINITY;
+                size_t best_index;
+
+                for (size_t i = 0; i < candidates.size (); ++i) {
+                    auto const p = candidates[i];
+                    float d = INFINITY;
+
+                    // find the minimum distance from this candidate to the
+                    // selected points
+                    for (auto q: selected)
+                        d = util::min (d, distance (p, q));
+
+                    // record if it's the furthest away
+                    if (d > best_distance) {
+                        best_distance = d;
+                        best_index = i;
+                    }
+                }
+
+                // if we didn't find a suitable point then we give up and
+                // return the points we found, otherwise record the best point
+                if (best_distance < minimum_distance)
+                    break;
+
+                selected.push_back (candidates[best_index]);
+            }
+
+            return selected;
+        }
+    }
 }
 
 #endif

tools/poisson.cpp

@@ -1,20 +1,45 @@
 #include "cmdopt.hpp"
 #include "geom/sample.hpp"
+#include "geom/aabb.hpp"
 
 #include <cstdlib>
 #include <iostream>
+#include <random>
+
+
+template <size_t S>
+class util::geom::surface::sampler<util::extent<S,float>> {
+public:
+    sampler (util::extent<S,float> _target):
+        target (_target)
+    { ; }
+
+    template <typename GeneratorT>
+    util::point<S,float>
+    operator() (GeneratorT &&gen) const noexcept
+    {
+        util::point<S,float> p;
+        for (size_t i = 0; i < S; ++i)
+            p[i] = std::uniform_real_distribution<float> (0, target[i]) (gen);
+        return p;
+    }
+
+private:
+    util::extent<S,float> target;
+};
+
 
 ///////////////////////////////////////////////////////////////////////////////
 int
 main (int argc, char **argv)
 {
-    util::extent2i area {256, 256};
+    util::extent2f area {256, 256};
     float distance = 5.f;
     int samples = 15;
 
     util::cmdopt::parser opts;
-    opts.add<util::cmdopt::option::value<int>>   ('w', "width",   "width of the space to fill", area.w);
+    opts.add<util::cmdopt::option::value<float>> ('w', "width",   "width of the space to fill", area.w);
-    opts.add<util::cmdopt::option::value<int>>   ('h', "height",  "height of the space to fill", area.h);
+    opts.add<util::cmdopt::option::value<float>> ('h', "height",  "height of the space to fill", area.h);
     opts.add<util::cmdopt::option::value<float>> ('d', "distance", "minimum distance between samples", distance);
     opts.add<util::cmdopt::option::value<int>>   ('s', "samples", "number of samples per iteration", samples);
 
@@ -22,9 +47,15 @@ main (int argc, char **argv)
 
 
     std::cout << "<svg height='" << area.h << "' width='" << area.h << "'>";
-    for (auto p: util::geom::poisson_sample (area, distance, samples))
+    for (auto p: util::geom::surface::poisson (util::geom::surface::sampler (area),
+                                               std::default_random_engine (),
+                                               distance,
+                                               samples))
+    {
         std::cout << "<circle cx='" << p.x << "' cy='" << p.y << "' r='1' />";
-    std::cout << "</svg>";
+    }
+
+    std::cout << "</svg>\n";
 
     return EXIT_SUCCESS;
 }