Fitting Implicit Surfaces to Oriented Point Sets

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?