lindhard/
contour.rs

1//! Iso-dose contours of a gridded scalar field, with the area and centroid of
2//! each closed contour.
3//!
4//! The contours are found by marching squares (W. E. Lorensen and H. E. Cline,
5//! "Marching cubes: a high resolution 3D surface construction algorithm",
6//! Comput. Graph. 21, 163 (1987), the 2D case), with linear interpolation
7//! along the cell edges. A grid point is *inside* when its value is
8//! `>= threshold`.
9//!
10//! * Every contour is directed with the inside on its left. An outer
11//!   boundary of an exposed feature is therefore counter-clockwise (positive
12//!   [`Contour::signed_area_m2`]) and the boundary of a hole in a feature is
13//!   clockwise (negative). The area of a feature with holes is the sum of the
14//!   signed areas of its outer contour and its holes.
15//! * The saddle cells (two diagonal inside corners) are resolved by the mean
16//!   of the four corner values: inside connects the inside corners when the
17//!   mean is `>= threshold`, and separates them otherwise.
18//! * **Open contours.** A feature that touches the grid edge has a contour
19//!   that leaves the grid. It is returned with `closed = false`, in the
20//!   order walked, and has no area and no centroid (they are `None`): the
21//!   grid does not say where the feature ends. Enlarge the grid to close it.
22//! * The area of a closed contour is the shoelace area of its polygon, and the
23//!   centroid is the area centroid of that polygon (Bourke's polygon
24//!   formulas; textbook). The polygon is the piecewise-linear contour, so the
25//!   area differs from the smooth contour's by `O(h²)` in the grid spacing.
26//!
27//! Contours are produced in a fixed order (by the cell and edge where each
28//! starts), so the output is deterministic.
29
30use std::collections::BTreeMap;
31
32use crate::exposure::{DoseGrid, Grid};
33
34/// One contour.
35#[derive(Debug, Clone, PartialEq)]
36pub struct Contour {
37    /// Vertices `[x, y]`, m. For a closed contour the first vertex is not
38    /// repeated at the end.
39    pub points: Vec<[f64; 2]>,
40    /// True if the contour closes on itself inside the grid.
41    pub closed: bool,
42    /// Signed shoelace area, m² (positive counter-clockwise, i.e. an outer
43    /// boundary); `None` for an open contour.
44    pub signed_area_m2: Option<f64>,
45    /// Area centroid `[x, y]`, m; `None` for an open contour or one of zero
46    /// area.
47    pub centroid_m: Option<[f64; 2]>,
48}
49
50/// A grid edge: `(i, j, vertical)`; horizontal runs from `(i, j)` to
51/// `(i+1, j)`, vertical from `(i, j)` to `(i, j+1)`.
52type Edge = (usize, usize, bool);
53
54/// The contours of `dose` at `threshold`.
55pub fn contours(dose: &DoseGrid, threshold: f64) -> Vec<Contour> {
56    contours_of(&dose.grid, &dose.dose_ev_per_m2, threshold)
57}
58
59/// The contours of `values` (row-major on `grid`) at `threshold`.
60pub fn contours_of(grid: &Grid, values: &[f64], threshold: f64) -> Vec<Contour> {
61    assert_eq!(values.len(), grid.len(), "one value per grid point");
62    let (nx, ny) = (grid.nx, grid.ny);
63    let v = |i: usize, j: usize| values[j * nx + i];
64    // Directed segments, keyed by start edge. Marching-square contours never
65    // branch, so each edge starts at most one segment.
66    let mut next: BTreeMap<Edge, Edge> = BTreeMap::new();
67    for j in 0..ny.saturating_sub(1) {
68        for i in 0..nx.saturating_sub(1) {
69            let (bl, br, tr, tl) = (v(i, j), v(i + 1, j), v(i + 1, j + 1), v(i, j + 1));
70            let case = usize::from(bl >= threshold)
71                | usize::from(br >= threshold) << 1
72                | usize::from(tr >= threshold) << 2
73                | usize::from(tl >= threshold) << 3;
74            let b: Edge = (i, j, false);
75            let r: Edge = (i + 1, j, true);
76            let t: Edge = (i, j + 1, false);
77            let l: Edge = (i, j, true);
78            let centre_in = (bl + br + tr + tl) / 4.0 >= threshold;
79            let segs: &[(Edge, Edge)] = match case {
80                1 => &[(b, l)],
81                2 => &[(r, b)],
82                3 => &[(r, l)],
83                4 => &[(t, r)],
84                5 if centre_in => &[(b, r), (t, l)],
85                5 => &[(b, l), (t, r)],
86                6 => &[(t, b)],
87                7 => &[(t, l)],
88                8 => &[(l, t)],
89                9 => &[(b, t)],
90                10 if centre_in => &[(l, b), (r, t)],
91                10 => &[(r, b), (l, t)],
92                11 => &[(r, t)],
93                12 => &[(l, r)],
94                13 => &[(b, r)],
95                14 => &[(l, b)],
96                _ => &[],
97            };
98            for &(from, to) in segs {
99                next.insert(from, to);
100            }
101        }
102    }
103    let point = |e: Edge| -> [f64; 2] {
104        let (i, j, vertical) = e;
105        let (i2, j2) = if vertical { (i, j + 1) } else { (i + 1, j) };
106        let (a, b) = (v(i, j), v(i2, j2));
107        let t = (threshold - a) / (b - a);
108        let [xa, ya] = grid.point(i, j);
109        let [xb, yb] = grid.point(i2, j2);
110        [xa + t * (xb - xa), ya + t * (yb - ya)]
111    };
112    let has_incoming: std::collections::BTreeSet<Edge> = next.values().copied().collect();
113    let mut visited: std::collections::BTreeSet<Edge> = Default::default();
114    let mut out = Vec::new();
115    // Open contours start at an edge nothing leads into.
116    let starts: Vec<Edge> = next
117        .keys()
118        .copied()
119        .filter(|e| !has_incoming.contains(e))
120        .collect();
121    for s in starts {
122        let mut pts = vec![point(s)];
123        visited.insert(s);
124        let mut cur = s;
125        while let Some(&n) = next.get(&cur) {
126            pts.push(point(n));
127            visited.insert(n);
128            cur = n;
129        }
130        out.push(open(pts));
131    }
132    // What remains is closed loops.
133    for &s in next.keys() {
134        if visited.contains(&s) {
135            continue;
136        }
137        let mut pts = Vec::new();
138        let mut cur = s;
139        loop {
140            visited.insert(cur);
141            pts.push(point(cur));
142            cur = next[&cur];
143            if cur == s {
144                break;
145            }
146        }
147        out.push(closed(pts));
148    }
149    out
150}
151
152fn open(points: Vec<[f64; 2]>) -> Contour {
153    Contour {
154        points,
155        closed: false,
156        signed_area_m2: None,
157        centroid_m: None,
158    }
159}
160
161fn closed(points: Vec<[f64; 2]>) -> Contour {
162    let n = points.len();
163    let (mut a2, mut cx, mut cy) = (0.0, 0.0, 0.0);
164    for k in 0..n {
165        let [x0, y0] = points[k];
166        let [x1, y1] = points[(k + 1) % n];
167        let cross = x0 * y1 - x1 * y0;
168        a2 += cross;
169        cx += (x0 + x1) * cross;
170        cy += (y0 + y1) * cross;
171    }
172    let area = a2 / 2.0;
173    let centroid = (area != 0.0).then(|| [cx / (6.0 * area), cy / (6.0 * area)]);
174    Contour {
175        points,
176        closed: true,
177        signed_area_m2: Some(area),
178        centroid_m: centroid,
179    }
180}
181
182impl Contour {
183    /// CSV of a list of contours: `contour,closed,signed_area_m2,centroid_x_m,
184    /// centroid_y_m,x_m,y_m`, one row per vertex (area and centroid on every
185    /// row of the contour, empty for an open one).
186    pub fn csv(contours: &[Contour]) -> String {
187        let mut s =
188            String::from("contour,closed,signed_area_m2,centroid_x_m,centroid_y_m,x_m,y_m\n");
189        let f = |o: Option<f64>| o.map(|v| v.to_string()).unwrap_or_default();
190        for (k, c) in contours.iter().enumerate() {
191            for p in &c.points {
192                s.push_str(&format!(
193                    "{k},{},{},{},{},{},{}\n",
194                    u8::from(c.closed),
195                    f(c.signed_area_m2),
196                    f(c.centroid_m.map(|m| m[0])),
197                    f(c.centroid_m.map(|m| m[1])),
198                    p[0],
199                    p[1]
200                ));
201            }
202        }
203        s
204    }
205}
206
207#[cfg(test)]
208mod tests {
209    use super::*;
210
211    fn grid(n: usize) -> Grid {
212        Grid {
213            x0_m: -1.0,
214            y0_m: -1.0,
215            dx_m: 2.0 / (n - 1) as f64,
216            dy_m: 2.0 / (n - 1) as f64,
217            nx: n,
218            ny: n,
219        }
220    }
221
222    fn field(g: &Grid, f: impl Fn(f64, f64) -> f64) -> Vec<f64> {
223        (0..g.len())
224            .map(|k| {
225                let [x, y] = g.point(k % g.nx, k / g.nx);
226                f(x, y)
227            })
228            .collect()
229    }
230
231    #[test]
232    fn circle_area_and_centroid() {
233        let g = grid(201);
234        let v = field(&g, |x, y| {
235            0.3 - ((x - 0.1).powi(2) + (y + 0.2).powi(2)).sqrt()
236        });
237        let c = contours_of(&g, &v, 0.0);
238        assert_eq!(c.len(), 1);
239        let a = c[0].signed_area_m2.unwrap();
240        assert!(c[0].closed && a > 0.0);
241        assert!((a - std::f64::consts::PI * 0.09).abs() / (std::f64::consts::PI * 0.09) < 1e-3);
242        let m = c[0].centroid_m.unwrap();
243        assert!((m[0] - 0.1).abs() < 1e-3 && (m[1] + 0.2).abs() < 1e-3);
244    }
245
246    #[test]
247    fn two_features_and_a_hole() {
248        let g = grid(161);
249        // Two discs, one with a hole (an annulus), via a radial field.
250        let v = field(&g, |x, y| {
251            let d1 = ((x + 0.5).powi(2) + y * y).sqrt();
252            let d2 = ((x - 0.5).powi(2) + y * y).sqrt();
253            (0.3 - d1).max(0.1 - (d2 - 0.2).abs())
254        });
255        let c = contours_of(&g, &v, 0.0);
256        let mut areas: Vec<f64> = c.iter().map(|c| c.signed_area_m2.unwrap()).collect();
257        areas.sort_by(f64::total_cmp);
258        // Annulus: outer radius 0.3 (area π 0.09), inner 0.1 (hole, negative);
259        // plus the 0.3 disc.
260        assert_eq!(areas.len(), 3, "{areas:?}");
261        assert!(areas[0] < 0.0);
262        let pi = std::f64::consts::PI;
263        assert!((areas[0] + pi * 0.01).abs() < 2e-3);
264        assert!((areas[1] - pi * 0.09).abs() < 2e-3 && (areas[2] - pi * 0.09).abs() < 2e-3);
265    }
266
267    #[test]
268    fn no_contour_when_threshold_above_or_below_all() {
269        let g = grid(11);
270        let v = field(&g, |x, y| x + y);
271        assert!(contours_of(&g, &v, 10.0).is_empty());
272        assert!(contours_of(&g, &v, -10.0).is_empty());
273    }
274
275    #[test]
276    fn feature_touching_the_edge_is_open_without_area() {
277        let g = grid(41);
278        let v = field(&g, |x, _| 0.2 - x); // inside is x <= 0.2, reaching the edge
279        let c = contours_of(&g, &v, 0.0);
280        assert_eq!(c.len(), 1);
281        assert!(!c[0].closed && c[0].signed_area_m2.is_none() && c[0].centroid_m.is_none());
282    }
283
284    #[test]
285    fn saddle_is_resolved_by_the_cell_mean() {
286        let g = Grid {
287            x0_m: 0.0,
288            y0_m: 0.0,
289            dx_m: 1.0,
290            dy_m: 1.0,
291            nx: 2,
292            ny: 2,
293        };
294        // BL and TR inside; mean 0.5 >= 0.5 joins them (one open contour per
295        // outside corner), mean below separates them (one per inside corner).
296        let joined = contours_of(&g, &[1.0, 0.0, 0.0, 1.0], 0.4);
297        let split = contours_of(&g, &[1.0, 0.0, 0.0, 1.0], 0.6);
298        assert_eq!((joined.len(), split.len()), (2, 2));
299        assert_ne!(joined, split);
300    }
301}
lindhard/contour.rs

lindhard/
contour.rs
