Comparative investigation of parallel spatial interpolation algorithms for building large-scale digital elevation models

Original MLS Interpolation Algorithm

The MLS approximation is used to approximate field variables and their derivatives. In a domain Ω, the MLS approximation $f^h\left(x\right)$ of the field variable $f\left(x\right)$ in the vicinity of a point $\bar{x}$ is given as (1)${\displaystyle f^h\left(x\right)=\sum _{j=1}^m{p}_j\left(x\right)\cdot a_j\left(\bar{x}\right)=P^T\left(x\right)\cdot a\left(\bar{x}\right)}$ where $p_j\left(x\right),j=1$ , 2, …, m is a complete basis function with coefficients $a_j\left(\bar{x}\right)$. At each point $\bar{x}$, $a_j\left(\bar{x}\right)$ is chosen to minimize the weighted residual L₂ − norm (L₂ − norm refers to ${∥x∥}_2$, where $x={\left[x_1,x_2,\dots ,x_n\right]}^T$, and ${∥x∥}_2=\sqrt{\left({\left|x_1\right|}^2+{\left|x_2\right|}^2+{\left|x_3\right|}^2+\cdots +{\left|x_n\right|}^2\right)}$): (2)${\displaystyle J=\sum _{I=1}^{N}w\left(\bar{x}-x_I\right){\left[P^T\left(x_I\right)a\left(\bar{x}\right)-f_I\right]}^2}$ where N is the number of nodes in the compact-supported neighborhood of $\bar{x}$ and f_I refers to the nodal parameter of f at x = x_I. Nodes refer to data points in the compact-supported neighborhood of $\bar{x}$. Compact-supported, i.e., point $\bar{x}$ is only related to the nodes of its neighborhood, x_I is one of the nodes in the compact-supported neighborhood. And $w\left(x-x_k\right)$ is the compact-supported weight function. The most commonly used weight functions are the spline functions, for example, the cubic spline weight function (Eq. (3)): (3)${\displaystyle w\left(\bar{s}\right)=\left\{\begin{array}{c}\frac{2}{3}-4{\bar{s}}^2+4{\bar{s}}^3,\hfill \\ \frac{4}{3}-4\bar{s}+4{\bar{s}}^2-\frac{4}{3}{\bar{s}}^3,\hfill \\ 0,\hfill \\ \hfill \end{array}\right.\begin{array}{c}\bar{s}\le \frac{1}{2}\hfill \\ \frac{1}{2}<\bar{s}\le 1\hfill \\ \bar{s}>1\hfill \end{array}}$ where $\bar{s}=\frac{s}{s_{\max }}$ and $s=\bar{x}-x_I$.

The minimum of J with respect to $a\left(\bar{x}\right)$ gives the standard form of MLS approximation: (4)${\displaystyle f^h\left(x\right)=\sum _{I=1}^N{\varphi }_I\left(x\right)f_I=\Phi \left(x\right)F.}$

Orthogonal MLS interpolation algorithm

For a given polynomial basis function $p_i\left(x\right),\begin{array}{c}\hfill \hfill \\ \hfill \hfill \end{array}$i = 1 , 2, ⋅⋅⋅, m, there is an orthonormal basis function $q_i\left(x,\bar{x}\right)$ that satisfies: ${\displaystyle q_1\left(x,\bar{x}\right)=p_1\left(x\right)}$ (5)${\displaystyle q_i\left(x,\bar{x}\right)=p_i\left(x\right)-\sum _{j=1}^{i-1}{\alpha }_{ij}\left(x,\bar{x}\right)q_j\left(x,\bar{x}\right),i=2,3,\cdot \cdot \cdot ,m}$ where ${\alpha }_{ij}\left(x,\bar{x}\right)$ is the coefficient that makes $q_i\left(x,\bar{x}\right)$ perpendicular to $q_j\left(x,\bar{x}\right)$. (6)${\displaystyle {\alpha }_{ij}\left(\bar{x}\right)=\frac{\sum _{k=1}^N{w}_k\left(\bar{x}\right)p_i\left(x_k\right)q_j\left(x_k,\bar{x}\right)}{{\sum }_{k=1}^N{w}_k\left(\bar{x}\right)q_j^2\left(x_k,\bar{x}\right)}}$

Because the coefficient matrix is a diagonal matrix, the solution for $a_i\left(x\right)$ does not require matrix inversion, i.e., (7)${\displaystyle a_i\left(\bar{x}\right)=\frac{{\sum }_{k=1}^N{w}_k\left(\bar{x}\right)q_i\left(x_k,\bar{x}\right)f_k}{{\sum }_{k=1}^N{w}_k\left(\bar{x}\right)q_i^2\left(x_k,\bar{x}\right)}}$ where a_i and $a_j\left(\bar{x}\right)$ (Eq. (1)) have the same definition. f_k and f_I (Eq. (2)) have the same definition, i.e., the nodal parameter of f at x = x_k. Finally, a_i and the orthonormal basis function $q_i\left(x,\bar{x}\right)$ are fitted into Eq. (1) to obtain the orthogonal MLS approximation $f^h\left(x\right)$.

When the number or order of basis functions increases, only a_m+1 and α_m+1 need to be calculated in Gram–Schmidt orthogonalization (Steve, 2011); recalculation of all entries in the coefficient matrix is not needed. This could reduce the computational cost and the computational error.

Lancaster’s MLS interpolation algorithm

A singular weight function is adopted to make the approximation function $f^h\left(x\right)$ constructed by the interpolation type MLS method satisfy the properties of the Kronecker δ function: (8)${\displaystyle \omega \left(x,x_k\right)=\left\{\begin{array}{c}{∥\left(x-x_k\right)∕{\rho }_k∥}^{\text{-}\alpha },\hfill \\ 0,\hfill \\ \hfill \end{array}\right.\begin{array}{c}∥x-x_k∥\le {\rho }_k\hfill \\ ∥x-x_k∥>{\rho }_k\hfill \end{array}}$

Let $p_0\left(x\right)\equiv 1,p_1\left(x\right),\dots ,p_{\bar{m}}\left(x\right)$ denote the basis function used to construct the approximation function, where the number of basis functions is $\bar{m}+1$. To implement the interpolation properties, a new set of basis functions is constructed for a given basis function. First, p₀(x) are standardized, i.e., (9)${\displaystyle {\tilde{p}}_0\left(x,\bar{x}\right)=\frac{1}{{\left[\sum _{k=1}^N\omega \left(x,x_k\right)\right]}^{1∕2}}}$

Then, we construct a new basis function of the following form: (10)${\displaystyle {\tilde{p}}_i\left(x,\bar{x}\right)=p_i\left(\bar{x}\right)-\sum _{k=1}^N\frac{\omega \left(x,x_k\right)}{\sum _{l=1}^N\omega \left(x,x_l\right)}{P}_i\left(x_k\right),i=1,2,\dots ,\bar{m}.}$

Variant A of Shepard’s interpolation algorithm

First, select the influence radius R > 0 and let the weight function be (14)${\displaystyle w\left(r\right)=\left\{\begin{array}{c}\frac{1}{r},\hfill \\ \frac{27}{4}{\left(\frac{r}{R}-1\right)}^2,\hfill \\ 0,\hfill \\ \hfill \end{array}\right.\begin{array}{c}0<r\le \frac{R}{3}\hfill \\ \frac{R}{3}<r\le R\hfill \\ r>R\hfill \end{array}}$

Then, a variation of Shepard’s interpolation will be obtained.

Variant B of Shepard’s Interpolation Algorithm

When employing the following weight function (Eq. (15)), a new variation of Shepard’s interpolation will be obtained. (15)${\displaystyle w\left(\bar{s}\right)=\left\{\begin{array}{c}\frac{2}{3}-4{\bar{s}}^2+4{\bar{s}}^3,\hfill \\ \frac{4}{3}-4\bar{s}+4{\bar{s}}^2-\frac{4}{3}{\bar{s}}^3,\hfill \\ 0,\hfill \\ \hfill \end{array}\right.\begin{array}{c}\bar{s}\le \frac{1}{2}\hfill \\ \frac{1}{2}<\bar{s}\le 1\hfill \\ \bar{s}>1\hfill \end{array}}$

Inverse Distance Weighted (IDW) interpolation algorithm

If the weight function is selected as (16)${\displaystyle w_i\left(x\right)=\frac{1}{d{\left(x,x_i\right)}^{\alpha }}}$ the IDW interpolation is obtained. Typically, α = 2 in the standard IDW. Where $d\left(x,x_i\right)$ is the distance between the interpolation point x_i and the nearby discrete point x.

AIDW interpolation algorithm

The Adaptive Inverse Distance Weighted (AIDW) is an improved version of the standard IDW (Shepard, 1968) originated by Lu & Wong (2008). The distance-decay parameter α is no longer a prespecified constant value but is adaptively adjusted for a specific unknown interpolated point according to the distribution of the nearest neighboring data points.

The parameter α is taken as (17)${\displaystyle \alpha \left({\mu }_R\right)=\left\{\begin{array}{c}{\alpha }_1,\hfill \\ {\alpha }_1\left[1-5\left({\mu }_R-0.1\right)\right]+5{\alpha }_2\left({\mu }_R-0.1\right),\hfill \\ 5{\alpha }_3\left({\mu }_R-0.3\right)+{\alpha }_2\left[1-5\left({\mu }_R-0.3\right)\right],\hfill \\ {\alpha }_3\left[1-5\left({\mu }_R-0.5\right)\right]+5{\alpha }_4\left({\mu }_R-0.5\right),\hfill \\ 5{\alpha }_5\left({\mu }_R-0.7\right)+{\alpha }_4\left[1-5\left({\mu }_R-0.7\right)\right],\hfill \\ {\alpha }_5,\hfill \\ \hfill \end{array}\right.\begin{array}{c}0.0\le {\mu }_R\le 0.1\hfill \\ 0.1\le {\mu }_R\le 0.3\hfill \\ 0.3\le {\mu }_R\le 0.5\hfill \\ 0.5\le {\mu }_R\le 0.7\hfill \\ 0.7\le {\mu }_R\le 0.9\hfill \\ 0.9\le {\mu }_R\le 1.0\hfill \end{array}}$ (18)${\displaystyle {\mu }_R=\left\{\begin{array}{c}0,\hfill \\ 0.5-0.5cos\left[\pi \left(R\left(S_0\right)-R_{\min }\right)∕R_{\max }\right],\hfill \\ 1,\hfill \\ \hfill \end{array}\right.\begin{array}{c}R\left(S_0\right)\le R_{\min }\hfill \\ R_{\min }\le R\left(S_0\right)\le R_{\max }\hfill \\ R\left(S_0\right)\ge R_{\max }\hfill \end{array}}$ where the α₁, α₂, α₃, α₄, α₅ are the to-be-assigned five levels or categories of distance decay value. R_min or R_max refer to a local nearest neighbor statistic value, and R_min and R_max can generally be set to 0.0 and 2.0, respectively. Then, (19)${\displaystyle R\left(S_0\right)=\frac{2\sqrt{N∕A}}{k}\sum _{i=1}^k{d}_i}$ where N is the number of points in the study area, A is the area of the study region, k is the number of nearest neighbor points, d_i is the nearest neighbor distances and S₀ is the location of an interpolated point.

Zone 1 (Flat Zone)

The first selected region is located in Hengshui City, Hebei Province. The DEM of this region has the identifier ASTGTM_N37E115 and is derived from the Geospatial Data Cloud (http://www.gscloud.cn/). The location and elevation of this region is illustrated in Fig. 2. In the region, the highest elevation is 48 m and the lowest is 8 m. We translated the X coordinate by 348,000 and the Y coordinate by 4,130,000 to obtain a 20 km ×20 km square area centered on the origin. Five sets of benchmark test data were generated in this region; see Table 1.

Zone 2 (Rugged Zone)

The second selected region is located in Ganzi Tibetan Autonomous Prefecture, Sichuan Province. The DEM of this region has the identifier ASTGTM_N29E099 and is derived from the Geospatial Data Cloud (http://www.gscloud.cn/). The location and elevation of this region is illustrated in Fig. 2. In the region, the highest elevation is 5,722 m and the lowest is 3,498 m. We translated the X coordinate by 570,000 and the Y coordinate by 3,300,000 to obtain a 20 km ×20 km square area centered on the origin. Five sets of benchmark test data are generated in this region; see Table 2.

Computational efficiency of sequential implementations

As listed in Table 6, for the sequential version, when giving the same sets of data points and interpolation points, the order of computational time from fastest to slowest is: the Shepard’s interpolation method, the MLS interpolation, and the RBF interpolation. The computational time of Shepard’s interpolation method is approximately 20% less than the MLS interpolation method, and it is approximately 70% less than the computational time of the computational time of RBF interpolation method.

Computational efficiency of parallel implementations

As shown in Figs. 6– 11, the parallel version developed on multi-GPUs has the highest speedup in the three parallel versions. Except for the RBF interpolation method, the maximum speedups of other interpolation algorithms are greater than 45.

As shown in Figs. 12 and 13, for the parallel version developed on multi-GPUs, the order of the computational time from fastest to slowest is: the Shepard’s interpolation, the MLS interpolation, the RBF interpolation method. The computational time of Shepard’s interpolation method is 3%–30% less than the computational time of the MLS interpolation method, and it is 70%–85% less than the computational time of the RBF interpolation method.

In the flat zone

As shown in Fig. 14, for the flat region, when the data points are evenly distributed, the frequency statistical curve of the MLS is the highest when it is close to zero, the lowest when it is far away from zero, and the relative error distribution range is smaller, which means that the error of MLS method is small. The characteristics of the frequency statistical curve of Shepard’s method are completely opposite to those of MLS, which means that the error of MLS method is large. For the RBF interpolation algorithm, the characteristic of the frequency statistics curve is a transitional phase between those for the MLS and those for Shepard’s method. The above curve features and E (Table 8) illustrate that the interpolation accuracy is from high to low in this condition: the MLS interpolation algorithm, RBF, and Shepard’s interpolation method.

When the data points are normally distributed, the relative error distribution ranges of all the interpolation methods are larger than that for the uniformly distributed data points. As shown in Fig. 14, the characteristics of the frequency statistics curve of RBF are obvious, the frequency statistical curve of RBF is above the frequency statistical curves of MLS and Shepard’s method, which means that the error of RBF method is larger. The characteristics of frequency statistical curves of MLS and Shepard’s method are very similar, and the relative error distribution range of MLS is the largest. As listed in Table 8, in the flat zone, the accuracy of MLS is slightly higher than Shepard’s method when the data points are normally distributed.

In the rugged zone

As shown in Fig. 15, for the rugged region, regardless of whether the data points are uniformly distributed or normally distributed, the characteristics of frequency statistical curves of MLS, RBF and Shepard’s method are similar to those illustrated in Fig. 14. However, in Fig. 15B, it is a little different in that most of the frequency statistical curve of Shepard’s method is higher than the RBF’s. As listed in Table 9, the interpolation accuracy is from high to low: the MLS interpolation algorithm, RBF, and Shepard’s interpolation method.

According to the above Figures and Tables, some summary conclusions are obtained as follows:

For the same region, when the density of data points is almost the same, the interpolation accuracy when the data points are evenly distributed is higher than the interpolation accuracy when the data points are nonuniformly distributed.

As listed in Tables 5 and 7, when the data points are evenly distributed, the gap of the accuracy between the three variations of the MLS method, RBF, and Shepard’s interpolation methods increases with the decrease of point density.

As shown in Figs. 14 and 15, when the data points are nonuniformly distributed, the maximum relative errors of MLS is larger than other algorithms’, however, MLS method has lower NRMSE and E. A small number of larger relative errors has little effect on the overall interpolation accuracy. A large number of small and medium relative errors are the key to determine the interpolation accuracy of the algorithm.

As listed in Table 5, compared with the uniform distribution, when the points are nonuniformly distributed the difference in the accuracy of the interpolation algorithms is not as sensitive to the changes of point density.

Compared with the three variations of the MLS method and the RBF method, Shepard’s interpolation method is quite suitable for cases where the data points have a smooth trend. When interpolating for the data points with an undulating trend, the accuracy of Shepard’s interpolation method will be poor. When the density of data points is small, this rule becomes more obvious.

In the flat zone

As illustrated in Fig. 12, for the flat region, except for the kNNRBF, when the number of data points is not much different, the nonuniformly distributed data point set requires significantly more interpolation time than the uniformly distributed data point set, and with the increase of the number of points, interpolation time does increase as well.

As illustrated in Fig. 10, the speedups achieved by the RBF interpolation method is generally small, and its speedups are not much different in various cases. However, when the size of data point set is Size 1 and the data point set is nonuniformly distributed, the speedup of the RBF interpolation method is larger than other methods, which means that the benefits of parallelism are lower in this case.

As indicated above, the distribution pattern of data points strongly influences the interpolation efficiency.

In the rugged zone

As illustrated in Figs. 11 and 13, the running time and the speedups in the rugged region are almost the same as those in the flat region. In other words, the characteristics of the terrain elevation of data points have a weak influence on computational efficiency.