ceph/src/boost/libs/geometry/test/algorithms/area.cpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2 // Unit Test
   3
   4 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
   5 // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
   6 // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
   7
   8 // This file was modified by Oracle on 2015, 2016.
   9 // Modifications copyright (c) 2015-2016, Oracle and/or its affiliates.
  10
  11 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  12
  13 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
  14 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
  15
  16 // Use, modification and distribution is subject to the Boost Software License,
  17 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  18 // http://www.boost.org/LICENSE_1_0.txt)
  19
  20
  21 #include <algorithms/test_area.hpp>
  22
  23 #include <boost/geometry/geometries/point_xy.hpp>
  24 #include <boost/geometry/geometries/point.hpp>
  25 #include <boost/geometry/geometries/box.hpp>
  26 #include <boost/geometry/geometries/ring.hpp>
  27 #include <boost/geometry/geometries/polygon.hpp>
  28
  29 #include <test_geometries/all_custom_ring.hpp>
  30 #include <test_geometries/all_custom_polygon.hpp>
  31 //#define BOOST_GEOMETRY_TEST_DEBUG
  32
  33 #include <boost/variant/variant.hpp>
  34
  35 template <typename Polygon>
  36 void test_polygon()
  37 {
  38     // Rotated square, length=sqrt(2) -> area=2
  39     test_geometry<Polygon>("POLYGON((1 1,2 2,3 1,2 0,1 1))", 2.0);
  40     test_geometry<Polygon>("POLYGON((1 1,2 2,3 1,2 0,1 1))", 2.0);
  41     test_geometry<Polygon>("POLYGON((0 0,0 7,4 2,2 0,0 0))", 16.0);
  42     test_geometry<Polygon>("POLYGON((1 1,2 1,2 2,1 2,1 1))", -1.0);
  43     test_geometry<Polygon>("POLYGON((0 0,0 7,4 2,2 0,0 0), (1 1,2 1,2 2,1 2,1 1))", 15.0);
  44 }
  45
  46
  47 template <typename P>
  48 void test_all()
  49 {
  50     test_geometry<bg::model::box<P> >("POLYGON((0 0,2 2))", 4.0);
  51     test_geometry<bg::model::box<P> >("POLYGON((2 2,0 0))", 4.0);
  52
  53     test_polygon<bg::model::polygon<P> >();
  54     test_polygon<all_custom_polygon<P> >();
  55
  56     // clockwise rings (second is wrongly ordered)
  57     test_geometry<bg::model::ring<P> >("POLYGON((0 0,0 7,4 2,2 0,0 0))", 16.0);
  58     test_geometry<bg::model::ring<P> >("POLYGON((0 0,2 0,4 2,0 7,0 0))", -16.0);
  59
  60     test_geometry<all_custom_ring<P> >("POLYGON((0 0,0 7,4 2,2 0,0 0))", 16.0);
  61
  62     // ccw
  63     test_geometry<bg::model::polygon<P, false> >
  64         ("POLYGON((0 0,0 7,4 2,2 0,0 0), (1 1,2 1,2 2,1 2,1 1))", -15.0);
  65
  66     test_geometry<bg::model::polygon<P, false> >
  67         ("POLYGON((1 0,0 1,-1 0,0 -1,1 0))", 2);
  68
  69     typedef typename bg::coordinate_type<P>::type coord_type;
  70     if (BOOST_GEOMETRY_CONDITION((boost::is_same<coord_type, double>::value)))
  71     {
  72         test_geometry<bg::model::polygon<P, false, false> >
  73             ("POLYGON((100000001 100000000, 100000000 100000001, \
  74                        99999999 100000000, 100000000  99999999))", 2);
  75     }
  76     else if (BOOST_GEOMETRY_CONDITION((boost::is_same<coord_type, float>::value)))
  77     {
  78         test_geometry<bg::model::polygon<P, false, false> >
  79             ("POLYGON((100001 100000, 100000 100001, \
  80                        99999 100000, 100000  99999))", 2);
  81     }
  82 }
  83
  84 template <typename Point>
  85 void test_spherical(bool polar = false)
  86 {
  87     typedef typename bg::coordinate_type<Point>::type ct;
  88     bg::model::polygon<Point> geometry;
  89
  90     // unit-sphere has area of 4-PI. Polygon covering 1/8 of it:
  91     // calculations splitted for ttmath
  92     ct const four = 4.0;
  93     ct const eight = 8.0;
  94     ct expected = four * boost::geometry::math::pi<ct>() / eight;
  95     bg::read_wkt("POLYGON((0 0,0 90,90 0,0 0))", geometry);
  96
  97     ct area = bg::area(geometry);
  98     BOOST_CHECK_CLOSE(area, expected, 0.0001);
  99
 100     // With strategy, radius 2 -> 4 pi r^2
 101     bg::strategy::area::huiller
 102         <
 103             typename bg::point_type<Point>::type
 104         > strategy(2.0);
 105
 106     area = bg::area(geometry, strategy);
 107     ct const two = 2.0;
 108     BOOST_CHECK_CLOSE(area, two * two * expected, 0.0001);
 109
 110     // Wrangel Island (dateline crossing)
 111     // With (spherical) Earth strategy
 112     bg::strategy::area::huiller
 113         <
 114             typename bg::point_type<Point>::type
 115         > spherical_earth(6373);
 116     bg::read_wkt("POLYGON((-178.7858 70.7852, 177.4758 71.2333, 179.7436 71.5733, -178.7858 70.7852))", geometry);
 117     area = bg::area(geometry, spherical_earth);
 118     // SQL Server gives: 4537.9654419375
 119     // PostGIS gives: 4537.9311668307
 120     // Note: those are Geographic, this test is Spherical
 121     BOOST_CHECK_CLOSE(area, 4506.6389, 0.001);
 122
 123     // Wrangel, more in detail
 124     bg::read_wkt("POLYGON((-178.568604 71.564148,-178.017548 71.449692,-177.833313 71.3461,-177.502838 71.277466 ,-177.439453 71.226929,-177.620026 71.116638,-177.9389 71.037491,-178.8186 70.979965,-179.274445 70.907761,-180 70.9972,179.678314 70.895538,179.272766 70.888596,178.791016 70.7964,178.617737 71.035538,178.872192 71.217484,179.530273 71.4383 ,-180 71.535843 ,-179.628601 71.577194,-179.305298 71.551361,-179.03421 71.597748,-178.568604 71.564148))", geometry);
 125     area = bg::area(geometry, spherical_earth);
 126     // SQL Server gives: 7669.10402181435
 127     // PostGIS gives: 7669.55565459832
 128     BOOST_CHECK_CLOSE(area, 7616.523769, 0.001);
 129
 130     // Check more at the equator
 131     /*
 132     select 1,geography::STGeomFromText('POLYGON((-178.7858 10.7852 , 179.7436 11.5733 , 177.4758 11.2333 , -178.7858 10.7852))',4326) .STArea()/1000000.0
 133     union select 2,geography::STGeomFromText('POLYGON((-178.7858 20.7852 , 179.7436 21.5733 , 177.4758 21.2333 , -178.7858 20.7852))',4326) .STArea()/1000000.0
 134     union select 3,geography::STGeomFromText('POLYGON((-178.7858 30.7852 , 179.7436 31.5733 , 177.4758 31.2333 , -178.7858 30.7852))',4326) .STArea()/1000000.0
 135     union select 0,geography::STGeomFromText('POLYGON((-178.7858 0.7852 , 179.7436 1.5733 , 177.4758 1.2333 , -178.7858 0.7852))',4326) .STArea()/1000000.0
 136     union select 4,geography::STGeomFromText('POLYGON((-178.7858 40.7852 , 179.7436 41.5733 , 177.4758 41.2333 , -178.7858 40.7852))',4326) .STArea()/1000000.0
 137     */
 138
 139     bg::read_wkt("POLYGON((-178.7858 0.7852, 177.4758 1.2333, 179.7436 1.5733, -178.7858 0.7852))", geometry);
 140     area = bg::area(geometry, spherical_earth);
 141     BOOST_CHECK_CLOSE(area, 14136.09946, 0.001); // SQL Server gives: 14064.1902284513
 142
 143
 144     bg::read_wkt("POLYGON((-178.7858 10.7852, 177.4758 11.2333, 179.7436 11.5733, -178.7858 10.7852))", geometry);
 145     area = bg::area(geometry, spherical_earth);
 146     BOOST_CHECK_CLOSE(area, 13760.2456, 0.001); // SQL Server gives: 13697.0941155193
 147
 148     bg::read_wkt("POLYGON((-178.7858 20.7852, 177.4758 21.2333, 179.7436 21.5733, -178.7858 20.7852))", geometry);
 149     area = bg::area(geometry, spherical_earth);
 150     BOOST_CHECK_CLOSE(area, 12987.8682, 0.001); // SQL Server gives: 12944.3970990317 -> -39m^2
 151
 152     bg::read_wkt("POLYGON((-178.7858 30.7852, 177.4758 31.2333, 179.7436 31.5733, -178.7858 30.7852))", geometry);
 153     area = bg::area(geometry, spherical_earth);
 154     BOOST_CHECK_CLOSE(area, 11856.3935, 0.001); // SQL Server gives: 11838.5338423574 -> -18m^2
 155
 156     bg::read_wkt("POLYGON((-178.7858 40.7852, 177.4758 41.2333, 179.7436 41.5733, -178.7858 40.7852))", geometry);
 157     area = bg::area(geometry, spherical_earth);
 158     BOOST_CHECK_CLOSE(area, 10404.627685523914, 0.001); // SQL Server gives: 10412.0607137119, -> +8m^2
 159
 160     // Concave
 161     bg::read_wkt("POLYGON((0 40,1 42,0 44,2 43,4 44,3 42,4 40,2 41,0 40))", geometry);
 162     area = bg::area(geometry, spherical_earth);
 163     BOOST_CHECK_CLOSE(area, 73538.2958, 0.001); // SQL Server gives: 73604.2047689719
 164
 165     // With hole POLYGON((0 40,4 40,4 44,0 44,0 40),(1 41,2 43,3 42,1 41))
 166     bg::read_wkt("POLYGON((0 40,0 44,4 44,4 40,0 40),(1 41,3 42,2 43,1 41))", geometry);
 167     area = bg::area(geometry, spherical_earth);
 168     BOOST_CHECK_CLOSE(area, 133233.844876, 0.001); // SQL Server gives: 133353.335
 169
 170     // around 0 meridian
 171     {
 172         bg::read_wkt("POLYGON((-10 0,-10 10,0 10,0 0,-10 0))", geometry);
 173         ct area1 = bg::area(geometry);
 174         bg::read_wkt("POLYGON((0 0,0 10,10 10,10 0,0 0))", geometry);
 175         ct area2 = bg::area(geometry);
 176         bg::read_wkt("POLYGON((-5 0,-5 10,5 10,5 0,-5 0))", geometry);
 177         ct area3 = bg::area(geometry);
 178         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 179         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 180         BOOST_CHECK_CLOSE(area1, 0.0303822, 0.001);
 181     }
 182     {
 183         bg::read_wkt("POLYGON((-10 -5,-10 5,0 5,0 -5,-10 -5))", geometry);
 184         ct area1 = bg::area(geometry);
 185         bg::read_wkt("POLYGON((0 -5,0 5,10 5,10 -5,0 -5))", geometry);
 186         ct area2 = bg::area(geometry);
 187         bg::read_wkt("POLYGON((-5 -5,-5 5,5 5,5 -5,-5 -5))", geometry);
 188         ct area3 = bg::area(geometry);
 189         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 190         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 191         BOOST_CHECK_CLOSE(area1, 0.0305, 0.001);
 192     }
 193     // around 180 meridian
 194     {
 195         bg::read_wkt("POLYGON((-180 0,-180 10,-170 10,-170 0,-180 0))", geometry);
 196         ct area1 = bg::area(geometry);
 197         bg::read_wkt("POLYGON((175 0,175 10,-175 10,-175 0,175 0))", geometry);
 198         ct area2 = bg::area(geometry);
 199         bg::read_wkt("POLYGON((170 0,170 10,180 10,180 0,170 0))", geometry);
 200         ct area3 = bg::area(geometry);
 201         bg::read_wkt("POLYGON((170 0,170 10,-180 10,-180 0,170 0))", geometry);
 202         ct area4 = bg::area(geometry);
 203         bg::read_wkt("POLYGON((180 0,180 10,-170 10,-170 0,180 0))", geometry);
 204         ct area5 = bg::area(geometry);
 205         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 206         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 207         BOOST_CHECK_CLOSE(area3, area4, 0.001);
 208         BOOST_CHECK_CLOSE(area4, area5, 0.001);
 209         BOOST_CHECK_CLOSE(area1, 0.0303822, 0.001);
 210     }
 211     {
 212         bg::read_wkt("POLYGON((-180 -5,-180 5,-170 5,-170 -5,-180 -5))", geometry);
 213         ct area1 = bg::area(geometry);
 214         bg::read_wkt("POLYGON((175 -5,175 5,-175 5,-175 -5,175 -5))", geometry);
 215         ct area2 = bg::area(geometry);
 216         bg::read_wkt("POLYGON((170 -5,170 5,180 5,180 -5,170 -5))", geometry);
 217         ct area3 = bg::area(geometry);
 218         bg::read_wkt("POLYGON((170 -5,170 5,-180 5,-180 -5,170 -5))", geometry);
 219         ct area4 = bg::area(geometry);
 220         bg::read_wkt("POLYGON((180 -5,180 5,-170 5,-170 -5,180 -5))", geometry);
 221         ct area5 = bg::area(geometry);
 222         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 223         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 224         BOOST_CHECK_CLOSE(area3, area4, 0.001);
 225         BOOST_CHECK_CLOSE(area4, area5, 0.001);
 226         BOOST_CHECK_CLOSE(area1, 0.0305, 0.001);
 227     }
 228     // around poles
 229 #ifdef BOOST_GEOMETRY_ENABLE_FAILING_TESTS
 230     {
 231         bg::read_wkt("POLYGON((0 80,-90 80,-180 80,90 80,0 80))", geometry);
 232         ct area1 = bg::area(geometry);
 233         bg::read_wkt("POLYGON((0 80,-90 80,180 80,90 80,0 80))", geometry);
 234         ct area2 = bg::area(geometry);
 235         bg::read_wkt("POLYGON((0 -80,90 -80,-180 -80,-90 -80,0 -80))", geometry);
 236         ct area3 = bg::area(geometry);
 237         bg::read_wkt("POLYGON((0 -80,90 -80,180 -80,-90 -80,0 -80))", geometry);
 238         ct area4 = bg::area(geometry);
 239         BOOST_CHECK_CLOSE(area1, area2, 0.001);
 240         BOOST_CHECK_CLOSE(area2, area3, 0.001);
 241         BOOST_CHECK_CLOSE(area3, area4, 0.001);
 242     }
 243 #endif
 244
 245     {
 246         bg::model::ring<Point> aurha; // a'dam-utr-rott.-den haag-a'dam
 247         bg::read_wkt("POLYGON((4.892 52.373,5.119 52.093,4.479 51.930,4.23 52.08,4.892 52.373))", aurha);
 248         if (polar)
 249         {
 250             // Create colatitudes (measured from pole)
 251             BOOST_FOREACH(Point& p, aurha)
 252             {
 253                 bg::set<1>(p, ct(90) - bg::get<1>(p));
 254             }
 255             bg::correct(aurha);
 256         }
 257         bg::strategy::area::huiller
 258             <
 259                 typename bg::point_type<Point>::type
 260             > huiller(6372.795);
 261         area = bg::area(aurha, huiller);
 262         BOOST_CHECK_CLOSE(area, 1476.645675, 0.0001);
 263
 264         // SQL Server gives: 1481.55595960659
 265         // for select geography::STGeomFromText('POLYGON((4.892 52.373,4.23 52.08,4.479 51.930,5.119 52.093,4.892 52.373))',4326).STArea()/1000000.0
 266     }
 267 }
 268
 269 template <typename P>
 270 void test_ccw()
 271 {
 272     typedef bg::model::polygon<P, false> ccw_polygon;
 273     // counterclockwise rings (second is wrongly ordered)
 274     test_geometry<ccw_polygon>("POLYGON((1 1,2 2,3 1,2 0,1 1))", -2.0);
 275     test_geometry<ccw_polygon>("POLYGON((1 1,2 0,3 1,2 2,1 1))", +2.0);
 276     test_geometry<ccw_polygon>("POLYGON((0 0,0 7,4 2,2 0,0 0))", -16.0);
 277     test_geometry<ccw_polygon>("POLYGON((0 0,2 0,4 2,0 7,0 0))", +16.0);
 278 }
 279
 280 template <typename P>
 281 void test_open()
 282 {
 283     typedef bg::model::polygon<P, true, false> open_polygon;
 284     test_geometry<open_polygon>("POLYGON((1 1,2 2,3 1,2 0))", 2.0);
 285     // Note the triangular testcase used in CCW is not sensible for open/close
 286 }
 287
 288 template <typename P>
 289 void test_open_ccw()
 290 {
 291     typedef bg::model::polygon<P, false, false> open_polygon;
 292     test_geometry<open_polygon>("POLYGON((1 1,2 0,3 1,2 2))", 2.0);
 293     // Note the triangular testcase used in CCW is not sensible for open/close
 294 }
 295
 296 template <typename P>
 297 void test_empty_input()
 298 {
 299     bg::model::polygon<P> poly_empty;
 300     bg::model::ring<P> ring_empty;
 301
 302     test_empty_input(poly_empty);
 303     test_empty_input(ring_empty);
 304 }
 305
 306 void test_large_integers()
 307 {
 308     typedef bg::model::point<int, 2, bg::cs::cartesian> int_point_type;
 309     typedef bg::model::point<double, 2, bg::cs::cartesian> double_point_type;
 310
 311     bg::model::polygon<int_point_type> int_poly;
 312     bg::model::polygon<double_point_type> double_poly;
 313
 314     std::string const polygon_li = "POLYGON((1872000 528000,1872000 192000,1536119 192000,1536000 528000,1200000 528000,1200000 863880,1536000 863880,1872000 863880,1872000 528000))";
 315     bg::read_wkt(polygon_li, int_poly);
 316     bg::read_wkt(polygon_li, double_poly);
 317
 318     double int_area = bg::area(int_poly);
 319     double double_area = bg::area(double_poly);
 320
 321     BOOST_CHECK_CLOSE(int_area, double_area, 0.0001);
 322 }
 323
 324 void test_variant()
 325 {
 326     typedef bg::model::point<double, 2, bg::cs::cartesian> double_point_type;
 327     typedef bg::model::polygon<double_point_type> polygon_type;
 328     typedef bg::model::box<double_point_type> box_type;
 329
 330     polygon_type poly;
 331     std::string const polygon_li = "POLYGON((18 5,18 1,15 1,15 5,12 5,12 8,15 8,18 8,18 5))";
 332     bg::read_wkt(polygon_li, poly);
 333
 334     box_type box;
 335     std::string const box_li = "BOX(0 0,2 2)";
 336     bg::read_wkt(box_li, box);
 337
 338     boost::variant<polygon_type, box_type> v;
 339
 340     v = poly;
 341     BOOST_CHECK_CLOSE(bg::area(v), bg::area(poly), 0.0001);
 342     v = box;
 343     BOOST_CHECK_CLOSE(bg::area(v), bg::area(box), 0.0001);
 344 }
 345
 346 int test_main(int, char* [])
 347 {
 348     test_all<bg::model::point<int, 2, bg::cs::cartesian> >();
 349     test_all<bg::model::point<float, 2, bg::cs::cartesian> >();
 350     test_all<bg::model::point<double, 2, bg::cs::cartesian> >();
 351
 352     test_spherical<bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> > >();
 353     //test_spherical<bg::model::point<double, 2, bg::cs::spherical<bg::degree> > >(true);
 354
 355     test_ccw<bg::model::point<double, 2, bg::cs::cartesian> >();
 356     test_open<bg::model::point<double, 2, bg::cs::cartesian> >();
 357     test_open_ccw<bg::model::point<double, 2, bg::cs::cartesian> >();
 358
 359 #ifdef HAVE_TTMATH
 360     test_all<bg::model::d2::point_xy<ttmath_big> >();
 361     test_spherical<bg::model::point<ttmath_big, 2, bg::cs::spherical_equatorial<bg::degree> > >();
 362 #endif
 363
 364     test_large_integers();
 365
 366     test_variant();
 367
 368     // test_empty_input<bg::model::d2::point_xy<int> >();
 369
 370     return 0;
 371 }