Is there a function to compute the exact (or approximate) shared volume between two ellipsoids with arbitrary rotation and translation?
-
Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. – Community Aug 29 '21 at 14:48
-
What kind of information I should provide? The question seems self-explanatory – i_am_juan Aug 29 '21 at 14:51
-
Depending on your context and requirements I might have had the same question. See my answer below how I dealt with it. – rohlemax Jun 28 '22 at 18:43
-
The analogue for ellipses (2 D) is still complicated. If the principal axes are parallel and proportional then it's like the problem for spheres. Otherwise, you can at least reduce to the case when one of the ellipsoids is a sphere. I am not sure it makes it much easier. Interesting question! – orangeskid Jun 28 '22 at 19:04
2 Answers
I don't know if this problem can be solved analytically. However it is quite feasible for numerical solution.
Suppose the first ellipsoid has the equation
$( r - C_1)^T Q_1 ( r- C_1) = 1 $
where $r$ is any point (vector) on the surface of the ellipsoid, $C_1$ is the center of the ellipsoid and $Q_1$ is a symmetric positive definite matrix. Thus, $Q_1$ can be decomposed into $Q_1 = R_1 D_1 R_1^T $
The matrix $D_1$ is diagonal, and its diagonal elements are the reciprocals of the square the semi-axes lengths of the ellipsoid. That is,
$D_1 = \begin{bmatrix} \dfrac{1}{a_1^2} && 0 && 0 \\ 0 && \dfrac{1}{b_1^2} && 0 \\ 0 && 0 && \dfrac{1}{c_1}^2 \end{bmatrix} $
Matrix $R_1$ is an orthogonal matrix whose columns are the unit vectors that are parallel to the axes of the ellipsoid.
A similar structure applies for the second ellipsoid. That is, its equation is given by:
$(r - r_2)^T Q_2 (r - r_2) = 1$
with $Q_2 = R_2 D_2 R_2^T $
and with $D_2$ containing the reciprocals of the semi-axes lengths for the second ellipsoid, and $R_2$ is orthogonal and has it columns parallel to the axes of the second ellipsoid.
The above described the algebraic equation of any two ellipsoid. An equivalent representation of the surface of the ellipsoids is given the parametric vector equations:
First ellipsoid: $r_1 = C_1 + R_1 {D_1}^{-1/2} u $
Note that:
${D_1}^{-1/2} = \begin{bmatrix} a_1 && 0 && 0 \\ 0 && b_1 && 0 \\ 0 && 0 && c_1 \end{bmatrix}$
The unit vector $u$ is parameterized as follows
$u = \begin{bmatrix} \sin \theta_1 \cos \phi_1 \\ \sin \theta_1 \sin \phi_1 \\ \cos \theta_1 \end{bmatrix} $
Exactly the same structure of the parametric vector equation carries over to the second ellipsoid, thus we have
$r_2 = C_2 + R_2 D_2^{-1/2} v $
where,
${D_2}^{-1/2} = \begin{bmatrix} a_2 && 0 && 0 \\ 0 && b_2 && 0 \\ 0 && 0 && c_2 \end{bmatrix}$
and
$v = \begin{bmatrix} \sin \theta_2 \cos \phi_2 \\ \sin \theta_2 \sin \phi_2 \\ \cos \theta_2 \end{bmatrix} $
Finally, we are ready to carry out the integral that will calculate the volume common to both ellipsoids.
For this, we will scan the surface of both ellipsoids in sequence (the first ellipsoid surface, then the second ellipsoid surface). Let $P_0$ be any point in space, then we will define
$V_1 = \displaystyle \iint_{\theta_1, \phi_1} \dfrac{1}{3} f_1(\theta_1, \phi_1) (r_1 - P_0) \cdot n_1 d \theta_1 d\phi_1 $
where
$f_1( \theta_1, \phi_1)$ is equal to $1$ if $r_1$ (as a function of $\theta_1$ and $\phi_1$) is inside the second ellipsoid, and is zero otherwise. The vector $n_1 $ is defined by
$n_1 = \dfrac{\partial r_1}{\partial \theta_1} \times \dfrac{\partial r_1}{\partial \phi_1} $
And in exactly the same way, we'll define
$V_2 = \displaystyle \iint_{\theta_2, \phi_2} \dfrac{1}{3} f_2(\theta_2, \phi_2) (r_2 - P_0) \cdot n_2 d \theta_2 d\phi_2 $
where
$f_2( \theta_2, \phi_2)$ is equal to $1$ if $r_2$ (as a function of $\theta_2$ and $\phi_2$) is inside the first ellipsoid, and is zero otherwise. The vector $n_2 $ is defined by
$n_2 = \dfrac{\partial r_2}{\partial \theta_2} \times \dfrac{\partial r_2}{\partial \phi_2} $
These integrations are two-dimensional integrations which can be evaluated quite easily, using for example the two dimensional version of Simpson's rule for numerical integration.
The numerical estimate of the volume of the region common to both ellipsoids is then
$V_{\text{Common} } = V_1 + V_2 $
Approximated volume overlap of two ellipsoids
An easy (yet not so ellegant) way to approximate the size of the shared volume of two ellipsoids is sampling their volume.
Suppose we have two Ellipsoids defined as quadrics $Q_{a}, Q_{b}$ along with known volumes $V_a, V_b$. A quadric is a compact representation of an Ellipsoid of arbitrary position and orientation. See this post how to get such a quadric representation from a Centroid and Covariance Matrix. These can be easily obtained from the translation and 3D-rotations your referring to.
Sample ellipsoid volume
Uniformly (random or grid) sample the volume defined by the Quadric $Q_{a}$ to obtain a set of points $X_a$:
$$X_a := \{x_{i}:x^{\top}Q_{a}x < 1\},\; i \in [1, ... , N]$$
Create overlap subset
Now, for all samples within the first ellipsoid, check if their also included within the second ellipsoid $Q_{2}$ to define the subset $X_{a \subset b}$ of $X_a$:
$$X_{a \cap b} := \{x \in X_{a} \wedge x^{\top}Q_{b}x < 1 \}$$
Calculate Overlap Volume from sample ratio
Now the overlap volume can be estimated as the ratio of overlapping points to all points sampled inside the ellipsoid.
$$V_{a \cap b} = \frac{|X_{a \cap b}|}{|X_{a}|} V_{a} = \frac{|X_{a \cap b}|}{N} V_{a}$$
This can be implemented rather efficiently on a GPU. Testing a sample is just two matrix-vector multiplications and one comparison. So depending on your use case, if an estimate is enough you could look at the approximation error with increasing $N$.
Code Example
See this notebook for a detailed code example. Here I also show how to construct a quadric from a rotation and translation.
- 59
-
The sampling method looks very interesting! Can we estimate the volume of just about any sub-region ( of some cube) given by some conditions ( inequalities). Is it easier if the restrictions are quadratic? Thank you for your attention! – orangeskid Jun 28 '22 at 18:58
-
This method requires at least one bounded shape with a finite volume; You cant sample an infinite space. But sure, you can estimate any sub-volume of a cube given some conditions! You only need to be able to check points for this condition for the second step. – rohlemax Jun 29 '22 at 12:25