TopologyException on operations on valid geometries [JTS works!]

[yabo]

I get a TopologyException for operations between these two geometries (reduced as much as I could from original geometry) :

wkb1 = 01060000000100000001030000000100000007000000000000005C3613C1000000000028A54000000000208D14C1000000000028A5C000000000A83F15C100000000C0A6D040981124BD63B414C1B82DC4EB7C969940CA5EE746E76814C120462E0E5B4A8A40323BFF843F4B14C1C77B5B2D590C9AC0000000005C3613C1000000000028A540

wkb2 = 01060000000100000001030000000100000005000000EF2CF4FC0ED315C1F4792E801F82B6C0477F50A83BA313C12242DCE0BE81B2C0203301B73FCE13C14FB2744AB1EBD24029062D840ECE15C161DDB50E8D6BD140EF2CF4FC0ED315C1F4792E801F82B6C0

isValid() returns true for both geometries while almost any operation fails with geos-3.2.1 and svn-trunk (r2999).

[yabo]

Minimal C++ program reproducing the TopologyException.

[strk]

Confirmed. XML test added in r3270. JTS works.

[strk]

JTS, like GEOS, fails before engaging snaprounding: found non-noded intersection between LINESTRING ( -322378.6881546058 830.7344192302965, -336712.0 -2708.0 ) and LINESTRING ( -334393.81924198254 841.2944606413985, -332495.8798798799 -1667.0870870870874 ) [ (-332495.8798798799, -1667.0870870870874, NaN) ]

[strk]

First spotted difference is on common bits computation:

GEOS:
Computed snap tolerance: 1.9759000000000001741e-05
Computed common bits: -327680 0

JTS:
Snap tol = 1.9759000000000002E-5
Computed common bits: (-262144.0, 0.0, NaN)

[strk]

I've found that even using the geometries resulting from JTS snapping (which happens) GEOS still fails, so the problem is not with snapping but after that.

Next difference is the number of intersections found between edges of the two rings. JTS finds 2 intersections, GEOS only 1:

GEOS:
Intersections:
-60234.68815460579, 830.7344192302965 seg # = 1 dist = 5568.480037827466

JTS:
Intersections:
-60234.68815460579, 830.7344192302963 seg # = 1 dist = 5568.480037827466
-60234.68815460579, 830.7344192302967 seg # = 1 dist = 5568.480037827467

[strk]

On a closer look GEOS also finds 2 intersection points, but only regsiter them on the first set of edges, not on the second:

GEOS:

Self edges:
edge  LINESTRING(-52631 2708, -74568 -2708, -85994 17051, -77080.934707903769 1637.6219931271462, -72249.819241982535 841.29446064139847, -70351.879879879882 -1667.0870870870874, -52631 2708)  A:ebi B:---  0
Intersections:
-60234.688154605792 830.73441923029645 seg # = 0 dist = 7603.6881546057921
-60234.688154605792 830.73441923029657 seg # = 5 dist = 10117.19172527409

Other edges:
edge  LINESTRING(-95427.747025205696 -5762.1230496452081, -59598.914369572361 -4737.7456185971696, -62351.928715514019 19374.77016942418, -95107.62907800317 17838.204022852125, -95427.747025205696 -5762.1230496452081)  A:--- B:ibe  0
Intersections:
-60234.688154605792 830.73441923029645 seg # = 1 dist = 5568.4800378274658


JTS:

Self Edges:
edge null: LINESTRING (-52631.0 2708.0,-74568.0 -2708.0,-85994.0 17051.0,-77080.93470790377 1637.6219931271462,-72249.81924198254 841.2944606413985,-70351.87987987988 -1667.0870870870874,-52631.0 2708.0)  A:ebi B:--- 0
Intersections:
(-60234.68815460579, 830.7344192302963, NaN) seg # = 0 dist = 7603.688154605792
(-60234.68815460579, 830.7344192302967, NaN) seg # = 5 dist = 10117.19172527409

Other Edges:
edge null: LINESTRING (-95427.7470252057 -5762.123049645208,-59598.91436957236 -4737.74561859717,-62351.92871551402 19374.77016942418,-95107.62907800317 17838.204022852125,-95427.7470252057 -5762.123049645208)  A:--- B:ibe 0
Intersections:
(-60234.68815460579, 830.7344192302963, NaN) seg # = 1 dist = 5568.480037827466
(-60234.68815460579, 830.7344192302967, NaN) seg # = 1 dist = 5568.480037827467

[strk]

"Other edges" is likely the edge of the rectangle-lookin polygon, which is crossed by the tiny spike, as suggested by the very close distance of the two intersections in "Other Edges". GEOS fails to create two nodes in there. Maybe collapsing the two intersection point as a single one.

Not even sure this would be a problem on itself, but the NodingValidator evidently does find two distinct intersections, one of which isn't found to be equal to an existing node...

[strk]

So, GEOS segment intersector properly finds both intersection points on that single segment of the rectangle. Both edges of the spike are found an intersection for. These two points are _very_ close to each other, but slitly closer in GEOS than in JTS:

GEOS:
 Adding EI -60234.688154605792 830.73441923029645 segIdx:1 dist:5568.4800378274658
 Adding EI -60234.688154605792 830.73441923029657 segIdx:1 dist:5568.4800378274658
JTS:
 Adding EI -60234.68815460579  830.7344192302963  segIdx:1 dist:5568.480037827466
 Adding EI -60234.68815460579  830.7344192302967  segIdx:1 dist:5568.480037827467

Note that this is a vertical edge, with the two intersection points being at a vertical distance of about 1E-13 (GEOS) to 4E-13 (JTS) units. The last "dist" is the distance of the intersection point from the start of the segment, which is much more far away (over 5568 units and counting). In GEOS this latest distance becomes the same for the two points, whereas in JTS you can still distinguish it (1E-12).

I'm going to check the distance computer next.

[strk]

To complete the analisys, I confirm that EdgeIntersectionList is a set usin EdgeINtersection comparator to distinguish intersections, and the comparator currently only uses segmentIndex and distance, which means in the GEOS case the second intersection would NOT be added.

I guess a possible solution could be using a long double for the distance.

[strk]

I confirm using a "long double" for the distance of an INtersection along a segment makes the operation succeed. Still the lower intersection point found by GEOS on the vertical edge of the rectangle is slighly closer to the higher one (intersections with the spike edges).

It might be worth tracking that difference out.

[strk]

On a side note, using the long double doesn't even fail with the original data, so no need to snapround.

[strk]

Another approach is, in case of same segment and same distance, look at the coordinate to compare two intersections. Problem is how to tell which comes before the other, as I guess it matters (altough no test in GEOS testsuite fails usin either comparison side for Coordinates).

[strk]

Fixed with r3334 by using long double for edge distance. Regress-tested with r3335.