Change Shewchuk Orientation filter to Ozaki et al. filter.

[tinko92]

This commit changes the current Shewchuk-based (notably, the original code doesn't use Shewchuk's constant which is significantly narrower) to the filter published by Ozaki et al. in https://doi.org/10.1007/s10543-015-0574-9 . It is tighter, so it has strictly fewer misses and it has fewer and better predictable branches and fewer instructions so it is theoretically and practically faster (the referenced paper has benchmarks).

It does not come with new tests because it should not change any observable behavior by design except possibly speed, which is not tested, so all tests testing the current predicate are tests for the patched implementation.

[dr-jts]

Where does the value 3.3306690621773714e-16 come from? I don't see it in the Ozaki paper (source from here).

[dr-jts]

See also libgeos/geos#1184

[tinko92]

Where does the value 3.3306690621773714e-16 come from? I don't see it in the Ozaki paper (source from here).

The algorithm is algorithm 3. u_N is an underflow guard to make the errorbound no zero in case of underflow, but I omitted it since I don't think the exact fallback computation handles underflow cases correctly, would need an extended exponent for that. The constant is theta, which is given in Theorem 3.1, page 10, the phi is given for double precision in Lemma 3.1 page 8. u is the roundoff unit for double, so 2^-53, given on page 4. Putting it all together should yield the constant with 1516 digits of precision for faithful double representation. The constant may not be sharp (never saw an input with an error that came very close to it and had the wrong sign) but it is almost sharp for one of the intermediate computations.

Edit: I only ever computed this value with C++, it's 3u-(phi-22)u^2 with phi=94906264 and u=2^(-53), assuming the computation is exact. But now that you mention it, 22 is not divisible by 4 but phi is, so it does get rounded to even and the nearest even number seems actually lower, so it has to be 3.3306690621773724E-16, not 3.3306690621773714E-16. I'll change the commit and thanks for double-checking.

[dr-jts]

I wonder if this solves #750?

Specifically, what is the orientation index computed for

LINESTRING (0 100, 3 106.3246)

POINT (1 102.1082)

Currently this is computed as -1, but by inspection it should be 0 (collinear).

If so, this case can be added as a unit test.

[tinko92]

I wonder if this solves #750?

Specifically, what is the orientation index computed for

LINESTRING (0 100, 3 106.3246)

POINT (1 102.1082)

Currently this is computed as -1, but by inspection it should be 0 (collinear).

If so, this case can be added as a unit test.

-1 is correct for this test case if the numbers are read as "nearest double to 106.3246" and "nearest double to 102.1082", which is how a compiler or number parser will read them, here is a visualization of the double values close to that:

This is like a 2D version of testing 0.2 + 0.1 == 0.3. The black circle marks where the point is.

Edit: Here is the tool for references, I know the UI isn't labelled well: https://www.tinkobartels.de/failure_viewer/ . It is for looking at different predicate filter failures. Selecting the last predicate in the list will guarantee correct results, when zoomed in close enough a grid of neighbouring double values will be displayed.

Blue represents -1, Red represents 1, Green represents 0. Notably, there is no green in this picture because the exact colinear point at x=1 is not representable in double precision. When picking two points with random decimal fractions that don't happen to be a sum of low powers of 1/2 and rounding them to double, almost none of the collinear points exist in the set of double numbers. The pixels are not quadratic because at 1 and 100, the exponents of the x and y coordinates are different. The other pattern in the grid is just a visual artifact of me not exactly aligning pixels in the canvas with actual screen pixels.

I can see how this is an unexpected result for library users, though. I see three possible solutions to solve this:

• Use this exact predicate for construction to guarantee algorithmic correctness and anywhere else to guarantee consistency at the cost of unintuitive results like this.
• Use another rounding predicate just for simple point to linear tests to be more in line with user expectation where inconsistent individual predicate results cannot mess up the consistency of a more complicated datastructure like a computed intersection polynomial or a triangulation. In this case in the filter failure case 0 should be returned and the constant should be changed from ~3e-16 to roughly ~5e-16 (can try to compute the exact constant later) and the expression needs to be changed slightly to guarantee that including the initial rounding error from a decimal input, "collinear up to rounding" points are recognised as such after the accumulated error of the determinant computation. This would also find some points collinear whose decimal representations are not collinear and will produce subtly inconsistent results between different algorithm.
• Use a decimal number type and extended precision arithmetic at the cost of a severe performance penalty (I'd expect at least 1 order of magnitude slower operations for everything where predicate computation cost is not negligible like triangulations).

[dr-jts]

@tinko92 Great analysis, thanks. It sounds like there is no easy/fast solution to this problem. That's not a surprise, given the well-known issues of geometric robustness when using floating-point. We'll just have document this situation, and ideally provide tolerance-based algorithms as a work-around.