Minimum bounding sphere of a Frustum

If you cant get it work you could always take all the corners and find the min and max points, and get the center using those (minpoint+maxpoint)/2

Then get the radius from the distance between center and either of those points...

You could try calculating the circumcircle from any 3 of the points. Because the shape is symmetrical the left out point should also be on the circle, satisfying your requirement for a minimum bounding sphere.

@Waterlimon: yes I know that, that will do the work, but as long as I'm using this for shadow mapping, I really need the minimum sphere.
@ Postie: ok, I'll try that and keep you posted, ...

I'm interested in the solution to this as well.

Let E be the eyepoint. Let the view direction be D, the up direction be U, and the right direction be R. D, U, and R are unit length and mutually perpendicular. The frustum values are dmin (near), dmax (far), rmax (right), -rmax (left), umax (up), and -umax (down). The sphere center is C = E + s*D for some s. A vertex of the near plane is P0 = E + dmin*D + umax*U + rmax*R. A vertex of the far plane is P1 = E + dmax*D + (dmax*umax/dmin)*U + (dmax*rmax/dmin)*R. The squared distances from P0 to C and P1 to C must be equal, so (s-dmin)^2 + rmax^2 + umax^2 = (s - dmax)^2 + (dmax*rmax/dmin)^2 + (dmax*umax/dmin)^2. This leads to s = ((dmin + dmax)/2)*(1 + (rmax/dmin)^2 + (umax/dmin)^2).

check this.

This might be better.

1) Dot(P - A, P - A) = Dot(P - B, P - B)
2) P = S + t*N

by putting P into the equation (1), I could solve t and find sphere center

Have you tried writing that equation down on paper? It becomes quite a long equation, I doubt these two lines solve it:

[source lang=cpp]
fl32 d = 2.0f*(math::Vec3Dot(N, A) - math::Vec3Dot(N, B));
fl32 t = ( math::Vec3Dot(A, A) - math::Vec3Dot(B, B) - 2.0f*S.x*(A.x+B.x) - 2.0f*S.y*(A.y+B.y) - 2.0f*S.z*(A.z+B.z) ) / d;[/source]
[font="Arial"]There should be a sqrt in the final equation.[/font]

Mussi, you're wrong:

||P-A||^2 = ||P-B||^2
P = S + tN

||S-A||^2 + t^2||N||^2 - 2t.dot(N,A) = ||S-B||^2 + t^2||N||^2 - 2t.dot(N,B)
2t.dot(N,A-B) = ||S-A||^2 - ||S-B||^2
t = (||S-A||^2 - ||S-B||^2) / (2.dot(N,A-B))

the error then, is that the OP has incorrectly equivalated:

||S-A||^2 - ||S-B||^2

with

||A||^2 - ||B||^2 - 2.dot(S,A+B)

whereas it would correctly be:

||A||^2 - ||B||^2 - 2.dot(S,A-B)

noting the sign change in the dot product from A+B to A-B

My vector algebra is a bit weak, what I was thinking was:

(P-A).(P-A) = (P-B).(P-B) =
(S+tN-A).(S+tN-A) = (S+tN-B).(S+tN-B) =
(Sx+t*Nx - Ax)[sup]2[/sup] + (Sy+t*Ny - Ay)[sup]2[/sup] + (Sz+t*Nz - Az)[sup]2[/sup] = (Sx+t*Nx - Bx)[sup]2[/sup] + (Sy+t*Ny - By)[sup]2[/sup] + (Sz+t*Nz - Bz)[sup]2[/sup] =

Hmm I notice now that (t*Nx)[sup]2[/sup], (t*Ny)[sup]2[/sup], (t*Nz)[sup]2[/sup] cancel out on both side of the equation, along with some other stuff. I guess you don't need a sqrt after all. Moving on:

Sx*Ax + t*Nx*Ax + Sy*Ay + t*Ny*Ay + Sz*Az + t*Nz*Az = Sx*Bx + t*Nx*Bx + Sy*By + t*Ny*By + Sz*Bz + t*Nz*Bz =
t = (Sx*Bx + Sy*By + Sz*Bz - Sx*Ax - Sy*Ay - Sz*Az) / (Nx*Ax + Ny*Ay + Nz*Az - Nx*Bx - Ny*By - Nz*Bz)

My head just started spinning, think I've missed something but it's getting too late .

I think I have it figured out.

First of all A and B need to be in the camera space, meaning that the line that goes out the center of the camera (the line formed by S + N * t) will fall directly on the z axis. This means that the X and Y of this line will always be 0 or.

Px = 0
Py = 0

This leaves us with only one variable to solve for or Pz.

We know that the distance between P and A is the same as P and B so

Ax^2 + Ay^2 + (Az - Pz)^2 = Bx^2 + By^2 + (Bz - Pz)^2

at this point we just solve for Pz.

Pz = ((Bx^2 + By^2 + Bz^2) - (Ax^2 + Ay^2 + Az^2)) / (2(Bz - Az))

