Area and center location of an ellipse generated by the intersection of an ellipsoid and plane

Question

I am working on a model that requires that I know the area and center coordinates of the ellipse that is created by the intersection of an ellipsoid and a plane.

Specifically, for the location of the center of the ellipse I want to know the coordinates of this ellipse in Cartesian coordinates.

The general equations for the ellipsoid and plane I have started with are: $$\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1$$

$$m(x-x_0) + n(y-y_0) + k(z-z_0) = 0$$

I am writing the general forms of these equations because actually I need to be able to solve this for a number of different ellipses and planes with different orientations.

One specific case that I would like to use this solution for first is the case of an ellipsoid that is defined by:

$$a = 7, b = 5, c = 6$$

and a plane defined by:

$$(y+b)tan(\theta)-z + \frac{1}{2} = 0$$

Where $\theta$ is the desired angle of the plane, in this case $\theta=30^o$ is perfectly fine (arbitrary example).

Please note that "b" is the same b used in the equation of the ellipsoid.

I looked at some of the other threads that have asked about plane ellipsoid intersections. However, since I specifically need to calculate the area of the ellipse generated by this intersection and the location of its center, when I tried to use a parametric solution, I had difficulty doing this once I had the parametric equations.

I would really love to learn how to solve this problem so I can include it into my model. Help is greatly appreciated.

Thank you!

-Christian

Answer 1

Without loss of generality, you can translate the ellipsoid to origin, rotate it so that its semiaxes are parallel to the Cartesian coordinate axes, and reorder the axes so that the $z$ component of the intersection plane normal in the rotated coordinates has the largest magnitude. Then, the points on the ellipsoid fulfill $$\frac{x^2}{r_x^2} + \frac{y^2}{r_y^2} + \frac{z^2}{r_z^2} = 1 \tag{1}\label{AC1}$$ where $r_x$, $r_y$, and $r_z$ are the ellipsoid semiaxes.

Define the plane using its normal $(n_x, n_y, n_z)$ and signed distance from origin $n_d$, i.e. point $(x, y, z)$ is on the plane if and only if $$x n_x + y n_y + z n_z = n_d \tag{2a}\label{AC2a}$$ Solving this for $z$ yields $$z = \frac{n_d - x n_x - y n_y}{n_z} \tag{2b}\label{AC2b}$$ (We reorder the axes so that $\lvert n_z \rvert \ge \lvert n_y \rvert$, $\lvert n_z \rvert \ge \lvert n_x \rvert$, to avoid division by zero, and in numerical computations, for best numerical stability.)

Substituting $\eqref{AC2b}$ into $\eqref{AC1}$ yields general quadratic form $$a x^2 + b x y + c y^2 + d x + e y + f = 0$$ where (assuming $n_z r_z \ne 0$, both sides multiplied by $n_z^2 r_z^2$) $$\left\lbrace ~ \begin{aligned} a &= n_x^2 + n_z^2 r_z^2 / r_x^2 \\ b &= 2 n_x n_y \\ c &= n_y^2 + n_z^2 r_z^2 / r_y^2 \\ d &= -2 n_x n_d \\ e &= -2 n_y n_d \\ f &= n_d^2 - n_z^2 r_z^2 \\ \end{aligned} \right . \tag{3}\label{AC3}$$ Using this answer by Osmund Francis and the Wikipedia Ellipse article, we can express the following:

The discriminant $D$ (note sign convention!) is $$D = b^2 - 4 a c$$ where the intersection is an ellipse if and only if $D \lt 0$.

Center of the ellipse is at $(x_0, y_0)$, $$\left\lbrace ~ \begin{aligned} \displaystyle x_0 &= \frac{2 c d - b e}{D} \\ \displaystyle y_0 &= \frac{2 a e - b d}{D} \\ \end{aligned} \right.$$ Semimajor axis length $r_+$ and semiminor axis length $r_-$ are $$\begin{aligned} r_+ &= \frac{1}{-D}\sqrt{ 2 ( a e^2 + c d^2 - b d e + f D ) ( a + c + \sqrt{ b^2 + ( a - c )^2 } ) } \\ r_- &= \frac{1}{-D}\sqrt{ 2 ( a e^2 + c d^2 - b d e + f D ) ( a + c - \sqrt{ b^2 + ( a - c )^2 } ) } \\ \end{aligned}$$ and the area $A$ of the ellipse is $$A = \pi ~ r_+ ~ r_-$$ The angle $\theta$ between $x$ axis and major axis is $$\theta = \begin{cases} \operatorname{atan2}\left(c - a - \sqrt{b^2 + (a-c)^2}, b\right), & b \ne 0 \\ 0, & b = 0, \quad a \lt c \\ 90^o, & b = 0, \quad a \gt c \\ 0, & b = 0, \quad a = c, \quad \lvert d \rvert \ge \lvert f \rvert \\ 90^o, & b = 0, \quad a = c, \quad \lvert d \rvert \lt \lvert f \rvert \\ \end{cases}$$ Note that the major axis intersects the ellipse at $(x_{+1}, y_{+1})$ and $(x_{+2}, y_{+2})$, $$\left\lbrace ~ \begin{aligned} x_{+1} &= x_0 + r_+ ~ \cos\theta \\ y_{+1} &= y_0 + r_+ ~ \sin\theta \\ \end{aligned} \right ., \quad \left\lbrace ~ \begin{aligned} x_{+2} &= x_0 - r_+ ~ \cos\theta \\ y_{+2} &= y_0 - r_+ ~ \sin\theta \\ \end{aligned} \right .$$ and similarly the minor axis intersects the ellipse at $(x_{-1}, y_{-1})$ and $(x_{-2}, y_{-2})$, $$\left\lbrace ~ \begin{aligned} x_{-1} &= x_0 + r_- ~ \cos\theta \\ y_{-1} &= y_0 + r_- ~ \sin\theta \\ \end{aligned} \right ., \quad \left\lbrace ~ \begin{aligned} x_{-2} &= x_0 - r_- ~ \cos\theta \\ y_{-2} &= y_0 - r_- ~ \sin\theta \\ \end{aligned} \right .$$ Finally, using angle parameter $\varphi$, we can parametrize the ellipse as $\bigr(x(\varphi), y(\varphi)\bigr)$, $0 \le \varphi \le 360^o$, $$\left\lbrace ~ \begin{aligned} x(\varphi) &= x_0 + r_+ \cos\theta \cos\varphi - r_- \sin\theta \sin\varphi \\ y(\varphi) &= y_0 + r_+ \sin\theta \cos\varphi + r_- \cos\theta \sin\varphi \\ \end{aligned} \right .$$ which is just $x = r_+ \cos \varphi$, $y = r_- \sin\varphi$, rotated by $\theta$, and translated to $(x_0, y_0)$.

Answer 2

This is pretty straightforward if you work in homogeneous coordinates.

The general equation of a quadric surface can be written in the form $\mathbf X^TQ\mathbf X=0$, where $Q$ is a symmetric $4\times4$ matrix. Given a coordinate system (i.e., a parameterization) $\lambda\mathbf u+\mu\mathbf v+\tau\mathbf w$ of the plane, let $M = [\mathbf u\;\mathbf v\;\mathbf w]$ so that $\mathbf X=M\mathbf x$. Then we have $$\mathbf X^TQ\mathbf X = (M\mathbf x)^TQ(M\mathbf x) = \mathbf x^T(M^TQM)\mathbf x=0,$$ that is, the intersection is a conic with matrix $M^TQM$ relative to this coordinate system of the plane.

