Fitting of Analytic Surfaces to Noisy Point Clouds

Fitting C-continuous or superior surfaces to a set S of points sampled on a 2-manifold is central to reverse engineering, computer aided geometric modeling, entertaining, modeling of art heritage, etc. This article addresses the fitting of analytic (ellipsoid, cones, cylinders) surfaces in general position in R3. Currently, the state of the art presents limitations in (i) automatically finding an initial guess for the analytic surface F sought, and (ii) economically estimating the geometric distance between a point of S and the analytic surface F . These issues are central in estimating an analytic surface which minimizes its accumulated distances to the point set. In response to this situation, this article presents and tests novel userindependent strategies for addressing aspects (i) and (ii) above, for cylinders, cones and ellipsoids. A conjecture for the calculation of the distance point-ellipsoid is also proposed. Our strategies produce good initial guesses for F and fast fitting error estimation for F , leading to an agile and robust optimization algorithm. Ongoing work addresses the fitting of free-form parametric surfaces to S.


S
: {p 0 , p 1 ..., p n } Noisy point sample F : best fit surface to S P CA : Principal Component Analysis LM : Lenvenberg-Marquardt k : norm degree f : objective function d i : minimum distance between the i -th point of S and its corresponding point on F

Introduction
Surface reconstruction is a widely studied field because of its importance in CAD-CAM applications, virtual reality, medical imaging and movie industry. Particularly, the reconstruction of analytic surfaces is important since they are frequently used in mechanical parts ( [1]). Surface reconstruction process consists of obtaining a surface that minimizes the distance between each point p i of a point sample S and its corresponding point on surface F . It is assumed that S fulfills the Nyquist-Shannon criteria ( [2], [3]).

Optimization Approach
The optimization problem of fitting F to S is described by the objective function f shown in Eq. (1).
Where the residual d i represents the minimum distance between the i -th point of S and its corresponding point on F and w indicates the order of the residual. Then d i is given by: Where k is the norm-degree to calculate the distance.
To minimize f and find the best surface fit, some variables are tunned. These variables are specific for each situation. Table 1 shows the decision variables for each surface addressed in this paper. On the other hand, norm k remains constant in the optimization process and is considered as a parameter of the problem.
The number of decision variables V and the number of equality constrains E in a optimization problem, allow to know the degrees of freedom with the equiation G = V − E. Table 1 presents the degrees of freedom for each analytic surface addressed in this paper. Notice

Optimization Method
The Gauss-Newton iterative method for solving nonlinear optimization problems uses the Hessian approximation H = J * J T to calculate the next iteration, as is shown in Eq. 3.
Where x is the decision variables vector and r is residuals vector.
Notice that in the case in which f is not strictly convex, J can be singular at some iteration possibly causing the algorithm to diverge. This problem can be overcome by using the Levenberg-Marquardt (LM) Method ( [4], [5]): Where µ is the LM parameter and I is the identity matrix.

Function and Region Convexity
The convexity condition of an objective function and its feasible region determines if a local extrema corresponds to a global extrema or not. In order to evaluate this condition in f , it is required to examine the eigenvalues e of its Hessian matrix H f by solving Eq. 5 If all eigenvalues of H f are positive, then f is strictlyconvex, but if at least one eigenvalue is equal to zero, f is convex ( [6]). In the case studies on this research, an exact calculation of H f is not possible. Thus, a numerical calculation is required by approximating partial derivatives numerically.

Optimality Conditions
Just like Zhou and Salvado ( [8]),Jiang and Cheng ( [9]) classify the surface fitting problem as a non-convex one. The authors do not discuss the convexity analysis neither for the objective function nor for the optimization region. Other references reviewed do not report any classification of the surface fitting problem in terms of the convexity.

Optimization Methods
Yan, Liu and Wang ([10]) use Lloyd iterations to reconstruct quadric surfaces from a 3D point cloud.
Ying, Yang and Zha ( [11]) fit ellipsoids to data using semidefinite programming obtaining low runtime.
Jiang and Cheng ( [9]) apply a decomposition technique to reduce the dimensions of the optimization space. This implies that the possibilities of dropping into local minima decrease. However, as the approach is out of the geometrical field, coming up with an initial guess of the parameters is not an easy task.
Li and Griffiths ( [12]) fit ellipsoids by a least squares method using quadrics. This method does not require an initial estimation but, as it is not based on real geometrical distances, the results do not provide the best geometric fitting ellipsoid.
References [13] and [8] report the use of LM method to fit analytic surfaces without mentioning the selection of the LM parameter and its influence on the optimization process.

