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.