#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct ContourLine {
pub paths: Vec<String>,
}
pub fn marching_squares(
grid: &[Vec<f64>],
x_coords: &[f64],
y_coords: &[f64],
threshold: f64,
) -> ContourLine {
let rows = grid.len();
if rows < 2 {
return ContourLine { paths: Vec::new() };
}
let cols = grid[0].len();
if cols < 2 || x_coords.len() < cols || y_coords.len() < rows {
return ContourLine { paths: Vec::new() };
}
let mut segments: Vec<((f64, f64), (f64, f64))> = Vec::new();
for row in 0..rows - 1 {
for col in 0..cols - 1 {
let v00 = grid[row][col];
let v10 = grid[row][col + 1];
let v01 = grid[row + 1][col];
let v11 = grid[row + 1][col + 1];
let case = ((v00 >= threshold) as u8) << 3
| ((v10 >= threshold) as u8) << 2
| ((v11 >= threshold) as u8) << 1
| ((v01 >= threshold) as u8);
if case == 0 || case == 15 {
continue; }
let x0 = x_coords[col];
let x1 = x_coords[col + 1];
let y0 = y_coords[row];
let y1 = y_coords[row + 1];
let top = lerp_x(x0, x1, v00, v10, threshold);
let bottom = lerp_x(x0, x1, v01, v11, threshold);
let left = lerp_y(y0, y1, v00, v01, threshold);
let right = lerp_y(y0, y1, v10, v11, threshold);
let top_pt = (top, y0);
let bottom_pt = (bottom, y1);
let left_pt = (x0, left);
let right_pt = (x1, right);
match case {
1 | 14 => segments.push((left_pt, bottom_pt)),
2 | 13 => segments.push((bottom_pt, right_pt)),
3 | 12 => segments.push((left_pt, right_pt)),
4 | 11 => segments.push((top_pt, right_pt)),
5 => {
let avg = (v00 + v10 + v01 + v11) / 4.0;
if avg >= threshold {
segments.push((left_pt, top_pt));
segments.push((bottom_pt, right_pt));
} else {
segments.push((left_pt, bottom_pt));
segments.push((top_pt, right_pt));
}
}
6 | 9 => segments.push((top_pt, bottom_pt)),
7 | 8 => segments.push((left_pt, top_pt)),
10 => {
let avg = (v00 + v10 + v01 + v11) / 4.0;
if avg >= threshold {
segments.push((left_pt, bottom_pt));
segments.push((top_pt, right_pt));
} else {
segments.push((left_pt, top_pt));
segments.push((bottom_pt, right_pt));
}
}
_ => {} }
}
}
let paths = chain_segments_to_svg(segments);
ContourLine { paths }
}
fn lerp_x(x0: f64, x1: f64, v0: f64, v1: f64, threshold: f64) -> f64 {
if (v1 - v0).abs() < 1e-12 {
(x0 + x1) / 2.0
} else {
let t = (threshold - v0) / (v1 - v0);
x0 + t * (x1 - x0)
}
}
fn lerp_y(y0: f64, y1: f64, v0: f64, v1: f64, threshold: f64) -> f64 {
if (v1 - v0).abs() < 1e-12 {
(y0 + y1) / 2.0
} else {
let t = (threshold - v0) / (v1 - v0);
y0 + t * (y1 - y0)
}
}
fn chain_segments_to_svg(segments: Vec<((f64, f64), (f64, f64))>) -> Vec<String> {
if segments.is_empty() {
return Vec::new();
}
let eps = 0.5; let mut used = vec![false; segments.len()];
let mut paths = Vec::new();
for start_idx in 0..segments.len() {
if used[start_idx] {
continue;
}
used[start_idx] = true;
let mut chain: Vec<(f64, f64)> = vec![segments[start_idx].0, segments[start_idx].1];
loop {
let tail = *chain.last().unwrap();
let mut found = false;
for i in 0..segments.len() {
if used[i] {
continue;
}
let (a, b) = segments[i];
if dist(tail, a) < eps {
used[i] = true;
chain.push(b);
found = true;
break;
} else if dist(tail, b) < eps {
used[i] = true;
chain.push(a);
found = true;
break;
}
}
if !found {
break;
}
}
loop {
let head = chain[0];
let mut found = false;
for i in 0..segments.len() {
if used[i] {
continue;
}
let (a, b) = segments[i];
if dist(head, b) < eps {
used[i] = true;
chain.insert(0, a);
found = true;
break;
} else if dist(head, a) < eps {
used[i] = true;
chain.insert(0, b);
found = true;
break;
}
}
if !found {
break;
}
}
if chain.len() >= 2 {
let mut path = format!("M {:.1} {:.1}", chain[0].0, chain[0].1);
for pt in &chain[1..] {
path.push_str(&format!(" L {:.1} {:.1}", pt.0, pt.1));
}
paths.push(path);
}
}
paths
}
fn dist(a: (f64, f64), b: (f64, f64)) -> f64 {
((a.0 - b.0).powi(2) + (a.1 - b.1).powi(2)).sqrt()
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_simple_contour() {
let grid = vec![
vec![0.0, 0.0, 0.0],
vec![0.0, 1.0, 0.0],
vec![0.0, 0.0, 0.0],
];
let x = vec![0.0, 50.0, 100.0];
let y = vec![0.0, 50.0, 100.0];
let contour = marching_squares(&grid, &x, &y, 0.5);
assert!(
!contour.paths.is_empty(),
"Should generate at least one contour path"
);
for path in &contour.paths {
assert!(path.starts_with("M "), "Path should start with M");
}
}
#[test]
fn test_no_contour_all_below() {
let grid = vec![vec![0.0, 0.0], vec![0.0, 0.0]];
let x = vec![0.0, 100.0];
let y = vec![0.0, 100.0];
let contour = marching_squares(&grid, &x, &y, 0.5);
assert!(contour.paths.is_empty(), "No contour when all below");
}
#[test]
fn test_no_contour_all_above() {
let grid = vec![vec![1.0, 1.0], vec![1.0, 1.0]];
let x = vec![0.0, 100.0];
let y = vec![0.0, 100.0];
let contour = marching_squares(&grid, &x, &y, 0.5);
assert!(contour.paths.is_empty(), "No contour when all above");
}
#[test]
fn test_horizontal_contour() {
let grid = vec![vec![1.0, 1.0], vec![0.0, 0.0]];
let x = vec![0.0, 100.0];
let y = vec![0.0, 100.0];
let contour = marching_squares(&grid, &x, &y, 0.5);
assert_eq!(contour.paths.len(), 1);
}
}