axiolid_overlay/
circle.rs

1//! The minimum enclosing circle of a point set (#118).
2//!
3//! # Exact choice, enclosed output
4//!
5//! Welzl's algorithm, in its iterative form: points are visited in a fixed
6//! pseudo-random order, and a point outside the current circle is put on
7//! the boundary of the next one. The one geometric decision is whether a
8//! point lies in the circle spanned by one, two or three support points,
9//! and it is exact:
10//!
11//! - one support point: equality;
12//! - two, the circle on them as diameter: the sign of `(p - a) . (p - b)`;
13//! - three, their circumcircle: the incircle determinant against the
14//!   triangle's orientation.
15//!
16//! Each is a polynomial in differences of the input `f64`s, decided in
17//! intervals and else in dyadics, so the support set is the exact minimum
18//! circle's. The visiting order is fixed, so the result does not depend on
19//! chance; the circle does not depend on the input order either (the
20//! minimum circle is unique), though which support points are reported can
21//! when four or more lie on it.
22//!
23//! Only the output is rounded. The exact centre is enclosed from the
24//! support points (a midpoint, or a circumcentre as a quotient of exact
25//! dyadic polynomials), and the returned radius is rounded up so the
26//! returned circle contains the exact one. [`CircleEvidence::error`] bounds
27//! both the centre's distance from the exact centre and the radius's excess
28//! over the exact radius.
29
30use axiolid_core::Point2;
31use axiolid_exact::{certify, Arith, Dyadic, Interval, SignExpr};
32use axiolid_guarantees::Sign;
33
34/// A circle by its centre and radius.
35#[derive(Debug, Clone, Copy, PartialEq)]
36pub struct EnclosingCircle {
37    /// Centre.
38    pub centre: Point2,
39    /// Radius; zero for a single distinct point.
40    pub radius: f64,
41}
42
43/// Which points determine the circle and how exact the output is.
44#[derive(Debug, Clone, PartialEq)]
45#[non_exhaustive]
46pub struct CircleEvidence {
47    /// Indices into the input of the one, two or three points the exact
48    /// minimum circle passes through and is determined by, ascending.
49    pub support: Vec<usize>,
50    /// A bound on the distance between the returned centre and the exact
51    /// one, and on the returned radius's excess over the exact radius. The
52    /// returned radius is never below the exact radius, so the returned
53    /// circle contains every input point. Zero for a single point.
54    pub error: f64,
55}
56
57/// The circle and its evidence.
58#[derive(Debug, Clone, PartialEq)]
59#[non_exhaustive]
60pub struct MinimumCircle {
61    /// The circle.
62    pub circle: EnclosingCircle,
63    /// Its support and error bound.
64    pub evidence: CircleEvidence,
65}
66
67/// Why no circle was built.
68#[derive(Debug, Clone, Copy, PartialEq, Eq)]
69#[non_exhaustive]
70pub enum CircleError {
71    /// No points.
72    Empty,
73    /// A coordinate was not finite.
74    NonFinite,
75}
76
77/// `(p - a) . (p - b)`: not positive exactly when `p` lies in the circle
78/// on `a b` as diameter.
79struct Diametral {
80    a: Point2,
81    b: Point2,
82    p: Point2,
83}
84
85impl SignExpr for Diametral {
86    fn sign_in<T: Arith>(&self) -> Option<Sign> {
87        let f = T::from_f64;
88        let (px, py) = (f(self.p.x), f(self.p.y));
89        let ax = px.sub(&f(self.a.x));
90        let ay = py.sub(&f(self.a.y));
91        let bx = px.sub(&f(self.b.x));
92        let by = py.sub(&f(self.b.y));
93        ax.mul(&bx).add(&ay.mul(&by)).sign()
94    }
95}
96
97/// Orientation of `c` against `a -> b`.
98struct Orient {
99    a: Point2,
100    b: Point2,
101    c: Point2,
102}
103
104impl SignExpr for Orient {
105    fn sign_in<T: Arith>(&self) -> Option<Sign> {
106        let f = T::from_f64;
107        let (ux, uy) = (f(self.b.x).sub(&f(self.a.x)), f(self.b.y).sub(&f(self.a.y)));
108        let (vx, vy) = (f(self.c.x).sub(&f(self.a.x)), f(self.c.y).sub(&f(self.a.y)));
109        ux.mul(&vy).sub(&uy.mul(&vx)).sign()
110    }
111}
112
113/// The incircle determinant: positive when `p` lies inside the circle
114/// through `a b c` taken counter-clockwise.
115struct InCircle {
116    a: Point2,
117    b: Point2,
118    c: Point2,
119    p: Point2,
120}
121
122impl SignExpr for InCircle {
123    fn sign_in<T: Arith>(&self) -> Option<Sign> {
124        let f = T::from_f64;
125        let row = |q: Point2| {
126            let x = f(q.x).sub(&f(self.p.x));
127            let y = f(q.y).sub(&f(self.p.y));
128            let l = x.mul(&x).add(&y.mul(&y));
129            (x, y, l)
130        };
131        let (ax, ay, al) = row(self.a);
132        let (bx, by, bl) = row(self.b);
133        let (cx, cy, cl) = row(self.c);
134        let minor = |x1: &T, y1: &T, x2: &T, y2: &T| x1.mul(y2).sub(&y1.mul(x2));
135        al.mul(&minor(&bx, &by, &cx, &cy))
136            .sub(&bl.mul(&minor(&ax, &ay, &cx, &cy)))
137            .add(&cl.mul(&minor(&ax, &ay, &bx, &by)))
138            .sign()
139    }
140}
141
142/// The inputs are finite, so the exact tier always decides.
143fn sign<E: SignExpr>(e: &E) -> Sign {
144    certify(e).unwrap_or(Sign::Zero)
145}
146
147/// Whether `p` lies in (or on) the circle determined by `support`.
148fn inside(points: &[Point2], support: &[usize], p: Point2) -> bool {
149    match *support {
150        [a] => points[a] == p,
151        [a, b] => {
152            sign(&Diametral {
153                a: points[a],
154                b: points[b],
155                p,
156            }) != Sign::Positive
157        }
158        [a, b, c] => {
159            let (a, b, c) = (points[a], points[b], points[c]);
160            let turn = sign(&Orient { a, b, c });
161            let side = sign(&InCircle { a, b, c, p });
162            side == Sign::Zero || side == turn
163        }
164        _ => unreachable!("a circle has one to three support points"),
165    }
166}
167
168/// A fixed pseudo-random permutation of `0..n` (Fisher-Yates driven by
169/// splitmix64 from a constant seed): the expected linear running time of
170/// a random order, with a result that never changes between runs.
171fn visiting_order(n: usize) -> Vec<usize> {
172    let mut order: Vec<usize> = (0..n).collect();
173    let mut state: u64 = 0x9e37_79b9_7f4a_7c15;
174    let mut next = || {
175        state = state.wrapping_add(0x9e37_79b9_7f4a_7c15);
176        let mut z = state;
177        z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9);
178        z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb);
179        z ^ (z >> 31)
180    };
181    for i in (1..n).rev() {
182        let j = (next() % (i as u64 + 1)) as usize;
183        order.swap(i, j);
184    }
185    order
186}
187
188/// The support of the minimum circle, by Welzl's iterative algorithm.
189fn welzl(points: &[Point2]) -> Vec<usize> {
190    let order = visiting_order(points.len());
191    let mut support = vec![order[0]];
192    for i in 1..order.len() {
193        let pi = order[i];
194        if inside(points, &support, points[pi]) {
195            continue;
196        }
197        // `pi` lies on the minimum circle of the points visited so far.
198        support = vec![pi];
199        for j in 0..i {
200            let pj = order[j];
201            if inside(points, &support, points[pj]) {
202                continue;
203            }
204            // So do `pi` and `pj`.
205            support = vec![pi, pj];
206            for &pk in &order[..j] {
207                if !inside(points, &support, points[pk]) {
208                    support = vec![pi, pj, pk];
209                }
210            }
211        }
212    }
213    support
214}
215
216fn exact(x: f64) -> Dyadic {
217    Dyadic::from_f64(x)
218}
219
220/// Enclosures of the exact centre of the circle on `support`.
221fn centre_enclosure(points: &[Point2], support: &[usize]) -> [Interval; 2] {
222    match *support {
223        [a] => [Interval::point(points[a].x), Interval::point(points[a].y)],
224        [a, b] => {
225            let half = exact(0.5);
226            let mid = |s: f64, t: f64| exact(s).add(&exact(t)).mul(&half).enclosure();
227            [mid(points[a].x, points[b].x), mid(points[a].y, points[b].y)]
228        }
229        [a, b, c] => {
230            let (a, b, c) = (points[a], points[b], points[c]);
231            let (ux, uy) = (exact(b.x).sub(&exact(a.x)), exact(b.y).sub(&exact(a.y)));
232            let (vx, vy) = (exact(c.x).sub(&exact(a.x)), exact(c.y).sub(&exact(a.y)));
233            let uu = ux.mul(&ux).add(&uy.mul(&uy));
234            let vv = vx.mul(&vx).add(&vy.mul(&vy));
235            let d = ux.mul(&vy).sub(&uy.mul(&vx)).mul(&exact(2.0)).enclosure();
236            let nx = uu.mul(&vy).sub(&vv.mul(&uy)).enclosure();
237            let ny = vv.mul(&ux).sub(&uu.mul(&vx)).enclosure();
238            [
239                Interval::point(a.x).add(&nx.quotient(d)),
240                Interval::point(a.y).add(&ny.quotient(d)),
241            ]
242        }
243        _ => unreachable!("a circle has one to three support points"),
244    }
245}
246
247/// An upper bound on the distance from `p` to any point of the box
248/// `centre`, and a lower bound on the distance to the nearest.
249fn distance_bounds(p: Point2, centre: &[Interval; 2]) -> (f64, f64) {
250    let dx = Interval::point(p.x).sub(&centre[0]);
251    let dy = Interval::point(p.y).sub(&centre[1]);
252    let squared = dx.mul(&dx).add(&dy.mul(&dy));
253    // `sqrt` is correctly rounded, so one step outward covers it.
254    let low = squared.lo().max(0.0).sqrt().next_down().max(0.0);
255    (low, squared.hi().sqrt().next_up())
256}
257
258/// The minimum circle enclosing `points`.
259///
260/// # Errors
261///
262/// [`CircleError::Empty`] for no points, [`CircleError::NonFinite`] for a
263/// coordinate that is not finite.
264pub fn minimum_enclosing_circle(points: &[Point2]) -> Result<MinimumCircle, CircleError> {
265    if points.is_empty() {
266        return Err(CircleError::Empty);
267    }
268    if !points.iter().all(|p| p.is_finite()) {
269        return Err(CircleError::NonFinite);
270    }
271    let mut support = welzl(points);
272    support.sort_unstable();
273    if let [a] = *support {
274        return Ok(MinimumCircle {
275            circle: EnclosingCircle {
276                centre: points[a],
277                radius: 0.0,
278            },
279            evidence: CircleEvidence {
280                support,
281                error: 0.0,
282            },
283        });
284    }
285    let enclosure = centre_enclosure(points, &support);
286    let mid = |i: Interval| i.lo() + 0.5 * (i.hi() - i.lo());
287    let centre = Point2::new(mid(enclosure[0]), mid(enclosure[1]));
288    // How far the returned centre can be from the exact one: the sum of
289    // the per-axis gaps is at least their Euclidean length.
290    let gap = |i: Interval, m: f64| (i.hi() - m).max(m - i.lo());
291    let centre_error = (gap(enclosure[0], centre.x) + gap(enclosure[1], centre.y)).next_up();
292    let (radius_low, radius_high) = distance_bounds(points[support[0]], &enclosure);
293    let mut radius = (radius_high + centre_error).next_up();
294    // The circle of that radius about the returned centre contains the
295    // exact circle; confirm every point with enclosures all the same.
296    let at = [Interval::point(centre.x), Interval::point(centre.y)];
297    for &p in points {
298        radius = radius.max(distance_bounds(p, &at).1);
299    }
300    Ok(MinimumCircle {
301        circle: EnclosingCircle { centre, radius },
302        evidence: CircleEvidence {
303            support,
304            error: (radius - radius_low).next_up(),
305        },
306    })
307}
axiolid_overlay/circle.rs

axiolid_overlay/
circle.rs
