Holger Bettag (hobold@vectorizer.org)
I grant permission to redistribute this
document, as long as you give credit
where credit is due.
ABSTRACT: I will present an analytical and practical algorithm that finds the maximal and minimal values of an arbitrary quadric function restricted to a tetrahedral subvolume of 3-space. Mathematical proofs are only sketched, but pseudocode is included. All special cases of quadrics are handled naturally. Generalization to other convex polyhedra is straight forward.
The method is based on the observation that quadric functions, when restricted to compact domains, will assume extremal values either at the domain's border or at a single apex point. Degenerate quadrics may have zero or more than one apex points, but that implies that the extremal values will be assumed somewhere on the border.
The border of a polyhedron consists of polygons, the border of a polygon consists of edges, and the border of an edge are its two end points. I will therefore present the method by building up from one dimensional to higher dimensional cases.
Without loss of generality, we regard an arbitrary parabola:
f(x) := a x² + b x + c
on the interval [-1,+1], that is defined by the three values it assumes at the borders and the center:
f(-1) = L,
f(0) = M,
f(1) = R.
Solving for a, b, and c yields:
a = (L + R)/2 - M,
b = (R - L)/2,
c = M.
The apex of that parabola is the point p where the derivative is zero, which ultimately works out to:
p = -b/(2 a).
Case I: If a happens to equal zero, then the parabola has degenerated into a straight line. In that case, the extremal values are L and R, and we are done.
Case II: If a doesn't equal zero, then we check if the apex's position p is contained in the interval [-1,+1]. If p is outside, then the respective portion of the parabola is strictly increasing or strictly decreasing. So again the extremal values are L and R and we are done.
Case III: If the apex position p is located inside the interval [-1,+1], then only one of L and R is extremal, while the other extreme is the function value at the apex:
v = f(p).
To conclude the one dimensional case, I note that we had to check the borders and possibly one apex point inside the interval.
Without loss of generality, we regard an arbitrary quadric in two variables x and y:
q(x,y) = a x² + b y² + c x y + d x + e y + f
on the triangle with vertices (0,0), (1,0), and (0,1). Let the quadric be defined by the six values it assumes at the corners and the edges' midpoints:
q(0,0) = V1,
q(1,0) = V2,
q(0,1) = V3,
q(0.5,0) = M1,
q(0,0.5) = M2,
q(0.5,0.5) = M3.
Solving for the quadric coefficients yields:
a = 2 V1 + 2 V2 - 4 M1,
b = 2 V1 + 2 V3 - 4 M2,
c = 4 V1 - 4 M1 - 4 M2 + 4 M3,
d = -3 V1 - V2 + 4 M1,
e = -3 V1 - V3 + 4 M2,
f = V1.
The apex of the quadric is located at the single point where the gradient is zero, thus we arrive at the position (px,py), where
px = (c e - 2 b d)/(4 a b - c²), and
py = (c d - 2 a e)/(4 a b - c²).
Case I: If the denominator (4 a b - c²) is zero, then either there is no apex point, and we will find the extremal values on some edge, or
Case II: the denominator is zero and there are infinitely many "apexes". Those points must all lie on a straight line (or even the whole plane!), which (if it touches the triangle at all) will intersect at least one of the edges, so we will find the extremal value when we look at the edges.
Case III: If the denominator does not equal zero, we have to check the position of the apex to see if it is relevant. All of the conditions
px > 0,
py > 0, and
px + py < 1
must be met for the value q(px,py) to be taken into account. The apex could actually be a saddle point (in which case the true extremes are on the border), but there is nothing to be gained here by further analysis of special cases.
Case IV: If the apex is outside the triangle, the extremal values must be located on the edges.
To conclude the two dimensional case, I note that exactly like in the one dimensional case, we had to look at the borders and possibly one apex point inside the triangle.
Analogous to above, we have a quadric in three variables:
q(x,y,z) = a x² + b y² + c z² + d y z + e x z + f x y + g x + h y + i z + j
on the tetrahedron with corners (0,0,0), (1,0,0), (0,1,0), and (0,0,1). We specify the values at vertices and edge midpoints to define the quadric:
q(0,0,0) = V1,
q(1,0,0) = V2,
q(0,1,0) = V3,
q(0,0,1) = V4,
q(0.5,0,0) = M1,
q(0,0.5,0) = M2,
q(0,0,0.5) = M3,
q(0.5,0.5,0) = M4,
q(0.5,0,0.5) = M5,
q(0,0.5,0.5) = M6.
The quadric coefficients then are:
a = 2 V1 + 2 V2 - 4 M1,
b = 2 V1 + 2 V3 - 4 M2,
c = 2 V1 + 2 V4 - 4 M3,
d = 4 V1 - 4 M2 - 4 M3 + 4 M6,
e = 4 V1 - 4 M1 - 4 M3 + 4 M5,
f = 4 V1 - 4 M1 - 4 M2 + 4 M4,
g = -3 V1 - V2 + 4 M1,
h = -3 V1 - V3 + 4 M2,
i = -3 V1 - V4 + 4 M3,
j = V1.
Solving for a vanishing gradient again, we get:
denominator = 8 a b c + 2 d e f
- 2 a d² - 2 b e² - 2 c f²,
px = (2 b e i + d² g + 2 c f h - 4 b c g - d f i - d e h)/denominator,
py = (e² h + 2 a d i + 2 c f g - 4 a c h - d e g - e f i)/denominator,
and
pz = (2 b e g + 2 a d h + f² i - 4 a b i - e f h - d f g)/denominator.
Case I: If denominator is zero, we will find the extremes on the faces of the tetrahedron. Just like in two dimensions, the sub-cases here don't need to be explicitly distinguished by the algorithm, and for the same reasons.
Case II: If (px,py,pz) is not inside the tetrahedron, the extremal values will be on the border.
Case III: Finally, if (px,py,pz) lies inside the tetrahedron, the quadric's function value at that point has to be taken into account.
And that concludes my formal explanation along with the sketched proof. I admit it gets repetitive towards the end. :-)
You might wonder why I specified the quadrics through explicit function values in the formal chapters above. Doing that made the implementation a bit nicer, and removed the need to specify any vertex positions as function arguments. The topology of the tetrahedron is actually irrelevant (so long as it is not degenerate).
// auxiliary function for visibility test, one dimensional
// internal use only
static inline void refineMinMax1D (float& min, float& max,
float l, float m, float r) {
float a, b, c, x; // coefficients of parabola, position/value of apex
a = 0.5*(l + r) - m;
b = 0.5*(r - l);
c = m;
if (a != 0.0f) { // otherwise the parabola is really a line
x = -b/(2.0f*a);
if (x*x < 1) { // otherwise the apex is outside the interval
x = (a*x + b)*x + c;
if (x < min) {
min = x;
}
if (x > max) {
max = x;
}
}
}
return;
}
// auxiliary function for visibility test, two dimensional
// internal use only
static inline void refineMinMax2D (float& min, float& max,
float K1, float K2, float K3,
float K4, float K5, float K6) {
float a, b, c, d, e, f; // coefficients of quadric
float x, y, v; // position/value of apex
a = 2.0f*(K1 + K2 - 2.0f*K4);
b = 2.0f*(K1 + K3 - 2.0f*K5);
c = 4.0f*(K1 - K4 - K5 + K6);
d = -3.0f*K1 - K2 + 4.0f*K4;
e = -3.0f*K1 - K3 + 4.0f*K5;
f = K1;
if (4.0f*a*b - c*c != 0.0f) { // otherwise extremum is not a point
// position of apex
x = (c*e - 2.0f*b*d)/(4.0f*a*b - c*c);
y = (c*d - 2.0f*a*e)/(4.0f*a*b - c*c);
if ((x > 0.0f) && (y > 0.0f) && (x + y < 1.0f)) { // is the apex inside?
v = (a*x + c*y + d)*x + (b*y + e)*y + f;
if (v < min) {
min = v;
}
if (v > max) {
max = v;
}
}
}
return;
}
// test if tetrahedron contains quadric surface (isovalue zero)
//
// K1 to K4 are implicit function values at the four vertices,
// K5 is the value at the midpoint of edge K1-K2,
// K6 at edge K1-K3, K7 at edge K1-K4, then K2-K3 and so on.
//
bool isVisible(float K1, float K2, float K3,
float K4, float K5, float K6, float K7,
float K8, float K9, float K10) {
float min, max;
// upper bound of minimum from vertices (zero dimensional borders)
min = K1;
if (K2 < min) { min = K2; }
if (K3 < min) { min = K3; }
if (K4 < min) { min = K4; }
// lower bound of maximum from vertices (zero dimensional borders)
max = K1;
if (K2 > max) { max = K2; }
if (K3 > max) { max = K3; }
if (K4 > max) { max = K4; }
if (min*max <= 0.0f) {
return true; // early exit if we are bracketing zero
}
// refine bounds from edges (one dimensional borders)
refineMinMax1D(min, max, K1, K5, K2);
refineMinMax1D(min, max, K1, K6, K3);
refineMinMax1D(min, max, K1, K7, K4);
refineMinMax1D(min, max, K2, K8, K3);
refineMinMax1D(min, max, K2, K9, K4);
refineMinMax1D(min, max, K3, K10, K4);
if (min*max <= 0.0f) {
return true; // early exit if we are bracketing zero
}
// refine bounds from faces (two dimensional borders)
refineMinMax2D(min, max, K2, K3, K4, K8, K9, K10);
refineMinMax2D(min, max, K1, K3, K4, K6, K7, K10);
refineMinMax2D(min, max, K1, K2, K4, K5, K7, K9);
refineMinMax2D(min, max, K1, K2, K3, K5, K6, K8);
if (min*max <= 0.0f) {
return true; // early exit if we are bracketing zero
}
// refine bounds from volume (three dimensional)
float A, B, C, D, E, F, G, H, I, J; // quadric coefficients
float den; // denominator
float x, y, z; // position of apex
float v; // value at apex
A = 2.0f*(K1 + K2 - 2.0f*K5);
B = 2.0f*(K1 + K3 - 2.0f*K6);
C = 2.0f*(K1 + K4 - 2.0f*K7);
D = 4.0f*(K1 - K6 - K7 + K10);
E = 4.0f*(K1 - K5 - K7 + K9);
F = 4.0f*(K1 - K5 - K6 + K8);
G = -3.0f*K1 - K2 + 4.0f*K5;
H = -3.0f*K1 - K3 + 4.0f*K6;
I = -3.0f*K1 - K4 + 4.0f*K7;
J = K1;
den = 8.0f*A*B*C + 2.0f*(D*E*F - A*D*D - B*E*E - C*F*F);
if (den != 0.0f) {
// if apex is a single point, compute its location
x = (-4.0f*B*C*G - D*F*I - D*E*H + 2.0f*(B*E*I + C*F*H) + D*D*G) / den;
y = (-4.0f*A*C*H - D*E*G - E*F*I + E*E*H + 2.0f*(A*D*I + C*F*G)) / den;
z = (-4.0f*A*B*I - E*F*H - D*F*G + 2.0f*(B*E*G + A*D*H) + F*F*I) / den;
// is the apex inside?
if ((x > 0.0f) && (y > 0.0f) && (z > 0.0f) && (x + y + z < 1.0f)) {
v = (A*x + G + F*y)*x
+ (B*y + H + D*z)*y
+ (C*z + I + E*x)*z
+ J;
if (v < min) {
min = v;
}
if (v > max) {
max = v;
}
}
}
if (min*max <= 0.0f) {
return true; // min and max are indeed bracketing zero
}
return false;
}
As you can see, quite a bit of computation is necessary in general. However, only basic arithmetic is required, and the method arrives at the definitive answer (numerical issues notwithstanding) in constant time (early exit conditions notwithstanding). Maybe a smarter person than me can find some link between my method and Linear Programming, which would drastically reduce the number of candidate points that need to be evaluated.
Beware, the above code is only proven, not extensively tested. Permission is granted to use the above code, but you do so at your own risk. For the benefit of the poor guy who has to maintain it, you might want to put a reference to this paper into the comments somewhere.
The method immediately generalizes to axis aligned bounding boxes and other convex polyhedra. You just have to process more vertices and more faces, and the inside/outside tests for the apex points have to be adapted. Efficiency decreases with more complex polyhedra, but is always better than repeated application of the tetrahedral routine to a tessellated polyhedron.
An implementation for the more standard hexahedral spatial tiles is left as an exercise for the reader. Really. :-) All the boring formulas in sections one to three were written down so that you won't have to derive them again. Not only will you gain a deeper understanding of the method, but you can also tell your boss that there are no potential copyright issues. (Alternatively, you could look at this code fragment for inspiration.)
I presented a method that can analytically determine if an arbitrary quadric surface touches a given tetrahedron. No prior knowledge of the quadric surface's geometry is required, and all special cases are handled consistently and uniformly. Generalization to hexahedral tiles is straight forward. The method is not exceptionally efficient, but should be fast enough for practical applications such as the construction of a spatial subdivision acceleration structure of a ray tracer.
13-May-2012: added link to reference implementation.
3-Dec-2010: renamed auxiliary functions to better reflect their purpose.
added note to section 3, case I.
2-Dec-2010: initial public release.
1-Dec-2010: a few clarifications and minor edits.
30-Nov-2010: first HTML draft.