From 650babf9c5c0ab094aa61319a97b4c6cfaafaaaa Mon Sep 17 00:00:00 2001 From: Dominik Strebel Date: Tue, 3 Nov 2026 08:29:16 +0100 Subject: [PATCH] perf(analysis): slope, aspect and hillshade rows on several cores (#1252) The terrain row pass (#1130) runs its rows on parallel::for_each_row: each output row is made from the three input rows around it and written where it lies, so the output is the same bit for bit as on one thread. Hillshade of a 9667 x 11947 DEM (the one #1111 measured): 12.6 s on one thread, 1.4 s on 16. --- crates/analysis/src/rasterops.rs | 81 +++++++++++++++++++++++--------- 1 file changed, 58 insertions(+), 23 deletions(-) diff --git a/crates/analysis/src/rasterops.rs b/crates/analysis/src/rasterops.rs index 22d98a3c1..6b558b566 100644 --- a/crates/analysis/src/rasterops.rs +++ b/crates/analysis/src/rasterops.rs @@ -66,7 +66,13 @@ fn per_cell(r: &Raster, f: impl Fn(i64, i64) -> Option) -> Result, } /// Draw the interior from three borrowed rows; edges and no-data keep the scalar path. -fn terrain(r: &Raster, data: &[f32], zf: f64, f: impl Fn(f64, f64) -> f64) -> Result, String> { +fn terrain(r: &Raster, data: &[f32], zf: f64, f: impl Fn(f64, f64) -> f64 + Sync) -> Result, String> { + terrain_on(None, r, data, zf, f) +} + +/// [`terrain`] on `threads` threads (none: as many as pay off). Each output row is made from the +/// three input rows around it and written where it lies, so rows go on several threads (#1252). +fn terrain_on(threads: Option, r: &Raster, data: &[f32], zf: f64, f: impl Fn(f64, f64) -> f64 + Sync) -> Result, String> { let (w, h) = (r.width as usize, r.height as usize); let cell = r.cell_size_m(); let scalar = |c, row| { @@ -79,28 +85,48 @@ fn terrain(r: &Raster, data: &[f32], zf: f64, f: impl Fn(f64, f64) -> f64) -> Re return per_cell(r, scalar); } let value = |v: Option| v.filter(|v| v.is_finite()).map(|v| v as f32).unwrap_or(NODATA); - let mut out = grid(w, h)?; - out.extend((0..w).map(|c| value(scalar(c as i64, 0)))); - for (row, ((top, mid), bottom)) in data.chunks_exact(w).zip(data.chunks_exact(w).skip(1)).zip(data.chunks_exact(w).skip(2)).enumerate() - { - let row = row + 1; - out.push(value(scalar(0, row as i64))); - for (col, ((a, b), c)) in top.windows(3).zip(mid.windows(3)).zip(bottom.windows(3)).enumerate() { - let ([a, b, c], [d, center, f0], [g, h, i]) = (a, b, c) else { continue }; - if !r.is_valid(*center) { - out.push(NODATA); - } else if [*a, *b, *c, *d, *f0, *g, *h, *i].iter().all(|v| r.is_valid(*v)) { - let [a, b, c, d, f0, g, h, i] = [*a, *b, *c, *d, *f0, *g, *h, *i].map(|v| f64::from(v) * zf); - let x = ((c + 2.0 * f0 + i) - (a + 2.0 * d + g)) / (8.0 * cell.0) * r.transform.pixel_width.signum(); - let y = ((a + 2.0 * b + c) - (g + 2.0 * h + i)) / (8.0 * cell.1) * -r.transform.pixel_height.signum(); - out.push(value(Some(f(x, y)))); - } else { - out.push(value(scalar((col + 1) as i64, row as i64))); + let mut out = new_band(w, h, NODATA)?; + // A cell's gradient and its function cost about one unit of the threshold's (some 30 ns). + let threads = threads.unwrap_or_else(|| crate::parallel::threads_for_rows(h, len, 0)); + crate::control::begin(h); + crate::parallel::for_each_row( + out.chunks_mut(w), + threads, + || (), + |_, row, line: &mut [f32]| { + // The first and last rows, and each row's first and last cells, take the edge path. + let edge = |col: usize| value(scalar(col as i64, row as i64)); + let (Some(top), Some(mid), Some(bottom)) = ( + row.checked_sub(1).and_then(|t| data.get(t * w..(t + 1) * w)), + data.get(row * w..(row + 1) * w), + data.get((row + 1) * w..(row + 2) * w), + ) else { + for (col, o) in line.iter_mut().enumerate() { + *o = edge(col); + } + return; + }; + for (col, o) in line.iter_mut().enumerate() { + *o = if col == 0 || col + 1 == w { + edge(col) + } else { + let at = |s: &[f32], dc: usize| s.get(col + dc - 1).copied().unwrap_or(f32::NAN); + let ([a, b, c], [d, center, f0], [g, hh, i]) = + ([0, 1, 2].map(|k| at(top, k)), [0, 1, 2].map(|k| at(mid, k)), [0, 1, 2].map(|k| at(bottom, k))); + if !r.is_valid(center) { + NODATA + } else if [a, b, c, d, f0, g, hh, i].iter().all(|v| r.is_valid(*v)) { + let [a, b, c, d, f0, g, hh, i] = [a, b, c, d, f0, g, hh, i].map(|v| f64::from(v) * zf); + let x = ((c + 2.0 * f0 + i) - (a + 2.0 * d + g)) / (8.0 * cell.0) * r.transform.pixel_width.signum(); + let y = ((a + 2.0 * b + c) - (g + 2.0 * hh + i)) / (8.0 * cell.1) * -r.transform.pixel_height.signum(); + value(Some(f(x, y))) + } else { + edge(col) + } + }; } - } - out.push(value(scalar((w - 1) as i64, row as i64))); - } - out.extend((0..w).map(|c| value(scalar(c as i64, (h - 1) as i64)))); + }, + )?; Ok(out) } @@ -1199,7 +1225,7 @@ mod tests { #[test] fn terrain_rows_match_scalar_edges_nodata_and_orientations() { - for (w, h) in [(1, 1), (1, 5), (2, 3), (3, 2), (3, 3), (9, 7)] { + for (w, h) in [(1, 1), (1, 5), (2, 3), (3, 2), (3, 3), (9, 7), (60, 41)] { for (pixel_width, pixel_height) in [(0.01, -0.01), (0.01, 0.01), (-0.01, -0.01), (-0.01, 0.01)] { for with_holes in [false, true] { let values = (0..w * h) @@ -1229,6 +1255,15 @@ mod tests { expected, "{w}x{h}, {pixel_width}, {pixel_height}, holes={with_holes}, zf={zf}" ); + // Rows made on several threads give the same, bit for bit (#1252). + let bits = |v: Result, String>| v.map(|v| v.iter().map(|x| x.to_bits()).collect::>()); + for threads in [1, 2, 7, 16] { + assert_eq!( + bits(terrain_on(Some(threads), &r, data, zf, f)), + bits(expected.clone()), + "{w}x{h}, {threads} threads" + ); + } } } }