Port DoubleDouble and CGAlgorithmsDD support #40

Closed
wants to merge 5 commits into
from

Projects

None yet

3 participants

@mloskot
Member
mloskot commented Dec 1, 2014

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:

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.

mloskot added some commits Dec 1, 2014
@mloskot mloskot Add ttmath 0.9.3 from http://www.ttmath.org
Configure ttmath with CMake.
ttmath is header-only library, but for Visual C++ 64-bit compiler
requires compilation of ttmathuint_x86_64_msvc.asm file (Microsoft
doesn't provide the inline assembler in its 64bit compiler).
ttmath will be used as Java/JTS DoubleDouble analog of
quadruple precision float-point.
59a02a5
@mloskot mloskot Add doubledouble 2.2 from http://www.boutell.com/fracster-src/doubled…
…ouble/doubledouble.html

NOTE: doubledouble is GPL
Apply number of changes and fixes to doubledouble 2.2:
1. Rename .cc file extension to .cpp
Add drand48() and srand48(long) for non-POSIX platforms from
https://gist.github.com/mortennobel/8665258, based on Freebsd
http://fxr.watson.org/fxr/ident?v=FREEBSD-LIBC;im=bigexcerpts;i=_dorand48
2. Correct names of C++ std library headers.
3. Add missing std:: namespace qualifier.
4. Visual C++ does not provide nan.h header.
4883a60
@mloskot mloskot Port CGAlgorithmsDD from JTS 1.14 (r789).
Two methods have been added orientationIndex and orientationIndexFilter.
The implementation uses either ttmath or doubledouble libraries as the
JTS DoubleDouble implementation analog.
In CGAlgorithmsDD.cpp, two defines provided to control which one is used:
GEOS_ENABLE_DD_TMNATH and GEOS_ENABLE_DD_DOUBLEDOUBLE.
(Uncomment one or the other and rebuild GEOS.)
Change CGalgorithms::orientationIndex to call CGAlgorithmsDD
(Note, previous calculations have been kept for comparison and testing
purposes, see TODO comments.)
5e5cff8
@mloskot mloskot Port curious JTS r626 test case for CGAlgorithms{DD}::orientationIndex.
The test case comes from comment in JTS r626 file
jts/java/src/com/vividsolutions/jts/algorithm/RobustDeterminant.java
(see http://sourceforge.net/p/jts-topo-suite/code/626/)

IMPORTANT: Both, JTS as well as GEOS with CGAlgorithmsDD support (through
either  ttmath or doubledouble), pass this orientationIndex test case.
See test<3> in computeOrientationTest.cpp.
836e701
@mloskot mloskot Add two test cases: point on segment and point on vertex.
Both tests have been added in two forms:
1) C++ version based on geometries in binary (WKB) format.
2) XML version based on geometry in textual (WKT) format.

IMPORTANT: Both, JTS 1.14 and GEOS with DD support (either through ttmath,
or doubledouble), fail on the point on segment test case, but pass on
the point on vertex test case.
The quadruple precision float-point support (DD) does not help in the
point on segment test case, so both, JTS and GEOS, seem to
have robustness issues algorithmically.
c3ec049
@mloskot mloskot self-assigned this Dec 1, 2014
@mloskot mloskot added a commit to mloskot/libgeos that referenced this pull request Dec 1, 2014
@mloskot mloskot Improve points co-linearity test in LineIntersector::hasIntersection.
Add simple test based on area of triangle created by three points
of the given segment and test point.
The point intersects the segment, if area of the triangle is
zero or close to zero: insignificant relatively to length of the segment.

Two interesting observations:
1) The fix allows to pass the "point on segment" test discussed in
http://trac.osgeo.org/geos/ticket/591 and PR: libgeos/libgeos#40
2) All GEOS tests pass apart from (testA, test<13>):
tests/unit/algorithm/RobustLineIntersectorTest.cpp:247
which is opposite to the "point on segment" test:
segment located as close and parallel as possible to X axis,
so it seems on X axis (so 0,0 intersects), but
JTS/GEOS are expected to detect it's not on X axis.
65b281b
@dr-jts
dr-jts commented Dec 1, 2014

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
Member
mloskot commented Dec 1, 2014

@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
dr-jts commented Dec 1, 2014

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
Member
mloskot commented Dec 2, 2014

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
dr-jts commented Dec 2, 2014

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
dr-jts commented Dec 2, 2014

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
Member
mloskot commented Dec 3, 2014

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
Member
mloskot commented Dec 3, 2014

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
Member
strk commented Dec 3, 2014

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

@mloskot
Member
mloskot commented Dec 3, 2014

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

@strk
Member
strk commented Dec 3, 2014

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
Member
mloskot commented Dec 3, 2014

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
Member
strk commented Dec 3, 2014

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

@strk
Member
strk commented Dec 3, 2014

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

@mloskot
Member
mloskot commented Dec 3, 2014

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
dr-jts commented Dec 3, 2014

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
Member
mloskot commented Dec 3, 2014

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 mloskot added a commit that referenced this pull request Dec 5, 2014
@mloskot mloskot Add two test cases, point-on-segment and point-on-vertex.
Curious detail of the tests is that points of interest have nearly exact X of tested points, the values differ after 14th decimal place.
The tests provided test geometries for intersection with and without coordinates trimming after the 14th place (as per Martin Davis suggestion).
It has been extensively discussed in Ticket #591 and #40 with Martin Davis' input.


git-svn-id: http://svn.osgeo.org/geos/trunk@4038 5242fede-7e19-0410-aef8-94bd7d2200fb
7196b9a
@mloskot
Member
mloskot commented Dec 5, 2014

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
Member
strk commented Jan 16, 2015

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

@mloskot mloskot closed this Aug 5, 2015
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment