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

A quadric surface is defined as the zero set of this rather odd function:
f(P) = PT 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.  PT, 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've been rendering, "z2 + k = x2 + y2", could be represented as a quadric by rewriting the hyperboloid equation as "x2 + y2 - z2 - k", and hence setting the diagonal terms of A to +1, +1, -1, and -k.

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) + t2 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 + t2*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 = PT BT  (transpose of a product is product of transposes in reverse order).   Hence:
f(B P) = PT BT A (B P)

Or, moving the parenthesis around:
f(B P) = PT (BT A B) P

Note that this is just a quadric with a new matrix (BT A B).  That is, we can evaluate f with points transformed  by B simply by adjusting the quadric matrix to (BT 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".

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: