Repository navigation
OverlayNG area check heuristic causing issues when snapping noder is used #951
Description
After upgrading JTS from 1.18.2 to 1.19.0 ran into a number of tests that use difference operation failing with TopologyException: Result area inconsistent with overlay operation.
I tracked it down to OverlayUtil#isResultAreaConsistent (which was added in #812):
case OverlayNG.DIFFERENCE:
isConsistent = isLess(areaResult, areaA, AREA_HEURISTIC_TOLERANCE) // true; args are: 0, 6.512
&& isGreater(areaResult, areaA - areaB, AREA_HEURISTIC_TOLERANCE); // false; args are 0; 8.8818E-16
break;
In our case we are performing a difference operation on two almost-exactly-the-same geometries and expect an empty result.
Note: we are using snapping noder like in OverlayNGSnappingFunctions code.
Test case 1
This case fails with TopologyException:
@Test
void test2() throws ParseException {
Geometry g1 = new WKTReader().read("MULTIPOLYGON (((0.37676311 2.57570853, 7.28652472 0.00028375, 7.60034931 0.81686059, 0.50229292 3.4551325, 0.37676311 2.57570853)))");
Geometry g2 = new WKTReader().read("MULTIPOLYGON (((0.50229292 3.4551325, 7.60034931 0.81686059, 7.28652472 0.00028375, 0.37676311 2.57570853, 0.50229292 3.4551325)))");
// Note: in our code we have an utility function that wraps this logic
Geometry overlayResult = OverlayNG.overlay(
g1,
g2,
OverlayNG.DIFFERENCE,
g1.getPrecisionModel(),
new ValidatingNoder(new SnappingNoder(1e-4))
);
assertTrue(overlayResult.isEmpty());
}
In this case OverlayUtil#isResultAreaConsistent is reached once, and the expression in isGreater check works out to be:
v1 >= v2 * (1 - tol);0 >= 8.8819E-16 * (1 - 0.1)==>false
Test case 2
This case passes with the same data:
@Test // Passes
void test() throws ParseException {
Geometry p1 = new WKTReader().read("MULTIPOLYGON (((0.37676311 2.57570853, 7.28652472 0.00028375, 7.60034931 0.81686059, 0.50229292 3.4551325, 0.37676311 2.57570853)))");
Geometry p2 = new WKTReader().read("MULTIPOLYGON (((0.50229292 3.4551325, 7.60034931 0.81686059, 7.28652472 0.00028375, 0.37676311 2.57570853, 0.50229292 3.4551325)))");
Geometry result = p1.difference(p2);
assertEquals(result.getArea(), 0);
}
In this case it seems OverlayNGRobust is used instead, which results in multiple calls to OverlayUtil#isResultAreaConsistent (with slightly different values) and eventually it succeeds in passing the heuristic check.
Question
To me this feels like a bug in the heuristic check, but maybe our utility function for the difference operation is implemented in a suboptimal way?
If it's indeed a bug, and the fix is relatively simple (for someone who does not know JTS internals well), we'd be happy to contribute a PR with the fix.
Activity
dr-jts commentedon Oct 3, 2023
Thanks for the detailed bug report. Sorry for the delay in looking into this. It looks like a design flaw in the heuristic area check.
In this case the result is an empty geometry with zero area. The input geometries are essentially identical, but have rotated vertex sequences. Due to numeric precision issues (the eternal bugbear) the areas evaluate to be slightly different. This causes the area heuristic check to fail, incorrectly.
The heuristic needs to be modified to be less stringent when evaluating very small differences in area. I'm not sure what a suitable fix is, but will think about it. Suggestions welcome.
grimsa commentedon Oct 4, 2023
In case of
DIFFERENCEoperation this is the current code:AREA_HEURISTIC_TOLERANCEis0.1and the comparisons are defined asIf I understand correctly, the intention is to check that:
areaOfResult) is smaller or equal to the 110% of the area of the first input geometry (as when we remove something from the input geometry, it cannot become bigger)areaOfResult) is>=to 90% of the difference of areas between two input geometries.isGreater(positiveAlmostZero, 0, tolerance), which meanspositiveAlmostZero>= 0 * 90%, i.e. the check also works.7.
isGreater(0, negativeAlmostZero, tolerance)would result inreturn 0 >= negativeAlmostZero * 90%, which istrue, so the check works.8.
isGreater(0, positiveAlmostZero, tolerance)would result inreturn 0 >= positiveAlmostZero * 90%, which isfalse(like in the reported case), so the check does not work.I tried to think of possible solutions, but did not come up with any definitive answer quickly (but I can think more of it if needed).
One line of thinking was to add special handling for 0 or almost-zero values. But I did not come up with a way how to come up with "what is almost zero" without hardcoding a constant, e.g.
0.000001. Maybe that could be expressed as some % of the area of A or B (whichever is smaller and non-empty)? But that would likely require passing an extra arg toisGreater.Another idea was to base it on relative difference between sizes (one such example function is provided at the bottom of this page: https://c-faq.com/fp/fpequal.html), but then I think it does not work for non-near-zero cases as if there is no overlap between two similarly-sized polygons, the relative difference between area of the result (area A) and difference between areas (close to 0) can be very large.
dr-jts commentedon Oct 11, 2023
See #1005 for a fix for this issue.
@grimsa it would be great if you can confirm the fix works.
grimsa commentedon Oct 13, 2023
@dr-jts #1005 seems to fix the issue we had - thank you for the fix! It will let us benefit from JTS 1.20.0 once it is released.
However, while testing I noticed an unrelated minor change in the output of negative buffer operation since 1.19.0 - reported it in #1007