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 $\mathbf{x}^TQ\mathbf{x}=0$, where: $$ Q = \begin{bmatrix} A & B & C & D \\ B & E & F & G \\ C & F & H & I \\ D & G & 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?


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

See also this paper

Hermite Interpolation of Implicit Surfaces with Radial Basis Functions

and ist extension

Hermite Radial Basis Functions Implicits


Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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