source: NonGTP/OpenEXR/include/Imath/ImathBoxAlgo.h @ 855

Revision 855, 14.0 KB checked in by igarcia, 19 years ago (diff)
Line 
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> &line,
66//                      Vec3<T> result)
67//
68//      bool intersects(const Box<Vec3<T>> &box, const Line3<T> &line)
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 clip(const T& in, const Box<T>& box)
82{
83    //
84    //  Clip a point so that it lies inside the given bbox
85    //
86
87    T out;
88
89    for (int i=0; i<(int)box.min.dimensions(); i++)
90    {
91        if (in[i] < box.min[i]) out[i] = box.min[i];
92        else if (in[i] > box.max[i]) out[i] = box.max[i];
93        else out[i] = in[i];
94    }
95
96    return out;
97}
98
99
100//
101// Return p if p is inside the box.
102//
103 
104template <class T>
105Vec3<T>
106closestPointInBox(const Vec3<T>& p, const Box< Vec3<T> >& box )
107{
108    Imath::V3f b;
109
110    if (p.x < box.min.x)
111        b.x = box.min.x;
112    else if (p.x > box.max.x)
113        b.x = box.max.x;
114    else
115        b.x = p.x;
116
117    if (p.y < box.min.y)
118        b.y = box.min.y;
119    else if (p.y > box.max.y)
120        b.y = box.max.y;
121    else
122        b.y = p.y;
123
124    if (p.z < box.min.z)
125        b.z = box.min.z;
126    else if (p.z > box.max.z)
127        b.z = box.max.z;
128    else
129        b.z = p.z;
130
131    return b;
132}
133
134template <class T>
135Vec3<T> closestPointOnBox(const Vec3<T>& pt, const Box< Vec3<T> >& box )
136{
137    //
138    //  This sucker is specialized to work with a Vec3f and a box
139    //  made of Vec3fs.
140    //
141
142    Vec3<T> result;
143   
144    // trivial cases first
145    if (box.isEmpty())
146        return pt;
147    else if (pt == box.center())
148    {
149        // middle of z side
150        result[0] = (box.max[0] + box.min[0])/2.0;
151        result[1] = (box.max[1] + box.min[1])/2.0;
152        result[2] = box.max[2];
153    }
154    else
155    {
156        // Find the closest point on a unit box (from -1 to 1),
157        // then scale up.
158
159        // Find the vector from center to the point, then scale
160        // to a unit box.
161        Vec3<T> vec = pt - box.center();
162        T sizeX = box.max[0]-box.min[0];
163        T sizeY = box.max[1]-box.min[1];
164        T sizeZ = box.max[2]-box.min[2];
165
166        T halfX = sizeX/2.0;
167        T halfY = sizeY/2.0;
168        T halfZ = sizeZ/2.0;
169        if (halfX > 0.0)
170            vec[0] /= halfX;
171        if (halfY > 0.0)
172            vec[1] /= halfY;
173        if (halfZ > 0.0)
174            vec[2] /= halfZ;
175
176        // Side to snap side that has greatest magnitude in the vector.
177        Vec3<T> mag;
178        mag[0] = fabs(vec[0]);
179        mag[1] = fabs(vec[1]);
180        mag[2] = fabs(vec[2]);
181
182        result = mag;
183
184        // Check if beyond corners
185        if (result[0] > 1.0)
186            result[0] = 1.0;
187        if (result[1] > 1.0)
188            result[1] = 1.0;
189        if (result[2] > 1.0)
190            result[2] = 1.0;
191
192        // snap to appropriate side         
193        if ((mag[0] > mag[1]) && (mag[0] >  mag[2]))
194        {
195            result[0] = 1.0;
196        }
197        else if ((mag[1] > mag[0]) && (mag[1] >  mag[2]))
198        {
199            result[1] = 1.0;
200        }
201        else if ((mag[2] > mag[0]) && (mag[2] >  mag[1]))
202        {
203            result[2] = 1.0;
204        }
205        else if ((mag[0] == mag[1]) && (mag[0] == mag[2]))
206        {
207            // corner
208            result = Vec3<T>(1,1,1);
209        }
210        else if (mag[0] == mag[1])
211        {
212            // edge parallel with z
213            result[0] = 1.0;
214            result[1] = 1.0;
215        }
216        else if (mag[0] == mag[2])
217        {
218            // edge parallel with y
219            result[0] = 1.0;
220            result[2] = 1.0;
221        }
222        else if (mag[1] == mag[2])
223        {
224            // edge parallel with x
225            result[1] = 1.0;
226            result[2] = 1.0;
227        }
228
229        // Now make everything point the right way
230        for (int i=0; i < 3; i++)
231        {
232            if (vec[i] < 0.0)
233                result[i] = -result[i];
234        }
235
236        // scale back up and move to center
237        result[0] *= halfX;
238        result[1] *= halfY;
239        result[2] *= halfZ;
240
241        result += box.center();
242    }
243    return result;
244}
245
246template <class S, class T>
247Box< Vec3<S> >
248transform(const Box< Vec3<S> >& box, const Matrix44<T>& m)
249{
250    // Transforms Box3f by matrix, enlarging Box3f to contain result.
251    // Clever method courtesy of Graphics Gems, pp. 548-550
252    //
253    // This works for projection matrices as well as simple affine
254    // transformations.  Coordinates of the box are rehomogenized if there
255    // is a projection matrix
256
257    // a transformed empty box is still empty
258    if (box.isEmpty())
259        return box;
260
261    // If the last column is close enuf to ( 0 0 0 1 ) then we use the
262    // fast, affine version.  The tricky affine method could maybe be
263    // extended to deal with the projection case as well, but its not
264    // worth it right now.
265
266    if (m[0][3] * m[0][3] + m[1][3] * m[1][3] + m[2][3] * m[2][3]
267        + (1.0 - m[3][3]) * (1.0 - m[3][3]) < 0.00001)
268    {
269        // Affine version, use the Graphics Gems hack
270        int             i, j;
271        Box< Vec3<S> >  newBox;
272
273        for (i = 0; i < 3; i++)
274        {
275            newBox.min[i] = newBox.max[i] = (S) m[3][i];
276
277            for (j = 0; j < 3; j++)
278            {
279                float a, b;
280
281                a = (S) m[j][i] * box.min[j];
282                b = (S) m[j][i] * box.max[j];
283
284                if (a < b)
285                {
286                    newBox.min[i] += a;
287                    newBox.max[i] += b;
288                }
289                else
290                {
291                    newBox.min[i] += b;
292                    newBox.max[i] += a;
293                }
294            }
295        }
296
297        return newBox;
298    }
299
300    // This is a projection matrix.  Do things the naive way.
301    Vec3<S> points[8];
302
303    /* Set up the eight points at the corners of the extent */
304    points[0][0] = points[1][0] = points[2][0] = points[3][0] = box.min[0];
305    points[4][0] = points[5][0] = points[6][0] = points[7][0] = box.max[0];
306
307    points[0][1] = points[1][1] = points[4][1] = points[5][1] = box.min[1];
308    points[2][1] = points[3][1] = points[6][1] = points[7][1] = box.max[1];
309
310    points[0][2] = points[2][2] = points[4][2] = points[6][2] = box.min[2];
311    points[1][2] = points[3][2] = points[5][2] = points[7][2] = box.max[2];
312
313    Box< Vec3<S> > newBox;
314    for (int i = 0; i < 8; i++)
315        newBox.extendBy(points[i] * m);
316
317    return newBox;
318}
319
320template <class T>
321Box< Vec3<T> >
322affineTransform(const Box< Vec3<T> > &bbox, const Matrix44<T> &M)
323{
324    float       min0, max0, min1, max1, min2, max2, a, b;
325    float       min0new, max0new, min1new, max1new, min2new, max2new;
326
327    min0 = bbox.min[0];
328    max0 = bbox.max[0];
329    min1 = bbox.min[1];
330    max1 = bbox.max[1];
331    min2 = bbox.min[2];
332    max2 = bbox.max[2];
333
334    min0new = max0new = M[3][0];
335    a = M[0][0] * min0;
336    b = M[0][0] * max0;
337    if (a < b) {
338        min0new += a;
339        max0new += b;
340    } else {
341        min0new += b;
342        max0new += a;
343    }
344    a = M[1][0] * min1;
345    b = M[1][0] * max1;
346    if (a < b) {
347        min0new += a;
348        max0new += b;
349    } else {
350        min0new += b;
351        max0new += a;
352    }
353    a = M[2][0] * min2;
354    b = M[2][0] * max2;
355    if (a < b) {
356        min0new += a;
357        max0new += b;
358    } else {
359        min0new += b;
360        max0new += a;
361    }
362
363    min1new = max1new = M[3][1];
364    a = M[0][1] * min0;
365    b = M[0][1] * max0;
366    if (a < b) {
367        min1new += a;
368        max1new += b;
369    } else {
370        min1new += b;
371        max1new += a;
372    }
373    a = M[1][1] * min1;
374    b = M[1][1] * max1;
375    if (a < b) {
376        min1new += a;
377        max1new += b;
378    } else {
379        min1new += b;
380        max1new += a;
381    }
382    a = M[2][1] * min2;
383    b = M[2][1] * max2;
384    if (a < b) {
385        min1new += a;
386        max1new += b;
387    } else {
388        min1new += b;
389        max1new += a;
390    }
391
392    min2new = max2new = M[3][2];
393    a = M[0][2] * min0;
394    b = M[0][2] * max0;
395    if (a < b) {
396        min2new += a;
397        max2new += b;
398    } else {
399        min2new += b;
400        max2new += a;
401    }
402    a = M[1][2] * min1;
403    b = M[1][2] * max1;
404    if (a < b) {
405        min2new += a;
406        max2new += b;
407    } else {
408        min2new += b;
409        max2new += a;
410    }
411    a = M[2][2] * min2;
412    b = M[2][2] * max2;
413    if (a < b) {
414        min2new += a;
415        max2new += b;
416    } else {
417        min2new += b;
418        max2new += a;
419    }
420
421    Box< Vec3<T> > xbbox;
422
423    xbbox.min[0] = min0new;
424    xbbox.max[0] = max0new;
425    xbbox.min[1] = min1new;
426    xbbox.max[1] = max1new;
427    xbbox.min[2] = min2new;
428    xbbox.max[2] = max2new;
429
430    return xbbox;
431}
432
433
434template <class T>
435bool findEntryAndExitPoints(const Line3<T>& line,
436                            const Box<Vec3<T> >& box,
437                            Vec3<T> &enterPoint,
438                            Vec3<T> &exitPoint)
439{
440    if ( box.isEmpty() ) return false;
441    if ( line.distanceTo(box.center()) > box.size().length()/2. ) return false;
442
443    Vec3<T>     points[8], inter, bary;
444    Plane3<T>   plane;
445    int         i, v0, v1, v2;
446    bool        front = false, valid, validIntersection = false;
447
448    // set up the eight coords of the corners of the box
449    for(i = 0; i < 8; i++)
450    {
451        points[i].setValue( i & 01 ? box.min[0] : box.max[0],
452                            i & 02 ? box.min[1] : box.max[1],
453                            i & 04 ? box.min[2] : box.max[2]);
454    }
455
456    // intersect the 12 triangles.
457    for(i = 0; i < 12; i++)
458    {
459        switch(i)
460        {
461        case  0: v0 = 2; v1 = 1; v2 = 0; break;         // +z
462        case  1: v0 = 2; v1 = 3; v2 = 1; break;
463
464        case  2: v0 = 4; v1 = 5; v2 = 6; break;         // -z
465        case  3: v0 = 6; v1 = 5; v2 = 7; break;
466
467        case  4: v0 = 0; v1 = 6; v2 = 2; break;         // -x
468        case  5: v0 = 0; v1 = 4; v2 = 6; break;
469
470        case  6: v0 = 1; v1 = 3; v2 = 7; break;         // +x
471        case  7: v0 = 1; v1 = 7; v2 = 5; break;
472
473        case  8: v0 = 1; v1 = 4; v2 = 0; break;         // -y
474        case  9: v0 = 1; v1 = 5; v2 = 4; break;
475
476        case 10: v0 = 2; v1 = 7; v2 = 3; break;         // +y
477        case 11: v0 = 2; v1 = 6; v2 = 7; break;
478        }
479        if((valid=intersect (line, points[v0], points[v1], points[v2],
480                             inter, bary, front)) == true)
481        {
482            if(front == true)
483            {
484                enterPoint = inter;
485                validIntersection = valid;
486            }
487            else
488            {
489                exitPoint = inter;
490                validIntersection = valid;
491            }
492        }
493    }
494    return validIntersection;
495}
496
497template<class T>
498bool intersects(const Box< Vec3<T> > &box,
499                const Line3<T> &line,
500                Vec3<T> &result)
501{
502    /*
503       Fast Ray-Box Intersection
504       by Andrew Woo
505       from "Graphics Gems", Academic Press, 1990
506    */
507
508    const int right     = 0;
509    const int left      = 1;
510    const int middle    = 2;
511
512    const Vec3<T> &minB = box.min;
513    const Vec3<T> &maxB = box.max;
514    const Vec3<T> &origin = line.pos;
515    const Vec3<T> &dir = line.dir;
516
517    bool inside = true;
518    char quadrant[3];
519    int whichPlane;
520    float maxT[3];
521    float candidatePlane[3];
522
523    /* Find candidate planes; this loop can be avoided if
524        rays cast all from the eye(assume perpsective view) */
525    for (int i=0; i<3; i++)
526    {
527        if(origin[i] < minB[i])
528        {
529            quadrant[i] = left;
530            candidatePlane[i] = minB[i];
531            inside = false;
532        }
533        else if (origin[i] > maxB[i])
534        {
535            quadrant[i] = right;
536            candidatePlane[i] = maxB[i];
537            inside = false;
538        }
539        else   
540        {
541            quadrant[i] = middle;
542        }
543    }
544
545    /* Ray origin inside bounding box */
546    if ( inside )       
547    {
548        result = origin;
549        return true;
550    }
551
552
553        /* Calculate T distances to candidate planes */
554    for (int i = 0; i < 3; i++)
555    {
556        if (quadrant[i] != middle && dir[i] !=0.)
557        {
558            maxT[i] = (candidatePlane[i]-origin[i]) / dir[i];
559        }
560        else
561        {
562            maxT[i] = -1.;
563        }
564    }
565
566    /* Get largest of the maxT's for final choice of intersection */
567    whichPlane = 0;
568
569    for (int i = 1; i < 3; i++)
570    {
571        if (maxT[whichPlane] < maxT[i])
572        {
573            whichPlane = i;
574        }
575    }
576
577    /* Check final candidate actually inside box */
578    if (maxT[whichPlane] < 0.) return false;
579
580    for (int i = 0; i < 3; i++)
581    {
582        if (whichPlane != i)
583        {
584            result[i] = origin[i] + maxT[whichPlane] *dir[i];
585
586            if ((quadrant[i] == right && result[i] < minB[i]) ||
587                (quadrant[i] == left && result[i] > maxB[i]))
588            {
589                return false;   /* outside box */
590            }
591        }
592        else
593        {
594            result[i] = candidatePlane[i];
595        }
596    }
597
598    return true;
599}
600
601template<class T>
602bool intersects(const Box< Vec3<T> > &box, const Line3<T> &line)
603{
604    Vec3<T> ignored;
605    return intersects(box,line,ignored);
606}
607
608
609} // namespace Imath
610
611#endif
Note: See TracBrowser for help on using the repository browser.