Feature/uniform point distribution

extensions/example/Jamfile.v2

@@ -13,3 +13,4 @@ project boost-geometry-examples-extensions
     ;
 
 build-project gis ;
+build-project random ;

extensions/example/random/Jamfile.v2

@@ -0,0 +1,16 @@
+# Boost.Geometry (aka GGL, Generic Geometry Library)
+#
+# Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
+# Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
+# Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
+
+# Use, modification and distribution is subject to the Boost Software License,
+# Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+# http://www.boost.org/LICENSE_1_0.txt)
+
+
+project boost-geometry-example-extensions-random
+    : # requirements
+    ;
+
+exe random_example : random_example.cpp ;

extensions/example/random/random_example.cpp

@@ -0,0 +1,82 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2019 Tinko Bartels, Berlin, Germany.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+// Random Example - showing random sampling in various geometries
+
+#include <fstream>
+#include <random>
+#include <cmath>
+
+#include <boost/geometry.hpp>
+#include <boost/geometry/extensions/random/uniform_point_distribution.hpp>
+
+const int samples = 100;
+
+template<typename Distribution, typename Generator, typename Mapper>
+void draw_samples(Distribution& d, Generator& g, Mapper& m) {
+    for (int i = 0 ; i < samples ; ++i)
+    {
+        m.map(d(g), "fill-opacity:1;fill:rgb(255,0,0);stroke:rgb(255,0,0);stroke-width:1", 1);
+    }
+}
+
+int main()
+{
+    namespace bg = boost::geometry;
+    typedef bg::model::point<double, 2, bg::cs::cartesian> point;
+    typedef bg::model::segment<point> segment;
+    typedef bg::model::box<point> box;
+    typedef bg::model::multi_point<point> multi_point;
+    typedef bg::model::polygon<point> polygon;
+
+    std::default_random_engine generator(1);
+    polygon poly;
+    boost::geometry::read_wkt("POLYGON((16 21,17.1226 17.5451,20.7553 17.5451,17.8164 15.4098,18.9389 11.9549,16 14.0902,13.0611 11.9549,14.1836 15.4098,11.2447 17.5451,14.8774 17.5451,16 21))", poly);
+    box b(point(0, 0), point(10, 10));
+    segment s(point(11, 0), point(21, 10));
+    multi_point mp;
+    for (double y = 11.0 ; y <= 21.0 ; y += 0.5)
+    {
+        for(double x = 0.0 ; x <= 10.0 ; x += 0.5)
+        {
+            bg::append(mp, point(x, y));
+        }
+    }
+    bg::random::uniform_point_distribution<box> box_dist(b);
+    bg::random::uniform_point_distribution<multi_point> mp_dist(mp);
+    bg::random::uniform_point_distribution<segment> seg_dist(s);
+    bg::random::uniform_point_distribution<polygon> poly_dist(poly);
+    std::ofstream svg("random.svg");
+    bg::svg_mapper<point> mapper(svg, 720, 720);
+    mapper.add(poly);
+    mapper.add(b);
+    mapper.add(mp);
+    mapper.add(s);
+    mapper.add(box(point(0, 22), point(22, 35)));
+    mapper.map(b, "fill-opacity:0.3;fill:rgb(51,51,153);stroke:rgb(51,51,153);stroke-width:2");
+    mapper.map(s, "fill-opacity:0.3;fill:rgb(51,51,153);stroke:rgb(51,51,153);stroke-width:2");
+    mapper.map(mp, "fill-opacity:0.3;fill:rgb(51,51,153);stroke:rgb(51,51,153);stroke-width:2", 2);
+    mapper.map(poly, "fill-opacity:0.3;fill:rgb(51,51,153);stroke:rgb(51,51,153);stroke-width:2");
+    draw_samples(box_dist, generator, mapper);
+    draw_samples(mp_dist, generator, mapper);
+    draw_samples(seg_dist, generator, mapper);
+    draw_samples(poly_dist, generator, mapper);
+
+    typedef bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree>> point_spherical;
+    typedef bg::model::box<point_spherical> box_spherical;
+    box_spherical sb(point_spherical(0, 0), point_spherical(90, 90));
+    bg::random::uniform_point_distribution<box_spherical> sb_dist(sb);
+    for (int i = 0 ; i < 10*samples ; ++i)
+    {
+        point_spherical p = sb_dist(generator);
+        double x = bg::get<0>(p) / 90 * 22;
+        double y = 32 - bg::get<1>(p) / 90 * 10;
+        mapper.map(point(x,y), "fill-opacity:1;fill:rgb(255,0,0);stroke:rgb(255,0,0);stroke-width:1", 1);
+    }
+    return 0;
+}

extensions/test/Jamfile.v2

@@ -27,4 +27,4 @@ build-project algorithms ;
 build-project gis ;
 build-project iterators ;
 build-project nsphere ;
-
+build-project random ;

extensions/test/random/Jamfile.v2

@@ -0,0 +1,11 @@
+# Boost.Geometry (aka GGL, Generic Geometry Library)
+#
+# Use, modification and distribution is subject to the Boost Software License,
+# Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+# http://www.boost.org/LICENSE_1_0.txt)
+
+test-suite boost-geometry-extensions-random
+    :
+    [ run random.cpp ]
+    ;
+

