packed_spatial_index/ray.rs
1//! Finite ray segments for ray/AABB traversal over the packed indexes.
2//!
3//! A ray is defined by an origin, an unnormalized direction, and a maximum ray
4//! parameter `max_distance` (the segment covers `origin + t * dir` for
5//! `t in [0, max_distance]`). Box intersections use the standard slab test with
6//! precomputed reciprocal directions; axis-parallel rays (a direction component
7//! that is exactly zero) are handled explicitly so a ray lying exactly on a box
8//! face still hits.
9
10use crate::geometry::{Box2D, Box3D, Point2D, Point3D};
11use crate::triangle::{Triangle3, TriangleHit};
12
13/// Finite 2D ray segment used by `raycast` searches.
14///
15/// # Example
16///
17/// ```
18/// use packed_spatial_index::{Box2D, Point2D, Ray2D};
19///
20/// let ray = Ray2D::new(Point2D::new(-1.0, 0.5), 1.0, 0.0, 10.0);
21/// assert!(ray.intersects_box(Box2D::new(0.0, 0.0, 1.0, 1.0)));
22/// assert_eq!(ray.enter_t(Box2D::new(0.0, 0.0, 1.0, 1.0)), Some(1.0));
23/// ```
24#[derive(Clone, Copy, Debug, PartialEq)]
25pub struct Ray2D {
26 /// Ray origin.
27 pub origin: Point2D,
28 /// X component of the ray direction.
29 pub dir_x: f64,
30 /// Y component of the ray direction.
31 pub dir_y: f64,
32 // Precomputed to remove one division from every slab test.
33 pub(crate) inv_dir_x: f64,
34 pub(crate) inv_dir_y: f64,
35 /// Maximum ray parameter to consider.
36 pub max_distance: f64,
37}
38
39impl Ray2D {
40 /// Create a finite ray segment covering `origin + t * dir` for
41 /// `t in [0, max_distance]`.
42 ///
43 /// The direction does not need to be normalized. `max_distance` and every
44 /// returned entry `t` are in units of the direction length, so the
45 /// Euclidean distance to a hit is `t * hypot(dir_x, dir_y)`; normalize the
46 /// direction (length 1) if you want `t` and `max_distance` in world units.
47 ///
48 /// A fully zero direction (`dir_x == dir_y == 0.0`) is a point probe: it
49 /// hits only boxes that contain `origin`, all at `t == 0.0`. Direction
50 /// components should be finite; a `NaN` direction produces unspecified
51 /// results.
52 #[inline]
53 pub const fn new(origin: Point2D, dir_x: f64, dir_y: f64, max_distance: f64) -> Self {
54 Self {
55 origin,
56 dir_x,
57 dir_y,
58 inv_dir_x: 1.0 / dir_x,
59 inv_dir_y: 1.0 / dir_y,
60 max_distance,
61 }
62 }
63
64 /// `true` if any direction component is exactly zero (an axis-parallel ray). The
65 /// vectorized slab test uses `1/dir = inf` and is not NaN-safe at a box face, so
66 /// such rays take a masked path.
67 #[inline]
68 #[allow(dead_code)] // used by the SIMD raycast paths only
69 pub(crate) fn has_zero_direction(self) -> bool {
70 self.dir_x == 0.0 || self.dir_y == 0.0
71 }
72
73 /// `true` when the ray segment touches `bounds` (edges inclusive).
74 #[inline]
75 pub fn intersects_box(self, bounds: Box2D) -> bool {
76 if self.max_distance < 0.0 || self.max_distance.is_nan() {
77 return false;
78 }
79 let mut t_min: f64 = 0.0;
80 let mut t_max = self.max_distance;
81 slab(
82 self.origin.x,
83 self.dir_x,
84 self.inv_dir_x,
85 bounds.min_x,
86 bounds.max_x,
87 &mut t_min,
88 &mut t_max,
89 ) && slab(
90 self.origin.y,
91 self.dir_y,
92 self.inv_dir_y,
93 bounds.min_y,
94 bounds.max_y,
95 &mut t_min,
96 &mut t_max,
97 )
98 }
99
100 /// Entry parameter `t` where the ray segment first touches `bounds` (`0.0` if the
101 /// origin is inside), or `None` if the segment misses. Used by ordered closest-hit
102 /// traversal.
103 #[inline]
104 pub fn enter_t(self, bounds: Box2D) -> Option<f64> {
105 if self.max_distance < 0.0 || self.max_distance.is_nan() {
106 return None;
107 }
108 let mut t_min: f64 = 0.0;
109 let mut t_max = self.max_distance;
110 let hit = slab(
111 self.origin.x,
112 self.dir_x,
113 self.inv_dir_x,
114 bounds.min_x,
115 bounds.max_x,
116 &mut t_min,
117 &mut t_max,
118 ) && slab(
119 self.origin.y,
120 self.dir_y,
121 self.inv_dir_y,
122 bounds.min_y,
123 bounds.max_y,
124 &mut t_min,
125 &mut t_max,
126 );
127 hit.then_some(t_min)
128 }
129}
130
131/// Finite 3D ray segment used by `raycast` searches.
132///
133/// # Example
134///
135/// ```
136/// use packed_spatial_index::{Box3D, Point3D, Ray3D};
137///
138/// let ray = Ray3D::new(Point3D::new(-1.0, 0.5, 0.5), 1.0, 0.0, 0.0, 10.0);
139/// assert!(ray.intersects_box(Box3D::new(0.0, 0.0, 0.0, 1.0, 1.0, 1.0)));
140/// ```
141#[derive(Clone, Copy, Debug, PartialEq)]
142pub struct Ray3D {
143 /// Ray origin.
144 pub origin: Point3D,
145 /// X component of the ray direction.
146 pub dir_x: f64,
147 /// Y component of the ray direction.
148 pub dir_y: f64,
149 /// Z component of the ray direction.
150 pub dir_z: f64,
151 // Precomputed to remove one division from every slab test.
152 pub(crate) inv_dir_x: f64,
153 pub(crate) inv_dir_y: f64,
154 pub(crate) inv_dir_z: f64,
155 /// Maximum ray parameter to consider.
156 pub max_distance: f64,
157}
158
159impl Ray3D {
160 /// Create a finite ray segment covering `origin + t * dir` for
161 /// `t in [0, max_distance]`.
162 ///
163 /// The direction does not need to be normalized. `max_distance` and every
164 /// returned entry `t` are in units of the direction length, so the
165 /// Euclidean distance to a hit is `t * (dir_x.hypot(dir_y).hypot(dir_z))`;
166 /// normalize the direction (length 1) if you want `t` and `max_distance` in
167 /// world units.
168 ///
169 /// A fully zero direction (`dir_x == dir_y == dir_z == 0.0`) is a point
170 /// probe: it hits only boxes that contain `origin`, all at `t == 0.0`.
171 /// Direction components should be finite; a `NaN` direction produces
172 /// unspecified results.
173 #[inline]
174 pub const fn new(
175 origin: Point3D,
176 dir_x: f64,
177 dir_y: f64,
178 dir_z: f64,
179 max_distance: f64,
180 ) -> Self {
181 Self {
182 origin,
183 dir_x,
184 dir_y,
185 dir_z,
186 inv_dir_x: 1.0 / dir_x,
187 inv_dir_y: 1.0 / dir_y,
188 inv_dir_z: 1.0 / dir_z,
189 max_distance,
190 }
191 }
192
193 /// `true` if any direction component is exactly zero (an axis-parallel ray). The
194 /// vectorized slab test uses `1/dir = inf` and is not NaN-safe at a box face, so
195 /// such rays take a masked path.
196 #[inline]
197 #[allow(dead_code)] // used by the SIMD raycast paths only
198 pub(crate) fn has_zero_direction(self) -> bool {
199 self.dir_x == 0.0 || self.dir_y == 0.0 || self.dir_z == 0.0
200 }
201
202 /// `true` when the ray segment touches `bounds` (faces inclusive).
203 #[inline]
204 pub fn intersects_box(self, bounds: Box3D) -> bool {
205 if self.max_distance < 0.0 || self.max_distance.is_nan() {
206 return false;
207 }
208 let mut t_min: f64 = 0.0;
209 let mut t_max = self.max_distance;
210 slab(
211 self.origin.x,
212 self.dir_x,
213 self.inv_dir_x,
214 bounds.min_x,
215 bounds.max_x,
216 &mut t_min,
217 &mut t_max,
218 ) && slab(
219 self.origin.y,
220 self.dir_y,
221 self.inv_dir_y,
222 bounds.min_y,
223 bounds.max_y,
224 &mut t_min,
225 &mut t_max,
226 ) && slab(
227 self.origin.z,
228 self.dir_z,
229 self.inv_dir_z,
230 bounds.min_z,
231 bounds.max_z,
232 &mut t_min,
233 &mut t_max,
234 )
235 }
236
237 /// Entry parameter `t` where the ray segment first touches `bounds` (`0.0` if the
238 /// origin is inside), or `None` if the segment misses. Used by ordered closest-hit
239 /// traversal.
240 #[inline]
241 pub fn enter_t(self, bounds: Box3D) -> Option<f64> {
242 if self.max_distance < 0.0 || self.max_distance.is_nan() {
243 return None;
244 }
245 let mut t_min: f64 = 0.0;
246 let mut t_max = self.max_distance;
247 let hit = slab(
248 self.origin.x,
249 self.dir_x,
250 self.inv_dir_x,
251 bounds.min_x,
252 bounds.max_x,
253 &mut t_min,
254 &mut t_max,
255 ) && slab(
256 self.origin.y,
257 self.dir_y,
258 self.inv_dir_y,
259 bounds.min_y,
260 bounds.max_y,
261 &mut t_min,
262 &mut t_max,
263 ) && slab(
264 self.origin.z,
265 self.dir_z,
266 self.inv_dir_z,
267 bounds.min_z,
268 bounds.max_z,
269 &mut t_min,
270 &mut t_max,
271 );
272 hit.then_some(t_min)
273 }
274
275 /// The closest triangle in `triangles` hit by this ray segment, as a
276 /// [`TriangleHit`] (`index` into the slice and `t` in direction-length
277 /// units), or `None` if the ray misses them all, by the Moller-Trumbore test.
278 ///
279 /// Works over either record type: [`Triangle3D`](crate::Triangle3D) tests in
280 /// `f64`, [`Triangle3DF32`](crate::Triangle3DF32) in `f32` (8 at a time with
281 /// the `simd` feature). Pair this with a query that yields candidate triangles
282 /// — an [`Index3DView`](crate::Index3DView)'s
283 /// [`triangles`](crate::Index3DView::triangles) payload, or the items a
284 /// `raycast` / `search` returns — for exact ray-mesh intersection.
285 #[inline]
286 pub fn closest_triangle<T: Triangle3>(&self, triangles: &[T]) -> Option<TriangleHit> {
287 let o = [self.origin.x, self.origin.y, self.origin.z];
288 let d = [self.dir_x, self.dir_y, self.dir_z];
289 T::closest_hit(o, d, self.max_distance, triangles)
290 }
291}
292
293#[inline]
294fn slab(
295 origin: f64,
296 direction: f64,
297 inverse: f64,
298 min: f64,
299 max: f64,
300 t_min: &mut f64,
301 t_max: &mut f64,
302) -> bool {
303 if direction == 0.0 {
304 return origin >= min && origin <= max;
305 }
306
307 let mut near = (min - origin) * inverse;
308 let mut far = (max - origin) * inverse;
309 if near > far {
310 core::mem::swap(&mut near, &mut far);
311 }
312
313 *t_min = (*t_min).max(near);
314 *t_max = (*t_max).min(far);
315 *t_min <= *t_max
316}
317
318#[cfg(test)]
319mod triangle_tests {
320 use super::*;
321 use crate::{Point3D, Triangle3D};
322
323 fn ray(o: [f64; 3], d: [f64; 3]) -> Ray3D {
324 Ray3D::new(Point3D::new(o[0], o[1], o[2]), d[0], d[1], d[2], 100.0)
325 }
326
327 #[test]
328 fn closest_triangle_hits_nearest_and_misses() {
329 // Two triangles facing the ray at z = 2 and z = 5; a +z ray should pick
330 // the nearer one regardless of slice order.
331 let near = Triangle3D::new([0.0, 0.0, 2.0], [2.0, 0.0, 2.0], [0.0, 2.0, 2.0]);
332 let far = Triangle3D::new([0.0, 0.0, 5.0], [2.0, 0.0, 5.0], [0.0, 2.0, 5.0]);
333 let tris = [far, near];
334 let hit = ray([0.25, 0.25, 0.0], [0.0, 0.0, 1.0])
335 .closest_triangle(&tris)
336 .unwrap();
337 assert_eq!(hit.index, 1);
338 assert!((hit.t - 2.0).abs() < 1e-4, "t={}", hit.t);
339
340 // Pointing away, and an empty slice, both miss.
341 assert!(
342 ray([0.25, 0.25, 0.0], [0.0, 0.0, -1.0])
343 .closest_triangle(&tris)
344 .is_none()
345 );
346 assert!(
347 ray([0.25, 0.25, 0.0], [0.0, 0.0, 1.0])
348 .closest_triangle::<Triangle3D>(&[])
349 .is_none()
350 );
351 }
352
353 #[test]
354 fn closest_triangle_f32_records() {
355 use crate::Triangle3DF32;
356 let near = Triangle3DF32::new([0.0, 0.0, 2.0], [2.0, 0.0, 2.0], [0.0, 2.0, 2.0]);
357 let far = Triangle3DF32::new([0.0, 0.0, 5.0], [2.0, 0.0, 5.0], [0.0, 2.0, 5.0]);
358 let hit = ray([0.25, 0.25, 0.0], [0.0, 0.0, 1.0])
359 .closest_triangle(&[far, near])
360 .unwrap();
361 assert_eq!(hit.index, 1);
362 assert!((hit.t - 2.0).abs() < 1e-4, "t={}", hit.t);
363 }
364}