ggplot_rs/stat/
marching_squares.rs

1use crate::data::Value;
2
3/// Extract contour line segments from a regular grid at specified level thresholds.
4/// Returns vectors of (x, y, level, group_id) for each segment endpoint.
5#[allow(clippy::too_many_arguments)]
6pub fn extract_contours(
7    grid: &[f64],
8    nx: usize,
9    ny: usize,
10    x_min: f64,
11    y_min: f64,
12    dx: f64,
13    dy: f64,
14    levels: &[f64],
15) -> (Vec<Value>, Vec<Value>, Vec<Value>, Vec<Value>) {
16    let mut x_vals = Vec::new();
17    let mut y_vals = Vec::new();
18    let mut level_vals = Vec::new();
19    let mut group_vals = Vec::new();
20    let mut group_id = 0;
21
22    for &threshold in levels {
23        for iy in 0..ny {
24            for ix in 0..nx {
25                let tl = grid[iy * (nx + 1) + ix];
26                let tr = grid[iy * (nx + 1) + ix + 1];
27                let br = grid[(iy + 1) * (nx + 1) + ix + 1];
28                let bl = grid[(iy + 1) * (nx + 1) + ix];
29
30                let case = ((tl >= threshold) as u8)
31                    | (((tr >= threshold) as u8) << 1)
32                    | (((br >= threshold) as u8) << 2)
33                    | (((bl >= threshold) as u8) << 3);
34
35                if case == 0 || case == 15 {
36                    continue;
37                }
38
39                let cell_x = x_min + ix as f64 * dx;
40                let cell_y = y_min + iy as f64 * dy;
41
42                let segments = marching_squares_segments(
43                    case, tl, tr, br, bl, threshold, cell_x, cell_y, dx, dy,
44                );
45
46                for (p1, p2) in segments {
47                    x_vals.push(Value::Float(p1.0));
48                    y_vals.push(Value::Float(p1.1));
49                    level_vals.push(Value::Float(threshold));
50                    group_vals.push(Value::Integer(group_id));
51
52                    x_vals.push(Value::Float(p2.0));
53                    y_vals.push(Value::Float(p2.1));
54                    level_vals.push(Value::Float(threshold));
55                    group_vals.push(Value::Integer(group_id));
56
57                    group_id += 1;
58                }
59            }
60        }
61    }
62
63    (x_vals, y_vals, level_vals, group_vals)
64}
65
66/// Linear interpolation between two values.
67pub fn lerp(a: f64, b: f64, t: f64) -> f64 {
68    a + t * (b - a)
69}
70
71/// Interpolation fraction for threshold crossing between v0 and v1.
72pub fn frac(v0: f64, v1: f64, threshold: f64) -> f64 {
73    let denom = v1 - v0;
74    if denom.abs() < f64::EPSILON {
75        0.5
76    } else {
77        (threshold - v0) / denom
78    }
79}
80
81/// Return line segments for a marching squares cell.
82/// Corners: tl=top-left, tr=top-right, br=bottom-right, bl=bottom-left.
83#[allow(clippy::too_many_arguments)]
84pub fn marching_squares_segments(
85    case: u8,
86    tl: f64,
87    tr: f64,
88    br: f64,
89    bl: f64,
90    threshold: f64,
91    cx: f64,
92    cy: f64,
93    dx: f64,
94    dy: f64,
95) -> Vec<((f64, f64), (f64, f64))> {
96    let top = || {
97        let t = frac(tl, tr, threshold);
98        (lerp(cx, cx + dx, t), cy)
99    };
100    let bottom = || {
101        let t = frac(bl, br, threshold);
102        (lerp(cx, cx + dx, t), cy + dy)
103    };
104    let left = || {
105        let t = frac(tl, bl, threshold);
106        (cx, lerp(cy, cy + dy, t))
107    };
108    let right = || {
109        let t = frac(tr, br, threshold);
110        (cx + dx, lerp(cy, cy + dy, t))
111    };
112
113    match case {
114        1 => vec![(top(), left())],
115        2 => vec![(top(), right())],
116        3 => vec![(left(), right())],
117        4 => vec![(right(), bottom())],
118        5 => {
119            let avg = (tl + tr + br + bl) / 4.0;
120            if avg >= threshold {
121                vec![(top(), right()), (left(), bottom())]
122            } else {
123                vec![(top(), left()), (right(), bottom())]
124            }
125        }
126        6 => vec![(top(), bottom())],
127        7 => vec![(left(), bottom())],
128        8 => vec![(left(), bottom())],
129        9 => vec![(top(), bottom())],
130        10 => {
131            let avg = (tl + tr + br + bl) / 4.0;
132            if avg >= threshold {
133                vec![(top(), left()), (right(), bottom())]
134            } else {
135                vec![(top(), right()), (left(), bottom())]
136            }
137        }
138        11 => vec![(right(), bottom())],
139        12 => vec![(left(), right())],
140        13 => vec![(top(), right())],
141        14 => vec![(top(), left())],
142        _ => vec![],
143    }
144}
ggplot_rs/stat/marching_squares.rs

ggplot_rs/stat/
marching_squares.rs
