Port DoubleDouble and CGAlgorithmsDD support

[mloskot]

This is a prototype of DoubleDouble support for GEOS, as an analog to JTS DoubleDouble.

The PR includes:

• ttmath 0.9.3 library
• doubledouble 2.2 library (GPL!)
• Minimal port of JTS CGAlgorithmsDD class
• DoubleDouble (DD) support can be switched to either ttmath or doubledouble at build-time, see detailed logs of the commits below.
• A bunch of curious test cases:
  • Test case of orientationIndex method for orientation vs inverse (see comment in
    RobustDeterminant.java/RobustDeterminant.orienationIndex; JTS r931 added this test in jts/algorithm/CGAlgorithmsTest.java)
  • Test cases of test on segment vs test on vertex of segment based on geometries in both formats, WKT and WKB. See commit logs for explanation of the curious results obtained using both JTS and GEOS (with DD support), for WKT vs WKB. Here is corresponding patch with those tests for JTS jts-algorithm-test-point-on-segment.patch
    .

This PR is directly related to the problems discussed in this GEOS ticket http://trac.osgeo.org/geos/ticket/591

Next, I'm going to contact @dr-jts if he would be willing to comment on the JTS results obtained for the test on segment case.

[dr-jts]

It's hard to see what the precise issue is. I assume there is some problem with the testPointOnSegment and/or the testPointOnVertex cases?

First of all, the tests use LineStrings with many points - is it possible to reduce the cases to a point against a single segment?

[mloskot]

@dr-jts I meant to write to you directly with clearer explanation.

First thing I'd like to confirm is how to interpret the point on segment test case results I've got:

1. The test passes if run by the XML tester which sources WKT-based geometries. The test files are TestIntersectsPL.xml and TestPreparedIntersectsPL.xml from the commit c3ec049
2. The test fails if run as JUnit test in Java which sources WKB-based geometries. See patch for RobustLineIntersectorTest.java.

So, the JTS test performs the same operations on the same geometries (different though at very far right precision places, after 13th digit), but depending on input format, WKT or WKB, the results are very different.

Martin, thanks for getting in touch about this.

[dr-jts]

In the XML test the 3rd segment of the LineString is vertical (identical X ordinates). The test point has the same Y ordinate, so it obviously lies on the segment, as the test reports.

It looks to me like the problem is that in the WKB geometry the 3rd segment is NOT vertical. The X ordinates are slightly different:

LINESTRING (-23.122057005538952 50.520197677479445,
-23.11534769669951 50.51334048151994,
-23.109468960005515 50.522337645220134, <--------
-23.1094689600055 50.51691776295595,
-23.096196792094247 50.53304648480942,
-23.088799100603424 50.52585152131849,
-23.08523026223618 50.52645822384092)

Obviously the point POINT(-23.1094689600055 50.5195368635957) can't lie on that segment, as the test reports.

So the question is: where did the WKB come from, and why do you think it exactly represents the WKT in the comment?

[mloskot]

That is indeed what I tried to point out in my previous comment.
On my machine, the WKB parsed with GEOS gives more differences comparing to the one you get:

POINT (-23.1094689600055080 50.5195368635957180)
--------------------^^^^^^^------------^^^^^^^^
LINESTRING 3rd and 4th points
       -23.1094689600055150 50.5223376452201340, 
       -23.1094689600055010 50.5169177629559480,
--------------------^^^^^^^------------^^^^^^^^

Use of the textual representation does not preserve the complete bits patters and doesn't help to make the tests predictable either, does it?
That is,
23.1094689600055
becomes
23.1094689600055080 or
23.1094689600055000 or
23.1094689600055001 or...

I'd argue that the coordinates in WKT are equally precise as the ones in WKB, so the 3 points are equally collinear in both.

The differences you point out seem significant(ish), but I was surprised the algorithm doesn't trim at 15th and later decimal place really.

It seems, that GEOS with DoubleDouble detects the point/segment above intersect is by a flux.

[dr-jts]

There's quite a few issues being raised here, so I'm going to try and dissect them.

First of all, the discussion of whether the point-on-segment algorithm is correct is not the issue here, as far as I can see. In the JTS WKT test, the test reports correctly that the point lies exactly on a vertical segment with the same Y ordinates (in FP representation).

In general, JTS algorithms work with FP representations of numbers. These are assumed to be stable during copying, so there's no need to round the numbers internally, because they produce consistent results.

So the real issue is the question of how stable the conversion from WKT and WKB to FP is. In JTS at least I believe both the WKT parser and the WKB parsers are stable, in the sense that identical text numbers produce identical FP bit patterns. Perhaps this isn't the case in GEOS (at least on some machines) ?

If you are faced with parsers that aren't stable, or if you would like input ordinates to be rounded to a lower precision, then you need to implement this by processing the geometries read from parsers. JTS provides the GeometryPrecisionReducer class to do this. Or, you can set the PrecisionModel on the parsers.

[dr-jts]

To put this issue another way, in the WKB for the LineString in the testPointOnSegment method has the 3rd and 4th X ordinate values are slightly different in the 15th least significant hex digit:

3: 9c266328061c37c0 56d8bff5db424940
4: 98266328061c37c0 034f7b5c2a424940
^

So the line segment is not vertical, and thus the point with WKB:

9c266328061c37c0 56d8bff5db424940

obviously does not lie exactly on the segment. So no issue there - everything is working as it is supposed to.

So not sure what the actual issue is, if there is one?

[mloskot]

On 2 December 2014 at 19:39, Martin Davis notifications@github.com wrote:

To put this issue another way, in the WKB for the LineString in the
testPointOnSegment method has the 3rd and 4th X ordinate values are
slightly different in the 15th least significant hex digit:

3: 9c266328061c37c0 56d8bff5db424940
4: 98266328061c37c0 034f7b5c2a424940
^

So the line segment is not vertical, and thus the point with WKB:

9c266328061c37c0 56d8bff5db424940

obviously does not lie exactly on the segment.

That is the difference I was pointing at as well and explained that I'd
expect trimming after 14th decimal place, otherwise even exact values may
never be float-point exact.
Nevertheless, I'm not suggesting there is a bug in JTS. just trying to
understand how it works to validate my expectation.

[mloskot]

On 2 December 2014 at 18:31, Martin Davis notifications@github.com wrote:

First of all, the discussion of whether the point-on-segment algorithm is
correct is not the issue here, as far as I can see. In the JTS WKT test,
the test reports correctly that the point lies exactly on a vertical
segment with the same Y ordinates (in FP representation).

OK, good.

So the real issue is the question of how stable the conversion from WKT
and WKB to FP is.

Yes, this is an important aspect of the issue as well.

If you are faced with parsers that aren't stable, or if you would like
input ordinates to be rounded to a lower precision, then you need to
implement this by processing the geometries read from parsers. JTS provides
the GeometryPrecisionReducer class to do this. Or, you can set the
PrecisionModel on the parsers.

Sounds like a viable alternative. I will give it a try.

Martin, thank you for your feedback, much appreciated.

Regardless all the precision issues we've discussed here, this patch with
the DoubleDouble support for GEOS is still valid and perhaps interesting to
some GEOS users. So, I'd like to keep it open. @strk what's your comment on
that, any interest to submit it (with some necessary refactoring and
updates in the build configurations, obviously).

[strk]

Does the DD code help fixing any case that would fail otherwise ?

[mloskot]

@strk Yes. The orientation index test I copied from JTS in 836e701 is one case. Not sure if enough to roll the cannon.

[strk]

I haven't checked the code, why were external libraries needed ? Doesn't JTS implement the "doubledouble" concept internally ? How can we ensure JTS and GEOS implementations are the same ? Does JTS have unit tests for the DoubleDouble class ?

About the double dependency, is that a compile-time option to use one or the other ?
Why should one be choosen over the other ?
Note that in one case the resulting binary would need to be GPL.

[mloskot]

On 3 December 2014 at 17:31, Sandro Santilli notifications@github.com
wrote:

I haven't checked the code, why were external libraries needed ? Doesn't
JTS implement the "doubledouble" concept internally ? How can we ensure JTS
and GEOS implementations are the same ? Does JTS have unit tests for the
DoubleDouble class ?

AFAIU, JTS DoubleDouble is port of C++ doubledouble I linked above.
I added ttmath for comparisons, experiments and it was just easy to use.

About the double dependency, is that a compile-time option to use one or

the other ?

Yes.

Why should one be choosen over the other ?

I haven't evaluated them to compare.
Both are very easy to use.

Note that in one case the resulting binary would need to be GPL.

That's one of the reason to choose ttmath over doubledouble.
The ttmath was easy to use and it's also adapted to Boost.Geometry, so it
was an obvious choice for me.

Perhaps there are other libraries that could be tried out and even simpler
to use.
So, wider feedback would be helpful.

[strk]

I'm happy to see this in trunk if you're up to take care of the automake part.

[strk]

@dr-jts: would it make sense to re-test existing known robustness cases filed on trac with the new DD code ?

[mloskot]

On 3 December 2014 at 17:50, Sandro Santilli notifications@github.com
wrote:

@dr-jts https://github.com/dr-jts: would it make sense to re-test
existing known robustness cases filed on trac with the new DD code ?

I was going to ask about that too.

[dr-jts]

Well, more testing is always good. That said, I'm pretty sure that DD doesn't help all that much with large-scale robustness (e.g.robustness due to issues in overlaying polygons). That's because overlay robustness failures are often due to the limited precision of computed points, which the DD tests don't help with.

The motivation in introducing DD is to finally eliminate robustness issues from core low-level predicates. This is obviously necessary to attain robustness in higher-order algorithms - but unfortunately it's not sufficient for many algorithms (such as overlay).

[mloskot]

The motivation in introducing DD is to finally eliminate robustness issues from core
low-level predicates. This is obviously necessary to attain robustness in higher-order
algorithms - but unfortunately it's not sufficient for many algorithms (such as overlay).

This is actually a very good summary and useful documentation from users' point.

Apart from testing and curiosity drive, my main motivation to for trying out the DD support was: it's in JTS, it should be in GEOS, for completeness of the port.

[mloskot]

The (in)famous point-on-segment and point-on-vertext test cases have been corrected according to @dr-jts suggestions and added to the GEOS tests suite, see http://trac.osgeo.org/geos/changeset/4038

[strk]

See also #42 as it contains a testcase that might be fixed by this work