Messing around with using halton sequences for quasi-random point filling of polygons (applied to mapnik @ ff8f3f1d4597eed89e04aa1dc7cb4755f6d4a744)

dot-density-points.diff

diff --git a/src/agg/process_point_symbolizer.cpp b/src/agg/process_point_symbolizer.cpp
index 6854903..abf1972 100644
--- a/src/agg/process_point_symbolizer.cpp
+++ b/src/agg/process_point_symbolizer.cpp
@@ -40,8 +40,58 @@
 // boost
 #include <boost/make_shared.hpp>
 
+#include <boost/utility.hpp>
+
 namespace mapnik {
 
+inline double halton_term(unsigned i, int base=2)
+{
+    double h=0;
+    double fac = 1.0/base;
+
+    while (i != 0)
+    {
+        double digit = i % base;
+        h = h + digit*fac;
+        i = (i-digit)/base;
+        fac = fac/base;
+    }
+    return h;
+};
+
+class quasi_random_sequence : private::boost::noncopyable
+{
+public:
+    quasi_random_sequence(unsigned terms, int base1=2,int base2=3)
+    : terms_(terms),
+      base1_(base1),
+      base2_(base2),
+      idx_(0) {}
+
+    void rewind(unsigned idx) const
+    {
+        idx_=idx;
+    }
+
+    bool vertex(double* x, double* y) const
+    {
+        while (idx_ < terms_)
+        {
+            *x = halton_term(idx_,base1_);
+            *y = halton_term(idx_,base2_);
+            ++idx_;
+            return true;
+        }
+        return false;
+    }
+
+private:
+    unsigned terms_;
+    int base1_;
+    int base2_;
+    mutable unsigned idx_;
+};
+
 template <typename T>
 void agg_renderer<T>::process(point_symbolizer const& sym,
                               mapnik::feature_impl & feature,
@@ -73,35 +123,71 @@ void agg_renderer<T>::process(point_symbolizer const& sym,
         for (unsigned i=0; i<feature.num_geometries(); ++i)
         {
             geometry_type const& geom = feature.get_geometry(i);
-            double x;
-            double y;
-            double z=0;
-            if (sym.get_point_placement() == CENTROID_POINT_PLACEMENT)
+            if (true) // test using halton sequences to fill polygon with points
             {
-                if (!label::centroid(geom, x, y))
-                    return;
+                double x;
+                double y;
+                double z=0;
+                //draw_geo_extent(geom.envelope(),mapnik::color("red"));
+                box2d<double> box = t_.forward(geom.envelope(),prj_trans);
+                double b_w = box.width();
+                double b_h = box.height();
+                unsigned term = 100; // density - needs to be user-configurable
+                // reduce the density of points for smaller  shapes
+                double res = ((t_.width()/b_w*.5) + (t_.height()/b_h*.5));
+                term /= res;
+                quasi_random_sequence seq(term);
+                while (seq.vertex(&x,&y))
+                {
+                    //std::clog << "x: " << x << " y: " << y << "\n";
+                    double dx = box.minx() + (x * b_w);
+                    double dy = box.miny() + (y * b_h);
+                    //std::clog << "dx: " << dx << " dy: " << dy << "\n";
+                    // make sure the point is within the poly
+                    double geo_x = dx;
+                    double geo_y = dy;
+                    t_.backward(&geo_x,&geo_y);
+                    prj_trans.forward(geo_x,geo_y,z);
+                    if (label::hit_test(geom,geo_x,geo_y,0))
+                    {
+                        render_marker(pixel_position(dx, dy),
+                                      **marker,
+                                      tr,
+                                      sym.get_opacity(),
+                                      sym.comp_op());
+                    }
+                }
             }
             else
             {
-                if (!label::interior_position(geom ,x, y))
-                    return;
-            }
-
-            prj_trans.backward(x,y,z);
-            t_.forward(&x,&y);
-            label_ext.re_center(x,y);
-            if (sym.get_allow_overlap() ||
-                detector_->has_placement(label_ext))
-            {
-
-                render_marker(pixel_position(x, y),
-                              **marker,
-                              tr,
-                              sym.get_opacity(),
-                              sym.comp_op());
-
-                if (!sym.get_ignore_placement())
-                    detector_->insert(label_ext);
+                double x;
+                double y;
+                double z=0;
+                if (sym.get_point_placement() == CENTROID_POINT_PLACEMENT)
+                {
+                    if (!label::centroid(geom, x, y))
+                        return;
+                }
+                else
+                {
+                    if (!label::interior_position(geom ,x, y))
+                        return;
+                }
+
+                prj_trans.backward(x,y,z);
+                t_.forward(&x,&y);
+                label_ext.re_center(x,y);
+                if (sym.get_allow_overlap() ||
+                    detector_->has_placement(label_ext))
+                {
+                    render_marker(pixel_position(x, y),
+                                  **marker,
+                                  tr,
+                                  sym.get_opacity(),
+                                  sym.comp_op());
+                    if (!sym.get_ignore_placement())
+                        detector_->insert(label_ext);
+                }
             }
         }
     }

example.md

<img src="http://f.cl.ly/items/0w27350J383s1f130704/halton.png />

test-case.xml

<!DOCTYPE Map>
<Map background-color="#b5d0d0" srs="+init=epsg:4326" minimum-version="0.7.2">
    <Style name="1">
        <Rule>
            <PolygonSymbolizer/>
        </Rule>
        <Rule>
            <PointSymbolizer />
        </Rule>
    </Style>

    <Layer name="point" srs="+init=epsg:4326">
        <StyleName>1</StyleName>
        <Datasource>
            <Parameter name="inline">
i|wkt
1|MULTIPOLYGON(((90 40,50 0,10 40,50 80,90 40)),((190 40,150 0,110 40,150 80,190 40)),((190 140,150 100,110 140,150 180,190 140)))
2|MULTIPOLYGON(((48 130,40 122,32 130,40 138,48 130)),((40 140,20 120,0 140,20 160,40 140)),((48 150,40 142,32 150,40 158,48 150)))
            </Parameter>
            <Parameter name="type">csv</Parameter>
        </Datasource>
    </Layer>

</Map>