Path: utzoo!attcan!uunet!mailrus!cs.utexas.edu!rutgers!njin!princeton!phoenix!gauss!markv
From: markv@gauss.Princeton.EDU (Mark VandeWettering)
Newsgroups: comp.graphics
Subject: Quadric Rendering (Re: Ray Tracing in 3D space)
Message-ID: <14638@phoenix.Princeton.EDU>
Date: 17 Mar 90 17:40:56 GMT
References: <14440@phoenix.Princeton.EDU> <20490@versatc.versatec.COM>
Sender: news@phoenix.Princeton.EDU
Reply-To: markv@gauss.Princeton.EDU (Mark VandeWettering)
Organization: Princeton University
Lines: 111

In article <20490@versatc.versatec.COM> ritter@versatc.versatec.COM (Jack Ritter) writes:

[ Lots of stuff about how he renders quadrics with raytracing ] 

I would be very interested in how people organize/optimize CSG trees for
things like quadrics.  Paper references would be welcome, and I will 
summarize to the net.

Recently, spurred on by some papers given to me by Pat Hanrahan, I have
taken an interest in quadric surfaces.  They really are quite interesting
beasties, and there are some really subtle tricks associated with them.
The math is even tractable for people who can only remember the quadratic
formula :-)

For instance, Jack used a formula like

	ax^2 + bxy +cxz + dy^2 + eyz + fz^2 + g = 0 

And that's fine for raytracing, because you would prefer to do all your
work in world space (at least MY raytracer does this (and I hate it ;-)).
It is also not exactly obvious how you might apply transformation matrices
to this form.

But you could also represent the quadric using homogeneous coordinates.
If points are represented by row vectors, then a quadric can be represented
by a 4x4 matrix (wow, just like homogeneous transformations!).  The quadric
satisfies the equation:

				T
	[ x y z w ] Q [x y z w]

where Q is of the form

	A 	B/2 	C/2 	D/2
	B/2	E 	F/2	G/2
	C/2	F/2	H	I/2
	D/2	G/2	I/2	J

The resulting quadric is

	Ax^2 + Bxy + Cxz + Dxw + Ey^2 + Fyz + Gyw + Hz^2 + Izw + Jw^2 == 0

Note that the Q matrix is symmetric.  

[ From now on, vectors followed by a small t will be transposes ]

Now, how to transform them:  Well, if a point P' = P M where M is a 
transformation matrix, we would like to have a new matrix Q ' such that
P' Q' P't == 0 as well. 	So:
	
	P Q Pt == 0 == P' Q' P't

	P Q Pt == (P M) Q' (P M)t
	P Q Pt == P M Q' Mt Pt
	  Q    == M Q' Mt
	  Q'   == Inv(M) Q Inv(M)t

Where Inv is the inverse.  Actually, since all multiples of a quadric are
equivalent, you can use the adjoint (matrix of cofactors) in place of the
inverse.  This has the advantage that it is always defined.

Now, how can you raytrace these?  Well, it may be much simpler if you do
the math with vectors and matrices rather than expanding it all out.

A ray is of the form:
	R = O + s D

If we substitute into our quadric equation:
	
	(O + s D) Q (O + s D)t == 0

	(OQ + sDQ)(Ot + (sD)t) == 0

	O Q Ot + O Q s Dt + s D Q Ot + s D Q s Dt == 0

[ and since Q is symmetric ]

	O Q Ot + 2 s O Q Dt + s^2 D Q Dt == 0

Which we can solve by the quadratic equation:

	a = D Q Dt
	b = 2 O Q Dt 
	c = O Q Ot 

And there you go.  I think that this is the way you should view ray/quadric
intersection, whether your implementation optimizes this or not.

The other big calculation is to determine the normal to the quadric at a 
particular point.  The normal to the surface is merely the partial derivatives
of the surface equation with respect to x, y, z, and w.  Well, golly, 
use the formula we wrote above:
   
Ax^2 + Bxy + Cxz + Dxw + Ey^2 + Fyz + Gyw + Hz^2 + Izw + Jw^2 == 0

Nt = 		2A x + By + Cz + Dw
		Bx + 2Ey + Fz + Gw
		Cx + Fy + 2Hz + Iw
		Dx + Gy + Iz + J

But if you look carefully, this is just (P Q)!!!!  (remember that 
all multiples of Q are equivalent)  So, all we have to do is evaluate
the normal with repect to the quadric matrix, and normalize the result.
(Pretty snazzy if you ask me)  This is a nice generalization of the 
special case of a sphere, where the point on the unit sphere IS the 
unit normal.

Anyway, I thought this was all quite interesting.  Anyone else have 
some nifty insights?

Mark