Best-fit cuboid to 8 points

Question

Given 8 target points $p_0,...,p_7$ in $\mathbb{R}^3$, what is the cuboid (as defined on that wikipedia page-- not necessarily rectangular) whose ordered vertices are as close as possible to the respective target points, in the sense of minimizing the sums-of-squares-of-distances?

I know this can be solved straightforwardly using numerical optimization, by parametrizing the cuboid by 18 scalar parameters in either of the following ways:

• the plane equations (3 scalar params each) of the 6 cuboid face planes
• any 6 independent vertices (3 scalar params each) out of the 8 cuboid vertices (there are 4 choices for the 6 independent vertices)

But I'd like a direct closed-form formula for the solution points (or solution face planes), in terms of the target points, instead.

(More generally, I'd like to be able to solve this for the vertices of any polyhedral graph, such as that of a square pyramid or triangular prism. But this question focuses on a cuboid for concreteness; I suspect when the cuboid solution is understood, the generalization to other polyhedra may become obvious.)

Some observations

Similar problems:

• Best-fit rigid transformation (i.e. translation and rotation) of unit cube (or any given rigid cuboid) to the target points. This is a case of Wahba's problem or the orthogonal Procrustes problem, solvable by aligning the centroids and then using SVD to find the best-fit rotation.
• Best-fit parallelopiped. This is linear least squares (aka non-orthogonal Procrustes), solvable by aligning the centroid of a unit cube to that of the target points and then finding the best-fit linear transformation, by solving the normal equations (or any of many other methods).

It's not too hard to see that the solution will have the following properties (and solutions to the two similar problems mentioned above also have these properties):

• The centroid (average) of the solution vertices will be the same that of as the target vertices (why? because, if not, rigidly translating the solution so that its centroid is the same as that of the target would yield an even better solution)
• If the spatial differences between the solution points and the target points are interpreted as forces, then:
  • there is zero net force (this is the same as saying the solution verts' centroid is the target verts' centroid, i.e. the property already mentioned above)
  • there is zero net torque, about the centroid or anywhere else (why? because, if not, rigidly rotating the solution a bit in the direction of the torque would yield an even better solution).

One obvious way of approximating a solution would be to simply set the 6 cuboid face planes to be the 6 best-fit planes through each respective set of 4 of the 8 target points. But the cuboid formed in this way will not be optimal, in general, in the required sense of minimizing the sums-of-squares-of-errors of the vertices. In particular, it's easy to come up with well-behaved solvable examples for which this method produces vertices arbitrarily far away, or fails to produce vertices at all (i.e. a triple of best-fit face planes whose common intersection "should" be a given cuboid vertex may not have a unique intersection point).

Answer 1

PARTIAL SOLUTION

First consider the case where the directions of the edges of the cuboid are fixed. Without loss of generality assume that the edges of the cuboid are parallel to the $x, y, z$ axes. Then the sum of squares of distance can be minimized indenpendently for each face. For example. if one of the faces that are perpendicular to the $x$ axis corresponds to points $p_0, p_1, p_2, p_3$, and we want to find the optimal $x$ coordinate $a$ to put the face. This reduces to minimizing $\sum_{i=0}^3 (x_i-a)^2$ where $x_i$ is the $x$ coordinate of $p_i$, so the optimal solution is $a = \frac{1}{4}(x_0+x_1+x_2+x_3)$. This means that if $q_0 := \frac{1}{4}(p_0+p_1+p_2+p_3) \in \mathbb{R}^3$, the plane that corresponds to points $p_0, p_1, p_2, p_3$ will always contain the point $q_0$. If $p'_i := p_i-q_0$ for $i \in \{0,1,2,3\}$, the objective function for this face is $\sum_{i=0}^3 (e_x \cdot p'_i)^2$ where $e_x = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}$.

