PBRT
|
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