search

locationtech / jts

The JTS Topology Suite is a Java library for creating and manipulating vector geometry.
Other
1.99k stars 443 forks source link

Change Shewchuk Orientation filter to Ozaki et al. filter. #1093

Open tinko92 opened 3 weeks ago

tinko92 commented 3 weeks ago

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 commented 3 weeks ago

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

dr-jts commented 3 weeks ago

See also https://github.com/libgeos/geos/pull/1184

tinko92 commented 3 weeks ago

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 ~15~16 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.