Merge pull request #434 from awulkiew/fix/geographic_union

Fix for union in geographic CS not generating result if it's too big
This commit is contained in:
Adam Wulkiewicz
2017-11-15 00:25:03 +01:00
committed by GitHub
5 changed files with 118 additions and 34 deletions

View File

@@ -74,7 +74,8 @@ template
inline OutputIterator add_rings(SelectionMap const& map,
Geometry1 const& geometry1, Geometry2 const& geometry2,
RingCollection const& collection,
OutputIterator out)
OutputIterator out,
bool require_positive_area = true)
{
typedef typename SelectionMap::const_iterator iterator;
typedef typename SelectionMap::mapped_type property_type;
@@ -123,7 +124,8 @@ inline OutputIterator add_rings(SelectionMap const& map,
// This cannot be done earlier (during traversal), not
// everything is figured out yet (sum of positive/negative rings)
if (geometry::num_points(result) >= min_num_points
&& math::larger(geometry::area(result), zero))
&& (! require_positive_area
|| math::larger(geometry::area(result), zero)))
{
*out++ = result;
}

View File

@@ -387,7 +387,13 @@ std::cout << "traverse" << std::endl;
assign_parents(geometry1, geometry2, rings, selected_ring_properties, strategy);
return add_rings<GeometryOut>(selected_ring_properties, geometry1, geometry2, rings, out);
// NOTE: There is no need to check result area for union because
// as long as the polygons in the input are valid the resulting
// polygons should be valid as well.
// This also solves the issue with non-cartesian CSes, the result may
// be too big so the area is negative but it's returned anyway.
return add_rings<GeometryOut>(selected_ring_properties, geometry1, geometry2, rings, out,
OverlayType != overlay_union);
}
template <typename RobustPolicy, typename OutputIterator, typename Strategy>

View File

@@ -4,8 +4,8 @@
# Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
# Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
#
# This file was modified by Oracle on 2014, 2015.
# Modifications copyright (c) 2014-2015, Oracle and/or its affiliates.
# This file was modified by Oracle on 2014, 2015, 2017.
# Modifications copyright (c) 2014-2017, Oracle and/or its affiliates.
#
# Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
# Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -18,6 +18,7 @@ test-suite boost-geometry-algorithms-union
:
[ run union.cpp : : : <define>BOOST_GEOMETRY_TEST_ONLY_ONE_TYPE
: algorithms_union ]
[ run union_aa_geo.cpp : : : : algorithms_union_aa_geo ]
[ run union_linear_linear.cpp : : : : algorithms_union_linear_linear ]
[ run union_multi.cpp : : : <define>BOOST_GEOMETRY_TEST_ONLY_ONE_TYPE
: algorithms_union_multi ]

View File

@@ -482,33 +482,6 @@ void test_areal()
1, 1, -1, 220.5);
}
void test_geographic()
{
typedef bg::model::point<double, 2, bg::cs::geographic<bg::degree> > point;
typedef bg::model::polygon<point> polygon;
typedef bg::model::multi_polygon<polygon> multipolygon;
bg::srs::spheroid<double> sph(6378137.0000000000, 6356752.3142451793);
bg::strategy::intersection::geographic_segments<> is(sph);
bg::strategy::area::geographic<point> as(sph);
polygon p1, p2;
boost::geometry::read_wkt("POLYGON((16 15,-132 10,-56 89,67 5,16 15))", p1);
boost::geometry::read_wkt("POLYGON((101 49,12 40,-164 10,117 0,101 49))", p2);
multipolygon result;
boost::geometry::union_(p1, p2, result, is);
double result_area = bg::area(result, as);
BOOST_CHECK(boost::size(result) == 1
&& boost::size(bg::exterior_ring(bg::range::at(result, 0))) == 9
&& boost::size(bg::interior_rings(bg::range::at(result, 0))) == 0);
BOOST_CHECK_CLOSE(result_area, 144265751613509.06, 0.001);
}
template <typename P>
void test_all()
{
@@ -584,7 +557,5 @@ int test_main(int, char* [])
#endif
#endif
test_geographic();
return 0;
}

View File

@@ -0,0 +1,104 @@
// Boost.Geometry
// Unit Test
// Copyright (c) 2017, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// 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 <iostream>
#include <string>
#include "test_union.hpp"
#include <boost/geometry/geometries/point_xy.hpp>
struct exterior_points_counter
{
exterior_points_counter() : count(0) {}
template <typename Polygon>
void operator()(Polygon const& poly)
{
count += boost::size(bg::exterior_ring(poly));
}
std::size_t count;
};
struct interiors_counter
: exterior_points_counter
{
template <typename Polygon>
void operator()(Polygon const& poly)
{
count += boost::size(bg::interior_rings(poly));
}
};
void test_geographic_one(std::string const& wkt1, std::string const& wkt2,
std::size_t count, std::size_t exterior_points_count, std::size_t interiors_count,
double expected_area)
{
typedef bg::model::point<double, 2, bg::cs::geographic<bg::degree> > point;
typedef bg::model::polygon<point> polygon;
typedef bg::model::multi_polygon<polygon> multipolygon;
bg::srs::spheroid<double> sph(6378137.0000000000, 6356752.3142451793);
bg::strategy::intersection::geographic_segments<> is(sph);
bg::strategy::area::geographic<point> as(sph);
polygon p1, p2;
boost::geometry::read_wkt(wkt1, p1);
boost::geometry::read_wkt(wkt2, p2);
multipolygon result;
boost::geometry::union_(p1, p2, result, is);
double result_area = bg::area(result, as);
std::size_t result_count = boost::size(result);
std::size_t result_exterior_points = std::for_each(boost::begin(result),
boost::end(result),
exterior_points_counter()).count;
std::size_t result_interiors = std::for_each(boost::begin(result),
boost::end(result),
interiors_counter()).count;
BOOST_CHECK_EQUAL(result_count, count);
BOOST_CHECK_EQUAL(result_exterior_points, exterior_points_count);
BOOST_CHECK_EQUAL(result_interiors, interiors_count);
BOOST_CHECK_CLOSE(result_area, expected_area, 0.001);
}
void test_geographic()
{
// input ok and the result is ok
test_geographic_one("POLYGON((16 15,-132 10,-56 89,67 5,16 15))",
"POLYGON((101 49,12 40,-164 10,117 0,101 49))",
1, 9, 0, 144265751613509.06);
// input ok but the result is too big
test_geographic_one("POLYGON((16 -15,-132 -22,-56 89,67 -29,16 -15))",
"POLYGON((101 49,12 40,-164 -21,117 -61,101 49))",
1, 9, 0, -163427005620080.0);
// the second polygon is reversed i.e. it covers more than half of the globe
// so the result is also too big
test_geographic_one("POLYGON((16 -15,-132 -22,-56 89,67 -29,16 -15))",
"POLYGON((101 49,117 -61,-164 -21,12 40,101 49))",
1, 7, 0, -125258931656228.08);
}
int test_main(int, char* [])
{
test_geographic();
return 0;
}