Ellipsoid and Plane Intersection Equation

Date: 04/20/2000 at 09:07:09
From: Martin Sheen
Subject: Solid Geometry (Ellipsoids)

I need the equation for the intersection of a plane surface with an 
ellipsoid. For example, if the ellipsoid has half-axes a, b, and c, 
and the plane is normal to the vector [i,j,k] and also passes through 
the point (i,j,k), what are the half-axes and orientation of the 
ellipse of intersection?

If (i,j,k) lies on an axis the answer is easy enough, but I'm stuck 
for the general case.

I'm not a student, but a physicist doing research in health care. I 
would be most grateful for any pointers.

Date: 04/20/2000 at 14:04:17
From: Doctor Rob
Subject: Re: Solid Geometry (Ellipsoids)

Thanks for writing to Ask Dr. Math, Martin.

The intersection you describe will be a closed curve in space. As 
such, it doesn't have a single equation. Single equations describe 
surfaces, not curves. You will either need two equations in x, y, and 
z, or else three equations for x, y, and z in terms of some parameter. 
Of course, you already have two equations: those of the ellipsoid and 
the plane:

     x^2/a^2 + y^2/b^2 + z^2/c^2 = 1
     (x-i)*i + (y-j)*j + (z-k)*k = 0

Apparently that is not what you want, however.

You can find the center of the ellipse using the following technique. 
Find the projection of the ellipse on the xy-plane by solving the 
second equation for z and substituting that into the first equation. 
Bring this equation into general form

     A*x^2 + 2*B*x*y + C*y^2 +2*D*x + 2*E*y + F = 0

by expanding and collecting like terms. Then if the center of the 
projection has coordinates (H,K), you can find them by solving the 
simultaneous equations

     A*H + B*K = -D
     B*H + C*K = -E

These will be the x- and y-coordinates of the center of the ellipse 
before projection. Now repeat this procedure with the projection of 
the ellipse on the xz-plane, by solving for y and substituting, and 
finding the center of that projection. It will have coordinates 
(x,z) = (H,L). Then the center of the ellipse in space will be at 
(H,K,L). If I did the calculations correctly,

     H = a^2*i*(i^2+j^2+k^2)/(a^2*i^2+b^2*j^2+c^2*k^2)
     K = b^2*j*(i^2+j^2+k^2)/(a^2*i^2+b^2*j^2+c^2*k^2)
     L = c^2*k*(i^2+j^2+k^2)/(a^2*i^2+b^2*j^2+c^2*k^2)

(a rather pleasing set of formulae, I think).

To find the direction and magnitude of the major axis, you can use the 
following. Let the square distance from a point (x,y,z) on the ellipse 
from the center (H,K,L) be

     f(x,y,z) = (x-H)^2 + (y-K)^2 + (z-L)^2

You want to maximize this value subject to the constraints of the two 
equations of the ellipsoid and plane. To do this, solve the ellipsoid 
equation for z and substitute that in the last equation above, and the 
plane equation. Then solve the plane equation for y, and substitute 
that into the above equation. This will eliminate both y and z, so 
that f is a function of x alone (but a rather complicated one). Now 
find its derivative df/dx, set that equal to zero, and solve for x. 
You should get four values of x. Substituting these into earlier 
equations, you can solve for the corresponding values of y and z. 
These four points will be the endpoints of the major and minor axes of 
the ellipse in space. Then it is a simple matter to obtain the lengths 
and directions of the half-axes.

To do all this symbolically is quite a mess, but with actual numbers 
this should be quite feasible. Try it and see.

I hope this satisfies your request. If not, write again.

- Doctor Rob, The Math Forum
  http://mathforum.org/dr.math/   

Date: 05/02/2000 at 12:41:55
From: Martin Sheen
Subject: Re: Solid Geometry (Ellipsoids)

Very many thanks for your extremely prompt reply.

You did invite me to write again; may I trespass on your time further?

Sorry; I wasn't clear quite about what I wanted - I was trying to be 
brief. I was actually thinking of the equation of the two-dimensional 
ellipse in the coordinate frame of the intersecting plane. Easy enough 
to transform, of course.

What I'm doing is viewing a section through an ellipsoidal eye from a 
viewpoint normal to the intersecting plane, and displaying the 
intersection on that plane along with a projection of the eye 
structures and isodose lines for planning radio-therapy treatments. 
From your help I've got the ellipse center (H,K,L) okay (in eye 
coordinates), but I'm a bit stuck with the half-axes.

