What is the enclosed volume of an irregular cube given the x,y,z coordinates of the 8 corners?

Question

I have the xyz coordinates of 8 points that forms an irregular-shaped cube. This is an animation, so the cube is undergoing periodic or cyclical shape-change. The co-planarism of any group or set of any combination of 4 of these points can't be assumed (it's a virtual certainty it practically never happens).

This question is very close to my situation: Volume of irregular solid

but (a) has some very specific constraints, and (b) doesn't appear to have been answered despite showing very simple geometry.

Is there a simple formula for computing the internal volume from the xyz coordinates of the corners of this irregular cube?

Answer 1

As pointed out by others in comment, if you don't place any constraint on the vertices/faces, there won't be a formula at all.

When the shape is close to a cube, one reasonable approximation has been suggested by vonbrand in the comment. Namely, linear interpolate between the 4 points delimiting a face.

More precisely, given any quadlateral face $f$, let $\vec{A}_f$, $\vec{B}_f$, $\vec{C}_f$, $\vec{D}_f$ be the $4$ vertices of it, arranged in counterclockwise orientation when viewed from outside of the "cube". We will approximate the face by the one with following parametrization:

$$[0,1]\times[0,1] \ni (s,t) \quad\mapsto\quad \vec{X}(s,t) = (1-t)((1-s)\vec{A}_f + s\vec{B}_f) + t((1-s)\vec{D}_f + s\vec{C}_f)\in \mathbb{R}^3$$

Let $F$ be the collection of the six linear interpolated faces and $\Omega$ be the region bounded by them. As long as these faces doesn't do anything strange like self-intersecting, the volume bounded by them is given by

$$\verb/Volume/ = \int_\Omega dV = \frac13 \int_\Omega \vec{\nabla}\cdot \vec{x} dV = \frac13 \int_{\partial \Omega} \vec{x}\cdot d\vec{S} = \frac13 \sum_{f\in F}\int_{[0,1]^2} \vec{X} \cdot \left(\partial_s \vec{X} \times \partial_t \vec{X}\right) ds dt $$ If one struggle through the algebra, one find

$$\verb/Volume/ = \frac{1}{12}\sum_{f\in F} \bigg( [A_fB_fC_f] + [B_fC_fD_f] + [C_fD_fA_f] + [D_fA_fB_f]\bigg)\tag{*1}$$

where $[uvw]$ is the shorthand for the triple product $\vec{u}\cdot (\vec{v} \times \vec{w})$ of any three vectors $\vec{u},\vec{v},\vec{w}$.

Above formula has an interesting geometric interpretation. Given any face $ABCD$, there are two possible ways to triangular them, either as $\{ ABC, CDA \}$ or $\{ BCD, DAB \}$.

If one construct two tetrahedron with apex at origin and $ABC$, $CDA$ as base, the sum of their volume equals to $\frac16 ( [ABC] + [CDA] )$. Similarly, if one construct two tetrahedron with apex at origin and $BCD$, $DAB$ as base, the sum of their volume equals to $\frac16 ( [BCD] + [DAB] )$.

This means for each face $f$, one can interpret the expression $$\frac{1}{12}\bigg( [A_fB_fC_f] + [B_fC_fD_f] + [C_fD_fA_f] + [D_fA_fB_f]\bigg)$$ as the "average" volume of the two possible ways of forming tetrahedron with apex at origin and base at $A_fB_fC_fD_f$. The expression in $(*1)$ can be interpreted as the average of the volume under all possible triangulation of the faces.

Answer 2

Adding to achille hui's answer, there is a variety of nice formulations.

Let the vertex coordinates be $v_i = (x_i,y_i,z_i)$ for $i\in\{0,\ldots,7\}$ such that $v_i$ and $v_j$ are connected by an edge if and only if the binary representations of $i$ and $j$ differ in exactly one bit.

