PBRT
/home/felix/UBC/projects/AdaptiveLightfieldSampling/pbrt_v2/src/3rdparty/ilmbase-1.0.2/ImathLineAlgo.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_IMATHLINEALGO_H
00038 #define INCLUDED_IMATHLINEALGO_H
00039 
00040 //------------------------------------------------------------------
00041 //
00042 //      This file contains algorithms applied to or in conjunction
00043 //      with lines (Imath::Line). These algorithms may require
00044 //      more headers to compile. The assumption made is that these
00045 //      functions are called much less often than the basic line
00046 //      functions or these functions require more support classes
00047 //
00048 //      Contains:
00049 //
00050 //      bool closestPoints(const Line<T>& line1,
00051 //                         const Line<T>& line2,
00052 //                         Vec3<T>& point1,
00053 //                         Vec3<T>& point2)
00054 //
00055 //      bool intersect( const Line3<T> &line,
00056 //                      const Vec3<T> &v0,
00057 //                      const Vec3<T> &v1,
00058 //                      const Vec3<T> &v2,
00059 //                      Vec3<T> &pt,
00060 //                      Vec3<T> &barycentric,
00061 //                      bool &front)
00062 //
00063 //      V3f
00064 //      closestVertex(const Vec3<T> &v0,
00065 //                    const Vec3<T> &v1,
00066 //                    const Vec3<T> &v2,
00067 //                    const Line3<T> &l)
00068 //
00069 //      V3f
00070 //      rotatePoint(const Vec3<T> p, Line3<T> l, float angle)
00071 //
00072 //------------------------------------------------------------------
00073 
00074 #include "ImathLine.h"
00075 #include "ImathVecAlgo.h"
00076 #include "ImathFun.h"
00077 
00078 namespace Imath {
00079 
00080 
00081 template <class T>
00082 bool
00083 closestPoints
00084     (const Line3<T>& line1,
00085      const Line3<T>& line2,
00086      Vec3<T>& point1,
00087      Vec3<T>& point2)
00088 {
00089     //
00090     // Compute point1 and point2 such that point1 is on line1, point2
00091     // is on line2 and the distance between point1 and point2 is minimal.
00092     // This function returns true if point1 and point2 can be computed,
00093     // or false if line1 and line2 are parallel or nearly parallel.
00094     // This function assumes that line1.dir and line2.dir are normalized.
00095     //
00096 
00097     Vec3<T> w = line1.pos - line2.pos;
00098     T d1w = line1.dir ^ w;
00099     T d2w = line2.dir ^ w;
00100     T d1d2 = line1.dir ^ line2.dir;
00101     T n1 = d1d2 * d2w - d1w;
00102     T n2 = d2w - d1d2 * d1w;
00103     T d = 1 - d1d2 * d1d2;
00104     T absD = abs (d);
00105 
00106     if ((absD > 1) ||
00107         (abs (n1) < limits<T>::max() * absD &&
00108          abs (n2) < limits<T>::max() * absD))
00109     {
00110         point1 = line1 (n1 / d);
00111         point2 = line2 (n2 / d);
00112         return true;
00113     }
00114     else
00115     {
00116         return false;
00117     }
00118 }
00119 
00120 
00121 template <class T>
00122 bool
00123 intersect
00124     (const Line3<T> &line,
00125      const Vec3<T> &v0,
00126      const Vec3<T> &v1,
00127      const Vec3<T> &v2,
00128      Vec3<T> &pt,
00129      Vec3<T> &barycentric,
00130      bool &front)
00131 {
00132     //
00133     // Given a line and a triangle (v0, v1, v2), the intersect() function
00134     // finds the intersection of the line and the plane that contains the
00135     // triangle.
00136     //
00137     // If the intersection point cannot be computed, either because the
00138     // line and the triangle's plane are nearly parallel or because the
00139     // triangle's area is very small, intersect() returns false.
00140     //
00141     // If the intersection point is outside the triangle, intersect
00142     // returns false.
00143     //
00144     // If the intersection point, pt, is inside the triangle, intersect()
00145     // computes a front-facing flag and the barycentric coordinates of
00146     // the intersection point, and returns true.
00147     //
00148     // The front-facing flag is true if the dot product of the triangle's
00149     // normal, (v2-v1)%(v1-v0), and the line's direction is negative.
00150     //
00151     // The barycentric coordinates have the following property:
00152     //
00153     //     pt = v0 * barycentric.x + v1 * barycentric.y + v2 * barycentric.z
00154     //
00155 
00156     Vec3<T> edge0 = v1 - v0;
00157     Vec3<T> edge1 = v2 - v1;
00158     Vec3<T> normal = edge1 % edge0;
00159 
00160     T l = normal.length();
00161 
00162     if (l != 0)
00163         normal /= l;
00164     else
00165         return false;   // zero-area triangle
00166 
00167     //
00168     // d is the distance of line.pos from the plane that contains the triangle.
00169     // The intersection point is at line.pos + (d/nd) * line.dir.
00170     //
00171 
00172     T d = normal ^ (v0 - line.pos);
00173     T nd = normal ^ line.dir;
00174 
00175     if (abs (nd) > 1 || abs (d) < limits<T>::max() * abs (nd))
00176         pt = line (d / nd);
00177     else
00178         return false;  // line and plane are nearly parallel
00179 
00180     //
00181     // Compute the barycentric coordinates of the intersection point.
00182     // The intersection is inside the triangle if all three barycentric
00183     // coordinates are between zero and one.
00184     //
00185 
00186     {
00187         Vec3<T> en = edge0.normalized();
00188         Vec3<T> a = pt - v0;
00189         Vec3<T> b = v2 - v0;
00190         Vec3<T> c = (a - en * (en ^ a));
00191         Vec3<T> d = (b - en * (en ^ b));
00192         T e = c ^ d;
00193         T f = d ^ d;
00194 
00195         if (e >= 0 && e <= f)
00196             barycentric.z = e / f;
00197         else
00198             return false; // outside
00199     }
00200 
00201     {
00202         Vec3<T> en = edge1.normalized();
00203         Vec3<T> a = pt - v1;
00204         Vec3<T> b = v0 - v1;
00205         Vec3<T> c = (a - en * (en ^ a));
00206         Vec3<T> d = (b - en * (en ^ b));
00207         T e = c ^ d;
00208         T f = d ^ d;
00209 
00210         if (e >= 0 && e <= f)
00211             barycentric.x = e / f;
00212         else
00213             return false; // outside
00214     }
00215 
00216     barycentric.y = 1 - barycentric.x - barycentric.z;
00217 
00218     if (barycentric.y < 0)
00219         return false; // outside
00220 
00221     front = ((line.dir ^ normal) < 0);
00222     return true;
00223 }
00224 
00225 
00226 template <class T>
00227 Vec3<T>
00228 closestVertex
00229     (const Vec3<T> &v0,
00230      const Vec3<T> &v1,
00231      const Vec3<T> &v2,
00232      const Line3<T> &l)
00233 {
00234     Vec3<T> nearest = v0;
00235     T neardot       = (v0 - l.closestPointTo(v0)).length2();
00236     
00237     T tmp           = (v1 - l.closestPointTo(v1)).length2();
00238 
00239     if (tmp < neardot)
00240     {
00241         neardot = tmp;
00242         nearest = v1;
00243     }
00244 
00245     tmp = (v2 - l.closestPointTo(v2)).length2();
00246     if (tmp < neardot)
00247     {
00248         neardot = tmp;
00249         nearest = v2;
00250     }
00251 
00252     return nearest;
00253 }
00254 
00255 
00256 template <class T>
00257 Vec3<T>
00258 rotatePoint (const Vec3<T> p, Line3<T> l, T angle)
00259 {
00260     //
00261     // Rotate the point p around the line l by the given angle.
00262     //
00263 
00264     //
00265     // Form a coordinate frame with <x,y,a>. The rotation is the in xy
00266     // plane.
00267     //
00268 
00269     Vec3<T> q = l.closestPointTo(p);
00270     Vec3<T> x = p - q;
00271     T radius = x.length();
00272 
00273     x.normalize();
00274     Vec3<T> y = (x % l.dir).normalize();
00275 
00276     T cosangle = Math<T>::cos(angle);
00277     T sinangle = Math<T>::sin(angle);
00278 
00279     Vec3<T> r = q + x * radius * cosangle + y * radius * sinangle; 
00280 
00281     return r;
00282 }
00283 
00284 
00285 } // namespace Imath
00286 
00287 #endif