I have a point cloud sampled from the surface of an elliptic paraboloid, whose axis is not the coordinate axes x, y, and z. I want to determine the equation of the axis.
Method1: I used the equation: $a x^{2}+b x y+c y^{2}+d x+e y+f=z$ and determined the coefficients using scipy. However, I do not know how to determine the equation of axis from it.
Method2: I wanted to find the general equation for an elliptic paraboloid with random axis. After a lot of searching, I failed to find such equation.
How can I determine the axis of an elliptic paraboloid with random axis.
The equation of the basic elliptic paraboloid with vertex at the origin and opening in the direction of the positive $z$ axis, while its $z$-cross sections are ellipses with major/minor axes along the $x, y$ coordinate axes is given by
$ z = a x^2 + b y^2 $
This paraboloid equation can be written in compact form as
$ r^T Q_0 r + b_0^T r = 0 $
where $r = [x, y, z]^T $ is the position coordinate vector of a point on the elliptic paraboloid. And
$ Q_0 = \begin{bmatrix} a && 0 && 0 \\ 0 && b && 0 \\ 0 && 0 && 0 \end{bmatrix} $
and
$ b_0 = \begin{bmatrix} 0 \\ 0 \\ -1 \end{bmatrix} $
Now if you shift the vertex of the paraboloid to a new location $r_0$, then image of $r$ is $r'$ with
$r' = r_0 + r $
and then you apply a rotation about the point $r_0$ whose rotation matrix is $R$, then the image of a point $r'$ due to the rotation, is
$ r'' = r_0 + R (r' - r_0 ) = r_0 + R r $
Therefore, to get the equation of the shifted/rotated paraboloid, we solve for $r$ from the last equation:
$ r = R^T (r'' - r_0) $
And we plug this into the equation of the original paraboloid, to get
$ (r'' - r_0)^T R Q_0 R^T (r'' - r_0) + b_0^T R^T (r'' - r_0) = 0 $
Now, we can drop the primes, and replace the variable name $r''$ with simple $r$:
$ (r - r_0)^T R Q_0 R^T (r - r_0) + b_0^T R^T (r - r_0) = 0 $
Written compactly, the last equation becomes
$ (r - r_0)^T Q (r - r_0) + b^T (r - r_0) = 0 \hspace{20pt} (*) $
where $ Q = R Q_0 R^T $ and $ b = R b_0 $
Expanding equation $(*)$ in terms of $x, y, z$ yields
$ A x^2 + B y^2 + C z^2 + D xy + E xz + F yz \\ + G x + H y + I z + J = 0 \hspace{20pt} (**)$
Now, given the point cloud, as suggested by @Narasimham, you can use the method of least squares to identify the quadratic equation $(**)$. Suppose you have $N$ points where $N \ge 9$, we will define
$ \theta = [ A, B, C, D, E, F, G, H, I, J]^T $
as the parameter vector, and the matrix $\Phi$ whose $i$-th row is
$[ x_i^2 , y_i^2 , z_i^2, x_i y_i, x_i z_i, y_i, z_i , x_i, y_i, z_i, 1 ] $
The equation for all the data is
$ \Phi \theta = \epsilon $
Since we want the minimum of $\epsilon^T \epsilon$, then $\theta$ will be the eigenvector of $ \Phi^T \Phi$ corresponding to the minimum eigenvalue.
That is one way to do it, the other way that I can think of is to pre-multiply the above equation through by $\Phi^T$ to get
$ \Phi^T \Phi \theta = \Phi^T \epsilon $
Assuming that $\epsilon = 0 $, and truncating the matrix $\Phi^T \Phi $ (which is $10 \times 10 $ to the first $9$ rows, a direction for $\theta$ can be found using the standard Gauss-Jordan elimination procedure.
Once this is done, then we now have the matrix $Q$ up to a scalar multiple, so we can diagonalize it, into
$ Q = R^T D R $
where $D_{33} = 0 $, then the axis of the paraboloid will be the third column of $R$.
Actually, we can identify all the features of the paraboloid from the parameter vector $\theta$. Let
$Q_1 = \begin{bmatrix} A && \dfrac{D}{2} && \dfrac{E}{2} \\ \dfrac{D}{2} && B && \dfrac{F}{2} \\ \dfrac{ E}{2 } && \dfrac{F}{2} && C \end{bmatrix} $
$b_1 = \begin{bmatrix} G \\ H \\ I \end{bmatrix} $
$c_1 = J $
The equation obtained from these three matrices after scaling by a scalar $t$ are
$ t r^T Q_1 r + t r^T b_1 + t c_1 = 0 \hspace{25pt} (1) $
We will identify $Q_0, r_0$ starting from equation $(*)$ repeated below
$ (r - r_0)^T Q (r - r_0) + b^T (r - r_0) = 0 \hspace{20pt} (*) $
which when expanded gives
$ r^T Q r + r^T (- 2 Q r_0 + b ) + (r_0^T Q r_0 - b^T r_0) = 0 \hspace{25pt} (2) $
Then, we directly have
$ Q = t Q_1 $
however $t$ is unknown yet. Since $Q$ is a scalar multiple of $Q_1$ which is known, then by diagonalizing $Q_1$ we get
$Q_1 = R D R^T $
And therefore
$ Q = R Q_0 R^T = R (t D) R^T $
i.e. Q_0 = t D
Moving on to the linear term, we have
$- 2 Q r_0 + b = t b_1 $
upon substituting $Q = t R D R^T $ and $ b = R b_0 $ this becomes
$ - 2 t R D R^T r_0 + R b_0 = t b_1 \hspace{25pt} (3)$
Pre-multiply through by $R^T$
$ - 2 t D R^T r_0 + b_0 = t R^T b_1 \hspace{25pt} (4)$
Looking at the third row of this vector equation, and since the third row of $D$ is zero (because $D_{33} = 0 $), then the third row says
$ -1 = t \bigg[ R^T b_1 \bigg]_z $
And this gives us the value of $t$, and this means we now have the value of matrix $Q = t Q_1$ and $ b = R b_0 $
Using the top two rows of equation $(4)$ gives
$ -2 D' R^t r_0 = R^T b_1 $
where
$ D' = \begin{bmatrix} D_{11} && 0 && 0 \\ 0 && D_{22} && 0 \end{bmatrix} $
Let vector $u = R^T r_0$, and vector $v = R^T b_1$, then
$ - 2D' u = v $
has infinite solutions, $u = [ -\dfrac{1}{2 D_{11}} v_1 , - \dfrac{1}{2 D_{22}} v_2 , s ]^T $
where $s$ is to be determined.
This means that
$ r_0 = R u = V_0 + s V_1 $
where $V_1$ is along the axis of the paraboloid which is the third column vector of $R$.
Moving on to the third condition,
$ (r_0^T Q r_0 - b^T r_0) = t c_1 \hspace{25pt} (5)$
$ b^T r_0 = (R b_0)^T = b_0^T R^T r_0 = b_0^T u = - s$
$ r_0^T Q r_0 = u^T D u $ which because of the structure of $D$ is not a function of $s$. Therefore, we have a single solution of $(5)$.
Thus, we have identified all the parameters of the elliptical paraboloid, which are the matrix $Q_0$, the vertex $r_0$, and the rotation matrix $R$.