Ray-Object Intersection for Planes, Spheres, and Quadrics

CS 481 Lecture, Dr. Lawlor

In a uniform transparent medium, light travels in straight lines.

Straight lines have a very simple equation:

(1)	position_along_line = point_on_line + some_float * line_direction;
or	P = C + t * D;

This is called a parametric equation, because "some_float" is a free parameter. Raytracing is a way to draw arbitrary objects by solving for this floating-point parameter.

Ray-Plane Intersection

For example, consider a plane. If we specify the plane using a surface normal vector "plane_normal", the distance along this normal from the plane to the origin, then points on a plane satisfy this equation:

(2)	dot(point_in_plane,plane_normal) = distance_to_origin
or	dot(P,N) = k

(Rationale: moving in the plane is motion perpendicular to the normal. The dot product of perpendicular vectors is zero.)

If we want to find when the plane and line intersect, we just set:

	point_in_plane = position_along_line

This lets us substitute equation (1) into equation (2), giving:

	dot(point_on_line + some_float * line_direction,plane_normal) = distance_to_origin
or	dot(C+t*D,N) = k

which we have to solve for "some_float" (t). It's not immediately clear how to do this, because t is trapped inside the dot product. But linearity to the rescue! It turns out that dot product distributes over vector addition and scalar multiplication, so

	dot(C+t*D,N) = dot(C,N) + t*dot(D,N) = k

This is now just a linear equation in t, with solution:

	t=(k-dot(C,N))/dot(D,N)

We can then plug this parameter t back into the ray equation (1) to get the ray/line intersection point.

Note that if dot(D,N)==0, then the ray is parallel to the plane, and there is no solution for t. Also, depending on the orientation of the camera and plane, the camera ray may hit the plane behind the viewer, at negative t values. So in practice, you compute t as above, then check to see if it's a reasonable intersection. If so, we draw the object.

To actually draw the object, we also need a surface normal. This is really easy for a plane, because a plane only has one fixed value for its surface normal, and we've already had to specify it just to define the plane!

Ray-Sphere Intersection

Points on a sphere satisfy this equation:

(3)	length(point_on_sphere) = radius

Annoyingly, computing length takes a square root, which makes this equation difficult to solve. However, if we square both sides of this equation (radius is positive, so this will always work), we can express the length-squared as a dot product:

	radius^2 = length(point_on_sphere)^2 = dot(point_on_sphere,point_on_sphere)

Now we've reduced the square root business to just dot products. With shorter symbol names:

	r*r = dot(P,P)

Substitute in the ray equation (1) to find the ray/line intersection point:

	r*r = dot(C+t*D,C+t*D) = dot(C,C) + 2*dot(C,t*D) + dot(t*D,t*D)

or using dot product linearity again:

	r*r = dot(C+t*D,C+t*D) = dot(C,C) + 2*t*dot(C,D) + t*t*dot(D,D)

This is a quadratic equation in t, with constants c=dot(C,C)-r*r, b=2*dot(C,D), and a=dot(D,D). The t values are then a problem from high school algebra:

	t = (-b +- sqrt(b*b-4.0*a*c))/(2.0*a)

Note that there can be:

Two solutions, corresponding to the - and + versions of the square root, which are the ray entering and leaving the sphere respectively. Some of those solutions may be behind the camera!
Exactly one solution, if b*b-4.0*a*c==0. This corresponds to the ray just barely grazing the surface of the sphere (tangent).
No solutions, if b*b-4.0*a*c<0. This corresponds to the ray missing the sphere entirely.

A sphere's normal is very simple--draw a line from the center point (often the origin) to the intersection point you just computed. That's the normal vector.

Ray-Hyperboloid Intersections

Let's say we're looking for 3D points that satisfy the following odd equation (a hyperboloid)

	z^2 + k = x^2 + y^2

Substituting in the ray equation, we get:

	(C.z+t*D.z)^2 + k = (C.x + t*D.x)^2 + (C.y + t*D.y)^2
so	(C.z+t*D.z)^2 + k - (C.x + t*D.x)^2 - (C.y + t*D.y)^2 = 0

We've got to solve this for t. Each of the (C+t*D) terms expands out to C*C+2*t*C*D + t*t*D*D, so we have a quadratic with:

	c = C.z*C.z + k - C.x*C.x - C.y*C.y
	b = 2.0*(C.z*D.z - C.x*D.x - C.y*D.y)
	a = D.z*D.z-D.x*D.x-D.y*D.y

We can now apply the quadratic equation, exactly as before.

The surface normal of this shape is a little trickier to compute. One way to compute the surface normal is to write the defining equation as a space-filling function:

	f(P) = P.z^2 + k - P.x^2 - P.y^2

The shape itself consists of the set of points where f(P)=0. But the gradient of f points along the surface normal:

	N = +- normalize(vec3( df / dP.x, df / dP.y, df/dP.z ))
or	N = +- normalize(vec3( -2*P.x, -2*P.y, + 2*P.z ))

