Raytracing Quadric Surfaces

2010, Dr. Lawlor, CS 481/681, CS, UAF

A quadric surface is defined as the zero set 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've been rendering, "z² + k = x² + y²", could be represented as a quadric by rewriting the hyperboloid equation as "x² + y² - z² - k", 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.