## Two-dimensional interpolation using a cell-based searching procedure

## Roberto Cavoretto

Department of Mathematics 'G. Peano', University of Torino, Via Carlo Alberto 10, 10123 Torino, Italy roberto.cavoretto@unito.it

## Abstract

In this paper we present an efficient algorithm for bivariate interpolation, which is based on the use of the partition of unity method for constructing a global interpolant. It is obtained by combining local radial basis function interpolants with locally supported weight functions. In particular, this interpolation scheme is characterized by the construction of a suitable partition of the domain in cells so that the cell structure strictly depends on the dimension of its subdomains. This fact allows us to construct an efficient cell-based searching procedure, which provides a significant reduction of CPU times. Complexity analysis and numerical results show such improvements on the algorithm performances.

Keywords : Partition of unity, Local methods, Searching techniques, Fast algorithms, Scattered data

## 1. Introduction

Let { ( x i , f i ) , i = 1 2 , , . . . , n } be a finite set of discrete data, with x i ∈ Ω ⊆ R 2 , and f i ∈ R . The x i are called the nodes , while the (corresponding) f i are the data values . The latter are obtained by sampling some (unknown) function f : Ω → R at the nodes, i.e. f i = f ( x i ), i = 1 2 , , . . . , n .

Therefore, the scattered data interpolation problem consists in finding a continuous function R : Ω → R such that

$$\mathcal { R } ( \mathbf x _ { i } ) = f _ { i }, \quad \ i = 1, 2, \dots, n.$$

Here, we consider the problem of constructing an efficient algorithm for bivariate interpolation of (large) scattered data sets. It is based on the partition of unity method for constructing a global interpolant by blending radial basis functions (RBFs) as local approximants and using locally supported weight functions.

Now, starting from the results of previous researches (see [1, 3, 4, 5, 6]) where efficient searching procedures based on the partition of the domain in strips or spherical zones are considered, we extend the previous ideas replacing the strip-based partition structure with a cell-based one. The latter leads to the creation of a cell-based searching procedure, whose origin comes from the repeated use of a quicksort routine with respect to different directions, enabling

us to pass from not ordered to ordered data structures. In particular, this process is strictly related to the construction of a partition of the domain Ω in square cells, which consists in generating two orthogonal families of parallel strips, where the original data set is suitably split up in ordered data subsets.

Then, exploiting the ordered data structure and the domain partition, the cell algorithm is efficiently implemented and optimized by connecting the method itself with the effective cell-based searching procedure. More precisely, the considered technique is characterized by the construction of a double structure of crossed strips, which partitions the domain in square cells and strictly depends on the dimension of its subdomains, providing a significant improvement in the searching procedures of the nearest neighbour points compared to the searching techniques in [1, 3, 4]. The final result is an efficient algorithm for bivariate interpolation of scattered data. Finally, complexity analysis and numerical tests show the high efficiency of the proposed algorithm.

The paper is organized as follows. In Section 2 we recall some theoretical results: firstly, we introduce the radial basis functions referring to existence and uniqueness of RBF interpolants, then we give a general description of the partition of unity method, which makes use of local RBF approximants. In Section 3, we present in detail the cell algorithm, which is efficiently implemented and optimized by using a cell-based searching procedure. Then, in Section 4 complexity of this algorithm is analyzed as well. Finally, Section 5 shows numerical results concerning efficiency and accuracy of the cell-based partition algorithm.

## 2. Local interpolation scheme

## 2.1. Radial basis functions

A suitable approach to solving the scattered data interpolation problem is to make the assumption that the interpolating function R is expressed as a linear combination of radial basis functions φ : [0 , ∞ → ) R , i.e.,

$$\mathcal { R } ( \mathbf x ) = \sum _ { j = 1 } ^ { n } c _ { j } \phi ( | | \mathbf x - \mathbf x _ { j } | | _ { 2 } ), \quad \mathbf x \in \mathbb { R } ^ { 2 },$$

where ||·|| 2 is the Euclidean distance, and R satisfies the interpolation conditions (1).

Thus, solving the interpolation problem under this assumption leads to a system of linear equations of the form

$$A c = f,$$

where the entries of the interpolation matrix A are given by

$$a _ { i j } = \phi ( | | \mathbf x _ { i } - \mathbf x _ { j } | | _ { 2 } ), \quad \ i, j = 1, 2, \dots, n,$$

c = [ c 1 , c 2 , . . . , c n ] T , and f = [ f , f 1 2 , . . . , f n ] T . Moreover, we know that the interpolation problem is well-posed , that is a solution to a problem exists and

is unique, if and only if the matrix A is non-singular. A sufficient condition to have non-singularity is that the corresponding matrix is positive definite (see, e.g., [7]).

Now, we remind that a real-valued continuous even function φ is called positive definite on R 2 if

$$\sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { n } c _ { i } c _ { j } \phi \left ( { \mathbf x } _ { i } - { \mathbf x } _ { j } \right ) \geq 0$$

for any n pairwise distinct nodes x 1 , x 2 , . . . , x n ∈ R 2 , and c = [ c 1 , c 2 , . . . , c n ] T ∈ R n . The function φ is called strictly positive definite on R 2 if the quadratic form (4) is zero only for c ≡ 0 .

Therefore, if φ is strictly positive definite, the interpolant (2) is unique, since the corresponding interpolation matrix (3) is positive definite and hence non-singular.

Some of the most popular strictly positive definite RBFs are

$$\begin{smallmatrix} \phi _ { G } ( r ) & = & e ^ { - \alpha ^ { 2 } r ^ { 2 } }, \\ \phi _ { W 2 } ( r ) & = & ( 1 - c r ) _ { + } ^ { 4 } ( 4 c r + 1 ), \end{smallmatrix} \quad ( \text{Wendand} ^ { \prime } s \ C ^ { 2 } \text{ function} )$$

where α, c ∈ R + are the shape parameters , r = || x -x i || 2 , and ( ) · + denotes the truncated power function. Note that Gaussian (G) is a globally supported function of infinite smoothness, whereas Wendland's function (W2) is compactly supported one with degree of smoothness 2. For further details, see [14].

## 2.2. Partition of unity interpolant

The partition of unity method was firstly suggested in [2, 9] in the context of meshfree Galerkin methods for the solution of partial differential equations (PDEs), but now it is also commonly used in the field of approximation theory (see [14]). In particular, this approach enables us to decompose a large problem into many small problems, and at the same time ensures that the accuracy obtained for the local fits is carried over to the global one. In fact, the partition of unity method we here consider can be thought as a modified Shepard's method with higher-order data, where local approximations R j are RBFs. Similar local approaches involving Shepard's type methods were considered in [1, 8, 10, 11, 12].

Thus, the partition of unity method consists in partitioning the open and bounded domain Ω ⊆ R 2 into d subdomains Ω j such that Ω ⊆ ⋃ d j =1 Ω j with some mild overlap among the subdomains. At first, we choose a partition of unity, i.e. a family of compactly supported, non-negative, continuous functions W j with supp( W j ) ⊆ Ω j such that

$$\sum _ { j = 1 } ^ { d } W _ { j } ( { \mathbf x } ) = 1, \quad \mathbf x \in \Omega.$$

Then, we can consider the global interpolant

$$\mathcal { I } ( \mathbf x ) = \sum _ { j = 1 } ^ { d } \mathcal { R } _ { j } ( \mathbf x ) W _ { j } ( \mathbf x ),$$

where the local radial basis function

$$\mathcal { R } _ { j } ( \mathbf x ) = \sum _ { k = 1 } ^ { m _ { j } } c _ { k } \phi ( | | \mathbf x - \mathbf x _ { k } | | _ { 2 } )$$

is obtained by solving a local interpolation problem, which is constructed using the m j nodes belonging to each subdomain Ω . j Note that if the local approximants satisfy the interpolation conditions at node x i , i.e. R j ( x i ) = f ( x i ), then the global approximant also interpolates at this node, i.e. I ( x i ) = f ( x i ), for i = 1 2 , , . . . , n .

