locationtech / jts

The JTS Topology Suite is a Java library for creating and manipulating vector geometry.
Other
1.95k stars 440 forks source link

DistanceOp.nearestPoints doesn't always return points covered by the input geometries #964

Closed joebarton3301 closed 1 year ago

joebarton3301 commented 1 year ago

I would expect coordinates returned by this method to always be covered by their respective input geometries, but I've found that's not true.

Here's an example:

GeometryFactory factory = new GeometryFactory();
Geometry inputCoordinate = factory.createPoint(new Coordinate(147.5666425807729, 6.587929382181272));
Geometry inputPolygon = factory.createPolygon(Arrays.asList(
                new Coordinate(148.9112, 5.5221),
                new Coordinate(146.7683, 5.5222),
                new Coordinate(147.5174, 2.2012),
                new Coordinate(148.9112, 5.5221)
).toArray(new Coordinate[0]));

Coordinate nearestCoordinateInPolygon = DistanceOp.nearestPoints(inputPolygon, inputCoordinate)[0]; // I'm getting 147.5665928459913, 5.522162747078912
Geometry nearestCoordinateAsGeometry = factory.createPoint(nearestCoordinateInPolygon);

assert inputPolygon.covers(nearestCoordinateAsGeometry); // this assertion fails

I assume this is just a floating point math precision issues as I see that nearestPoints is attempting to project the input point onto the input geometry's line segments. I've been dealing with it by reducing the input geometry by a small amount to guarantee that any nearest point calculations will fall within the geometry.

I don't see anywhere that nearestPoints make such a guarantee, but the behavior seems off and the workaround is a bit awkward. Is there a better way to approach this, or is it worth considering as a change to nearestPoints?

dr-jts commented 1 year ago

You're right, this is a numerical precision issue. The computation of a point along a line segment perpendicular to another offset point is unavoidable subject to numerical roundoff. So there's basically no way to guarantee that the result of nearest points always lies within the input geometries.

Really the issue lies with the intersects predicate. It would be nice if the predicate accepted a tolerance value, below which geometries would be considered to intersect. This is planned for future release (but no timeline yet).

As a workaround, you can compute the distance between geometries, and if it is less than the desired tolerance value consider that the geometries interesect.