1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2002-2012, Industrial Light & Magic, a division of Lucas
4 // Digital Ltd. LLC
5 //
6 // All rights reserved.
7 //
8 // Redistribution and use in source and binary forms, with or without
9 // modification, are permitted provided that the following conditions are
10 // met:
11 // * Redistributions of source code must retain the above copyright
12 // notice, this list of conditions and the following disclaimer.
13 // * Redistributions in binary form must reproduce the above
14 // copyright notice, this list of conditions and the following disclaimer
15 // in the documentation and/or other materials provided with the
16 // distribution.
17 // * Neither the name of Industrial Light & Magic nor the names of
18 // its contributors may be used to endorse or promote products derived
19 // from this software without specific prior written permission.
20 //
21 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 //
33 ///////////////////////////////////////////////////////////////////////////
34
35
36
37 #ifndef INCLUDED_IMATHBOXALGO_H
38 #define INCLUDED_IMATHBOXALGO_H
39
40
41 //---------------------------------------------------------------------------
42 //
43 // This file contains algorithms applied to or in conjunction
44 // with bounding boxes (Imath::Box). These algorithms require
45 // more headers to compile. The assumption made is that these
46 // functions are called much less often than the basic box
47 // functions or these functions require more support classes.
48 //
49 // Contains:
50 //
51 // T clip<T>(const T& in, const Box<T>& box)
52 //
53 // Vec3<T> closestPointOnBox(const Vec3<T>&, const Box<Vec3<T>>& )
54 //
55 // Vec3<T> closestPointInBox(const Vec3<T>&, const Box<Vec3<T>>& )
56 //
57 // Box< Vec3<S> > transform(const Box<Vec3<S>>&, const Matrix44<T>&)
58 // Box< Vec3<S> > affineTransform(const Box<Vec3<S>>&, const Matrix44<T>&)
59 //
60 // void transform(const Box<Vec3<S>>&, const Matrix44<T>&, Box<V3ec3<S>>&)
61 // void affineTransform(const Box<Vec3<S>>&,
62 // const Matrix44<T>&,
63 // Box<V3ec3<S>>&)
64 //
65 // bool findEntryAndExitPoints(const Line<T> &line,
66 // const Box< Vec3<T> > &box,
67 // Vec3<T> &enterPoint,
68 // Vec3<T> &exitPoint)
69 //
70 // bool intersects(const Box<Vec3<T>> &box,
71 // const Line3<T> &ray,
72 // Vec3<T> intersectionPoint)
73 //
74 // bool intersects(const Box<Vec3<T>> &box, const Line3<T> &ray)
75 //
76 //---------------------------------------------------------------------------
77
78 #include "ImathBox.h"
79 #include "ImathMatrix.h"
80 #include "ImathLineAlgo.h"
81 #include "ImathPlane.h"
82 #include "ImathNamespace.h"
83
84 IMATH_INTERNAL_NAMESPACE_HEADER_ENTER
85
86
87 template <class T>
88 inline T
clip(const T & p,const Box<T> & box)89 clip (const T &p, const Box<T> &box)
90 {
91 //
92 // Clip the coordinates of a point, p, against a box.
93 // The result, q, is the closest point to p that is inside the box.
94 //
95
96 T q;
97
98 for (int i = 0; i < int (box.min.dimensions()); i++)
99 {
100 if (p[i] < box.min[i])
101 q[i] = box.min[i];
102 else if (p[i] > box.max[i])
103 q[i] = box.max[i];
104 else
105 q[i] = p[i];
106 }
107
108 return q;
109 }
110
111
112 template <class T>
113 inline T
closestPointInBox(const T & p,const Box<T> & box)114 closestPointInBox (const T &p, const Box<T> &box)
115 {
116 return clip (p, box);
117 }
118
119
120 template <class T>
121 Vec3<T>
closestPointOnBox(const Vec3<T> & p,const Box<Vec3<T>> & box)122 closestPointOnBox (const Vec3<T> &p, const Box< Vec3<T> > &box)
123 {
124 //
125 // Find the point, q, on the surface of
126 // the box, that is closest to point p.
127 //
128 // If the box is empty, return p.
129 //
130
131 if (box.isEmpty())
132 return p;
133
134 Vec3<T> q = closestPointInBox (p, box);
135
136 if (q == p)
137 {
138 Vec3<T> d1 = p - box.min;
139 Vec3<T> d2 = box.max - p;
140
141 Vec3<T> d ((d1.x < d2.x)? d1.x: d2.x,
142 (d1.y < d2.y)? d1.y: d2.y,
143 (d1.z < d2.z)? d1.z: d2.z);
144
145 if (d.x < d.y && d.x < d.z)
146 {
147 q.x = (d1.x < d2.x)? box.min.x: box.max.x;
148 }
149 else if (d.y < d.z)
150 {
151 q.y = (d1.y < d2.y)? box.min.y: box.max.y;
152 }
153 else
154 {
155 q.z = (d1.z < d2.z)? box.min.z: box.max.z;
156 }
157 }
158
159 return q;
160 }
161
162
163 template <class S, class T>
164 Box< Vec3<S> >
transform(const Box<Vec3<S>> & box,const Matrix44<T> & m)165 transform (const Box< Vec3<S> > &box, const Matrix44<T> &m)
166 {
167 //
168 // Transform a 3D box by a matrix, and compute a new box that
169 // tightly encloses the transformed box.
170 //
171 // If m is an affine transform, then we use James Arvo's fast
172 // method as described in "Graphics Gems", Academic Press, 1990,
173 // pp. 548-550.
174 //
175
176 //
177 // A transformed empty box is still empty, and a transformed infinite box
178 // is still infinite
179 //
180
181 if (box.isEmpty() || box.isInfinite())
182 return box;
183
184 //
185 // If the last column of m is (0 0 0 1) then m is an affine
186 // transform, and we use the fast Graphics Gems trick.
187 //
188
189 if (m[0][3] == 0 && m[1][3] == 0 && m[2][3] == 0 && m[3][3] == 1)
190 {
191 Box< Vec3<S> > newBox;
192
193 for (int i = 0; i < 3; i++)
194 {
195 newBox.min[i] = newBox.max[i] = (S) m[3][i];
196
197 for (int j = 0; j < 3; j++)
198 {
199 S a, b;
200
201 a = (S) m[j][i] * box.min[j];
202 b = (S) m[j][i] * box.max[j];
203
204 if (a < b)
205 {
206 newBox.min[i] += a;
207 newBox.max[i] += b;
208 }
209 else
210 {
211 newBox.min[i] += b;
212 newBox.max[i] += a;
213 }
214 }
215 }
216
217 return newBox;
218 }
219
220 //
221 // M is a projection matrix. Do things the naive way:
222 // Transform the eight corners of the box, and find an
223 // axis-parallel box that encloses the transformed corners.
224 //
225
226 Vec3<S> points[8];
227
228 points[0][0] = points[1][0] = points[2][0] = points[3][0] = box.min[0];
229 points[4][0] = points[5][0] = points[6][0] = points[7][0] = box.max[0];
230
231 points[0][1] = points[1][1] = points[4][1] = points[5][1] = box.min[1];
232 points[2][1] = points[3][1] = points[6][1] = points[7][1] = box.max[1];
233
234 points[0][2] = points[2][2] = points[4][2] = points[6][2] = box.min[2];
235 points[1][2] = points[3][2] = points[5][2] = points[7][2] = box.max[2];
236
237 Box< Vec3<S> > newBox;
238
239 for (int i = 0; i < 8; i++)
240 newBox.extendBy (points[i] * m);
241
242 return newBox;
243 }
244
245 template <class S, class T>
246 void
transform(const Box<Vec3<S>> & box,const Matrix44<T> & m,Box<Vec3<S>> & result)247 transform (const Box< Vec3<S> > &box,
248 const Matrix44<T> &m,
249 Box< Vec3<S> > &result)
250 {
251 //
252 // Transform a 3D box by a matrix, and compute a new box that
253 // tightly encloses the transformed box.
254 //
255 // If m is an affine transform, then we use James Arvo's fast
256 // method as described in "Graphics Gems", Academic Press, 1990,
257 // pp. 548-550.
258 //
259
260 //
261 // A transformed empty box is still empty, and a transformed infinite
262 // box is still infinite
263 //
264
265 if (box.isEmpty() || box.isInfinite())
266 {
267 return;
268 }
269
270 //
271 // If the last column of m is (0 0 0 1) then m is an affine
272 // transform, and we use the fast Graphics Gems trick.
273 //
274
275 if (m[0][3] == 0 && m[1][3] == 0 && m[2][3] == 0 && m[3][3] == 1)
276 {
277 for (int i = 0; i < 3; i++)
278 {
279 result.min[i] = result.max[i] = (S) m[3][i];
280
281 for (int j = 0; j < 3; j++)
282 {
283 S a, b;
284
285 a = (S) m[j][i] * box.min[j];
286 b = (S) m[j][i] * box.max[j];
287
288 if (a < b)
289 {
290 result.min[i] += a;
291 result.max[i] += b;
292 }
293 else
294 {
295 result.min[i] += b;
296 result.max[i] += a;
297 }
298 }
299 }
300
301 return;
302 }
303
304 //
305 // M is a projection matrix. Do things the naive way:
306 // Transform the eight corners of the box, and find an
307 // axis-parallel box that encloses the transformed corners.
308 //
309
310 Vec3<S> points[8];
311
312 points[0][0] = points[1][0] = points[2][0] = points[3][0] = box.min[0];
313 points[4][0] = points[5][0] = points[6][0] = points[7][0] = box.max[0];
314
315 points[0][1] = points[1][1] = points[4][1] = points[5][1] = box.min[1];
316 points[2][1] = points[3][1] = points[6][1] = points[7][1] = box.max[1];
317
318 points[0][2] = points[2][2] = points[4][2] = points[6][2] = box.min[2];
319 points[1][2] = points[3][2] = points[5][2] = points[7][2] = box.max[2];
320
321 for (int i = 0; i < 8; i++)
322 result.extendBy (points[i] * m);
323 }
324
325
326 template <class S, class T>
327 Box< Vec3<S> >
affineTransform(const Box<Vec3<S>> & box,const Matrix44<T> & m)328 affineTransform (const Box< Vec3<S> > &box, const Matrix44<T> &m)
329 {
330 //
331 // Transform a 3D box by a matrix whose rightmost column
332 // is (0 0 0 1), and compute a new box that tightly encloses
333 // the transformed box.
334 //
335 // As in the transform() function, above, we use James Arvo's
336 // fast method.
337 //
338
339 if (box.isEmpty() || box.isInfinite())
340 {
341 //
342 // A transformed empty or infinite box is still empty or infinite
343 //
344
345 return box;
346 }
347
348 Box< Vec3<S> > newBox;
349
350 for (int i = 0; i < 3; i++)
351 {
352 newBox.min[i] = newBox.max[i] = (S) m[3][i];
353
354 for (int j = 0; j < 3; j++)
355 {
356 S a, b;
357
358 a = (S) m[j][i] * box.min[j];
359 b = (S) m[j][i] * box.max[j];
360
361 if (a < b)
362 {
363 newBox.min[i] += a;
364 newBox.max[i] += b;
365 }
366 else
367 {
368 newBox.min[i] += b;
369 newBox.max[i] += a;
370 }
371 }
372 }
373
374 return newBox;
375 }
376
377 template <class S, class T>
378 void
affineTransform(const Box<Vec3<S>> & box,const Matrix44<T> & m,Box<Vec3<S>> & result)379 affineTransform (const Box< Vec3<S> > &box,
380 const Matrix44<T> &m,
381 Box<Vec3<S> > &result)
382 {
383 //
384 // Transform a 3D box by a matrix whose rightmost column
385 // is (0 0 0 1), and compute a new box that tightly encloses
386 // the transformed box.
387 //
388 // As in the transform() function, above, we use James Arvo's
389 // fast method.
390 //
391
392 if (box.isEmpty())
393 {
394 //
395 // A transformed empty box is still empty
396 //
397 result.makeEmpty();
398 return;
399 }
400
401 if (box.isInfinite())
402 {
403 //
404 // A transformed infinite box is still infinite
405 //
406 result.makeInfinite();
407 return;
408 }
409
410 for (int i = 0; i < 3; i++)
411 {
412 result.min[i] = result.max[i] = (S) m[3][i];
413
414 for (int j = 0; j < 3; j++)
415 {
416 S a, b;
417
418 a = (S) m[j][i] * box.min[j];
419 b = (S) m[j][i] * box.max[j];
420
421 if (a < b)
422 {
423 result.min[i] += a;
424 result.max[i] += b;
425 }
426 else
427 {
428 result.min[i] += b;
429 result.max[i] += a;
430 }
431 }
432 }
433 }
434
435
436 template <class T>
437 bool
findEntryAndExitPoints(const Line3<T> & r,const Box<Vec3<T>> & b,Vec3<T> & entry,Vec3<T> & exit)438 findEntryAndExitPoints (const Line3<T> &r,
439 const Box<Vec3<T> > &b,
440 Vec3<T> &entry,
441 Vec3<T> &exit)
442 {
443 //
444 // Compute the points where a ray, r, enters and exits a box, b:
445 //
446 // findEntryAndExitPoints() returns
447 //
448 // - true if the ray starts inside the box or if the
449 // ray starts outside and intersects the box
450 //
451 // - false otherwise (that is, if the ray does not
452 // intersect the box)
453 //
454 // The entry and exit points are
455 //
456 // - points on two of the faces of the box when
457 // findEntryAndExitPoints() returns true
458 // (The entry end exit points may be on either
459 // side of the ray's origin)
460 //
461 // - undefined when findEntryAndExitPoints()
462 // returns false
463 //
464
465 if (b.isEmpty())
466 {
467 //
468 // No ray intersects an empty box
469 //
470
471 return false;
472 }
473
474 //
475 // The following description assumes that the ray's origin is outside
476 // the box, but the code below works even if the origin is inside the
477 // box:
478 //
479 // Between one and three "frontfacing" sides of the box are oriented
480 // towards the ray's origin, and between one and three "backfacing"
481 // sides are oriented away from the ray's origin.
482 // We intersect the ray with the planes that contain the sides of the
483 // box, and compare the distances between the ray's origin and the
484 // ray-plane intersections. The ray intersects the box if the most
485 // distant frontfacing intersection is nearer than the nearest
486 // backfacing intersection. If the ray does intersect the box, then
487 // the most distant frontfacing ray-plane intersection is the entry
488 // point and the nearest backfacing ray-plane intersection is the
489 // exit point.
490 //
491
492 const T TMAX = limits<T>::max();
493
494 T tFrontMax = -TMAX;
495 T tBackMin = TMAX;
496
497 //
498 // Minimum and maximum X sides.
499 //
500
501 if (r.dir.x >= 0)
502 {
503 T d1 = b.max.x - r.pos.x;
504 T d2 = b.min.x - r.pos.x;
505
506 if (r.dir.x > 1 ||
507 (abs (d1) < TMAX * r.dir.x &&
508 abs (d2) < TMAX * r.dir.x))
509 {
510 T t1 = d1 / r.dir.x;
511 T t2 = d2 / r.dir.x;
512
513 if (tBackMin > t1)
514 {
515 tBackMin = t1;
516
517 exit.x = b.max.x;
518 exit.y = clamp (r.pos.y + t1 * r.dir.y, b.min.y, b.max.y);
519 exit.z = clamp (r.pos.z + t1 * r.dir.z, b.min.z, b.max.z);
520 }
521
522 if (tFrontMax < t2)
523 {
524 tFrontMax = t2;
525
526 entry.x = b.min.x;
527 entry.y = clamp (r.pos.y + t2 * r.dir.y, b.min.y, b.max.y);
528 entry.z = clamp (r.pos.z + t2 * r.dir.z, b.min.z, b.max.z);
529 }
530 }
531 else if (r.pos.x < b.min.x || r.pos.x > b.max.x)
532 {
533 return false;
534 }
535 }
536 else // r.dir.x < 0
537 {
538 T d1 = b.min.x - r.pos.x;
539 T d2 = b.max.x - r.pos.x;
540
541 if (r.dir.x < -1 ||
542 (abs (d1) < -TMAX * r.dir.x &&
543 abs (d2) < -TMAX * r.dir.x))
544 {
545 T t1 = d1 / r.dir.x;
546 T t2 = d2 / r.dir.x;
547
548 if (tBackMin > t1)
549 {
550 tBackMin = t1;
551
552 exit.x = b.min.x;
553 exit.y = clamp (r.pos.y + t1 * r.dir.y, b.min.y, b.max.y);
554 exit.z = clamp (r.pos.z + t1 * r.dir.z, b.min.z, b.max.z);
555 }
556
557 if (tFrontMax < t2)
558 {
559 tFrontMax = t2;
560
561 entry.x = b.max.x;
562 entry.y = clamp (r.pos.y + t2 * r.dir.y, b.min.y, b.max.y);
563 entry.z = clamp (r.pos.z + t2 * r.dir.z, b.min.z, b.max.z);
564 }
565 }
566 else if (r.pos.x < b.min.x || r.pos.x > b.max.x)
567 {
568 return false;
569 }
570 }
571
572 //
573 // Minimum and maximum Y sides.
574 //
575
576 if (r.dir.y >= 0)
577 {
578 T d1 = b.max.y - r.pos.y;
579 T d2 = b.min.y - r.pos.y;
580
581 if (r.dir.y > 1 ||
582 (abs (d1) < TMAX * r.dir.y &&
583 abs (d2) < TMAX * r.dir.y))
584 {
585 T t1 = d1 / r.dir.y;
586 T t2 = d2 / r.dir.y;
587
588 if (tBackMin > t1)
589 {
590 tBackMin = t1;
591
592 exit.x = clamp (r.pos.x + t1 * r.dir.x, b.min.x, b.max.x);
593 exit.y = b.max.y;
594 exit.z = clamp (r.pos.z + t1 * r.dir.z, b.min.z, b.max.z);
595 }
596
597 if (tFrontMax < t2)
598 {
599 tFrontMax = t2;
600
601 entry.x = clamp (r.pos.x + t2 * r.dir.x, b.min.x, b.max.x);
602 entry.y = b.min.y;
603 entry.z = clamp (r.pos.z + t2 * r.dir.z, b.min.z, b.max.z);
604 }
605 }
606 else if (r.pos.y < b.min.y || r.pos.y > b.max.y)
607 {
608 return false;
609 }
610 }
611 else // r.dir.y < 0
612 {
613 T d1 = b.min.y - r.pos.y;
614 T d2 = b.max.y - r.pos.y;
615
616 if (r.dir.y < -1 ||
617 (abs (d1) < -TMAX * r.dir.y &&
618 abs (d2) < -TMAX * r.dir.y))
619 {
620 T t1 = d1 / r.dir.y;
621 T t2 = d2 / r.dir.y;
622
623 if (tBackMin > t1)
624 {
625 tBackMin = t1;
626
627 exit.x = clamp (r.pos.x + t1 * r.dir.x, b.min.x, b.max.x);
628 exit.y = b.min.y;
629 exit.z = clamp (r.pos.z + t1 * r.dir.z, b.min.z, b.max.z);
630 }
631
632 if (tFrontMax < t2)
633 {
634 tFrontMax = t2;
635
636 entry.x = clamp (r.pos.x + t2 * r.dir.x, b.min.x, b.max.x);
637 entry.y = b.max.y;
638 entry.z = clamp (r.pos.z + t2 * r.dir.z, b.min.z, b.max.z);
639 }
640 }
641 else if (r.pos.y < b.min.y || r.pos.y > b.max.y)
642 {
643 return false;
644 }
645 }
646
647 //
648 // Minimum and maximum Z sides.
649 //
650
651 if (r.dir.z >= 0)
652 {
653 T d1 = b.max.z - r.pos.z;
654 T d2 = b.min.z - r.pos.z;
655
656 if (r.dir.z > 1 ||
657 (abs (d1) < TMAX * r.dir.z &&
658 abs (d2) < TMAX * r.dir.z))
659 {
660 T t1 = d1 / r.dir.z;
661 T t2 = d2 / r.dir.z;
662
663 if (tBackMin > t1)
664 {
665 tBackMin = t1;
666
667 exit.x = clamp (r.pos.x + t1 * r.dir.x, b.min.x, b.max.x);
668 exit.y = clamp (r.pos.y + t1 * r.dir.y, b.min.y, b.max.y);
669 exit.z = b.max.z;
670 }
671
672 if (tFrontMax < t2)
673 {
674 tFrontMax = t2;
675
676 entry.x = clamp (r.pos.x + t2 * r.dir.x, b.min.x, b.max.x);
677 entry.y = clamp (r.pos.y + t2 * r.dir.y, b.min.y, b.max.y);
678 entry.z = b.min.z;
679 }
680 }
681 else if (r.pos.z < b.min.z || r.pos.z > b.max.z)
682 {
683 return false;
684 }
685 }
686 else // r.dir.z < 0
687 {
688 T d1 = b.min.z - r.pos.z;
689 T d2 = b.max.z - r.pos.z;
690
691 if (r.dir.z < -1 ||
692 (abs (d1) < -TMAX * r.dir.z &&
693 abs (d2) < -TMAX * r.dir.z))
694 {
695 T t1 = d1 / r.dir.z;
696 T t2 = d2 / r.dir.z;
697
698 if (tBackMin > t1)
699 {
700 tBackMin = t1;
701
702 exit.x = clamp (r.pos.x + t1 * r.dir.x, b.min.x, b.max.x);
703 exit.y = clamp (r.pos.y + t1 * r.dir.y, b.min.y, b.max.y);
704 exit.z = b.min.z;
705 }
706
707 if (tFrontMax < t2)
708 {
709 tFrontMax = t2;
710
711 entry.x = clamp (r.pos.x + t2 * r.dir.x, b.min.x, b.max.x);
712 entry.y = clamp (r.pos.y + t2 * r.dir.y, b.min.y, b.max.y);
713 entry.z = b.max.z;
714 }
715 }
716 else if (r.pos.z < b.min.z || r.pos.z > b.max.z)
717 {
718 return false;
719 }
720 }
721
722 return tFrontMax <= tBackMin;
723 }
724
725
726 template<class T>
727 bool
intersects(const Box<Vec3<T>> & b,const Line3<T> & r,Vec3<T> & ip)728 intersects (const Box< Vec3<T> > &b, const Line3<T> &r, Vec3<T> &ip)
729 {
730 //
731 // Intersect a ray, r, with a box, b, and compute the intersection
732 // point, ip:
733 //
734 // intersect() returns
735 //
736 // - true if the ray starts inside the box or if the
737 // ray starts outside and intersects the box
738 //
739 // - false if the ray starts outside the box and intersects it,
740 // but the intersection is behind the ray's origin.
741 //
742 // - false if the ray starts outside and does not intersect it
743 //
744 // The intersection point is
745 //
746 // - the ray's origin if the ray starts inside the box
747 //
748 // - a point on one of the faces of the box if the ray
749 // starts outside the box
750 //
751 // - undefined when intersect() returns false
752 //
753
754 if (b.isEmpty())
755 {
756 //
757 // No ray intersects an empty box
758 //
759
760 return false;
761 }
762
763 if (b.intersects (r.pos))
764 {
765 //
766 // The ray starts inside the box
767 //
768
769 ip = r.pos;
770 return true;
771 }
772
773 //
774 // The ray starts outside the box. Between one and three "frontfacing"
775 // sides of the box are oriented towards the ray, and between one and
776 // three "backfacing" sides are oriented away from the ray.
777 // We intersect the ray with the planes that contain the sides of the
778 // box, and compare the distances between ray's origin and the ray-plane
779 // intersections.
780 // The ray intersects the box if the most distant frontfacing intersection
781 // is nearer than the nearest backfacing intersection. If the ray does
782 // intersect the box, then the most distant frontfacing ray-plane
783 // intersection is the ray-box intersection.
784 //
785
786 const T TMAX = limits<T>::max();
787
788 T tFrontMax = -1;
789 T tBackMin = TMAX;
790
791 //
792 // Minimum and maximum X sides.
793 //
794
795 if (r.dir.x > 0)
796 {
797 if (r.pos.x > b.max.x)
798 return false;
799
800 T d = b.max.x - r.pos.x;
801
802 if (r.dir.x > 1 || d < TMAX * r.dir.x)
803 {
804 T t = d / r.dir.x;
805
806 if (tBackMin > t)
807 tBackMin = t;
808 }
809
810 if (r.pos.x <= b.min.x)
811 {
812 T d = b.min.x - r.pos.x;
813 T t = (r.dir.x > 1 || d < TMAX * r.dir.x)? d / r.dir.x: TMAX;
814
815 if (tFrontMax < t)
816 {
817 tFrontMax = t;
818
819 ip.x = b.min.x;
820 ip.y = clamp (r.pos.y + t * r.dir.y, b.min.y, b.max.y);
821 ip.z = clamp (r.pos.z + t * r.dir.z, b.min.z, b.max.z);
822 }
823 }
824 }
825 else if (r.dir.x < 0)
826 {
827 if (r.pos.x < b.min.x)
828 return false;
829
830 T d = b.min.x - r.pos.x;
831
832 if (r.dir.x < -1 || d > TMAX * r.dir.x)
833 {
834 T t = d / r.dir.x;
835
836 if (tBackMin > t)
837 tBackMin = t;
838 }
839
840 if (r.pos.x >= b.max.x)
841 {
842 T d = b.max.x - r.pos.x;
843 T t = (r.dir.x < -1 || d > TMAX * r.dir.x)? d / r.dir.x: TMAX;
844
845 if (tFrontMax < t)
846 {
847 tFrontMax = t;
848
849 ip.x = b.max.x;
850 ip.y = clamp (r.pos.y + t * r.dir.y, b.min.y, b.max.y);
851 ip.z = clamp (r.pos.z + t * r.dir.z, b.min.z, b.max.z);
852 }
853 }
854 }
855 else // r.dir.x == 0
856 {
857 if (r.pos.x < b.min.x || r.pos.x > b.max.x)
858 return false;
859 }
860
861 //
862 // Minimum and maximum Y sides.
863 //
864
865 if (r.dir.y > 0)
866 {
867 if (r.pos.y > b.max.y)
868 return false;
869
870 T d = b.max.y - r.pos.y;
871
872 if (r.dir.y > 1 || d < TMAX * r.dir.y)
873 {
874 T t = d / r.dir.y;
875
876 if (tBackMin > t)
877 tBackMin = t;
878 }
879
880 if (r.pos.y <= b.min.y)
881 {
882 T d = b.min.y - r.pos.y;
883 T t = (r.dir.y > 1 || d < TMAX * r.dir.y)? d / r.dir.y: TMAX;
884
885 if (tFrontMax < t)
886 {
887 tFrontMax = t;
888
889 ip.x = clamp (r.pos.x + t * r.dir.x, b.min.x, b.max.x);
890 ip.y = b.min.y;
891 ip.z = clamp (r.pos.z + t * r.dir.z, b.min.z, b.max.z);
892 }
893 }
894 }
895 else if (r.dir.y < 0)
896 {
897 if (r.pos.y < b.min.y)
898 return false;
899
900 T d = b.min.y - r.pos.y;
901
902 if (r.dir.y < -1 || d > TMAX * r.dir.y)
903 {
904 T t = d / r.dir.y;
905
906 if (tBackMin > t)
907 tBackMin = t;
908 }
909
910 if (r.pos.y >= b.max.y)
911 {
912 T d = b.max.y - r.pos.y;
913 T t = (r.dir.y < -1 || d > TMAX * r.dir.y)? d / r.dir.y: TMAX;
914
915 if (tFrontMax < t)
916 {
917 tFrontMax = t;
918
919 ip.x = clamp (r.pos.x + t * r.dir.x, b.min.x, b.max.x);
920 ip.y = b.max.y;
921 ip.z = clamp (r.pos.z + t * r.dir.z, b.min.z, b.max.z);
922 }
923 }
924 }
925 else // r.dir.y == 0
926 {
927 if (r.pos.y < b.min.y || r.pos.y > b.max.y)
928 return false;
929 }
930
931 //
932 // Minimum and maximum Z sides.
933 //
934
935 if (r.dir.z > 0)
936 {
937 if (r.pos.z > b.max.z)
938 return false;
939
940 T d = b.max.z - r.pos.z;
941
942 if (r.dir.z > 1 || d < TMAX * r.dir.z)
943 {
944 T t = d / r.dir.z;
945
946 if (tBackMin > t)
947 tBackMin = t;
948 }
949
950 if (r.pos.z <= b.min.z)
951 {
952 T d = b.min.z - r.pos.z;
953 T t = (r.dir.z > 1 || d < TMAX * r.dir.z)? d / r.dir.z: TMAX;
954
955 if (tFrontMax < t)
956 {
957 tFrontMax = t;
958
959 ip.x = clamp (r.pos.x + t * r.dir.x, b.min.x, b.max.x);
960 ip.y = clamp (r.pos.y + t * r.dir.y, b.min.y, b.max.y);
961 ip.z = b.min.z;
962 }
963 }
964 }
965 else if (r.dir.z < 0)
966 {
967 if (r.pos.z < b.min.z)
968 return false;
969
970 T d = b.min.z - r.pos.z;
971
972 if (r.dir.z < -1 || d > TMAX * r.dir.z)
973 {
974 T t = d / r.dir.z;
975
976 if (tBackMin > t)
977 tBackMin = t;
978 }
979
980 if (r.pos.z >= b.max.z)
981 {
982 T d = b.max.z - r.pos.z;
983 T t = (r.dir.z < -1 || d > TMAX * r.dir.z)? d / r.dir.z: TMAX;
984
985 if (tFrontMax < t)
986 {
987 tFrontMax = t;
988
989 ip.x = clamp (r.pos.x + t * r.dir.x, b.min.x, b.max.x);
990 ip.y = clamp (r.pos.y + t * r.dir.y, b.min.y, b.max.y);
991 ip.z = b.max.z;
992 }
993 }
994 }
995 else // r.dir.z == 0
996 {
997 if (r.pos.z < b.min.z || r.pos.z > b.max.z)
998 return false;
999 }
1000
1001 return tFrontMax <= tBackMin;
1002 }
1003
1004
1005 template<class T>
1006 bool
intersects(const Box<Vec3<T>> & box,const Line3<T> & ray)1007 intersects (const Box< Vec3<T> > &box, const Line3<T> &ray)
1008 {
1009 Vec3<T> ignored;
1010 return intersects (box, ray, ignored);
1011 }
1012
1013
1014 IMATH_INTERNAL_NAMESPACE_HEADER_EXIT
1015
1016 #endif // INCLUDED_IMATHBOXALGO_H
1017