I would like to be able to find two ellipsoid for collision. I have found a paper which describes an algorithm in which they use Lagrange multipliers and Cramer's rule to find the collision. However I end up with a 6th degree equation which has to be solved numerically.
I would like to optimize this by first checking for collision and if there is a collision I do that expensive calculation. My idea was to project the ellipsoids onto x, y and z axis and if there is an overlap between all projection pairs - there is a collision between the ellipsoids as well.
How computationally expensive would be such an approach and how I could do it?
Thanks a lot!
This is not an answer, but a rough outline on a possible iterative numerical method to find the minimum distance between two ellipsoids. The idea is to use two vectors, each on the surface of an ellipsoid, and connect them with an attractive "force".
Note that if $r_i$ and $R_i$ are the axes of the two ellipsoids, $i = 1, 2, 3$, then if the distance between the ellipsoid centers is $\left(\min(r_i) + \min(R_i)\right)$ or less, the two ellipsoids must intersect. (This is because a sphere of radius $\min(r_i)$ at the center of the first ellipsoid is contained within the first ellipsoid, a sphere of radius $\min(R_i)$ at the center of the second ellipsoid is contained within the second ellipsoid, and since those two spheres intersect, the ellipsoids must also intersect.)
Also note that if the distance between the ellipsoid centers is greater than $\left(\max(r_i) + \max(R_i)\right)$, the two ellipsoids cannot intersect. (This is because a sphere of radius $\max(r_i)$ contains the first ellipsoid, a sphere of radius $\max(R_i)$ the second ellipsoid, and since the spheres do not intersect, the ellipsoids cannot intersect either.)
The interesting case I explore here is when we have two non-degenerate (all $r_i \gt 0$, $R_i \gt 0$) ellipsoids with the distance between their centers between $\left(\min(r_i) + \min(R_i)\right)$ and $\left(\max(r_i) + \max(R_i)\right)$.
The reason I approach the issue with a function that finds the minimum distance between two ellipsoids, is that when that distance is known at some point, we know that the two ellipsoids cannot collide unless their orientation changes, or the sum of the distances their centers have moved since the last check is at least the distance found during that last check. If the ellipsoids typically stay in the same orientation, and only their centers move, this can potentially save a lot of unnecessary computation.
Let's assume each ellipsoid is defined by its center, $\vec{c}$, and an orthogonal matrix $\mathbf{r}$ whose column vectors are the three axis vectors $\vec{r}_1$, $\vec{r}_2$, and $\vec{r}_3$.
Note that the axis vectors must be orthogonal, $$\vec{r}_1 \cdot \vec{r}_2 = 0 , \; \; \vec{r}_1 \cdot \vec{r}_3 = 0 , \; \; \vec{r}_2 \cdot \vec{r}_3 = 0$$ and have nonzero length, $$\vec{r}_1 \cdot \vec{r}_1 \gt 0, \; \; \vec{r}_2 \cdot \vec{r}_2 \gt 0, \; \; \vec{r}_3 \cdot \vec{r}_3 \gt 0$$
The function that describes the ellipsoid with respect to its center is $$f(\vec{s}) = \left(\frac{\vec{s}\cdot\vec{r}_1}{\vec{r}_1\cdot\vec{r}_1}\right)^2 + \left(\frac{\vec{s}\cdot\vec{r}_2}{\vec{r}_2\cdot\vec{r}_2}\right)^2 + \left(\frac{\vec{s}\cdot\vec{r}_3}{\vec{r}_3\cdot\vec{r}_3}\right)^2 = 1$$ so the normal at any surface point $\vec{s}$ (relative to its center) is $$\nabla f(\vec{s}) = \vec{n}(\vec{s}) = 2 \vec{r}_1 \frac{\vec{s}\cdot\vec{r}_1}{\left(\vec{r}_1\cdot\vec{r}_1\right)^2} + 2 \vec{r}_2 \frac{\vec{s}\cdot\vec{r}_2}{\left(\vec{r}_2\cdot\vec{r}_2\right)^2} + 2 \vec{r}_3 \frac{\vec{s}\cdot\vec{r}_3}{\left(\vec{r}_3\cdot\vec{r}_3\right)^2}$$
To speed up the test whether point $\vec{s}$ relative to the center of the ellipsoid is contained within the ellipsoid ($f\left(\vec{s}\right) \le 1$), we can precalculate a matrix $\mathbf{b}$, whose three row vectors $\vec{b}^T_1$, $\vec{b}^T_2$, $\vec{b}^T_3$ only depend on $\mathbf{r}$: $$\vec{b}^T_1 = \frac{\vec{r}_1}{\vec{r}_1\cdot\vec{r}_1} , \; \; \vec{b}^T_2 = \frac{\vec{r}_2}{\vec{r}_2\cdot\vec{r}_2} , \; \; \vec{b}^T_3 = \frac{\vec{r}_3}{\vec{r}_3\cdot\vec{r}_3}$$ Essentially, each $\vec{b}^T_i$ is parallel to $\vec{r}_i$, but with inverted norm: $\lvert\vec{b}^T_i\rvert = 1 / \lvert\vec{r}_i\rvert$. The test is then $$\vec{v} \cdot \vec{v} \le 1 \; \text{where} \; \vec{v} = \mathbf{b} \vec{s} = \left( \vec{s} \cdot \vec{b}^T_1 ,\, \vec{s} \cdot \vec{b}^T_2 ,\, \vec{s} \cdot \vec{b}^T_3 \right)$$ i.e. only costs four dot products and one (scalar) comparison; or twelve scalar multiplications, eight scalar additions, and one (scalar) comparison.
If we have some point $\vec{p}$ with respect to the center of the ellipsoid (and not too close to the center of ellipsoid, say $\lvert\vec{p}\rvert \ge \min(\lvert\vec{r}_1\rvert, \lvert\vec{r}_2\rvert, \lvert\vec{r}_3\rvert)$), we can find the point $\vec{s}$ on the surface of the ellipsoid that is on the line between the center of the ellipsoid and $\vec{p}$ (on the same side as $\vec{p}$), using an unit vector $\hat{u}$, $\hat{u}\cdot\hat{u} = 1$: $$\hat{u} = \frac{\vec{u}}{\sqrt{\vec{u}\cdot\vec{u}}} , \; \; \vec{u} = \mathbf{b} \vec{p}$$
Side note: I believe $\mathbf{b}$ is an affine transform that deforms (reforms?) the ellipsoid back into a sphere, but I'm not a mathematician, so I might be wrong. However, I do have verified numerically that it works (to quite a high precision) -- specifically, that $\mathbf{r}\hat{u}$ is on the surface of the ellipsoid and parallel to $\vec{p}$ --, with a few thousand random ellipsoids and a few hundred thousand random test points in each.
Initially, start at the ellipse surface points on the line connecting the two centers. Check if the points are contained within the other ellipse.
On each iteration step, let $\vec{s}$ is the surface point on the current ellipse, relative to the center of the current ellipse, and $\vec{a}$ is the surface point on the other ellipse, relative to the center of the current ellipse. Calculate the vector difference between $\vec{a}$ (scaled to unit length) and the ellipsoid surface unit normal at $\vec{s}$: $$\vec{d} = \frac{\vec{a}}{\sqrt{\vec{a}\cdot\vec{a}}} - \frac{\vec{n}(\vec{s})}{\sqrt{\vec{n}(\vec{s})\cdot\vec{n}(\vec{s})}}$$
If $\vec{d}$ is very short for both ellipsoids, $\vec{s}$ should be the nearest points between the two ellipsoids, and you can return their absolute distance (as the distance between the ellipsoids).
The "attractive force" "pulls" the current surface point $\vec{s}$ in the direction of $\vec{d}$. Calculate $$\vec{u} = \mathbf{b}(\vec{s} + \lambda \vec{d})$$ where $\lambda$ is a small positive factor (that probably should decrease during iteration), the amount of movement in direction $\vec{d}$; then, $$\hat{u} = \frac{\vec{u}}{\sqrt{\vec{u}\cdot\vec{u}}}$$ and finally $$\vec{s} = \mathbf{r}\hat{u}$$
If $\vec{s}$ is contained in the other ellipsoid, the two ellipsoids intersect. (Obviously, you should check the $\vec{s}$ of both ellipsoids.)
Otherwise, I would calculate the squared distance between the $\vec{s}$'s in absolute coordinates, and remember the minimum value, so that if some predefined iteration limit is reached, the function can return an estimate (the square root of the remembered minimum).
Does this work? I don't know; I haven't actually coded it yet.
Is it faster than numerically solving a sixth-degree function? I don't know; I don't even know what that sixth-degree function looks like. (I suspect that computing the coefficients is actually more computationally expensive than finding the roots of a sixth-degree polynom.)
The above can also be written using two unit vectors, $\hat{u}_A$ and $\hat{u}_B$, with the intent of minimizing $\lvert \vec{c}_{AB} + \mathbf{b}_A \hat{u}_A - \mathbf{b}_B \hat{u}_B \rvert$, where $\vec{c}_{AB} = \vec{c}_B - \vec{c}_A$.
One could even precalculate some $N$ unit vectors $\hat{u}_i$, and the $N$ points on an ellipsoid $\vec{s}_i = \mathbf{b}\hat{u}_i$ whenever the axes of the ellipsoid change (i.e. size or orientation changes; requiring $9N$ multiplications and $6N$ additions per ellipsoid); then the $O(N^2)$ search for $\min\lvert \vec{c}_{AB} + \vec{s}^A_i - \vec{s}^B_j \rvert$ would only require one square root, $3N^2$ multiplications, and $4N^2$ additions or subtractions.