Fix intersection computation of slightly overlapping polygons

[petrbel]

Problem: The polygon intersection computation doesn't work well when two polygons almost share an edge overlap very slightly. One of the symptoms is that intersection operation isn't commutative (as described in #1405)

Content: This PR consists of two commits: one adding a trivial not-working example (a unit test), two a proposed fix using equals2D instead of trivial independent coord-by-coord comparison. What i propose is to actually compare the dst point instead of float comparison of subtraction. That, should the points be very close, can produce 0 even for non-equivalent edges due to float arithmetics.

---

Example description: Below I enclose the problem description. It's similar to #1405 but, I believe, more straightforward and debuggable.

• Polygon A: POLYGON ((1 1, 0 0.5, 0 0, 1 1)) - a triangle with area 0.25
• Polygon B: POLYGON ((1 1, 0 0.49999999999999994, 0 1, 1 1)) - another triangle with area (almost) 0.25 which shares a common vertex 1 1 and a tiny overlap close to one of the edges with A and almost shares a second vertex 0 0.49999999999999994 with A (0 0.5)
• Apart from point (1 1) these polygons don't share any point. Their intersection is a single point and the intersection area is mathematically 0. Their union area is (almost) 0.5
• The intersection of these polygons a tiny neighborhood of edge (1 1) -> (0 0.5) with area close to 0
• Prior this PR, the intersection area of these polygons was 0.25 as well as the union area. The newly proposed unit tests demonstrate this fact.
• With the second commit applied, the behavior is fixed

Example visualization: I used dbeaver with postgis only for the visualization purposes. I know this postgis isn't up-to-date and linked geos is only 3.13.1, however, the same behavior is replicated directly in the tests.
postgis_full_version()

: POSTGIS="3.6.4 94d984b" [EXTENSION] PGSQL="180" GEOS="3.13.1-CAPI-1.19.2" SFCGAL="SFCGAL 2.0.0, CGAL 6.0, BOOST 1.83.0" PROJ="9.6.0 NETWORK_ENABLED=OFF URL_ENDPOINT=https://cdn.proj.org USER_WRITABLE_DIRECTORY=/var/lib/postgresql/.local/share/proj DATABASE_PATH=/usr/share/proj/proj.db" (compiled against PROJ 9.6.0) LIBXML="2.9.14" LIBJSON="0.18" LIBPROTOBUF="1.5.1" WAGYU="0.5.0 (Internal)" TOPOLOGY

The source sql in case anybody would like to replicate:

with cte as (
select
	'POLYGON ((1 1, 0 0.5, 0 0, 1 1))'::geometry as a,
    'POLYGON ((1 1, 0 0.49999999999999994, 0 1, 1 1))'::geometry as b
)
select
    a,
    b,
	st_area(a) as area_a,
	st_area(b) as area_b,
	st_area(st_intersection(a, b)) as area_ab,
	st_area(st_intersection(b, a)) as area_ba,
	postgis_full_version() as postgis_version
from cte

---

Disclaimer: while I found out about the problem myself and basically crafted the examples, I admit using a LLM for some hints and coding assistence. In case it's not allowed feel free to close the PR, however, please consider the revealed/described bug. I believe the code sheds some light to it.

---

Closes #1405

[dr-jts]

Why do you say "Apart from point (1 1) these polygons don't share any point"?

In fact they properly intersect, as can be seen in this image with magnified topology:

[petrbel]

@dr-jts you are right. I got confused by the y-axis.

However, the behavior is still flawed, isn't it? The intersection area is most definitely not 0.25. What do you think?

edit: I updated the initial description so that it's more clear

[dr-jts]

@dr-jts you are right. I got confused by the y-axis.

However, the behavior is still flawed, isn't it? The intersection area is most definitely not 0.25. What do you think?

Looks like a bug in GEOS. JTS gives the expected intersection result (with area = 2.7755575615628914E-1 ):

POLYGON ((0 0.5, 1 1, 0 0.4999999999999999, 0 0.5))

[dbaston]

Looks like a bug in GEOS

Yes and no, I think. When I pull this up in the debugger, JTS and GEOS return the same result for all invocations of HalfEdge::compareAngularDirection. In particular, both return 0 for when comparing the edges (1 1, 0 0.5) and (1 1, 0 0.49999999999999994). But we know the angular direction is not equal here. I suspect that treating this as equal sets us up for an insertion order dependency later on, and it currently happens to work out correctly in JTS and not in GEOS.

I think the fix here is OK, though I would be more comfortable if we also asserted that orig().equals2D(e->orig())

[petrbel]

@dbaston I think your exaplanation is correct. My explanation is that if we use the coordinate-wise float comparison and the difference is less than what floats can handle, the points are considered the same. Then, the quadrands tie-breaker never runs.
This, I believe, leads to HalfEdge::isEdgesSorted (by using compareTo which finally uses compareAngularDirection) considering an edge unsorted (because the condition checks for > 0). In such case the order of the nodes is determined by the order of their insertion. And that, in my understanding, might swap the interior and exterior of the intersection, leading to the whole input polygon is considered the intersection.

If that makes sense.

---

I think the fix here is OK, though I would be more comfortable if we also asserted that orig().equals2D(e->orig())

Good idea. I added the asserts

Please let me know if you have any further remarks

[dr-jts]

Agreed, it's always better to use more robust logic that avoids arithmetic operations.

[dr-jts]

Mpre notes:

• compareAngularDirection should use Quadrant.quadrant(Coordinate, Coordinate), since that's more robust by avoiding the arithmetic to compute direction vector components (see Make HalfEdge.compareAngularDirection more robust locationtech/jts#1224 and Improve Quadrant computation robustness by using vector endpoints locationtech/jts#1226)
• the methods directionX and directionY are badly named. They are not the X and Y ordinates of the direction point, but the direction vector. They should probably be called directionVectorX and directionVectorY. And with the above change they can be removed.