Non-linear least squares fitting of Bézier surfaces to unstructured point clouds

1.   Introduction

Bézier curves and surfaces are widely used in computer graphics and computer aided design for representing complex geometries as a set of smooth analytically driven curves. For example, they have been used for aerofoil geometry optimisation ^[1] for representing the outline of textual characters ^[2] and are used extensively for designing and reconstructing complex free form objects in computer aided geometric design ^[3]. A Bézier curve is a parametric curve based on Bernstein polynomials. The shape of the curve is modified by a set of control points, consecutive control points are termed a control polygon, note that the resulting curve does not necessarily pass through these control points. The number of control points defines the possible complexity of the curve, for example, with two control points only a straight line can be evaluated, this being a first degree Bézier curve, with three control points (a second degree Bézier curve) quadratic curvature is possible, and so on. Some example fourth degree (5 by 5 control points) Bézier surfaces are shown in Figure 1 alongside the control points.

In this work we are interested in fitting Bézier surfaces to a set of noisy, unstructured, scattered points, the particular application we are interested in is the removal of geometric form in order to characterise the surface texture of highly complex engineering surfaces ^[4]. It is straightforward to fit a polynomial surface to a set of measured points using the least squares method, many commercial software packages are available to do this fitting task, however Bézier surfaces offer the desirable property of passing through the control points at the corners of its control polygon, this means that Bézier patches can be stitched together to form a patchwork which can be used to approximate geometries of even greater complexity, such as the surface geometries that can now be fabricated using additive manufacturing.

Having reviewed the literature, the linear least squares (LLS) solution for fitting Bézier surfaces is well documented, however, we are unable to find the non-linear least squares (NLLS) Bézier surface fitting algorithm that is presented in this work. Fitting Bézier curves (not surfaces) via LLS and NLLS is considered in references ^[5] and ^[6] and a NLLS spline curve fitting algorithm is presented in ^[7]. In references ^[8,9,10] the LLS Bézier surface fitting algorithm is given, but iterative refinement is achieved using a genetic algorithm, the fireworks algorithm and the firefly algorithm, respectively, rather than NLLS, which makes the fitting task more complicated than it need be for many applications such as form removal for surface texture characterisation. The LLS fitting algorithm for Bézier surfaces is given in references ^[11] and ^[12] but the authors do not present the NLLS fitting algorithm. In ^[13] free form surface fitting is discussed from a reverse engineering perspective but neither the LLS or NLLS Bézier surface fitting algorithms are presented. We therefore adopt the general LLS and NLLS approach presented in references ^[5] and ^[6] but extend to the case of a fourth degree Bézier surface.

Fourth degree Bézier surfaces are considered here to allow for the fitting of highly complex surfaces such as the engineered surface considered at the end of the paper. Linear, quadratic and cubic surfaces are deemed too simple to represent the complex surfaces that are of interest in this work. It is worth noting that the presented algorithm can be modified to consider any degree of Bézier surface.

In this work the partial derivatives that are required for both the LLS and NLLS fitting are evaluated explicitly (see Appendix 2 and 3) as opposed to numerically as per the work in ^[5], this leads to a more computationally efficient implementation.

The structure or order of the data points to which we fit Bézier surfaces is not required apriori, the data points can be randomly scattered. Furthermore, no initial estimate of the control points is required, this is provided by the LLS algorithm. Therefore, the developed method is highly practical and should be of great use for a number of surface fitting tasks.

In order to validate the developed method and benchmark it, we compare the algorithm to the fitting of polynomial surfaces for 3 different surfaces, two simulated surfaces with superimposed noise and one measured surface. The algorithm is benchmarked against polynomials because these are the default fitting functions used by engineers, especially in the field of surface metrology ^[4]. The convergence of the iterative NLLS algorithm is also plotted for each of the considered surfaces to demonstrate the stability of the algorithm.

2.   Method

The general formula for a Bézier surface of degree $n, m$ is:

where $u$ and $v$ are parametric coordinates with $0\le u\le 1$ and $0\le v\le 1$. The term $t$ is the index of the parametric coordinates. The terms ${k}_{ij}$ are the control point coordinates in $x, y, z$ and these are the variables that need to be estimated such that the Bézier surface $P\left({u}_{t}, {v}_{t}\right)$ fits the measured data. We will consider a fourth degree Bézier surface where $m = n = 4$. The fourth degree Bézier surface is written explicitly as per Appendix 1.

Let the measured data we wish to fit the Bézier surface to be denoted ${d}_{t}$ with $t = 1, \dots, N$. The least squares solution requires we minimise the sum of the squared errors. Therefore, define error ${e}_{t}$ as:

Summing the squared errors gives:

We want to minimise $M$, therefore take partial derivatives with respect to ${k}_{ij}$ and set these to zero:

The partial derivatives $\partial {e}_{t}/\partial {k}_{i, j}$ are given in Appendix 2. Taking partial derivatives leads to 25 equations that must be solved simultaneously for the coefficients ${k}_{ij}$. These equations can be written in the form $Ax = b$ and solved in the usual way. First, consider the partial derivatives set equal to zero and then rearranged:

Let the coefficients of the ${k}_{ij}$ terms in $P\left({u}_{t}, {v}_{t}\right)$ be denoted $c1\; to\; c25$:

which takes the form $Ax = b$ where:

We therefore solve for $x$, where $x = {A}^{-1}b$. This gives best fit control points ${k}_{i, j}$ for the initial nodal points ${u}_{t}, {v}_{t}$. Now we need to minimise ${e}_{t}$ by finding the best fit nodal points ${u}_{t}, {v}_{t}$. Given that ${e}_{t}$ is non-linear with respect to ${u}_{t}, {v}_{t}$ we must use non-linear least squares, for this we will use the Newton Raphson method, which states the following:

where ${x}_{n}$ is an initial guess of the root of $f\left(x\right)$ and ${x}_{n+1}$ is an improved guess of the root. For our purpose we want to find the root of ${e}_{t}$ with respect to $u$ and $v$, we therefore write:

for ${u}_{t}$ and ${v}_{t}$ respectively. These equations require that the partial derivatives of ${e}_{t}$ are evaluated with respect to ${u}_{t}$ and ${v}_{t}$, these are given in Appendix 3. The partial derivatives are known as the Jacobian matrices, $J$, and are of size $N$ by $N$ and are filled with zeros apart from the diagonal. Rewriting $\partial {e}_{t}/\partial {u}_{t}$ and $\partial {e}_{t}/\partial {v}_{t}$ as ${J}_{u}$ and ${J}_{v}$ respectively and rearranging to solve for ${u}_{2, t}$ and ${v}_{2, t}$ gives:

where $\alpha$ is a relaxation parameter, $\alpha = 0.5$ is used throughout this work and is chosen empirically. These new estimates of ${u}_{t}$ and ${v}_{t}$ are then used to recalculate the control points ${k}_{i, j}$ using the linear least squares solution previously derived. The above can be implement as an iterative algorithm as follows:

1. Evaluate the linear least squares solution (Eq 7) to calculate the control points ${k}_{i, j}^{x}, {k}_{i, j}^{y}, {k}_{i, j}^{z}$, use an initial estimate of ${u}_{t}, {v}_{t}$ which can be uniformly spaced values from 0 to 1.

2. Evaluate the non-linear least squares solution (Eq 11) to calculate a better estimate of ${u}_{t}, {v}_{t}$.

3. Use the linear least squares solution to calculate new control points ${k}_{i, j}^{x}, {k}_{i, j}^{y}, {k}_{i, j}^{z}$ using the updated estimate of ${u}_{t}, {v}_{t}$.