In order to be able to formulate error bounds we need some technical conditions. Then, we require the partition of unity functions W j to be k-stable , i.e. each W j ∈ C k ( R 2 ), j = 1 2 , , . . . , d , and for every multi-index β ∈ N m 0 with | β | ≤ k there exists a constant C β &gt; 0 such that

$$| | D ^ { \beta } W _ { j } | | _ { L _ { \infty } ( \Omega _ { j } ) } \leq C _ { \beta } / \delta _ { j } ^ { | \beta | },$$

where δ j = diam(Ω ). j

In accordance with the statements in [13] we require some additional regularity assumptions on the covering { Ω j } d j =1 . Therefore, setting X n = { x i , i = 1 2 , , . . . , n } ⊆ Ω, an open and bounded covering { Ω j } d j =1 is called regular for (Ω , X n ) if the following properties are satisfied:

- (a) for each x ∈ Ω, the number of subdomains Ω j with x ∈ Ω j is bounded by a global constant K ;
- (b) each subdomain Ω j satisfies an interior cone condition;
- (c) the local fill distances h X j , Ω j are uniformly bounded by the global fill distance h X n , Ω , where X j = X n ∩ Ω . j

Therefore, assuming that:

- · φ ∈ C k ν ( R 2 ) is a strictly positive definite function;
- · { Ω j } d j =1 is a regular covering for (Ω , X n );
- · { W j } d j =1 is k -stable for { Ω j } d j =1 ;

we have the following convergence result (see, e.g., [7, 14]), i.e., the error between f ∈ N φ (Ω), where N φ is the native space of φ , and its partition of unity interpolant (5) can be bounded by

$$| D ^ { \beta } f ( { \mathbf x } ) - D ^ { \beta } \mathcal { I } ( { \mathbf x } ) | \leq C h _ { \mathcal { X } _ { n }, \Omega } ^ { ( k + \nu ) / 2 - | \beta | } | f | _ { \mathcal { N } _ { \phi } ( \Omega ) },$$

for all x ∈ Ω and all | β | ≤ k/ 2.

## 3. Cell algorithm

In this section we present an efficient algorithm for bivariate interpolation of scattered data sets lying on the domain Ω = [0 1] , 2 ⊂ R 2 , which is based on the partition of unity method for constructing a global interpolant by blending RBFs as local approximants and using locally supported weight functions. It is efficiently implemented and optimized by connecting the interpolation method itself with an effective cell-based searching procedure. In particular, the implementation of this algorithm is based on the construction of a cell structure , which is obtained by partitioning the domain Ω in square cells, whose sizes strictly depend on the dimension of its subdomains. Such approach leads to important improvements in the searching processes of the nearest neighbour points compared to the searching techniques presented in [1, 3, 4].

## 3.1. Input and output

INPUT:

- · X n = ( { x , y i i ) , i = 1 2 , , . . . , n } , set of nodes;
- · F n = { f , i i = 1 2 , , . . . , n } , set of data values;
- · C d = (¯ { x , y i ¯ ) i , i = 1 2 , , . . . , d } , set of subdomain points (centres);
- · E s = (˜ { x , y i ˜ ) i , i = 1 2 , , . . . , s } , set of evaluation points.

## OUTPUT:

- · A s = {I (˜ x , y i ˜ ) i , i = 1 2 , , . . . , s } , set of approximated values.

## 3.2. Data partition phase

Stage 1. The set X n of nodes and the set E s of evaluation points are ordered with respect to a common direction (e.g. the y -axis), by applying a quicksort y procedure .

Stage 2. For each subdomain point (¯ x , y i ¯ ), i i = 1 2 , , . . . , d , a local circular subdomain is constructed, whose half-size (the radius) depends on the subdomain number d , that is

$$\delta _ { s u b d o m } = \sqrt { \frac { 2 } { d } }.$$

This value is suitably chosen, supposing to have a nearly uniform node distribution and assuming that the ratio n/d ≈ 4.