Solving the ellipsoid equation for z and substituting into f(x,y,z) 
gives an f(x,y) with terms not only in x^2 and y^2 but also a 
Sqrt[1 - (x/a)^2 - (y/b)^2], and substituting into the plane equation 
will give terms in x and y and the same Sqrt[...]. If we get rid of 
the square roots we'll get simultaneous equations with fourth-order 
and second-order terms respectively, and I don't know how to cope with 
that. If instead we solve the plane equation for z and substitute that 
into f(x,y,z) and the ellipsoid equations, we get two simultaneous 
equations of the form A*x^2 + 2*B* x*y + C*y^2 + ..., which may be 
simpler but I'm still stuck.

So; sorry to bother you again, but is there an alternative? The 
difference between the ellipsoid axes is not going to be extreme, and 
absolute accuracy is not required, so if there's a simpler 
approximation assuming only first-order difference from a sphere, that 
would probably do.

Very many thanks again for your time.
Martin Sheen

Date: 05/02/2000 at 16:32:39
From: Doctor Rob
Subject: RE: Solid Geometry (Ellipsoids)

Thanks for the further explanation. You can rotate coordinates so that 
the x'y'-plane is parallel to the plane in question. This you can do 
as follows.

Let (i,j,k) be direction numbers for the z' axis in the xyz-system. 
Then (-k,0,i) is a perpendicular direction (rather arbitrarily chosen 
among all such), which we'll make the direction numbers for the x' 
axis. Then we can normalize these to direction cosines

     lambda3 = i/sqrt(i^2+j^2+k^2)
     mu3     = j/sqrt(i^2+j^2+k^2)
     nu3     = k/sqrt(i^2+j^2+k^2)

     lambda1 = -k/sqrt(i^2+k^2)
     mu1     = 0
     nu1     = i/sqrt(i^2+k^2)

To find the direction cosines of the y' axis, they are

     lambda2 = i*j/[sqrt(i^2+k^2)*sqrt{i^2+j^2+k^2)]
     mu2     = -(i^2+k^2)/[sqrt(i^2+k^2)*sqrt{i^2+j^2+k^2)]
     nu2     = j*k/[sqrt(i^2+k^2)*sqrt{i^2+j^2+k^2)]

You can verify that the three direction vectors are mutually 
orthogonal unit vectors, so will represent a rotation of coordinates 
which will not change distances (an *isometry* or *orthogonal 
transformation*). Now make the substitution

     x = lambda1*x' + lambda2*y' + lambda3*z'
     y =     mu1*x' +     mu2*y' +     mu3*z'
     z =     nu1*x' +     nu2*y' +     nu3*z'

This will transform the equation of the ellipsoid into a rather 
general second-degree equation in x, y, and z, and the equation of the 
plane will be transformed into

     z' = sqrt(i^2+y^2+z^2)

Then you can substitute this value for z' in the equation of the 
ellipsoid. This will yield the x'y'-equation of the ellipse in the 
general form

     A*x'^2 + 2*B*x'*y' + C*y'^2 + 2*D*x' + 2*E*y' + F = 0

Next you'll want to rotate these coordinates in order to remove the 
x'*y' term. This is done by determining any angle t such that

     tan(2*t) = 2*B/(A-C)

Then the transformation is

     x' = cos(t)*x" - sin(t)*y"
     y' = sin(t)*x" + cos(t)*y"
     z' = z"

Now it is just a matter of completing the square on x" and y", moving 
the constant term to the other side of the equation, and dividing by 
it to bring the equation of the ellipse into the standard form

     (x"-h')^2/a'^2 + (y"-k')^2/b'^2 = 1

Then you can read off the semi-major and semi-minor axis lengths a' 
and b' from this form, and the center (h',k',sqrt(i^2+j^2+k^2)) in the 
x"y"z"-coordinate system. Finally you can transform this back into the 
old coordinate system using the substitution equations given above, if 
desired. Alternatively, you can translate the coordinate system via 
the substitution

     X = x" - h'
     Y = y" - k'
     Z = z" - sqrt(i^2+j^2+k^2)

and the equation of the ellipse and the plane become

     X^2/a'^2 + Y^2/b'^2 = 1
     Z = 0

I am not sure which of the coordinate systems above may have been what 
you had in mind when you wrote, "the equation of the two-dimensional 
ellipse in the coordinate frame of the intersecting plane." You pick.

If you would like me to work an example, I can oblige.

- Doctor Rob, The Math Forum
  http://mathforum.org/dr.math/   

Date: 05/08/2000 at 09:20:27
From: Martin Sheen
Subject: RE: Solid Geometry (Ellipsoids)

Again, I'm most grateful for your time and effort. This is just what I 
need.

Coding it without any mistakes is another matter, of course, but 
that's a familiar problem.

All the best,

Martin Sheen