terra_texture_rs/
common.rs

1//! Pieces shared across kernel modules: the serial/parallel dispatch
2//! threshold, and the numpy-compatible 2-D gradient used by both
3//! `curvature` and `hillshade`.
4
5use ndarray::{Array2, ArrayView2};
6
7/// Element count at or above which kernels switch from a serial loop to
8/// rayon's parallel iteration.
9///
10/// Type: `usize`, counted in array elements. For
11/// [`luminosity_blend_core`](crate::luminosity_blend_core) it is counted
12/// in pixels (H × W) rather than elements (H × W × 3).
13///
14/// **Provisional placeholder, not a measured value.** 256 × 256 = 65,536
15/// was picked only as a plausible order of magnitude for rayon's
16/// dispatch overhead to be amortized. Run `cargo bench` and replace it
17/// with the crossover point `benches/blend_bench.rs` actually measures
18/// on the target hardware.
19// TODO(benchmark): replace with a measured value.
20pub const PARALLEL_THRESHOLD: usize = 65_536;
21
22/// numpy-compatible `np.gradient(arr, spacing)` along both axes.
23///
24/// Uses numpy's default `edge_order=1`: a one-sided forward/backward
25/// difference at the first/last index of each axis, and a second-order
26/// central difference everywhere in between.
27///
28/// # Arguments
29///
30/// * `arr` - `ArrayView2<f32>`, shape (H, W). Any memory layout.
31/// * `spacing` - `f32`, grid spacing (the DEM cell size), applied to
32///   both axes.
33///
34/// # Returns
35///
36/// `(d_axis0, d_axis1)`: two newly allocated, C-contiguous
37/// `Array2<f32>` of shape (H, W), the derivative down the rows and across
38/// the columns respectively. Same order as numpy's
39/// `zy, zx = np.gradient(dem, cellsize)`.
40///
41/// # Differences from numpy
42///
43/// If an axis has length 1, this returns zeros for that axis. numpy
44/// instead raises `ValueError` (it needs at least `edge_order + 1 = 2`
45/// elements per axis).
46///
47/// # Panics
48///
49/// If either axis has length 0 while the other does not (out-of-bounds
50/// index on the empty axis).
51pub(crate) fn gradient2d(arr: ArrayView2<f32>, spacing: f32) -> (Array2<f32>, Array2<f32>) {
52    let (h, w) = arr.dim();
53    let mut d_axis0 = Array2::<f32>::zeros((h, w));
54    let mut d_axis1 = Array2::<f32>::zeros((h, w));
55
56    // axis 0 (down rows), column by column
57    if h == 1 {
58        // np.gradient on a length-1 axis returns zeros
59        d_axis0.fill(0.0);
60    } else {
61        for j in 0..w {
62            d_axis0[[0, j]] = (arr[[1, j]] - arr[[0, j]]) / spacing;
63            d_axis0[[h - 1, j]] = (arr[[h - 1, j]] - arr[[h - 2, j]]) / spacing;
64        }
65        for i in 1..h - 1 {
66            for j in 0..w {
67                d_axis0[[i, j]] = (arr[[i + 1, j]] - arr[[i - 1, j]]) / (2.0 * spacing);
68            }
69        }
70    }
71
72    // axis 1 (across columns), row by row
73    if w == 1 {
74        d_axis1.fill(0.0);
75    } else {
76        for i in 0..h {
77            d_axis1[[i, 0]] = (arr[[i, 1]] - arr[[i, 0]]) / spacing;
78            d_axis1[[i, w - 1]] = (arr[[i, w - 1]] - arr[[i, w - 2]]) / spacing;
79            for j in 1..w - 1 {
80                d_axis1[[i, j]] = (arr[[i, j + 1]] - arr[[i, j - 1]]) / (2.0 * spacing);
81            }
82        }
83    }
84
85    (d_axis0, d_axis1)
86}
terra_texture_rs/common.rs

terra_texture_rs/
common.rs
