One of the key concepts behind realtime raytracing is very effective spatial subdivision to limit the number of surfaces that have to be tested against each ray. I am not sure if there is a standard, but kD-Trees and variants thereof are being used widely.
To make any of these methods work for quadrics, the computer has to determine if a given axis aligned bounding box overlaps a given quadric surface. A practical algorithm must not miss any quadric surface that does in fact overlap the box. In traditional raytracing, it was tolerable if the algorithm occasionally claimed a hit when a surface did in fact miss the box. In realtime raytracing, we cannot afford such false positives, because they cause unnecessary computation, costing valuable time.
Thanks to the relative simplicity of quadrics, the overlap problem can be solved analytically. This idea has been looming in my mind for quite a while before I could finally prove it (thinking in implicit surfaces is an unusual thing to do for computer graphics). Applying the idea to axis aligned boxes took me more than one attempt, but eventually I arrived at a reasonably concise procedure. The code has survived extensive testing against a brute force implementation, and is currently part of Boar. It has been slightly cleaned up to serve as an illustration here:
/*
Intersecting an Arbitrary Quadric with an Axis Aligned Box
(an illustrative code fragment)
An explanation of the geometry and math behind this algorithm can be found at
http://www.vectorizer.org/IntersectQuadricTetrahedron.html
Feel free to send comments, improvements and questions to
"hobold@vectorizer.org". To my knowledge this algorithm is new, and
unencumbered by intellectual property rights.
This code fragment is Copyright 2012 by Holger Bettag, licensed under GPLv2.
If the GPLv2 is too restrictive for you, contact me. I just want to learn
about improvements you might find. I am sure we can agree on something that
won't cost you a dime, and doesn't force all the obligations of the GPLv2 on
you. The GPLv2 is merely the default.
*/
class Quadric {
public:
double A, B, C; // A x^2 + B y^2 + C z^2
double D, E, F; // + D yz + E zx + F xy
double G, H, I; // + G x + H y + I z
double J; // + J
}
// handle apex of a parabola
// internal use
static
inline void refineMinMax1D(double& minval, double& maxval,
double a, double b, double c,
double tmin, double tmax) {
if(a != 0.0) {
double pt = -b/(2.0*a);
// clamp to limits
pt = qmin(pt, tmax);
pt = qmax(pt, tmin);
double value = (a*pt + b)*pt + c;
minval = qmin(minval, value);
maxval = qmax(maxval, value);
}
return;
}
// handle apex of a 2D quadric and the two edges u = umin, u = umax and
// the one corner at u = umax, v = vmin
// internal use
static
void refineMinMax2D(double& minval, double& maxval,
double a, double b, double c, double d, double e, double f,
double umin, double umax, double vmin, double vmax) {
// check 2D apex point
double denominator = 4.0*a*b - c*c;
double value;
if (denominator != 0.0) {
double pu = (c*e - 2.0*b*d)/denominator;
double pv = (c*d - 2.0*a*e)/denominator;
// clamp to limits
pu = qmin(pu, umax);
pu = qmax(pu, umin);
pv = qmin(pv, vmax);
pv = qmax(pv, vmin);
value = (a*pu + c*pv + d)*pu + (b*pv + e)*pv + f;
minval = qmin(minval, value);
maxval = qmax(maxval, value);
}
// handle border at u = umin
refineMinMax1D(minval, maxval,
b, c*umin + e, (a*umin + d)*umin + f,
vmin, vmax);
// ... and border at u = umax
// extract coefficients of parabola (for reuse)
c = c*umax + e;
a = (a*umax + d)*umax + f;
refineMinMax1D(minval, maxval, b, c, a, vmin, vmax);
// look at corner at u = umax, v = vmin
value = (b*vmin + c)*vmin + a;
minval = qmin(minval, value);
maxval = qmax(maxval, value);
return;
}
// intersect a quadric with an axis aligned box
//
// returns 2 (VoxInside) if the box is contained in the quadric,
// 1 (VoxOutside) if the box is completely outside the quadric
// 0 (VoxSurface) if the box overlaps the quadric surface
VoxClass Quadric::intersectVoxel(const Vec3& Min, const Vec3& Max) const {
// handle the coordinate minimum corner and the coordinate maximum corner
double value = evaluatePoint(Min);
double minval = evaluatePoint(Max);
double maxval = qmax(minval, value);
minval = qmin(minval, value);
// process the faces and edges of the box
// (each of the six face calls handles two of the 12 edges and
// one of the 8 corners; the remaining two corners are already done)
// extract face at z = zmin as 2 dimensional quadric and process it
// q(u,v) := Q(u,v,zmin)
refineMinMax2D(minval, maxval,
A, B, F, E*Min.z + G, D*Min.z + H, (C*Min.z + I)*Min.z + J,
Min.x, Max.x, Min.y, Max.y);
// q(u,v) := Q(u,v,zmax)
refineMinMax2D(minval, maxval,
A, B, F, E*Max.z + G, D*Max.z + H, (C*Max.z + I)*Max.z + J,
Min.x, Max.x, Min.y, Max.y);
// y = ymin (switched u and v to q(v,ymin,u) for symmetry)
// q(u,v) := Q(v,ymin,u)
refineMinMax2D(minval, maxval,
C, A, E, D*Min.y + I, F*Min.y + G, (B*Min.y + H)*Min.y + J,
Min.z, Max.z, Min.x, Max.x);
// q(u,v) := Q(v,ymax,u)
refineMinMax2D(minval, maxval,
C, A, E, D*Max.y + I, F*Max.y + G, (B*Max.y + H)*Max.y + J,
Min.z, Max.z, Min.x, Max.x);
// x = xmin
// q(u,v) := Q(xmin,u,v)
refineMinMax2D(minval, maxval,
B, C, D, F*Min.x + H, E*Min.x + I, (A*Min.x + G)*Min.x + J,
Min.y, Max.y, Min.z, Max.z);
// q(u,v) := Q(xmax,u,v)
refineMinMax2D(minval, maxval,
B, C, D, F*Max.x + H, E*Max.x + I, (A*Max.x + G)*Max.x + J,
Min.y, Max.y, Min.z, Max.z);
// find the 3D apex point
double denominator = 8.0*A*B*C + 2.0*(D*E*F - A*D*D - B*E*E - C*F*F);
if (denominator != 0.0) {
double px = (2.0*B*E*I + D*D*G + 2.0*C*F*H - 4.0*B*C*G - D*F*I - D*E*H);
double py = (E*E*H + 2.0*A*D*I + 2.0*C*F*G - 4.0*A*C*H - D*E*G - E*F*I);
double pz = (2.0*B*E*G + 2.0*A*D*H + F*F*I - 4.0*A*B*I - E*F*H - D*F*G);
px /= denominator;
py /= denominator;
pz /= denominator;
// clamp position to limits
px = qmin(px, Max.x);
px = qmax(px, Min.x);
py = qmin(py, Max.y);
py = qmax(py, Min.y);
pz = qmin(pz, Max.z);
pz = qmax(pz, Min.z);
value = evaluatePoint(Vec3(px, py, pz));
minval = qmin(minval, value);
maxval = qmax(maxval, value);
}
return VoxClass(((maxval < 0.0) << 1) + (minval > 0.0));
}
// brute force method, for debugging & test reference
VoxClass Quadric::intersectVoxelDBG(const Vec3& Min, const Vec3& Max) const {
const int size = 100;
int x,y,z;
double min = evaluatePoint(Min);
double max = min;
double val;
double fx, fy, fz;
for (z = size; z >= 0; z--) {
fz = Min.z + double(z)*(Max.z - Min.z)/double(size);
for (y = size; y >= 0; y--) {
fy = Min.y + double(y)*(Max.y - Min.y)/double(size);
for (x = size; x >= 0; x--) {
fx = Min.x + double(x)*(Max.x - Min.x)/double(size);
val = evaluatePoint(Vec3(fx,fy,fz));
min = qmin(min, val);
max = qmax(max, val);
}
}
}
if (min > 0.0)
return voxOutside; // positive minimum value => voxel is outside
if (max < 0.0)
return voxInside; // negative maximum value => voxel is inside
return voxSurface; // otherwise voxel must straddle quadric surface
}
I admit it is rather lengthy, and quite a bit of computation is required. But for now I am pretending that the construction of the spatial subdivision structure is not in the critical path.