From bogart@utah-gr.UUCP Thu Oct 20 02:42:04 1988
Path: leah!bingvaxu!sunybcs!boulder!ncar!mailrus!wasatch!utah-gr!bogart
From: bogart@utah-gr.UUCP (Rod G. Bogart)
Newsgroups: comp.graphics
Subject: Ray/Triangle Intersection with Barycentric Coordinates
Keywords: How we do it.
Message-ID: <2925@utah-gr.UUCP>
Date: 20 Oct 88 06:42:04 GMT
Organization: University of Utah CS Dept
Lines: 63


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)