This gradient trick is handy way to extract normals for any surface that you have an equation for!

Ray-Object Intersection in General

Given a function f(P), we can raytrace the 3D surface f(P)=0 by plugging in the ray equation for P, and then solving f(C+t*D)=0 for t.

Once we find a point on the surface, we can compute surface normals from the gradient of f.

The only thing this *doesn't* work for is:

Objects whose function is too complex to solve analytically. Unfortunately, even simple functions like sin(x)-x=k have no elementary inverse function.
Objects that don't have a simple equation. For example, I can trivially write the equation for a sphere, a cube, etc; but I can't write one equation for the shape of my cat.

For these more complex objects, we'll need another approach--typically, we break the object down into simpler objects.

Quadric Surfaces

A quadric surface is defined as the f(P)=0 solution of this rather odd function:
f(P) = P^T A P

Here I'm treating the 3D vector P as a homogenous 4-vector, with the extra "w" coordinate set to 1.0 in the usual homogeneous coordinates fashion also used for perspective projection. To work properly with the 4x4 parameter matrix A, I've got to treat P as a column vector. P^T, the transpose of P, is a row vector.

If you expand out the rules of matrix multiplication, the above is exactly equivalent to:
f(P) = dot(P,A P)

That is, we first apply the matrix A to the point P, then dot the resulting vector with P. This gives us a scalar.

Here's a list of the terms of P that get applied to each entry of the matrix A. The corresponding entry of the matrix goes into the "__" spot.

f(x,y,z) =

__xx	__xy	__xz	__*x
__yx	__yy	__yz	__*y
__zx	__zy	__zz	__*z
__*x	__*y	__*z	__*1.0

So, for example, the hyperboloid we rendered above, "z² + k = x² + y²", could be represented as a quadric by rewriting the hyperboloid equation as "x² + y² - z² - k = 0", and hence setting the diagonal terms of A to +1, +1, -1, and -k.

Ray-Quadric Intersection

Following the usual technique for raytracing, we plug the ray equation P = C + t D into the object's function:
    f(P) = dot(P,A P)
    f(C + t D) = dot(C + t D, A (C + t D) )
    f(C + t D) = dot(C + t D, A C + t A D )
    f(C + t D) = dot(C,AC) + t dot(C,AD) + t dot(D,AC) + t² dot(D,AD)

This is just a quadratic equation, with coefficients:
    float a = dot(D,A*D);
    float b = dot(C,A*D) + dot(D,A*C);
    float c = dot(C,A*C)

You can then apply the quadratic equation normally. However, one thing that happens with quadrics is that the quadratic term "a" can be zero, which causes a divide by zero in the quadratic equation; or "a" can be very close to zero, causing serious roundoff problems. For these cases, it's better to switch to a linear solver, where instead of solving:
   0 = c + t*b + t²*a
you solve this linear equation instead:
   0 = c + t*b
This has only one solution, t=-c/b.

Transforming a Quadric by a Matrix

One very cool aspect of quadric objects is that a quadric transformed by a matrix is just another quadric:
     f(B P) = (B P)^T A (B P)

If you look up the definition of transpose, you find that (B P)^T = P^T B^T (transpose of a product is product of transposes in reverse order).   Hence:
     f(B P) = P^T B^T A (B P)

Or, moving the parenthesis around:
     f(B P) = P^T (B^T A B) P

Note that this is just a quadric with a new matrix (B^T A B). That is, we can evaluate f with points transformed by B simply by adjusting the quadric matrix to (B^T A B)!

In OpenGL, the modelview matrix transforms vertices on the object into world coordinates. Sadly, for the transformation above, we'd like to take world-coordinates points P, and transform them into object coordinates. This means B = inverse(gl_ModelViewMatrix). I've written some easy-to-call C++ code to invert a 4x4 matrix in "osl/mat4_inverse.h".

Quadric Normals

The surface normals for quadric can be found by taking partial derivatives of the object function:

f(x,y,z) =

__xx	__xy	__xz	__*x
__yx	__yy	__yz	__*y
__zx	__zy	__zz	__*z
__*x	__*y	__*z	__*1.0

Taking derivatives of f along each axis gives the gradient:

grad(f(x,y,z)) =vec3(df/dx,df/dy,df/dz) =

vec3(

__2x	__*y	__*z	__*1.0
__*y
__*z
__*1.0

	__*x
__*x	__2y	__*z	__*1.0
	__*z
	__*1.0

		__*x
		__*y
__*x	__*y	__2z	__*1.0
		__*1.0

)

The easy way to sum up these terms is to assemble a "gradient matrix" with the correct terms from the matrix summed up. You can then apply the gradient matrix to your hit point to get the gradient.

641 Students: Ray-Object Intersection Literature

Charles Loop and Jim Blinn. "Real-Time GPU Rendering of Piecewise Algebraic Surfaces". Siggraph Proceedings, Boston, 2006.
Presents an analytic method for raytracing algebraic surfaces, such as low-order Bezier patches.