extensions/test/random/random.cpp

@@ -0,0 +1,144 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+// Unit Test
+
+// Copyright (c) 2019 Tinko Bartels, Berlin, Germany.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#include <random>
+#include <vector>
+#include <sstream>
+
+#include <geometry_test_common.hpp>
+
+#include <boost/geometry.hpp>
+#include <boost/geometry/extensions/random/uniform_point_distribution.hpp>
+
+typedef bg::model::point<double, 2, bg::cs::cartesian> point2d_cart;
+typedef bg::model::point<double, 3, bg::cs::cartesian> point3d_cart;
+typedef bg::model::point<double, 2, bg::cs::geographic<bg::degree>> point2d_geog;
+
+void test_geographic()
+{
+    //We check whether the generated points lie roughly along
+    //the great circle segment (<2 km distance on a sphere with
+    //the radius of earth) and are distributed uniformly with respect
+    //to great circle arc length.
+    typedef bg::model::linestring<point2d_geog> linestring;
+    linestring ls {{ 0.0, 0.0 }, { 45.0, 45.0 }, { 60.0, 60.0 }};
+    bg::random::uniform_point_distribution<linestring> l_dist(ls);
+    std::mt19937 generator(0);
+    int sample_count = 2000;
+    int count_below_45 = 0;
+    for (int i = 0 ; i < sample_count ; ++i)
+    {
+        point2d_geog sample = l_dist(generator);
+        BOOST_CHECK( bg::distance(sample, ls) < 2000 );
+        if(bg::get<0>(sample) < 45.0) count_below_45++;
+    }
+    double length_ratio = bg::distance(ls[0], ls[1]) / 
+        ( bg::distance(ls[0], ls[1]) + bg::distance(ls[1], ls[2]) );
+    double sample_ratio = ((double) count_below_45) / sample_count;
+    bool in_range = sample_ratio * 0.95 < length_ratio 
+        && sample_ratio * 1.05 > length_ratio;
+    BOOST_CHECK( in_range );
+
+    //We check whether the generated points lie in the spherical box
+    //(which is actually a triangle in this case) and whether the latitude
+    //is distributed as expected for uniform spherical distribution, using
+    //known area ratios of spherical caps.
+    typedef bg::model::box<point2d_geog> box;
+    box b {{ 0.0, 0.0 }, { 90.0, 90.0 }};
+    bg::random::uniform_point_distribution<box> b_dist(b);
+    int under_60 = 0;
+    for (int i = 0 ; i < sample_count ; ++i)
+    {
+        point2d_geog sample = b_dist(generator);
+        BOOST_CHECK( bg::within(sample, b) );
+        if(bg::get<1>(sample) < 60.0) ++under_60;
+    }
+    BOOST_CHECK_GT(under_60, 0.5 * 0.95 * sample_count);
+    BOOST_CHECK_LT(under_60, 0.5 * 1.05 * sample_count);
+}
+
+void test_polygon()
+{
+    //This test will test uniform sampling in polygon, which also checks
+    //uniform sampling in boxes. We check whether two equal distributions
+    //(copied using operator<< and operator>>) generate the same sequence
+    //of points and whether those points are uniformly distributed with
+    //respect to cartesian area.
+    typedef bg::model::polygon<point2d_cart> polygon;
+    polygon poly;
+    bg::read_wkt(
+        "POLYGON((16 21,17.1226 17.5451,20.7553 17.5451, 17.8164 15.4098,18.9389 11.9549,16 14.0902,13.0611 11.9549, 14.1836 15.4098,11.2447 17.5451,14.8774 17.5451,16 21))",
+        poly);
+    bg::random::uniform_point_distribution<polygon> poly_dist(poly);
+    bg::random::uniform_point_distribution<polygon> poly_dist2;
+    BOOST_CHECK( !(poly_dist == poly_dist2) );
+    std::stringstream ss;
+    ss << poly_dist;
+    ss >> poly_dist2;
+    BOOST_CHECK( poly_dist == poly_dist2 );
+    std::mt19937 generator(0), generator2(0);
+    for (int i = 0 ; i < 100 ; ++i)
+    {
+        point2d_cart sample1 = poly_dist(generator);
+        BOOST_CHECK( bg::equals(sample1, poly_dist2(generator2)) );
+        BOOST_CHECK( bg::within(sample1, poly) );
+    }
+    std::vector<point2d_cart> randoms;
+    const int uniformity_test_samples = 2000;
+    for (int i = 0 ; i < uniformity_test_samples ; ++i)
+    {
+        randoms.push_back(poly_dist(generator));
+    }
+    typedef bg::model::box<point2d_cart> box;
+    box env, lhalf;
+    bg::envelope(poly, env);
+    bg::set<bg::min_corner, 0>(lhalf, bg::get<bg::min_corner, 0>(env));
+    bg::set<bg::min_corner, 1>(lhalf, bg::get<bg::min_corner, 1>(env));
+    bg::set<bg::max_corner, 0>(lhalf, bg::get<bg::max_corner, 0>(env));
+    bg::set<bg::max_corner, 1>(lhalf,
+        (bg::get<bg::max_corner, 1>(env) + bg::get<bg::min_corner, 1>(env)) / 2);
+    std::vector<polygon> lower;
+    bg::intersection(lhalf, poly, lower);
+    double area_ratio = bg::area(lower[0])/bg::area(poly);
+    int in_lower = 0;
+    for (int i = 0 ; i < uniformity_test_samples ; ++i)
+    {
+        if(bg::within(randoms[i], lhalf))
+            ++in_lower;
+    }
+    double sample_ratio = ((double) in_lower ) / uniformity_test_samples;
+    BOOST_CHECK_GT( sample_ratio * 1.05, area_ratio );
+    BOOST_CHECK_LT( sample_ratio * 0.95, area_ratio );
+}
+
+void test_multipoint()
+{
+    typedef bg::model::multi_point<point3d_cart> multipoint;
+    multipoint mp {{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}};
+    int first = 0;
+    bg::random::uniform_point_distribution<multipoint> mp_dist(mp);
+    std::mt19937 generator(0);
+    int sample_count = 1000;
+    for (int i = 0 ; i < sample_count ; ++i)
+    {
+        point3d_cart sample = mp_dist(generator);
+        BOOST_CHECK( bg::within(sample, mp) );
+        if(bg::equals(sample, mp[0])) ++first;
+    }
+    BOOST_CHECK_GT(first * 1.05, sample_count / 3);
+    BOOST_CHECK_LT(first * 0.95, sample_count / 3);
+}
+
+int test_main(int, char* [])
+{
+    test_polygon();
+    test_geographic();
+    test_multipoint();
+    return 0;
+}