Initial Guess
Numerical optimization strategies are sensitive to the initial guess. The closer to the ideal solution is the initial guess, the less number of iterations in the optimization process. On the other hand, a bad initial guess could make the algorithm to diverge.
Just like Ruiz and Cadavid ([14]), Kwon et al. ([15]) use PCA for finding an approximation to axis orientation, center coordinates and radius of point clouds belonging to circular cylindrical surfaces. Because this technique reduces the dimensionality of the data giving the direction of largest dispersion ( [16]), it is limited to cylinders with aspect ratio lower than 5.0 ( [14]) and to cylindrical caps. Similarly, Zhou and Salvado ( [8]) use the eigenvectors of the covariance matrix of the whole data as the axes of the ellipsoid coordinate system. As noted, this technique gives axes that probably will not correspond to the ellipsoid coordinate system in the case of ellipsoidal caps.
Simari and Singh ( [17]) estimate the ellipsoid's center as the geometric centroid of the data set. This proposal works in the cases of ellipsoids completely sampled but it loses validity when the cloud is only a subsample of the whole ellipsoid.
References [18] and [13] report the implementation of an algebraic method based on a least squares solution of the general quadric equation to find an initial guess of the parameters of analytic surfaces.

Literature Review Conclusions and Contribution of this Paper
As was shown in the taxonomy conducted in the literature review, there are several issues that remain open in optimized analytic surfaces fitting which are studied in this work: (1) Estimation of the real geometric distance between a point and an ellipsoid, (2) Identification of the effect of the parameters such as the norm k in the distance measurement, (3) Analysis of the optimality conditions effect on the convergence of the algorithm.

Circular Cylinder Fitting
A circular cylinder is defined by a radius R, an axis vector and its pivot pointv and pv respectively. For purposes of this research, no assumption is made on the orientation or position in space of the cylinder from which the data set belongs.

Initial Guess for Circular Cylinder
The initial guess of the cylinder's parameters is obtained with a statistical and geometrical procedure as explained below and shown in Fig. 1: 1. Random seed points are selected from S and a local neighborhood Ln found for each seed. See Fig. 1(a).
2. Crossing each other the line segments defined by a seed point and the normal vectorn of the best plane fit to Ln, a set of points Cp = {q 1 , q 2 , ..., q n } passing nearv are found. See Fig. 1(b).
3. A PCA is executed over Cp for finding an initial approximation of the cylinder axis and its pivot point, v * and pv * respectively. See Fig. 1(c).
4. An approximation to the cylinder radius is calculated as the average of the minimum distances between S andv * .
This method allows processing both complete cylindrical surfaces and cylindrical caps.

Estimation of Point-Cylinder Distance
The minimum distance d i between the point p i and a circular cylindrical surface is calculated as: Where R is the radius of the cylinder and pr i is the orthogonal projection of p i ontov k as is shown in Fig. 2.

Circular Cone Fitting
A circular right cone can be defined by an axis vectorv, an apex Ap and a semi-opening angle ψ.

Initial Guess for Circular Cone
The initial approximation of a circular conical surface to a point cloud is obtained by an statistical and geometrical procedure as is depicted in Fig. 3 and explained as follows: 1. A set of seed points and local neighborhoods Ln are taken from S. See Fig. 4(a).
2. The minimum curvature directionK min of Ln, being collinear with a generatrix of the cone, is found by fitting a paraboloid p(x, y) = a + bx + cy + dxy + ex 2 + f y 2 and calculating the eigenvectors of it Hes- . See Fig. 4 (a) Seed points and them local neighborhoods Ap * being the first approximation to Ap, is obtained by averaging the crossing points of all lines defined by a seed point and its correspondingK min . Notice that Ap * represents an statistical apex. See Fig.  4(c).
5. The initial estimation of ψ is taken as the average of the angles between the vectorsK min andv * . See Fig. 4(e).

Point-Cone Distance Estimation
The distance d i from a point p i and a cone is calculated by solving the Eq. (7) for α, β and γ ⊂ R.
Whereû is a rotation ofv aroundn 1 an angle −ψ, n 1 =ŵ ×v ŵ×v andn 2 =û ×n1 û×n1 as is shown in the Fig. 5. Finally d i = β is the signed distance between p i and the cone.

Ellipsoid Fitting
As is shown in Table 1, an ellipsoid in general position is defined by the center coordinates (C x , C y , C z ), the semi axes R x , R y , R z and the Euler angles θ x , θ y , θ z .