Stage 3. A double structure of crossed strips is constructed as follows:

i) a first family of q strips, parallel to the x -axis, is considered taking

$$q = \left \lceil \frac { 1 } { \delta _ { s u b d o m } } \right \rceil,$$

and a quicksort x procedure is applied to order the nodes belonging to each strip;

- ii) a second family of q strips, parallel to the y -axis, is considered.

Note that each of the two strip structures are ordered and numbered from 1 to q ; moreover, the choice in (7) follows directly from the side length of the domain Ω (unit square), that here is 1, and the subdomain radius δ subdom .

Stage 4. The domain (unit square) is partitioned by a cell-based structure consisted of q 2 square cells, whose length of the sides is given by δ cell ≡ δ subdom . Then, the following structure is considered:

- · the sets X n , C d and E s are partitioned by the cell structure into q 2 subsets X n k , C d k and E p k , k = 1 2 , , . . . , q 2 ,

where n k , d k and p k are the number of points in the k -th cell.

## 3.3. Localization phase

Stage 5. In order to identify the cells to be examined in the searching procedure, we consider two steps as follows:

- (A) since δ cell ≡ δ subdom , the ratio between these quantities is denoted by i ∗ = δ subdom /δ cell = 1. So the number j ∗ = (2 i ∗ + 1) 2 of cells to be examined for each node is 9.
- (B) for each cell k = [ v, w ], u, v = 1 2 , , . . . , q , a cell-based searching procedure is considered, examining the points from the cell [ v -i ∗ , w -i ∗ ] to the cell [ v + i ∗ , w + i ∗ ]. Note that if v -i ∗ &lt; 1 and/or w -i ∗ &lt; 1, or v + i ∗ &gt; q and/or w + i ∗ &gt; q , then we set v -i ∗ = 1 and/or w -i ∗ = 1, and v + i ∗ = q and/or w + i ∗ = . q

Then, after defining which and how many cells are to be examined, the cell-based searching procedure is applied:

- · for each subdomain point of C d k , k = 1 2 , , . . . , q 2 , to determine all nodes belonging to a subdomain. The number of nodes of the subdomain centred at (¯ x , y i ¯ ) is counted and stored in i m i , i = 1 2 , , . . . , d ;
- · for each evaluation point of E p k , k = 1 2 , , . . . , q 2 , in order to find all those belonging to a subdomain of centre (¯ x , y i ¯ ) and radius i δ subdom . The number of subdomains containing the i -th evaluation point is counted and stored in r i , i = 1 2 , , . . . , s .

## 3.4. Evaluation phase

Stage 7. A local approximant R j ( x, y ) and a weight function W x,y j ( ), j = 1 2 , , . . . , d , is found for each evaluation point.

Stage 8. Applying the global fit (5), the surface can be approximated at any evaluation point (˜ x, y ˜) ∈ E s .

## 4. Complexity

The partition of unity algorithm involves the use of the standard sorting routine quicksort, which requires on average a time complexity O ( M log M ), where M is the number of nodes to be sorted. Specifically, we have a data partition phase consisting of building the data structure, where the computational cost is:

- · O ( n log n ) for sorting all n nodes;
- · O ( s log s ) for sorting all s evaluation points.

Moreover, in order to compute the local RBF interpolants, we have to solve d linear systems and the cost is:

- · O ( m 3 i ), i = 1 2 , , . . . , d , where m i is the number of nodes in the i -th subdomain.

Then, for the k -th evaluation point of E s the cost is:

$$\bullet \ r _ { k } \cdot \mathcal { O } ( m _ { i } ), \, i = 1, 2, \dots, d, \, k = 1, 2, \dots, s.$$

Finally, the algorithm requires 3 n , 3 d and 3 s storage requirements for the data, and m i , i = 1 2 , , . . . , d , locations for the coefficients of each local RBF interpolant.

## 5. Numerical results