Furthermore, let us gather all control point coordinates in the matrix $$C = \begin{pmatrix} x_0 & y_0 & z_0 \\ x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_3 & y_3 & z_3 \\ x_4 & y_4 & z_4 \\ x_5 & y_5 & z_5 \\ x_6 & y_6 & z_6 \\ x_7 & y_7 & z_7 \end{pmatrix}$$ We assume trilinear interpolation; this is an example of the general Ansatz $$\begin{align} v: [0,1]^3 &\mapsto \mathbb{R}^3 && \text{(deformation from unit cube)} \\ v(t) &= C^\top\,F(t) && \text{(interpolation between control points)} \\ F: [0,1]^3 &\mapsto \mathbb{R}^8 && \text{(vector of fixed shape functions)} \end{align}$$ I won't go into integration details, but as a consequence of the Cauchy-Binet formula, the oriented volume can be written as a linear combination of all $3\times3$ minors of $C$ with fixed coefficients. The formula can then be expanded to $$\begin{align} V(C) &= \frac{1}{12}\sum_{i,j,k=0,\ldots,7}a_{ijk} x_i y_j z_k \\ a_{ijk} &= \begin{cases} 0 & \text{if there is no unique face with vertices $i,j,k$;} \\ & \text{i.e., if $i,j,k,7-i,7-j,7-k$ are not pairwise distinct} \\ & \text{or none or all of $i,j,k$ are in $\{0,3,5,6\}$ (same bit parity)} \\ +1 & \text{otherwise, if vertices $i,j,k$ of the undeformed cube appear} \\ & \text{in counterclockwise order when their face is viewed} \\ & \text{straight from the outside} \\ -1 & \text{otherwise.} \end{cases} \end{align}$$ Note that $a_{ijk}$ is totally anti-symmetric. This is generally the case with oriented volume formulae.

In our particular case, there are exactly $144$ nonzero $a_{ijk}$, corresponding to $6\times4=24$ size-$3$ determinants that contribute to a surface integral like the one in Achille's answer, whereby each determinant expands to $6$ monomials.

Can we do with less, but perhaps larger, determinants? Sure. Another formulation that still reminds us of the surface-integral approach is $$\begin{align} V(C) &= \frac{Q_{0123}-Q_{4567}-Q_{0145}+Q_{2367}+Q_{0246}-Q_{1357}}{12} \\ \text{where}\quad Q_{hijk} &= \begin{vmatrix} +1 & x_h & y_h & z_h \\ -1 & x_i & y_i & z_i \\ -1 & x_j & y_j & z_j \\ +1 & x_k & y_k & z_k \end{vmatrix} = \begin{vmatrix} x_i + x_h & y_i + y_h & z_i + z_h \\ x_j - x_i & y_j - y_i & z_j - z_i \\ x_k - x_h & y_k - y_h & z_k - z_h \end{vmatrix} \end{align}$$ which uses $6$ determinants. Another intriguing formulation is $$\begin{align} V(C) &= \frac{\Delta_{\text{diag}}+\Delta_{0356}-\Delta_{1247}}{12} \\ \text{where}\quad \Delta_{\text{diag}} &= \begin{vmatrix} 1 & 0 & 0 & 0 & 1 & x_0 & y_0 & z_0 \\ 0 & 1 & 0 & 0 & 0 & x_1 & y_1 & z_1 \\ 0 & 0 & 1 & 0 & 0 & x_2 & y_2 & z_2 \\ 0 & 0 & 0 & 1 & 1 & x_3 & y_3 & z_3 \\ 0 & 0 & 0 & 1 & 0 & x_4 & y_4 & z_4 \\ 0 & 0 & 1 & 0 & 1 & x_5 & y_5 & z_5 \\ 0 & 1 & 0 & 0 & 1 & x_6 & y_6 & z_6 \\ 1 & 0 & 0 & 0 & 0 & x_7 & y_7 & z_7 \end{vmatrix} \\& = \begin{vmatrix} 1 & x_3 - x_4 & y_3 - y_4 & z_3 - z_4 \\ 1 & x_5 - x_2 & y_5 - y_2 & z_5 - z_2 \\ 1 & x_6 - x_1 & y_6 - y_1 & z_6 - z_1 \\ 1 & x_0 - x_7 & y_0 - y_7 & z_0 - z_7 \end{vmatrix} \\ \text{and}\quad \Delta_{hijk} &= \begin{vmatrix} 1 & x_h & y_h & z_h \\ 1 & x_i & y_i & z_i \\ 1 & x_j & y_j & z_j \\ 1 & x_k & y_k & z_k \end{vmatrix} = \begin{vmatrix} x_i - x_h & y_i - y_h & z_i - z_h \\ x_j - x_h & y_j - y_h & z_j - z_h \\ x_k - x_h & y_k - y_h & z_k - z_h \end{vmatrix} \end{align}$$ which uses three determinants. $\Delta_{\text{diag}}$ uses only the spatial diagonal vectors, and $\Delta_{0356}$ and $\Delta_{1247}$ contribute corrections for the influence of the shape of the two interior tetrahedra. Herein, an interior tetrahedron is formed by four vertices that are two edges (one face diagonal) apart from each other.

