The new `--de9im` reports duplicate or even contradictory relations for identical geometries with ≥ 500 points

[hannahbast]

For two byte-identical geometries with at least DUPLICATE_REMOVAL_MIN_SIZE = 500 points, the --de9im output contains every relation twice. With fewer points, or with the classic relations (equals etc.), the same input is reported once. For polygons the two reports are identical. For linestrings they can disagree: the same directed pair is reported once with matrix M and once with M transposed (e.g. 10FF0FFF2 and 1FF00FFF2), and for byte-identical geometries the true matrix must be transpose-symmetric, so at least one of the two is wrong.

Here is a minimal example, two identical 1000-vertex ellipses (generator below):

$ ./spatialjoin --de9im < ellipse.txt
a	2FFF1FFF2	b
a	2FFF1FFF2	b
b	2FFF1FFF2	a
b	2FFF1FFF2	a

With 400 instead of 1000 vertices, each line appears once. The probable cause and a suggested (verified) fix are in the first comment below.

NOTE: I found this via the DE-9IM spatial join in QLever (ad-freiburg/qlever#3183). Joining the 162 Luxembourg OSM admin boundaries with themselves returns 180 instead of 162 pairs. The 18 extra rows are exactly the boundaries with ≥ 500 points.

python3 -c "
import math
n = 1000  # with 400, the output is correct
pts = ','.join(f'{5+3*math.cos(2*math.pi*i/n):.7f} {50+2*math.sin(2*math.pi*i/n):.7f}' for i in range(n+1))
print(f'a\tPOLYGON(({pts}))\nb\tPOLYGON(({pts}))')" > ellipse.txt

[hannahbast]

Probable cause: geometries with ≥ 500 points that are byte-identical to an earlier one are removed from the sweep, and their relations are later reconstructed via _refs. For a self-check pair (a == b), the two reference-expansion loops at the end of Sweeper::writeDE9IM (Sweeper.cpp:1994 at d1170ee) both emit the same relations. The classic write* functions guard the second loop with if (a != b || aSub != bSub) ("no need to check exactly the same direction again"). writeDE9IM lacks exactly this guard. Adding it fixes the problem:

   // handle references
   ...
+  // no need to check exactly the same direction again
+  if (a != b || aSub != bSub) {
     if (referersA != _refs.end()) {
       const auto& subs = referersA->second.find(aSub);
       ...
     }
+  }

Verified on a local build at d1170ee: the minimal example reports each relation once, the OSM Luxembourg join returns exactly the 162 expected pairs (matching the classic equals output), and outputs for small geometries are unchanged. For the linestring case, the patch keeps one of the two contradictory matrices. I did not verify whether the surviving one is the correct one for byte-identical closed ways, that may warrant a separate look.

[patrickbr]

For linestrings they can disagree: the same directed pair is reported once with matrix M and once with M transposed (e.g. 10FF0FFF2 and 1FF00FFF2)

This strongly hints at another small typo/seitendreher somewhere. The correct side may be reported purely by sheer luck. Will check this.

[patrickbr]

Should be fixed with #24

[patrickbr]

For linestrings they can disagree: the same directed pair is reported once with matrix M and once with M transposed (e.g. 10FF0FFF2 and 1FF00FFF2)

This strongly hints at another small typo/seitendreher somewhere. The correct side may be reported purely by sheer luck. Will check this.

Both were wrong, the hardcoded DE9IM matrix for equal lines simply had a small typo (boundary / interior intersection does not exist between equal lines but was hardcoded as 0, hence no transpose symmetry)