In this section we present some tests to verify performance and effectiveness of the cell-based partition algorithm on scattered data sets. The code is implemented in C/C++ language, while numerical results are carried out on a Intel Core 2 Duo Computer (2.1 GHz). In the experiments we consider a node distribution with n = 4225 16641 66049 uniformly random Halton nodes generated , , by using the program given in [15]. The partition of unity algorithm is run considering d = 1024 4096 16384 subdomain points and , , s = 33 × 33 evaluation (or grid) points, which are contained in the unit square Ω = [0 , 1] 2

The performance of the interpolation algorithm is verified taking the data values by Franke's test function

$$& f ( x, y ) = \frac { 3 } { 4 } \exp \left [ - \frac { ( 9 x - 2 ) ^ { 2 } + ( 9 y - 2 ) ^ { 2 } } { 4 } \right ] + \frac { 3 } { 4 } \exp \left [ - \frac { ( 9 x + 1 ) ^ { 2 } } { 4 9 } - \frac { 9 y + 1 } { 1 0 } \right ] \\ & \quad + \frac { 1 } { 2 } \exp \left [ - \frac { ( 9 x - 7 ) ^ { 2 } + ( 9 y - 3 ) ^ { 2 } } { 4 } \right ] - \frac { 1 } { 5 } \exp \left [ - ( 9 x - 4 ) ^ { 2 } - ( 9 y - 7 ) ^ { 2 } \right ]. \\ & \text{Moreover, since we are concerned to noint, out the effectiveness of the nro-}$$

Moreover, since we are concerned to point out the effectiveness of the proposed algorithm, in Table 1 we compare CPU times (in seconds) obtained by running the cell algorithm described in Section 3, and the strip algorithm proposed in [6]. This comparison highlights the high efficiency of the cell algorithm, which gives us a considerable saving of time.

Table 1: CPU times (in seconds) obtained by running the cell algorithm ( t cell ) and the strip algorithm ( t strip ).

|     n |     d |   t cell | t strip   |
|-------|-------|----------|-----------|
|  4225 |  1024 |      0.3 | 0 . 4     |
| 16641 |  4096 |      0.8 | 1 . 3     |
| 66049 | 16384 |      2.6 | 6 . 5     |

Then, in order to test accuracy of the local algorithm, in Table 2 we report the Root Mean Square Errors (RMSEs), i.e.

$$R M S E = \sqrt { \frac { 1 } { s } \sum _ { i = 1 } ^ { s } | f ( { \mathbf x } _ { i } ) - \mathcal { I } ( { \mathbf x } _ { i } ) | ^ { 2 } }. \\ \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \ dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots.$$

The error computation is achieved by considering both globally and locally supported RBFs for suitable values of the shape parameters, i.e., α 2 = 50 for φ G , and c = 1 for φ W 2 . We observe that the local scheme turns out to be accurate, even if we do not consider the optimal values for the shape parameters, namely those values for which we get the best possible results. However, these choices give a good compromise among accuracy and stability.

Table 2: RMSEs with α 2 = 50 and c = 1 for Franke's function.

| n     | 4225          | 16641         | 66049         |
|-------|---------------|---------------|---------------|
| φ G   | 3 . 0045E - 4 | 2 . 8214E - 5 | 1 . 5200E - 6 |
| φ W 2 | 2 . 2145E - 4 | 5 . 3127E - 5 | 9 . 3027E - 6 |

Then, in Figure 1 we plot the behavior of the RMSEs by varying values of the shape parameters for Franke's function. These graphs (and other ones obtained considering different test functions we omit for brievity) point out that, if an optimal search of the shape parameters was performed, in some cases the results of accuracy reported in this section could be improved of one or even two orders of magnitude. Note that each evaluation is carried out by choosing equispaced values of the shape parameter with α 2 ∈ [1 , 100], and c ∈ [0 1 . , 2].

By analyzing numerical tests and the related pictures, we observe that Wendland's function φ W 2 has greater stability than Gaussian φ G and good accuracy. However, these graphs give an idea on stability and enable us to choose 'sure' values for the shape parameters. These tests confirm theoretical results and suggest to use basis functions with a moderate order of smoothness, thus avoiding the well-known ill-conditioning problems of infinitely smooth RBFs, in particular if we are dealing with a very large number of nodes.

