Add uniform_point_distribution.

boostorg · tinko92 · Aug 12, 2019 · Aug 12, 2019 · Aug 12, 2019 · Aug 13, 2019
commit 0704d7e660f9faa57f4c68eb1bd95cebb921b704
diff --git a/extensions/example/Jamfile.v2 b/extensions/example/Jamfile.v2
@@ -13,3  13,4 @@ project boost-geometry-examples-extensions
     ;
 
 build-project gis ;
 build-project random ;
diff --git a/extensions/example/random/Jamfile.v2 b/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 ;
diff --git a/extensions/example/random/random_example.cpp b/extensions/example/random/random_example.cpp
@@ -0,0  1,78 @@
 // 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));
         }
     }
     auto box_dist = bg::random::uniform_point_distribution(b);
     auto mp_dist = bg::random::uniform_point_distribution(mp);
     auto seg_dist = bg::random::uniform_point_distribution(s);
     auto poly_dist = bg::random::uniform_point_distribution(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));
     auto sb_dist = bg::random::uniform_point_distribution(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;
 }
diff --git a/extensions/test/Jamfile.v2 b/extensions/test/Jamfile.v2
@@ -27,4  27,4 @@ build-project algorithms ;
 build-project gis ;
 build-project iterators ;
 build-project nsphere ;
-
 build-project random ;
diff --git a/extensions/test/random/Jamfile.v2 b/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 ]
     ;
 
diff --git a/extensions/test/random/random.cpp b/extensions/test/random/random.cpp
@@ -0,0  1,141 @@
 // 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>
 
 namespace bg = boost::geometry;
 
 using point2d_cart = bg::model::point<double, 2, bg::cs::cartesian>;
 using point3d_cart = bg::model::point<double, 3, bg::cs::cartesian>;
 using point2d_geog =
     bg::model::point<double, 2, bg::cs::geographic<bg::degree>>;
 
 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.
     using linestring = bg::model::linestring<point2d_geog>;
     linestring ls {{ 0.0, 0.0 }, { 45.0, 45.0 }, { 60.0, 60.0 }};
     auto l_dist = bg::random::uniform_point_distribution(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.
     using box = bg::model::box<point2d_geog>;
     box b {{ 0.0, 0.0 }, { 90.0, 90.0 }};
     auto b_dist = bg::random::uniform_point_distribution(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.
     using polygon = bg::model::polygon<point2d_cart>;
     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);
     auto poly_dist = bg::random::uniform_point_distribution(poly);
     decltype(poly_dist) 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));
     }
     using box = bg::model::box<point2d_cart>;
     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()
 {
     using multipoint = bg::model::multi_point<point3d_cart>;
     multipoint mp {{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}};
     int first = 0;
     auto mp_dist = bg::random::uniform_point_distribution(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;
 }
Original file line number	Diff line number	Diff line change
Expand Up		@@ -13,3 13,4 @@ project boost-geometry-examples-extensions
		;

		build-project gis ;
		build-project random ;