/src/FreeImage/Source/OpenEXR/Imath/ImathBoxAlgo.h

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