Closer inspection of the above formula reveals that $V(C)$ is unchanged when the vertices of one (or the other) interior tetrahedron are parallel-shifted, i. e. moved by the same vector.

Can we do with less than three determinants?

For such investigations, the following orthogonal transformation will be useful: $$\begin{align} \underbrace{\begin{pmatrix} x_0 & y_0 & z_0 \\ x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_3 & y_3 & z_3 \\ x_4 & y_4 & z_4 \\ x_5 & y_5 & z_5 \\ x_6 & y_6 & z_6 \\ x_7 & y_7 & z_7 \end{pmatrix}}_{C} &= \underbrace{\frac{1}{2}\begin{pmatrix} -1 & -1 & -1 & -1 & -1 & -1 & -1 & +1 \\ +1 & +1 & -1 & +1 & -1 & -1 & +1 & +1 \\ +1 & -1 & +1 & -1 & +1 & -1 & +1 & +1 \\ -1 & +1 & +1 & +1 & +1 & -1 & -1 & +1 \\ -1 & +1 & +1 & -1 & -1 & +1 & +1 & +1 \\ +1 & -1 & +1 & +1 & -1 & +1 & -1 & +1 \\ +1 & +1 & -1 & -1 & +1 & +1 & -1 & +1 \\ -1 & -1 & -1 & +1 & +1 & +1 & +1 & +1 \end{pmatrix}}_{Q} \underbrace{\begin{pmatrix} u_0 & v_0 & w_0 \\ u_1 & v_1 & w_1 \\ u_2 & v_2 & w_2 \\ u_3 & v_3 & w_3 \\ u_4 & v_4 & w_4 \\ u_5 & v_5 & w_5 \\ u_6 & v_6 & w_6 \\ u_7 & v_7 & w_7 \end{pmatrix}}_{\widetilde{C}} \end{align}$$ Note that $$\begin{align} Q^{-1} &= \frac{1}{2}Q^\top & \det Q &= 16 \end{align}$$ and that the new coordinates indexed $6$ and $7$ measure those parallel shifts that we have identified above as irrelevant.

In the $(u,v,w)$ coordinate system, the volume formula can be written as $$V(C) = \widetilde{V}(\widetilde{C}) = \frac{1}{3}\left( \begin{vmatrix} u_0 & v_0 & w_0 \\ u_1 & v_1 & w_1 \\ u_3 & v_3 & w_3 \end{vmatrix} -\begin{vmatrix} u_0 & v_0 & w_0 \\ u_2 & v_2 & w_2 \\ u_4 & v_4 & w_4 \end{vmatrix} +\begin{vmatrix} u_1 & v_1 & w_1 \\ u_2 & v_2 & w_2 \\ u_5 & v_5 & w_5 \end{vmatrix} \right) +\begin{vmatrix} u_3 & v_3 & w_3 \\ u_4 & v_4 & w_4 \\ u_5 & v_5 & w_5 \end{vmatrix}$$ so the coordinates indexed $6$ and $7$ are unused, which should come as no surprise.

