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?
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@lists.boost.org> on behalf of Heriberto Delgado via Boost-users <boost-users@lists.boost.org> Sent: Monday, March 11, 2019 6:42:36 PM To: boost-users@lists.boost.org 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?
Yes, I’m familiar with the myriad of issues that mishandling floating-point numbers can bring to the table. Notice, though, that the 1E-15 numbers you see there are essentially almost-zero. If you replace those with zero in both the SQL statement and the program, you will notice that the results don’t change at all. I tried. That’s why I don’t think these particular numbers are the origin of the issue here. El 11 mar. 2019, a la(s) 13:12, Stian Zeljko Vrba <vrba@quine.no> escribió:
Hi Heriberto, Op 11-3-2019 om 18:42 schreef Heriberto Delgado via Boost-users:
We're working on it. I didn't try your polygons yet, they look quite challenging. You can try to define BOOST_GEOMETRY_NO_ROBUSTNESS, as this will turn off internal rescaling and might help you. In the future, this will be the default. Alternatively, or on top of this, you might try to use long double. Please let me know your results. Regards, Barend
Thank you, Barend! Regarding the changes you mentioned, I found the following: - Long double precision for the polygons didn't seem to change the results. - When I added "#define BOOST_GEOMETRY_NO_ROBUSTNESS" to the top of the program, things DID change: Now the result list is empty (== the function did not detect any intersections). Hope this helps. El lun., 11 mar. 2019 a las 15:12, Barend Gehrels via Boost-users (< boost-users@lists.boost.org>) escribió:
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 [cid:fba3592a-a7ad-42fa-b2db-3299deabc7c1] ________________________________ From: Boost-users <boost-users-bounces@lists.boost.org> on behalf of Heriberto Delgado via Boost-users <boost-users@lists.boost.org> Sent: Monday, March 11, 2019 6:42:36 PM To: boost-users@lists.boost.org 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?
Hi Stian, Op 13-3-2019 om 18:14 schreef Stian Zeljko Vrba via Boost-users:
Indeed the lines were nearly collinear. Anyway, it is fixed in 1.70, see ticket https://github.com/boostorg/geometry/issues/566 Regards, Barend
participants (3)
-
Barend Gehrels -
Heriberto Delgado -
Stian Zeljko Vrba