terra_texture_rs/
hillshade.rs

1//! Hillshade of a DEM, as a fused kernel.
2//!
3//! Mirrors `TerraTexture.derivatives`' `hillshade()` exactly, including
4//! its `az = 360 - azimuth + 90` convention and `np.gradient`'s
5//! `edge_order=1` boundary handling. Verified against `derivatives.py`
6//! in `tests/test_derivatives_rust.py`.
7//!
8//! # NaN handling
9//!
10//! Same as [`curvature`](crate::curvature): the DEM must be NaN-free,
11//! and the Python caller nan-fills before and re-masks after.
12//!
13//! # Formula
14//!
15//! The textbook form is
16//!
17//! ```text
18//! slope  = π/2 - atan(hypot(zx, zy))
19//! aspect = atan2(-zx, zy)
20//! shaded = sin(alt)·sin(slope) + cos(alt)·cos(slope)·cos(az - aspect)
21//! ```
22//!
23//! Expanding it with `sin(atan g) = g/√(1+g²)`, `cos(atan g) = 1/√(1+g²)`
24//! and the angle-difference identity, the `g = hypot(zx, zy)` factor
25//! cancels, leaving
26//!
27//! ```text
28//! shaded = (sin(alt) + cos(alt)·(cos(az)·zy - sin(az)·zx)) / √(1 + zx² + zy²)
29//! ```
30//!
31//! That is one square root per pixel and no inverse trig. Flat cells
32//! (zx = zy = 0) need no special case, since the denominator is 1 rather
33//! than the 0/0 that `hypot(0, 0)` would produce in the aspect term.
34
35use ndarray::{Array2, ArrayView2, Zip};
36
37use crate::common::{gradient2d, PARALLEL_THRESHOLD};
38
39/// Hillshade for one cell.
40///
41/// Takes the cell's gradients `zx`, `zy` and the precomputed sines and
42/// cosines of the (converted) azimuth and altitude, all `f32`. Returns
43/// the illumination as an `f32` clamped to \[0, 1\]. See the module docs
44/// for the formula.
45#[inline]
46fn hillshade_pixel(zx: f32, zy: f32, sin_az: f32, cos_az: f32, sin_alt: f32, cos_alt: f32) -> f32 {
47    let denom = (1.0 + zx * zx + zy * zy).sqrt();
48    let shaded = (sin_alt + cos_alt * (cos_az * zy - sin_az * zx)) / denom;
49    shaded.clamp(0.0, 1.0)
50}
51
52/// Compute a hillshade (simulated illumination) of a DEM.
53///
54/// # Arguments
55///
56/// * `dem` - `ArrayView2<f32>`, shape (H, W): elevations. **Must be
57///   NaN-free** (see module docs). Any memory layout. H and W should
58///   each be at least 2; an axis of length 1 is treated as flat, where
59///   numpy would raise.
60/// * `cellsize` - `f32`: grid spacing, in the same units as the
61///   elevations.
62/// * `azimuth` - `f32`, degrees: direction the light comes from,
63///   clockwise from north (e.g. `315.0` for north-west). Converted
64///   internally with `derivatives.py`'s `az = 360 - azimuth + 90`.
65/// * `altitude` - `f32`, degrees: height of the light above the horizon
66///   (`0.0` = horizon, `90.0` = directly overhead).
67/// * `out` - `&mut Array2<f32>`, shape (H, W): overwritten with the
68///   hillshade, every element in \[0, 1\] (0 = fully shaded, 1 = fully lit).
69///
70/// The four trig values are computed once per call, not per pixel.
71/// Runs serially below [`PARALLEL_THRESHOLD`]
72/// elements and in parallel at or above it.
73///
74/// # Panics
75///
76/// * If `out` does not have the same shape as `dem`.
77/// * If exactly one of H and W is 0.
78///
79/// # Example
80///
81/// ```
82/// use ndarray::Array2;
83/// use terra_texture_rs::hillshade_core;
84///
85/// // flat ground lit from directly overhead is fully lit
86/// let dem = Array2::<f32>::zeros((4, 4));
87/// let mut out = Array2::<f32>::zeros(dem.raw_dim());
88/// hillshade_core(dem.view(), 1.0, 315.0, 90.0, &mut out);
89///
90/// assert!(out.iter().all(|&v| (v - 1.0).abs() < 1e-6));
91/// ```
92pub fn hillshade_core(dem: ArrayView2<f32>, cellsize: f32, azimuth: f32, altitude: f32, out: &mut Array2<f32>) {
93    let (zy, zx) = gradient2d(dem, cellsize);
94    let az = (360.0 - azimuth + 90.0).to_radians();
95    let alt = altitude.to_radians();
96    // Computed ONCE for the whole DEM, not per pixel: 4 trig calls total
97    // instead of up to 2*H*W.
98    let (sin_az, cos_az) = az.sin_cos();
99    let (sin_alt, cos_alt) = alt.sin_cos();
100
101    let n = dem.len();
102    let combine = |o: &mut f32, &zx: &f32, &zy: &f32| {
103        *o = hillshade_pixel(zx, zy, sin_az, cos_az, sin_alt, cos_alt);
104    };
105    let z = Zip::from(out).and(&zx).and(&zy);
106    if n >= PARALLEL_THRESHOLD {
107        z.par_for_each(combine);
108    } else {
109        z.for_each(combine);
110    }
111}
terra_texture_rs/hillshade.rs

terra_texture_rs/
hillshade.rs
