steel_worldgen/density/
spline_eval.rs

1//! Public spline evaluation helpers for generated density function code.
2//!
3//! These operate on raw data (slices, scalars) rather than the `CubicSpline`
4//! struct, so generated code can call them with statically embedded spline data.
5
6/// Binary search to find the interval start: largest `i` where `locations[i] <= input`.
7///
8/// Returns -1 if `input` is before all locations.
9#[inline]
10#[must_use]
11pub fn find_interval(locations: &[f32], input: f32) -> i32 {
12    let mut lo = 0i32;
13    let mut hi = locations.len() as i32;
14    while lo < hi {
15        let mid = i32::midpoint(lo, hi);
16        if input < locations[mid as usize] {
17            hi = mid;
18        } else {
19            lo = mid + 1;
20        }
21    }
22    lo - 1
23}
24
25/// Hermite cubic interpolation between two adjacent spline points.
26///
27/// Given two points `(x1, y1)` and `(x2, y2)` with derivatives `d1` and `d2`,
28/// evaluates the hermite cubic at `input`.
29///
30/// Matches vanilla's formula: `lerp(t, y1, y2) + t * (1 - t) * lerp(t, a, b)`.
31#[inline]
32#[must_use]
33pub fn hermite_interpolate(
34    x1: f32,
35    x2: f32,
36    y1: f32,
37    y2: f32,
38    d1: f32,
39    d2: f32,
40    input: f32,
41) -> f32 {
42    let t = (input - x1) / (x2 - x1);
43    let h = x2 - x1;
44    let a = d1 * h - (y2 - y1);
45    let b = -d2 * h + (y2 - y1);
46    let lerp_y = y1 + t * (y2 - y1);
47    let lerp_ab = a + t * (b - a);
48    lerp_y + t * (1.0 - t) * lerp_ab
49}
50
51/// Evaluate a spline defined by static data arrays.
52///
53/// `locations`, `derivatives` must have the same length.
54/// `value_at(index)` returns the value at a given point index
55/// (can be a constant or a nested spline evaluation).
56///
57/// Uses binary search + hermite cubic interpolation, matching vanilla's
58/// `CubicSpline.Multipoint.apply()`.
59#[inline]
60pub fn evaluate_spline(
61    locations: &[f32],
62    derivatives: &[f32],
63    input: f32,
64    value_at: impl Fn(usize) -> f32,
65) -> f32 {
66    if locations.is_empty() {
67        return 0.0;
68    }
69
70    let last = locations.len() - 1;
71    let start = find_interval(locations, input);
72
73    if start < 0 {
74        let value = value_at(0);
75        return value + derivatives[0] * (input - locations[0]);
76    }
77
78    let start = start as usize;
79    if start == last {
80        let value = value_at(last);
81        return value + derivatives[last] * (input - locations[last]);
82    }
83
84    let y1 = value_at(start);
85    let y2 = value_at(start + 1);
86    hermite_interpolate(
87        locations[start],
88        locations[start + 1],
89        y1,
90        y2,
91        derivatives[start],
92        derivatives[start + 1],
93        input,
94    )
95}
96
97#[cfg(test)]
98mod tests {
99    use super::*;
100
101    #[test]
102    fn find_interval_before_all() {
103        assert_eq!(find_interval(&[0.0, 1.0, 2.0], -1.0), -1);
104    }
105
106    #[test]
107    fn find_interval_exact_match() {
108        assert_eq!(find_interval(&[0.0, 1.0, 2.0], 1.0), 1);
109    }
110
111    #[test]
112    fn find_interval_between() {
113        assert_eq!(find_interval(&[0.0, 1.0, 2.0], 0.5), 0);
114    }
115
116    #[test]
117    fn find_interval_after_all() {
118        assert_eq!(find_interval(&[0.0, 1.0, 2.0], 3.0), 2);
119    }
120
121    #[test]
122    fn hermite_linear() {
123        // With zero derivatives, hermite reduces to linear interpolation
124        let result = hermite_interpolate(0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.5);
125        assert!((result - 0.5).abs() < 1e-6);
126    }
127
128    #[test]
129    fn evaluate_spline_extrapolate_before() {
130        let locs = [0.0_f32, 1.0];
131        let derivs = [2.0_f32, 0.0];
132        let result = evaluate_spline(&locs, &derivs, -1.0, |i| [0.0, 1.0][i]);
133        // value_at(0) + derivative[0] * (input - location[0]) = 0.0 + 2.0 * (-1.0) = -2.0
134        assert!((result - (-2.0)).abs() < 1e-6);
135    }
136
137    #[test]
138    fn evaluate_spline_extrapolate_after() {
139        let locs = [0.0_f32, 1.0];
140        let derivs = [0.0_f32, 3.0];
141        let result = evaluate_spline(&locs, &derivs, 2.0, |i| [0.0, 1.0][i]);
142        // value_at(1) + derivative[1] * (input - location[1]) = 1.0 + 3.0 * 1.0 = 4.0
143        assert!((result - 4.0).abs() < 1e-6);
144    }
145
146    #[test]
147    fn evaluate_spline_constant_values() {
148        let locs = [0.0_f32, 1.0, 2.0];
149        let derivs = [0.0_f32, 0.0, 0.0];
150        let vals = [1.0_f32, 1.0, 1.0];
151        // All constant values with zero derivatives = flat spline
152        let result = evaluate_spline(&locs, &derivs, 0.5, |i| vals[i]);
153        assert!((result - 1.0).abs() < 1e-6);
154    }
155}
steel_worldgen/density/spline_eval.rs

steel_worldgen/density/
spline_eval.rs
