View source for Progressive-iterative approximation method

{{Short description|Computer-aided geometric design}}

{{Multiple issues|{{Technical|date=August 2024}}
{{Orphan|date=August 2024}}}}
In mathematics, the '''progressive-iterative approximation method''' is an [[iterative method]] of [[Curve fitting|data fitting]] with geometric meanings.<ref name=":0">{{Cite journal |last1=Lin |first1=Hong-Wei |last2=Bao |first2=Hu-Jun |last3=Wang |first3=Guo-Jin |title=Totally positive bases and progressive iteration approximation |journal=Computers & Mathematics with Applications |date=2005 |volume=50 |issue=3–4 |pages=575–586 |doi=10.1016/j.camwa.2005.01.023 |issn=0898-1221}}</ref> 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 (mathematics)|limit]] curve (surface) can [[Interpolation|interpolate]] or [[Approximation theory|approximate]] the given data points.<ref name=":9" /> It avoids solving a linear system of equations directly and allows flexibility in adding constraints during the iterative process.<ref name=":1" /> Therefore, it has been widely used in geometric design and related fields.<ref name=":9" />

The study of the iterative method with geometric meaning can be traced back to the work of scholars such as Dongxu Qi and [[Carl R. de Boor|Carl de Boor]] in the 1970s.<ref name=":10" /><ref name=":11" /> In 1975, Qi et al. developed and proved the "profit and loss" algorithm for uniform cubic [[B-spline]] curves,<ref name=":10">{{Cite journal |last1=Qi |first1=Dongxu |last2=Tian |first2=Zixian |last3=Zhang |first3=Auxin |last4=Feng |first4=Jiabin |year=1975 |title=The method of numeric polish in curve fitting |journal=Acta Math Sinica |volume=18 |issue=3 |pages=173–184}}</ref> and in 1979, de Boor independently proposed this algorithm.<ref name=":11">{{Cite journal |last=Carl |first=de Boor |year=1979 |title=How does Agee's smoothing method work? |journal=Proceedings of the 1979 Army Numerical Analysis and Computers Conference, ARO Report.}}</ref> In 2004, Hongwei Lin and coauthors proved that non-uniform cubic B-spline curves and surfaces have the "profit and loss" property.<ref name=":1">{{Cite journal |last1=Lin |first1=Hongwei |last2=Wang |first2=Guojin |last3=Dong |first3=Chenshi |date=2004 |title=Constructing iterative non-uniform B-spline curve and surface to fit data points |journal=Science in China Series F |volume=47 |issue=3 |pages=315 |doi=10.1360/02yf0529 |s2cid=966980 |issn=1009-2757}}</ref> 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).<ref name=":0" /> In 2007, Maekawa et al. changed the algebraic distance in PIA to geometric distance and named it ''geometric interpolation'' (GI).<ref>{{Cite journal |last1=Maekawa |first1=Takashi |last2=Yasunori |first2=Matsumoto |last3=Ken |first3=Namiki |year=2007 |title=Interpolation by geometric algorithm |journal=Computer-Aided Design |volume=39 |issue=4 |pages=313–323|doi=10.1016/j.cad.2006.12.008 }}</ref> In 2008, Cheng et al. extended it to [[Subdivision surface|subdivision surfaces]] and named the method ''progressive interpolation'' (PI).<ref>{{Cite book |last1=Cheng |first1=Fuhua |last2=Fan |first2=Fengtao |last3=Lai |first3=Shuhua |last4=Huang |first4=Conglin |last5=Wang |first5=Jiaxi |last6=Yong |first6=Junhai |title=Advances in Geometric Modeling and Processing |chapter=Progressive Interpolation Using Loop Subdivision Surfaces  |series=Lecture Notes in Computer Science |date=2008 |volume=4975 |pages=526–533|doi=10.1007/978-3-540-79246-8_43 |isbn=978-3-540-79245-1 }}</ref> 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).<ref name=":9">{{Cite journal |last1=Lin |first1=Hongwei |last2=Maekawa |first2=Takashi |last3=Deng |first3=Chongyang |title=Survey on geometric iterative methods and their applications |journal=Computer-Aided Design |date=2018 |volume=95 |pages=40–51 |doi=10.1016/j.cad.2017.10.002 |issn=0010-4485}}</ref>

PIA is now extended to several common curves and surfaces in the [[geometric design]] field,<ref name=":13">{{Cite book |last=Hoschek |first=Josef |url=https://dl.acm.org/doi/abs/10.5555/174506 |title=Fundamentals of computer aided geometric design |date=February 1993 |publisher=A. K. Peters, Ltd. |isbn=978-1-56881-007-2 |location=USA}}</ref> including [[Non-uniform rational B-spline|NURBS]] curves and surfaces,<ref name=":2">{{Cite journal |last1=Shi |first1=Limin |last2=Wang |first2=Renhong |year=2006 |title=An iterative algorithm of NURBS interpolation and approximation |journal=Journal of Mathematical Research with Applications |volume=26 |issue=4 |pages=735–743}}</ref> [[T-spline]] surfaces,<ref name=":3">{{Cite journal |last1=Lin |first1=Hongwei |last2=Zhang |first2=Zhiyu |title=An Efficient Method for Fitting Large Data Sets Using T-Splines |journal=SIAM Journal on Scientific Computing |date=2013 |volume=35 |issue=6 |pages=A3052–A3068 |doi=10.1137/120888569 |bibcode=2013SJSC...35A3052L |issn=1064-8275}}</ref> and [[implicit curve]]s and [[implicit surface|surfaces]].<ref name=":5">{{Cite journal |last1=Hamza |first1=Yusuf Fatihu |last2=Lin |first2=Hongwei |last3=Li |first3=Zhao |year=2020 |title=Implicit progressive-iterative approximation for curve and surface reconstruction |journal=Computer Aided Geometric Design |volume=77 |pages=101817|doi=10.1016/j.cagd.2020.101817 |arxiv=1909.00551 |s2cid=202540812 }}</ref>

== Iteration methods ==
Generally, progressive-iterative approximation (PIA) can be divided into interpolation and approximation schemes.<ref name=":9" /> 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,<ref name=":6">{{Cite journal |last=Lin |first=Hongwei |date=2010 |title=Local progressive-iterative approximation format for blending curves and patches |journal=Computer Aided Geometric Design |volume=27 |issue=4 |pages=322–339 |doi=10.1016/j.cagd.2010.01.003 |issn=0167-8396}}</ref> implicit-PIA,<ref name=":5" /> fairing-PIA,<ref name=":7">{{Cite journal |last1=Jiang |first1=Yini |last2=Lin |first2=Hongwei |last3=Huang |first3=Weixian |date=2023-05-16 |title=Fairing-PIA: progressive-iterative approximation for fairing curve and surface generation |journal=The Visual Computer |volume=40 |issue=3 |pages=1467–1484 |arxiv=2211.11416 |doi=10.1007/s00371-023-02861-7 |issn=0178-2789}}</ref> and isogeometric least-squares progressive-iterative approximation (IG-LSPIA)<ref name=":8" />—that are specialized for solving the [[isogeometric analysis]] problem.<ref name=":14">{{Cite journal |last1=Hughes |first1=T. J. R. |last2=Cottrell |first2=J. A. |last3=Bazilevs |first3=Y. |date=2005-10-01 |title=Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement |url=https://www.sciencedirect.com/science/article/pii/S0045782504005171 |journal=Computer Methods in Applied Mechanics and Engineering |volume=194 |issue=39 |pages=4135–4195 |doi=10.1016/j.cma.2004.10.008 |bibcode=2005CMAME.194.4135H |issn=0045-7825}}</ref>

=== Interpolation scheme: PIA ===
[[File:截屏2024-08-10 09.23.08.png|thumb|439x439px|Interpolation scheme of PIA<br/>Top left: The data points and the initial control polygon (Here, the initial control points are taken as the data points). Top right: The initial curve and the difference vectors. Bottom left: New control polygon is generated by adding the difference vectors to the old control points. Bottom right: The new control polygon and new curve (in purple). 

]]
In interpolation algorithms of PIA,<ref name=":0" /><ref name=":1" /><ref name=":2" /><ref>{{Cite journal |last1=Chen |first1=Jie |last2=Wang |first2=Guo-Jin |date=2011 |title=Progressive iterative approximation for triangular Bézier surfaces |journal=Computer-Aided Design |volume=43 |issue=8 |pages=889–895 |doi=10.1016/j.cad.2011.03.012 |issn=0010-4485}}</ref> 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:
<math display="block">\mathbf{P}(\mathbf{t})=\sum_{i=1}^n\mathbf{P}_iB_i(\mathbf{t}).</math>
For example:
* If <math>\mathbf{P}(\mathbf{t})</math> is a B-spline curve, then <math>\mathbf{t}</math> is a scalar, <math>B_i(t)</math> is a B-spline basis function, and <math>\mathbf{P}_i</math> denotes the control point;<ref name=":13" />
* If <math>\mathbf{P}(\mathbf{t})</math> is a B-spline patch with <math>n_u\times n_v</math> control points, then <math>\mathbf{t}=(u,v)</math> and <math>B_i(\mathbf{t})=N_i(u)N_i(v)</math>, where <math>N_i(u)</math> and <math>N_i(v)</math> are B-spline basis functions;<ref name=":13" />
* If <math>\mathbf{P}(\mathbf{t})</math> is a trivariate B-spline solid with <math>n_u \times n_v \times n_w</math> control points, then <math>\mathbf{t}=(u,v,w)</math> and <math>B_i(\mathbf{t})=N_i(u)N_i(v)N_i(w)</math>, where <math>N_i(u)</math>, <math>N_i(v)</math>, and <math>N_i(w)</math> are B-spline basis functions.<ref>{{Cite journal |last1=Lin |first1=Hongwei |last2=Jin |first2=Sinan |last3=Hu |first3=Qianqian |last4=Liu |first4=Zhenbao |date=2015 |title=Constructing B-spline solids from tetrahedral meshes for isogeometric analysis |url=http://dx.doi.org/10.1016/j.cagd.2015.03.013 |journal=Computer Aided Geometric Design |volume=35-36 |pages=109–120 |doi=10.1016/j.cagd.2015.03.013 |issn=0167-8396}}</ref>

Additionally, this can be applied to NURBS curves and surfaces, T-spline surfaces, and triangular Bernstein–Bézier surfaces.<ref name=":12" />

Given an ordered data set <math>\mathbf{Q}_i</math> with parameters <math>t_i</math> satisfying <math>t_1<t_2<\cdots</math> for <math>i=1,2,\cdots,n</math>, the initial fitting curve is:<ref name=":0" />
<math display="block">\mathbf{P}^{(0)}(t)=\sum_{i=1}^n\mathbf{P}_i^{(0)}B_i(t)</math>
where the initial control points of the initial fitting curve <math>\mathbf{P}_i^{(0)}</math> can be randomly selected. Suppose that after the <math>k</math>th iteration, the <math>k</math>th fitting curve <math>\mathbf{P}^{(k)}(t)</math> is generated by

{{NumBlk|:|<math> \mathbf{P}^{(k)}(t)=\sum_{i=1}^n\mathbf{P}_i^{(k)}B_i(t).</math>|{{EquationRef|1}}}}

To construct the <math>(k+1)</math>st curve, we first calculate the ''difference vectors'',
<math display="block">\mathbf{\Delta}^{(k)}_i=\mathbf{Q}_i-\mathbf{P}^{(k)}(t_i), \quad 
i=1,2,\cdots,n</math>
and use them to update the control points by
<math display="block">\mathbf{P}_i^{(k+1)}=\mathbf{P}_i^{(k)}+\mathbf{\Delta}_i^{(k)} </math>
which leads to the <math>(k+1)</math>st fitting curve:
<math display="block">
\mathbf{P}^{(k+1)}(t)=\sum_{i=1}^n\mathbf{P}_i^{(k+1)}B_i(t).</math>
In this way, we obtain a sequence of curves <math display="inline">
\mathbf{P}^{(\alpha)}(t),\alpha=0,1,2,\cdots
</math>, which converges to a limit curve that interpolates the give data points,<ref name=":0" /><ref name=":2" /> i.e.,
<math display="block">
\lim \limits_{\alpha\rightarrow\infty}\mathbf{P}^{(\alpha)}(t_i)=\mathbf{Q}_i, \quad 
i=1,2,\cdots,n.
</math>

=== Approximation scheme: LSPIA ===
[[File:截屏2024-08-13 09.51.05.png|thumb|641x641px|Approximation scheme: LSPIA<br />Top left: Data points '''''<math>\mathbf{Q}_i</math>''''' (blue circles), initial control polygon (green lines) constructed from a subset of '''''<math>\mathbf{Q}</math>''''', and initial fitting curve <math>\mathbf{P}^{(0)} (t)</math>. Top right: Difference vectors '''''<math>\boldsymbol{\delta}^{(k)}_i</math>''''' for data points and difference vectors '''''<math>\mathbf{\Delta}^{(k)}_j</math>''''' for control points. Bottom: A new control polygon (purple lines) is generated by adding '''''<math>\mathbf{\Delta}^{(k)}_j</math>''''' to the old control points; it then creates the next fitting curve <math>\mathbf{P}^{(1)} (t)</math> (purple curve). 

]]
For the B-spline curve and surface fitting problem, Deng and Lin proposed a '''least-squares progressive–iterative approximation''' (LSPIA),<ref name=":3" /><ref name=":4">{{Cite journal |last1=Deng |first1=Chongyang |last2=Lin |first2=Hongwei |date=2014 |title=Progressive and iterative approximation for least squares B-spline curve and surface fitting |journal=Computer-Aided Design |volume=47 |pages=32–44 |doi=10.1016/j.cad.2013.08.012 |issn=0010-4485}}</ref> 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.<ref name=":3" />

Assume there exists <math>m</math> data points and <math>n</math> control points, where <math>n\le m</math>. Start with equation ({{EquationNote|1}}), which gives the <math>k</math>th fitting curve as
<math display="block">
\mathbf{P}^{(k)}(t)=\sum_{j=1}^n\mathbf{P}_j^{(k)}B_j(t).
</math>
To generate the <math>(k+1)</math>th fitting curve, first compute the difference vectors for the data points<ref name=":3" /><ref name=":4" />
<math display="block">
\boldsymbol{\delta}^{(k)}_i=\mathbf{Q}_i-\mathbf{P}^{(k)}(t_i), \quad 
i=1,2,\cdots,m
</math>
and then the difference vectors for the control points
<math display="block">\mathbf{\Delta}^{(k)}_j=\frac{
\sum_{i \in I_j}{c_i B_j(t_i) \boldsymbol{\delta}_i^{(k)}}}{\sum_{i \in I_j}c_i B_j(t_i)}, \quad 
j = 1,2,\cdots,n</math>
where <math>I_j</math> is the index set of the data points in the <math>j</math>th group, whose parameters fall in the local support of the <math>j</math>th basis function, i.e., <math>B_j(t_i)\ne0</math>. The <math>c_i</math> are weights that guarantee the convergence of the algorithm, usually taken as <math>c_i = 1, i \in I_j</math>.

Finally, the control points of the <math>(k+1)</math>th curve are updated by <math> \mathbf{P}_j^{(k+1)}=\mathbf{P}_j^{(k)}+\mathbf{\Delta}_j^{(k)}, </math> leading to the <math>(k+1)</math>th fitting curve <math>\mathbf{P}^{(k+1)}(t)</math>. In this way, we obtain a sequence of curve, and the limit curve converges to the [[Least squares|least-squares fitting]] result to the given data points.<ref name=":3" /><ref name=":4" />

=== Local-PIA ===
[[File:截屏2024-08-13 13.45.43.png|thumb|573x573px|Local PIA: If only one control point is adjusted, the Bézier curve just interpolates the data point (in red) corresponding to the adjusted control point. ]]
In the local-PIA method,<ref name=":6" /> the control points are divided into active and fixed control points, whose subscripts are denoted as <math display="inline">I=\left\{i_1,i_2,\cdots,i_I\right\}</math> and <math display="inline">J=\left\{j_1,j_2,\cdots,j_J\right\}</math>, respectively. Assume that, the <math display="inline">k</math>th fitting curve is <math display="inline">\mathbf{P}^{(k)}(t)=\sum_{j=1}^n\mathbf{P}_j^{(k)}B_j(t)</math>, where the fixed control points satisfy
<math display="block">
\mathbf{P}_j^{(k)}=\mathbf{P}_j^{(0)},\quad j\in J,\quad k=0,1,2,\cdots.
</math>
Then, on the one hand, the iterative formula of the difference vector <math display="inline">\mathbf{\Delta}_h^{(k+1)}</math> corresponding to the fixed control points is
<math display="block">\begin{aligned}
\mathbf{\Delta}_h^{(k+1)}&=\mathbf{Q}_h-\sum_{j=1}^n\mathbf{P}_j^{(k+1)}B_j(t_h)\\
&=\mathbf{Q}_h-\sum_{j\in J}\mathbf{P}_j^{(k+1)}B_j(t_h)-\sum_{i\in I}\left(\mathbf{P}_i^{(k)}+\mathbf{\Delta}_i^{(k)}\right)B_i(t_h)\\
&=\mathbf{Q}_h-\sum_{j=1}^n\mathbf{P}_j^{(k)}B_j(t_h)-\sum_{i\in I}\mathbf{\Delta}_i^{(k)}B_i(t_h)\\
&=\mathbf{\Delta}_h^{(k)}-\sum_{i\in I}\mathbf{\Delta}_i^{(k)}B_i(t_h), \quad h\in J.
\end{aligned}</math>
On the other hand, the iterative formula of the difference vector <math display="inline">\mathbf{D}_l^{(k+1)}</math> corresponding to the active control points is
<math display="block">
\begin{aligned}
\mathbf{\Delta}_l^{(k+1)}&=\mathbf{Q}_l-\sum_{j=1}^n\mathbf{P}_j^{(k+1)}B_j(t_l)\\
&=\mathbf{Q}_l-\sum_{j=1}^n\mathbf{P}_j^{(k)}B_j(t_l)-\sum_{i\in I}\mathbf{\Delta}_i^{(k)}B_i(t_l)\\
&=\mathbf{\Delta}_l^{(k)}-\sum_{i\in I}\mathbf{\Delta}_i^{(k)}B_i(t_l)\\
&=-\mathbf{\Delta}_{i_1}^{(k)}B_{i_1}(t_l)-\mathbf{\Delta}_{i_2}^{(k)}B_{i_2}(t_l)-\cdots+\left(1-B_l(t_l)\right)\mathbf{\Delta}_l ^{(k)}-\cdots-\mathbf{\Delta}_{i_I}^{(k)}B_{i_I}(t_l),\quad l\in I.
\end{aligned}
</math>
Arranging the above difference vectors into a one-dimensional sequence,
<math display="block">
\mathbf{D}^{(k+1)}=\left[\mathbf{\Delta}_{j_1}^{(k+1)} ,\mathbf{\Delta}_{j_2}^{(k+1)},\cdots,\mathbf{\Delta}_{j_J}^{(k+1)},\mathbf{\Delta}_{i_1}^{(k+1)},\mathbf{\Delta}_{i_2}^{(k+1)},\cdots,\mathbf{\Delta}_{i_I}^{(k+1)}\right]^T,\quad k=0,1,2,\cdots,
</math>
the local iteration format in matrix form is,
<math display="block">
\mathbf{D}^{(k+1)}=\mathbf{T}\mathbf{D}^{(k)},\quad k=0,1,2,\cdots,
</math>
where <math display="inline">\mathbf{T}</math> is the iteration matrix:
<math display="block">
\mathbf{T}=
\begin{bmatrix}
\mathbf{E}_J  & -\mathbf{B}_1\\
0             & \mathbf{E}_I-\mathbf{B}_2
\end{bmatrix},
</math>
where <math display="inline">\mathbf{E}_J</math> and <math display="inline">\mathbf{E}_I</math> are the identity matrices and
<math display="block">
\mathbf{B}_1=
\begin{bmatrix}
B_{i_1}\left(t_{j_1} \right) & B_{i_2}\left(t_{j_1} \right) & \cdots &B_{i_I}\left(t_{j_1} \right) \\
B_{i_1}\left(t_{j_2} \right) & B_{i_2}\left(t_{j_2} \right) & \cdots &B_{i_I}\left(t_{j_2} \right) \\
\vdots & \vdots &\vdots & \vdots \\
B_{i_1}\left(t_{j_J} \right) & B_{i_2}\left(t_{j_J} \right) & \cdots &B_{i_I}\left(t_{j_J} \right) \\
\end{bmatrix},

\mathbf{B}_2=
\begin{bmatrix}
B_{i_1}\left(t_{i_1} \right) & B_{i_2}\left(t_{i_1} \right) & \cdots &B_{i_I}\left(t_{i_1} \right) \\
B_{i_1}\left(t_{i_2} \right) & B_{i_2}\left(t_{i_2} \right) & \cdots &B_{i_I}\left(t_{i_2} \right) \\
\vdots & \vdots &\vdots & \vdots \\
B_{i_1}\left(t_{i_I} \right) & B_{i_2}\left(t_{i_I} \right) & \cdots &B_{i_I}\left(t_{i_I} \right) \\
\end{bmatrix}.
</math>
The above local iteration format converges and can be extended to blending surfaces<ref name=":6" /> and subdivision surfaces.<ref>{{Cite journal |last1=Zhao |first1=Yu |last2=Lin |first2=Hongwei |last3=Bao |first3=Hujun |date=2012 |title=Local progressive interpolation for subdivision surface fitting |journal=Journal of Computer Research and Development |volume=49 |issue=8 |pages=1699–1707}}</ref>

=== Implicit-PIA ===

The PIA format for implicit curve and surface reconstruction is presented in the following.<ref name=":5" /> Given an ordered point cloud <math display="inline">\left\{\mathbf{Q}_i\right\}_{i=1}^n</math> and a unit normal vector <math display="inline">\left\{\mathbf{n}_i\right\}_{i=1}^n</math> on the data points, we want to reconstruct an implicit curve from the given point cloud. To avoid a trivial solution, some offset points <math display="inline">\left\{\mathbf{Q}_l\right\}_{l=n+1}^{2n}</math> are added to the point cloud.<ref name=":5" /> They are offset by a distance <math display="inline">\sigma</math> along the unit normal vector of each point
<math display="block">
\mathbf{Q}_l=\mathbf{Q}_i+\sigma\mathbf{n}_i,\quad l=n+i,\quad i=1,2,\cdots,n.
</math>
Assume that <math display="inline">\epsilon</math> is the value of the implicit function at the offset point
<math display="block">
f\left(\mathbf{Q}_l\right)=\epsilon,\quad l=n+1,n+2,\cdots,2n.
</math>
Let the implicit curve after the <math display="inline">\alpha</math>th iteration be
<math display="block">
f^{(\alpha)}(x,y)=\sum_{i=1}^{N_u}\sum_{j=1}^{N_v}C_{ij}^{(\alpha)}B_i(x)B_j(y),
</math>
where <math display="inline">C_{ij}^{(\alpha)}</math> is the control point.

Define the difference vector of data points as<ref name=":5" />
<math display="block">
\begin{aligned}
\boldsymbol{\delta}_k^{(\alpha)}&=0-f^{(\alpha)}(x_k,y_k),\quad k=1,2,\cdots,n,\\
\boldsymbol{\delta}_l^{(\alpha)}&=\epsilon-f^{(\alpha)}(x_l,y_l),\quad l=n+1,n+2,\cdots, 2n.
\end{aligned}
</math>
Next, calculate the difference vector of control coefficients
<math display="block">
\boldsymbol{\Delta}_{ij}^{(\alpha)}=\mu\sum_{k=1}^{2n} B_i(x_k)B_j(y_k) \boldsymbol{\delta}_k^{(\alpha)},\quad i=1,2,\cdots,N_u,\quad j=1,2,\cdots,N_v,
</math>
where <math display="inline">\mu</math> is the convergence coefficient. As a result, the new control coefficients are
<math display="block">
C_{ij}^{(\alpha+1)}=C_{ij}^{(\alpha)}+\boldsymbol{\Delta}_{ij}^{(\alpha)},
</math>
leading to the new algebraic B-spline curve
<math display="block">
f^{(\alpha+1)}(x,y)=\sum_{i=1}^{N_u}\sum_{j=1}^{N_v}C_{ij}^{(\alpha+1)}B_i(x)B_j(y).
</math>
The above procedure is carried out iteratively to generate a sequence of algebraic B-spline functions <math display="inline">\left\{f^{(\alpha)}(x,y), \quad \alpha=0,1,2,\cdots\right\}</math>. The sequence converges to a minimization problem with constraints when the initial control coefficients <math display="inline">C_{ij}^{(0)}=0</math>.<ref name=":5" />

Assume that the implicit surface generated after the <math display="inline">\alpha</math>th iteration is
<math display="block">
f^{(\alpha)}(x,y,z)=\sum_{i=1}^{N_u}\sum_{j=1}^{N_v}\sum_{k=1}^{N_w}C_{ijk}^{(\alpha)}B_i(x)B_j(y)B_k(z),
</math>
the iteration format is similar to that of the curve case.<ref name=":5" /><ref>{{Cite journal |last1=Liu |first1=Shengjun |last2=Liu |first2=Tao |last3=Hu |first3=Ling |last4=Shang |first4=Yuanyuan |last5=Liu |first5=Xinru |date=2021-09-01 |title=Variational progressive-iterative approximation for RBF-based surface reconstruction |url=https://doi.org/10.1007/s00371-021-02213-3 |journal=The Visual Computer |language=en |volume=37 |issue=9 |pages=2485–2497 |doi=10.1007/s00371-021-02213-3 |issn=1432-2315}}</ref>

=== Fairing-PIA ===

To develop fairing-PIA, we first define the functionals as follows:<ref name=":7" />
<math display="block">
\mathcal{F}_{r,j}(f) = \int_{t_1}^{t_m}B_{r,j}(t)fdt,\quad j=1,2,\cdots,n,\quad r=1,2,3,
</math>
where <math display="inline">B_{r,j}(t)</math> represents the <math display="inline">r</math>th derivative of the basis function <math display="inline">B_j(t)</math>,<ref name=":13" /> (e.g. [[B-spline|B-spline basis function]]).

Let the curve after the <math display="inline">k</math>th iteration be
<math display="block">
\mathbf{P}^{[k]}(t)=\sum_{j=1}^nB_j(t)\mathbf{P}_j^{[k]},\quad t\in[t_1,t_m].
</math>
To construct the new curve <math display="inline">\mathbf{P}^{[k+1]}(t)</math>, we first calculate the <math display="inline">(k + 1)</math>st difference vectors for data points,<ref name=":7" />
<math display="block">
\mathbf{d}_i^{[k]} = \mathbf{Q}_i - \mathbf{P}^{[k]}(t_i),\quad i=1,2,\cdots,m.
</math>
Then, the fitting difference vectors and the fairing vectors for control points are calculated by<ref name=":7" />
<math display="block">
\begin{align} 
\boldsymbol{\delta}_j^{[k]} &= \sum_{h\in I_j}B_j(t_h)\mathbf{d}_h^{[k]},\quad j=1,2,\cdots,n \\
\boldsymbol{\eta}_{j}^{[k]} &= \sum_{l=1}^n \mathcal{F}_{r,l}\left(B_{r,j}(t)\right)\mathbf{P}_l^{[k]},\quad j=1,2,\cdots,n \\ 
\end{align}
</math>
Finally, the control points of the <math>(k+1)</math>st curve are produced by<ref name=":7" />
<math display="block">
\mathbf{P}_j^{[k+1]} = \mathbf{P}_j^{[k]} + \mu_j
\left[
\left(1-\omega_j\right)\boldsymbol{\delta}_j^{[k]} - \omega_j\boldsymbol{\eta}_{j}^{[k]}
\right],\quad j=1,2,\cdots,n,
</math>
where <math>\mu_j</math> is a normalization weight, and <math>\omega_j</math> is a smoothing weight corresponding to the <math>j</math>th control point. The smoothing weights can be employed to adjust the smoothness individually, thus bringing great flexibility for smoothness.<ref name=":7" /> The larger the smoothing weight is, the smoother the generated curve is. The new curve is obtained as follows
<math display="block">
\mathbf{P}^{[k+1]}(t)=\sum_{j=1}^nB_j(t)\mathbf{P}_j^{[k+1]},\quad t\in[t_1,t_m].
</math>
In this way, we obtain a sequence of curves <math display="inline">\left\{\mathbf{P}^{[k]}(t),\;k=1,2,3,\cdots\right\}</math>. The sequence converges to the solution of the conventional fairing method based on energy minimization when all smoothing weights are equal (<math display="inline">\omega_j=\omega</math>).<ref name=":7" /> Similarly, the fairing-PIA can be extended to the surface case.

=== IG-LSPIA ===

Isogeometric least-squares progressive-iterative approximation (IG-LSPIA).<ref name=":8">{{Cite journal |last1=Jiang |first1=Yini |last2=Lin |first2=Hongwei |date=2023-02-10 |title=IG-LSPIA: Least Squares Progressive Iterative Approximation for Isogeometric Collocation Method |journal=Mathematics |volume=11 |issue=4 |pages=898 |doi=10.3390/math11040898 |issn=2227-7390 |doi-access=free }}</ref> Given a [[boundary value problem]]<ref name=":14" />
<math display="block">
\left\{
\begin{aligned}
\mathcal{L}u=f,&\quad \text{in}\;\Omega,\\
\mathcal{G}u=g,&\quad \text{on}\;\partial\Omega,
\end{aligned}
\right.
</math>
where <math display="inline">u:\Omega\to\mathbb{R}</math> is the unknown solution, <math display="inline">\mathcal{L}</math> is the [[differential operator]], <math display="inline">\mathcal{G}</math> is the boundary operator, and <math display="inline">f</math> and <math display="inline">g</math> are the continuous functions. In the [[Isogeometric analysis|isogeometric analysis method]], NURBS basis functions<ref name=":13" /> are used as shape functions to solve the numerical solution of this boundary value problem.<ref name=":14" /> The same basis functions are applied to represent the numerical solution <math display="inline">u_h</math> and the geometric mapping <math display="inline">G</math>:
<math display="block">
\begin{aligned}
u_h\left(\hat{\tau}\right) &= \sum_{j=1}^nR_{j}(\hat\tau )u_j,\\
G({\hat \tau })            &= \sum_{j=1}^nR_{j}(\hat\tau )P_j,
\end{aligned}
</math>
where <math display="inline">R_j(\hat{\tau})</math> denotes the NURBS basis function, <math display="inline">u_j</math> is the control coefficient. After substituting the collocation points<ref name=":15">{{Cite journal |last1=Lin |first1=Hongwei |last2=Hu |first2=Qianqian |last3=Xiong |first3=Yunyang |date=2013-12-01 |title=Consistency and convergence properties of the isogeometric collocation method |url=https://www.sciencedirect.com/science/article/pii/S0045782513002521 |journal=Computer Methods in Applied Mechanics and Engineering |volume=267 |pages=471–486 |doi=10.1016/j.cma.2013.09.025 |bibcode=2013CMAME.267..471L |issn=0045-7825}}</ref> <math display="inline">\hat\tau_{i} ,i = 1,2,...,{m}</math> into the strong form of [[Partial differential equation|PDE]], we obtain a discretized problem<ref name=":15" />
<math display="block">
\left\{
\begin{aligned}
\mathcal{L}u_{h}(\hat\tau_{i})=f(G(\hat\tau_{i})),&\quad i\in\mathcal{I_L},\\
\mathcal{G}u_{h}(\hat\tau_{j})=g(G(\hat\tau_{j})),&\quad j\in\mathcal{I_G},
\end{aligned}
\right.
</math>
where <math display="inline">\mathcal{I_L}</math> and <math display="inline">\mathcal{I_G}</math> denote the subscripts of internal and boundary collocation points, respectively.

Arranging the control coefficients <math display="inline">u_j</math> of the numerical solution <math display="inline">u_h(\hat\tau)</math> into an <math display="inline">1</math>-dimensional column vector <math display="inline">\mathbf{U}=[u_1,u_2,...,u_n]^T</math>, the discretized problem can be reformulated in matrix form as
<math display="block">
\mathbf{AU}=\mathbf{b}
</math>
where <math display="inline">\mathbf{A}</math> is the collocation matrix and <math display="inline">\mathbf{b}</math> is the load vector.

Assume that the discretized load values are data points <math display="inline">\left\{b_i\right\}_{i=1}^m</math> to be fitted. Given the initial guess of the control coefficients <math display="inline">\left\{u_j^{(0)}\right\}_{j=1}^n, n<m </math>, we obtain an initial blending function<ref name=":8" />
<math display="block">
U^{(0)}(\hat\tau) = \sum_{j=1}^nA_j(\hat\tau)u_j^{(0)},\quad\hat\tau\in[\hat\tau_1,\hat\tau_m],
</math>
where <math display="inline">A_j(\hat\tau)</math>, <math display="inline">j=1,2,\cdots,n</math>, represents the combination of different order derivatives of the NURBS basis functions determined using the operators <math display="inline">\mathcal{L}</math> and <math display="inline">\mathcal{G}</math>
<math display="block">
A_j(\hat\tau) = \left\{
\begin{aligned}
\mathcal{L}R_j(\hat\tau), &\quad \hat{\tau}\ \text{in}\  \Omega_p^{in},\\
\mathcal{G}R_j(\hat\tau), &\quad \hat{\tau}\ \text{in}\ \Omega_p^{bd}, \quad j=1,2,\cdots,n,
\end{aligned}
\right.
</math>
where <math display="inline">\Omega_p^{in}</math> and <math display="inline">\Omega_p^{bd}</math> indicate the interior and boundary of the parameter domain, respectively. Each <math display="inline">A_j(\hat\tau)</math> corresponds to the <math display="inline">j</math>th control coefficient. Assume that <math display="inline">J_{in}</math> and <math display="inline">J_{bd}</math> 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., 
<math display="block">
u_{j}^{(k)}=u_{j}^{*},\quad j\in J_{bd},\quad k=0,1,2,\cdots.
</math>
The <math display="inline">k</math>th blending function, generated after the <math display="inline">k</math>th iteration of IG-LSPIA,<ref name=":8" /> is assumed to be as follows:
<math display="block">
U^{(k)}(\hat\tau) = \sum_{j=1}^nA_j(\hat\tau)u_j^{(k)},\quad\hat\tau\in[\hat\tau_1,\hat\tau_m].
</math>
Then, the difference vectors for collocation points (DCP) in the <math display="inline">(k + 1)</math>st iteration are obtained using
<math display="block">
\begin{align}
\boldsymbol{\delta}_i^{(k)}
&= b_i-\sum_{j=1}^{n}A_j(\hat\tau_i)u_j^{(k)}\\
&= b_i-\sum_{j\in J_{bd}}A_j(\hat\tau_i)u_j^{(k)}
-\sum_{j\in J_{in}}A_j(\hat\tau_i)u_j^{(k)}
,\quad i=1,2,...,m.
\end{align}
</math>
Moreover, group all load values whose parameters fall in the local support of the <math display="inline">j</math>th derivatives function, i.e., <math display="inline">A_j(\hat\tau_i)\ne 0</math>, into the <math display="inline">j</math>th group corresponding to the <math display="inline">j</math>th control coefficient, and denote the index set of the <math display="inline">j</math>th group of load values as <math display="inline">I_j</math>. Lastly, the differences for control coefficients (DCC) can be constructed as follows:<ref name=":8" />
<math display="block">
d_j^{(k)}=\mu\sum_{h\in I_j}A_j(\hat\tau_h)\boldsymbol{\delta}_h^{(k)},\quad j=1,2,...,n,
</math>
where <math display="inline">\mu</math> is a normalization weight to guarantee the convergence of the algorithm.

Thus, the new control coefficients are updated via the following formula, 
<math display="block">
u_j^{(k+1)}=u_j^{(k)}+d_j^{(k)},\quad j=1,2,...,n,
</math>
Consequently, the <math display="inline">(k + 1)</math>st blending function is generated as follows:
<math display="block">
U^{(k+1)}(\hat\tau) = \sum_{j=1}^nA_j(\hat\tau)u_j^{(k+1)}.
</math>
The above iteration process is performed until the desired fitting precision is reached and a sequence of blending functions is obtained
<math display="block">
\left \{ U^{(k)}(\hat\tau),k=0,1,\dots  \right \}.
</math>
The IG-LSPIA converges to the solution of a constrained least-squares collocation problem.<ref name=":8" />

== Proof of convergence ==

=== Non-singular case ===
Let {{Mvar|n}} be the number of control points and {{Mvar|m}} be the number of data points.

If <math display="inline">n=m</math>, the PIA iterative format in matrix form is
<ref name=":0" /><math display="block">
\begin{align}
\mathbf{P^{(\alpha+1)}} &=\mathbf{P^{(\alpha)}}+\mathbf{\Delta}^{(\alpha)} \\
&=\mathbf{P}^{(\alpha)}+\mathbf{Q}-\mathbf{B}\mathbf{P}^{(\alpha)} \\
&=\left(\mathbf{I}-\mathbf{B}\right)\mathbf{P}^{(\alpha)}+\mathbf{Q}
\end{align}
</math>
where
<math display="block">
\begin{align}
\mathbf{Q} &= \left[\mathbf{Q}_1,\mathbf{Q}_2,\cdots,\mathbf{Q}_m\right]^T \\
\mathbf{P^{(\alpha)}} &= \left[\mathbf{P}_1^{(\alpha)},\mathbf{P}_2^{(\alpha)},\cdots,\mathbf{P}_n^{(\alpha)}\right]^T \\
\mathbf{\Delta}^{(\alpha)} &= \left[\mathbf{\Delta}_1^{(\alpha)},\mathbf{\Delta}^{(\alpha)}_2,\cdots,\mathbf{\Delta}^{(\alpha)}_n\right]^T \\
\mathbf{B} &= \begin{bmatrix}
B_1(t_1) & B_2(t_1) &\cdots &B_n(t_1)\\
B_1(t_2) & B_2(t_2) &\cdots &B_n(t_2)\\
\vdots   & \vdots   &\ddots & \vdots \\
B_1(t_m) & B_2(t_m) &\cdots &B_n(t_m)\\
\end{bmatrix}.
\end{align}
</math>
The convergence of the PIA is related to the properties of the collocation matrix. If the [[spectral radius]] of the iteration matrix <math>\mathbf{I}-\mathbf{B}</math> is less than <math>1</math>, 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).<ref name=":9" />

If <math display="inline">n<m</math>, the LSPIA in matrix form is<ref name=":3" /><ref name=":4" />
<math display="block">
\begin{align}
\mathbf{P^{(\alpha+1)}}&=\mathbf{P^{(\alpha)}}+\mu\mathbf{B}^T\mathbf{\Delta}^{(\alpha)} \\
&=\mathbf{P}^{(\alpha)}+\mu\mathbf{B}^T\left(\mathbf{Q}-\mathbf{B}\mathbf{P}^{(\alpha)}\right) \\
&=\left(\mathbf{I}-\mu\mathbf{B}^T\mathbf{B}\right)\mathbf{P}^{(\alpha)}+\mu\mathbf{B}^T\mathbf{Q}.
\end{align}
</math>
When the matrix <math display="inline">\mathbf{B}^T\mathbf{B}</math> is [[Invertible matrix|nonsingular]], the following results can be obtained:<ref>{{Cite book |last1=Horn |first1=Roger A. |url=http://dx.doi.org/10.1017/cbo9781139020411 |title=Matrix Analysis |last2=Johnson |first2=Charles R. |date=2012-10-22 |publisher=Cambridge University Press |doi=10.1017/cbo9781139020411 |isbn=978-0-521-83940-2}}</ref>

{{Math theorem |If <math display="inline">0<\mu<\frac{2}{\lambda_0}</math> , where <math display="inline">\lambda_0</math> is the largest [[Eigenvalues and eigenvectors|eigenvalue]] of the matrix <math display="inline">\mathbf{B}^T\mathbf{B}</math>, then the eigenvalues of <math display="inline">\mu\mathbf{B}^T\mathbf{B}</math> are real numbers and satisfy <math display="inline">0<\lambda(\mu\mathbf{B}^T\mathbf{B})<2</math>. |name=Lemma}}

'''Proof''' Since <math display="inline">\mathbf{B}^T\mathbf{B}</math> is nonsingular, and <math display="inline">\mu>0</math>, then <math display="inline">\lambda(\mu\mathbf{B}^T\mathbf{B})>0</math>. Moreover,
<math display="block">
\lambda(\mu\mathbf{B}^T\mathbf{B}) =\mu\lambda(\mathbf{B}^T\mathbf{B})<2\frac{\lambda(\mathbf{B}^T\mathbf{B})}{\lambda_0}<2.
</math>
In summary, <math display="inline">0<\lambda(\mu\mathbf{B}^T\mathbf{B})<2</math>.

{{Math theorem |If <math display="inline">0<\mu<\frac{2}{\lambda_0}</math> , LSPIA is convergent, and converges to the least-squares fitting result to the given data points.<ref name=":3" /><ref name=":4" /> |name=Theorem}}

'''Proof''' From the matrix form of iterative format, we obtain the following:
<math display="block">
\begin{align}
\mathbf{P^{(\alpha+1)}}&=\left(\mathbf{I}-\mu\mathbf{B}^T\mathbf{B}\right)\mathbf{P}^{(\alpha)}+\mu\mathbf{B}^T\mathbf{Q},\\
&=\left(\mathbf{I}-\mu\mathbf{B}^T\mathbf{B}\right)\left[\left(\mathbf{I}-\mu\mathbf{B}^T\mathbf{B}\right)\mathbf{P}^{(\alpha-1)}+\mu\mathbf{B}^T\mathbf{Q}\right]+\mu\mathbf{B}^T\mathbf{Q},\\
&=\left(\mathbf{I}-\mu\mathbf{B}^T\mathbf{B}\right)^2\mathbf{P}^{(\alpha-1)}+\sum_{i=0}^1\left( \mathbf{I}-\mu\mathbf{B}^T\mathbf{B}\right)\mu\mathbf{B}^T\mathbf{Q},\\
&=\cdots\\
&=\left(\mathbf{I}-\mu\mathbf{B}^T\mathbf{B}\right)^{\alpha+1}\mathbf{P}^{(0)}+\sum_{i=0}^{\alpha}\left( \mathbf{I}-\mu\mathbf{B}^T\mathbf{B}\right)^{\alpha}\mu\mathbf{B}^T\mathbf{Q}.\\
\end{align}
</math>
According to above Lemma, the spectral radius of the matrix <math display="inline">\mu\mathbf{B}^T\mathbf{B}</math> satisfies
<math display="block">
0<\rho\left({\mu\mathbf{B}^T\mathbf{B}}\right)<2
</math>
and thus the spectral radius of the iteration matrix satisfies
<math display="block">
0<\rho\left({\mathbf{I}-\mu\mathbf{B}^T\mathbf{B}}\right)<1.
</math>
When <math display="inline">\alpha\rightarrow \infty</math>
<math display="block">
\left(\mathbf{I}-\mu\mathbf{B}^T\mathbf{B}\right)^{\infty}=0,\ \sum_{i=0}^{\infty}\left( \mathbf{I}-\mu\mathbf{B}^T\mathbf{B}\right)^{\alpha}=\frac{1}{\mu}\left(\mathbf{B}^T\mathbf{B}\right)^{-1}.
</math>
As a result,
<math display="block">  \mathbf{P}^{(\infty)}=\left(\mathbf{B}^T\mathbf{B}\right)^{-1}\mathbf{B}^T\mathbf{Q},
</math>
i.e., <math display="inline">\mathbf{B}^T\mathbf{B}\mathbf{P}^{(\infty)}=\mathbf{B}^T\mathbf{Q}</math>, 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.<ref name=":12">{{Cite journal |last1=Lin |first1=Hongwei |last2=Cao |first2=Qi |last3=Zhang |first3=Xiaoting |title=The Convergence of Least-Squares Progressive Iterative Approximation for Singular Least-Squares Fitting System |journal=Journal of Systems Science and Complexity |date=2018 |volume=31 |issue=6 |pages=1618–1632 |doi=10.1007/s11424-018-7443-y |s2cid=255157830 |issn=1009-6124}}</ref>

== 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.<ref>{{Cite journal |last1=Liu |first1=Chengzhi |last2=Han |first2=Xuli |last3=Li |first3=Juncheng |date=2020 |title=Preconditioned progressive iterative approximation for triangular Bézier patches and its application |journal=Journal of Computational and Applied Mathematics |volume=366 |pages=112389 |doi=10.1016/j.cam.2019.112389 |issn=0377-0427 |s2cid=202942809}}</ref>
* '''Iteration matrix inverse approximation''': Sajavičius improved the LSPIA based on the matrix approximate inverse method. In each iteration step, the approximate inverse of the coefficient matrix of the least-squares fitting problem is first computed and then used as the weight to adjust the control points.<ref>{{Cite journal |last=Sajavičius |first=Svajūnas |date=2023 |title=Hyperpower least squares progressive iterative approximation |journal=Journal of Computational and Applied Mathematics |volume=422 |pages=114888 |doi=10.1016/j.cam.2022.114888 |issn=0377-0427 |s2cid=252965212}}</ref>
* '''Optimal weight''': Lu initially presented a weighted progressive-iterative approximation (WPIA) that introduces the optimal weight of difference vectors for control points to accelerate the convergence.<ref>{{Cite journal |last=Lu |first=Lizheng |title=Weighted progressive iteration approximation and convergence analysis |journal=Computer Aided Geometric Design |date=2010 |volume=27 |issue=2 |pages=129–137 |doi=10.1016/j.cagd.2009.11.001 |issn=0167-8396}}</ref> Moreover, Zhang et al. proposed a weighted local PIA format for tensor Bézier surfaces.<ref>{{Cite journal |last=Zhang |first=Li |date=2014-05-01 |title=Weighted Local Progressive Iterative Approximation for Tensor-product B\'{e}zier Surfaces |journal=Journal of Information and Computational Science |volume=11 |issue=7 |pages=2117–2124 |doi=10.12733/jics20103359 |issn=1548-7741}}</ref> Li et al. assigned initial weights to each data point, and the weights of the interpolated points are determined adaptively during the iterative process.<ref>{{Cite journal |last1=Li |first1=Shasha |last2=Xu |first2=Huixia |last3=Deng |first3=Chongyang |year=2019 |title=Data-weighted least square progressive and iterative approximation and related B-Spline curve fitting |journal=Journal of Computer-Aided Design & Computer Graphics |volume=31 |issue=9 |pages=1574–1580}}</ref>
* '''Acceleration with memory''': In 2020, Huang et al. proposed a PIA method with memory for least square fitting (MLSPIA), which has a similar format to the momentum method. MLSPIA generates a series of fitting curves with three weights by iteratively adjusting the control points. With appropriate parameter selection, these curves converge to the least squares fit results for a given data point and are more efficient than LSPIA.<ref>{{Cite journal |last1=Huang |first1=Zheng-Da |last2=Wang |first2=Hui-Di |date=2020 |title=On a progressive and iterative approximation method with memory for least square fitting |journal=Computer Aided Geometric Design |volume=82 |pages=101931 |arxiv=1908.06417 |doi=10.1016/j.cagd.2020.101931 |issn=0167-8396 |s2cid=201070122}}</ref>
* '''Stochastic descent strategy''': Rios and Jüttle explored the relationship between LSPIA and gradient descent method and proposed a stochastic LSPIA algorithm with parameter correction.<ref>{{Cite journal |last1=Rios |first1=Dany |last2=Jüttler |first2=Bert |date=2022 |title=LSPIA, (stochastic) gradient descent, and parameter correction |journal=Journal of Computational and Applied Mathematics |volume=406 |pages=113921 |doi=10.1016/j.cam.2021.113921 |issn=0377-0427 |s2cid=244018717}}</ref>

== Applications ==

Since PIA has obvious geometric meaning, constraints can be easily integrated in the iterations. Currently, PIA has been widely applied in many fields, such as data fitting, reverse engineering, geometric design, mesh generation, data compression, fairing curve and surface generation, and isogeometric analysis.

=== Data fitting ===
* Adaptive data fitting: The control points are divided into ''active control points'' and ''fixed control points''. In each round of iteration, if the fitting error of a data point reaches a given precision, its corresponding control point is fixed and not updated. This iterative process is repeated until all control points are fixed. The algorithm performs well on large-scale data fitting by adaptively reducing the number of active control points.<ref>{{Cite journal |last=Lin |first=Hongwei |date=2012 |title=Adaptive data fitting by the progressive-iterative approximation |journal=Computer Aided Geometric Design |volume=29 |issue=7 |pages=463–473 |doi=10.1016/j.cagd.2012.03.005 |issn=0167-8396}}</ref>
* Large-scale data fitting: By combining T-spline with PIA, an incremental fitting algorithm suitable for fitting large-scale data sets is proposed. During the incremental iteration, each new round of iterations reuses information from the last round of iterations to save computation. While the convergence speed of the traditional point-by-point iterative algorithm decreases as the number of control points increases, in PIA the computation of each iteration step is unrelated to the number of control points; this gives PIA a powerful capability for data fitting.<ref name=":3" />
* Local fitting: Based on the local property of PIA, a series of local PIA formats have been proposed.<ref name=":6" /><ref>{{Cite journal |last1=Zhao |first1=Yu |last2=Lin |first2=Hongwei |last3=Bao |first3=Hujun |year=2012 |title=Local progressive interpolation for subdivision surface fitting |journal=Computer Research and Development |volume=49 |issue=8 |pages=1699–1707}}</ref>

=== Implicit reconstruction ===
For implicit curve and surface reconstruction, PIA avoids the additional zero level set and regularization term, which greatly improves the speed of the reconstruction algorithm.<ref name=":5" />

=== Offset curve approximation ===
Firstly, the data points are sampled on the original curve. Then, the initial polynomial approximation curve or rational approximation curve of the offset curve is generated from these sampled points. Finally, the offset curve is approximated iteratively using the PIA method.<ref>{{Cite journal |last1=Zhang |first1=Li |last2=Wang |first2=Huan |last3=Li |first3=Yuanyuan |last4=Tan |first4=Jieqing |year=2014 |title=A progressive iterative approximation method in offset approximation |journal=Journal of Computer Aided Design and Computer Graphics |volume=26 |issue=10 |pages=1646–1653}}</ref>

=== Mesh generation ===
Given a triangular mesh model as input, the algorithm first constructs the initial hexahedral mesh, then extracts the quadrilateral mesh of the surface as the initial boundary mesh. During the iterations, the movement of each mesh vertex is constrained to ensure the validity of the mesh. Finally, the hexahedral model is fitted to the given input model. The algorithm can guarantee the validity of the generated hexahedral mesh, i.e., the Jacobi value at each mesh vertex is greater than zero.<ref>{{Cite journal |last1=Lin |first1=Hongwei |last2=Jin |first2=Sinan |last3=Liao |first3=Hongwei |last4=Jian |first4=Qun |year=2015 |title=Quality guaranteed all-hex mesh generation by a constrained volume iterative fitting algorithm |journal=Computer-Aided Design |volume=67-68 |pages=107–117 |doi=10.1016/j.cad.2015.05.004 |issn=0010-4485}}</ref>

=== Data compression ===
First, the image data are converted into a one-dimensional sequence by Hilbert scan. Then, these data points are fitted by LSPIA to generate a Hilbert curve. Finally, the Hilbert curve is sampled, and the compressed image can be reconstructed. This method can well preserve the neighborhood information of pixels.<ref>{{Cite journal |last1=Hu |first1=Lijuan |last2=Yi |first2=Yeqing |last3=Liu |first3=Chengzhi |last4=Li |first4=Juncheng |year=2020 |title=Iterative method for image compression by using LSPIA |journal=IAENG International Journal of Computer Science |volume=47 |issue=4 |pages=1–7}}</ref>

=== Fairing curve and surface generation ===
Given a data point set, we first define the fairing functional, and calculate the fitting difference vector and the fairing vector of the control point; then, adjust the control points with fairing weights. According to the above steps, the fairing curve and surface can be generated iteratively. Due to the sufficient fairing parameters, the method can achieve global or local fairing. It is also flexible to adjust knot vectors, fairing weights, or data parameterization after each round of iteration. The traditional energy-minimization method is a special case of this method, i.e., when the smooth weights are all the same.<ref name=":7" />

=== Isogeometric analysis ===
The discretized load values are regarded as the set of data points, and the combination of the basis functions and their derivative functions is used as the blending function for fitting. The method automatically adjusts the degrees of freedom of the numerical solution of the partial differential equation according to the fitting result of the blending function to the load values. In addition, the average iteration time per step is only related to the number of data points (i.e., collocation points) and unrelated to the number of control coefficients.<ref name=":8" />

== References ==
*{{Creative Commons text attribution notice|cc=by4|url=https://www.mdpi.com/2227-7390/11/4/898 |author(s)=zju_cagd}}
{{reflist}}

{{Improve categories|date=August 2024}}
[[Category:Computer-aided design]]
[[Category:Computer graphics]]
[[Category:Computational geometry]]
[[Category:Geometric algorithms]]
[[Category:Curve fitting]]