Hi,

as I suspected, your example is highly degenerate; sides are almost overlapping and are ALMOST parallel but not quite; see the screenshot.

Red is the 1st rectangle, yellow is the 2nd rectangle and as well the result of SQL server query.

Now to the degeneracy, taking the top sides of both rectangles as vectors, Mathematica gives:

Dot[v1, v2]/(Norm[v1]*Norm[v2]) -> -1 (i.e., perfectly parallel in opposite directions) in machine precision. On the other hand, SetPrecision[Dot[v1, v2]/(Norm[v1]*Norm[v2]), 64] -> -0.999999999999999333866185224906075745820999145507812500000000000

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?