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

Revision 855, 9.2 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_IMATHLINEALGO_H
38#define INCLUDED_IMATHLINEALGO_H
39
40//------------------------------------------------------------------
41//
42//      This file contains algorithms applied to or in conjunction
43//      with lines (Imath::Line). These algorithms may require
44//      more headers to compile. The assumption made is that these
45//      functions are called much less often than the basic line
46//      functions or these functions require more support classes
47//
48//      Contains:
49//
50//      bool closestPoints(const Line<T>& line1,
51//                         const Line<T>& line2,
52//                         Vec3<T>& point1,
53//                         Vec3<T>& point2)
54//
55//      bool intersect( const Line3<T> &line,
56//                      const Vec3<T> &v0,
57//                      const Vec3<T> &v1,
58//                      const Vec3<T> &v2,
59//                      Vec3<T> &pt,
60//                      Vec3<T> &barycentric,
61//                      bool &front)
62//
63//      V3f
64//      closestVertex(const Vec3<T> &v0,
65//                    const Vec3<T> &v1,
66//                    const Vec3<T> &v2,
67//                    const Line3<T> &l)
68//
69//      V3f
70//      nearestPointOnTriangle(const Vec3<T> &v0,
71//                             const Vec3<T> &v1,
72//                             const Vec3<T> &v2,
73//                             const Line3<T> &l)
74//
75//      V3f
76//      rotatePoint(const Vec3<T> p, Line3<T> l, float angle)
77//
78//------------------------------------------------------------------
79
80#include <ImathLine.h>
81#include <ImathVecAlgo.h>
82
83namespace Imath {
84
85
86template <class T>
87bool closestPoints(const Line3<T>& line1,
88                   const Line3<T>& line2,
89                   Vec3<T>& point1,
90                   Vec3<T>& point2)
91{
92    //
93    //  Compute the closest points on two lines. This was originally
94    //  lifted from inventor. This function assumes that the line
95    //  directions are normalized. The original math has been collapsed.
96    //
97
98    T A = line1.dir ^ line2.dir;
99
100    if ( A == 1 ) return false;
101
102    T denom = A * A - 1;
103
104    T B = (line1.dir ^ line1.pos) - (line1.dir ^ line2.pos);
105    T C = (line2.dir ^ line1.pos) - (line2.dir ^ line2.pos);
106
107    point1 = line1(( B - A * C ) / denom);
108    point2 = line2(( B * A - C ) / denom);
109
110    return true;
111}
112
113
114
115template <class T>
116bool intersect( const Line3<T> &line,
117                const Vec3<T> &v0,
118                const Vec3<T> &v1,
119                const Vec3<T> &v2,
120                Vec3<T> &pt,
121                Vec3<T> &barycentric,
122                bool &front)
123{
124    //    Intersect the line with a triangle.
125    //    1. find plane of triangle
126    //    2. find intersection point of ray and plane
127    //    3. pick plane to project point and triangle into
128    //    4. check each edge of triangle to see if point is inside it
129
130    //
131    // XXX TODO - this routine is way too long
132    //          - the value of EPSILON is dubious
133    //          - there should be versions of this
134    //            routine that do not calculate the
135    //            barycentric coordinates or the
136    //            front flag
137
138    const float EPSILON = 1e-6;
139
140    T   d, t, d01, d12, d20, vd0, vd1, vd2, ax, ay, az, sense;
141    Vec3<T>     v01, v12, v20, c;
142    int         axis0, axis1;
143
144    // calculate plane for polygon
145    v01 = v1 - v0;
146    v12 = v2 - v1;
147
148    // c is un-normalized normal
149    c = v12.cross(v01);
150
151    d = c.length();
152    if(d < EPSILON)
153        return false;   // cant hit a triangle with no area
154    c = c * (1. / d);
155
156    // calculate distance to plane along ray
157
158    d = line.dir.dot(c);
159    if (d < EPSILON && d > -EPSILON)
160        return false;   // line is parallel to plane containing triangle
161
162    t = (v0 - line.pos).dot(c) / d;
163
164    if(t < 0)
165        return false;
166
167    // calculate intersection point
168    pt = line.pos + t * line.dir;
169
170    // is point inside triangle? Project to 2d to find out
171    // use the plane that has the largest absolute value
172    // component in the normal
173    ax = c[0] < 0 ? -c[0] : c[0];
174    ay = c[1] < 0 ? -c[1] : c[1];
175    az = c[2] < 0 ? -c[2] : c[2];
176
177    if(ax > ay && ax > az)
178    {
179        // project on x=0 plane
180
181        axis0 = 1;
182        axis1 = 2;
183        sense = c[0] < 0 ? -1 : 1;
184    }
185    else if(ay > az)
186    {
187        axis0 = 2;
188        axis1 = 0;
189        sense = c[1] < 0 ? -1 : 1;
190    }
191    else
192    {
193        axis0 = 0;
194        axis1 = 1;
195        sense = c[2] < 0 ? -1 : 1;
196    }
197
198    // distance from v0-v1 must be less than distance from v2 to v0-v1
199    d01 = sense * ((pt[axis0] - v0[axis0]) * v01[axis1]
200                 - (pt[axis1] - v0[axis1]) * v01[axis0]);
201
202    if(d01 < 0) return false;
203
204    vd2 = sense * ((v2[axis0] - v0[axis0]) * v01[axis1]
205                 - (v2[axis1] - v0[axis1]) * v01[axis0]);
206
207    if(d01 > vd2) return false;
208
209    // distance from v1-v2 must be less than distance from v1 to v2-v0
210    d12 = sense * ((pt[axis0] - v1[axis0]) * v12[axis1]
211                 - (pt[axis1] - v1[axis1]) * v12[axis0]);
212
213    if(d12 < 0) return false;
214
215    vd0 = sense * ((v0[axis0] - v1[axis0]) * v12[axis1]
216                 - (v0[axis1] - v1[axis1]) * v12[axis0]);
217
218    if(d12 > vd0) return false;
219
220    // calculate v20, and do check on final side of triangle
221    v20 = v0 - v2;
222    d20 = sense * ((pt[axis0] - v2[axis0]) * v20[axis1]
223                 - (pt[axis1] - v2[axis1]) * v20[axis0]);
224
225    if(d20 < 0) return false;
226
227    vd1 = sense * ((v1[axis0] - v2[axis0]) * v20[axis1]
228                 - (v1[axis1] - v2[axis1]) * v20[axis0]);
229
230    if(d20 > vd1) return false;
231
232    // vd0, vd1, and vd2 will always be non-zero for a triangle
233    // that has non-zero area (we return before this for
234    // zero area triangles)
235    barycentric = Vec3<T>(d12 / vd0, d20 / vd1, d01 / vd2);
236    front = line.dir.dot(c) < 0;
237
238    return true;
239}
240
241template <class T>
242Vec3<T>
243closestVertex(const Vec3<T> &v0,
244              const Vec3<T> &v1,
245              const Vec3<T> &v2,
246              const Line3<T> &l)
247{
248    Vec3<T> nearest = v0;
249    T neardot       = (v0 - l.closestPointTo(v0)).length2();
250   
251    T tmp           = (v1 - l.closestPointTo(v1)).length2();
252
253    if (tmp < neardot)
254    {
255        neardot = tmp;
256        nearest = v1;
257    }
258
259    tmp = (v2 - l.closestPointTo(v2)).length2();
260    if (tmp < neardot)
261    {
262        neardot = tmp;
263        nearest = v2;
264    }
265
266    return nearest;
267}
268
269template <class T>
270Vec3<T>
271nearestPointOnTriangle(const Vec3<T> &v0,
272                       const Vec3<T> &v1,
273                       const Vec3<T> &v2,
274                       const Line3<T> &l)
275{
276    Vec3<T> pt, barycentric;
277    bool front;
278
279    if (intersect (l, v0, v1, v2, pt, barycentric, front))
280        return pt;
281
282    //
283    // The line did not intersect the triangle, so to be picky, you should
284    // find the closest edge that it passed over/under, but chances are that
285    // 1) another triangle will be closer
286    // 2) the app does not need this much precision for a ray that does not
287    //    intersect the triangle
288    // 3) the expense of the calculation is not worth it since this is the
289    //    common case
290    //
291    // XXX TODO  This is bogus -- nearestPointOnTriangle() should do
292    //           what its name implies; it should return a point
293    //           on an edge if some edge is closer to the line than
294    //           any vertex.  If the application does not want the
295    //           extra calculations, it should be possible to specify
296    //           that; it is not up to this nearestPointOnTriangle()
297    //           to make the decision.
298
299    return closestVertex(v0, v1, v2, l);
300}
301
302template <class T>
303Vec3<T>
304rotatePoint(const Vec3<T> p, Line3<T> l, T angle)
305{
306    //
307    // Rotate the point p around the line l by the given angle.
308    //
309
310    //
311    // Form a coordinate frame with <x,y,a>. The rotation is the in xy
312    // plane.
313    //
314
315    Vec3<T> q = l.closestPointTo(p);
316    Vec3<T> x = p - q;
317    T radius = x.length();
318
319    x.normalize();
320    Vec3<T> y = (x % l.dir).normalize();
321
322    T cosangle = Math<T>::cos(angle);
323    T sinangle = Math<T>::sin(angle);
324
325    Vec3<T> r = q + x * radius * cosangle + y * radius * sinangle;
326
327    return r;
328}
329
330
331} // namespace Imath
332
333#endif
Note: See TracBrowser for help on using the repository browser.