Now we want to find the optimal rotation for the cuboid. If the cuboid needs to be rotated with an orthogonal matrix $R$, the objective function for the face $p_0, p_1, p_2, p_3$ is $\sum_{i=0}^3 (e_x \cdot Rp'_i)^2$. Letting $e_i := e_x$ this can be rewritten as $\sum_{i=0}^3 (e_i \cdot Rp'_i)^2$. By similar arguments for the other faces the problem reduces to minimizing $\sum_{i=0}^{23} (e_i \cdot Rp'_i)^2$ with the constraint that $R$ is an orthogonal matrix, where $e_0, e_1, \dots, e_{23}, p'_0, p'_1, \dots, p'_{23}$ are known vectors.

Next we express rotations in quaternions. Let $q$ be the unit quaternion corresponding to $R$, and $a_i, b_i$ be the quaternion for the vectors $e_i, p'_i$. The objective function becomes $\sum_{i=0}^{23} -\operatorname{Re}(a_iqb_i\overline{q})^2$. I could not obtain a closed expression for the optimal solution, but this form should be easier to optimize numerically than the original problem.

Answer 2

PARTIAL SOLUTION

It seems useful to organize the problem using the method of Lagrange multipliers with multiple constraints, which transforms it into a system of 30 degree-3 polynomial equations in 30 scalar unknowns, as follows.

Call the target points $\mathbf{\tilde{p_0}},...,\mathbf{\tilde{p_7}}$.

We want to find the values of the 8 cuboid vertices (24 scalar variables): $$\begin{aligned} \mathbf{p_0} &= (x_0,y_0,z_0)^T \\ &... \\ \mathbf{p_7} &= (x_7,y_7,z_7)^T \\ \end{aligned}$$ that minimizes the scalar function: $$\begin{aligned} f(\mathbf{p_0},...\mathbf{p_7}) = \sum{{\left\lVert\mathbf{p_i}-\mathbf{\tilde{p_i}}\right\rVert}^2} \end{aligned}$$ subject to the following 6 scalar constraints (which say that the 4 vertices of each of the 6 faces of the cuboid must be coplanar): $$\begin{aligned} g_0(\mathbf{p_0},...,\mathbf{p_7}) &= ((\mathbf{p_2}-\mathbf{p_0})\times(\mathbf{p_3}-\mathbf{p_0}))\cdot(\mathbf{p_1}-\mathbf{p_0}) = 0 \\ g_1(\mathbf{p_0},...,\mathbf{p_7}) &= ((\mathbf{p_1}-\mathbf{p_0})\times(\mathbf{p_5}-\mathbf{p_0}))\cdot(\mathbf{p_4}-\mathbf{p_0}) = 0 \\ g_2(\mathbf{p_0},...,\mathbf{p_7}) &= ((\mathbf{p_4}-\mathbf{p_0})\times(\mathbf{p_6}-\mathbf{p_0}))\cdot(\mathbf{p_2}-\mathbf{p_0}) = 0 \\ g_3(\mathbf{p_0},...,\mathbf{p_7}) &= ((\mathbf{p_5}-\mathbf{p_7})\times(\mathbf{p_1}-\mathbf{p_7}))\cdot(\mathbf{p_3}-\mathbf{p_7}) = 0 \\ g_4(\mathbf{p_0},...,\mathbf{p_7}) &= ((\mathbf{p_3}-\mathbf{p_7})\times(\mathbf{p_2}-\mathbf{p_7}))\cdot(\mathbf{p_6}-\mathbf{p_7}) = 0 \\ g_5(\mathbf{p_0},...,\mathbf{p_7}) &= ((\mathbf{p_6}-\mathbf{p_7})\times(\mathbf{p_4}-\mathbf{p_7}))\cdot(\mathbf{p_5}-\mathbf{p_7}) = 0 \end{aligned}$$

Notice that the expressions for the $g_i$'s are $3{\times}3$ determinants which are essentially tetrahedron volumes (times 6), so each of them can be rearranged in various ways by (evenly) permuting its 4 points, when convenient. E.g. when we want to differentiate $g_0$ with respect to $\mathbf{p_0}$, it helps to think of it as $g_0 = ((\mathbf{p_3}-\mathbf{p_1})\times(\mathbf{p_2}-\mathbf{p_1}))\cdot(\mathbf{p_0}-\mathbf{p_1})$, whose partial derivative with respect to $\mathbf{p_0}$ can then be clearly seen to be $(\mathbf{p_3}-\mathbf{p_1})\times(\mathbf{p_2}-\mathbf{p_1})$ (a row vector).

Following the method as described on the wikipedia page, the Lagrangian function is: $$ \mathcal{L}(\mathbf{p_0},...,\mathbf{p_7},\lambda_0,...,\lambda_5) = f(\mathbf{p_0},...\mathbf{p_7}) - \sum{\lambda_j g_j(\mathbf{p_0},...,\mathbf{p_7})} $$ The solution will be a stationary point of $\mathcal{L}$, i.e. a point at which all 30 scalar partial derivatives in the gradient of $\mathcal{L}$ are zero. This is a system of 30 equations in 30 scalar unknowns.

Expanding the first three of these partial derivatives which we want to be zero, $\frac{\partial \mathcal{L}}{\partial \mathbf{p_0}} = \left(\frac{\partial \mathcal{L}}{\partial x_0}, \frac{\partial \mathcal{L}}{\partial y_0}, \frac{\partial \mathcal{L}}{\partial z_0}\right)$, as an example of how it looks: $$\begin{aligned} \frac{\partial \mathcal{L}}{\partial \mathbf{p_0}} &= \frac{\partial f}{\partial \mathbf{p_0}} - \sum{\lambda_j \frac{\partial g_j}{\partial \mathbf{p_0}}} \\ &= \frac{\partial f}{\partial \mathbf{p_0}} - \lambda_0 \frac{\partial g_0}{\partial \mathbf{p_0}} - \lambda_1 \frac{\partial g_1}{\partial \mathbf{p_0}} - \lambda_2 \frac{\partial g_2}{\partial \mathbf{p_0}} \text{ (since }g_3,g_4,g_5\text{ do not depend on }\mathbf{p_0}\text{)} \\ &= 2\left(\mathbf{p_0} {-} \mathbf{\tilde{p_0}}\right) - \lambda_0 (\mathbf{p_3}{-}\mathbf{p_1}){\times}(\mathbf{p_2}{-}\mathbf{p_1}) - \lambda_1 (\mathbf{p_5}{-}\mathbf{p_4}){\times}(\mathbf{p_1}{-}\mathbf{p_4}) - \lambda_2 (\mathbf{p_6}{-}\mathbf{p_2}){\times}(\mathbf{p_4}{-}\mathbf{p_2}) \end{aligned}$$

And expanding $\frac{\partial \mathcal{L}}{\partial \lambda_5}$ as another representative example: $$\begin{aligned} \frac{\partial \mathcal{L}}{\partial \lambda_0} &= -g_0(\mathbf{p_0},...,\mathbf{p_7}) \\ &= - ((\mathbf{p_2}-\mathbf{p_0})\times(\mathbf{p_3}-\mathbf{p_0}))\cdot(\mathbf{p_1}-\mathbf{p_0}) \end{aligned}$$

Here is the full system of $8{\times}3{+}6{=}30$ polynomial equations in the 30 unknowns, expanded(ish). (There are probably typos.) $$\begin{aligned} \frac{\partial \mathcal{L}}{\partial \mathbf{p_0}} &{=} 2\left(\mathbf{p_0} {-} \mathbf{\tilde{p_0}}\right) {-} \lambda_0 (\mathbf{p_3}{-}\mathbf{p_1}){\times}(\mathbf{p_2}{-}\mathbf{p_1}) {-} \lambda_1 (\mathbf{p_5}{-}\mathbf{p_4}){\times}(\mathbf{p_1}{-}\mathbf{p_4}) {-} \lambda_2 (\mathbf{p_6}{-}\mathbf{p_2}){\times}(\mathbf{p_4}{-}\mathbf{p_2}) {=} (0,0,0) \\ \frac{\partial \mathcal{L}}{\partial \mathbf{p_1}} &{=} 2\left(\mathbf{p_1} {-} \mathbf{\tilde{p_1}}\right) {-} \lambda_0 (\mathbf{p_2}{-}\mathbf{p_0}){\times}(\mathbf{p_3}{-}\mathbf{p_0}) {-} \lambda_1 (\mathbf{p_4}{-}\mathbf{p_5}){\times}(\mathbf{p_0}{-}\mathbf{p_5}) {-} \lambda_3 (\mathbf{p_7}{-}\mathbf{p_5}){\times}(\mathbf{p_3}{-}\mathbf{p_5}) {=} (0,0,0) \\ \frac{\partial \mathcal{L}}{\partial \mathbf{p_2}} &{=} 2\left(\mathbf{p_2} {-} \mathbf{\tilde{p_2}}\right) {-} \lambda_0 (\mathbf{p_1}{-}\mathbf{p_3}){\times}(\mathbf{p_0}{-}\mathbf{p_3}) {-} \lambda_2 (\mathbf{p_4}{-}\mathbf{p_0}){\times}(\mathbf{p_6}{-}\mathbf{p_0}) {-} \lambda_4 (\mathbf{p_7}{-}\mathbf{p_3}){\times}(\mathbf{p_6}{-}\mathbf{p_3}) {=} (0,0,0) \\ \frac{\partial \mathcal{L}}{\partial \mathbf{p_3}} &{=} 2\left(\mathbf{p_3} {-} \mathbf{\tilde{p_3}}\right) {-} \lambda_0 (\mathbf{p_0}{-}\mathbf{p_2}){\times}(\mathbf{p_1}{-}\mathbf{p_2}) {-} \lambda_3 (\mathbf{p_5}{-}\mathbf{p_7}){\times}(\mathbf{p_1}{-}\mathbf{p_7}) {-} \lambda_4 (\mathbf{p_6}{-}\mathbf{p_2}){\times}(\mathbf{p_7}{-}\mathbf{p_2}) {=} (0,0,0) \\ \frac{\partial \mathcal{L}}{\partial \mathbf{p_4}} &{=} 2\left(\mathbf{p_4} {-} \mathbf{\tilde{p_4}}\right) {-} \lambda_1 (\mathbf{p_1}{-}\mathbf{p_0}){\times}(\mathbf{p_5}{-}\mathbf{p_0}) {-} \lambda_2 (\mathbf{p_2}{-}\mathbf{p_6}){\times}(\mathbf{p_0}{-}\mathbf{p_6}) {-} \lambda_5 (\mathbf{p_7}{-}\mathbf{p_6}){\times}(\mathbf{p_5}{-}\mathbf{p_6}) {=} (0,0,0) \\ \frac{\partial \mathcal{L}}{\partial \mathbf{p_5}} &{=} 2\left(\mathbf{p_5} {-} \mathbf{\tilde{p_5}}\right) {-} \lambda_1 (\mathbf{p_0}{-}\mathbf{p_1}){\times}(\mathbf{p_4}{-}\mathbf{p_1}) {-} \lambda_3 (\mathbf{p_3}{-}\mathbf{p_1}){\times}(\mathbf{p_7}{-}\mathbf{p_1}) {-} \lambda_5 (\mathbf{p_6}{-}\mathbf{p_7}){\times}(\mathbf{p_4}{-}\mathbf{p_7}) {=} (0,0,0) \\ \frac{\partial \mathcal{L}}{\partial \mathbf{p_6}} &{=} 2\left(\mathbf{p_6} {-} \mathbf{\tilde{p_6}}\right) {-} \lambda_2 (\mathbf{p_0}{-}\mathbf{p_4}){\times}(\mathbf{p_2}{-}\mathbf{p_4}) {-} \lambda_4 (\mathbf{p_3}{-}\mathbf{p_7}){\times}(\mathbf{p_2}{-}\mathbf{p_7}) {-} \lambda_5 (\mathbf{p_5}{-}\mathbf{p_4}){\times}(\mathbf{p_7}{-}\mathbf{p_4}) {=} (0,0,0) \\ \frac{\partial \mathcal{L}}{\partial \mathbf{p_7}} &{=} 2\left(\mathbf{p_7} {-} \mathbf{\tilde{p_7}}\right) {-} \lambda_3 (\mathbf{p_1}{-}\mathbf{p_3}){\times}(\mathbf{p_5}{-}\mathbf{p_3}) {-} \lambda_4 (\mathbf{p_2}{-}\mathbf{p_6}){\times}(\mathbf{p_3}{-}\mathbf{p_6}) {-} \lambda_5 (\mathbf{p_4}{-}\mathbf{p_5}){\times}(\mathbf{p_6}{-}\mathbf{p_5}) {=} (0,0,0) \\ \frac{\partial \mathcal{L}}{\partial \lambda_0} &= - ((\mathbf{p_2}-\mathbf{p_0})\times(\mathbf{p_3}-\mathbf{p_0}))\cdot(\mathbf{p_1}-\mathbf{p_0}) = 0 \\ \frac{\partial \mathcal{L}}{\partial \lambda_1} &= - ((\mathbf{p_1}-\mathbf{p_0})\times(\mathbf{p_5}-\mathbf{p_0}))\cdot(\mathbf{p_4}-\mathbf{p_0}) = 0 \\ \frac{\partial \mathcal{L}}{\partial \lambda_2} &= - ((\mathbf{p_4}-\mathbf{p_0})\times(\mathbf{p_6}-\mathbf{p_0}))\cdot(\mathbf{p_2}-\mathbf{p_0}) = 0 \\ \frac{\partial \mathcal{L}}{\partial \lambda_3} &= - ((\mathbf{p_5}-\mathbf{p_7})\times(\mathbf{p_1}-\mathbf{p_7}))\cdot(\mathbf{p_3}-\mathbf{p_7}) = 0 \\ \frac{\partial \mathcal{L}}{\partial \lambda_4} &= - ((\mathbf{p_3}-\mathbf{p_7})\times(\mathbf{p_2}-\mathbf{p_7}))\cdot(\mathbf{p_6}-\mathbf{p_7}) = 0 \\ \frac{\partial \mathcal{L}}{\partial \lambda_5} &= - ((\mathbf{p_6}-\mathbf{p_7})\times(\mathbf{p_4}-\mathbf{p_7}))\cdot(\mathbf{p_5}-\mathbf{p_7}) = 0 \\ \end{aligned}$$ (Notice that the last six of them just re-state the six coplanarity constraints.)

At this point, I'm not sure how to proceed. It seems like we've probably made some progress, compared to the original minimization problem, but I don't see how to solve this system of 30 degree-3 polynomials in 30 unknowns (other than by using a numeric solver; Newton's method would probably work well for this).