In mathematics, the progressive-iterative approximation method is an iterative method of data fitting with geometric meanings. Given a set of data points to be fitted, the method obtains a series of fitting curves (or surfaces) by iteratively updating the control points, and the limit curve (surface) can interpolate or approximate the given data points. It avoids solving a linear system of equations directly and allows flexibility in adding constraints during the iterative process. Therefore, it has been widely used in geometric design and related fields.
The study of the iterative method with geometric meaning can be traced back to the work of scholars such as Dongxu Qi and Carl de Boor in the 1970s. In 1975, Qi et al. developed and proved the "profit and loss" algorithm for uniform cubic B-spline curves, and in 1979, de Boor independently proposed this algorithm. In 2004, Hongwei Lin and coauthors proved that non-uniform cubic B-spline curves and surfaces have the "profit and loss" property. Later, in 2005, Lin et al. proved that the curves and surfaces with normalized and totally positive basis all have this property and named it progressive iterative approximation (PIA). In 2007, Maekawa et al. changed the algebraic distance in PIA to geometric distance and named it geometric interpolation (GI). In 2008, Cheng et al. extended it to subdivision surfaces and named the method progressive interpolation (PI). Since the iteration steps of the PIA, GI, and PI algorithms are similar and all have geometric meanings, they are collectively referred to as geometric iterative methods (GIM).
PIA is now extended to several common curves and surfaces in the geometric design field, including NURBS curves and surfaces, T-spline surfaces, and implicit curves and surfaces.
Iteration methods Generally, progressive-iterative approximation (PIA) can be divided into interpolation and approximation schemes. In interpolation algorithms, the number of control points is equal to that of the data points; in approximation algorithms, the number of control points can be less than that of the data points. Specifically, there are some representative iteration methods—such as local-PIA, implicit-PIA, fairing-PIA, and isogeometric least-squares progressive-iterative approximation (IG-LSPIA)—that are specialized for solving the isogeometric analysis problem.
Interpolation scheme: PIA
In interpolation algorithms of PIA, every data point is used as a control point. To facilitate the description of the PIA iteration format for different forms of curves and surfaces, the following formula is uniformly used:
For example: If is a B-spline curve, then is a scalar, is a B-spline basis function, and denotes the control point; If is a B-spline patch with control points, then and , where and are B-spline basis functions; If is a trivariate B-spline solid with control points, then and , where , , and are B-spline basis functions.
Additionally, this can be applied to NURBS curves and surfaces, T-spline surfaces, and triangular Bernstein–Bézier surfaces.
Given an ordered data set with parameters satisfying for , the initial fitting curve is:
where the initial control points of the initial fitting curve can be randomly selected. Suppose that after the th iteration, the th fitting curve is generated by
To construct the st curve, we first calculate the difference vectors,
and use them to update the control points by
which leads to the st fitting curve:
In this way, we obtain a sequence of curves , which converges to a limit curve that interpolates the give data points, i.e.,
Approximation scheme: LSPIA
For the B-spline curve and surface fitting problem, Deng and Lin proposed a least-squares progressive–iterative approximation (LSPIA), which allows the number of control points to be less than the number of the data points and is more suitable for large-scale data fitting problems.
Assume there exists data points and control points, where . Start with equation (), which gives the th fitting curve as
To generate the th fitting curve, first compute the difference vectors for the data points
and then the difference vectors for the control points
where is the index set of the data points in the th group, whose parameters fall in the local support of the th basis function, i.e., . The are weights that guarantee the convergence of the algorithm, usually taken as .
Finally, the control points of the th curve are updated by leading to the th fitting curve . In this way, we obtain a sequence of curve, and the limit curve converges to the least-squares fitting result to the given data points.
Local-PIA
In the local-PIA method, the control points are divided into active and fixed control points, whose subscripts are denoted as and , respectively. Assume that, the th fitting curve is , where the fixed control points satisfy
Then, on the one hand, the iterative formula of the difference vector corresponding to the fixed control points is
On the other hand, the iterative formula of the difference vector corresponding to the active control points is
Arranging the above difference vectors into a one-dimensional sequence,
the local iteration format in matrix form is,
where is the iteration matrix:
where and are the identity matrices and
The above local iteration format converges and can be extended to blending surfaces and subdivision surfaces.
Implicit-PIA
The PIA format for implicit curve and surface reconstruction is presented in the following. Given an ordered point cloud and a unit normal vector on the data points, we want to reconstruct an implicit curve from the given point cloud. To avoid a trivial solution, some offset points are added to the point cloud. They are offset by a distance along the unit normal vector of each point
Assume that is the value of the implicit function at the offset point
Let the implicit curve after the th iteration be
where is the control point.
Define the difference vector of data points as
Next, calculate the difference vector of control coefficients
where is the convergence coefficient. As a result, the new control coefficients are
leading to the new algebraic B-spline curve
The above procedure is carried out iteratively to generate a sequence of algebraic B-spline functions . The sequence converges to a minimization problem with constraints when the initial control coefficients .
Assume that the implicit surface generated after the th iteration is
the iteration format is similar to that of the curve case.
Fairing-PIA
To develop fairing-PIA, we first define the functionals as follows:
where represents the th derivative of the basis function , (e.g. B-spline basis function).
Let the curve after the th iteration be
To construct the new curve , we first calculate the st difference vectors for data points,
Then, the fitting difference vectors and the fairing vectors for control points are calculated by
Finally, the control points of the st curve are produced by
where is a normalization weight, and is a smoothing weight corresponding to the th control point. The smoothing weights can be employed to adjust the smoothness individually, thus bringing great flexibility for smoothness. The larger the smoothing weight is, the smoother the generated curve is. The new curve is obtained as follows
In this way, we obtain a sequence of curves . The sequence converges to the solution of the conventional fairing method based on energy minimization when all smoothing weights are equal (). Similarly, the fairing-PIA can be extended to the surface case.
IG-LSPIA
Isogeometric least-squares progressive-iterative approximation (IG-LSPIA). Given a boundary value problem
where is the unknown solution, is the differential operator, is the boundary operator, and and are the continuous functions. In the isogeometric analysis method, NURBS basis functions are used as shape functions to solve the numerical solution of this boundary value problem. The same basis functions are applied to represent the numerical solution and the geometric mapping :
where denotes the NURBS basis function, is the control coefficient. After substituting the collocation points into the strong form of PDE, we obtain a discretized problem
where and denote the subscripts of internal and boundary collocation points, respectively.
Arranging the control coefficients of the numerical solution into an -dimensional column vector , the discretized problem can be reformulated in matrix form as
where is the collocation matrix and is the load vector.
Assume that the discretized load values are data points to be fitted. Given the initial guess of the control coefficients , we obtain an initial blending function
where , , represents the combination of different order derivatives of the NURBS basis functions determined using the operators and
where and indicate the interior and boundary of the parameter domain, respectively. Each corresponds to the th control coefficient. Assume that and are the index sets of the internal and boundary control coefficients, respectively. Without loss of generality, we further assume that the boundary control coefficients have been obtained using strong or weak imposition and are fixed, i.e.,
The th blending function, generated after the th iteration of IG-LSPIA, is assumed to be as follows:
Then, the difference vectors for collocation points (DCP) in the st iteration are obtained using
Moreover, group all load values whose parameters fall in the local support of the th derivatives function, i.e., , into the th group corresponding to the th control coefficient, and denote the index set of the th group of load values as . Lastly, the differences for control coefficients (DCC) can be constructed as follows:
where is a normalization weight to guarantee the convergence of the algorithm.
Thus, the new control coefficients are updated via the following formula,
Consequently, the st blending function is generated as follows:
The above iteration process is performed until the desired fitting precision is reached and a sequence of blending functions is obtained
The IG-LSPIA converges to the solution of a constrained least-squares collocation problem.
Proof of convergence
Non-singular case Let be the number of control points and be the number of data points.
If , the PIA iterative format in matrix form is
where
The convergence of the PIA is related to the properties of the collocation matrix. If the spectral radius of the iteration matrix is less than , then the PIA is convergent. It has been shown that the PIA methods are convergent for Bézier curves and surfaces, B-spline curves and surfaces, NURBS curves and surfaces, triangular Bernstein–Bézier surfaces, and subdivision surfaces (Loop, Catmull-Clark, Doo-Sabin).
If , the LSPIA in matrix form is
When the matrix is nonsingular, the following results can be obtained:
Proof Since is nonsingular, and , then . Moreover,
In summary, .
Proof From the matrix form of iterative format, we obtain the following:
According to above Lemma, the spectral radius of the matrix satisfies
and thus the spectral radius of the iteration matrix satisfies
When
As a result,
i.e., , which is equivalent to the normal equation of the fitting problem. Hence, the LSPIA algorithm converges to the least squares result for a given sequence of points.
Singular case
Lin et al. showed that LSPIA converges even when the iteration matrix is singular.
Acceleration algorithms and others
Precondition: Liu et al. proposed a preconditioned PIA for Bézier surfaces via the diagonally compensated reduction method, effectively improving the accuracy and efficiency of the classical algorithm.