PBRT
/home/felix/UBC/projects/AdaptiveLightfieldSampling/pbrt_v2/src/3rdparty/ilmbase-1.0.2/ImathBoxAlgo.h
00001 
00002 //
00003 // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
00004 // Digital Ltd. LLC
00005 // 
00006 // All rights reserved.
00007 // 
00008 // Redistribution and use in source and binary forms, with or without
00009 // modification, are permitted provided that the following conditions are
00010 // met:
00011 // *       Redistributions of source code must retain the above copyright
00012 // notice, this list of conditions and the following disclaimer.
00013 // *       Redistributions in binary form must reproduce the above
00014 // copyright notice, this list of conditions and the following disclaimer
00015 // in the documentation and/or other materials provided with the
00016 // distribution.
00017 // *       Neither the name of Industrial Light & Magic nor the names of
00018 // its contributors may be used to endorse or promote products derived
00019 // from this software without specific prior written permission. 
00020 // 
00021 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00022 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00023 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
00024 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
00025 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
00027 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
00028 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
00029 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
00030 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
00031 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00032 //
00034 
00035 
00036 
00037 #ifndef INCLUDED_IMATHBOXALGO_H
00038 #define INCLUDED_IMATHBOXALGO_H
00039 
00040 
00041 //---------------------------------------------------------------------------
00042 //
00043 //      This file contains algorithms applied to or in conjunction
00044 //      with bounding boxes (Imath::Box). These algorithms require
00045 //      more headers to compile. The assumption made is that these
00046 //      functions are called much less often than the basic box
00047 //      functions or these functions require more support classes.
00048 //
00049 //      Contains:
00050 //
00051 //      T clip<T>(const T& in, const Box<T>& box)
00052 //
00053 //      Vec3<T> closestPointOnBox(const Vec3<T>&, const Box<Vec3<T>>& )
00054 //
00055 //      Vec3<T> closestPointInBox(const Vec3<T>&, const Box<Vec3<T>>& )
00056 //
00057 //      void transform(Box<Vec3<T>>&, const Matrix44<T>&)
00058 //
00059 //      bool findEntryAndExitPoints(const Line<T> &line,
00060 //                                  const Box< Vec3<T> > &box,
00061 //                                  Vec3<T> &enterPoint,
00062 //                                  Vec3<T> &exitPoint)
00063 //
00064 //      bool intersects(const Box<Vec3<T>> &box, 
00065 //                      const Line3<T> &ray, 
00066 //                      Vec3<T> intersectionPoint)
00067 //
00068 //      bool intersects(const Box<Vec3<T>> &box, const Line3<T> &ray)
00069 //
00070 //---------------------------------------------------------------------------
00071 
00072 #include "ImathBox.h"
00073 #include "ImathMatrix.h"
00074 #include "ImathLineAlgo.h"
00075 #include "ImathPlane.h"
00076 
00077 namespace Imath {
00078 
00079 
00080 template <class T>
00081 inline T
00082 clip (const T &p, const Box<T> &box)
00083 {
00084     //
00085     // Clip the coordinates of a point, p, against a box.
00086     // The result, q, is the closest point to p that is inside the box.
00087     //
00088 
00089     T q;
00090 
00091     for (int i = 0; i < int (box.min.dimensions()); i++)
00092     {
00093         if (p[i] < box.min[i])
00094             q[i] = box.min[i];
00095         else if (p[i] > box.max[i])
00096             q[i] = box.max[i];
00097         else
00098             q[i] = p[i];
00099     }
00100 
00101     return q;
00102 }
00103 
00104 
00105 template <class T>
00106 inline T
00107 closestPointInBox (const T &p, const Box<T> &box)
00108 {
00109     return clip (p, box);
00110 }
00111 
00112 
00113 template <class T>
00114 Vec3<T>
00115 closestPointOnBox (const Vec3<T> &p, const Box< Vec3<T> > &box)
00116 {
00117     //
00118     // Find the point, q, on the surface of
00119     // the box, that is closest to point p.
00120     //
00121     // If the box is empty, return p.
00122     //
00123 
00124     if (box.isEmpty())
00125         return p;
00126 
00127     Vec3<T> q = closestPointInBox (p, box);
00128 
00129     if (q == p)
00130     {
00131         Vec3<T> d1 = p - box.min;
00132         Vec3<T> d2 = box.max - p;
00133 
00134         Vec3<T> d ((d1.x < d2.x)? d1.x: d2.x,
00135                    (d1.y < d2.y)? d1.y: d2.y,
00136                    (d1.z < d2.z)? d1.z: d2.z);
00137 
00138         if (d.x < d.y && d.x < d.z)
00139         {
00140             q.x = (d1.x < d2.x)? box.min.x: box.max.x;
00141         }
00142         else if (d.y < d.z)
00143         {
00144             q.y = (d1.y < d2.y)? box.min.y: box.max.y;
00145         }
00146         else
00147         {
00148             q.z = (d1.z < d2.z)? box.min.z: box.max.z;
00149         }
00150     }
00151 
00152     return q;
00153 }
00154 
00155 
00156 template <class S, class T>
00157 Box< Vec3<S> >
00158 transform (const Box< Vec3<S> > &box, const Matrix44<T> &m)
00159 {
00160     //
00161     // Transform a 3D box by a matrix, and compute a new box that
00162     // tightly encloses the transformed box.
00163     //
00164     // If m is an affine transform, then we use James Arvo's fast
00165     // method as described in "Graphics Gems", Academic Press, 1990,
00166     // pp. 548-550.
00167     //
00168 
00169     //
00170     // A transformed empty box is still empty
00171     //
00172 
00173     if (box.isEmpty())
00174         return box;
00175 
00176     //
00177     // If the last column of m is (0 0 0 1) then m is an affine
00178     // transform, and we use the fast Graphics Gems trick.
00179     //
00180 
00181     if (m[0][3] == 0 && m[1][3] == 0 && m[2][3] == 0 && m[3][3] == 1)
00182     {
00183         Box< Vec3<S> > newBox;
00184 
00185         for (int i = 0; i < 3; i++) 
00186         {
00187             newBox.min[i] = newBox.max[i] = (S) m[3][i];
00188 
00189             for (int j = 0; j < 3; j++) 
00190             {
00191                 float a, b;
00192 
00193                 a = (S) m[j][i] * box.min[j];
00194                 b = (S) m[j][i] * box.max[j];
00195 
00196                 if (a < b) 
00197                 {
00198                     newBox.min[i] += a;
00199                     newBox.max[i] += b;
00200                 }
00201                 else 
00202                 {
00203                     newBox.min[i] += b;
00204                     newBox.max[i] += a;
00205                 }
00206             }
00207         }
00208 
00209         return newBox;
00210     }
00211 
00212     //
00213     // M is a projection matrix.  Do things the naive way:
00214     // Transform the eight corners of the box, and find an
00215     // axis-parallel box that encloses the transformed corners.
00216     //
00217 
00218     Vec3<S> points[8];
00219 
00220     points[0][0] = points[1][0] = points[2][0] = points[3][0] = box.min[0];
00221     points[4][0] = points[5][0] = points[6][0] = points[7][0] = box.max[0];
00222 
00223     points[0][1] = points[1][1] = points[4][1] = points[5][1] = box.min[1];
00224     points[2][1] = points[3][1] = points[6][1] = points[7][1] = box.max[1];
00225 
00226     points[0][2] = points[2][2] = points[4][2] = points[6][2] = box.min[2];
00227     points[1][2] = points[3][2] = points[5][2] = points[7][2] = box.max[2];
00228 
00229     Box< Vec3<S> > newBox;
00230 
00231     for (int i = 0; i < 8; i++) 
00232         newBox.extendBy (points[i] * m);
00233 
00234     return newBox;
00235 }
00236 
00237 
00238 template <class S, class T>
00239 Box< Vec3<S> >
00240 affineTransform (const Box< Vec3<S> > &box, const Matrix44<T> &m)
00241 {
00242     //
00243     // Transform a 3D box by a matrix whose rightmost column
00244     // is (0 0 0 1), and compute a new box that tightly encloses
00245     // the transformed box.
00246     //
00247     // As in the transform() function, above, we use James Arvo's
00248     // fast method.
00249     //
00250 
00251     if (box.isEmpty())
00252     {
00253         //
00254         // A transformed empty box is still empty
00255         //
00256 
00257         return box;
00258     }
00259 
00260     Box< Vec3<S> > newBox;
00261 
00262     for (int i = 0; i < 3; i++) 
00263     {
00264         newBox.min[i] = newBox.max[i] = (S) m[3][i];
00265 
00266         for (int j = 0; j < 3; j++) 
00267         {
00268             float a, b;
00269 
00270             a = (S) m[j][i] * box.min[j];
00271             b = (S) m[j][i] * box.max[j];
00272 
00273             if (a < b) 
00274             {
00275                 newBox.min[i] += a;
00276                 newBox.max[i] += b;
00277             }
00278             else 
00279             {
00280                 newBox.min[i] += b;
00281                 newBox.max[i] += a;
00282             }
00283         }
00284     }
00285 
00286     return newBox;
00287 }
00288 
00289 
00290 template <class T>
00291 bool
00292 findEntryAndExitPoints (const Line3<T> &r,
00293                         const Box<Vec3<T> > &b,
00294                         Vec3<T> &entry,
00295                         Vec3<T> &exit)
00296 {
00297     //
00298     // Compute the points where a ray, r, enters and exits a box, b:
00299     //
00300     // findEntryAndExitPoints() returns
00301     //
00302     //     - true if the ray starts inside the box or if the
00303     //       ray starts outside and intersects the box
00304     //
00305     //     - false otherwise (that is, if the ray does not
00306     //       intersect the box)
00307     //
00308     // The entry and exit points are
00309     //
00310     //     - points on two of the faces of the box when
00311     //       findEntryAndExitPoints() returns true
00312     //       (The entry end exit points may be on either
00313     //       side of the ray's origin)
00314     //
00315     //     - undefined when findEntryAndExitPoints()
00316     //       returns false
00317     //
00318 
00319     if (b.isEmpty())
00320     {
00321         //
00322         // No ray intersects an empty box
00323         //
00324 
00325         return false;
00326     }
00327 
00328     //
00329     // The following description assumes that the ray's origin is outside
00330     // the box, but the code below works even if the origin is inside the
00331     // box:
00332     //
00333     // Between one and three "frontfacing" sides of the box are oriented
00334     // towards the ray's origin, and between one and three "backfacing"
00335     // sides are oriented away from the ray's origin.
00336     // We intersect the ray with the planes that contain the sides of the
00337     // box, and compare the distances between the ray's origin and the
00338     // ray-plane intersections.  The ray intersects the box if the most
00339     // distant frontfacing intersection is nearer than the nearest
00340     // backfacing intersection.  If the ray does intersect the box, then
00341     // the most distant frontfacing ray-plane intersection is the entry
00342     // point and the nearest backfacing ray-plane intersection is the
00343     // exit point.
00344     //
00345 
00346     const T TMAX = limits<T>::max();
00347 
00348     T tFrontMax = -TMAX;
00349     T tBackMin = TMAX;
00350 
00351     //
00352     // Minimum and maximum X sides.
00353     //
00354 
00355     if (r.dir.x >= 0)
00356     {
00357         T d1 = b.max.x - r.pos.x;
00358         T d2 = b.min.x - r.pos.x;
00359 
00360         if (r.dir.x > 1 ||
00361             (abs (d1) < TMAX * r.dir.x &&
00362              abs (d2) < TMAX * r.dir.x))
00363         {
00364             T t1 = d1 / r.dir.x;
00365             T t2 = d2 / r.dir.x;
00366 
00367             if (tBackMin > t1)
00368             {
00369                 tBackMin = t1;
00370 
00371                 exit.x = b.max.x; 
00372                 exit.y = clamp (r.pos.y + t1 * r.dir.y, b.min.y, b.max.y);
00373                 exit.z = clamp (r.pos.z + t1 * r.dir.z, b.min.z, b.max.z);
00374             }
00375 
00376             if (tFrontMax < t2)
00377             {
00378                 tFrontMax = t2;
00379 
00380                 entry.x = b.min.x; 
00381                 entry.y = clamp (r.pos.y + t2 * r.dir.y, b.min.y, b.max.y);
00382                 entry.z = clamp (r.pos.z + t2 * r.dir.z, b.min.z, b.max.z);
00383             }
00384         }
00385         else if (r.pos.x < b.min.x || r.pos.x > b.max.x)
00386         {
00387             return false;
00388         }
00389     }
00390     else // r.dir.x < 0
00391     {
00392         T d1 = b.min.x - r.pos.x;
00393         T d2 = b.max.x - r.pos.x;
00394 
00395         if (r.dir.x < -1 ||
00396             (abs (d1) < -TMAX * r.dir.x &&
00397              abs (d2) < -TMAX * r.dir.x))
00398         {
00399             T t1 = d1 / r.dir.x;
00400             T t2 = d2 / r.dir.x;
00401 
00402             if (tBackMin > t1)
00403             {
00404                 tBackMin = t1;
00405 
00406                 exit.x = b.min.x; 
00407                 exit.y = clamp (r.pos.y + t1 * r.dir.y, b.min.y, b.max.y);
00408                 exit.z = clamp (r.pos.z + t1 * r.dir.z, b.min.z, b.max.z);
00409             }
00410 
00411             if (tFrontMax < t2)
00412             {
00413                 tFrontMax = t2;
00414 
00415                 entry.x = b.max.x; 
00416                 entry.y = clamp (r.pos.y + t2 * r.dir.y, b.min.y, b.max.y);
00417                 entry.z = clamp (r.pos.z + t2 * r.dir.z, b.min.z, b.max.z);
00418             }
00419         }
00420         else if (r.pos.x < b.min.x || r.pos.x > b.max.x)
00421         {
00422             return false;
00423         }
00424     }
00425 
00426     //
00427     // Minimum and maximum Y sides.
00428     //
00429 
00430     if (r.dir.y >= 0)
00431     {
00432         T d1 = b.max.y - r.pos.y;
00433         T d2 = b.min.y - r.pos.y;
00434 
00435         if (r.dir.y > 1 ||
00436             (abs (d1) < TMAX * r.dir.y &&
00437              abs (d2) < TMAX * r.dir.y))
00438         {
00439             T t1 = d1 / r.dir.y;
00440             T t2 = d2 / r.dir.y;
00441 
00442             if (tBackMin > t1)
00443             {
00444                 tBackMin = t1;
00445 
00446                 exit.x = clamp (r.pos.x + t1 * r.dir.x, b.min.x, b.max.x);
00447                 exit.y = b.max.y; 
00448                 exit.z = clamp (r.pos.z + t1 * r.dir.z, b.min.z, b.max.z);
00449             }
00450 
00451             if (tFrontMax < t2)
00452             {
00453                 tFrontMax = t2;
00454 
00455                 entry.x = clamp (r.pos.x + t2 * r.dir.x, b.min.x, b.max.x);
00456                 entry.y = b.min.y; 
00457                 entry.z = clamp (r.pos.z + t2 * r.dir.z, b.min.z, b.max.z);
00458             }
00459         }
00460         else if (r.pos.y < b.min.y || r.pos.y > b.max.y)
00461         {
00462             return false;
00463         }
00464     }
00465     else // r.dir.y < 0
00466     {
00467         T d1 = b.min.y - r.pos.y;
00468         T d2 = b.max.y - r.pos.y;
00469 
00470         if (r.dir.y < -1 ||
00471             (abs (d1) < -TMAX * r.dir.y &&
00472              abs (d2) < -TMAX * r.dir.y))
00473         {
00474             T t1 = d1 / r.dir.y;
00475             T t2 = d2 / r.dir.y;
00476 
00477             if (tBackMin > t1)
00478             {
00479                 tBackMin = t1;
00480 
00481                 exit.x = clamp (r.pos.x + t1 * r.dir.x, b.min.x, b.max.x);
00482                 exit.y = b.min.y; 
00483                 exit.z = clamp (r.pos.z + t1 * r.dir.z, b.min.z, b.max.z);
00484             }
00485 
00486             if (tFrontMax < t2)
00487             {
00488                 tFrontMax = t2;
00489 
00490                 entry.x = clamp (r.pos.x + t2 * r.dir.x, b.min.x, b.max.x);
00491                 entry.y = b.max.y; 
00492                 entry.z = clamp (r.pos.z + t2 * r.dir.z, b.min.z, b.max.z);
00493             }
00494         }
00495         else if (r.pos.y < b.min.y || r.pos.y > b.max.y)
00496         {
00497             return false;
00498         }
00499     }
00500 
00501     //
00502     // Minimum and maximum Z sides.
00503     //
00504 
00505     if (r.dir.z >= 0)
00506     {
00507         T d1 = b.max.z - r.pos.z;
00508         T d2 = b.min.z - r.pos.z;
00509 
00510         if (r.dir.z > 1 ||
00511             (abs (d1) < TMAX * r.dir.z &&
00512              abs (d2) < TMAX * r.dir.z))
00513         {
00514             T t1 = d1 / r.dir.z;
00515             T t2 = d2 / r.dir.z;
00516 
00517             if (tBackMin > t1)
00518             {
00519                 tBackMin = t1;
00520 
00521                 exit.x = clamp (r.pos.x + t1 * r.dir.x, b.min.x, b.max.x);
00522                 exit.y = clamp (r.pos.y + t1 * r.dir.y, b.min.y, b.max.y);
00523                 exit.z = b.max.z; 
00524             }
00525 
00526             if (tFrontMax < t2)
00527             {
00528                 tFrontMax = t2;
00529 
00530                 entry.x = clamp (r.pos.x + t2 * r.dir.x, b.min.x, b.max.x);
00531                 entry.y = clamp (r.pos.y + t2 * r.dir.y, b.min.y, b.max.y);
00532                 entry.z = b.min.z; 
00533             }
00534         }
00535         else if (r.pos.z < b.min.z || r.pos.z > b.max.z)
00536         {
00537             return false;
00538         }
00539     }
00540     else // r.dir.z < 0
00541     {
00542         T d1 = b.min.z - r.pos.z;
00543         T d2 = b.max.z - r.pos.z;
00544 
00545         if (r.dir.z < -1 ||
00546             (abs (d1) < -TMAX * r.dir.z &&
00547              abs (d2) < -TMAX * r.dir.z))
00548         {
00549             T t1 = d1 / r.dir.z;
00550             T t2 = d2 / r.dir.z;
00551 
00552             if (tBackMin > t1)
00553             {
00554                 tBackMin = t1;
00555 
00556                 exit.x = clamp (r.pos.x + t1 * r.dir.x, b.min.x, b.max.x);
00557                 exit.y = clamp (r.pos.y + t1 * r.dir.y, b.min.y, b.max.y);
00558                 exit.z = b.min.z; 
00559             }
00560 
00561             if (tFrontMax < t2)
00562             {
00563                 tFrontMax = t2;
00564 
00565                 entry.x = clamp (r.pos.x + t2 * r.dir.x, b.min.x, b.max.x);
00566                 entry.y = clamp (r.pos.y + t2 * r.dir.y, b.min.y, b.max.y);
00567                 entry.z = b.max.z; 
00568             }
00569         }
00570         else if (r.pos.z < b.min.z || r.pos.z > b.max.z)
00571         {
00572             return false;
00573         }
00574     }
00575 
00576     return tFrontMax <= tBackMin;
00577 }
00578 
00579 
00580 template<class T>
00581 bool
00582 intersects (const Box< Vec3<T> > &b, const Line3<T> &r, Vec3<T> &ip)
00583 {
00584     //
00585     // Intersect a ray, r, with a box, b, and compute the intersection
00586     // point, ip:
00587     //
00588     // intersect() returns
00589     //
00590     //     - true if the ray starts inside the box or if the
00591     //       ray starts outside and intersects the box
00592     //
00593     //     - false if the ray starts outside the box and intersects it,
00594     //       but the intersection is behind the ray's origin.
00595     //
00596     //     - false if the ray starts outside and does not intersect it
00597     //
00598     // The intersection point is
00599     //
00600     //     - the ray's origin if the ray starts inside the box
00601     //
00602     //     - a point on one of the faces of the box if the ray
00603     //       starts outside the box
00604     //
00605     //     - undefined when intersect() returns false
00606     //
00607 
00608     if (b.isEmpty())
00609     {
00610         //
00611         // No ray intersects an empty box
00612         //
00613 
00614         return false;
00615     }
00616 
00617     if (b.intersects (r.pos))
00618     {
00619         //
00620         // The ray starts inside the box
00621         //
00622 
00623         ip = r.pos;
00624         return true;
00625     }
00626 
00627     //
00628     // The ray starts outside the box.  Between one and three "frontfacing"
00629     // sides of the box are oriented towards the ray, and between one and
00630     // three "backfacing" sides are oriented away from the ray.
00631     // We intersect the ray with the planes that contain the sides of the
00632     // box, and compare the distances between ray's origin and the ray-plane
00633     // intersections.
00634     // The ray intersects the box if the most distant frontfacing intersection
00635     // is nearer than the nearest backfacing intersection.  If the ray does
00636     // intersect the box, then the most distant frontfacing ray-plane
00637     // intersection is the ray-box intersection.
00638     //
00639 
00640     const T TMAX = limits<T>::max();
00641 
00642     T tFrontMax = -1;
00643     T tBackMin = TMAX;
00644 
00645     //
00646     // Minimum and maximum X sides.
00647     //
00648 
00649     if (r.dir.x > 0)
00650     {
00651         if (r.pos.x > b.max.x)
00652             return false;
00653 
00654         T d = b.max.x - r.pos.x;
00655 
00656         if (r.dir.x > 1 || d < TMAX * r.dir.x)
00657         {
00658             T t = d / r.dir.x;
00659 
00660             if (tBackMin > t)
00661                 tBackMin = t;
00662         }
00663 
00664         if (r.pos.x <= b.min.x)
00665         {
00666             T d = b.min.x - r.pos.x;
00667             T t = (r.dir.x > 1 || d < TMAX * r.dir.x)? d / r.dir.x: TMAX;
00668 
00669             if (tFrontMax < t)
00670             {
00671                 tFrontMax = t;
00672 
00673                 ip.x = b.min.x; 
00674                 ip.y = clamp (r.pos.y + t * r.dir.y, b.min.y, b.max.y);
00675                 ip.z = clamp (r.pos.z + t * r.dir.z, b.min.z, b.max.z);
00676             }
00677         }
00678     }
00679     else if (r.dir.x < 0)
00680     {
00681         if (r.pos.x < b.min.x)
00682             return false;
00683 
00684         T d = b.min.x - r.pos.x;
00685 
00686         if (r.dir.x < -1 || d > TMAX * r.dir.x)
00687         {
00688             T t = d / r.dir.x;
00689 
00690             if (tBackMin > t)
00691                 tBackMin = t;
00692         }
00693 
00694         if (r.pos.x >= b.max.x)
00695         {
00696             T d = b.max.x - r.pos.x;
00697             T t = (r.dir.x < -1 || d > TMAX * r.dir.x)? d / r.dir.x: TMAX;
00698 
00699             if (tFrontMax < t)
00700             {
00701                 tFrontMax = t;
00702 
00703                 ip.x = b.max.x; 
00704                 ip.y = clamp (r.pos.y + t * r.dir.y, b.min.y, b.max.y);
00705                 ip.z = clamp (r.pos.z + t * r.dir.z, b.min.z, b.max.z);
00706             }
00707         }
00708     }
00709     else // r.dir.x == 0
00710     {
00711         if (r.pos.x < b.min.x || r.pos.x > b.max.x)
00712             return false;
00713     }
00714 
00715     //
00716     // Minimum and maximum Y sides.
00717     //
00718 
00719     if (r.dir.y > 0)
00720     {
00721         if (r.pos.y > b.max.y)
00722             return false;
00723 
00724         T d = b.max.y - r.pos.y;
00725 
00726         if (r.dir.y > 1 || d < TMAX * r.dir.y)
00727         {
00728             T t = d / r.dir.y;
00729 
00730             if (tBackMin > t)
00731                 tBackMin = t;
00732         }
00733 
00734         if (r.pos.y <= b.min.y)
00735         {
00736             T d = b.min.y - r.pos.y;
00737             T t = (r.dir.y > 1 || d < TMAX * r.dir.y)? d / r.dir.y: TMAX;
00738 
00739             if (tFrontMax < t)
00740             {
00741                 tFrontMax = t;
00742 
00743                 ip.x = clamp (r.pos.x + t * r.dir.x, b.min.x, b.max.x);
00744                 ip.y = b.min.y; 
00745                 ip.z = clamp (r.pos.z + t * r.dir.z, b.min.z, b.max.z);
00746             }
00747         }
00748     }
00749     else if (r.dir.y < 0)
00750     {
00751         if (r.pos.y < b.min.y)
00752             return false;
00753 
00754         T d = b.min.y - r.pos.y;
00755 
00756         if (r.dir.y < -1 || d > TMAX * r.dir.y)
00757         {
00758             T t = d / r.dir.y;
00759 
00760             if (tBackMin > t)
00761                 tBackMin = t;
00762         }
00763 
00764         if (r.pos.y >= b.max.y)
00765         {
00766             T d = b.max.y - r.pos.y;
00767             T t = (r.dir.y < -1 || d > TMAX * r.dir.y)? d / r.dir.y: TMAX;
00768             
00769             if (tFrontMax < t)
00770             {
00771                 tFrontMax = t;
00772 
00773                 ip.x = clamp (r.pos.x + t * r.dir.x, b.min.x, b.max.x);
00774                 ip.y = b.max.y; 
00775                 ip.z = clamp (r.pos.z + t * r.dir.z, b.min.z, b.max.z);
00776             }
00777         }
00778     }
00779     else // r.dir.y == 0
00780     {
00781         if (r.pos.y < b.min.y || r.pos.y > b.max.y)
00782             return false;
00783     }
00784 
00785     //
00786     // Minimum and maximum Z sides.
00787     //
00788 
00789     if (r.dir.z > 0)
00790     {
00791         if (r.pos.z > b.max.z)
00792             return false;
00793 
00794         T d = b.max.z - r.pos.z;
00795 
00796         if (r.dir.z > 1 || d < TMAX * r.dir.z)
00797         {
00798             T t = d / r.dir.z;
00799 
00800             if (tBackMin > t)
00801                 tBackMin = t;
00802         }
00803 
00804         if (r.pos.z <= b.min.z)
00805         {
00806             T d = b.min.z - r.pos.z;
00807             T t = (r.dir.z > 1 || d < TMAX * r.dir.z)? d / r.dir.z: TMAX;
00808             
00809             if (tFrontMax < t)
00810             {
00811                 tFrontMax = t;
00812 
00813                 ip.x = clamp (r.pos.x + t * r.dir.x, b.min.x, b.max.x);
00814                 ip.y = clamp (r.pos.y + t * r.dir.y, b.min.y, b.max.y);
00815                 ip.z = b.min.z; 
00816             }
00817         }
00818     }
00819     else if (r.dir.z < 0)
00820     {
00821         if (r.pos.z < b.min.z)
00822             return false;
00823 
00824         T d = b.min.z - r.pos.z;
00825 
00826         if (r.dir.z < -1 || d > TMAX * r.dir.z)
00827         {
00828             T t = d / r.dir.z;
00829 
00830             if (tBackMin > t)
00831                 tBackMin = t;
00832         }
00833 
00834         if (r.pos.z >= b.max.z)
00835         {
00836             T d = b.max.z - r.pos.z;
00837             T t = (r.dir.z < -1 || d > TMAX * r.dir.z)? d / r.dir.z: TMAX;
00838             
00839             if (tFrontMax < t)
00840             {
00841                 tFrontMax = t;
00842 
00843                 ip.x = clamp (r.pos.x + t * r.dir.x, b.min.x, b.max.x);
00844                 ip.y = clamp (r.pos.y + t * r.dir.y, b.min.y, b.max.y);
00845                 ip.z = b.max.z; 
00846             }
00847         }
00848     }
00849     else // r.dir.z == 0
00850     {
00851         if (r.pos.z < b.min.z || r.pos.z > b.max.z)
00852             return false;
00853     }
00854 
00855     return tFrontMax <= tBackMin;
00856 }
00857 
00858 
00859 template<class T>
00860 bool
00861 intersects (const Box< Vec3<T> > &box, const Line3<T> &ray)
00862 {
00863     Vec3<T> ignored;
00864     return intersects (box, ray, ignored);
00865 }
00866 
00867 
00868 } // namespace Imath
00869 
00870 #endif