The only thing left to do is the transform P into world space and calculate the radius of the circle by finding the distance between P and A or P and B.

If I didn't explain clearly enough let me know and I can give some detail.

Let E be the eyepoint. Let the view direction be D, the up direction be U, and the right direction be R. D, U, and R are unit length and mutually perpendicular. The frustum values are dmin (near), dmax (far), rmax (right), -rmax (left), umax (up), and -umax (down). The sphere center is C = E + s*D for some s. A vertex of the near plane is P0 = E + dmin*D + umax*U + rmax*R. A vertex of the far plane is P1 = E + dmax*D + (dmax*umax/dmin)*U + (dmax*rmax/dmin)*R. The squared distances from P0 to C and P1 to C must be equal, so (s-dmin)^2 + rmax^2 + umax^2 = (s - dmax)^2 + (dmax*rmax/dmin)^2 + (dmax*umax/dmin)^2. This leads to s = ((dmin + dmax)/2)*(1 + (rmax/dmin)^2 + (umax/dmin)^2).

This is actually the circumscribed sphere, which I believe is the minimum-volume sphere as long as center C is inside the frustum, which means as long as s <= dmax. When s > dmax, the minimum-volume sphere is the one whose center is E+dmax*D. The radius is distance from C to a vertex on the far plane. In the case s < dmax, the radius is the distance from C to any of the 8 vertices, but you can just substitute the value of s in either of the squared distances I mentioned and then take the square root.

@Mussi: I solved them using Maple, you should give it a try, nice app

@Owl: I have tried checking miniball algorithm out, Dave Eberly also has an implementation in his website, which is very complicated and as long as frustum has 8 verts and is symmetric I don't think I need to implement that.


@Dave, @HappyCoder
thanks guys, as long as I got HappyCoder's method more clearly, I'll try to implement that first and check the results

another method I'm currently using is that, I use barycentric minimum bounding spheres method (described here) for two tetrahedrons inside of the frustum. if indices 0~3 is vertices of the frustum's near plane, and 4~7 for far plane, I calculate one sphere for 0-2-5-7, and one for 1-3-4-6 (opposite tetrahedrons) , and merge two spheres. I think get acceptable results from this, but it's still pretty big especially for frustums with high aspect ratio (fov is high and (far-near) is small), although I don't think I can get small bounding spheres for those, because of the nature of the shape.

Did the circumcircle idea I suggested not work out?



//----------------------------------------------------------------------------
void ComputeMinimumVolumeSphereOfFrustum (const float eye[3],
const float direction[3], const float dmin, const float dmax,
const float rmax, const float umax, float center[3], float& radius)
{
float invDMin = 1.0f/dmin;
float u = umax*invDMin;
float r = rmax*invDMin;
float u2pr2 = u*u + r*r;
float s = 0.5f*(dmin + dmax)*(1.0f + u2pr2);
if (s >= dmax)
{
s = dmax;
radius = dmax*sqrt(u2pr2);
}
else
{
float diff = 1.0f - s/dmax;
radius = dmax*sqrt(diff*diff + u2pr2);
}
for (int i = 0; i < 3; ++i)
{
center = eye + s*direction;
}
}
//----------------------------------------------------------------------------
template <int n, typename Real>
void ComputeMiniballFrustum (const Real eye[n], const Real direction[n],
const Real min0, const Real max[n], Real center[n], Real& radius)
{
Real invMin0 = ((Real)1)/min0;
float ratio[n], sqrLength = (Real)0;
for (int i = 1; i < n; ++i)
{
ratio = max*invMin0;
sqrLength += ratio*ratio;
}
Real s = ((Real)0.5)*(min0 + max[0])*((Real)1 + sqrLength);
if (s >= max[0])
{
s = max[0];
radius = max[0]*sqrt(sqrLength);
}
else
{
Real diff = (Real)1 - s/max[0];
radius = max[0]*sqrt(diff*diff + sqrLength);
}
for (int i = 0; i < n; ++i)
{
center = eye + s*direction;
}
}
//----------------------------------------------------------------------------
int main ()
{
float eye[3] = { 0.0f, 0.0f, 0.0f };
float direction[3] = { 0.0f, 0.0f, 1.0f };
float dmin = 10.0f, dmax = 1000.0f, rmax = 2.0f, umax = 1.0f;
float center0[3], radius0;
ComputeMinimumVolumeSphereOfFrustum(eye, direction, dmin,
dmax, rmax, umax, center0, radius0);

float min0 = dmin, max[3] = { dmax, rmax, umax };
float center1[3], radius1;
ComputeMiniballFrustum<3,float>(eye, direction, min0, max, center1,
radius1);
return 0;
}
//----------------------------------------------------------------------------



Methods described by HappyCoder and David Eberly as well as the method I described in my previous post (slower), worked and gave identical results
thanks

@Postie: no I couldn't get that to work, I used barycentric coordinates to calculate circum circle of three points, and chose different set of points, but none of them gave minimum sphere of the frustum