Ray Tracing: Graphics for the Masses
Color trace_ray( Ray original_ray )
{- Color point_color, reflect_color, refract_color
Object obj
obj = get_first_intersection( original_ray )
point_color = get_point_color( obj )
if ( object is reflective )- reflect_color = trace_ray( get_reflected_ray( original_ray, obj ) )
if ( object is refractive )
- refract_color = trace_ray( get_refracted_ray( original_ray, obj ) )
return ( combine_colors( point_color, reflect_color, refract_color ))
}
Intersection of Rays/Segments with Triangles
Ray/Segment-Triangle Intersection
Consider a ray R (or a segment S) from
P0 to P1,
and a triangle T with vertices V0,
V1 and V2.
The triangle T lies in the plane p
through V0 with normal vector
n = (V1-V0)×(V2-V0).
To get the intersection of R (or S) and T, one first
determines the intersection of R (or S) and
p. If they do not
intersect, then the ray (or segment) also does not intersect T and we are
done. However, if they intersect in the point PI=P(rI),
we need to determine if this point is in the triangle T for there to be a
valid intersection.There are a number of ways to test for the
inclusion of a point inside a 3D planar triangle. The algorithms given by
[Badouel, 1990] and [O'Rourke, 1998] project the
point and triangle onto a 2D coordinate plane where inclusion is tested.
To implement these algorithms, one must choose a projection coordinate plane
which avoids a degenerate projection. This is done by excluding the
coordinate which has the largest component in the plane normal vector n [Synder
& Barr, 1987]. The intent is to reduce the 3D problem to a simpler 2D
problem which has an efficient solution. However, there is a small
overhead involved in selecting and applying the projection function. The
algorithm of [Moller-Trumbore, 1997] (MT) does not project into 2D, and
finds a solution using direct 3D computations. Testing with some complex
models shows that the MT algorithm is faster than the one by Badouel.
We present an alternate method which also uses direct 3D computations to
determine inclusion, avoiding projection onto a 2D coordinate plane. As a
result, the code is cleaner and more compact. Like [Moller-Trumbore,
1997], we use the parametric equation of
p relative
to T, but derive a different method of solution which computes the
parametric coordinates of the intersection point in the
plane. The parametric plane equation (see:
Planes in Algorithm
4) is given by:
where s and t are real
numbers, and u = (V1-V0)
and v = (V2-V0)
are edge vectors of T. A point P = V(s,t)
is in the triangle T when s>=0, t>=0, and
s+t<=1. So, given PI, one
just has to find the (sI,tI)-coordinate for
it, and then check these inequalities to verify inclusion in T. Further,
a point P = V(s,t) is on an edge of T if one of the conditions
s=0, t=0, or s+t=1 is true (each condition corresponds
to one edge). Also, the three vertices are given by: V0=V(0,0),
V1=V(1,0), and V2=V(0,1).
To solve for sI and tI,
we define a 3D generalization of Hill's "perp-dot" product
[Hill, 1994]. For an embedded 2D plane
p with a
normal vector n, and any vector a in the plane (that is, a
satisfies n·a = 0), define the "generalized perp operator
on p"
by: a^ = n×a.
Then, a^ is another
vector in the plane p
(since n·a^ = 0),
and it is perpendicular to a (since a·a^
= 0). Further, this perp operator is linear on vectors in
p;
that is, (Aa + Bb)^
= Aa^ + Bb^
where a and b are vectors in
p, and A
and B are scalars. Note that if
p is the 2D
xy-plane (z=0) with n = (0,0,1), then this is exactly the
same 2D perp operator given by [Hill, 1994].We can now use this generalized perp
operator to solve the plane's parametric equation for the intersection point PI.
Put w = (PI-V0) which is
another vector in p.
Then, we want to solve the equation: w = su + tv
for s and t. Take the dot product of both sides with v^
to get: w · v^ =
su · v^ + tv
· v^
= su · v^,
and solve for sI. Similarly, taking the dot
product with u^,
we get: w · u^
= su · u^ +
tv · u^
= tv · u^,
and solve for tI. We get:
and
The denominators are nonzero whenever the triangle T is
nondegenerate (that is, has a nonzero area). When T is degenerate,
it is either a segment or a point, and in either case does not uniquely define a
plane and the computed normal vector is (0,0,0). So, this case must be
treated separately as a 3D segment-segment intersection which was treated in Algorithm
5 (see: Lines and Segments).Altogether we have used 3 cross products which is a lot of
computation. However, for any three vectors a b c, the following
identity holds for left association of the cross product: (a × b)
× c = (a · c) b - (b · c) a.
[Note: the cross product is not associative, and so there is a different (but
similar) identity for right association]. Applying the left association
identity results in the simplifications:
And we can now compute sI
and tI using only dot products
as:
with 5 distinct dot products. We have arranged terms so
that the two denominators are the same and only need to be calculated once.// Assume that classes are
already given for the objects:
// Point and Vector with
// coordinates {float x, y, z;}
// operators for:
// == to test
equality
// != to test
inequality
// (Vector)0 =
(0,0,0) (null vector)
// Point
= Point ± Vector
// Vector =
Point - Point
// Vector =
Scalar * Vector (scalar product)
// Vector =
Vector * Vector (cross product)
// Line and Ray and Segment with defining
points {Point P0, P1;}
// (a Line is infinite, Rays and
Segments start at P0)
// (a Ray extends beyond P1, but a
Segment ends at P1)
// Plane with a point and a normal {Point V0; Vector
n;}
// Triangle with defining vertices {Point V0, V1, V2;}
// Polyline and Polygon with n vertices {int n;
Point *V;}
// (a Polygon has V[n]=V[0])
//===================================================================#define SMALL_NUM
0.00000001 // anything that avoids division overflow
// dot product (3D) which allows vector operations in arguments
#define dot(u,v) ((u).x * (v).x + (u).y * (v).y + (u).z * (v).z)//
intersect_RayTriangle(): intersect a ray with a 3D triangle
// Input: a ray R, and a triangle T
// Output: *I = intersection point (when it exists)
// Return: -1 = triangle is degenerate (a segment or point)
// 0 =
disjoint (no intersect)
// 1 =
intersect in unique point I1
// 2 =
are in the same plane
int
intersect_RayTriangle( Ray R, Triangle T, Point* I )
{
Vector u, v, n;
// triangle vectors
Vector dir, w0, w;
// ray vectors
float r, a, b;
// params to calc ray-plane intersect
// get triangle edge vectors and plane normal
u = T.V1 - T.V0;
v = T.V2 - T.V0;
n = u * v;
// cross product
if (n == (Vector)0)
// triangle is degenerate
return -1;
// do not deal with this case
dir = R.P1 - R.P0;
// ray direction vector
w0 = R.P0 - T.V0;
a = -dot(n,w0);
b = dot(n,dir);
if (fabs(b) < SMALL_NUM) { // ray is
parallel to triangle plane
if (a == 0)
// ray lies in triangle plane
return 2;
else return 0;
// ray disjoint from plane
}
// get intersect point of ray with triangle plane
r = a / b;
if (r < 0.0)
// ray goes away from triangle
return 0;
// => no intersect
// for a segment, also test if (r > 1.0) => no intersect
*I = R.P0 + r * dir;
// intersect point of ray and plane
// is I inside T?
float uu, uv, vv, wu, wv, D;
uu = dot(u,u);
uv = dot(u,v);
vv = dot(v,v);
w = *I - T.V0;
wu = dot(w,u);
wv = dot(w,v);
D = uv * uv - uu * vv;
// get and test parametric coords
float s, t;
s = (uv * wv - vv * wu) / D;
if (s < 0.0 || s > 1.0)
// I is outside T
return 0;
t = (uv * wu - uu * wv) / D;
if (t < 0.0 || (s + t) > 1.0) // I is outside T
return 0;
return 1;
// I is in T
}
Posted from Diigo. The rest of my favorite links are here.
