Shadow of an ellipsoid - the ultimate version

If $(x', y', z')$ are the coordinates of a point, on the ellipsoid, with respect to $\mathbf{u_1}, \mathbf{u_2}, \mathbf{u_3}$ , then we know that,

$\dfrac{ {x'}^2 }{15^2} + \dfrac{ {y'}^2 }{30^2} +\dfrac{ {z'}^2 }{10^2} = 1$

Thus if we define $\mathbf{r'} = (x', y', z')$ , and $\mathbf{r} = (x, y, z )$ ,then

${\mathbf{r'}}^T D_e \mathbf{r'} = 1$ ,

where $D_e$ is a diagonal matrix whose diagonal entries are, in order, $\dfrac{1}{15^2}$ , $\dfrac{1}{30^2}$ and $\dfrac{1}{10^2}$ . The relationship between $\mathbf{r'}$ and $\mathbf{r}$ is simply,

$\mathbf{r} = R_e \mathbf{r'}$ where $R_e = [ \mathbf{u_1}, \mathbf{u_2}, \mathbf{u_3} ]$

The three columns of $R_e$ are orthonormal, hence ${R_e}^{-1} = {R_e}^T$ , and thus, the equation of the ellipsoid is,

${\mathbf{r}}^T Q_e \mathbf{r} = 1$ , where $Q_e = R_e D_e {R_e}^T$

It is a simple matter to find the lowest point on the ellipsoid, because at that point the normal vector to the surface of the ellipsoid is pointing in the negative $\mathbf{k}$ direction, where $\mathbf{k} = (0, 0, 1)$ . Now, the normal vector is parallel to $Q_e r$ , hence, we have, $Q_e r = - \alpha \mathbf{k}$ for some positive constant $\alpha$ . Solving for $\mathbf{r}$ and plugging this into the ellipsoid equation, one finds that the lowest point $\mathbf{r}^*$ on the surface of the ellipsoid is given by

$\mathbf{r}^* = -\dfrac{1}{\sqrt{ \mathbf{k}^T {Q_e}^{-1} \mathbf{k} }} {Q_e}^{-1} \mathbf{k}$

Shifting the ellipsoid, such that the above found point becomes at the origin, one finds that the shift vector is $\mathbf{r}_0 = - \mathbf{r}^*$ . The shifted image of a point $\mathbf{r}$ is the point $\mathbf{r'} = \mathbf{r} + \mathbf{r_0}$ . Substituting for $\mathbf{r}$ into the ellispoid equation, we end up with,

$(\mathbf{r'}-\mathbf{r_0})^T Q_e (\mathbf{r'} - \mathbf{r_0}) = 1$

Dropping the primes off the primed variables, we end up with,

$(\mathbf{r}-\mathbf{r_0})^T Q_e (\mathbf{r} - \mathbf{r_0}) = 1$

And this is the final equation of ellipsoid that is now touching the $xy$ plane at the origin.

The shadow of the ellipsoid on the $xy$ plane is generated by the intersection of the $xy$ plane and the non-uniformly scaled cone that is tangent to the ellipsoid, and which has its apex at the point ligth source $p_0$ . Tangency points on the ellipsoid are characterized by the vector equation,

$(\mathbf{r} - \mathbf{p_0})^T Q_e (\mathbf{r} - \mathbf{r_0} ) = 0$

which can be manipulated as follows,

$( \mathbf{r} - \mathbf{r_0} + \mathbf{r_0} - \mathbf{p_0})^T Q_e (\mathbf{r} - \mathbf{r_0} ) = 0$

And this simplifies to,

$(\mathbf{r_0} - \mathbf{p_0})^T Q_e (\mathbf{r} - \mathbf{r_0}) = -1$

This is an equation of a plane, thus all tangent points lie in a plane. Intersecting a plane with an ellipsoid was addressed here . The resulting ellipse of intersection defines the base of the general (non-uniformly scaled) cone, which has its apex at $p_0$ . The next and final task is to intersect this cone with $xy$ plane, as was addressed in this problem . And this results in the shadow ellipse.

This is one way to proceed. The other way we can go is by transforming the ellipsoid into the unit sphere, then relating the original general cone to the regular cone resulting from the transformation. Recall that

$(\mathbf{r}-\mathbf{r_0})^T Q_e (\mathbf{r} - \mathbf{r_0}) = 1$

And that $Q_e = R_e D_e {R_e}^T = R_e \sqrt{D_e} \sqrt{D_e} {R_e}^T$

Thus, if define $\mathbf{r'} = \sqrt{D_e} {R_e}^T (\mathbf{r} - \mathbf{r_0}) = A (\mathbf{r} - \mathbf{r_0})$ , then

${\mathbf{r'}}^T {\mathbf{r'}} = 1$ , i.e. $\mathbf{r'}$ is on the unit sphere centered at the origin. The corresponding

position of the point ligth source is $\mathbf{p'_0} = A (\mathbf{p_0} - \mathbf{r_0} )$ . Our general cone is transformed

into a regular right circular cone, with apex $\mathbf{p'_0}$ , and whose axis is along the vector $\mathbf{p'_0}$ .

The semi-vertical angle of this right circular cone is simply $\theta_c = \sin^{-1} \dfrac{1}{| \mathbf{p'_0} | }$ . Therefore, the

equation of this cone is $(\mathbf{r'} - \mathbf{p'_0})^T R_c D_c {R_c}^T (\mathbf{r'} - \mathbf{p'_0} ) = 0$

where $D_c$ is given by,

$D_c = \begin{bmatrix} \cos^2 \theta_c && 0 && 0 \\ 0 && \cos^2 \theta_c && 0 \\ 0 && 0 && -\sin^2 \theta_c \end{bmatrix}$

And the rotation matrix is of the form

$R_c = [ \mathbf{w_1} , \mathbf{w_2}, \mathbf{w_3} ]$

such that $\mathbf{w_3}$ is a unit vector in the direction of $\mathbf{p'_0}$ (the axis of the right circular cone) and $\mathbf{w_1}$ , $\mathbf{w_2}$

are mutually orthogonal unit vectors that are also orthogonal to $\mathbf{w_3}$ . Define $Q'_c = R_c D_c {R_c}^T$ ,

then we can, in a straightforward manner, find the corresponding matrix, in the original coordinates $\mathbf{r}$

We start with,

$(\mathbf{r'} - \mathbf{p'_0})^T Q'_c (\mathbf{r'} - \mathbf{p'_0} ) = 0$

And substitute

$\mathbf{r'} = A (\mathbf{r} - \mathbf{r_0})$

Hence, the equation of the original cone is (recall that $\mathbf{p'_0} = A(\mathbf{p_0} - \mathbf{r_0} )$ ),

$(\mathbf{r} - \mathbf{p_0} )^T Q_c (\mathbf{r} - \mathbf{p_0} ) = 0$

where $Q_c = A^T Q'_c A$ . Now that we have the equation of the original (general, untransformed) cone, we intersect it with the $xy$ plane, using the procedure outlined in the solution of this problem , and is repeated below for the reader's convenience.

The vector equation of any plane is

$\mathbf{r} = \mathbf{r_1} + W \mathbf{u}$

where $\mathbf{r_1}$ is any arbitrary point on the plane, and $W = [ \mathbf{w_1} , \mathbf{w_2} ]$ is a $3 \times 2$ matrix whose columns are arbitrary unit vectors that together with the normal to the plane form an orthonormal basis for $\mathbb{R}^3$ .

In the case of the $xy$ -plane, it is natural to choose $\mathbf{r_1} = (0, 0, 0)$ , $\mathbf{w_1} = (1, 0, 0)$ and $\mathbf{w_2} = (0, 1, 0)$ .

Substituting this vector expression into the equation of the oblique cone, one obtains,

$( \mathbf{r_1} - \mathbf{p_0} + W \mathbf{u} )^T Q_c (\mathbf{r_1} - \mathbf{p_0} + W \mathbf{u}) = 0$

Expanding this quadratic form gives,

$\mathbf{u}^T (W^T Q_c W) \mathbf{u} + 2 \mathbf{u}^T W^T Q_c (\mathbf{r_1} - \mathbf{p_0}) + (\mathbf{r_1} - \mathbf{p_0})^T Q_c ( \mathbf{r_1} - \mathbf{p_0}) = 0$

Take $\mathbf{u_0} = - (W^T Q_c W)^{-1} W^T Q_c (\mathbf{r_1} - \mathbf{p_0} )$ , then the above equation becomes,

$(\mathbf{u} - \mathbf{u_0})^T (W^T Q_c W) (\mathbf{u} - \mathbf{u_0}) = C$

where the constant $C$ is given by $C = - (\mathbf{r_1} - \mathbf{p_0})^T Q_c (\mathbf{r_1} - \mathbf{p_0}) + {\mathbf{u_0}}^T (W^T Q_c W) \mathbf{u_0}$

Dividing by $C$ , we obtain,

$(\mathbf{u} - \mathbf{u_0})^T Q_1 (\mathbf{u} - \mathbf{u_0}) = 1$

where $Q_1 = \dfrac{1}{C} (W^T Q_c W)$ .

The final step is to diagonalize $Q_1$ into $Q_1 = R D R^T$ . The diagonal elements of the diagonal matrix $D$ are the eigenvalues of $Q_1$ which are the reciprocal of the square of the semi-major and semi-minorr axes lengths $a$ and $b$ of the resulting ellipse of intersection.

The area is given by $\pi a b$ .

Here's an outline of my solution. It was improvised, but ultimately effective. I'm sure there are more formal methods.

1) Find the point on the ellipsoid with the minimum z coordinate. Suppose I generate an arbitrary direction vector $\vec{v} = (v_x,v_y,v_z)$ . Starting from the center of the ellipsoid, a displacement of $\sigma \vec{v}$ will get us to the ellipsoid surface. What is the value of $\sigma$ ? To find out, first represent $\vec{v}$ as a linear combination of the axial vectors:

$\large{\vec{v} = \alpha \, \vec{u_1} + \beta \, \vec{u_2} + \gamma \, \vec{u_3} \\ v_x = \alpha \, u_{1x} + \beta \, u_{2x} + \gamma \, u_{3x} \\ v_y = \alpha \, u_{1y} + \beta \, u_{2y} + \gamma \, u_{3y} \\ v_z = \alpha \, u_{1z} + \beta \, u_{2z} + \gamma \, u_{3z}}$

We can solve for $\alpha, \beta, \gamma$ by solving the linear system. This solution is also very useful when transforming coordinates for the shadow calculation. Now consider the vector $\sigma \vec{v}$ .

$\large{\sigma \, \vec{v} = \sigma \, \alpha \, \vec{u_1} + \sigma \, \beta \, \vec{u_2} + \sigma \, \gamma \, \vec{u_3}}$

Suppose that $\sigma$ is such that we are displaced from the center to the ellipsoid surface. ( $a,b,c)$ are the semi-lengths:

$\large{\frac{\sigma^2 \alpha^2}{a^2} + \frac{\sigma^2 \beta^2}{b^2} + \frac{\sigma^2 \gamma^2}{c^2} = 1 \\ \sigma = \frac{1}{\sqrt{\frac{ \alpha^2}{a^2} + \frac{ \beta^2}{b^2} + \frac{ \gamma^2}{c^2}}}}$

Now we have determined the ellipsoid surface point corresponding to displacement from the center along an arbitrary direction vector. I used a hill-climbing algorithm, varying the components of $\vec{v}$ , finding the corresponding ellipsoid surface point, and storing the coordinates with the minimum z value. It is then a simple matter to determine the center of the ellipsoid, such that the point now lies at the origin

2) The next step is to determine the area of the shadow. Imagine sweeping over a rectangular grid in the $xy$ plane, and drawing a line from each evaluated $xy$ point to the light source. That line will intersect the ellipsoid zero, one, or two times. Zero times corresponds to bright (lit) area outside the shadow. One time corresponds to the boundary between light and shadow. Two times corresponds to shadow.

Let us mathematically examine the two-intersection case. To do this, take the center of the ellipsoid (found earlier) as the origin of the "prime" coordinate system, and use the semi-axial vectors as the basis vectors for the "prime" coordinate system. Represent the prime-referenced position of a point in the $xy$ plane (with z = 0) as $(x',y',z')$ . Represent the prime-referenced light source coordinates as $(P'x, P'y, P'z)$ . These transformations are accomplished in the manner described earlier (solution of the linear system for the three basis vector weights).

Suppose we start at $(x',y',z')$ and travel along a displacement vector $\vec{q}$ (also in "prime" coordinates) on the way to $(P'x, P'y, P'z)$ . The following equation describes the intersection. Note that the value $\alpha$ is used in a different context here than it was earlier.

$\large{\frac{(x' + \alpha \, q_x)^2}{a^2} + \frac{(y' + \alpha \, q_y)^2}{b^2} + \frac{(z' + \alpha \, q_z)^2}{c^2} = 1}$

Re-arranging into a form more suitable to the quadratic equation yields:

$\large{A \, \alpha^2 + B \, \alpha + C = 0 \\ A = \frac{q_x^2}{a^2} + \frac{q_y^2}{b^2} + \frac{q_z^2}{c^2} \\ B = \frac{2 \, x' \, q_x}{a^2} + \frac{2 \, y' \, q_y}{b^2} + \frac{2 \, z' \, q_z}{c^2} \\ C = \frac{x'^2}{a^2} + \frac{y'^2}{b^2} + \frac{z'^2}{c^2} - 1}$

Shadowed points are points with two solutions for $\alpha$ , meaning $B^2 > 4 A C$ . So we scan through the plane in normal $xyz$ space (with z = 0) in a high-resolution rectangular grid of sufficient size, transform the coordinates appropriately into the "prime" reference frame, calculate the number of ray intersections with the ellipsoid for each scan point, and add an incremental area $\Delta x \, \Delta y$ to a running area sum for each shadowed point. The total shadow area comes out to $1558.9$ .