Initial Guess for Ellipsoid
The first approximation to the parameters of the ellipsoid can be obtained with the following procedure:

A general quadric surface is defined by Eq 8.
a 1 x 2 + a 2 y 2 + a 3 z 2 + a 4 xy + a 5 xy + a 6 yz + a 7 x +a 8 y + a 9 z + a 10 = 0 (8) Rearranging Eq. 8 appropriately ( [18]), it can be written as:   In compact form, the above equation is: If the n points of S are taken into account, Λ is a (9 × n) matrix where its row i-th is: R is a (n × 1) vector with row i-th: x 2 i + y 2 i + z 2 i . B is the coefficients vector and it can be obtained by solving the linear system of Eq. 9. As the system is overspecified, a least squares solution is calculated the with the pseudo-inverse matrix of Λ, Λ + : The initial approximation of the center coordinates C x * , C y * , C z * , the semi axes R x * , R y * , R z * and the Euler angles θ * x , θ * y , θ * z , are obtained from the subdiscriminant A in the matrix notation of a general oriented quadric: The eigenvectors of A represent the axes of the ellipsoid, then θ * x , θ * y and θ * z can be calculated. The eigenvalues of A are proportional to 1 Rx * 2 , 1 Ry * 2 and 1 Rz * 2 ([19]), then R x * , R y * and R z * can be obtained. The initial guess of the ellipsoid center can be obtained as [C x * , C y * , C z * ] = −A −1 V . Figure 6: Calculation of point-ellipsoid distance

Point-Ellipsoid Distance Estimation
We present the following conjecture (Fig. 6) Conjecture. Let E be an ellipsoid centered in (C x , C y , C z ), with Euler angles θ x , θ y , θ z and semi axes R x , R y , R z . Let p i = (x i , y i , z i ) be a point in R 3 . Then, an ellipsoid E exists with the same center and orientation to E, but with semi axes R If E * and p i are translated to the origin and aligned with the principal axes by a rigid transformation, the following equation can be posed: Rearranging (11): The minimum absolute real root of the polynomial in Eq. (12) corresponds to the minimum signed distance from p i to E.

Results and Discussion
In order to test our fitting routines two study cases were proposed as follows.

Data Set 1. Frog
In order to prove the algorithm for fitting ellipsoids, a subset of the frog shown in Fig. 7(a) was taken. In Fig.  7(b) the highlighted points to which the ellipsoids were fitted can be seen. The result of the fitting process is presented in Figs. 7(c) and 7(d), with the history of the fitting error appearing in Fig. 8.
A good initial guess found by an algebraic approach, let to a fast convergence of the algorithm. In Fig. 8 it may be seen that the longest fitting process required of 12 iterations for finding the optimum according to the termination criteria. Fig. 9 shows the initial estimation of the ellipsoid and the best geometrical ellipsoid fit. Notice that the initial surface wraps most of the points, giving a good starting point for the LM algorithm. Table 2 presents a comparison between the the initial ellipsoids and the optimized ones.

Data Set 2. Fan
To test the cylinder and cone fitting algorithms some parts of the fan shown in Fig. 10(a) were processed with     Fig. 10 displays the results of the optimization process of two conical surfaces and one cylinder. As in the case of ellipsoids, the algorithm found the optimal surface after a few iterations. The history of the optimized function for the Fan Data Set is shown in Fig.  11. The cones (Fig. 12) required 4 iterations while the cylinders (Fig. 13) required 6 iterations. The good initial estimation of the surfaces allows the convergence of the algorithm and to reduce the number of iterations, therefore saving computing resources. Geometrical statistics for the Fan Data Set appear in Table 3.

Conclusions and Future Work
This article presents the fitting of analytic surfaces (such as cylinders, cones and ellipsoids) in the sense of mathematical optimization. The objective function for each surface was implemented in terms of the real geometric distance. In the case of cylinder and conical surfaces this metric is formulated and calculated easily. However, in the ellipsoid case the measurement of the distance between a point and the surface is not trivial. In response to this situation this work presented a novel methodology to calculate this distance. The addressed results allow to check that the proposed distance calculation works fine.
The routines for the initial guess of the surfaces were implemented using geometrical and statistical procedures. The study cases allow to prove that the iterative optimization algorithms converge fast with a good initial guess.
Future work includes the extension of the optimization strategies to other analytic and to free form parametric surfaces.