It is now instructive to write $\widetilde{V}(\widetilde{C})$ as a bilinear form acting on two columns of $\widetilde{C}$ and parameterized by the remaining column. This yields $$\begin{align} \widetilde{V}(\widetilde{C}) &= \frac{1}{3}\begin{pmatrix} v_0 \\ v_1 \\ v_2 \\ v_3 \\ v_4 \\ v_5 \end{pmatrix}^\top \underbrace{\begin{pmatrix} 0 & u_3 & -u_4 & -u_1 & u_2 & 0 \\ -u_3 & 0 & u_5 & u_0 & 0 & -u_2 \\ u_4 & -u_5 & 0 & 0 & -u_0 & u_1 \\ u_1 & -u_0 & 0 & 0 & 3u_5 & -3u_4 \\ -u_2 & 0 & u_0 & -3u_5 & 0 & 3u_3 \\ 0 & u_2 & -u_1 & 3u_4 & -3u_3 & 0 \end{pmatrix}}_{U} \begin{pmatrix} w_0 \\ w_1 \\ w_2 \\ w_3 \\ w_4 \\ w_5 \end{pmatrix} \end{align}$$ (If we had tried that with the original vertex coordinates, the expressions in the associated matrix would have been more complex.)

There is no single-det formula for $V(C)$ in the strict (real) sense. If there were, we could transform it to our $(u,v,w)$ system and write it as $$\begin{align} \widetilde{V}(\widetilde{C}) &= \begin{vmatrix} * & * & * & u_0 & v_0 & w_0 \\ * & * & * & u_1 & v_1 & w_1 \\ * & * & * & u_2 & v_2 & w_2 \\ * & * & * & u_3 & v_3 & w_3 \\ * & * & * & u_4 & v_4 & w_4 \\ * & * & * & u_5 & v_5 & w_5 \end{vmatrix} \end{align}$$ with the $*$ being unspecified constants. Then there would be a $3$-dimensional space for $(u_0,\ldots,u_5)$, namely the span of the unspecified columns of the above matrix, which would force volume zero regardless of $v$ and $w$. That is, the matrix $U$ corresponding to the bilinear form for $\widetilde{V}(\widetilde{C})$ would have to be all zeros in such cases. Formally, the linear map $\Phi$ from $(u_0,\ldots,u_5)$ to $U$ would need to have a $3$-dimensional kernel. However, we see that $u_0,\ldots,u_5$ occur as individual entries in $U$, so the kernel of $\Phi$ has dimension zero.

Does a two-det formula for $\widetilde{V}(\widetilde{C})$ exist? Or, roughly speaking, can we write $$\begin{align} \widetilde{V}(\widetilde{C}) &= \begin{vmatrix} * & * & * & u_0 & v_0 & w_0 \\ * & * & * & u_1 & v_1 & w_1 \\ * & * & * & u_2 & v_2 & w_2 \\ * & * & * & u_3 & v_3 & w_3 \\ * & * & * & u_4 & v_4 & w_4 \\ * & * & * & u_5 & v_5 & w_5 \end{vmatrix} + \begin{vmatrix} * & * & * & u_0 & v_0 & w_0 \\ * & * & * & u_1 & v_1 & w_1 \\ * & * & * & u_2 & v_2 & w_2 \\ * & * & * & u_3 & v_3 & w_3 \\ * & * & * & u_4 & v_4 & w_4 \\ * & * & * & u_5 & v_5 & w_5 \end{vmatrix} \end{align}$$ To answer this question, consider the case that $(u_0,\ldots,u_5)$ is in the span of the unspecified columns of the left matrix. Then the left determinant becomes zero. If furthermore $(w_0,\ldots,w_5)$ is in the span of the leftmost four columns of the right matrix, we get $V=0$ regardless of $(v_0,\ldots,v_5)$. This implies that the associated matrix $U$ must have rank less than four. Then all its $4\times4$ minors must be zero. It turns out that this happens if and only if $$\begin{pmatrix}u_0\\u_1\\u_2\end{pmatrix} = s\begin{pmatrix}u_5\\u_4\\u_3\end{pmatrix} \quad\text{where}\quad s^2 = -3$$ which is impossible for real nonzero $(u_0,\ldots,u_5)$, but possible if we allow complex entries. In the latter case, the two determinants summed in the desired volume formula are complex conjugates, so we can write the volume as real part of a single complex determinant: $$\begin{align} \widetilde{V}(\widetilde{C}) &= \operatorname{Re}\begin{vmatrix} 1 & 0 & 0 & u_0 & v_0 & w_0 \\ 0 & 1 & 0 & u_1 & v_1 & w_1 \\ 0 & 0 & 1 & u_2 & v_2 & w_2 \\ 0 & 0 & r & u_3 & v_3 & w_3 \\ 0 & r & 0 & u_4 & v_4 & w_4 \\ r & 0 & 0 & u_5 & v_5 & w_5 \end{vmatrix} = \operatorname{Re}\begin{vmatrix} u_3 - r u_2 & v_3 - r v_2 & w_3 - r w_2 \\ u_4 - r u_1 & v_4 - r v_1 & w_4 - r w_1 \\ u_5 - r u_0 & v_5 - r v_0 & w_5 - r w_0 \end{vmatrix} \\ r^2 &= \frac{1}{s^2} = -\frac{1}{3} \end{align}$$ Note the similarity to the octahedron's volume formula.

