axiolid_predicates/
sphere.rs

1//! `incircle` and `insphere`: is a point inside a circumscribed ball?
2//!
3//! These are the Delaunay predicates. `incircle(a, b, c, d)` asks whether `d`
4//! lies inside the circle through `a`, `b`, `c`; `insphere` is the 3D analogue.
5//! A zero means `d` lies exactly on the ball -- the cocircular/cospherical case
6//! that makes a Delaunay triangulation non-unique and, if misjudged, produces
7//! inverted or overlapping cells.
8//!
9//! Both are lifted determinants: adding a coordinate equal to the squared
10//! distance from the origin turns "inside a ball" into "below a hyperplane".
11//! That lift squares the operand magnitudes, so the error bound grows faster
12//! than `orient*`'s and the filter fails sooner -- which is why the escalation
13//! rate is measured rather than assumed.
14
15use axiolid_core::{Point2, Point3};
16use axiolid_guarantees::{Certified, Precision, Sign};
17
18use crate::arithmetic::{expansion_product, expansion_sign, expansion_sum, negate_expansion};
19use crate::expansion::two_diff;
20
21/// Machine epsilon for binary64.
22const EPSILON: f64 = f64::EPSILON / 2.0;
23
24/// Relative error bound for the lifted 3x3 `incircle` determinant.
25const INCIRCLE_ERROR_FACTOR: f64 = (10.0 + 96.0 * EPSILON) * EPSILON;
26
27/// Relative error bound for the lifted 4x4 `insphere` determinant.
28const INSPHERE_ERROR_FACTOR: f64 = (16.0 + 224.0 * EPSILON) * EPSILON;
29
30/// Is `d` inside the circle through `a`, `b`, `c`?
31///
32/// [`Sign::Positive`] means inside when `a, b, c` are counter-clockwise.
33/// Callers that cannot guarantee that orientation must normalise it first with
34/// `orient2d`, because the sign of this determinant flips with it.
35///
36/// Always [`Certified::Certain`].
37#[must_use]
38pub fn incircle(a: Point2, b: Point2, c: Point2, d: Point2) -> Certified {
39    match incircle_filter(a, b, c, d) {
40        Certified::Certain { sign, .. } => Certified::exact_sign(sign),
41        _ => Certified::exact_sign(incircle_exact(a, b, c, d)),
42    }
43}
44
45/// The fast filter alone, exposed so escalation can be measured.
46#[must_use]
47pub fn incircle_filter(a: Point2, b: Point2, c: Point2, d: Point2) -> Certified {
48    let (adx, ady) = (a.x - d.x, a.y - d.y);
49    let (bdx, bdy) = (b.x - d.x, b.y - d.y);
50    let (cdx, cdy) = (c.x - d.x, c.y - d.y);
51
52    let bdxcdy = bdx * cdy;
53    let cdxbdy = cdx * bdy;
54    let alift = adx * adx + ady * ady;
55
56    let cdxady = cdx * ady;
57    let adxcdy = adx * cdy;
58    let blift = bdx * bdx + bdy * bdy;
59
60    let adxbdy = adx * bdy;
61    let bdxady = bdx * ady;
62    let clift = cdx * cdx + cdy * cdy;
63
64    let determinant =
65        alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) + clift * (adxbdy - bdxady);
66
67    let permanent = (bdxcdy.abs() + cdxbdy.abs()) * alift
68        + (cdxady.abs() + adxcdy.abs()) * blift
69        + (adxbdy.abs() + bdxady.abs()) * clift;
70
71    Certified::from_filter(
72        determinant,
73        INCIRCLE_ERROR_FACTOR * permanent,
74        Precision::F64,
75    )
76}
77
78/// Exact sign of the lifted `incircle` determinant.
79///
80/// Every coordinate difference is kept exactly, as a two-term expansion
81/// (`two_diff`), and every product and sum after it is an expansion. The
82/// differences used to be rounded first, which made this "exact" fallback
83/// wrong exactly where the filter hands over -- nearly cocircular points
84/// whose differences do not fit an `f64` -- and Delaunay flips driven by it
85/// could cycle for ever (#190).
86#[must_use]
87fn incircle_exact(a: Point2, b: Point2, c: Point2, d: Point2) -> Sign {
88    let (adx, ady) = (diff(a.x, d.x), diff(a.y, d.y));
89    let (bdx, bdy) = (diff(b.x, d.x), diff(b.y, d.y));
90    let (cdx, cdy) = (diff(c.x, d.x), diff(c.y, d.y));
91
92    let bc = minor(&bdx, &cdy, &cdx, &bdy);
93    let ca = minor(&cdx, &ady, &adx, &cdy);
94    let ab = minor(&adx, &bdy, &bdx, &ady);
95
96    let total = expansion_sum(
97        &expansion_sum(&lift(&bc, &adx, &ady), &lift(&ca, &bdx, &bdy)),
98        &lift(&ab, &cdx, &cdy),
99    );
100    expansion_sign(&total)
101}
102
103/// `p - q`, exactly, as an expansion.
104#[must_use]
105fn diff(p: f64, q: f64) -> Vec<f64> {
106    let (d, err) = two_diff(p, q);
107    let mut e = Vec::with_capacity(2);
108    if err != 0.0 {
109        e.push(err);
110    }
111    if d != 0.0 || e.is_empty() {
112        e.push(d);
113    }
114    e
115}
116
117/// `p q - r s` over expansions, exactly.
118#[must_use]
119fn minor(p: &[f64], q: &[f64], r: &[f64], s: &[f64]) -> Vec<f64> {
120    expansion_sum(
121        &expansion_product(p, q),
122        &negate_expansion(&expansion_product(r, s)),
123    )
124}
125
126/// Multiply an expansion by `x*x + y*y`, exactly.
127#[must_use]
128fn lift(e: &[f64], x: &[f64], y: &[f64]) -> Vec<f64> {
129    let square = expansion_sum(&expansion_product(x, x), &expansion_product(y, y));
130    expansion_product(e, &square)
131}
132
133/// Is `e` inside the sphere through `a`, `b`, `c`, `d`?
134///
135/// [`Sign::Positive`] means inside when `a, b, c, d` are positively oriented
136/// (`orient3d(a, b, c, d) > 0`). As with [`incircle`], the sign flips with the
137/// base orientation, so a caller must normalise it.
138///
139/// Always [`Certified::Certain`].
140#[must_use]
141pub fn insphere(a: Point3, b: Point3, c: Point3, d: Point3, e: Point3) -> Certified {
142    match insphere_filter(a, b, c, d, e) {
143        Certified::Certain { sign, .. } => Certified::exact_sign(sign),
144        _ => Certified::exact_sign(insphere_exact(a, b, c, d, e)),
145    }
146}
147
148/// The fast filter alone, exposed so escalation can be measured.
149#[must_use]
150pub fn insphere_filter(a: Point3, b: Point3, c: Point3, d: Point3, e: Point3) -> Certified {
151    let v = |p: Point3| (p.x - e.x, p.y - e.y, p.z - e.z);
152    let (ax, ay, az) = v(a);
153    let (bx, by, bz) = v(b);
154    let (cx, cy, cz) = v(c);
155    let (dx, dy, dz) = v(d);
156
157    let ab = ax * by - bx * ay;
158    let bc = bx * cy - cx * by;
159    let cd = cx * dy - dx * cy;
160    let da = dx * ay - ax * dy;
161    let ac = ax * cy - cx * ay;
162    let bd = bx * dy - dx * by;
163
164    let abc = az * bc - bz * ac + cz * ab;
165    let bcd = bz * cd - cz * bd + dz * bc;
166    let cda = cz * da + dz * ac + az * cd;
167    let dab = dz * ab + az * bd + bz * da;
168
169    let alift = ax * ax + ay * ay + az * az;
170    let blift = bx * bx + by * by + bz * bz;
171    let clift = cx * cx + cy * cy + cz * cz;
172    let dlift = dx * dx + dy * dy + dz * dz;
173
174    let determinant = (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
175
176    let permanent =
177        (abc.abs() * dlift + dab.abs() * clift) + (cda.abs() * blift + bcd.abs() * alift);
178
179    Certified::from_filter(
180        determinant,
181        INSPHERE_ERROR_FACTOR * permanent,
182        Precision::F64,
183    )
184}
185
186/// Exact sign of the lifted 4x4 `insphere` determinant.
187///
188/// Expands along the lifted column: each 3x3 minor is built from exact 2x2
189/// cofactors, scaled by the remaining z difference, then by the squared
190/// distance. Nothing is rounded between those steps.
191#[must_use]
192fn insphere_exact(a: Point3, b: Point3, c: Point3, d: Point3, e: Point3) -> Sign {
193    // Differences exact, as for `incircle_exact` (#190).
194    let v = |p: Point3| [diff(p.x, e.x), diff(p.y, e.y), diff(p.z, e.z)];
195    let (a3, b3, c3, d3) = (v(a), v(b), v(c), v(d));
196
197    let minor3 = |p: &[Vec<f64>; 3], q: &[Vec<f64>; 3], r: &[Vec<f64>; 3]| {
198        let qr = minor(&q[0], &r[1], &r[0], &q[1]);
199        let rp = minor(&r[0], &p[1], &p[0], &r[1]);
200        let pq = minor(&p[0], &q[1], &q[0], &p[1]);
201        expansion_sum(
202            &expansion_sum(
203                &expansion_product(&qr, &p[2]),
204                &expansion_product(&rp, &q[2]),
205            ),
206            &expansion_product(&pq, &r[2]),
207        )
208    };
209
210    let bcd = minor3(&b3, &c3, &d3);
211    let cda = minor3(&c3, &d3, &a3);
212    let dab = minor3(&d3, &a3, &b3);
213    let abc = minor3(&a3, &b3, &c3);
214
215    // Cofactor expansion signs alternate: +d -c +b -a.
216    let total = expansion_sum(
217        &expansion_sum(&lift3(&abc, &d3), &negate_expansion(&lift3(&dab, &c3))),
218        &expansion_sum(&lift3(&cda, &b3), &negate_expansion(&lift3(&bcd, &a3))),
219    );
220    expansion_sign(&total)
221}
222
223/// Multiply an expansion by `x*x + y*y + z*z`, exactly.
224#[must_use]
225fn lift3(e: &[f64], p: &[Vec<f64>; 3]) -> Vec<f64> {
226    let square = expansion_sum(
227        &expansion_sum(
228            &expansion_product(&p[0], &p[0]),
229            &expansion_product(&p[1], &p[1]),
230        ),
231        &expansion_product(&p[2], &p[2]),
232    );
233    expansion_product(e, &square)
234}
axiolid_predicates/sphere.rs

axiolid_predicates/
sphere.rs