4. Repeat 2 and 3 until a convergence criterion is met. In this work convergence is deemed to have been met when the percentage difference between the sum of the squared residual (Eq 3) of two successive iterations is less than or equal to 0.5%.

3.   Results

To test the developed algorithm, Bézier surfaces are fitted to unstructured scattered data sets that have been generated by randomly sampling two analytical surfaces, these surfaces are then superimposed with uniformly distributed noise. The analytical surfaces are:

where the $x, y$coordinates are generated by randomly generating uniformly distributed values between -5 and 5. The term $\varepsilon (x, y)$ is a noise signal with uniformly distributed values that lie between -0.1 and 0.1. Both data sets are generated to have 5000 $x, y, z$ coordinates.

To benchmark the Bézier surface fitting algorithm, polynomial surfaces of 5th to 10th degree are fitted to the same scattered data points. The goodness of fit for both the Bézier and polynomial surfaces is calculated as the sum of the squared difference (error) between the $z$ coordinates of the scattered data and the fitted surfaces for all $x, y$ coordinates.

Figure 2A shows the unstructured scattered data points generated using Eq 12, these are represented as coloured dots, and the fitted Bézier surface is represented as a set of continuous red lines. The lines of the Bézier surface have been generated using uniformly spaced values of ${u}_{t}, {v}_{t}$, these are downsampled to show the curvature of the fitted surface more clearly. Figure 2B shows the same data but viewed from the top ($z$) direction, the curvature of the lines illustrates how the values of ${u}_{t}, {v}_{t}$ have been modified by the NLLS algorithm in order for the surface to fit the scattered data, this modification would not occur with the LLS algorithm alone.

Figure 3A shows the convergence of the NLLS fitting algorithm, the convergence is smooth and without oscillation. Also plotted in Figure 3A is a constant value that corresponds to the sum of the squared noise; this is the squared sum of $\varepsilon (x, y)$, the noise residual that would remain after subtracting the perfect analytical form in Eq 12. We would not expect the sum of the squared error of the surface fit to fall below this value, if it did, it would indicate overfitting. Figure 3B shows the sum of the squared error for fitting polynomials of varying degree, the plot shows that the 4th order (25 control points) Bézier surface achieves a fit similar in error to that of a 7th or 8th degree polynomial (36 and 45 fitting coefficients respectively). Also plotted in Figure 3B is the squared sum of $\varepsilon (x, y)$, notice that for polynomials of degree 9 and 10 the fit residual falls below this value, which suggests the higher order polynomials are overfitting.

The results for fitting a Bézier surface to Eq 13 are shown in Figure 4A. As before the unstructured scattered data points are represented as coloured dots, and the fitted Bézier surface is represented as a set of continuous red lines. Figure 4B shows the same data but viewed from the top ($z$) direction, the curvature of the lines illustrates how the values of ${u}_{t}, {v}_{t}$ have been modified by the NLLS algorithm in order for the surface to fit the scattered data.

Figure 5A shows the convergence of the NLLS fitting algorithm when fitting to data sampled from Eq 13, as before the convergence is smooth and without oscillation and the sum of the squared error does not fall below the sum of the squared noise. Fewer iterations are required for the algorithm to convergence when fitting to Eq 13 compared to fitting to Eq 12 (see Figure 3A). Figure 5B shows that the 4th order (25 control points) Bézier surface achieves a fit similar in error to that of a 6th or 7th degree polynomial (28 and 36 fitting coefficients respectively). The sum of the squared error for polynomials of 8th degree and higher converge to the sum of the squared noise, so little to no overfitting is observed.

