bogart@utah-gr.UUCP (Rod G. Bogart) (10/20/88)
A while back, there was a posting concerning ray/triangle
intersection. The goal was to determine if a ray intersects a
triangle, and if so, what are the barycentric coordinates. For the
uninitiated, barycentric coordinates are three values (r,s,t) all in
the range zero to one. Also, the sum of the values is one. These
values can be used as interpolation parameters for data which is known
at the triangle vertices (i.e. normals, colors, uv).
The algorithm presented previously involved a matrix inversion. The
math went something like this: Since (r,s,t) are interpolation
values, then the intersection point (P) must be a combination of the
triangle vertices scaled by (r,s,t).
[ x1 y1 z1 ] [ r ] [ Px ] [ r ] [ Px ]
[ x2 y2 z2 ] [ s ] = [ Py ] [ s ] = [ Py ] ~V
[ x3 y3 z3 ] [ t ] [ Pz ] [ t ] [ Pz ]
So, by inverting the vertex matrix (V -> ~V), and given any point in
the plane of the triangle, we can determine (r,s,t). If they are in
the range zero to one, the point is in the triangle.
The only problem with this method is numerical instability. If one
vertex is the origin, the matrix won't invert. If the triangle lies
in a coordinate plane, the matrix won't invert. In fact, for any
triangle which lies in a plane through the origin, the matrix won't
invert. (The vertex vectors don't span R3.) The reason this method
is so unstable is because it tries to solve a 2D problem in 3D. Once
the ray/plane intersection point is known, the barycentric coordinates
solution is a 2D issue.
Another way to think of barycentric coordinates is by the relative
areas of the subtriangles defined by the intersection point and the
triangle vertices.
1 If the area of triangle 123 is A, then the area of
/|\ P23 is rA. Area 12P is sA and area 1P3 is tA.
/ | \ With this image, it is obvious that r+s+t must equal
/ | \ one. If r, s, or t go outside the range zero to one,
/ t | s \ P will be outside the triangle.
/ _-P-_ \
/ _- -_ \
/_- r -_\
3---------------2
By using the above are relationships, the following equations define
r, s and t.
N = triangle normal = (vec(1 2) cross vec(1 3))
(vec(1 P) cross vec(1 3)) dot N
s = -------------------------------
length N
(vec(1 2) cross vec(1 P)) dot N
t = -------------------------------
length N
r = 1 - (s + t)
In actual code, it is better to avoid the divide and the square root.
So, you can set s equal to the numerator, and then test if s is less
than zero or greater than sqr(length N). For added efficiency,
preprocess the data and store sqr( length N) in the triangle data
structure. Even for extremely long thin triangles, this method is
accurate and numerically stable.
RGB Living life in the fast lane, eight items or less.
({ihnp4,decvax}!utah-cs!bogart, bogart@cs.utah.edu)arenberg@trwrb.UUCP (Jeff Arenberg) (10/21/88)
Ok, here is how I handle this calculation in my ray tracing program. I think it is quite effecient. Let a triangle be represented in the following manner : |\ | \ p1 | \ | \ O ------------> |________\ p0 p2 where p0 is the vector from the origin to one vertex and p1, p2 are the vectors from the first vertex to the other two vertices. Let N = p1 X p2 be the normal to the triangle. ------- | p1 X p2 | Construct the matrices b = | p1 | , bb = inv(b) = | bb[0] | | p2 | | bb[1] | | N | | bb[2] | and store away bb. Let the intersecting ray be parameterizes as r = t * D + P Now you can quickly intersect the ray with the triangle using the following psuedo code. ( . means vector dot product) Den = D . bb[2] if (Den == 0) then ray parallel to triangle plane, so return Num = (p0 - P) . bb[2] t = Num / Den if (t <= 0) then on or behind triangle, so return p = t * D + P - p0 a = p . bb[0] b = p . bb[1] if (a < 0.0 || b < 0.0 || a + b > 1.0) then not in triangle and return b1 = 1 - a - b /* barycentric coordinates */ b2 = a b3 = b The idea here is that the matrix bb transforms to a coordinate frame where the sides of the triangle form the X,Y axes and the normal the Z axis of the frame and the sides have been scaled to unit length. The variable Den represents the dZ component of the ray in this frame. If dZ is zero, then the ray must be parallel to the X,Y plane. Num is the Z location of the ray origin in the new frame and t is simply the parameter in both frames required to intersect the ray with the triangle's plane. Once t is known, the intersection point is found in the original frame, saved for latter use, and the X,Y coordinates of this point are found in the triangle's frame. A simple comparison is then made to determine if the point is inside the triangle. The barycenter coordinates are also easily found. I haven't seen this algorithm in any of the literature, but then I haven't really looked either. If anyone knows if this approach has been published before, I'd really like to know about it. Jeff Arenberg ------------------------------------------------------------- UUCP : ( ucbvax, ihnp4, uscvax ) !trwrb!csed-pyramid!arenberg GEnie: shifty -------------------------------------------------------------