https://orcid.org/0000-0001-9998-3016
Figures
Abstract
During the iterative process of the progressive iterative approximation, it is necessary to calculate the difference between the current interpolation curve and the corresponding data points, known as the adjustment vector. To achieve more precise adjustments of control points, this paper decomposes the adjustment vector into its coordinate components and introduces a weight for each component. By dynamically adjusting these weights, we can accelerate the convergence of iterations and enhance approximation accuracy. During the iteration, the weight coefficients are flexibly adjusted based on the error of the current iteration step, demonstrating the flexibility and precision of the geometric iterative method in addressing geometric approximation problems. Numerical experiment results indicate that this vector decomposition technique is a critical mathematical operation for improving algorithm efficiency and precisely adjusting the shapes of curves or surfaces to approximate the data set.
Citation: Liu Y, Wang Y, Liu C (2025) Adaptive weighted progressive iterative approximation based on coordinate decomposition. PLoS ONE 20(1): e0317225. https://doi.org/10.1371/journal.pone.0317225
Editor: Madhu Chetty, Federation University Australia, AUSTRALIA
Received: July 21, 2024; Accepted: December 23, 2024; Published: January 29, 2025
Copyright: © 2025 Liu et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the article and its Supporting information files.
Funding: This study was funded by National Natural Science Foundation of China (No. 12101225), Natural Science Foundation of Hunan Province (No. 2023JJ50080).
Competing interests: The authors have declared that no competing interests exist.
1 Introduction
Progressive iterative approximation (PIA) is an interpolation-based geometric iterative method. The basic idea of PIA is to continuously adjust the control points of curves and surfaces so that the resulting limit curves and surfaces interpolate the given set of data points [1, 2]. In the iteration of PIA, new control points are calculated based on the current control points and the data points to be interpolated. This process is repeated until the generated curves or surfaces meet specific accuracy requirements or reach a predetermined number of iterations. This method has clear geometric significance and is easy to understand and implement. More importantly, since it adopts an iterative approach, it avoids constructing complex linear systems, thereby reducing the complexity of the solution. In certain scenarios requiring high-precision fitting, interpolation-based geometric iterative methods have significant advantages and wide applications, see [2].
In recent years, researchers have made various improvements to the PIA method for data interpolation to enhance its computational efficiency and accuracy. For example, in [3], the introduction of weights led to the weighted PIA method, known as WPIA; the PIA with multiple weights [4, 5] uses different weights for different adjustment vectors, allowing more flexible handling of local detail features of curves/surfaces and enhance the accuracy; the Chebyshev polynomial accelerated iterative method [6], the preconditioned PIA method [7, 8], and the PIA method with shape parameter curves [9], the iterative construction methods based on the splitting of the collocation matrix, including Jacobi–PIA [10, 11], GS–PIA [12–14], and SOR–PIA [15], as well as other methods [16], etc.
In the aforementioned research findings on the acceleration of geometric iterative methods, some methods are only applicable to specific basis functions. For instance, the method presented in [16] is not ideal in terms of accuracy when approximating large-scale data with Bézier curves; the methods in [7, 8] can only accelerate the PIA for basis functions such as Bernstein and Said-Ball; the methods in [10–15] can only accelerate the PIA method for B-spline curves, etc. It is worth mentioning that research on accelerating convergence by introducing weights is of great significance in geometric iterative methods. Some weighted methods, such as the WPIA method [3] and the PIA with multiple weights [4, 5], are applicable to most PIA methods for curves and surfaces with normalized totally positive basis functions, exhibiting universality. By introducing weight coefficients, the convergence and computational efficiency of the PIA method is improved, as demonstrated in [3–5], and a lot literature therein. Reasonable weights not only increase the algorithm’s convergence rate but also provide more effective tools for handling large-scale data and complex shapes. Despite the significant achievements of weighted PIA methods, there are still some issues that merit further research. For example, how to automatically select and adjust weights to meet different approximation needs, and how to extend this method to more complex geometric shapes. The selection and adjustment strategies for weight coefficients remain key focuses and challenges for future research.
In order to achieve the desired interpolation effect, the control points of the interpolation curve/surface must be continuously adjusted during the iterative process of PIA. This process is not something that can be accomplished in one step, but rather requires consideration of the components of the control vertices along all coordinate axes. However, in-depth research and numerical examples have shown that not all coordinate components need to be continuously adjusted during the adjustment process. In fact, some coordinate components remain relatively stable after several iterations and do not require further adjustment. This finding is of great significance for improving the efficiency and accuracy of iterative methods, as it allows for a greater focus on those coordinate components that truly require optimization, thereby accelerating the entire interpolation process and enhancing the quality of the results.
To rapidly adjust the control points of the interpolating curve, this paper will delve into the adjustment vector weight coefficients in the PIA method. We will adopt a strategy of coordinate decomposition and introduce adjustable weight coefficients, aiming to enhance the flexibility and precision of the PIA method. The organization of this paper is as follows: Section 2 provides a brief overview of the PIA method and its related conclusions; Section 3 presents the adaptive weighted progressive iterative approximation based on coordinate decomposition, applicable to both curve and surface cases; Section 4 demonstrates the effectiveness of the proposed method through numerical experiments. Some conclusions are drawn in the last section.
2 Preliminary
2.1 Progressive iterative approximation
Let a set of data points in
or
to be interpolated, with the i-th point is associated with a parameter ti. Let
be a blending basis that satisfies bi(t) ≥ 0 and
for t ∈ [0, 1].
Initially, we use the interpolation points pi(i = 0 ⋯ n) as the initial control points, that is, . This allows us to derive an initial interpolation curve
given by
We then calculate the difference vector between the current curve and the interpolation points, for i = 0, 1, ⋯, n, and use
as the adjustment vector to adjust the i-th control point, yielding
. Subsequently, using the new set of control points
, a new interpolation curve γ(1)(t) can be generated.
Let the interpolation curve generated after k iterations be denoted as γ(k)(t). The adjustment vectors for the control points can be calculated by
(1)
The interpolation curve generated in the (k + 1)-th iteration is given by
where the updated control points are defined as
(2)
Substituting
from Eq (1) into Eq (2), we get
(3)
Thus, a sequence of sets of control points and a corresponding sequence of curves γ(k)(t), for k = 0, 1, ⋯ can be generated. If the limit
holds true for each i, then the sequence of curves γ(k)(t) is said to have the PIA property, and the method of generating this sequence of curves is referred to as the PIA method. As stated in [1], all normalized totally positive basis functions possess the PIA property.
Let , P = [p0, p1, ⋯, pn]T, and
, then the process of updating control points in Eq (3) can be written in matrix form as
(4)
where I is the identity matrix,
is the collocation matrix of the basis functions
evaluated at tj for j = 0, 1, ⋯, n.
2.2 Weighted progressive iterative approximation
To achieve faster convergence, Lu proposed the weighted progressive iterative approximation (WPIA) method [3]. This method is implemented by multiplying the difference vector by a weight ω where 0 < ω < 2, i.e.,
(5)
Lemma 1 [3] When (6) the WPIA achieves the fastest convergence. In such case, the spectral radius of WPIA is
where λn is the smallest eigenvalue of the collocation matrix B.
3 Adaptive PIA based on coordinate decomposition
3.1 The curve case
If is a set of points in the two-dimensional space
, the control point adjustment vector in Eq (3) can be decomposed along the x-coordinate and y-coordinate directions. If
is a set of points in the three-dimensional space
, the control point adjustment vector in Eq (3) can be decomposed along the x-coordinate, y-coordinate, and z-coordinate directions. For convenience, denote the x-coordinate, y-coordinate, and z-coordinate (if in
) of the point pi to be interpolated as
,
, and
(if in
), respectively. Similarly, denote the x-coordinate, y-coordinate, and z-coordinate (if in
) of the point
in the k-th iteration as
,
, and
. Denote the x-coordinate, y-coordinate, and z-coordinate (if in
) of the control point adjustment vector
as
,
, and
.
Next, we will use the interpolation curve in as an example to analyze the coordinate decomposition process of the PIA method in detail. For the case of
, the derivation process is similar to that of
.
From an algebraic perspective, the iterative process of updating the control points in the PIA method can be viewed as separate iterations in the x− and y− directions, i.e.,
(7)
From a geometric perspective, Fig 1 shows the decomposition diagram of the control point adjustment vectors along the x− and y− directions.
To achieve fast and precise adjustments of control points for the interpolation curve, we in this paper introduce an innovation upon the foundation of Eq (7). We assign unique weight coefficients to the adjustment vectors in each direction, thereby providing more flexible control. The improved iterative process can be expressed as follows:
(8)
where ω(k,x) and ω(k,y) are adjustable weight coefficients dependent on the iteration step k. We refer to the enhanced iterative process (8) as the weighted progressive iterative approximation based on coordinate decomposition (WPIA-CD). By introducing direction-specific weight coefficients, this method enables more precise and adaptable control of the interpolating curve.
Remark 1 It is straightforward to verify that the WPIA-CD method simplifies to the classical PIA method when ω(k,x) = ω(k,y) = 1. Furthermore, when , the WPIA-CD method reduces to the WPIA method.
3.2 Selection of adaptive weight coefficients
In this subsection, we provide an in-depth analysis of the methodology for selecting the weight coefficients ω(k,x) and ω(k,y). For the sake of brevity, this paper presents a detailed exposition of the derivation process for the selection of ω(k,x). Through this process, a more comprehensive understanding of how to reasonably set the weight coefficients to optimize the adjustment effect of the interpolating curve can be achieved.
Consider the matrix form of the iterative process in the x-direction
(9)
where
(10)
is the component of the control point adjustment vector Δ(k) along the x-coordinate, which is also the component of the difference vector between the current interpolating curve and the data points along the x-coordinate. From (9) and (10), the component of the difference vector in the next iteration along the x-coordinate can be expressed as
The weight coefficient ω(k,x) in (9) can be determined by minimizing the square of the l2 norm of Δ(k+1,x), that is
Let
, then we have
Setting the gradient of f(ω) to zero, we can obtain the optimal weight coefficient for updating the control points in the x−direction.
(11)
Similarly, the optimal weight coefficient for updating the control points in the y−direction, ω(k,y), can be expressed as
(12)
From (11) and (12), it is evident that the optimal weight coefficients for the (k + 1)-th iteration of control points depend on the norms of the adjustment vectors of each component in the k-th iteration. Dynamically adjusting the weight coefficients based on the convergence status during the iterative process allows the algorithm to better adapt to different iterative stages and data characteristics, thereby further improving the convergence speed.
Remark 2 In the existing weighted methods, to ensure the fastest convergence speed of the iterative method, it is usually necessary to compute the eigenvalues of the collocation matrix, as mentioned in [3, 11]. Although there are practical methods to solve for the weight coefficients, these methods inevitably affect the convergence speed of the iterative method. In [15], the authors employ optimization algorithms to find approximately optimal weight coefficients. Additionally, some studies, such as [4], do not provide specific theoretical optimal weight coefficients. However, for the method mentioned above, it is only necessary to calculate the norms of the difference vectors, which undoubtedly significantly reduces the computational load.
To ensure the efficiency of the iterative algorithm, a strategy is implemented in the iterations of the control points in the x− and y− directions: during the iterative process, if the iteration error in any direction meets the termination condition (i.e., or
, where tol is a pre-established threshold), the iteration in that direction will immediately stop. Under this strategy, when the weight coefficients are determined according to Eqs (11) and (12) respectively, the iterative process (8) is named the adaptive progressive iterative approximation with coordinate decomposition (APIA-CD) method. The APIA-CD method for curve interpolation in a 2D space can be summarized in the following algorithm.
Algorithm 1 Adaptive progressive iterative approximation with coordinate decomposition(APIA-CD)
1: Input: Set of interpolation points: , pre-established threshold tol, maximum number of iterations ITmax.
2: Set and generate the initial interpolation curve.
3: Set k = 1, calculate and
.
4: while k < ITmax do
5: If , stop the iteration of the x-coordinates of the control points.
6: If , stop the iteration of the y-coordinates of the control points.
7: Calculate ω(k,x) and ω(k,y) according to Eqs (11) and (12), respectively.
8: Update the x- and y-coordinates of the control points according to Eq (8).
9: Calculate and
.
10: k = k + 1.
11: end while
12: Output: The control points at iteration k, number of iterations, and interpolation error.
3.2.1 Convergence analysis.
Theorem 1 The APIA-CD method for curves corresponding to normalized totally positive basis functions is convergent.
Proof Consider the iterative process defined by the APIA-CD method. Specifically, the weight coefficients in the x-direction is chosen to minimize the squared l2 norm of Δ(k+1,x), which is a measure of the error in the x-direction, i.e.,
where Δ(k+1,x) represents the x-component of the difference vector between the current interpolating curve and the data points after k + 1 iterations.
According to Remark 1, the choice of ω(k,x) ensures that the l2 norm of Δ(k+1,x) obtained by the APIA-CD method is less than or equal to the l2-norm of Δ(k+1,x) obtained by the WPIA method, as described in Lemma 1. This implies that the APIA-CD method reduces the error in the x-direction at least as effectively as the WPIA method.
Furthermore, as established in [1], the WPIA method is convergent for curves corresponding to normalized totally positive basis functions. Since the APIA-CD method reduces the error at least as effectively as the WPIA method in each iteration, it follows that the APIA-CD method must also converge for curves corresponding to normalized totally positive basis functions.
3.3 The surface case
Given a set of data points to be interpolated, each point pij is associated with a pair of parameters (ui, vj). Using these points as control points, an initial interpolation surface can be defined as
where
and
denote two sets of basis functions in the u and v parametric directions, respectively.
Let S(k)(u, v) denote the interpolation surface that approximates the data points after k iterations (k = 0, 1, ⋯). The adjustment vector at the (i, j)-th data point for the (k + 1)-th iteration is given by
By decomposing
along the coordinate axes, we obtain the components
,
and
, respectively. Consequently, the coordinates of the control points
for the subsequent iteration of the interpolation surface are updated as
and
In this way, a sequence of surfaces S(k)(u, v) that successively approximate the data points
can be constructed for k = 0, 1, ⋯. This iterative methodology for generating sequences of interpolation surfaces is referred to as the weighted progressive iterative approximation based on coordinate decomposition (WPIA-CD) for surface cases.
Next, we discuss the selection of the weights ω(k,x), ω(k,y) and ω(k,z). Again, only the derivation of the choice of ω(k,x) is described.
First, we organize the interpolation points, control points, and the components of the adjustment vectors along the x-coordinate into vector forms
and
Note that Δ(k,x) can be equivalently expressed as
where ⊗ denotes the Kronecker product of matrices,
and
are the collocation matrices of the basis functions
and
evaluated at parameters uj for j = 0, 1, ⋯, m and vj for j = 0, 1, ⋯, n, respectively.
The iterative process of updating of control points along the x-direction can be expressed as a matrix form
(13)
where I is a identity matrix of the same order as B1 ⊗ B2.
The adjustment vector of the (k + 1)-th iteration along the x-direction can be expressed as
(14)
In analogy with the curve interpolation, the weight ω(k,x) in (14) is ascertained by minimizing the square of the l2-norm of Δ(k+1,x), i.e.,
(15)
Letting
, we have
(16)
Setting the gradient of g(ω) to 0, the optimal weight ω(k,x) for updating the control points in the x-direction can be given by
(17)
Similarly, the optimal weight coefficients for updating the control points in the y- and z- directions can be respectively expressed as
and
To guarantee the computational efficiency of the surface interpolation algorithm, the iterative processes along the x-, y-, and z- directions of updating the control points are stopped as soon as the error associated with the iteration in any given direction falls below a pre-established threshold. The optimal weights are dynamically adjusted depending on the norms of adjustment vectors at the current iteration. This adaptive approach is referred to as the adaptive progressive iterative approximation based on coordinate decomposition (AWPIA-CD) for surface interpolation.
4 Numerical results
In this section, a series of numerical examples are presented to assess the efficiency of APIA-CD in dealing with geometric approximation problems. In order to comprehensively evaluate the performance of APIA-CD in different application scenarios, the experiments cover both uniformly and non-uniformly sampled data sets in and
, including the application of APIA-CD to Bézier curves (surfaces) and quasi-uniform cubic B-spline curves (surfaces).
To evaluate the accuracy of the approximate interpolation algorithms for curve and surface cases, we respectively use
to measure the interpolation error at the k-the iteration.
For comparison, we also tested some other variants of PIA with weights, including the WPIA method [3] and the DWPIA method [4]. All numerical experiments were performed on a computer with an Intel(R) Core(TM) i7–8565U CPU @ 1.80GHz 1.99 GHz processor and 8.00 GB of RAM.
4.1 Numerical results of curve interpolation
Example 1 30 points are sampled from the function where 0 ≤ x ≤ 2π.
Example 2 49 data points are sampled uniformly from the four leaf clover function:
Example 3 53 data points from the parametric equations
For curve interpolation, we employed WPIA, DWPIA, and APIA-CD to iteratively interpolate the given data points in Examples 1–3. In our experiments, we individually assessed the convergence behavior for different coordinate components. To visually demonstrate the performance of the algorithms, Figs 2 to 4 show the l2-norm of each component of the difference vectors as a function of the number of iterations. The results showed that the norm of each component of the difference vector exhibited varying convergence behaviors. This variation in convergence rates is expected and can be attributed to several factors. Firstly, different coordinate components may correspond to different directions along the curve, and the distribution of control points in these directions may be uneven, leading to differential sensitivity to adjustments. Secondly, the interpolation algorithm may be more sensitive to certain components, necessitating more frequent updates for these parts. Furthermore, the strategy for selecting weights can significantly influence the convergence rates of the components. For instance, components with higher weights will undergo more significant adjustments during the iterations. The variation in convergence behavior is crucial for the design and optimization of the iterative method because it allows the algorithm to adjust weights flexibly according to the actual convergence of each component, thereby more effectively controlling the approximation process and ensuring that the final result meets accuracy requirements and converges within a reasonable number of iterations.
(a) Bézier surface. (b) Cubic B-spline surface.
(a) Bézier surface. (b) Cubic B-spline surface.
(a) Bézier surface. (b) Cubic B-spline surface.
The proposed algorithm takes into account both the overall approximation effect and the convergence characteristics of each component. By flexibly adjustment, the performance of the approximation algorithm can be optimized to achieve more efficient and accurate geometric modeling. Moreover, whether for Bézier curves or B-spline curves, the APIA-CD method has shown faster rapid convergence, which is due to the dynamic adjustment of weights during the iterative process of APIA-CD, thereby more effectively controlling the approximation process.
By employing the APIA-CD method, the point sets presented in Examples 1 to 3 are approximately interpolated using Bézier curves and cubic B-spline curves. These point sets and their corresponding interpolation curves are illustrated in Figs 5 to 10, showcasing the progression of curves after various iterations. Within Figs 5 to 10, the blue points denote the interpolation points, the black curves represent the interpolated curves, and the green lines indicate the control polygons. To maintain clarity and simplicity in the figures, a legend is included only in Fig 5, with the understanding that the same color-coding applies to the subsequent figures. To complement these visual representations, the interpolation errors are detailed below their respective figures. It is evident that the curves derived from APIA-CD progressively approximate the provided data points with increasing accuracy.
(a) . (b)
. (c)
.
(a) . (b)
. (c)
.
(a) . (b)
. (c)
.
(a) . (b)
. (c)
.
(a) . (b)
. (c)
.
(a) . (b)
. (c)
.
4.2 Numerical results of surface interpolation
Example 4 Consider interpolating 20 points: (1, 1, 1), (1, 2, 4), (1, 3, 2), (1, 4, 2), (1, 5, 2), (2, 1, 2), (2, 2, 2), (2, 3, 3), (2, 4, 6), (2, 5, 4), (3, 1, 3), (3, 2, 2), (3, 3, 4), (3, 4, 4), (3, 5, 3), (4, 1, 4), (4, 2, 6), (4, 3, 1), (4, 4, 1), (4, 5, 2).
Example 5 Consider interpolating 40 × 42 points sampled from the peaks function.
For surface interpolation, we we have also employed WPIA, DWPIA, and AWPIA-CD to iteratively interpolate the given data points in Examples 4 and 5. Similar to the curve case, we individually assessed the convergence behavior of different coordinate components in our experiments and visually demonstrated the performance of the algorithms in Figs 11 and 12. We observed that the convergence behaviors of the norms of different coordinate components also vary in surface interpolation. Specifically, the norms of the components of the difference vector along the x− and y− directions rapidly approach machine precision after the first iteration, indicating a swift convergence of the algorithm. During the iterative process, it is only necessary to adjust the z-coordinate of the control points, which will obviously reduce the computational load. Moreover, we listed in Tables 1 and 2 the number of iterations, the interpolation errors, and the runtime of WPIA, DWPIA, and AWPIA-CD. From Tables 1 and 2, we can see that under the requirement of the same accuracy, AWPIA-CD requires fewer iterations than WPIA and DWPIA, and the runtime of AWPIA-CD are less than those of WPIA and DWPIA. These results highlight the efficiency and effectiveness of AWPIA-CD in terms of iteration count and runtime.
(a) Bézier surface. (b) Bi-cubic B-spline surface.
(a) Bézier surface. (b) Bi-cubic B-spline surface.
By implementing APIA-CD, the point sets in Examples 4 and 5 are approximately interpolated by Bézier surfaces and bi-cubic B-spline surfaces. Figs 13 and 14 display the point sets and interpolation surfaces for Examples 4 and 5. It can be observed that both the Bézier and bi-cubic B-spline surfaces generated by AWPIA-CD interpolate the given point sets well.
(a) Bézier surface . (b) Bi-cubic B-spline surface
.
(a) Bézier surface . (b) Bi-cubic B-spline surface
.
5 Conclusions
This paper studies and optimizes the adjustment process of control vertices in the interpolated geometric iterative method. To improve control vertex adjustment, the strategy of coordinate decomposition is used. This involves projecting the adjustment vectors to each coordinate axis separately and obtaining the components using vector operations. This paper also introduces adjustable weight coefficients for each component. Adjusting these coefficients speeds up convergence and improves accuracy. The method is flexible and can be easily adapted to different situations. The method in this paper can be used for most all-positive basis functions. These features show how flexible the method is and how accurate it is.
Experiments show that the technique improves the efficiency of the algorithm and accurately adjusts curves or surfaces to approximate data points. This research helps develop and use the geometric iterative method more widely.
Supporting information
S1 Data. The data for Example 5 are sampled from the peaks function.
https://doi.org/10.1371/journal.pone.0317225.s001
(RAR)
References
- 1. Lin HW, Bao HJ, Wang GJ. Totally positive bases and progressive iteration approximation. Computers & Mathematics with Applications. 2005;50(3-4):575–586.
- 2. Lin H, Maekawa T, Deng C. Survey on geometric iterative methods and their applications. Computer-Aided Design. 2018;95:40–51.
- 3. Lu L. Weighted progressive iteration approximation and convergence analysis. Computer Aided Geometric Design. 2010;27(2):129–137.
- 4. Zhang L, Zhao L, Tan J. Progressive iterative approximation with different weights and its application. J Zhejiang Univ(Sci Edn). 2017;44(1):22–27.
- 5. Li Z, Zhonghua L, Lin Z, Xiangrong S, Jieqing T. Local interpolation type of geometric iterative method with multiple weights. Journal of Computer-Aided Design & Computer Graphics. 2018;30(9):1699–1704.
- 6. Liu C, Han X, Li J. The Chebyshev accelerating method for progressive iterative approximation. Communications in Information and Systems. 2017;17(1):25–43.
- 7. Liu C, Liu Z. Progressive iterative approximation with preconditioners. Mathematics. 2020;8(9):1503.
- 8. Liu C, Liu Z, Han X. Preconditioned progressive iterative approximation for tensor product Bézier patches. Mathematics and Computers in Simulation. 2021;185:372–383.
- 9. Chengzhi L, Xuli H, Juncheng L. Progressive-Iterative Approximation by Extension of Cubic Uniform B-spline Curves. Journal of Computer-Aided Design & Computer Graphics. 2019;31(6):899–910.
- 10. Xiaoyan L, Chongyang D. Jacobi-PIA algorithm for non-uniform cubic B-spline curve interpolation. Journal of Computer-Aided Design & Computer Graphics. 2015;27(3):485–491.
- 11. Liu C, Li J, Hu L. Jacobi–PIA algorithm for bi-cubic B-spline interpolation surfaces. Graphical Models. 2022;120:101134.
- 12. Liu X, Deng C. Non-uniform cubic B-spline curve interpolation algorithm of GS-PIA. J Hangzhou Dianzi Univ Nat Sci. 2015;35:79–82.
- 13. Zhihao W, Yajuan L, Chongyang D. Convergence proof of GS-PIA algorithm. Journal of Computer-Aided Design & Computer Graphics. 2018;30(11):2035–2041.
- 14. Xiang Y, Liu C. GS-PIA Algorithm for Bi-cubic B-spline Interpolation Surfaces. IAENG International Journal of Applied Mathematics. 2024;54(6).
- 15. Shou H, Hu L, Fang S. Progressive iterative approximation of non-uniform cubic B-spline curves and surfaces via successive over-relaxation iteration. Mathematics. 2022;10(20):3766.
- 16. Carnicer JM, Delgado J, Peña JM. On the progressive iteration approximation property and alternative iterations. Computer aided geometric design. 2011;28(9):523–526.