To test the proposed method with real measured data we consider the problem of fitting and subtracting the geometric form of an engineering surface in order to isolate the surface texture for subsequent surface texture analysis. A complex metallic surface is measured via X-ray computed tomography (XCT), Figure 6A shows a volume rendering of the surface. The surface is extracted from the XCT data as an unstructured point cloud with 14478 $x, yz$ coordinates (Figure 6B) and we fit a Bézier surface using the developed algorithm. The algorithm converges after 75 iterations. The fitted surface is shown in Figure 7A as a set of continuous red lines that are down sampled for clarity. Subtracting the fitted surface from the measured surface yields the waviness and roughness of the surface that can be analysed in the usual way ^[4], see Figure 7B.

4.   Discussion and conclusions

From the examples considered, it can be seen that the NLLS Bézier surface fitting algorithm is able to approximate surfaces of a high degree of complexity and is comparable to a 6th to 8th degree polynomial. The main advantage of using a Bézier surface rather than a polynomial is the ability to form piecewise Bézier patchworks to approximate surfaces of even higher complexity. Bézier surfaces have the highly desirable property of passing through the control points at the corners of the control polygon, see Figure 1. Adjacent Bézier patches can be joined by equating the control points at the joining edges, however, to ensure a smooth transition from one patch to another requires equating the tangents of the joining control points, this is no trivial task. Fitting a patchwork of Bézier surfaces is complicated by the fact that control points at the edge of a patch cannot be changed without influencing the rest of the fitted surface, hence it is desirable to consider the stitching and fitting steps simultaneously. To address this problem Cao et al. ^[15] proposed a method to automatically subdivide a surface into a patchwork, whilst Lin et al. ^[16] proposed an algorithm for adaptively fitting a Bézier patchwork and ensuring first order continuity. Addressing this challenge is beyond the scope of the present work, but will be considered in future work. A comprehensive review of constrained fitting methods can be found in reference ^[17].

In the field of dimensional and surface metrology LLS and NLLS fitting algorithms are well known and accepted for their robustness ^[14]. Although more exotic iterative Bézier surface fitting algorithms have been developed, their complexity is not required for the type of surface fitting considered here, the robustness of the fitting algorithm is of more importance, this being the motivation for developing the algorithm presented here. Furthermore, many of the methods presented in the literature require the coordinate points to be ordered in order to fit a Bézier surface. Our algorithm requires no such apriori information, only the $x, y, z$ coordinates of an unstructured point cloud are required as inputs, such as those generated by optical scanners and X-ray computed tomography.

To avoid overfitting or underfitting, the degree of the Bézier surface should be selected to match the complexity of the measured surface. This work has focused on an engineering surface with two inflection points (a peak and a trough), so at least 4 by 4 control points are required to approximate the surface. It was decided that the number of control points of the Bézier surface should be increased to 5 by 5 in order to better approximate the complexity of the measured surface, thus judgement was exercised in the choice of the degree of the Bézier surface. A more formal approach to selecting the degree of the Bézier surface is to conduct a convergence study, whereby the fit residual is plotted as a function of the degree of the surface, the point at which the fit residual converges is then selected as the appropriate degree of the Bézier surface. Alternatively, a Bézier surface of a fixed degree could be used and the patch size adjusted: if the fit residual is too large then the Bézier patch should be fitted to a smaller region of the measured surface and vice versa.

To conclude, linear and non-linear least squares Bézier surface fitting algorithms have been derived from first principles, the latter has not previously been published for Bézier surfaces. The analytical form of the partial derivatives required for the fitting process have been presented, these were previously evaluated numerically for the case of fitting Bézier curves, not surfaces ^[5]. The performance of the fitting algorithm has been evaluated for simulated and measured data and shown to be suitable for approximating complex surfaces. The developed algorithm has been shown to be stable and to smoothly converge to a solution.

Conflict of interest

All authors declare that there is no conflict of interest in this paper.

Appendix 1: Explicit form of a fourth degree Bézier surface

Appendix 2: Partial derivatives of ${e}_{t}$ with respect to ${k}_{ij}$

Appendix 3: Partial derivatives of ${e}_{t}$ with respect to ${u}_{t}$ and ${v}_{t}$