15
$\begingroup$

I have a question regarding quadric fit to a set of points and corresponding normals (or equivalently, tangents). Fitting quadric surfaces to point data is well explored. Some works are as follows:

Fitting to projective contours is also covered by some works, such as this one.

From all these works, I think Taubin's method for Quadric fitting is pretty popular:

Let me briefly summarize. A Quadric $Q$ can be written in the algebraic form: $$ f(\mathbf{c},\mathbf{x}) = A x^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz + J $$ where $\mathbf{c}$ is the coefficient vector and $\mathbf{x}$ are the 3D coordinates. Any point $\mathbf{x}$ lies on the quadric $Q$ if $\begin{pmatrix}\mathbf{x}^T & 1\end{pmatrix}Q\begin{pmatrix}\mathbf{x} \\ 1\end{pmatrix}=0$, where: $$ Q = \begin{bmatrix} A & D & E & G \\ D & B & F & H \\ E & F & C & I \\ G & H & I & J \\ \end{bmatrix} $$

Algebraic Fit In principle, we would like to solve for the parameters that minimize the sum of squared geometric distances between the points and the quadratic surface. Unfortunately, it turns out that this is a non-convex optimization problem with no known analytical solutions. Instead, a standard approach is to solve for an algebraic fit, that is to solve for the parameters $\mathbf{c}$ that minimize:

$$ \sum\limits_{i=1}^{n} f(\mathbf{c},\mathbf{x}^i)^2 = \mathbf{c}^T M \mathbf{c} $$ with $$ M = \sum\limits_{i=1}^{n} l(\mathbf{x}^i)l(\mathbf{x}^i)^T $$ where $\{\mathbf{x}^i\}$ are the points in the point cloud and $$ l = [x^2, y^2, z^2, xy, xz, yz, x,y,z, 1]^T $$

Notice that such direct minimization would yield the trivial solution with $\mathbf{c}$ at the origin. This question has been studied extensively in the literature. One resolution that has been found to work well in practice is Taubin’s method (cited above), introducing the constraint:

$$ \| \nabla_x f(\mathbf{c},\mathbf{x}^i) \|^2 = 1 $$

This can be solved as follows: Let : $$ N = \sum\limits_{i=1}^n l_x(\mathbf{x}^i)l_x(\mathbf{x}^i)^T+ l_y(\mathbf{x}^i)l_y(\mathbf{x}^i)^T + l_z(\mathbf{x}^i)l_z(\mathbf{x}^i)^T $$ where subscripts denote the derivatives. The solution is given by the generalized Eigen decomposition, $(M −\lambda N) \mathbf{c} = 0$. The best-fit parameter vector is equal to the Eigenvector corresponding to the smallest Eigenvalue.

Main Question In many applications, the normals of the point cloud are available (or computed). The normals of the quadric $\mathbf{N}(x)$ can also be calculated by differentiating and normalizing the implicit surface:

$$ \mathbf{N}(x) = \frac{\nabla f(\mathbf{c},\mathbf{x})}{\|\nabla f(\mathbf{c},\mathbf{x})\|} $$ where $$ \nabla f(\mathbf{c},\mathbf{x}) = 2 \begin{bmatrix} Ax + Dy + Fz + G \\ By + Dx + Ez + H \\ Cz + Ey + Fx + I \\ \end{bmatrix} $$

However, Taubin's method utilizes only the point geometry, and not the tangent space. And I am not aware of many methods, which are suitable for fitting quadrics such that the tangents of the quadric also match the tangents of underlying point cloud. I am looking for potential extensions of the method above, or any other to cover these first order derivatives.

What I would like to achieve is maybe addressed partially in lower dimensional spaces, with more primitive surface (curve) types. For example, fitting lines to image edges, taking into consideration the gradient information is covered here. Fitting planes (a simple type of quadric) to 3D clouds is very common (link 1) or fitting spheres or cylinders can be fit to oriented point sets (link 2). So what I'm wondering is something similar, but the fitted primitive is a quadric.

I would also welcome the analysis of the proposed method such as:

  • What is the minimum number of oriented points required?
  • What are the degenerate cases ?
  • Can anything be said about robustness?

Update: I would like to present a direction to follow. Formally, what I desire to achieve:

$$ \| \nabla f - \mathbf{n} \| = 0 $$ at the point $\mathbf{x}$. Maybe it might be possible to fuse it with Taubin's method to come up with an additional constraint and minimize using Lagrange multipliers?

Improve this question
$\endgroup$
2
  • $\begingroup$ Aren't many of Q's elements mis-positioned in Q? $\endgroup$
    – Museful
    Nov 27, 2019 at 16:07
  • $\begingroup$ You are right, and I have now fixed this. $\endgroup$
    – Tolga Birdal
    Nov 28, 2019 at 16:03

3 Answers 3

Reset to default
6
$\begingroup$

I was surprised for not receiving a satisfactory answer to the question above and my investigations showed me that this, indeed is an unexplored area. Hence, I put some effort in developing solutions to this problem, and published the following manuscripts:

T. Birdal, B. Busam, N. Navab, S. Ilic and P. Sturm. "A Minimalist Approach to Type-Agnostic Detection of Quadrics in Point Clouds." Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. 2018. http://openaccess.thecvf.com/content_cvpr_2018/html/Birdal_A_Minimalist_Approach_CVPR_2018_paper.html

T. Birdal, B. Busam, N. Navab, S. Ilic and P. Sturm, "Generic Primitive Detection in Point Clouds Using Novel Minimal Quadric Fits," in IEEE Transactions on Pattern Analysis and Machine Intelligence. https://arxiv.org/abs/1901.01255

I will briefly touch the main idea here:

This approach is similar to gradient-one fitting ($\nabla 1$). We align the gradient vector of the quadric $\nabla \mathbf{Q}(\mathbf{x}_i)$ with the normal of the point cloud $\mathbf{n}_i \in \mathbb{R}^3$. However, unlike $\nabla 1$-fits, we opt to use a linear constraint to increase the rank rather than regularizing the solution. This is seemingly non-trivial as the vector-vector alignment brings a non-linear constraint either of the form: $$ \begin{align} \frac{\nabla \mathbf{Q}(\mathbf{x}_i)}{\|\nabla \mathbf{Q}(\mathbf{x}_i)\|} - \mathbf{n}_i\, = \,0\,\quad\text{or}\,\quad \frac{\nabla \mathbf{Q}(\mathbf{x}_i)}{\|\nabla \mathbf{Q}(\mathbf{x}_i)\|} \,\cdot\, \mathbf{n}_i\, = \,1. \end{align} $$ The non-linearity is caused by the normalization as it is hard to know the magnitude and thus the homogeneous scale in advance. We solve this issue by introducing a per normal homogeneous scale $\alpha_i$ among the unknowns and write: $$ \begin{align} \nabla \mathbf{Q}( \mathbf{x}_i) = \nabla \mathbf{v}^T_i \mathbf{q} = \alpha_i \mathbf{n}_i \end{align} $$ where $$ \mathbf{v} = \begin{bmatrix} x^2 & y^2 & z^2 & 2xy & 2xz & 2yz & 2x & 2y & 2z & 1 \end{bmatrix}^T $$ Stacking this up for all $N$ points $\mathbf{x}_i$ and normals $\mathbf{n}_i$ leads to a system of the form $\mathbf{A}^{\prime}\mathbf{q}=\mathbf{0}$: $$ \begin{equation} \begin{bmatrix} \mathbf{v}^T_1 & 0 & 0 & \cdots & 0 \\ \mathbf{v}^T_2 & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathbf{v}^T_n & 0 & 0 & \cdots & 0 \\ {\nabla \mathbf{v}_1^T} & -\mathbf{n}_1 & \mathbf{0}_3 & \cdots & \mathbf{0}_3 \\ {\nabla \mathbf{v}_2^T} & \mathbf{0}_3 & -\mathbf{n}_2 & \cdots & \mathbf{0}_3 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ {\nabla \mathbf{v}_n^T} & \mathbf{0}_3 & \mathbf{0}_3 & \cdots & -\mathbf{n}_n \end{bmatrix} \begin{bmatrix} A \\ B \\ \vdots \\ I \\ J \\ \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_n \end{bmatrix} = \mathbf{0} \end{equation} $$ where ${\nabla \mathbf{v}_i^T}= {\nabla \mathbf{v}(\mathbf{x}_i)^T \in \mathbb{R}^{3\times 10}}$, $\mathbf{0}_3$ is a $3\times 1$ column vector of zeros, $\mathbf{A}^\prime$ is $4N\times(N+10)$ and $\boldsymbol{\alpha}=\{\alpha_i\}$ are the unknown homogeneous scales.

