1use super::Srgba;
2
3pub const SHADOW_STEPS: u32 = 8;
4pub const SHADOW_REACH: f32 = 4.0;
5
6#[derive(Clone, Copy, Debug, PartialEq)]
7pub struct CurvePoint {
8 pub elevation: f32,
9 pub offset: f32,
10 pub sigma: f32,
11 pub alpha: f32,
12}
13
14#[derive(Clone, Debug, PartialEq)]
15pub struct ShadowCurve {
16 pub colour: Srgba,
17 pub points: Vec<CurvePoint>,
18}
19
20impl ShadowCurve {
21 pub fn new(colour: Srgba, mut points: Vec<CurvePoint>) -> Result<Self, String> {
22 if points.iter().any(|point| {
23 !(point.elevation > 0.0
24 && point.elevation.is_finite()
25 && point.offset.is_finite()
26 && point.sigma >= 0.0
27 && point.sigma.is_finite()
28 && (0.0..=1.0).contains(&point.alpha))
29 }) {
30 return Err(
31 "curve points need a positive elevation, a finite offset, a nonnegative sigma and an alpha in 0..=1"
32 .into(),
33 );
34 }
35 points.sort_by(|a, b| a.elevation.total_cmp(&b.elevation));
36 if points
37 .windows(2)
38 .any(|pair| pair[0].elevation == pair[1].elevation)
39 {
40 return Err("two curve points share an elevation".into());
41 }
42 Ok(Self { colour, points })
43 }
44
45 pub fn none() -> Self {
46 Self {
47 colour: Srgba::TRANSPARENT,
48 points: Vec::new(),
49 }
50 }
51
52 pub fn at(&self, elevation: f32) -> Option<CurvePoint> {
53 if !super::positive(elevation) {
54 return None;
55 }
56 let first = *self.points.first()?;
57 let zero = CurvePoint {
58 elevation: 0.0,
59 offset: 0.0,
60 sigma: 0.0,
61 alpha: 0.0,
62 };
63 let mut low = zero;
64 for &point in &self.points {
65 if elevation <= point.elevation {
66 let t = (elevation - low.elevation) / (point.elevation - low.elevation);
67 let mix = |a: f32, b: f32| a + (b - a) * t;
68 return Some(CurvePoint {
69 elevation,
70 offset: mix(low.offset, point.offset),
71 sigma: mix(low.sigma, point.sigma),
72 alpha: mix(low.alpha, point.alpha),
73 });
74 }
75 low = point;
76 }
77 let last = *self.points.last().unwrap_or(&first);
78 Some(CurvePoint { elevation, ..last })
79 }
80}
81
82pub fn erf(x: f32) -> f32 {
83 let t = 1.0 / (1.0 + 0.327_591_1 * x.abs());
84 let y = 1.0
85 - (((((1.061_405_4 * t - 1.453_152_1) * t) + 1.421_413_8) * t - 0.284_496_74) * t
86 + 0.254_829_6)
87 * t
88 * (-x * x).exp();
89 if x >= 0.0 { y } else { -y }
90}
91
92fn reach(half: [f32; 2], radius: f32, y: f32) -> f32 {
93 let delta = (half[1] - radius - y.abs()).min(0.0);
94 half[0] - radius + (radius * radius - delta * delta).max(0.0).sqrt()
95}
96
97fn row(x: f32, y: f32, half: [f32; 2], radii: [f32; 4], sigma: f32) -> f32 {
98 let (left, right) = if y < 0.0 {
99 (radii[0], radii[1])
100 } else {
101 (radii[3], radii[2])
102 };
103 let k = std::f32::consts::FRAC_1_SQRT_2 / sigma;
104 0.5 * (erf((x + reach(half, left, y)) * k) - erf((x - reach(half, right, y)) * k))
105}
106
107pub fn rounded_box(p: [f32; 2], half: [f32; 2], radii: [f32; 4], sigma: [f32; 2]) -> f32 {
108 let k = std::f32::consts::FRAC_1_SQRT_2 / sigma[1];
109 let low = (-half[1]).max(p[1] - SHADOW_REACH * sigma[1]);
110 let high = half[1].min(p[1] + SHADOW_REACH * sigma[1]);
111 if high <= low {
112 return 0.0;
113 }
114 let first = (low / half[1]).clamp(-1.0, 1.0).asin();
115 let last = (high / half[1]).clamp(-1.0, 1.0).asin();
116 let step = (last - first) / SHADOW_STEPS as f32;
117 let mut total = 0.0;
118 let mut edge = erf((p[1] - low) * k);
119 for index in 0..SHADOW_STEPS {
120 let a = first + step * index as f32;
121 let b = a + step;
122 let next = erf((p[1] - half[1] * b.sin()) * k);
123 let weight = 0.5 * (edge - next);
124 edge = next;
125 total += weight * row(p[0], half[1] * ((a + b) * 0.5).sin(), half, radii, sigma[0]);
126 }
127 total
128}
129
130#[cfg(test)]
131mod tests {
132 use super::*;
133
134 fn exact(p: [f32; 2], half: [f32; 2], radii: [f32; 4], sigma: [f32; 2]) -> f32 {
135 let gauss = |x: f64, s: f64| {
136 (-(x * x) / (2.0 * s * s)).exp() / (s * (std::f64::consts::TAU).sqrt())
137 };
138 let n = 400;
139 let mut total = 0.0f64;
140 let (hx, hy) = (f64::from(half[0]), f64::from(half[1]));
141 let radii = radii.map(f64::from);
142 let reach = |r: f64, v: f64| {
143 let delta = (hy - r - v.abs()).min(0.0);
144 hx - r + (r * r - delta * delta).max(0.0).sqrt()
145 };
146 for i in 0..n {
147 let v = -hy + (i as f64 + 0.5) * 2.0 * hy / n as f64;
148 let (rl, rr) = if v < 0.0 {
149 (radii[0], radii[1])
150 } else {
151 (radii[3], radii[2])
152 };
153 let (left, right) = (reach(rl, v), reach(rr, v));
154 let width = left + right;
155 for j in 0..n {
156 let u = -left + (j as f64 + 0.5) * width / n as f64;
157 total += gauss(f64::from(p[0]) - u, f64::from(sigma[0]))
158 * gauss(f64::from(p[1]) - v, f64::from(sigma[1]))
159 * (width / n as f64)
160 * (2.0 * hy / n as f64);
161 }
162 }
163 total as f32
164 }
165
166 #[test]
167 fn erf_is_within_its_published_bound() {
168 for i in -40..=40 {
169 let x = i as f64 * 0.1;
170 let reference = libm_erf(x);
171 assert!((f64::from(erf(x as f32)) - reference).abs() < 1e-6, "{x}");
172 }
173 }
174
175 fn libm_erf(x: f64) -> f64 {
176 let n = 20000;
177 let step = x / n as f64;
178 let mut sum = 0.0;
179 for i in 0..n {
180 let t = (i as f64 + 0.5) * step;
181 sum += (-t * t).exp();
182 }
183 sum * step * 2.0 / std::f64::consts::PI.sqrt()
184 }
185
186 #[test]
187 fn the_shadow_matches_a_brute_force_blur() {
188 let cases = [
189 ([90.0, 30.0], [13.2; 4], [5.0, 5.0]),
190 ([170.0, 300.0], [30.0; 4], [16.0, 16.0]),
191 ([24.0, 24.0], [24.0; 4], [14.0, 14.0]),
192 ([60.0, 20.0], [0.0; 4], [5.0, 9.0]),
193 ([8.0, 8.0], [8.0; 4], [10.0, 10.0]),
194 ([40.0, 40.0], [40.0; 4], [4.0, 4.0]),
195 ([90.0, 30.0], [0.0, 24.0, 6.0, 12.0], [5.0, 5.0]),
196 ([60.0, 40.0], [30.0, 0.0, 30.0, 0.0], [8.0, 6.0]),
197 ([50.0, 50.0], [50.0, 0.0, 0.0, 0.0], [10.0, 10.0]),
198 ];
199 let mut worst = 0.0f32;
200 for (half, radii, sigma) in cases {
201 let mut case_worst = (0.0f32, [0.0f32; 2]);
202 for i in -6..=6 {
203 for j in -6..=6 {
204 let p = [
205 i as f32 * (half[0] + 3.0 * sigma[0]) / 6.0,
206 j as f32 * (half[1] + 3.0 * sigma[1]) / 6.0,
207 ];
208 let error =
209 (rounded_box(p, half, radii, sigma) - exact(p, half, radii, sigma)).abs();
210 worst = worst.max(error);
211 if error > case_worst.0 {
212 case_worst = (error, p);
213 }
214 }
215 }
216 eprintln!(
217 "{half:?} r {radii:?} sigma {sigma:?}: worst {} at {:?}",
218 case_worst.0, case_worst.1
219 );
220 }
221 assert!(worst < 8e-3, "worst {worst}");
222 }
223
224 #[test]
225 fn the_curve_is_piecewise_linear_from_zero() {
226 let curve = ShadowCurve::new(
227 Srgba::hex(0x2b3442),
228 vec![
229 CurvePoint {
230 elevation: 2.0,
231 offset: 12.0,
232 sigma: 16.0,
233 alpha: 0.26,
234 },
235 CurvePoint {
236 elevation: 1.0,
237 offset: 4.0,
238 sigma: 5.0,
239 alpha: 0.22,
240 },
241 ],
242 )
243 .unwrap();
244 assert_eq!(curve.at(0.0), None);
245 assert_eq!(curve.at(1.0).unwrap().offset, 4.0);
246 let mid = curve.at(1.5).unwrap();
247 assert!(
248 (mid.offset - 8.0).abs() < 1e-6
249 && (mid.sigma - 10.5).abs() < 1e-6
250 && (mid.alpha - 0.24).abs() < 1e-6
251 );
252 assert!((curve.at(0.5).unwrap().alpha - 0.11).abs() < 1e-6);
253 assert_eq!(curve.at(9.0).unwrap().sigma, 16.0);
254 assert!(
255 ShadowCurve::new(
256 Srgba::WHITE,
257 vec![CurvePoint {
258 elevation: 0.0,
259 offset: 0.0,
260 sigma: 1.0,
261 alpha: 0.5
262 }]
263 )
264 .is_err()
265 );
266 assert_eq!(ShadowCurve::none().at(3.0), None);
267 }
268}