Boost 1.59.0: geometry::intersection() asserts for spherical_equatorial coordinate system

[Dave Stacey <drstacey@…>]

For disjoint polygons, geometry::intersection() assert()s when used with the spherical_equatorial coordinate system. assert() comes from geometry::math::detail::normalize_spheroidal_coordinates<>::apply(), as shown by the following code:

#include <iostream>
#include <deque>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/adapted/c_array.hpp>
#include <boost/foreach.hpp>

typedef boost::geometry::cs::spherical_equatorial<boost::geometry::degree> coordinate_system;
BOOST_GEOMETRY_REGISTER_C_ARRAY_CS(coordinate_system)

int main()
{
    typedef boost::geometry::model::point<double, 2, coordinate_system> point_t;
    typedef boost::geometry::model::polygon<point_t> polygon_t;

    double green_pts[][2] = {
        {-4.5726431789237223,52.142932977753595},
        {-4.5743166242433153,52.143359442355219},
        {-4.5739141406075410,52.143957260988416},
        {-4.5722406991324354,52.143530796430468},
        {-4.5726431789237223,52.142932977753595}};
    polygon_t green;
    boost::geometry::append(green, green_pts);

    double blue_pts[][2] = {
        {-4.5714644516017975,52.143819810922480},
        {-4.5670821923630358,52.143819810922480},
        {-4.5670821923630358,52.143649055226163},
        {-4.5714644516017975,52.143649055226163},
        {-4.5714644516017975,52.143819810922480}};
    polygon_t blue;
    boost::geometry::append(blue, green_pts);

    std::deque<polygon_t> output;
    // assert() comes from the following call...
    boost::geometry::intersection(green, blue, output);

    int i = 0;
    std::cout << "green && blue:" << std::endl;
    BOOST_FOREACH(polygon_t const& p, output)
    {
        std::cout << i++ << ": " << boost::geometry::area(p) << std::endl;
    }

    return 0;
}

assert() is generated in the call to geometry::intersection(). Tested with both g++-4.9.3 and MSVC 2005 SP1. This is a regression on Boost 1.57.0, where the assert() is not produced.

assert() is only generated when using the cs::spherical_equatorial coordinate system; if cs::cartesian is used instead then the assert() is not generated.

Many thanks for your work on Boost and for looking at this ticket.

[Dave Stacey <drstacey@…>]

Looking at the normalize_spheroidal_coordinates<>::apply() code, I tried defining BOOST_GEOMETRY_NORMALIZE_LATITUDE, but this did not fix the assert().

I also tried calling geometry::correct() on the two polygons before calling geometry::intersection(), but likewise this had no effect and the assert() still occured.

[awulkiew]

This assertion fails because expand() function is called for Points containing very big integral coordinates. The reason is that currently, by default, floating-point coordinates are rescaled into integral coordinates internally in set operations and these integral coordinates are passed into expand() to calculate BoundingBoxes. This probably also means that since these coordinates are normalized, the result is not correct.

To fix it we could disable the rescaling entirely or for non-cartesian coordinate systems, or always use cartesian robust_point_type in rescale policy.

To work around it, for now you could disable the rescaling by defining BOOST_GEOMETRY_NO_ROBUSTNESS.

[awulkiew]

Fix: ​https://github.com/boostorg/geometry/pull/341