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) =
__*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, "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", 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^{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 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) =
__*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.