1
$\begingroup$

I am trying to determine whether a line and an ellipsoid (in general, inc. rotations) intersect.

Following a similar argument to the Line-sphere intersection using vector notation I have constructed the general equation for an ellipsoid:

$$(\mathbf{x} - \mathbf{c})^T \mathbf{R}^T \mathbf{A} \mathbf{R} (\mathbf{x}-\mathbf{c}) = 1$$

where:

  • $\mathbf{x}$ is a point on the ellipsoid

  • $\mathbf{c}$ is the centre of the ellipsoid

  • $\mathbf{A}$ is a diagonal matrix where the diagonal elements are the inverse of the square of the length of the principal semi-axes (e.g. $\frac{1}{a^2}$, $\frac{1}{b^2}$ and $\frac{1}{c^2}$ if defining the "standard ellipsoid")

  • $\mathbf{R}$ is a rotation matrix (built up from $x, y$ and $z$ rotations: e.g. $R=R_z(\alpha)R_y(\beta)R_z(\gamma)$)


We define a line by:

$\mathbf{x} = \mathbf{o} + \lambda\mathbf{l}$

where

$\mathbf{x}$ is a point on the line

$\mathbf{o}$ is the origin of the line

$\lambda$ is the distance along the line

$\mathbf{l}$ is the direction of the line (unit vector)


To find the intersection we substitute the equation of the line into the equation of the ellipsoid to get:

$(\mathbf{o} + \lambda\mathbf{l}-\mathbf{c})^T\mathbf{R}^T\mathbf{A}\mathbf{R}(\mathbf{o} + \lambda\mathbf{l}-\mathbf{c})=1$

which needs to be solved for $\lambda$. Unfortunately that is where I get stuck. I can solve an example numerically but I'm not sure how to rearrange it.

My guess is this can be rearranged into a quadratic in $\lambda$ (e.g. $a^2\lambda + b\lambda + c =0$) which means that I would only need to check the term under the square root sign in the quadratic formula for whether the line and ellipsoid intersect.

$\endgroup$
  • 2
    $\begingroup$ If you can solve the problem for a sphere, why not just transform the ellipsoid to a sphere? Translate it to the origin, and then transform it to a sphere with a liner transformation. These transformations take the line to a line, so you have the line-sphere intersection problem. $\endgroup$ – saulspatz Jul 31 '19 at 13:20
  • $\begingroup$ @saulspatz thanks, that sounds like an excellent idea! Any hints on how to do the transformation? I get the translation, but not sure how to apply the transformation. My main problem is that the sphere example has no matrix in it (just vectors) which I can rearrange. The matrix causes me the problems! I assume that remains even after the transformation? $\endgroup$ – Sean Elvidge Jul 31 '19 at 17:01
  • $\begingroup$ I'm not sure I understand exactly what the problem is. Are you trying to solve a particular problem? Are you trying to write a computer program? $\endgroup$ – saulspatz Jul 31 '19 at 17:29
  • $\begingroup$ @saulspatz yup you got it in one. Trying to write some code to determine if a (given) line intersects with an ellipsoid (given the centre, principal semi-axes and rotation angles). Code works fine checking a line with a translated ellipsoid (with no rotations). Adding the rotations is what causes the problems and made me try and go back to the underlying maths problem. $\endgroup$ – Sean Elvidge Jul 31 '19 at 19:06
  • $\begingroup$ Rotate the line at the same time as you rotated the ellipsoid then work with the new line and the new ellipsoid. Maybe it's rotating the line that's confusing you. Translate the line to the origin (by subtracting $\mathbf{o}$), apply the same rotation to the line as you did to the ellipse, and then translate back. (Add $\mathbf{o}$). I'm sure you'll find it easier to just transform the problem that to come up with some messy formula. $\endgroup$ – saulspatz Jul 31 '19 at 19:22
1
$\begingroup$

[Collecting the discussion from the comments into an answer.]

You can certainly expand the expression that you’ve got and collect terms to get a quadratic equation in $\lambda$: Setting $\mathbf p = \mathbf o-\mathbf c$ and $\mathbf Q=\mathbf R^T\mathbf A\mathbf R$, you have $$(\lambda\mathbf l+\mathbf p)^T\mathbf Q(\lambda\mathbf l+\mathbf p) = (\mathbf l^T\mathbf Q\mathbf l)\lambda^2+2(\mathbf l^T\mathbf Q\mathbf p)\lambda +\mathbf p^T\mathbf Q\mathbf p=1,$$ which has real solutions when $(\mathbf l^T\mathbf Q\mathbf p)^2\ge(\mathbf l^T\mathbf Q\mathbf l)(\mathbf p^T\mathbf Q\mathbf p-1)$.

On the other hand, you’re essentially constructing the ellipsoid by scaling, rotating and translating the unit sphere. Since you’re not interested in the actual points of intersection, you could apply the inverse transformation to the line and then use the point-line distance formula to determine the line’s distance from the origin. If this is no greater than one, then there is at least one intersection.

