Help with least squares ellipsoid specific fitting

Question

In following the method of Li and Griffiths paper which can be freely accessed here, I created the following toy problem to see if I am performing the calculations correctly. Below I describe the problem and I give my intermediate and final calculations. I'm hoping to understand where and what I am doing wrong.

The problem: Using an origin-centered ellipsoid defined by the equation $$\tag{1} 0.25x^2+1y^2+4z^2=1,$$ the following seven points $p_i(x_i,y_i,z_i)$ which satisfy Eqn(1) were selected:

A plot of these points along with the ellipsoid are as shown here:

My attempt to fit an ellipsoid to this set of points following the method of Li and Griffiths: With $k$ set equal to 4, the constraint matrix $C_1$ is:

The matrix $DD^T$ is (see Section 3 of paper, Eqn(11)):

The matrix $A=S_{11}-S_{12}S_{22}^{-1}S_{12}^T$ found in Eqn(15) of the paper, is:

The matrix $C_1^{-1}A$ is:

The largest eigenvalue $\lambda$ (the only positive eigenvalue) of matrix $C_1^{-1}A$ and associated eigenvector $u_1$ is:

Vector $u_2=-S_{22}^{-1}S_{12}^Tu_1$ is equal to 0.62. Thus, with the ellipsoid's parameter vector $u=(u_1,u_2)$, the best fit ellipsoid is found as $$\tag{2} ax^2+by^2+cz^2+2fyz+2gxz+2hxy+d=0$$ with the coefficients:

The value of $kJ-I^2$, where $I$ and $J$ are defined in Eqns(2) & (3) of the paper, is 0.3 while $k=4$. The value of $k$ needed to make $kJ-I^2=1$ is found to be 5.298, with the new associated parameter vector $u$ is

This is clearly wrong. Where have I gone wrong? What have I done wrong?

Answer 1

In the following, I will use a different method than the one suggested in your referenced paper.

Since it is given by your example problem that the ellipsoid to be identified is origin centered, then the algebraic equation of the ellipsoid is

$ \mathbf{r}^T Q \mathbf{r} = 1 $

where $\mathbf{r} = [x, y, z]^T $ and $Q$ is a symmetric positive definite $3 \times 3$ matrix.

Given a set of data points $\{ P_i = (x_i, y_i, z_i) , i = 1, N\}, N \ge 6 $, you want to fit these points to the model

$ A x^2 + B y^2 + C z^2 + D xy + E xz + F yz + G = 0 $

For that, you construct an error function $\mathbf{E}$ defined as follows

$ \mathbf{E} = \displaystyle \sum_{i=1}^N (A x^2 + B y^2 + C z^2 + D xy + E xz + F yz + G)^2 $

And you minimize this function. To ensure that you get a valid ellipsoid with the parameters $A$ through $G$ not all zero, you make the minimization conditioned on the constraint $A + B + C = 1 $

Substitute for $A$ in the above expression you get

$ \mathbf{E} = \displaystyle \sum_{i=1}^N ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G )^2 $

Differentiating $\mathbf{E}$ with respect to $B, C, D, E, F, G$ results in the following equations for the minimum

$ \displaystyle \sum_{i=1}^N ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G ) (y_i^2 - x_i^2) = 0 $

$ \displaystyle \sum_{i=1}^N ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G ) (z_i^2 - x_i^2) = 0 $

$ \displaystyle \sum_{i=1}^N ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G ) (x_i y_i) = 0 $

$ \displaystyle \sum_{i=1}^N ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G ) (x_i z_i) = 0 $

$ \displaystyle \sum_{i=1}^N ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G ) (y_i z_i) = 0 $

$ \displaystyle \sum_{i=1}^N ( ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G ) = 0 $

The above $6$ equations lead to the following linear system for finding $B, C, D, E, F, G$

Define the vector $X = [B, C, D, E, F, G]^T $, then the normal equations are

$ \mathbf{A} X = Y $

where

$ \mathbf{A} = \displaystyle \sum_{i=1}^N A_i $

with

$A_i = \begin{bmatrix} (y_i^2 - x_i^2)^2 && (z_i^2 - x_i^2)(y_i^2 - x_i^2) && x_i y_i (y_i^2 - x_i^2) && x_i z_i (y_i^2 - x_i^2) && y_i z_i (y_i^2 - x_i^2) && (y_i^2 - x_i^2 ) \\ (y_i^2 - x_i^2) (z_i^2 - x_i^2) && (z_i^2 - x_i^2)^2 && x_i y_i (z_i^2 - x_i^2) && x_i z_i (z_i^2 - x_i^2) && y_i z_i (z_i^2 - x_i^2) && (z_i^2 - x_i^2 )\\ (y_i^2 - x_i^2) (x_i y_i) && (z_i^2 - x_i^2) (x_i y_i) && (x_i y_i)^2 && x_i^2 y_i z_i && x_i y_i^2 z_i && x_i y_i \\ (y_i^2 - x_i^2) (x_i z_i) && (z_i^2 - x_i^2) (x_i z_i) && x_i^2 y_i z_i && x_i^2 z_i^2 && x_i y_i z_i^2 && x_i z_i \\ (y_i^2 - x_i^2) (y_i z_i) && (z_i^2 - x_i^2) (y_i z_i) && (x_i y_i)(y_i z_i) && (x_i z_i) (y_i z_i) && (y_i z_i)^2 && y_i z_i \\ (y_i^2 - x_i^2) && (z_i^2 - x_i^2) && (x_i y_i) && (x_i z_i) && (y_i z_i) && 1 \end{bmatrix}$

And

$ Y = \begin{bmatrix} - \sum x_i^2 (y_i^2 - x_i^2) \\ - \sum x_i^2 (z_i^2 - x_i^2) \\ - \sum x_i^2 (x_i y_i) \\ - \sum x_i^2 (x_i z_i) \\ - \sum x_i^2 (y_i z_i) \\ - \sum x_i^2 \end{bmatrix}$

Now we have our linear system which can be solve for the parameter vector $X$ using Gauss-Jordan elimination, which is a standard routine. Having found $X$ all the parameters of the ellipsoid are known, and it can be written in quadratic form as

$ r^T Q r = 0 $

where $r = [x, y, z, 1]^T$ , and

$ Q = \begin{bmatrix} A && D/2 && E/2 && 0 \\ D/2 && B && F/2 && 0 \\ E/2 && F/2 && C && 0 \\ 0&&0&&0 && G \end{bmatrix} $

Now define the $3 \times 3$ matrix

$ Q_0 = \begin{bmatrix} A && D / 2 && E / 2 \\ D / 2 && B && F/2 \\ E/2 && F/2 && C \end{bmatrix}$

Which enables us to write the equation of the fitted ellipsoid in the form

$ \mathbf{r}^T Q_0 \mathbf{r} = - G $

where $p = [x, y, z]^T $

Dividing both sides by $(-G ) $ gives us

$ \mathbf{r}^T Q_1 \mathbf{r} = 1$

where $Q_1 = Q_0 / ( - G) $

Using the data in the first table that you provided, I got

$Q_1 = \begin{bmatrix} 0.25 && -9 \times 10^{-7} && -2.7 \times 10^{-6} \\ -9\times 10^{-7} && 0.999999 && 1.56 \times 10^{-5} \\ -2.7 \times 10^{-6} && 1.56 \times 10^{-5} && 4 \end{bmatrix} $