Use Ozaki et al.'s error bound and single-branch evaluation in orientation index filter.

[tinko92]

This PR is separated out of #1162 . It contains two changes to the orientation predicate.

The first change changes the relative error bound from currently 1e-15 to ~3.33e-16, which is proven to be sufficient in https://doi.org/10.1007/s10543-015-0574-9 by Ozaki et al. A tighter error bound can reduce the calls to the expensive exact computation in some degenerate cases (Edit: Here is a gist demonstrating this for a sequence of 100k generated near-collinear inputs https://gist.github.com/tinko92/ce44f93d8d4e4d5909310c00cf19b391 ) and has no disadvantages.

The second changes is also motivated by https://doi.org/10.1007/s10543-015-0574-9 and reduces the number of branches and instructions in the evaluation. This is more concise and may result in faster evaluation (benchmarks are in the referenced paper) but I am not currently aware of any geos-benchmarks, beside the Orientation Index micro benchmark, for which the orientation predicate is the bottle neck. The original logic that is replaced is from the function body of REAL orient2d(pa, pb, pc) https://www.cs.cmu.edu/afs/cs/project/quake/public/code/predicates.c , where the early returns make use of the branching in the fabs-computation, discussed in the comment preceding #define Absolute(a) ((a) >= 0.0 ? (a) : -(a)). This reasoning is obsolete because modern compilers with enabled optimizations will compile fabs to a branchless masking of the sign bit.

Edit: I did not remove the orientation(double x) function, which seems unused now, because I don't know if that would be considered an API change (it's not in a namespace called detail or internal or sth. like that and I don't know the geos conventions). But I can edit the commit if it is no longer needed.

Edit2: Sorry for only now reading more about the jts-geos-relationship/-process. I will propose the same change to JTS, when I find time. Having said that, since patched and unpatched always return the correct orientation by deferring filter misses to the exact computation for inputs that do not go near overflow/underflow, this would not change observable behavior between JTS/Geos for the same inputs, except in computation time, so I think it classifies as an implementation detail/optimization.

[dr-jts]

Edit2: Sorry for only now reading more about the jts-geos-relationship/-process. I will propose the same change to JTS, when I find time. Having said that, since patched and unpatched always return the correct orientation by deferring filter misses to the exact computation for inputs that do not go near overflow/underflow, this would not change observable behavior between JTS/Geos for the same inputs, except in computation time, so I think it classifies as an implementation detail/optimization.

Thanks for providing a PR for JTS well. It's correct that the two libraries sometimes diverge in implementation, but it's nice to keep them in synch as far as possible, to reduce the mental complexity of future debugging.

[dr-jts]

@tinko92 it would be nice to have a stress test in GEOS (and JTS) for point-line orientation, to confirm that the implementation is working as expected. Do you know of anything like this?

The tricky part of course is to generate test cases with known correct answers. This could be done either a priori, by careful construction of the cases, or via some kind of higher-precision (but presumably slower algorithm).

[tinko92]

@dr-jts The implementation here changes the filter. The filter's correctness criterion is to return either the correct sign or failure so an error case would be returning the wrong sign, returning failure would never be incorrect (though, potentially a performance regression, should never occur with this filter without also occuring with the old one). The entire predicate runs an exact stage iff failure is returned, so any error case of the filter is an error case of the entire predicate so testing the entire predicate is sufficient (and strictly better coverage). The degenerate CCW tests currently in the test suite do some of that.

To stresstest the implementation I would run many degenerate cases (so maybe not as a unit test). Off the top of my head, I would suggest the following generators of degenerate cases:

1. Most straightforward: Two random, non-equal double-coordinate points pa, pb in [-1., 1.]×[-1., 1.] and a random scalar c in [-1., 1.], then set pc = pa + c * (pb - pa)

2. Most realistic (imho) representation of likely real-world degeneracies: Same thing with integers in [-1000, 1000]×[-1000, 1000], then divide all coordinates by 1000.0 to get double coordinates. Mirrors the degeneracies that come from rounded decimal input.
   Both test the usual rounding error with numbers at the same scale.

3. Random pa in [-1., 1.]×[-1., 1.], random scalars b and c in [-100., 100.], set pb = pa * exp(b), pc = pa * exp(c). Tests rounding-errors from catastrophic cancellation.

4. Same as above but c in [-700., 700.]. Will test overflows and underflows. The filter will be correct for overflows, some times incorrect for underflows (intentional omission because it costs some cycles and filter failures and is not covered by the exact stage in geos anyways as far as I can see).

Potential oracles to test the results against are:

1. Shewchuk's predicate.c (very mature, last revised in 1996, C++ port with minimal syntax updates available on github), can deal with 1. to 3.
2. The Orientation_2 method of Epick or Epeck kernel from CGAL, can deal with 1. to 4.

Stresstest without an oracle: The generators can be extended to return four points instead of three. The predicate would then not be tested for correctness but for consistency, i.e. something roughly like

for ( auto [p1, p2, p3, p4] : generator() ) for ( auto [pa, pb, pc, pd] : permutations(p1, p2, p3, p4) )
if ( (pa == pb || pa == pc || pb == pc) )
{
	if (orientation(pa, pb, pc) != 0) throw ...;
	else continue;
}
if ( orientation(pa, pb, pc) != -orientation(pa, pc, pb) ) throw ...;
if ( orientation(pa, pb, pc) == 0 )
{
	auto [pa_, pb_, pc_] = sort_lexicographically(pa, pb, pc);
	if ( orientation(pa_, pb_, pd) != orientation(pa_, pc_, pd) ) throw ...;
	else continue;
}
auto [pa_, pb_, pc_] = sort_to_ccw(pa, pb, pc);
if(    orientation(pa_, pb_, pd) == right
    && orientation(pb_, pc_, pd) == right
    && orientation(pc_, pa_, pd) == right ) throw ...;
}

I hope, I did not miss anything there. Any correct predicate will be consistent. I would be very surprised if any incorrect predicate (except something trivial like always returning the opposite of the correct sign or always 0) could realistically be consistent. This has the advantage of not needing external code and consistency is the property one actually cares about.

If you would like some or all of the above in geos (generators 1-4, test for correctness with external dependency or consistency without or both), I could write that as a PR later when I have time (maybe for the benchmark directory? Because I think to be useful this would be long-running and should not run by default through make check).

If test cases should be generated together with results from oracles and then stored in text, it would be very important to print the coordinates after std::cout << std::setprecision(std::numeric_limits::max_digits10);

[dr-jts]

Thanks for the detailed reply. I agree that for the orientation predicate consistency is at least as important as correctness. Your suggestion of a consistency checker is very interesting. A PR would be welcome.

On a related note, I had hoped that the relaxed filter bound would resolve this issue, where the DD orientation result is inconsistent with the FP result, and from a human perspective obviously incorrect:

bin/geosop -a "LINESTRING (2 0, 0 2)" -b "POINT(0.4 1.6)" intersection

But this still produces the visibly incorrect answer of POINT EMPTY - i.e. no intersection.

I've been wondering about always running the orientation computation first in FP, and accepting that answer if it is 0 (point lies on the line). This should avoid inconsistencies such as the above. @tinko92 any thoughts on this?

[tinko92]

@dr-jts For the given example, the filter returns failure (uncertainty) and the exact computation returns non-collinearity, which is correct for the closest floating-point approximations of (rd(0.4) rd(1.6)), where rd(x) is the function that rounds to nearest double.

Behaviour under different orientation filters and exact orientation predicate

A point that is collinear (green) with (2 0, 0 2) and has the y-coordinate rd(1.6) exists but its x-coordinate is slightly smaller than rd(0.4). In the following, I visualised the area of filter failures (yellow) respectively for Ozaki's filter, an interval arithmetic filter (slower) and exact computation (much slower), with Point(rd(0.4), rd(1.6)) marked by a black dot.

Filter as orientation predicate with tolerance

You could change the filter's return value for uncertainty cases from failure (trigger DD computation) to collinear (then DD computation would just never be called) to obtain a predicate with tolerance. This predicate will return collinear for many but not all inputs where decimal approximation would suggest collinearity. To obtain a tolerant filter that guarantees returning collinearity on inputs that are collinear in decimial, you also need to account for rounding of the input values (which Ozaki's filter does not do) and you obtain an error bound of the form ~6e-16 ( (abs(a.x) + abs(c.x)) * (abs(b.y) + abx(c.y)) + (abs(a.y) + abs(c.y)) * (abs(b.x) + abx(c.x)) ), iirc. This would be similar to Burnikel et al.'s filter, with a much larger tolerance region but still on the relative order of 10 times machine epsilon (canvas zoomed out compared to above ones).

New problems/side effects expected from that approach

This approach would definitely not be consistent for more than three points, e.g. pc might be judged not-collinear to line(pa, pb) under tolerance but for some pd that is far away from pa, pb and pc and thereby induces a much larger tolerance area, both pb and pc could be judged collinear to line (pa, pd) under tolerance. The approach of rounding determinants smaller than some epsilon-dependent bound to zero is discussed by Kettner et al. as "epsilon-tweaking" in Section 6 Non-Solutions.