If $\mathbf x_1$ and $\mathbf x_2$ are two points on a line, then the distance of the point $\mathbf x_0$ from this line is given by the formula (from Wolfram MathWorld) $${\lvert (\mathbf x_1-\mathbf x_0)\times(\mathbf x_2-\mathbf x_0)\rvert \over \lvert\mathbf x_2-\mathbf x_1\rvert}.$$ If $\mathbf x_0$ is the origin, this reduces to $${\lvert \mathbf x_2\times\mathbf x_1\rvert \over \lvert \mathbf x_2-\mathbf x_1\rvert}.\tag1$$ Now, since in your question you have $\mathbf R^T\mathbf A\mathbf R$ for the matrix of the ellipse, the matrix $\mathbf R$ already represents the inverse rotation. That is, it aligns the ellipsoid with the coordinate axes. (If instead you meant for it to represent a rotation that maps the standard basis onto the principal axes, then the ellipsoid’s matrix should be $\mathbf R\mathbf A\mathbf R^T$ instead and the inverse map below should use $\mathbf R^{-1}=\mathbf R^T$ instead of $\mathbf R$.) Setting $\mathbf S=\mathbf A^{1/2} = \operatorname{diag}(1/a,1/b,1/c)$, the required inverse map is therefore $\mathbf x\mapsto \mathbf S\mathbf R(\mathbf x-\mathbf c)$. A convenient choice for the two points on the line are the images of $\mathbf o$ and $\mathbf o+\mathbf l$. Plugging them into (1), we get $${\lvert \mathbf S\mathbf R(\mathbf o-\mathbf c)\times \mathbf S\mathbf R(\mathbf o+\mathbf l-\mathbf c) \rvert \over \lvert \mathbf S\mathbf R(\mathbf o-\mathbf c)-\mathbf S\mathbf R(\mathbf o+\mathbf l-\mathbf c) \rvert} = {\lvert \mathbf S\mathbf R(\mathbf o-\mathbf c)\times\mathbf S\mathbf R\,\mathbf l \rvert \over \lvert \mathbf S\mathbf R\,\mathbf l\rvert}.$$ Since we’re comparing a ratio of nonnegative numbers to $1$, we can square and rearrange to get $$\lvert\mathbf S\mathbf R(\mathbf o-\mathbf c)\times\mathbf S\mathbf R\,\mathbf l \rvert^2 \le \lvert \mathbf S\mathbf R\,\mathbf l\rvert^2.$$ This looks like it could be more efficient to compute that the quadratic equation’s discriminant.

We can even go a bit farther. Using $\lvert\mathbf u\times\mathbf v\rvert^2=\lvert\mathbf u\rvert^2\lvert\mathbf v\rvert^2-(\mathbf u\cdot\mathbf v)^2$, $$\lvert\mathbf S\mathbf R(\mathbf o-\mathbf c)\times\mathbf S\mathbf R\,\mathbf l\rvert^2=\lvert\mathbf S\mathbf R(\mathbf o-\mathbf c)\rvert^2\lvert\mathbf S\mathbf R\,\mathbf l\rvert^2-\left((\mathbf S\mathbf R(\mathbf o-\mathbf c))\cdot\mathbf S\mathbf R\,\mathbf l\right)^2$$ so $$\left(\lvert\mathbf S\mathbf R(\mathbf o-\mathbf c)\rvert^2-1\right) \lvert\mathbf S\mathbf R\,\mathbf l\rvert^2 \le \left((\mathbf S\mathbf R(\mathbf o-\mathbf c))\cdot\mathbf S\mathbf R\,\mathbf l\right)^2,$$ but $\mathbf Q=(\mathbf S\mathbf R)^T(\mathbf S\mathbf R)$, so this is just the original inequality derived from the quadratic equation, refactored.

Note that if you already have unit direction vectors for the principal axes, there’s no reason for the purposes of this calculation to decompose the rotation into three primitive rotations as in your question: $\mathbf R$ is simply the matrix with these vectors as its rows. Note also that since $\mathbf S$ is diagonal, you can implement multiplication of a vector by $\mathbf S$ with three multiplications.

| cite | improve this answer | |
$\endgroup$
1
$\begingroup$

Since $\bf o$ and $\bf c$ are known vectors, as to simplify the notation let's put $\bf o - \bf c = \bf v$.

Next consider that $\bf A$ and therefore $\bf B = {\bf R}^T \; \bf A \; \bf R$ are symmetric matrices, which means $$ {\bf a}^{\,T} \,{\bf B}\,{\bf b} = {\bf b}^{\,T} \,{\bf B}\,{\bf a} $$

So the equation you already achieved can simply be developed as: $$ \eqalign{ & 1 = \left( {{\bf v} + \lambda \,{\bf l}} \right)^{\,T} {\bf B}\left( {{\bf v} + \lambda \,{\bf l}} \right) = \cr & = \left( {{\bf v}^{\,T} + \lambda \,{\bf l}^{\,T} } \right)\left( {{\bf B}\,{\bf v} + \lambda \,{\bf B}\,{\bf l}} \right) = \cr & = \left( {{\bf v}^{\,T} {\bf B}\,{\bf v} + \lambda \,{\bf v}^{\,T} {\bf B}\,{\bf l}} \right) + \left( {\lambda \,{\bf l}^{\,T} {\bf B}\,{\bf v} + \lambda ^2 \,\,{\bf l}^{\,T} {\bf B}\,{\bf l}} \right) = \cr & \lambda ^2 \,\,\left( {{\bf l}^{\,T} {\bf B}\,{\bf l}} \right) + 2\lambda \,\left( {{\bf v}^{\,T} {\bf B}\,{\bf l}} \right) + \left( {{\bf v}^{\,T} {\bf B}\,{\bf v}} \right) \cr} $$

which is in fact a simple scalar quadratic equation in $\lambda$.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ @amd: thanks for signalling, amended $\endgroup$ – G Cab Aug 2 '19 at 20:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.