Geometry Boost List: [geometry] Expectation for robustness and rescale policy

Subject: [geometry] Expectation for robustness and rescale policy
From: Martinez, Max (max.martinez_at_[hidden])
Date: 2016-08-02 14:34:34

A recent upgrade from boost 1.55 to 1.61 resulted in a number of unit test failures in our development environment related to the robustness polices in boost geometry (specifically the rescale policy in the context of intersection) that seem to now be applied by default. The failures go away (previous behavior is restored) if we set BOOST_GEOMETRY_NO_ROBUSTNESS, but we were wondering what our expectations should be with respect to robustness (there doesn't seem to be any documentation on it). I'm not a COGO expert so after googling around a bit my impression is that a robust algorithm is one that behaves correctly (and fast) even for "degenerate" cases.

The failures are related to floating point differences but one that was surprising to us was that intersections generated on horizontal (and vertical?) line segments did not result in points that had the exact y (x) value as the points representing the line segment intersected.

The example below illustrates this. It is a test that has been massaged a bit so that it can be examined outside our environment. It passes if the coordinate type is float, but not if it is double (which our adaptation is based upon). I understand that the results will be slightly different due to the scaling to integer and back to double, but is it expected that the horizontal line intersection, which can easily (and will) be inspected by the geospatial end user, might not have a y value that matches the horizontal line (value ends up as 3647747.9999977136)? Is this also possible if robustness is disabled (i.e., were/are we just getting lucky with this set of values)?

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/geometries/ring.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/convert.hpp>
.
.
.
       typedef boost::geometry::model::d2::point_xy<double> point_2d;
       typedef boost::geometry::model::box<point_2d> box_2d;
       typedef boost::geometry::model::ring<point_2d> ring_2d;

       box_2d b = boost::geometry::make<box_2d>(598413.0, 3647748.0, 598641.0, 3647586.5);
       ring_2d r1;
       boost::geometry::convert(b, r1);

       b = boost::geometry::make<box_2d>(598427.84499100002, 3647759.164202, 598607.84499100002, 3647619.164202);
       ring_2d r2;
       boost::geometry::convert(b, r2);

       {
             using namespace boost::geometry;
             std::deque<ring_2d> outputs;
             ring_2d &rRing1 = r1;
             ring_2d &rRing2 = r2;

             ring_2d rRing1Corrected = rRing1;
             boost::geometry::correct(rRing1Corrected);
             ring_2d rRing2Corrected = rRing2;
             boost::geometry::correct(rRing2Corrected);

             intersection(rRing1Corrected, rRing2Corrected, outputs);

             box_2d bounds;
             envelope(outputs[0], bounds);

       CPPUNIT_ASSERT_EQUAL(boost::geometry::traits::coordinate_type<point_2d>::type(3647748.0), bounds.max_corner().y());
       }