Ray tracing to general second order surfaces

First, let's look at a special case: ray tracing to an unit radius cylindrical tube.
The implicit function for a unit diameter infinitely long cylindrical tube along the z axis is:
x² + y² - 1 ≤ 0
We can solve the ray equation manually, like we did for the sphere, as follows:

Substituting in V+tW for (x,y,z):
(V_x + t W_x)² + (V_y + t W_y)² - 1 ≤ 0

Expanding out:
(W_x² + W_y²) t² + (2 V_x W_x + 2 V_y W_y) t + (V_x² + V_y² - 1) ≤ 0

Which just equals:
(W•W) t² + 2(V•W) t + (V•V - 1) ≤ 0

So to solve the quadratic equation in t:
A = W•W
B = 2 (V•W)
C = V•V - 1

But this doesn't give us the ability to transform the cylinder, nor would it generalize to other quadratic shapes, including cones, paraboloids and hyperboloids.
So instead, let's work on the general case of being able to transform any implicit quadratic function by a 4×4 matrix, and then ray trace to the transformed surface described by that transformed function.
The general second order implicit surface:
Any quadratic implicit surface can be expressed as follows:
ax² + by² + cz² + dyz + ezx + fxy + gx + hy + iz + j ≤ 0
The ten coefficients of this equation can be stored as a vector:
[a, b, c, d, e, f, g, h, i, j]
The unit sphere, the infinite unit cylinder and the infinite unit cone, respectively, can all be expressed in such a form:

Sphere: x² + y² + z² - 1 ≤ 0 [1, 1, 1, 0, 0, 0, 0, 0, 0, -1]

Cylinder: x² + y² - 1 ≤ 0 [1, 1, 0, 0, 0, 0, 0, 0, 0, -1]

Cone: x² + y² - z² ≤ 0 [1, 1, -1, 0, 0, 0, 0, 0, 0, 0]

Raytracing to a general second order implicit surface:

Ray tracing to a general second order surface just requires substituting the ray into this ten-coefficient equation to get a quadratic equation in t, exactly as we did above in the case of the cylinder:

Substituting in V+tW for (x,y,z):
a(V_x+tW_x)² + b(V_y+tW_y)² + c(V_z+tW_z)² + d(V_y+tW_y)(V_z+tW_z) + e(V_z+tW_z)(V_x+tW_x) + f(V_x+tW_x)(V_y+tW_y) + g(V_x+tW_x) + h(V_y+tW_y) + i(V_z+tW_z) + j ≤ 0

Expanding out:
a(V_x+tW_x)(V_x+tW_x) + b(V_y+tW_y)(V_y+tW_y) + c(V_z+tW_z)(V_z+tW_z) + d(V_y+tW_y)(V_z+tW_z) + e(V_z+tW_z)(V_x+tW_x) + f(V_x+tW_x)(V_y+tW_y) + g(V_x+tW_x) + h(V_y+tW_y) + i(V_z+tW_z) + j ≤ 0

Rearranging terms:
(aW_x² + bW_y² + cW_z² + dW_yW_z + eW_zW_x + fW_xW_y) t² + (2aV_xW_x + 2bV_yW_y + 2cV_zW_z + d(V_yW_z+V_zW_y) + e(V_zW_x+V_xW_z) + f(V_xW_y+V_yW_x) + gW_x + hW_y + iW_z) t + (aV_x² + bV_y² + cV_z² + dV_yV_z + eV_zV_x + fV_xV_y + gV_x + hV_y + iV_z + j) ≤ 0

So to solve the quadratic equation in t:

A = aW_x² + bW_y² + cW_z² + dW_yW_z + eW_zW_x + fW_xW_y

B = 2(aV_xW_x + bV_yW_y + cV_zW_z) + d(V_yW_z+V_zW_y) + e(V_zW_x+V_xW_z) + f(V_xW_y+V_yW_x) + gW_x + hW_y + iW_z

C = aV_x² + bV_y² + cV_z² + dV_yV_z + eV_zV_x + fV_xV_y + gV_x + hV_y + iV_z + j

Applying a linear transform to a second order surface

Let us first review the simpler case of linear surfaces: Half space P = (a, b, c, d) contains point p = (x,y,z), if and only if ax + by + cz + d ≤ 0. This can also be expressed as the product:

a b c d • x
y
z
1 ≤ 0

or, in other words, L p^T ≤ 0. To find the plane that contains transformed points M•p^T, you need to replace L with L•M^-1.

Similarly, we can express quadratic equation

ax² + by² + cz² + dyz + ezx + fxy + gx + hy + iz + j ≤ 0

as the following double product:

x y z 1 •

a f e g
0 b d h
0 0 c i
0 0 0 j
• x
y
z
1 ≤ 0

or in other words, p Q p^T = 0, where Q denotes the 10 quadratic coefficients arranged into the above 4×4 matrix.

This means that if you want to find the three-variable quadratic equation that contains transformed points M•p^T, you need to replace Q by (M^-1)^T Q M^-1 because that will give you:

(M p^T)^T (M^-1)^T Q M^-1 (M p^T) =
p^T M^T (M^-1)^T Q M^-1 M p^T =
p Q p^T

when you create the product of

(M^-1)^T P M^-1

you can end up with non-zero values in all sixteen slots of your 4×4 matrix. This is because in addition to x*y, x*z, x*1, y*z, y*1, z*1 terms, you also end up with y*x, z*y, 1*z, z*y, 1*z, 1*z terms, as follows:

x*x x*y x*z x*1
y*x y*y y*z y*1
z*x z*y z*z z*1
1*x 1*y 1*z 1*1

Since multiplication of numbers is commutative, you can just add the six numbers in the lower left (shown in magenta above) into their mirror locations on the upper right, as follows:

A E I M
B F J N
C G K O
D H L P
→

A E+B I+C M+D
0 F J+G N+H
0 0 K O+L
0 0 0 P

This will leave you with just the ten unique coefficients of the quadratic in three variables, in the upper right of your transformed matrix.

Computing the surface normal:

Once you have computed your transformed values for [a,b,c,d,e,f,g,h,i,j], you can use high school differential calculus to compute its vector valued derivative:

f(x,y,z) = ax² + by² + cz² + dyz + ezx + fxy + gx + hy + iz + j
f'(x,y,z) = [ 2ax + ez + fy + g , 2by + dz + fx + h , 2cz + dy + ex + i ]

Just normalize this derivative to get the normal at the surface point:

N = normalize(f'(x,y,z))