robust-predicates

Exact-sign geometric predicates for Zig: orientation and incircle / insphere tests, ported from Jonathan Shewchuk's adaptive-precision arithmetic. No dependencies, no allocations, thread-safe.

Floating-point rounding can make a naive orientation test output a wrong sign. These predicates return the exact sign, while doing only as much work as the input requires. The easy cases stay cheap, the degenerate ones get more expensive as more work is required to keep the result exact.

Predicates

The sign convention is the one in Shewchuk's paper: positive means counterclockwise or inside.

Each predicate has a fast* variant (fastOrient2D, fastIncircle, …). Those are fast, non-exact, non-adaptive. They are faster but not robust, so the sign can be wrong on hard cases.

Available predicates are:

• orient2D(a, b, c): orientation of triangle a, b, c. Positive if counterclockwise, negative if clockwise, zero if collinear. Magnitude approximates twice the signed triangle area.
• orient3D(a, b, c, d): orientation of d relative to the plane through a, b, c. Positive if d is below the plane (with a, b, c counterclockwise seen from above), negative if above, zero if coplanar. Magnitude approximates six times the signed tetrahedron volume.
• incircle(a, b, c, d): tests if d is inside the circle through a, b, c. Positive if inside, negative if outside, zero if cocircular. a, b, c must be counterclockwise, or the sign is reversed.
• insphere(a, b, c, d, e): tests if e is inside the sphere through a, b, c, d. Positive if inside, negative if outside, zero if cospherical. a, b, c, d must be positively oriented (per orient3D), or the sign is reversed.

Install

zig fetch --save git+https://github.com/woldendans/zig-robust-predicates#v0.1.1

Then add the module in build.zig:

const rp = b.dependency("robust_predicates", .{
    .target = target,
    .optimize = optimize,
});
exe.root_module.addImport("robust_predicates", rp.module("robust_predicates"));

Usage

const rp = @import("robust_predicates");
const RP = rp.F64; // f64, use rp.F32 or rp.RobustPredicates(...) for another float

const a = RP.Point2D{ .x = 0, .y = 0 };
const b = RP.Point2D{ .x = 1, .y = 0 };
const c = RP.Point2D{ .x = 0, .y = 1 };

const s = RP.orient2D(a, b, c); // > 0 (CCW)

Requires Zig 0.16.0 or later.

Performance

To leverage modern CPU capabilities and simplify the original implementation this library's twoProduct tries to use @mulAdd (fused multiply add, one rounding). In case the library is compiled targeting non-FMA ready hardware it will fall back to a Dekker/Shewchuk split-based implementation. This is a comptime decision.

When using the split version, for simplicity's sake, the pre-split optimizations of the original paper have not been implemented. After some tests, it seems modern compilers are very capable of optimizing the code without them anyway, bringing the potential gain to a very low value.

Unless you require targeting 10+ years old hardware or a specific architecture, it is advised to build for a FMA-ready CPU target. Examples:

• x86-64: -Dcpu=baseline+fma or -Dcpu=x86_64_v3
• AArch64 / arm64: nothing to do

License

Released under the MIT license, see LICENSE.

Credits

This library is a Zig port of robust-predicates (https://github.com/mourner/robust-predicates) by Vladimir Agafonkin, released under the Unlicense (public domain), itself a port of the public-domain C routines by Jonathan Richard Shewchuk (https://www.cs.cmu.edu/~quake/robust.html).