1 /**
2  * This is an implementation of the <a
3  * href="http://www.cg.cs.tu-bs.de/media/publications/fast-rayaxis-aligned-bounding-box-overlap-tests-using-ray-slopes.pdf">Fast Ray/Axis-Aligned Bounding Box
4  * Overlap Tests using Ray Slopes</a> paper.
5  * <p>
6  * It is an efficient implementation when testing many axis-aligned boxes against the same ray.
7  * <p>
8  * This class is thread-safe and can be used in a multithreaded environment when testing many axis-aligned boxes against the same ray concurrently.
9  * 
10  * @author Kai Burjack
11  */
12 module doml.ray_aabb_intersection;
13 
14 import Math = doml.math;
15 
16 import std.math.traits: isNaN;
17 
18 /*
19  * The MIT License
20  *
21  * Copyright (c) 2015-2021 Kai Burjack
22  %#%# Translated by jordan4ibanez
23  *
24  * Permission is hereby granted, free of charge, to any person obtaining a copy
25  * of this software and associated documentation files (the "Software"), to deal
26  * in the Software without restriction, including without limitation the rights
27  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
28  * copies of the Software, and to permit persons to whom the Software is
29  * furnished to do so, subject to the following conditions:
30  *
31  * The above copyright notice and this permission notice shall be included in
32  * all copies or substantial portions of the Software.
33  *
34  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
35  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
36  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
37  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
38  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
39  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
40  * THE SOFTWARE.
41  */
42 
43 /**
44  * This is an implementation of the <a
45  * href="http://www.cg.cs.tu-bs.de/media/publications/fast-rayaxis-aligned-bounding-box-overlap-tests-using-ray-slopes.pdf">Fast Ray/Axis-Aligned Bounding Box
46  * Overlap Tests using Ray Slopes</a> paper.
47  * <p>
48  * It is an efficient implementation when testing many axis-aligned boxes against the same ray.
49  * <p>
50  * This class is thread-safe and can be used in a multithreaded environment when testing many axis-aligned boxes against the same ray concurrently.
51  * 
52  * @author Kai Burjack
53  */
54 struct RayAabIntersection {
55     private double originX, originY, originZ;
56     private double dirX, dirY, dirZ;
57 
58     /* Needed for ray slope intersection method */
59     private double c_xy, c_yx, c_zy, c_yz, c_xz, c_zx;
60     private double s_xy, s_yx, s_zy, s_yz, s_xz, s_zx;
61     private byte classification;
62 
63     /**
64      * Create a new {@link RayAabIntersection} and initialize it with a ray with origin <code>(originX, originY, originZ)</code>
65      * and direction <code>(dirX, dirY, dirZ)</code>.
66      * <p>
67      * In order to change the direction and/or origin of the ray later, use {@link #set(double, double, double, double, double, double) set()}.
68      * 
69      * @see #set(double, double, double, double, double, double)
70      * 
71      * @param originX
72      *          the x coordinate of the origin
73      * @param originY
74      *          the y coordinate of the origin
75      * @param originZ
76      *          the z coordinate of the origin
77      * @param dirX
78      *          the x coordinate of the direction
79      * @param dirY
80      *          the y coordinate of the direction
81      * @param dirZ
82      *          the z coordinate of the direction
83      */
84     this(double originX, double originY, double originZ, double dirX, double dirY, double dirZ) {
85         set(originX, originY, originZ, dirX, dirY, dirZ);
86     }
87 
88     /**
89      * Update the ray stored by this {@link RayAabIntersection} with the new origin <code>(originX, originY, originZ)</code>
90      * and direction <code>(dirX, dirY, dirZ)</code>.
91      * 
92      * @param originX
93      *          the x coordinate of the ray origin
94      * @param originY
95      *          the y coordinate of the ray origin
96      * @param originZ
97      *          the z coordinate of the ray origin
98      * @param dirX
99      *          the x coordinate of the ray direction
100      * @param dirY
101      *          the y coordinate of the ray direction
102      * @param dirZ
103      *          the z coordinate of the ray direction
104      */
105     public void set(double originX, double originY, double originZ, double dirX, double dirY, double dirZ) {
106         this.originX = originX;
107         this.originY = originY;
108         this.originZ = originZ;
109         this.dirX = dirX;
110         this.dirY = dirY;
111         this.dirZ = dirZ;
112         precomputeSlope();
113     }
114     /**
115      * Precompute the values necessary for the ray slope algorithm.
116      */
117     private void precomputeSlope() {
118         double invDirX = 1.0f / dirX;
119         double invDirY = 1.0f / dirY;
120         double invDirZ = 1.0f / dirZ;
121         s_yx = dirX * invDirY;
122         s_xy = dirY * invDirX;
123         s_zy = dirY * invDirZ;
124         s_yz = dirZ * invDirY;
125         s_xz = dirZ * invDirX;
126         s_zx = dirX * invDirZ;
127         c_xy = originY - s_xy * originX;
128         c_yx = originX - s_yx * originY;
129         c_zy = originY - s_zy * originZ;
130         c_yz = originZ - s_yz * originY;
131         c_xz = originZ - s_xz * originX; // <- original paper had a bug here. It switched originZ/originX
132         c_zx = originX - s_zx * originZ; // <- original paper had a bug here. It switched originZ/originX
133         int sgnX = cast(int)Math.math_signum(dirX);
134         int sgnY = cast(int)Math.math_signum(dirY);
135         int sgnZ = cast(int)Math.math_signum(dirZ);
136         classification = cast(byte) ((sgnZ+1) << 4 | (sgnY+1) << 2 | (sgnX+1));
137     }
138 
139     /**
140      * Test whether the ray stored in this {@link RayAabIntersection} intersect the axis-aligned box
141      * given via its minimum corner <code>(minX, minY, minZ)</code> and its maximum corner <code>(maxX, maxY, maxZ)</code>.
142      * <p>
143      * This implementation uses a tableswitch to dispatch to the correct intersection method.
144      * <p>
145      * This method is thread-safe and can be used to test many axis-aligned boxes concurrently.
146      * 
147      * @param minX
148      *          the x coordinate of the minimum corner
149      * @param minY
150      *          the y coordinate of the minimum corner
151      * @param minZ
152      *          the z coordinate of the minimum corner
153      * @param maxX
154      *          the x coordinate of the maximum corner
155      * @param maxY
156      *          the y coordinate of the maximum corner
157      * @param maxZ
158      *          the z coordinate of the maximum corner
159      * @return <code>true</code> iff the ray intersects the given axis-aligned box; <code>false</code> otherwise
160      */
161     public bool test(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
162         // tableswitch with dense and consecutive cases (will be a simple jump based on the switch argument)
163         switch (classification) {
164         case 0: // 0b000000: // MMM
165             return MMM(minX, minY, minZ, maxX, maxY, maxZ);
166         case 1: // 0b000001: // OMM
167             return OMM(minX, minY, minZ, maxX, maxY, maxZ);
168         case 2: // 0b000010: // PMM
169             return PMM(minX, minY, minZ, maxX, maxY, maxZ);
170         case 3: // 0b000011: // not used
171             return false;
172         case 4: // 0b000100: // MOM 
173             return MOM(minX, minY, minZ, maxX, maxY, maxZ);
174         case 5: // 0b000101: // OOM
175             return OOM(minX, minY, minZ, maxX, maxY);
176         case 6: // 0b000110: // POM
177             return POM(minX, minY, minZ, maxX, maxY, maxZ);
178         case 7: // 0b000111: // not used
179             return false;
180         case 8: // 0b001000: // MPM
181             return MPM(minX, minY, minZ, maxX, maxY, maxZ);
182         case 9: // 0b001001: // OPM
183             return OPM(minX, minY, minZ, maxX, maxY, maxZ);
184         case 10: // 0b001010: // PPM
185             return PPM(minX, minY, minZ, maxX, maxY, maxZ);
186         case 11: // 0b001011: // not used
187         case 12: // 0b001100: // not used
188         case 13: // 0b001101: // not used
189         case 14: // 0b001110: // not used
190         case 15: // 0b001111: // not used
191             return false;
192         case 16: // 0b010000: // MMO
193             return MMO(minX, minY, minZ, maxX, maxY, maxZ);
194         case 17: // 0b010001: // OMO
195             return OMO(minX, minY, minZ, maxX, maxZ);
196         case 18: // 0b010010: // PMO
197             return PMO(minX, minY, minZ, maxX, maxY, maxZ);
198         case 19: // 0b010011: // not used
199             return false;
200         case 20: // 0b010100: // MOO
201             return MOO(minX, minY, minZ, maxY, maxZ);
202         case 21: // 0b010101: // OOO
203             return false; // <- degenerate case
204         case 22: // 0b010110: // POO
205             return POO(minY, minZ, maxX, maxY, maxZ);
206         case 23: // 0b010111: // not used
207             return false;
208         case 24: // 0b011000: // MPO
209             return MPO(minX, minY, minZ, maxX, maxY, maxZ);
210         case 25: // 0b011001: // OPO
211             return OPO(minX, minZ, maxX, maxY, maxZ);
212         case 26: // 0b011010: // PPO
213             return PPO(minX, minY, minZ, maxX, maxY, maxZ);
214         case 27: // 0b011011: // not used
215         case 28: // 0b011100: // not used
216         case 29: // 0b011101: // not used
217         case 30: // 0b011110: // not used
218         case 31: // 0b011111: // not used
219             return false;
220         case 32: // 0b100000: // MMP
221             return MMP(minX, minY, minZ, maxX, maxY, maxZ);
222         case 33: // 0b100001: // OMP
223             return OMP(minX, minY, minZ, maxX, maxY, maxZ);
224         case 34: // 0b100010: // PMP
225             return PMP(minX, minY, minZ, maxX, maxY, maxZ);
226         case 35: // 0b100011: // not used
227             return false;
228         case 36: // 0b100100: // MOP
229             return MOP(minX, minY, minZ, maxX, maxY, maxZ);
230         case 37: // 0b100101: // OOP
231             return OOP(minX, minY, maxX, maxY, maxZ);
232         case 38: // 0b100110: // POP
233             return POP(minX, minY, minZ, maxX, maxY, maxZ);
234         case 39: // 0b100111: // not used
235             return false;
236         case 40: // 0b101000: // MPP
237             return MPP(minX, minY, minZ, maxX, maxY, maxZ);
238         case 41: // 0b101001: // OPP
239             return OPP(minX, minY, minZ, maxX, maxY, maxZ);
240         case 42: // 0b101010: // PPP
241             return PPP(minX, minY, minZ, maxX, maxY, maxZ);
242         default:
243             return false;
244         }
245     }
246 
247     /* Intersection tests for all possible ray direction cases */
248 
249     private bool MMM(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
250         return originX >= minX && originY >= minY && originZ >= minZ
251             && s_xy * minX - maxY + c_xy <= 0.0f
252             && s_yx * minY - maxX + c_yx <= 0.0f
253             && s_zy * minZ - maxY + c_zy <= 0.0f
254             && s_yz * minY - maxZ + c_yz <= 0.0f
255             && s_xz * minX - maxZ + c_xz <= 0.0f
256             && s_zx * minZ - maxX + c_zx <= 0.0f;
257     }
258     private bool OMM(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
259         return originX >= minX && originX <= maxX && originY >= minY && originZ >= minZ
260             && s_zy * minZ - maxY + c_zy <= 0.0f
261             && s_yz * minY - maxZ + c_yz <= 0.0f;
262     }
263     private bool PMM(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
264         return originX <= maxX && originY >= minY && originZ >= minZ
265             && s_xy * maxX - maxY + c_xy <= 0.0f
266             && s_yx * minY - minX + c_yx >= 0.0f
267             && s_zy * minZ - maxY + c_zy <= 0.0f
268             && s_yz * minY - maxZ + c_yz <= 0.0f
269             && s_xz * maxX - maxZ + c_xz <= 0.0f
270             && s_zx * minZ - minX + c_zx >= 0.0f;
271     }
272     private bool MOM(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
273         return originY >= minY && originY <= maxY && originX >= minX && originZ >= minZ
274             && s_xz * minX - maxZ + c_xz <= 0.0f
275             && s_zx * minZ - maxX + c_zx <= 0.0f;
276     }
277     private bool OOM(double minX, double minY, double minZ, double maxX, double maxY) {
278         return originZ >= minZ && originX >= minX && originX <= maxX && originY >= minY && originY <= maxY;
279     }
280     private bool POM(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
281         return originY >= minY && originY <= maxY && originX <= maxX && originZ >= minZ
282             && s_xz * maxX - maxZ + c_xz <= 0.0f
283             && s_zx * minZ - minX + c_zx >= 0.0f;
284     }
285     private bool MPM(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
286         return originX >= minX && originY <= maxY && originZ >= minZ
287             && s_xy * minX - minY + c_xy >= 0.0f
288             && s_yx * maxY - maxX + c_yx <= 0.0f
289             && s_zy * minZ - minY + c_zy >= 0.0f
290             && s_yz * maxY - maxZ + c_yz <= 0.0f
291             && s_xz * minX - maxZ + c_xz <= 0.0f
292             && s_zx * minZ - maxX + c_zx <= 0.0f;
293     }
294     private bool OPM(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
295         return originX >= minX && originX <= maxX && originY <= maxY && originZ >= minZ
296             && s_zy * minZ - minY + c_zy >= 0.0f
297             && s_yz * maxY - maxZ + c_yz <= 0.0f;
298     }
299     private bool PPM(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
300         return originX <= maxX && originY <= maxY && originZ >= minZ
301             && s_xy * maxX - minY + c_xy >= 0.0f
302             && s_yx * maxY - minX + c_yx >= 0.0f
303             && s_zy * minZ - minY + c_zy >= 0.0f
304             && s_yz * maxY - maxZ + c_yz <= 0.0f
305             && s_xz * maxX - maxZ + c_xz <= 0.0f
306             && s_zx * minZ - minX + c_zx >= 0.0f;
307     }
308     private bool MMO(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
309         return originZ >= minZ && originZ <= maxZ && originX >= minX && originY >= minY
310             && s_xy * minX - maxY + c_xy <= 0.0f
311             && s_yx * minY - maxX + c_yx <= 0.0f;
312     }
313     private bool OMO(double minX, double minY, double minZ, double maxX, double maxZ) {
314         return originY >= minY && originX >= minX && originX <= maxX && originZ >= minZ && originZ <= maxZ;
315     }
316     private bool PMO(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
317         return originZ >= minZ && originZ <= maxZ && originX <= maxX && originY >= minY
318             && s_xy * maxX - maxY + c_xy <= 0.0f
319             && s_yx * minY - minX + c_yx >= 0.0f;
320     }
321     private bool MOO(double minX, double minY, double minZ, double maxY, double maxZ) {
322         return originX >= minX && originY >= minY && originY <= maxY && originZ >= minZ && originZ <= maxZ;
323     }
324     private bool POO(double minY, double minZ, double maxX, double maxY, double maxZ) {
325         return originX <= maxX && originY >= minY && originY <= maxY && originZ >= minZ && originZ <= maxZ;
326     }
327     private bool MPO(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
328         return originZ >= minZ && originZ <= maxZ && originX >= minX && originY <= maxY
329             && s_xy * minX - minY + c_xy >= 0.0f
330             && s_yx * maxY - maxX + c_yx <= 0.0f;
331     }
332     private bool OPO(double minX, double minZ, double maxX, double maxY, double maxZ) {
333         return originY <= maxY && originX >= minX && originX <= maxX && originZ >= minZ && originZ <= maxZ;
334     }
335     private bool PPO(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
336         return originZ >= minZ && originZ <= maxZ && originX <= maxX && originY <= maxY
337             && s_xy * maxX - minY + c_xy >= 0.0f
338             && s_yx * maxY - minX + c_yx >= 0.0f;
339     }
340     private bool MMP(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
341         return originX >= minX && originY >= minY && originZ <= maxZ
342             && s_xy * minX - maxY + c_xy <= 0.0f
343             && s_yx * minY - maxX + c_yx <= 0.0f
344             && s_zy * maxZ - maxY + c_zy <= 0.0f
345             && s_yz * minY - minZ + c_yz >= 0.0f
346             && s_xz * minX - minZ + c_xz >= 0.0f
347             && s_zx * maxZ - maxX + c_zx <= 0.0f;
348     }
349     private bool OMP(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
350         return originX >= minX && originX <= maxX && originY >= minY && originZ <= maxZ
351             && s_zy * maxZ - maxY + c_zy <= 0.0f
352             && s_yz * minY - minZ + c_yz >= 0.0f;
353     }
354     private bool PMP(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
355         return originX <= maxX && originY >= minY && originZ <= maxZ
356             && s_xy * maxX - maxY + c_xy <= 0.0f
357             && s_yx * minY - minX + c_yx >= 0.0f
358             && s_zy * maxZ - maxY + c_zy <= 0.0f
359             && s_yz * minY - minZ + c_yz >= 0.0f
360             && s_xz * maxX - minZ + c_xz >= 0.0f
361             && s_zx * maxZ - minX + c_zx >= 0.0f;
362     }
363     private bool MOP(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
364         return originY >= minY && originY <= maxY && originX >= minX && originZ <= maxZ
365             && s_xz * minX - minZ + c_xz >= 0.0f
366             && s_zx * maxZ - maxX + c_zx <= 0.0f;
367     }
368     private bool OOP(double minX, double minY, double maxX, double maxY, double maxZ) {
369         return originZ <= maxZ && originX >= minX && originX <= maxX && originY >= minY && originY <= maxY;
370     }
371     private bool POP(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
372         return originY >= minY && originY <= maxY && originX <= maxX && originZ <= maxZ
373             && s_xz * maxX - minZ + c_xz >= 0.0f
374             && s_zx * maxZ - minX + c_zx <= 0.0f;
375     }
376     private bool MPP(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
377         return originX >= minX && originY <= maxY && originZ <= maxZ
378             && s_xy * minX - minY + c_xy >= 0.0f
379             && s_yx * maxY - maxX + c_yx <= 0.0f
380             && s_zy * maxZ - minY + c_zy >= 0.0f
381             && s_yz * maxY - minZ + c_yz >= 0.0f
382             && s_xz * minX - minZ + c_xz >= 0.0f
383             && s_zx * maxZ - maxX + c_zx <= 0.0f;
384     }
385     private bool OPP(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
386         return originX >= minX && originX <= maxX && originY <= maxY && originZ <= maxZ
387             && s_zy * maxZ - minY + c_zy <= 0.0f
388             && s_yz * maxY - minZ + c_yz <= 0.0f;
389     }
390     private bool PPP(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
391         return originX <= maxX && originY <= maxY && originZ <= maxZ
392             && s_xy * maxX - minY + c_xy >= 0.0f
393             && s_yx * maxY - minX + c_yx >= 0.0f
394             && s_zy * maxZ - minY + c_zy >= 0.0f
395             && s_yz * maxY - minZ + c_yz >= 0.0f
396             && s_xz * maxX - minZ + c_xz >= 0.0f
397             && s_zx * maxZ - minX + c_zx >= 0.0f;
398     }
399 }