shirley@uiucdcsm.cs.uiuc.edu (04/17/88)
The fastest way I know of to do this is to first find the point p
where the ray hits the plane. Then write down the following equation:
p = v1 + u * (v2 - v1) + v * (v3 - v1)
where v1, v2, v3 are the vertices of the triangle, and u,v are
parametric cooordinates on the plane containing the triangle.
The meaning of this equation is apparent in a diagram. Now rewrite
this equation to get three linear equations (from x/y/z coordinates)
with two unknowns (u/v). Use two equations and cramer's rule to
solve for u and v. Be sure to use the best two equations-- if the
triangle lies in the xy plane discard the z equation, etc.
Now that you have u and v, the point p is in the triangle iff
1. u > 0
2. v > 0
3. u+v < 1
Again, draw a diagram to see why this is true.
You can use the same equations for intersection with parallelogram,
except that condition three is replaced by
3. u < 1
4. v < 1
Hope this helps!
Peter Shirley
U of Illinois at UCroy@aero.ARPA (Roy Hashimoto) (04/19/88)
Here is my ray-triangle intersection algorithm. I won't claim
that it's unique, but I haven't seen it in any texts. On the
other hand, with thousands of people writing ray tracers, it's
hard to believe that anything is really new. Stop me if you've
heard this one...
A ray can be defined parametrically as
(x0,y0,z0) + t*(dx,dy,dz), t > 0
where (x0,y0,z0) is the initial point and (dx,dy,dz) is a direction
vector. The triangle primitive can similarly be partially defined
by
(xp,yp,zp) + t1*(v1x,v1y,v1z) + t2*(v2x,v2y,v2z), t1,t2 > 0
where (xp,yp,zp) is a vertex, and the two direction vectors are
obtained by subtracting that vertex from each of the remaining two.
Note that this is only a partial definition, as additional restrictions
must be placed on t1 and t2.
An intersection occurs where these vector equations are satisfied
simultaneously, and it is easy to write a set of equations in
t, t1, and t2. Only solutions with t, t1, and t2 all positive are
of interest, and this limitation is useful to pare the computation.
Solving the system with Cramer's formula requires the calculations
of four determinants. Careful ordering of the calculations involved
will reduce the amount of work for rays that miss.
Referring to the four determinants as D (coefficients only), D*t,
D*t1, and D*t2, by playing with the expressions for D and D*t, we
can see that:
1. They share the terms:
vy1*vz2 - vy2*vz1
vx1*vz2 - vx2*vz1
vx1*vy2 - vx2*vy1
2. These terms are independent of the ray.
3. These terms, with one sign change, form a normal vector of
the triangle.
Since a normal probably must be stored anyway for shading purposes,
why not store this one? It is true that usually a normalized vector
is more appropriate, so storing a scaling factor to convert to a
normalized vector would be helpful. Now D and D*t can be calculated
using 6 multiplies and 7 adds (I count subtracts=adds).
If D = 0, or D and D*t have opposite signs, there is no positive
t solution, and all remaining calculations can be eliminated. If
the signs are the same, compute D*t1. As far as I know, there are
no shortcuts, and D*t1 will take 9 multiplies and 5 adds.
If D*t1 has a different sign from D and D*t, there is no appropriate
solution, and this terminates the test. If not, go ahead with
D*t2, another 9 multiplies and 5 adds, and check the signs again.
Assuming all the signs match, the ray satisfies the conditions for
the partially defined triangle. The remaining constraint is:
t1 + t2 < 1
which is equivalent to
abs(D*t1) + abs(D*t2) < abs(D)
which requires 2 adds and maybe some sign changes. All the rays
that pass this stage have been determined to intersect the triangle.
The only remaining information required is the intersection point,
which can be found by finding t with 1 division and substituting
in our original ray equation, requiring 3 multiplies and 3 adds.
Making the following assumptions:
1. Comparison of D and D*t eliminates half the rays (they point
away from the plane containing the triangle).
2. Comparison with D*t1 eliminates half the remaining rays (wrong
half-plane).
3. Comparison with D*t2 eliminates two-thirds of the remaining
rays (a series of randomly chosen angles of random triangles
should average 60 degrees).
4. No rays are eliminated in the t1 + t2 < 1 constraint.
then the expected calculation per test is 13.0 multiplies, 11.2
adds, and 0.1 divisions.
It is very important to note that these assumptions work for
random rays and random triangles, which is not generally correct
for real ray tracers. The use of bounding boxes or other methods
of object or ray coherence will also invalidate these results.
However, I have used this algorithm for a number of different
approaches, and it is quite fast and doesn't require a lot of
additional storage.
I would appreciate hearing any comments on this technique, or,
better still, on any faster algorithms.
Roy T. Hashimoto
The Aerospace Corporation
roy@aerospace.aero.orgshirley@uiucdcsm.cs.uiuc.edu (04/21/88)
Roy Hashimoto's note reminded me that I forgot to mention a sidenote--
if you can constrain your ray to start at the origin (easy for first
generation rays-- otherwise it's no-go) you can take a shortcut:
The ray can now be described 'tv' instead of 'o + tv' (o = origin point,
v = view direction, t = parametric distance along ray). As before, the
parametric equation for the plane containing the triangle is:
p = p0 + k1*(p1 - p0) + k2*(p2 - p0)
Here p is a point on the plane that varies with k1 and k2, and p1/p2/p3
are the vertices. The conditions for inside/outside are the same as in
my last response (I apologize if I have changed notation).
We can generalize this to:
p = k0*p0 + k1*p1 + k2*p2
The point p is inside now if k0/k1/k2 are all positive (the condition
'k1 + k2 < 1' becomes '1 - (k1+k2) > 0' is 'k0 > 0').
Now assume that the ray hits the plane:
tv = k0*p0 + k1*p1 + k2*p2
since k0/k1/k2 are unknowns we can rewrite:
v = k0'*p0 + k1'*p1 + k2'*p2
or
|p0x p1x p2x| |k0'| |vx|
|p0y p1y p2y| * |k1'| = |vy|
|p0z p1z p2z| |k2'| |vz|
To solve :
|k0'| |p0x p1x p2x|-1 |vx|
|k1'| = |p0y p1y p2y| * |vy|
|k2'| |p0z p1z p2z| |vz|
If t is positive then the ray hits inside iff all k's are positive, otherwise
it hits only if they are all negative (a case we probably aren't
interested in for normal ray tracing applications).
Since we can precompute the inverted matrix, that gives a maximum of
9 multiplies, 6 adds, and 3 compares for each intersection. Usually we
can get away with fewer operations since all k's must be above 0.
Another thing I forgot to mention in my last response-- k1 and k2 are texture
coordinates for interpolation if colors or normal vectors are defined
at the vertices (if using the method here renormalize k's using condition
that k0 + k1 + k2 = 1. This also will give t). The interpolation
formula is
(1 - k1 - k2)*c0 + k1*c1 + k2*c2
As a final note, this method was shown to me by Hiroshi Degushi, whose
group implemented this on LINKS-1 hardware in Japan.
Peter Shirley