In the original $(x,y,z)$ coordinate system, the volume formula can be written as $$\begin{align} \text{Let}\quad t &= \frac{\zeta_{12}}{\sqrt{3}} = \frac{1}{2} + \frac{\mathrm{i}}{2\sqrt{3}} = \frac{1 + \zeta_6}{3},\quad \zeta_n = \exp\frac{2\pi\mathrm{i}}{n} \\ \text{Then}\quad V(C) &= \frac{1}{2}\operatorname{Re}\begin{vmatrix} 1 & 0 & 0 & 0 & 1 & t x_0 & t y_0 & t z_0 \\ 0 & 1 & 0 & 0 & 0 & \bar{t}x_1 & \bar{t}y_1 & \bar{t}z_1 \\ 0 & 0 & 1 & 0 & 0 & \bar{t}x_2 & \bar{t}y_2 & \bar{t}z_2 \\ 0 & 0 & 0 & 1 & 1 & t x_3 & t y_3 & t z_3 \\ 0 & 0 & 0 & 1 & 0 & \bar{t}x_4 & \bar{t}y_4 & \bar{t}z_4 \\ 0 & 0 & 1 & 0 & 1 & t x_5 & t y_5 & t z_5 \\ 0 & 1 & 0 & 0 & 1 & t x_6 & t y_6 & t z_6 \\ 1 & 0 & 0 & 0 & 0 & \bar{t}x_7 & \bar{t}y_7 & \bar{t}z_7 \end{vmatrix} \\& = \frac{1}{2}\operatorname{Re}\begin{vmatrix} 1 & t x_3 - \bar{t}x_4 & t y_3 - \bar{t}y_4 & t z_3 - \bar{t}z_4 \\ 1 & t x_5 - \bar{t}x_2 & t y_5 - \bar{t}y_2 & t z_5 - \bar{t}z_2 \\ 1 & t x_6 - \bar{t}x_1 & t y_6 - \bar{t}y_1 & t z_6 - \bar{t}z_1 \\ 1 & t x_0 - \bar{t}x_7 & t y_0 - \bar{t}y_7 & t z_0 - \bar{t}z_7 \end{vmatrix} \end{align}$$ which reminds us of $\Delta_{\text{diag}}$. In terms of $3\times3$ determinants, the volume formula is $$\begin{align} V(C) &= \frac{1}{2}\operatorname{Re}\begin{vmatrix} t (x_6 - x_0) + \bar{t}(x_7 - x_1) & t (y_6 - y_0) + \bar{t}(y_7 - y_1) & t (z_6 - z_0) + \bar{t}(z_7 - z_1) \\ t (x_5 - x_0) + \bar{t}(x_7 - x_2) & t (y_5 - y_0) + \bar{t}(y_7 - y_2) & t (z_5 - z_0) + \bar{t}(z_7 - z_2) \\ t (x_3 - x_0) + \bar{t}(x_7 - x_4) & t (y_3 - y_0) + \bar{t}(y_7 - y_4) & t (z_3 - z_0) + \bar{t}(z_7 - z_4) \end{vmatrix} \end{align}$$ You might regard that as a single-det formula, at the cost of complex-number trickery.