motion_core/noise/
advanced.rs

1//! Cellular, simplex, evolving three-dimensional and warped fractal noise.
2
3use super::*;
4
5/// Euclidean distances to the nearest two seeded cell features, in increasing order.
6/// One feature per unit cell; a 5x5 search includes both nearest neighbors.
7/// Unsupported coordinates return zero distances.
8pub fn worley2(point: [f64; 2], seed: u32) -> (f64, f64) {
9    if !lattice_supported(point) {
10        return (0.0, 0.0);
11    }
12    let cell = point.map(f64::floor);
13    let fraction = [point[0] - cell[0], point[1] - cell[1]];
14    let mut nearest = [f64::INFINITY; 2];
15    for vertical in -2i32..=2 {
16        for horizontal in -2i32..=2 {
17            let hash = hash2(
18                (cell[0] as i32).wrapping_add(horizontal),
19                (cell[1] as i32).wrapping_add(vertical),
20                seed,
21            );
22            let feature = [
23                f64::from(hash32(hash ^ 0xa511_e9b3)) / 4294967296.0,
24                f64::from(hash32(hash ^ 0x63d8_3595)) / 4294967296.0,
25            ];
26            let distance = (f64::from(horizontal) + feature[0] - fraction[0]).powi(2)
27                + (f64::from(vertical) + feature[1] - fraction[1]).powi(2);
28            if distance < nearest[0] {
29                nearest = [distance, nearest[0]];
30            } else if distance < nearest[1] {
31                nearest[1] = distance;
32            }
33        }
34    }
35    (libm::sqrt(nearest[0]), libm::sqrt(nearest[1]))
36}
37
38/// Seeded triangular-lattice simplex noise with compact fourth-power kernels.
39/// Unsupported coordinates or skewed cells return zero.
40pub fn simplex2(point: [f64; 2], seed: u32) -> f64 {
41    if !lattice_supported(point) {
42        return 0.0;
43    }
44    let skew = (point[0] + point[1]) * 0.3660254037844386;
45    let cell = [(point[0] + skew).floor(), (point[1] + skew).floor()];
46    if !lattice_supported(cell) {
47        return 0.0;
48    }
49    let unskew = (cell[0] + cell[1]) * 0.21132486540518713;
50    let local = [point[0] - (cell[0] - unskew), point[1] - (cell[1] - unskew)];
51    let middle = if local[0] > local[1] {
52        [1i32, 0]
53    } else {
54        [0, 1]
55    };
56    let corner = |offset: [i32; 2], shift: f64| {
57        let delta = [
58            local[0] - f64::from(offset[0]) + shift,
59            local[1] - f64::from(offset[1]) + shift,
60        ];
61        let weight = (0.5 - delta[0] * delta[0] - delta[1] * delta[1]).max(0.0);
62        let gradient = GRADIENTS[(hash2(
63            (cell[0] as i32).wrapping_add(offset[0]),
64            (cell[1] as i32).wrapping_add(offset[1]),
65            seed,
66        ) % 12) as usize];
67        weight.powi(4) * (gradient[0] * delta[0] + gradient[1] * delta[1])
68    };
69    70.0 * (corner([0, 0], 0.0)
70        + corner(middle, 0.21132486540518713)
71        + corner([1, 1], 0.42264973081037427))
72}
73
74/// Quintic three-dimensional Perlin, with the third coordinate usable as time.
75/// Twelve normalized edge gradients; unsupported lattice coordinates return zero.
76pub fn perlin3(point: [f64; 3], seed: u32) -> f64 {
77    if point.into_iter().any(|value| {
78        !value.is_finite() || value < f64::from(i32::MIN) || value >= f64::from(i32::MAX)
79    }) {
80        return 0.0;
81    }
82    let cell = point.map(f64::floor);
83    let fraction: [f64; 3] = std::array::from_fn(|axis| point[axis] - cell[axis]);
84    let amount = fraction.map(|value| fade(value).0);
85    let gradients = [
86        [1.0, 1.0, 0.0],
87        [-1.0, 1.0, 0.0],
88        [1.0, -1.0, 0.0],
89        [-1.0, -1.0, 0.0],
90        [1.0, 0.0, 1.0],
91        [-1.0, 0.0, 1.0],
92        [1.0, 0.0, -1.0],
93        [-1.0, 0.0, -1.0],
94        [0.0, 1.0, 1.0],
95        [0.0, -1.0, 1.0],
96        [0.0, 1.0, -1.0],
97        [0.0, -1.0, -1.0],
98    ];
99    let mut result = 0.0;
100    for depth in 0i32..=1 {
101        for vertical in 0i32..=1 {
102            for horizontal in 0i32..=1 {
103                let offset = [horizontal, vertical, depth];
104                let gradient = gradients[(hash3(
105                    (cell[0] as i32).wrapping_add(horizontal),
106                    (cell[1] as i32).wrapping_add(vertical),
107                    (cell[2] as i32).wrapping_add(depth),
108                    seed,
109                ) % 12) as usize];
110                let dot = (0..3)
111                    .map(|axis| gradient[axis] * (fraction[axis] - f64::from(offset[axis])))
112                    .sum::<f64>()
113                    * std::f64::consts::FRAC_1_SQRT_2;
114                let weight = (0..3)
115                    .map(|axis| {
116                        if offset[axis] == 0 {
117                            1.0 - amount[axis]
118                        } else {
119                            amount[axis]
120                        }
121                    })
122                    .product::<f64>();
123                result += dot * weight;
124            }
125        }
126    }
127    result
128}
129
130/// Bounded displacement by two independent shared Perlin fractals.
131#[derive(Clone, Copy, Debug, PartialEq)]
132pub struct DomainWarp {
133    /// Displacement in input-coordinate units, in 0..=1000000.
134    pub strength: f64,
135    /// Frequency applied only to the displacement potential, in (0, 1000000].
136    pub frequency: f64,
137    /// Explicit displacement seed; the second axis uses a decorrelated seed.
138    pub seed: u32,
139    /// Bounded displacement octave recipe.
140    pub fractal: Fractal,
141}
142
143impl Default for DomainWarp {
144    fn default() -> Self {
145        Self {
146            strength: 0.5,
147            frequency: 1.0,
148            seed: 0,
149            fractal: Fractal::default(),
150        }
151    }
152}
153
154impl DomainWarp {
155    /// Validates finite bounded displacement parameters and its octave recipe.
156    pub fn validate(self) -> Result<(), crate::ValidationError> {
157        self.fractal.validate()?;
158        if !self.strength.is_finite()
159            || !(0.0..=1e6).contains(&self.strength)
160            || !self.frequency.is_finite()
161            || self.frequency <= 0.0
162            || self.frequency > 1e6
163        {
164            return Err(crate::ValidationError::invalid(
165                "noise",
166                "warp",
167                "Require strength 0..=1000000 and frequency (0,1000000]",
168            ));
169        }
170        Ok(())
171    }
172}
173
174/// Warps a point without changing the sampling algorithm applied afterwards.
175/// Rejects unsupported input, octave coordinates and displaced output.
176pub fn domain_warp(point: [f64; 2], warp: DomainWarp) -> Result<[f64; 2], crate::ValidationError> {
177    warp.validate()?;
178    if !lattice_supported(point) {
179        return Err(crate::ValidationError::invalid(
180            "noise",
181            "point",
182            "Warp coordinates exceed lattice support",
183        ));
184    }
185    if warp.strength == 0.0 {
186        return Ok(point);
187    }
188    let source = point.map(|value| value * warp.frequency);
189    let horizontal = fbm(|point| perlin2(point, warp.seed), source, warp.fractal)?;
190    let vertical = fbm(
191        |point| perlin2(point, warp.seed ^ 0xa511_e9b3),
192        source,
193        warp.fractal,
194    )?;
195    let output = [
196        point[0] + horizontal * warp.strength,
197        point[1] + vertical * warp.strength,
198    ];
199    if !lattice_supported(output) {
200        return Err(crate::ValidationError::invalid(
201            "noise",
202            "point",
203            "Displaced coordinates exceed lattice support",
204        ));
205    }
206    Ok(output)
207}
208
209/// Amplitude-normalized absolute Perlin fractal in `[0, 1]`.
210pub fn turbulence2(
211    point: [f64; 2],
212    seed: u32,
213    fractal: Fractal,
214) -> Result<f64, crate::ValidationError> {
215    fbm(|point| perlin2(point, seed).abs(), point, fractal)
216}
217
218/// Amplitude-normalized squared ridge fractal in `[0, 1]`.
219pub fn ridged2(
220    point: [f64; 2],
221    seed: u32,
222    fractal: Fractal,
223) -> Result<f64, crate::ValidationError> {
224    fbm(
225        |point| (1.0 - perlin2(point, seed).abs()).max(0.0).powi(2),
226        point,
227        fractal,
228    )
229}
230
231#[cfg(test)]
232mod tests {
233    use super::*;
234
235    #[test]
236    fn worley_nearest_pair_matches_larger_reference_search() {
237        for index in 0..257 {
238            let point = [
239                f64::from(index % 17) * 0.173 - 1.2,
240                f64::from(index / 17) * 0.137 - 0.9,
241            ];
242            let cell = point.map(f64::floor);
243            let mut distances = Vec::new();
244            for vertical in -4i32..=4 {
245                for horizontal in -4i32..=4 {
246                    let hash = hash2(cell[0] as i32 + horizontal, cell[1] as i32 + vertical, 42);
247                    let feature = [
248                        cell[0]
249                            + f64::from(horizontal)
250                            + f64::from(hash32(hash ^ 0xa511_e9b3)) / 4294967296.0,
251                        cell[1]
252                            + f64::from(vertical)
253                            + f64::from(hash32(hash ^ 0x63d8_3595)) / 4294967296.0,
254                    ];
255                    distances.push(libm::hypot(feature[0] - point[0], feature[1] - point[1]));
256                }
257            }
258            distances.sort_by(f64::total_cmp);
259            let actual = worley2(point, 42);
260            assert!((actual.0 - distances[0]).abs() < 1e-12);
261            assert!((actual.1 - distances[1]).abs() < 1e-12);
262            assert!(actual.0 >= 0.0 && actual.0 <= actual.1);
263        }
264    }
265
266    #[test]
267    fn simplex_perlin3_and_warp_are_finite_continuous_and_seeded() {
268        assert_eq!(simplex2([0.0; 2], 42), 0.0);
269        assert_eq!(perlin3([1.0, -2.0, 3.0], 42), 0.0);
270        let mut sum = 0.0;
271        let mut energy = 0.0;
272        for index in 0..4096 {
273            let point = [
274                f64::from(index % 64) * 0.173 - 4.0,
275                f64::from(index / 64) * 0.137 - 3.0,
276            ];
277            let value = simplex2(point, 42);
278            assert!(value.is_finite() && value.abs() <= 1.0);
279            sum += value;
280            energy += value * value;
281            assert!((value - simplex2([point[0] + 1e-7, point[1]], 42)).abs() < 1e-5);
282            assert!(perlin3([point[0], point[1], 0.37], 42).abs() <= 1.0);
283        }
284        assert!((sum / 4096.0).abs() < 0.06 && energy / 4096.0 > 0.01);
285        let point = [0.37, -1.23];
286        assert_ne!(
287            perlin3([point[0], point[1], 0.3], 42),
288            perlin3([point[0], point[1], 0.4], 42)
289        );
290        assert_ne!(simplex2(point, 42), simplex2(point, 43));
291        let warp = DomainWarp {
292            seed: 42,
293            ..DomainWarp::default()
294        };
295        assert_eq!(
296            domain_warp(
297                point,
298                DomainWarp {
299                    strength: 0.0,
300                    ..warp
301                }
302            )
303            .unwrap(),
304            point
305        );
306        let output = domain_warp(point, warp).unwrap();
307        assert_eq!(output, domain_warp(point, warp).unwrap());
308        assert_ne!(output, point);
309        assert_ne!(
310            output,
311            domain_warp(point, DomainWarp { seed: 43, ..warp }).unwrap()
312        );
313        assert!(
314            domain_warp(
315                point,
316                DomainWarp {
317                    strength: f64::NAN,
318                    ..warp
319                }
320            )
321            .is_err()
322        );
323        assert!(domain_warp([f64::INFINITY, 0.0], warp).is_err());
324        for value in [
325            turbulence2(point, 42, Fractal::default()).unwrap(),
326            ridged2(point, 42, Fractal::default()).unwrap(),
327        ] {
328            assert!((0.0..=1.0).contains(&value));
329        }
330    }
331
332    #[test]
333    fn curl_is_divergence_free_and_depends_on_both_coordinates() {
334        let point = [0.37, -1.23];
335        let delta = 1e-5;
336        let horizontal = (curl2([point[0] + delta, point[1]], 42)[0]
337            - curl2([point[0] - delta, point[1]], 42)[0])
338            / (2.0 * delta);
339        let vertical = (curl2([point[0], point[1] + delta], 42)[1]
340            - curl2([point[0], point[1] - delta], 42)[1])
341            / (2.0 * delta);
342        assert!((horizontal + vertical).abs() < 1e-7);
343        assert_ne!(
344            curl2(point, 42)[0],
345            curl2([point[0] + 0.1, point[1]], 42)[0]
346        );
347        assert_ne!(
348            curl2(point, 42)[1],
349            curl2([point[0], point[1] + 0.1], 42)[1]
350        );
351    }
352}
motion_core/noise/advanced.rs

motion_core/noise/
advanced.rs
