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