While the solution of this formulation lying in the nullspace of $\mathbf{A}^{\prime}$ produces acceptable results, the system is quite unconstrained in what it could do (the scale factors are too free). It is better to find an appropriate regularizer that is also not too complicated to implement. In practice, once again analogous to gradient-one fitting, we could prefer unit-norm polynomial gradients, and thus, can write $\alpha_i=1$ or equivalently $\alpha_i \gets \bar{\alpha}$, one common scale factor. This soft constraint will try to force zero set of the polynomial respect the local continuity of the data. Such regularization also saves us from solving the sensitive homogeneous system, and lets us re-write the system in a more compact form $\mathbf{A}\mathbf{q}=\mathbf{n}$:

$$ \small \begin{bmatrix} x_1^2 & y_1^2 & z_1^2 & 2x_1y_1 & 2x_1z_1 & 2y_1z_1 & 2x_1 & 2y_1 & 2z_1 & 1\\ x_2^2 & y_2^2 & z_2^2 & 2x_2y_2 & 2x_2z_2 & 2y_2z_2 & 2x_2 & 2y_2 & 2z_2 & 1\\ &&&& \vdots &&&&&\\ 2x_1 & 0 & 0 & 2y_1 & 2z_1 & 0 & 2 & 0 & 0 & 0\\ 0 & 2y_1 & 0 & 2x_1 & 0 & 2z_1 & 0 & 2 & 0 & 0\\ 0 & 0 & 2z_1 & 0 & 2x_1 & 2y_1 & 0 & 0 & 2 & 0\\ 2x_2 & 0 & 0 & 2y_2 & 2z_2 & 0 & 2 & 0 & 0 & 0\\ 0 & 2y_2 & 0 & 2x_2 & 0 & 2z_2 & 0 & 2 & 0 & 0\\ 0 & 0 & 2z_2 & 0 & 2x_2 & 2y_2 & 0 & 0 & 2 & 0\\ &&&& \vdots &&&&& \end{bmatrix} \begin{bmatrix} A \\ B \\ C \\ D \\ E \\ F \\ G \\ H \\ I \\ J \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \dots \\ {n}^1_x \\ {n}^1_y \\ {n}^1_z \\ {n}^2_x \\ {n}^2_y \\ {n}^2_z \\ \dots \end{bmatrix} $$

All in all, solving this system of equations will simultaneously guide the quadric to be incident to the point cloud while aligning its gradients towards the normals. It is also possible to weigh the contributions of points and normals differently. In certain cases, to obtain a type-specific fit, a minor redesign of $\mathbf{A}$ tailored to the desired primitive suffices. For all these details as well as some theoretical analysis and pseudo-code, I refer you to the aforementioned publications.

Improve this answer
$\endgroup$
3
  • $\begingroup$ This is great! How would one modify A to weight the relative contributions of points and normals differently? $\endgroup$
    – Museful
    Nov 27, 2019 at 21:32
  • $\begingroup$ Just multiply the first rows that are the point-equations, with the desired weight. Optionally, to scale the rows corresponding to the normals, one would also need to scale the right hand side of the equation: $\mathbf{n}$. $\endgroup$
    – Tolga Birdal
    Nov 28, 2019 at 15:59
  • $\begingroup$ Thanks. Shouldn't the transpose symbol be removed from q and n in the last equation? $\endgroup$
    – Museful
    Nov 28, 2019 at 20:18
1
$\begingroup$

I know of one example where the normals have been included in the fitting procedure. It is not a direct quadric fitting though. A locally parametrized patch is fitted to the points and normals. Using normals gives more equations in the fitting problem, allowing higher order polynomials to be used.

  1. A novel cubic-order algorithm for approximating principal direction vectors
Improve this answer
$\endgroup$
1
$\begingroup$

See also this paper

Hermite Interpolation of Implicit Surfaces with Radial Basis Functions

and ist extension

Hermite Radial Basis Functions Implicits

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.