Modify

Opened 2 years ago

Closed 21 months ago

#11789 closed Bugs (fixed)

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

Reported by: Dave Stacey <drstacey@…> Owned by: Barend Gehrels
Milestone: Boost 1.61.0 Component: geometry
Version: Boost 1.59.0 Severity: Problem
Keywords: Cc:

Description

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.

Attachments (0)

Change History (3)

comment:1 Changed 2 years ago by 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.

comment:2 Changed 22 months ago by 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.

Last edited 22 months ago by awulkiew (previous) (diff)

comment:3 Changed 21 months ago by awulkiew

Milestone: To Be DeterminedBoost 1.61.0
Resolution: fixed
Status: newclosed

Modify Ticket

Change Properties
Set your email in Preferences
Action
as closed The owner will remain Barend Gehrels.
The resolution will be deleted.

Add Comment


E-mail address and name can be saved in the Preferences.

 
Note: See TracTickets for help on using tickets.