RayObject 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 floatingpoint parameter.
RayPlane 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=(kdot(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!
RaySphere 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 lengthsquared 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*b4.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*b4.0*a*c==0. This corresponds to the ray just
barely grazing the surface of the sphere (tangent).
 No solutions, if b*b4.0*a*c<0. This corresponds to the ray missing the sphere entirely.
A sphere's normal is very simpledraw a line from the center point
(often the origin) to the intersection point you just computed.
That's the normal vector.
RayHyperboloid 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.zD.x*D.xD.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 spacefilling 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!
RayObject 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 approachtypically, 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 4vector, 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) = 
__*x*x

__*x*y

__*x*z

__*x

__*y*x

__*y*y

__*y*z

__*y

__*z*x

__*z*y

__*z*z

__*z

__*x

__*y

__*z

__*1.0


So, for example, the hyperboloid we rendered above, "z^{2} + k = x^{2} + y^{2}", could be represented as a quadric by rewriting the hyperboloid equation as "x^{2} + y^{2}  z^{2}  k = 0", and hence setting the diagonal terms of A to +1, +1, 1, and k.
RayQuadric 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^{2} 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^{2}*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 worldcoordinates points P, and transform them into object
coordinates. This means B = inverse(gl_ModelViewMatrix).
I've written some easytocall 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) = 
__*x*x

__*x*y

__*x*z

__*x

__*y*x

__*y*y

__*y*z

__*y

__*z*x

__*z*y

__*z*z

__*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(

__*2*x

__*y

__*z

__*1.0

__*y




__*z




__*1.0





, 

__*x



__*x

__*2*y

__*z

__*1.0


__*z




__*1.0




, 


__*x




__*y


__*x

__*y

__*2*z

__*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: RayObject Intersection Literature
Charles Loop and Jim Blinn. "RealTime GPU Rendering of Piecewise Algebraic Surfaces". Siggraph Proceedings, Boston, 2006.
Presents an analytic method for raytracing algebraic surfaces, such as loworder Bezier patches.