Recommendation and Warning

If you want to use such an orientation test with tolerance I would exclusively use it for simple predicate algorithms (like point-in-linestring predicate, etc.) but never for construction algorithms that rely on consistency in non-trival ways (like triangulation, convex hull, polygon clipping etc.), due to the discussion by Kettner et al., and still expect some non-intuitive situations, e.g. a user computes an intersection point of line a and b, which has a larger accumulated error than you would get from just rounded decimal inputs of an actual intersection point, the predicate with tolerance judges the computed intersection point non-collinear with lines a and or b. The user then prints the computed intersection point to stdout with default precision (like 6 digits?) that masks the accumulated error from the intersection computation, tests the point directly against the line again and then the predicate suddenly returns collinearity.

Potential rigorous (but difficult) alternative without above problems

I believe (I am not aware of any implementation), a geometric predicate that is 1) consistent, 2) returns intuitive results for decimal inputs, and 3) not orders of magnitude slower than what you have now could work like this:

• take coordinate input as pair of strings, store double approximation and pointer to pair of decimals in each point (50% memory overhead in point for extra pointer + extra memory for decimal number types somewhere else)
• Use filter similar to Burnikel's (maybe ~10-20% slower than Ozaki's filter?)
• On filter failure, use exact decimal computation (very slow but rarely called).

Extending this to constructions (e.g. also computing intersection points as approximations and in decimal if needed for set operations) would be a major undertaking, though.

[dr-jts]

@tinko92 once again, thanks for your detailed and informative comment.

Agreed that any approximated orientation test should only be used for spatial predicates (where topology is only evaluated at points, so overall consistency is less (or non) critical.

I'm still unclear whether it is reasonable to use an approach of first evaluating in FP, and accepting the answer if it is exact (point on line)? This seems like it would take care of "obvious" situations such as the one above. As you say, it probably won't be consistent, and perhaps even "obviously wrong" for larger inputs (i.e. ones which are collinear using exact arithmetic on the decimal input, but for which the FP evaluation returns non-collinear due to round-off from either decimal-to-FP conversion or determinant computation).

But maybe this doesn't really matter. Users should understand the limitations of finite machine arithmetic, and not expect infinite accuracy. (This would be more palatable if the topological predicates supported a precision model - which is technically feasible, and will be implemented in the future, I hope).

[dr-jts]

For further inspiration for future work on this, this is a nice C++ implementation of Shewchuk. It has an Orient2D consistency test using perturbation (although GEOS has many more test cases).

[tinko92]

of first evaluating in FP, and accepting the answer if it is exact (point on line)? This seems like it would take care of "obvious" situations such as the one above

@dr-jts , If by that you mean if(det == 0) return 0; somewhere in orientationIndexFilter then I don't think that would be helpful. It does not solve the 'obvious' case 2 0 0 2 0.4 1.6. Shewchuk's code will also not return 0 for that. I would instead recommend that if you want 0 for this case and similar cases to do if(std::abs(det) < error) return 0. Reasoning below.

Demo

Consider:

#include <iostream>
#include <string>
double orientation(double ax, double ay, double bx, double by, double cx, double cy)
{
    double const detleft  = (ax - cx) * (by - cy);
    double const detright = (ay - cy) * (bx - cx);
    double det = detleft - detright;
    return det;
};
int main(int argc, char *argv[]) {
  std::cout << orientation(std::stod(argv[1]), std::stod(argv[2]), std::stod(argv[3]), std::stod(argv[4]), std::stod(argv[5]), std::stod(argv[6])) << '\n';
  return 0;
}

On my ARM machine (you can get similar chaos on x86, discussed below):

clang++ main.cpp -O2 && ./a.out 2 0 0 2 0.4 1.6
-2.22045e-16
g++ main.cpp -O2 && ./a.out 2 0 0 2 0.4 1.6
-2.30926e-16

In neither case det == 0. Also, note the different results from two compilers with the same (generally considered safe) flag -O2 on the same machine. The difference is

GCC

        fsub    d12, d12, d0
        fsub    d15, d15, d11
        fsub    d14, d14, d11
        fsub    d0, d13, d0
        //...
        fmul    d15, d15, d12
        fnmsub  d0, d0, d14, d15 //<-- basically std::fma( (ax - cx), (by - cy), -detright);

Clang

        fsub    d1, d8, d12
        fsub    d2, d11, d0
        fsub    d0, d9, d0
        fsub    d3, d10, d12
        fmul    d1, d1, d2
        fmul    d0, d3, d0
        fsub    d0, d1, d0 //<-- doesn't insert FMA on O2

Even for a simpler symmetric case:

g++ main.cpp -O2 && ./a.out -0.1 -0.1 0.1 0.1 0 0
8.32667e-19
g++ main.cpp -O2 && ./a.out 0.1 0.1 -0.1 -0.1 0 0
8.32667e-19
clang++ main.cpp -O2 && ./a.out 0.1 0.1 -0.1 -0.1 0 0
0

Note that this is even self-contradictory between the same three points on GCC, switching two points does not invert the sign.

You'll get the same discrepancies on x86-64 with GCC between march=x86-64-v2 (what RHEL 9 is using afaik, will get you the results I get for Clang) and -march=x86-64-v3 (what they're exploring for possibly using in RHEL 10 afaik, if used with GCC on x86 it should get you the results that I get for GCC on my machine). Worse, the orientation code is a good candidate for inlining, depending on context, I could see GCC possibly deciding to put in fma in some call sites and not in others, leading to different behaviours of the same code in the build.

