Re: [Boost-users] Issue with boost::geometry::intersection

Subject: Re: [Boost-users] Issue with boost::geometry::intersection
From: Stian Zeljko Vrba (vrba_at_[hidden])
Date: 2019-03-11 19:12:13

Having much experience with computational geometry, at a first glance I'm not convinced that SQLServer returns the correct result either: it seems you have a very degenerate case there (coordinates around 1e-15 with other coordinates in the "normal" range) that exhausts the floating point precision/dynamic range.

You should verify the result using arbitrary-precision or rational arithmetic using Mathematica or another tool.

In any case: you cannot handle correctly all possible inputs using standard floating-point arithmetic. It is always possible to construct a case where the intersection point(s) fall "in between" representable numbers in FP arithmetic and then it's an implementation detail whether an intersection will be reported or not. Or whether the intersection point will go haywire due to simple (numerically unstable) implementation.

2D intersection requires solving a linear system and the required precision for the output may be much greater (double?) than the precision of the inputs. The precision requirements multiply if you're about to do further computations with the computed points, so arbitrary precision is no silver bullet as memory requirements (and speed impact) can grow exponentially.

I did a comprehensive literature survey for the application I was working on and indeed I haven't found a silver bullet. You have to work with your clients and gain requirements in order to succeed here.

-- Stian

________________________________
From: Boost-users <boost-users-bounces_at_[hidden]> on behalf of Heriberto Delgado via Boost-users <boost-users_at_[hidden]>
Sent: Monday, March 11, 2019 6:42:36 PM
To: boost-users_at_[hidden]
Cc: Heriberto Delgado
Subject: [Boost-users] Issue with boost::geometry::intersection

I'm having an issue while using boost::geometry to find the intersection of two specific polygons, and looking at the evidence in front of me I can't help but think this might be a bug in the library.

The two polygons in question are:

POLYGON((16.8739 0.353458,8.80872e-14 5.50545e-15,0.266003 -12.6988,17.1399 -12.3454,16.8739 0.353458))

POLYGON((-4.53257e-19 2.16383e-17,9.99781 0.209424,10.1459 -6.86009,0.148086 -7.06952,-4.53257e-19 2.16383e-17))

Both of them look like two slightly rotated rectangles, first one being bigger than the other. They intersect at their upper left corner.

In order to verify that these polygons CAN be intersected, I ran the following SQL sentence through SQL Server Management Studio (17.9.1):

SELECT geometry::STGeomFromText('POLYGON((16.8739 0.353458,8.80872e-14 5.50545e-15,0.266003 -12.6988,17.1399 -12.3454,16.8739 0.353458))', 0)
.STIntersection(
geometry::STGeomFromText('POLYGON((-4.53257e-19 2.16383e-17,9.99781 0.209424,10.1459 -6.86009,0.148086 -7.06952,-4.53257e-19 2.16383e-17))', 0)
).STAsText()

The resulting WKT is:

POLYGON ((0.1480860710144043 -7.0695199966430664, 10.145899772644043 -6.8600900173187256, 9.997809886932373 0.20942401885986328, 0 0, 0.1480860710144043 -7.0695199966430664))

which looks about right for the requested intersection rectangle. SSMS itself displays a nice rectangle when you remove ".STAsText()" from that sentence.

However, when I tried to duplicate that result using the following program using Visual Studio 2017 (15.9.8) to create an x64 console app, using boost 1.69 as a NuGet:

-----------------------------------------------------------------------------------------------------------------
#include <boost/geometry.hpp>
#include <vector>
#include <iostream>
int main(int argc, char** argv)
{
 boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>> first;
 boost::geometry::read_wkt("POLYGON((16.8739 0.353458,8.80872e-14 5.50545e-15,0.266003 -12.6988,17.1399 -12.3454,16.8739 0.353458))", first);
 boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>> second;
 boost::geometry::read_wkt("POLYGON((-4.53257e-19 2.16383e-17,9.99781 0.209424,10.1459 -6.86009,0.148086 -7.06952,-4.53257e-19 2.16383e-17))", second);
 std::vector<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double>>> results;
 boost::geometry::intersection(first, second, results);
 for (auto& result : results)
 {
  std::cout << boost::geometry::wkt(result) << std::endl;
  std::cout << boost::geometry::is_valid(result) << std::endl;
 }
 return 0;
}
-----------------------------------------------------------------------------------------------------------------

the resulting polygon is:

POLYGON((0 -8.88178e-16,0.148086 -7.06952,0.148086 -7.06952,0 -8.88178e-16))

which you can easily see it's not a valid polygon, even though boost::geometry returns "1" when is_valid() is called upon it.

I was planning to use boost::geometry on a large scale server application that processes quite large amounts of geometry in a similar way, and this is seriously impairing my ability to do so.

Do you guys know what is going on here, and if there is a solution or workaround for it?