Figure 1: RMSEs obtained by varying α 2 and c for Franke's function.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0183v1-figure-9794f118-67cd-4680-b113-024307e01e5a.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0183v1-figure-b796abec-7000-4d42-95b1-3cd96ce3160c.jpg)

## References

- [1] G. Allasia, R. Besenghi, R. Cavoretto &amp; A. De Rossi, Scattered and track data interpolation using an efficient strip searching procedure, Appl. Math. Comput. 217 (2011), 5949-5966.
- [2] I. Babuˇ ska &amp; J. M. Melenk, The partition of unity method, Internat. J. Numer. Methods. Engrg. 40 (1997), 727-758.
- [3] R. Cavoretto &amp; A. De Rossi, Fast and accurate interpolation of large scattered data sets on the sphere, J. Comput. Appl. Math. 234 (2010), 15051521.
- [4] R. Cavoretto &amp; A. De Rossi, Spherical interpolation using the partition of unity method: an efficient and flexible algorithm, Appl. Math. Lett. 25 (2012), 1251-1256.
- [5] R. Cavoretto, A unified version of efficient partition of unity algorithms for meshless interpolation, in: T. E. Simos et al. (Eds.), Proceedings of the ICNAAM 2012 , AIP Conf. Proc., vol. 1479, Amer. Inst. Phys., Melville, New York, 2012, pp. 1054-1057.
- [6] R. Cavoretto, Partition of unity algorithm for two-dimensional interpolation using compactly supported radial basis functions, 2012, submitted for publication.
- [7] G. E. Fasshauer, Meshfree Approximation Methods with MATLAB , World Scientific Publishers, Singapore, 2007.
- [8] D. Lazzaro, L. B. Montefusco, Radial basis functions for the multivariate interpolation of large scattered data sets, J. Comput. Appl. Math. 140 (2002), 521-536.

- [9] J. M. Melenk &amp; I. Babuˇ ska, The partition of unity finite element method: basic theory and applications, Comput. Methods. Appl. Mech. Engrg. 139 (1996) 289-314.
- [10] R. J. Renka, Multivariate interpolation of large sets of scattered data, ACM Trans. Math. Software 14 (1988), 139-148.
- [11] R. J. Renka, Algorithm 660: QSHEP2D: Quadratic Shepard method for bivariate interpolation of scattered data, ACM Trans. Math. Software 14 (1988), 149-150.
- [12] W. I. Thacker, J. Zhang, L. T. Watson, J. B. Birch, M. A. Iyer &amp; M. W. Berry, Algorithm 905: SHEPPACK: Modified Shepard algorithm for interpolation of scattered multivariate data, ACM Trans. Math. Software 37 (2010), Art. 34, 1-20.
- [13] H. Wendland, Fast evaluation of radial basis functions: Methods based on partition of unity, in: C. K. Chui et al. (Eds.), Approximation Theory X: Wavelets, Splines, and Applications , Vanderbilt Univ. Press, Nashville, TN, 2002, pp. 473-483.
- [14] H. Wendland, Scattered Data Approximation , Cambridge Monogr. Appl. Comput. Math., vol. 17, Cambridge Univ. Press, Cambridge, 2005.
- [15] T.-T. Wong, W.-S. Luk &amp; P.-A. Heng, Sampling with Hammersley and Halton points, J. Graphics Tools 2 (1997), 9-24.
Roberto Cavoretto
Two-dimensional interpolation using a cell-based searching procedure
In this paper we present an efficient algorithm for bivariate interpolation,
which is based on the use of the partition of unity method for constructing a
global interpolant. It is obtained by combining local radial basis function
interpolants with locally supported weight functions. In particular, this
interpolation scheme is characterized by the construction of a suitable
partition of the domain in cells so that the cell structure strictly depends on
the dimension of its subdomains. This fact allows us to construct an efficient
cell-based searching procedure, which provides a significant reduction of CPU
times. Complexity analysis and numerical results show such improvements on the
algorithm performances.