Behaviour for all common cases of collinearities

For the sake of clarity, consider the following formalisation: $a, b, c \in R^2$ real points, $\tilde a, \tilde b, \tilde c \in F^2$ the corresponding floating point points, $f$ the true orientation function (real arithmetic), $\tilde f$ the FP orientation (function above, FP arithmetic), $f_\text{filter}$ the filtered orientation index (FP arithmetic with error bound) in this PR and $f_\text{Shewchuk}$ the orient2d you just linked.

Always: $f(\tilde a, \tilde b, \tilde c) = f_\text{Shewchuk}(\tilde a, \tilde b, \tilde c)$.

Case 1: Trivial

• Condition: a = b or a = c or b = c or a.x = b.x = c.x or a.y = b.y = c.y
• $f(a,b,c) = 0$, $f_\text{Shewchuk} (\tilde a, \tilde b, \tilde c) = 0$
• $\tilde f (\tilde a, \tilde b, \tilde c) = 0$, $f_\text{filter}(\tilde a, \tilde b, \tilde c) = 0$.
• => No problems.

Case 2: Non-Trivial but power of two-denominator coordinates

• Condition: Not trivial but $f(a,b,c) = 0$ and $a = \tilde a, b = \tilde b, c = \tilde c$. This is the case for numbers like integers or fractions with power of two denominator $0.5$, $0.375$, etc.
• $f_\text{Shewchuk}(\tilde a, \tilde b, \tilde c) = 0$
• probably $\tilde f(\tilde a, \tilde b, \tilde c) = 0$
• almost certainly (with another filter: certainly) $f_\text{filter}(\tilde a, \tilde b, \tilde c) = \text{failure}$.

Case 3: General Case (like 2 0 0 2 0.4 1.6)

• Condition: Not trivial and $f(a, b, c) = 0$ but $(a, b, c) \neq (\tilde a, \tilde b, \tilde c)$, so at least one point is not representable in double, e.g. -0.1 -0.1 0 0 0.1 0.1 or 2 0 0 2 0.4 1.6. This is the case for by far most decimal coordinates like $0.1, 0.2, 0.3, 0.4, 0.6, 0.7, ...$, etc.
• Very likely $f(\tilde a, \tilde b, \tilde c) \neq 0$, so $f_\text{Shewchuk}(\tilde a, \tilde b, \tilde c) \neq 0$.
• Sign of $\tilde f(\tilde a, \tilde b, \tilde c)$ is basically pseudorandom in that case (as shown in the above example)
• But almost certainly (with another filter: certainly) $f_\text{filter}(\tilde a, \tilde b, \tilde c) = \text{failure}$.

$a = (2\ 0), b = (0\ 2), c(0.4\ 1.6)$ is an example of the last kind of case $a = \tilde a, b = \tilde b$ but $c \neq \tilde c$ 0.4 is the infinite binary fraction $0.011011011..._2$ and it gets rounded after 54 digits to be stored in a double, thatswhy point(0.4 1.6) for some type with constructor point(double, double) is truly not collinear.

Rounding determinants to zero if their absolute values are smaller than some carefully constructed error bound (i.e. filter failure) is the only heuristic that will handle all three cases well if you want a predicate that returns 0 for the double points representing some real points a, b, c with $f(a,b,c) = 0$.

[dr-jts]

@tinko92 I added a stress test to check orientation index consistency in #1237. Looks to me like it's showing the expected results, and confirms the robustness of the filter and DD implementation.