Abstract
We propose a modified bivariate Shepard operator of least squares thinplate spline type based on an iterative method, introduced by A.V. Masjukov and V.V. Masjukov. This iterative method does not require an artificial setup of the parameters and it is based on successive scaling. For the numerical examples, using an idea of J. R. McMahon, we construct some representative sets of knot points for the initial sets of the interpolation nodes.
Authors
Keywords
Paper coordinates
A. Malina, Iterative Shepard operator of least squares thinplate spline type, Dolomites Res. Notes Approx., 16 (2023) no. 3, pp. 5762. http://doi.org/10.14658/PUPJDRNA202338
About this paper
Journal
Dolomites Research Notes on Approximation
Publisher Name
Padova University Press
DOI
10.14658/PUPJDRNA202338
Online ISSN
20356803
Google Scholar Profile
Paper (preprint) in HTML form
Iterative Shepard operator of least squares thinplate spline type
Abstract
We propose a modified bivariate Shepard operator of least squares thinplate spline type based on an iterative method, introduced by A.V. Masjukov and V.V. Masjukov. This iterative method does not require an artificial setup of the parameters and it is based on successive scaling. For the numerical examples, using an idea of J. R. McMahon, we construct some representative sets of knot points for the initial sets of the interpolation nodes.
Keywords: Shepard operator, least squares approximation, thinplate spline, knot points, iterative multiscale method.
2010 Mathematics Subject Classification: 41A05, 41A25.
1 Preliminaries
The Shepard method, introduced by D. Shepard in [16], is one of the best suited methods for interpolating scattered data. It has many computational advantages and many of its drawbacks have been overcome during the time by combining it with other interpolation operators or by modifying the basis functions. One way to improve the method is described by A. V. Masjukov and V. V. Masjukov in [11] and consists of constructing an iterative multiscale method. This idea has been used by other authors in order to obtain better approximation results (see, e.g., [5]).
The modified Shepard operator introduced by Franke and Nielson in [10] and later on improved in [9], [14], [15] requires some artificial parameters such as a number of closest nodes or a radius of influence. The iterative method [11] is free of these setup parameters and performs at each iteration a reduction of the current interpolation result’s residue. The accuracy of this method is comparable to the modified Shepard procedure as shown in [11], although for both methods there are cases when one is preferred to the other.
Consider $f:X\subset {\mathbb{R}}^{2}\to \mathbb{R}$ and the distinct nodes $({x}_{i},{y}_{i})\in X,i=1,\mathrm{\dots},N$. Denoting by ${d}_{i}(x,y)$ the distance between a given point $(x,y)\in X$ and the node $({x}_{i},{y}_{i}),i=1,\mathrm{\dots},N$, the Shepard operator in the bivariate case is defined as
$$\left({S}^{\mu}f\right)(x,y)=\sum _{i=1}^{N}{A}_{i,\mu}(x,y)f({x}_{i},{y}_{i}),$$ 
where
$${A}_{i,\mu}(x,y)=\frac{\prod _{\begin{array}{c}j=1\\ j\ne i\end{array}}^{N}{d}_{j}^{\mu}(x,y)}{\sum _{k=1}^{N}\prod _{\begin{array}{c}j=1\\ j\ne k\end{array}}^{N}{d}_{j}^{\mu}(x,y)},$$  (1) 
with the parameter $\mu >0$, whose choice could change the behaviour of the operator.
It is wellknown that method’s accuracy can be improved by considering the modified Shepard operator, defined as:
$$\left({S}^{w}f\right)(x,y)=\frac{\sum _{i=1}^{N}{W}_{i}(x,y)f({x}_{i},{y}_{i})}{\sum _{i=1}^{N}{W}_{i}(x,y)},$$ 
with
$${W}_{i}(x,y)={\left[\frac{{({R}^{w}{d}_{i})}_{+}}{{R}^{w}{d}_{i}}\right]}^{2},$$  (2) 
where
$${u}_{+}=\{\begin{array}{cc}u,\hfill & \text{if}u0,\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}$$ 
and ${R}^{w}$ is a radius of influence about the node $({x}_{i},{y}_{i})$ which depends on $i$. It is computed as the distance between the node $i$ and the $n$th nearest node to it for $n>{N}_{w}$, where ${N}_{w}$ is fixed. In addition, $n$ should be as small as possible within the constraint that the $n$th nearest node is significantly more distant than the $(n1)$st nearest node (see, e.g., [15]).
Starting with Shepard himself, many other authors have proposed to apply the operator not on $f({x}_{i},{y}_{i})$ directly, but on some interpolation operator $P(f)({x}_{i},{y}_{i})$ in order to obtain better results. So, the combined Shepard operator is given by
$$\left({S}_{P}^{\mu}f\right)(x,y)=\sum _{i=1}^{N}{A}_{i,\mu}(x,y)P(f)({x}_{i},{y}_{i}),$$ 
with ${A}_{i,\mu}$ defined in (1) and $\mu >0$.
Based on the least squares thinplate spline type function, the classical and the modified Shepard operators, it is introduced in [6] the Shepard operator of least squares thinplate spline type. An improvement of the method is obtained by replacing the initial set of nodes with a representative smaller set of knot points, using an idea described by McMahon and Franke in [12] and [13].
2 The Iterative Shepard operator of least squares thinplate spline type
Consider $f:X\subset {\mathbb{R}}^{2}\to \mathbb{R}$ and the distinct nodes $({x}_{i},{y}_{i})\in X,i=1,\mathrm{\dots},N$. The classical Shepard operator combined with a least squares thinplate spline is defined in [6] as
$$({S}_{1}f)(x,y)=\sum _{i=1}^{N}{A}_{i,\mu}(x,y){F}_{i}(x,y),$$  (3) 
where ${A}_{i,\mu},$ $i=1,\mathrm{\dots},N,$ are defined by (1), for a given parameter $\mu >0$.
The modified Shepard operator of least squares thinplate spline type has the form
$$({S}_{2}f)(x,y)=\frac{\sum _{i=1}^{N}{W}_{i}(x,y){F}_{i}(x,y)}{\sum _{i=1}^{N}{W}_{i}(x,y)},$$  (4) 
where ${W}_{i},i=1,\mathrm{\dots},N,$ are defined by (2).
In both cases the least squares thinplate spline ${F}_{i},i=1,\mathrm{\dots},N$, is given by
$${F}_{i}(x,y)=\sum _{j=1}^{i}{C}_{j}{d}_{j}^{2}\mathrm{log}({d}_{j})+ax+by+c,$$  (5) 
with ${d}_{j}=\sqrt{{(x{x}_{j})}^{2}+{(y{y}_{j})}^{2}}.$
The coefficients ${C}_{j},$ $a,$ $b,$ $c$ are found in order to minimize the expression (see, e.g., [6], [12]):
$$E=\sum _{i=1}^{N}{[{F}_{i}({x}_{i},{y}_{i})f({x}_{i},{y}_{i})]}^{2},$$ 
by solving a system that has the form
$$\left(\begin{array}{cc}M& X\\ {X}^{T}& {O}_{3}\end{array}\right)\cdot \left(\begin{array}{c}\bm{C}\\ \mathbf{v}\end{array}\right)=\left(\begin{array}{c}\mathbf{f}\\ \mathrm{\U0001d7ce}\end{array}\right),$$ 
considering the following notations:

•
$M={\left({m}_{ij}\right)}_{i,j=\overline{1,N}}$, ${m}_{ij}=\{\begin{array}{cc}{d}_{ij}^{2}\mathrm{log}({d}_{ij}),& i\ne j\\ 0,& i=j\end{array}$, with ${d}_{ij}=\sqrt{{({x}_{i}{x}_{j})}^{2}+{({y}_{i}{y}_{j})}^{2}}$;

•
$\mathbf{X}=\left(\begin{array}{ccc}{x}_{1}& {y}_{1}& 1\\ \mathrm{\vdots}& \mathrm{\vdots}& \mathrm{\vdots}\\ {x}_{N}& {y}_{N}& 1\end{array}\right)$;

•
$C={({C}_{1},\mathrm{\dots},{C}_{N})}^{T}$, $\mathbf{v}={(a,b,c)}^{T}$, $\mathbf{f}={({f}_{1},\mathrm{\dots},{f}_{N})}^{T}$ and $\mathrm{\U0001d7ce}={(0,0,0)}^{T}$.
To improve method’s results, the initial set of $N$ nodes has been replaced by a smaller subset of $k$ representative knot points $({x}_{i},{y}_{i}),i=1,\mathrm{\dots},k$, following an idea from [12]. The algorithm is also described in [6]. In short, using a randomly generated set of knot points, the procedure consists in replacing the set of knot points by the centroids of each subset of points that are closest to the same knot, until the set of knots has not changed for two consecutive iterations. Using these knot points, the classical Shepard operator and the modified Shepard operator are defined in a similar way.
The iterative Shepard operator introduced in [11] is of the following form
$$u(x,y)=\sum _{k=0}^{K}\sum _{j=1}^{{N}^{\prime}}\left[{u}_{j}^{(k)}w\left((x{x}_{j},y{y}_{j})/{\tau}_{k}\right)/\sum _{p=1}^{{N}^{\prime}}w\left(({x}_{p}{x}_{j},{y}_{p}{y}_{j})/{\tau}_{k}\right)\right],$$  (6) 
with $(x,y)\in X,({x}_{i},{y}_{i})\in X,i=1,\mathrm{\dots},{N}^{\prime}$, are the interpolation nodes and $w$ is a continuously differentiable weight function, having the properties that
$$w(x,y)\ge 0,\forall (x,y)\in {\mathbb{R}}^{2},w(0,0)>0\text{and}w(x,y)=0\text{if}\Vert (x,y)\Vert 1.$$ 
${u}_{j}^{(k)}$ denotes the interpolation residuals at the $k$th step, with ${u}_{j}^{(0)}\equiv {u}_{j}$.
Keeping the ideas from [11] and based on the method described in the paper, we define the iterative Shepard operator of least squares thinplate spline type, denoted by ${u}_{L}$, as
$${u}_{L}(x,y)=\sum _{k=0}^{K}\sum _{j=1}^{{N}^{\prime}}\left[{u}_{{F}_{j}}^{(k)}w\left((x{x}_{j},y{y}_{j})/{\tau}_{k}\right)/\sum _{p=1}^{{N}^{\prime}}w\left(({x}_{p}{x}_{j},{y}_{p}{y}_{j})/{\tau}_{k}\right)\right],$$  (7) 
where $(x,y)\in X,({x}_{i},{y}_{i})\in X,i=1,\mathrm{\dots},{N}^{\prime}$, and the interpolation residuals at the $k$th step ${u}_{{F}_{j}}^{(k)}$ are given by
$${u}_{{F}_{j}}^{(0)}={F}_{j}({x}_{j},{y}_{j}),({x}_{j},{y}_{j})\in X,j=1,\mathrm{\dots},{N}^{\prime}$$ 
and
$${u}_{{F}_{j}}^{(k+1)}={u}_{{F}_{j}}^{(k)}\sum _{q=1}^{{N}^{\prime}}\left[{u}_{{F}_{q}}^{(k)}w\left(({x}_{j}{x}_{q},{y}_{j}{y}_{q})/{\tau}_{k}\right)/\sum _{p=1}^{{N}^{\prime}}w\left(({x}_{p}{x}_{q},{y}_{p}{y}_{q})/{\tau}_{k}\right)\right].$$ 
The functions ${F}_{i}$ are the least squares thinplate splines, given in (5). The parameter ${\tau}_{k}$ is chosen as in [11] and it decreases from a given value ${\tau}_{0}$, which can be, for example
$${\tau}_{0}>\underset{(x,y)\in X}{sup}\underset{1\le j\le {N}^{\prime}}{\mathrm{max}}\Vert (x{x}_{j},y{y}_{j})\Vert $$ 
to
$$ 
The sequence $\{{\tau}_{k}\}$ of scale factors is given by
$${\tau}_{k}={\gamma}^{k}{\tau}_{0},\gamma \in (0,1).$$ 
It was also shown in [11] that by applying the idea of multiscale analysis, the behaviour of the interpolant does not change very much for values of $\gamma $ from $0.6$ to $0.95$. Smaller values of $\gamma $ could also be chosen for decreased computational time, in the case of sparsed interpolation nodes.
The weight function $w$ is defined as
$$w(x,y)=w(x)w(y)$$ 
with
$$ 
We construct the iterative Shepard operator of least squares thinplate spline type on both sets of interpolation nodes that were mentioned before: firstly for ${N}^{\prime}=N$ nodes and then for ${N}^{\prime}=k$ knot points that are representative for the initial set of $N$ nodes, with $$.
3 Test results
We consider the following test functions (see, e.g., [9], [14], [15]):
$$\begin{array}{cc}\text{Gentle:}\hfill & {f}_{1}(x,y)=\mathrm{exp}[\frac{81}{16}({(x0.5)}^{2}+{(y0.5)}^{2})]/3,\hfill \\ \text{Saddle:}\hfill & {f}_{2}(x,y)=\frac{(1.25+\mathrm{cos}5.4y)}{6+6{(3x1)}^{2}}.\hfill \end{array}$$  (8) 
Table 1 contains the maximum errors for approximating the functions given in (8) by the corresponding two iterative Shepard operators given in (7), for N=100 random generated points in $[0,1]\times [0,1]$ and k=25 knot points. We consider ${\tau}_{0}=3$, $K=20$ and $\gamma =0.66,\mathrm{\hspace{0.33em}0.84}$. We also compute the maximum errors for the classical and the modified Shepard operators of least squares type given in (3) and (4), respectively. We consider both sets of nodes, the 100 initial interpolation nodes and the 25 representative knot points, for $\mu =3$ and ${N}_{w}=19$.
In Figures 1  2 we plot the graphs of ${f}_{1},{f}_{2}$ and the corresponding iterative Shepard operator ${u}_{L}$ when ${N}^{\prime}=N=100$ and ${N}^{\prime}=k=25$.
We observe that in all the cases, the choice of the 25 representative knot points produces better approximation results, this illustrating the benefits of the algorithm of choosing the representative set of knot points. Also, the smallest error for the approximation of the Saddle function is attained using the iterative Shepard operator of least squares type and for the approximation of the Gentle function is obtained for the modified Shepard operator of least squares type. Because the differences are not large, the two methods are comparable, depending on the set of data that is used, both of them having important advantages, as previously mentioned.
${f}_{1}$  ${f}_{2}$  

${u}_{L}$ , $\gamma =0.66$  $k=25$  0.1635  0.2166 
$N=100$  0.1693  0.2496  
${u}_{L}$, $\gamma =0.84$  $k=25$  0.1336  0.2176 
$N=100$  0.1603  0.2233  
${S}_{1}f$  $k=25$  0.1578  0.3783 
$N=100$  0.1644  0.4001  
${S}_{2}f$  $k=25$  0.1212  0.2834 
$N=100$  0.1246  0.2858 
4 Conclusions
In this paper we have proposed an iterative modification of the Shepard operator of least squares thinplate spline type, introduced in [6]. For the given sets of interpolation nodes, smaller subsets of representative knot points are constructed and the numerical results demonstrate the benefits of this algorithm. The advantages of this method are also confirmed by the numerical results.
References
 [1] T. Cătinaş. The combined ShepardAbelGoncharov univariate operator. Rev. Anal. Numér. Théor. Approx., 32:11–20, 2003.
 [2] T. Cătinaş. The combined ShepardLidstone bivariate operator. In: de Bruin, M.G. et al. (eds.): Trends and Applications in Constructive Approximation. International Series of Numerical Mathematics, Springer GroupBirkhäuser Verlag, 151:77–89, 2005.
 [3] T. Cătinaş. Bivariate interpolation by combined Shepard operators. Proceedings of ${17}^{th}$ IMACS World Congress, Scientific Computation, Applied Mathematics and Simulation, 7pp., 2005.
 [4] T. Cătinaş. The bivariate Shepard operator of Bernoulli type. Calcolo, 44(4):189–202, 2007.
 [5] T. Cătinaş. An iterative modification of ShepardBernoulli Operator. Results Math., 69:387–395, 2016.
 [6] T. Cătinaş, A. Malina. Shepard operator of least squares thinplate spline type. Stud. Univ. BabeşBolyai Math. 66(2):257–265, 2021.
 [7] Gh. Coman. Hermitetype Shepard operators. Rev. Anal. Numér. Théor. Approx., 26:33–38, 1997.
 [8] Gh. Coman. Shepard operators of Birkhoff type. Calcolo, 35:197–203, 1998.
 [9] R. Franke. Scattered data interpolation: tests of some methods. Math. Comp., 38:181–200, 1982.
 [10] R. Franke, G. Nielson. Smooth interpolation of large sets of scattered data. Int. J. Numer. Meths. Engrg. 15:1691–1704, 1980.
 [11] A.V. Masjukov, V.V. Masjukov. Multiscale modification of Shepard’s method for interpolation of multivariate scattered data. Mathematical Modelling and Analysis, Proceeding of the 10th International Conference MMA2005 & CMAM2, 467472.
 [12] J.R. McMahon. Knot selection for least squares approximation using thin plate splines. M.S. Thesis, Naval Postgraduate School, 1986.
 [13] J.R. McMahon, R. Franke. Knot selection for least squares thin plate splines. Technical Report, Naval Postgraduate School, Monterey, 1987.
 [14] R.J. Renka, A.K. Cline. A trianglebased ${C}^{1}$ interpolation method. Rocky Mountain J. Math., 14:223–237, 1984.
 [15] R.J. Renka. Multivariate interpolation of large sets of scattered data. ACM Trans. Math. Software, 14:139–148, 1988.
 [16] D. Shepard. A two dimensional interpolation function for irregularly spaced data. Proceedings of the 23rd ACM National Conference, New York: ACM Press, 517–523, 1968.
[1] T. Cătinaș. The combined ShepardAbelGoncharov univariate operator. Rev. Anal. Numér. Théor. Approx., 32:11–20, 2003.
[2] T. Cătinaș. The combined ShepardLidstone bivariate operator. In: de Bruin, M.G. et al. (eds.): Trends and Applications in Constructive
Approximation. International Series of Numerical Mathematics, Springer GroupBirkhäuser Verlag, 151:77–89, 2005.
[3] T. Cătinaș. Bivariate interpolation by combined Shepard operators. Proceedings of 17th IMACS World Congress, Scientific Computation,
Applied Mathematics and Simulation, 7pp., 2005.
[4] T. Cătinaș. The bivariate Shepard operator of Bernoulli type. Calcolo, 44(4):189–202, 2007.
[5] T. Cătinaș. An iterative modification of ShepardBernoulli Operator. Results Math., 69:387–395, 2016.
[6] T. Cătinaș., A. Malina. Shepard operator of least squares thinplate spline type. Stud. Univ. Babe¸sBolyai Math. 66(2):257–265, 2021.
[7] Gh. Coman. Hermitetype Shepard operators. Rev. Anal. Numér. Théor. Approx., 26:33–38, 1997.
[8] Gh. Coman. Shepard operators of Birkhoff type. Calcolo, 35:197–203, 1998.
[9] R. Franke. Scattered data interpolation: tests of some methods. Math. Comp., 38:181–200, 1982.
[10] R. Franke, G. Nielson. Smooth interpolation of large sets of scattered data. Int. J. Numer. Meths. Engrg. 15:1691–1704, 1980.
[11] A.V. Masjukov, V.V. Masjukov. Multiscale modification of Shepard’s method for interpolation of multivariate scattered data. Mathematical
Modelling and Analysis, Proceeding of the 10th International Conference MMA2005 & CMAM2, 467472.
[12] J.R. McMahon. Knot selection for least squares approximation using thin plate splines. M.S. Thesis, Naval Postgraduate School, 1986.
[13] J.R. McMahon, R. Franke. Knot selection for least squares thin plate splines. Technical Report, Naval Postgraduate School, Monterey, 1987.
[14] R.J. Renka, A.K. Cline. A trianglebased C1 interpolation method. Rocky Mountain J. Math., 14:223–237, 1984.
[15] R.J. Renka. Multivariate interpolation of large sets of scattered data. ACM Trans. Math. Software, 14:139–148, 1988.
[16] D. Shepard. A two dimensional interpolation function for irregularly spaced data. Proceedings of the 23rd ACM National Conference, New
York: ACM Press, 517–523, 1968.