include/boost/geometry/extensions/random/strategies/agnostic/uniform_convex_fan.hpp

@@ -0,0 +1,95 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2019 Tinko Bartels, Berlin, Germany.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_EXTENSIONS_RANDOM_STRATEGIES_AGNOSTIC_UNIFORM_CONVEX_FAN_HPP
+#define BOOST_GEOMETRY_EXTENSIONS_RANDOM_STRATEGIES_AGNOSTIC_UNIFORM_CONVEX_FAN_HPP
+
+#include <random>
+#include <vector>
+#include <cstdlib>
+#include <cmath>
+
+#include <boost/geometry/algorithms/equals.hpp>
+#include <boost/geometry/algorithms/transform.hpp>
+#include <boost/geometry/algorithms/within.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace strategy { namespace uniform_point_distribution {
+
+// The following strategy is suitable for convex rings and polygons with
+// non-empty interior.
+template
+<
+    typename Point,
+    typename DomainGeometry,
+    typename TriangleStrategy,
+    typename SideStrategy   //Actually, we need a triangle area strategy here.
+>
+struct uniform_convex_fan
+{
+private:
+    std::vector<double> accumulated_areas;
+    // It is hard to see a reason not to use double here. If a triangles
+    // relative size is smaller than doubles epsilon, it is too unlikely to
+    // realistically occur in a random sample anyway.
+public:
+    uniform_convex_fan(DomainGeometry const& g)
+	{
+        accumulated_areas.push_back(0);
+        for (int i = 2 ; i < g.size() ; ++i) {
+            accumulated_areas.push_back(
+                accumulated_areas.back() +
+                std::abs(SideStrategy::template side_value<double, double>(
+                         *g.begin(),
+                         *(g.begin() + i - 1),
+                         *(g.begin() + i))));
+        }
+    }
+    bool equals(DomainGeometry const& l_domain,
+                DomainGeometry const& r_domain,
+                uniform_convex_fan const& r_strategy) const
+    {
+        if( l_domain.size() != r_domain.size() ) return false;
+        for (int i = 0; i < l_domain.size(); ++i) {
+            if( !boost::geometry::equals(*(l_domain.begin() + i),
+                                         *(r_domain.begin() + i)))
+                return false;
+        }
+        return true;
+    }
+    template<typename Gen>
+    Point apply(Gen& g, DomainGeometry const& d)
+    {
+        std::uniform_real_distribution<double> dist(0, 1);
+        double r = dist(g) * accumulated_areas.back(),
+               s = dist(g);
+        std::size_t i = std::distance(
+            accumulated_areas.begin(),
+            std::lower_bound(accumulated_areas.begin(),
+                             accumulated_areas.end(),
+                             r));
+        return TriangleStrategy::template map
+            <
+                double
+            >(* d.begin(),
+              *(d.begin() + i),
+              *(d.begin() + i + 1),
+              ( r - accumulated_areas[ i - 1 ]) /
+                ( accumulated_areas[ i ] - accumulated_areas[ i - 1 ] ),
+              s);
+    }
+    void reset(DomainGeometry const&) {};
+};
+
+}} // namespace strategy::uniform_point_distribution
+
+}} // namespace boost::geometry
+
+#endif //  BOOST_GEOMETRY_EXTENSIONS_RANDOM_STRATEGIES_AGNOSTIC_UNIFORM_CONVEX_FAN_HPP