fast_voxel_traversal/
raycast_3d.rs

1use glam::{IVec3, Vec3};
2
3#[derive(Debug, Default, Clone, Copy, PartialEq, Eq)]
4pub struct BoundingVolume3 {
5    pub size: (i32, i32, i32),
6}
7
8impl BoundingVolume3 {
9    pub fn traverse_ray(&self, ray: Ray3) -> VoxelRay3Iterator {
10        VoxelRay3Iterator::new(self.clone(), ray)
11    }
12
13    #[inline(always)]
14    pub(crate) fn contains_point(&self, point: IVec3) -> bool {
15        point.cmpge(IVec3::ZERO).all() && point.cmplt(self.size.into()).all()
16    }
17}
18
19#[derive(Debug, Clone, Copy)]
20pub struct Ray3 {
21    pub origin: (f32, f32, f32),
22    pub direction: (f32, f32, f32),
23    pub length: f32,
24}
25
26#[derive(Debug, Clone, Copy)]
27pub struct Ray3hit {
28    pub distance: f32,
29    pub voxel: (i32, i32, i32),
30    pub normal: Option<(i32, i32, i32)>,
31}
32
33#[derive(Debug, Default, Clone, Copy)]
34pub struct VoxelRay3Iterator {
35    volume: BoundingVolume3,
36    max_d: f32,
37    i: IVec3,
38    step: IVec3,
39    delta: Vec3,
40    dist: Vec3,
41    t_max: Vec3,
42    t: f32,
43    norm: Option<IVec3>,
44    done: bool,
45}
46
47// Based on https://github.com/fenomas/fast-voxel-raycast/blob/master/index.js
48impl VoxelRay3Iterator {
49    pub fn new(volume: BoundingVolume3, ray: Ray3) -> Self {
50        let mut p = Vec3::from(ray.origin);
51
52        // Normalize direction vector
53        let d = Vec3::from(ray.direction).normalize();
54
55        // How long we have traveled thus far (modified by initial 'jump to volume bounds').
56        let mut t = 0.0;
57
58        // If the point it outside the chunk, AABB test to 'jump ahead'.
59        if !volume.contains_point(p.floor().as_ivec3()) {
60            // First AABB test the chunk bounds
61            let aabb = test_aabb_of_chunk(volume, p, d, ray.length);
62
63            // Chunk AABB test failed, no way we could intersect a voxel.
64            if aabb.is_none() {
65                return Self {
66                    done: true,
67                    ..Default::default()
68                };
69            }
70
71            let aabb = aabb.unwrap();
72
73            // Back the hit off at least 1 voxel
74            p = aabb - d * 2.0;
75
76            // Set t to the already traveled distance.
77            t += (p - aabb).length() - 2.0;
78        }
79
80        // Max distance we can travel. This is either the ray length, or the current `t` plus the
81        // corner to corner length of the voxel volume.
82        let max_d = f32::min(
83            ray.length,
84            t + IVec3::from(volume.size).as_vec3().length() + 2.0,
85        );
86
87        // The starting voxel for the raycast.
88        let i = p.floor().as_ivec3();
89
90        // The directionVec we are stepping for each component.
91        let step = d.signum().as_ivec3();
92
93        // Just abs(Vec3::ONE / d) but acounts for zeros in the distance vector.
94        let delta = (Vec3::new(
95            if d.x.abs() < f32::EPSILON {
96                f32::INFINITY
97            } else {
98                1.0 / d.x
99            },
100            if d.y.abs() < f32::EPSILON {
101                f32::INFINITY
102            } else {
103                1.0 / d.y
104            },
105            if d.z.abs() < f32::EPSILON {
106                f32::INFINITY
107            } else {
108                1.0 / d.z
109            },
110        ))
111        .abs();
112
113        let dist = Vec3::new(
114            if step.x > 0 {
115                i.x as f32 + 1.0 - p.x
116            } else {
117                p.x - i.x as f32
118            },
119            if step.y > 0 {
120                i.y as f32 + 1.0 - p.y
121            } else {
122                p.y - i.y as f32
123            },
124            if step.z > 0 {
125                i.z as f32 + 1.0 - p.z
126            } else {
127                p.z - i.z as f32
128            },
129        );
130
131        // The nearest voxel boundary.
132        let t_max = Vec3::new(
133            if delta.x < f32::INFINITY {
134                delta.x * dist.x
135            } else {
136                f32::INFINITY
137            },
138            if delta.y < f32::INFINITY {
139                delta.y * dist.y
140            } else {
141                f32::INFINITY
142            },
143            if delta.z < f32::INFINITY {
144                delta.z * dist.z
145            } else {
146                f32::INFINITY
147            },
148        );
149
150        Self {
151            volume,
152            max_d,
153            i,
154            step,
155            delta,
156            dist,
157            t_max,
158            t,
159            norm: None,
160            done: false,
161        }
162    }
163}
164
165impl Iterator for VoxelRay3Iterator {
166    type Item = Ray3hit;
167
168    fn next(&mut self) -> Option<Self::Item> {
169        if self.done {
170            return None;
171        }
172
173        while self.t <= self.max_d {
174            // Test if the current traverse is within the volume.
175            let mut hit = None;
176            if self.volume.contains_point(self.i) {
177                hit = Some(Ray3hit {
178                    distance: self.t,
179                    voxel: self.i.into(),
180                    normal: self.norm.map(|n| n.into()),
181                });
182            }
183
184            // Select the smallest t_max
185            if self.t_max.x < self.t_max.y {
186                if self.t_max.x < self.t_max.z {
187                    self.i.x += self.step.x;
188                    self.t = self.t_max.x;
189                    self.t_max.x += self.delta.x;
190                    self.norm = Some(IVec3::new(-self.step.x, 0, 0));
191                } else {
192                    self.i.z += self.step.z;
193                    self.t = self.t_max.z;
194                    self.t_max.z += self.delta.z;
195                    self.norm = Some(IVec3::new(0, 0, -self.step.z));
196                }
197            } else {
198                if self.t_max.y < self.t_max.z {
199                    self.i.y += self.step.y;
200                    self.t = self.t_max.y;
201                    self.t_max.y += self.delta.y;
202                    self.norm = Some(IVec3::new(0, -self.step.y, 0));
203                } else {
204                    self.i.z += self.step.z;
205                    self.t = self.t_max.z;
206                    self.t_max.z += self.delta.z;
207                    self.norm = Some(IVec3::new(0, 0, -self.step.z));
208                }
209            }
210
211            if hit.is_some() {
212                return hit;
213            }
214        }
215
216        self.done = true;
217        return None;
218    }
219}
220
221fn test_aabb_of_chunk(
222    volume: BoundingVolume3,
223    from: Vec3,
224    direction: Vec3,
225    distance: f32,
226) -> Option<Vec3> {
227    let min = Vec3::ZERO;
228    let max = IVec3::from(volume.size).as_vec3();
229    let mut t = Vec3::ZERO;
230
231    for i in 0..3 {
232        if direction[i] > 0.0 {
233            t[i] = (min[i] - from[i]) / direction[i];
234        } else {
235            t[i] = (max[i] - from[i]) / direction[i];
236        }
237    }
238
239    let mi = if t[0] > t[1] {
240        if t[0] > t[2] {
241            0
242        } else {
243            2
244        }
245    } else {
246        if t[1] > t[2] {
247            1
248        } else {
249            2
250        }
251    };
252
253    if t[mi] >= 0.0 && t[mi] <= distance {
254        // The intersect point (distance along the ray).
255        let pt = from + direction * t[mi];
256
257        // The other two value that need to be checked
258        let o1 = (mi + 1) % 3;
259        let o2 = (mi + 2) % 3;
260
261        if pt[o1] >= min[o1] && pt[o1] <= max[o1] && pt[o2] >= min[o2] && pt[o2] <= max[o2] {
262            return Some(pt);
263        }
264    }
265
266    // AABB test failed.
267    return None;
268}