git.maemo.org Git - opencv/blob - 3rdparty/include/OpenEXR/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> &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
  77 namespace Imath {
  78
  79
  80 template <class T>
  81 inline 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
 104 template <class T>
 105 Vec3<T>
 106 closestPointInBox(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
 134 template <class T>
 135 Vec3<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
 246 template <class S, class T>
 247 Box< Vec3<S> >
 248 transform(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
 320 template <class T>
 321 Box< Vec3<T> >
 322 affineTransform(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
 434 template <class T>
 435 bool 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
 497 template<class T>
 498 bool 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
 601 template<class T>
 602 bool 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