Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Make line intersection computation more robust #937

Merged
merged 2 commits into from
Jul 10, 2023
Merged

Make line intersection computation more robust #937

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3a50712
Change to use DD line intersection
dr-jts Jul 10, 2023
6d23e3b
Make RobustLineIntersection unit test less brittle
dr-jts Jul 10, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 2 additions & 53 deletions src/algorithm/Intersection.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#include <vector>

#include <geos/algorithm/Intersection.h>
#include <geos/algorithm/CGAlgorithmsDD.h>

namespace geos {
namespace algorithm { // geos.algorithm
Expand All @@ -26,59 +27,7 @@ geom::CoordinateXY
Intersection::intersection(const geom::CoordinateXY& p1, const geom::CoordinateXY& p2,
const geom::CoordinateXY& q1, const geom::CoordinateXY& q2)
{
double minX0 = p1.x < p2.x ? p1.x : p2.x;
double minY0 = p1.y < p2.y ? p1.y : p2.y;
double maxX0 = p1.x > p2.x ? p1.x : p2.x;
double maxY0 = p1.y > p2.y ? p1.y : p2.y;

double minX1 = q1.x < q2.x ? q1.x : q2.x;
double minY1 = q1.y < q2.y ? q1.y : q2.y;
double maxX1 = q1.x > q2.x ? q1.x : q2.x;
double maxY1 = q1.y > q2.y ? q1.y : q2.y;

double intMinX = minX0 > minX1 ? minX0 : minX1;
double intMaxX = maxX0 < maxX1 ? maxX0 : maxX1;
double intMinY = minY0 > minY1 ? minY0 : minY1;
double intMaxY = maxY0 < maxY1 ? maxY0 : maxY1;

double midx = (intMinX + intMaxX) / 2.0;
double midy = (intMinY + intMaxY) / 2.0;

// condition ordinate values by subtracting midpoint
double p1x = p1.x - midx;
double p1y = p1.y - midy;
double p2x = p2.x - midx;
double p2y = p2.y - midy;
double q1x = q1.x - midx;
double q1y = q1.y - midy;
double q2x = q2.x - midx;
double q2y = q2.y - midy;

// unrolled computation using homogeneous coordinates eqn
double px = p1y - p2y;
double py = p2x - p1x;
double pw = p1x * p2y - p2x * p1y;

double qx = q1y - q2y;
double qy = q2x - q1x;
double qw = q1x * q2y - q2x * q1y;

double x = py * qw - qy * pw;
double y = qx * pw - px * qw;
double w = px * qy - qx * py;

double xInt = x/w;
double yInt = y/w;
geom::CoordinateXY rv;
// check for parallel lines
if (!std::isfinite(xInt) || !std::isfinite(yInt)) {
rv.setNull();
return rv;
}
// de-condition intersection point
rv.x = xInt + midx;
rv.y = yInt + midy;
return rv;
return CGAlgorithmsDD::intersection(p1, p2, q1, q2);
}


Expand Down
17 changes: 5 additions & 12 deletions tests/unit/algorithm/RobustLineIntersectionTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -305,43 +305,37 @@ void object::test<2>
// but apparently works now.
// Possibly normalization has fixed this?
//
#if 0 // fails: finds 1 intersection rather then two
template<>
template<>
void object::test<3>
()
{
std::vector<Coordinate> intPt;
intPt.push_back(Coordinate(2089426.5233462777, 1180182.3877339689));
intPt.push_back(Coordinate(2085646.6891757075, 1195618.7333999649));
intPt.push_back(Coordinate(2087600.4716727887, 1187639.7426241424));

checkIntersection(
"LINESTRING ( 2089426.5233462777 1180182.3877339689, 2085646.6891757075 1195618.7333999649 )",
"LINESTRING ( 1889281.8148903656 1997547.0560044837, 2259977.3672235999 483675.17050843034 )",
2,
1,
intPt, 0);
}
#endif // fails

#if 0 // fails: the intersection point doesn't match
// 4 - Outside envelope using HCoordinate method (testCmp5CaseWKT)
template<>
template<>
void object::test<4>
()
{
std::vector<Coordinate> intPt;
intPt.push_back(Coordinate(4348437.0557510145, 5552597.375203926));
intPt.push_back(Coordinate(4348440.8493873989, 5552599.2720221197));

checkIntersection(
"LINESTRING (4348433.262114629 5552595.478385733, 4348440.849387404 5552599.272022122 )",
"LINESTRING (4348433.26211463 5552595.47838573, 4348440.8493874 5552599.27202212 )",
1,
intPt, 0);
}
#endif // fails

#if 0 // fails: the intersection point doesn't match
// 5 - Result of this test should be the same as the WKT one!
// (testCmp5CaseRaw)
template<>
Expand All @@ -356,11 +350,10 @@ void object::test<5>
pt.push_back(Coordinate(4348440.8493874, 5552599.27202212));

std::vector<Coordinate> intPt;
intPt.push_back(Coordinate(4348437.0557510145, 5552597.375203926));
intPt.push_back(Coordinate(4348440.8493873989, 5552599.2720221197));

checkIntersection(pt, 1, intPt, 0);
}
#endif // fails

/**
* Test involving two non-almost-parallel lines.
Expand Down Expand Up @@ -453,7 +446,7 @@ void object::test<10>
"LINESTRING (-215.22279674875 -158.65425425385, -218.1208801283 -160.68343590235)",
1,
"POINT ( -215.22279674875 -158.65425425385 )",
0);
1.0e-10);
}

/**
Expand Down
23 changes: 23 additions & 0 deletions tests/unit/operation/relate/RelateOpTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,5 +120,28 @@ void object::test<3>()
}
}

// see https://github.com/locationtech/jts/issues/396
// testContainsNoding
template<>
template<>
void object::test<4> ()
{
checkRelate(
"LINESTRING (1 0, 0 2, 0 0, 2 2)",
"LINESTRING (0 0, 2 2)",
"101F00FF2" );
}

// see https://github.com/libgeos/geos/issues/933
// testContainsNoding
template<>
template<>
void object::test<5> ()
{
checkRelate(
"MULTILINESTRING ((0 0, 1 1), (0.5 0.5, 1 0.1, -1 0.1))",
"LINESTRING (0 0, 1 1)",
"1F1000FF2" );
}

} // namespace tut
49 changes: 47 additions & 2 deletions tests/xmltester/tests/general/TestRelateLL.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
<run>
<precisionModel scale="1.0" offsetx="0.0" offsety="0.0"/>

<case>
<desc>LL - disjoint, non-overlapping envelopes</desc>
Expand Down Expand Up @@ -308,7 +307,7 @@
</case>

<case>
<desc>LA - closed multiline / empty line</desc>
<desc>LL - closed multiline / empty line</desc>
<a>
MULTILINESTRING ((0 0, 0 1), (0 1, 1 1, 1 0, 0 0))
</a>
Expand All @@ -320,4 +319,50 @@
</test>
</case>

<case>
<desc>LL - test intersection node computation (see https://github.com/locationtech/jts/issues/396)</desc>
<a>
LINESTRING (1 0, 0 2, 0 0, 2 2)
</a>
<b>
LINESTRING (0 0, 2 2)
</b>
<test>
<op name="relate" arg1="A" arg2="B" arg3="101F00FF2">true</op>
</test>
<test><op name="contains" arg1="A" arg2="B"> true </op></test>
<test><op name="coveredBy" arg1="A" arg2="B"> false </op></test>
<test><op name="covers" arg1="A" arg2="B"> true </op></test>
<test><op name="crosses" arg1="A" arg2="B"> false </op></test>
<test><op name="disjoint" arg1="A" arg2="B"> false </op></test>
<test><op name="equalsTopo" arg1="A" arg2="B"> false </op></test>
<test><op name="intersects" arg1="A" arg2="B"> true </op></test>
<test><op name="overlaps" arg1="A" arg2="B"> false </op></test>
<test><op name="touches" arg1="A" arg2="B"> false </op></test>
<test><op name="within" arg1="A" arg2="B"> false </op></test>
</case>

<case>
<desc>LL - test intersection node computation (see https://github.com/libgeos/geos/issues/933)</desc>
<a>
MULTILINESTRING ((0 0, 1 1), (0.5 0.5, 1 0.1, -1 0.1))
</a>
<b>
LINESTRING (0 0, 1 1)
</b>
<test>
<op name="relate" arg1="A" arg2="B" arg3="1F1000FF2">true</op>
</test>
<test><op name="contains" arg1="A" arg2="B"> true </op></test>
<test><op name="coveredBy" arg1="A" arg2="B"> false </op></test>
<test><op name="covers" arg1="A" arg2="B"> true </op></test>
<test><op name="crosses" arg1="A" arg2="B"> false </op></test>
<test><op name="disjoint" arg1="A" arg2="B"> false </op></test>
<test><op name="equalsTopo" arg1="A" arg2="B"> false </op></test>
<test><op name="intersects" arg1="A" arg2="B"> true </op></test>
<test><op name="overlaps" arg1="A" arg2="B"> false </op></test>
<test><op name="touches" arg1="A" arg2="B"> false </op></test>
<test><op name="within" arg1="A" arg2="B"> false </op></test>
</case>

</run>
5 changes: 2 additions & 3 deletions tests/xmltester/tests/issue/issue-geos-398.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,9 +32,8 @@ POLYGON ((50.0044466166182886 0.0,
</b>
<test>
<op name="union" arg1="A" arg2="B">
MULTIPOLYGON (((0 25.950779066386122, 0 26.860827855779647, 0 29.83879230192533, 0 35.73921397193218, 60 13.183894681853793, 60 10.39852836763784, 60 7.802134559422377, 60 6.657099879646016, 60 6.510515132098641, 60 0, 55.31611974965083 0, 50.00444661661829 0, 43.31611974965083 0, 0 0, 0 14.034613342373582, 0 17.92266612923108, 0 21.587486526024364, 0 25.950779066386122),
(0 21.587486526024364, 34.025852439655786 6.898140262297274, 34.02585243965579 6.898140262297273, 0 21.587486526024364)),
((13.44557253233479 36, 60 36, 60 16.794451829809802, 60 16.36440115550932, 60 14.043996030454757, 2.9187843276549756 36, 11.8945390820011 36, 13.44557253233479 36)))
MULTIPOLYGON (((0 26.860827855779647, 0 29.83879230192533, 0 35.73921397193218, 60 13.183894681853793, 60 10.39852836763784, 60 7.802134559422377, 60 6.657099879646016, 60 6.510515132098641, 60 0, 55.31611974965083 0, 50.00444661661829 0, 43.31611974965083 0, 0 0, 0 14.034613342373582, 0 17.92266612923108, 0 21.587486526024364, 0 25.950779066386122, 0 26.860827855779647)),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

First coordinate in ring is different? Seems odd?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might just be a bit of non-determinism. Will investigate.

Copy link
Contributor Author

@dr-jts dr-jts Jul 10, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think there was a report that overlay would "rotate" results by a single coordinate over repeated operations. Perhaps this is the same phenomenon at work?

((60 36, 60 16.794451829809802, 60 16.36440115550932, 60 14.043996030454757, 2.9187843276549756 36, 11.8945390820011 36, 13.44557253233479 36, 60 36)))
</op>
</test>
</case>
Expand Down
2 changes: 1 addition & 1 deletion tests/xmltester/tests/robust/overlay/TestOverlay-pg-list.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ MULTIPOLYGON (((742591.4178018438 5087769.04526206, 742591.417801844 5087769.045

<test>
<op name="intersection" arg1="A" arg2="B">
MULTIPOLYGON (((434608.9999999999 2012571.9999999907, 434608.99999999977 2012571.9999999907, 434505.9999999997 2013031.9999999907, 434889.4374999999 2013081.8073453507, 434506 2013031.99999999, 434608.9999999999 2012571.9999999907)), ((434990.9999999999 2013094.99999999, 435264.9999999997 2013496.9999999898, 435311.2499999999 2013535.9783236892, 435265 2013496.99999999, 434991 2013094.99999999, 434990.9999999999 2013094.99999999)), ((436813.9999999997 2014014.99999999, 436802.9999999997 2014278.9999999893, 436803 2014278.9999999893, 436814 2014014.99999999, 436813.9999999997 2014014.99999999)), ((437427.6587183307 2014442.90461996, 437488.99999999977 2014458.99999999, 437372.9999999997 2016223.9999999898, 437686.7499999999 2016676.926282041, 437373 2016223.99999999, 437489 2014458.99999999, 437427.6587183307 2014442.90461996)))
MULTIPOLYGON (((437372.9999999997 2016223.9999999898, 437686.7499999997 2016676.9262820408, 437373 2016223.99999999, 437489 2014458.99999999, 437427.6587183306 2014442.90461996, 437488.99999999977 2014458.99999999, 437372.9999999997 2016223.9999999898)), ((434990.9999999999 2013094.99999999, 435264.9999999997 2013496.9999999898, 435311.2499999997 2013535.9783236892, 435265 2013496.99999999, 434991 2013094.99999999, 434990.9999999999 2013094.99999999)), ((434505.9999999997 2013031.9999999907, 434889.4374999997 2013081.8073453507, 434506 2013031.99999999, 434608.9999999999 2012571.9999999907, 434608.99999999977 2012571.9999999907, 434505.9999999997 2013031.9999999907)), ((436802.9999999997 2014278.9999999893, 436803 2014278.9999999893, 436814 2014014.99999999, 436813.9999999997 2014014.99999999, 436802.9999999997 2014278.9999999893)))
</op>
</test>

Expand Down
Loading