doryen_extra/noise/algorithms/
wavelet.rs

1/* BSD 3-Clause License
2 *
3 * Copyright © 2019, Alexander Krivács Schrøder <alexschrod@gmail.com>.
4 * Copyright © 2008-2019, Jice and the libtcod contributors.
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions are met:
9 *
10 * 1. Redistributions of source code must retain the above copyright notice,
11 *    this list of conditions and the following disclaimer.
12 *
13 * 2. Redistributions in binary form must reproduce the above copyright notice,
14 *    this list of conditions and the following disclaimer in the documentation
15 *    and/or other materials provided with the distribution.
16 *
17 * 3. Neither the name of the copyright holder nor the names of its
18 *    contributors may be used to endorse or promote products derived from
19 *    this software without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 */
33
34use crate::noise::algorithms::AlgorithmInitializer;
35use crate::noise::Algorithm;
36use crate::random::algorithms::Algorithm as RandomAlgorithm;
37use crate::random::{Random, Rng};
38use crate::util::FloorRem;
39use derivative::Derivative;
40use std::mem::MaybeUninit;
41
42#[allow(dead_code)]
43const WAVELET_TILE_SIZE: usize = 32;
44const WAVELET_TILE_SIZE_SQUARED: usize = WAVELET_TILE_SIZE * WAVELET_TILE_SIZE;
45const WAVELET_TILE_SIZE_CUBED: usize = WAVELET_TILE_SIZE_SQUARED * WAVELET_TILE_SIZE;
46const WAVELET_ARAD: usize = 16;
47
48const WAVELET_SCALE: f32 = 2.0;
49
50/* wavelet noise, adapted from Robert L. Cook and Tony Derose 'Wavelet noise' paper */
51
52/// Wavelet noise algorithm.
53#[derive(Clone, Derivative)]
54#[derivative(Debug)]
55pub struct Wavelet {
56    dimensions: usize,
57    #[derivative(Debug = "ignore")]
58    tile_data: Box<[f32; WAVELET_TILE_SIZE_CUBED]>,
59}
60
61impl Algorithm for Wavelet {
62    fn new<R: RandomAlgorithm>(dimensions: usize, initializer: AlgorithmInitializer<R>) -> Self {
63        let mut random = initializer.random;
64        let tile_data = WaveletTileData::initialize(&mut random);
65        Self {
66            dimensions,
67            tile_data,
68        }
69    }
70
71    #[allow(clippy::many_single_char_names)]
72    fn generate(&self, f: &[f32]) -> f32 {
73        if self.dimensions > 3 {
74            panic!("Wavelet noise only supports up to 3 dimensions");
75        }
76
77        let mut pf = [0.0; 3];
78        for (pfe, &fe) in Iterator::zip(
79            pf.iter_mut().take(self.dimensions),
80            f.iter().take(self.dimensions),
81        ) {
82            *pfe = fe * WAVELET_SCALE;
83        }
84
85        let mut mid = [0; 3];
86        let mut w = [[0.0; 3]; 3];
87        let mut t;
88        for i in 0..3 {
89            mid[i] = (pf[i] - 0.5).ceil() as i32;
90            t = mid[i] as f32 - (pf[i] - 0.5);
91            w[i][0] = t * t * 0.5;
92            w[i][2] = (1.0 - t) * (1.0 - t) * 0.5;
93            w[i][1] = 1.0 - w[i][0] - w[i][2];
94        }
95
96        let mut c = [0; 3];
97        let mut result = 0.0;
98        let mid = mid;
99        for p2 in -1..=1 {
100            for p1 in -1..=1 {
101                for p0 in -1..=1 {
102                    let mut weight = 1.0;
103                    for i in 0..3 {
104                        let p = match i {
105                            0 => p0,
106                            1 => p1,
107                            2 => p2,
108                            _ => unreachable!(),
109                        };
110
111                        c[i] = (mid[i].wrapping_add(p)).floor_modulo(WAVELET_TILE_SIZE as i32);
112                        weight *= w[i][(p + 1) as usize];
113                    }
114                    result += weight
115                        * self.tile_data[c[2] as usize * WAVELET_TILE_SIZE_SQUARED
116                            + c[1] as usize * WAVELET_TILE_SIZE
117                            + c[0] as usize];
118                }
119            }
120        }
121
122        result.max(-1.0).min(1.0)
123    }
124}
125
126pub(crate) struct WaveletTileData;
127
128impl WaveletTileData {
129    pub(crate) fn initialize<R: RandomAlgorithm>(
130        random: &mut Random<R>,
131    ) -> Box<[f32; WAVELET_TILE_SIZE_CUBED]> {
132        let mut noise = Self::generate_noise(random);
133        let mut temp1 = [0.0; WAVELET_TILE_SIZE_CUBED];
134        let mut temp2 = [0.0; WAVELET_TILE_SIZE_CUBED];
135
136        for iy in 0..WAVELET_TILE_SIZE {
137            for iz in 0..WAVELET_TILE_SIZE {
138                let i = iy * WAVELET_TILE_SIZE + iz * WAVELET_TILE_SIZE_SQUARED;
139                Self::downsample(&noise[i..], &mut temp1[i..], 1);
140                Self::upsample(&temp1[i..], &mut temp2[i..], 1);
141            }
142        }
143        for ix in 0..WAVELET_TILE_SIZE {
144            for iz in 0..WAVELET_TILE_SIZE {
145                let i = ix + iz * WAVELET_TILE_SIZE_SQUARED;
146                Self::downsample(&temp2[i..], &mut temp1[i..], WAVELET_TILE_SIZE);
147                Self::upsample(&temp1[i..], &mut temp2[i..], WAVELET_TILE_SIZE);
148            }
149        }
150        for ix in 0..WAVELET_TILE_SIZE {
151            for iy in 0..WAVELET_TILE_SIZE {
152                let i = ix + iy * WAVELET_TILE_SIZE;
153                Self::downsample(&temp2[i..], &mut temp1[i..], WAVELET_TILE_SIZE_SQUARED);
154                Self::upsample(&temp1[i..], &mut temp2[i..], WAVELET_TILE_SIZE_SQUARED);
155            }
156        }
157        for (n, &t) in Iterator::zip(noise.iter_mut(), temp2.iter()) {
158            *n -= t;
159        }
160        let mut offset = WAVELET_TILE_SIZE / 2;
161        if offset & 1 == 0 {
162            offset += 1
163        }
164        let mut i = 0;
165        for ix in 0..WAVELET_TILE_SIZE {
166            for iy in 0..WAVELET_TILE_SIZE {
167                for iz in 0..WAVELET_TILE_SIZE {
168                    temp1[i] = noise[((ix + offset) % (WAVELET_TILE_SIZE))
169                        + ((iy + offset) % (WAVELET_TILE_SIZE)) * WAVELET_TILE_SIZE
170                        + ((iz + offset) % (WAVELET_TILE_SIZE)) * WAVELET_TILE_SIZE_SQUARED];
171                    i += 1;
172                }
173            }
174        }
175        for (n, &t) in Iterator::zip(noise.iter_mut(), temp1.iter()) {
176            *n += t;
177        }
178
179        noise
180    }
181
182    #[allow(unsafe_code)]
183    fn generate_noise<R: RandomAlgorithm>(
184        random: &mut Random<R>,
185    ) -> Box<[f32; WAVELET_TILE_SIZE_CUBED]> {
186        let mut noise: Box<[MaybeUninit<f32>; WAVELET_TILE_SIZE_CUBED]> =
187            Box::new(unsafe { MaybeUninit::uninit().assume_init() });
188        for n in noise.iter_mut() {
189            unsafe {
190                n.as_mut_ptr().write(random.get_f32(-1.0, 1.0));
191            }
192        }
193
194        unsafe { Box::from_raw(Box::into_raw(noise).cast::<[f32; WAVELET_TILE_SIZE_CUBED]>()) }
195    }
196
197    fn downsample(from: &[f32], to: &mut [f32], stride: usize) {
198        const A_COEFFICIENTS: [f32; 2 * WAVELET_ARAD] = [
199            0.000_334, -0.001_528, 0.000_410, 0.003_545, -0.000_938, -0.008_233, 0.002_172,
200            0.019_120, -0.005_040, -0.044_412, 0.011_655, 0.103_311, -0.025_936, -0.243_780,
201            0.033_979, 0.655_340, 0.655_340, 0.033_979, -0.243_780, -0.025_936, 0.103_311,
202            0.011_655, -0.044_412, -0.005_040, 0.019_120, 0.002_172, -0.008_233, -0.000_938,
203            0.003_546, 0.000_410, -0.001_528, 0.000_334,
204        ];
205
206        for i in 0..WAVELET_TILE_SIZE as isize / 2 {
207            to[i as usize * stride] = 0.0;
208            for k in 2 * i - WAVELET_ARAD as isize..2 * i + WAVELET_ARAD as isize {
209                to[i as usize * stride] += A_COEFFICIENTS[(16 + k - 2 * i) as usize]
210                    * from[k.floor_modulo(WAVELET_TILE_SIZE as isize) as usize * stride];
211            }
212        }
213    }
214
215    fn upsample(from: &[f32], to: &mut [f32], stride: usize) {
216        const P_COEFFICIENTS: [f32; 4] = [0.25, 0.75, 0.75, 0.25];
217
218        for i in 0..WAVELET_TILE_SIZE as isize {
219            to[i as usize * stride] = 0.0;
220            for k in i / 2..=i / 2 {
221                to[i as usize * stride] += P_COEFFICIENTS[(2 + i - 2 * k) as usize]
222                    * from[k.floor_modulo(WAVELET_TILE_SIZE as isize / 2) as usize * stride];
223            }
224        }
225    }
226}