Applying this to your problem, we have $Q=\operatorname{diag}(1/a^2,1/b^2,1/c^2,-1)$. For the matrix $M$, we need a point on the plane and any two linearly independent vectors orthogonal to its normal, $(m,n,k)$. A pair of orthogonal unit vectors removes the need to scale the computed area at the end, but it’s not necessary. For a point on the plane, you have $(x_0,y_0,z_0)$, and for the two vectors parallel to the plane, you can take any two of $(0,k,-n)$, $(-k,0,m)$ and $(n,-m,0)$ that are linearly independent. For example, taking the last two of these vectors we would have $$M=\begin{bmatrix}-k&n&x_0\\0&-m&y_0\\m&0&z_0\\0&0&1\end{bmatrix}.$$

Once you have $C=M^TQM$, you can use standard formulas to find the center of the ellipse and compute its area. For instance, the last row of $C^{-1}$ or, equivalently, of $\operatorname{adj}(C)$ are the homogeneous coordinates of the conic’s center, if it has one. After dehomogenizing to obtain the center coordinates $(u,v)$, you can translate the origin to this point, which will give you a matrix of the form $$C'=\begin{bmatrix}a&b&0\\b&c&0\\0&0&f\end{bmatrix}.$$ Note that translation doesn’t affect the upper-left $2\times2$ submatrix, but it does replace the lower-right element of $C$ with $f=(u,v,1)C(u,v,1)^T$. If this represents an ellipse, which it should if you started with an ellipsoid, its area is equal to ${\pi\lvert f\rvert\over\sqrt{ac-b^2}}$. Since the coordinate system chosen for the plane is probably not orthonormal, to get the true area you need to multiply this by $\lVert\mathbf u\times\mathbf v\rVert$, where $\mathbf u$ and $\mathbf v$ are the two basis vectors that you chose when constructing $M$. To obtain the 3-D coordinates of the center point, just multiply its homogeneous coordinates by $M$.

You might also want to know if the plane even intersects the ellipsoid in more than one point in the first place. This can be accomplished at various stages in the process, but it’s fairly easy to check this early on before doing any of the other work. If the ellipsoid were a unit sphere, then we could just check that the distance of the plane to the origin was less than one, so we can apply to the plane the same transformation that maps our ellipsoid to the unit sphere and then check the distance of the transformed plane from the origin. Skipping ahead to the result, the plane will intersect the ellipsoid in a nontrivial ellipse when $$(mx_0+ny_0+kz_0)^2\le (am)^2+(bn)^2+(ck)^2.$$

Applying this to your example with $\theta=\pi/6$, we first check that $$\left(-5\cdot\frac1{\sqrt3}-\frac12\right)^2\le\left(5\cdot\frac1{\sqrt3}\right)^2+(-1\cdot 6)^2.$$ It is, so we have an ellipse. We can choose $$M = \begin{bmatrix}1&0&0\\0&1&-5\\0&\frac1{\sqrt3}&\frac12\\0&0&1\end{bmatrix}.$$ This turns out to be convenient because the basis vectors are orthogonal and the ellipse’s axes are aligned with them. We then obtain $$C = M^TQM = \begin{bmatrix}\frac1{49}&0&0\\0&\frac{133}{2700}&\frac1{72\sqrt3}-\frac15\\0&\frac1{72\sqrt3}-\frac15&\frac1{144}\end{bmatrix}.$$ The center of this ellipse works out to be $\left(0,-\frac5{266}(5\sqrt3-216)\right)$, which gives us $f = \frac1{532}(20\sqrt3-429)$. Since the basis vectors are orthogonal, the norm of their cross product is just the product of their norms, which is equal to $2/\sqrt3$. Multiplying everything out produces an area of $${15(429-20\sqrt3)\over19\sqrt{133}}\pi \approx 84.8112,$$ which matches the area computed by GeoGebra. Finally, the 3-D coordinates of the ellipse’s center are $\left(0,-\frac{25}{266}(10+\sqrt3),\frac{18}{133}(3+10\sqrt3)\right) \approx (0,-1.10264,2.75014)$, which also matches the center computed by GeoGebra.

We could’ve normalized both of the basis vectors up front when constructing $M$, but I don’t think that really saves much work in this example: doing so only changes the nonzero off-diagonal entries of $C$ and eliminates one multiplication at the end at the cost of one division up front, which is a wash.