[AM-04-071] Splines and Radial Basis Functions Interpolation

Mitasova, H. and Mitas, L. (2025). Splines and Radial Basis Functions Interpolation. The Geographic Information Science & Technology Body of Knowledge (2025 Edition), John P. Wilson (ed.). DOI: 10.22224/gistbok/2025.1.25.

1. Introduction
2. RBF Splines Foundations
3. Spline Interpolation Properties
4. Multivariate and Partial Splines
5. Bivariate Multiquadrics
6. Summary

1. Introduction

Spatial interpolation methods (Raposo 2025) use the measured values at given locations to estimate the values at unsampled locations (Figure 1). The measured data are often irregularly distributed points, with significant noise and potential discontinuities. The datasets, acquired by modern mapping technologies such as LiDAR can be very large, obtained from various sources with variable accuracies. The modeled fields, such as elevations, precipitation or chemical concentrations, are usually very complex. Many interpolation and approximation methods were developed to address these challenges and compute regular grid (raster) representation of modeled fields from scattered point data. Methods based on linear superposition of radial basis functions (RBF) centered at the data points include multivariate splines and multiquadrics.  Spline functions minimize the smoothness seminorm or roughness penalty leading to a smooth surface that passes through or close to the given data points.

It is important to note that the RBF class of splines described in this section are, in general, different from a rich class of piecewise polynomial splines (Wahba 1990) such as locally refined B-splines (LR B-splines, Skytt et al., 2015)  or non-uniform rational B-splines (NURBS, Fisher and Wales, 1992) commonly used in computer-aided design for interpolation of line features or creating smooth 3-D surfaces with meshes (3D Parametric Surfaces, forthcoming). 

2. Radial Basis Functions (RBF) Splines Foundations

Spline interpolation functions were derived using a variational approach based on the assumption that the interpolation function should pass through, or close to, the data points and, at the same time, should be as smooth as possible. These two requirements are combined into a single condition of minimizing the sum of the deviations from the N measured points and the smoothness seminorm (roughness penalty in Hutchinson 1989):     

Equation 1

where  is location of a given point   is the measured value in are positive weights (smoothing parameters), I(F) denotes the smoothness seminorm (Table 1) and F(r) is the interpolation function. The solution of Equation (1) can be expressed as a sum of two components (Wahba 1990):

Equation 2

where T(r) is a ‘trend’ function and is a radial basis function which has a form dependent on the choice of smoothness seminorm I(F). Smoothness seminorm minimization can be performed analytically (Mitas and Mitasova, 1988) leading to an explicit solution for  (Table 2). The coefficients 𝜆 are obtained by solving a system of linear equations which makes the computational cost proportional to N³. Various segmentation techniques were developed to reduce the computational cost and make the spline interpolation feasible for massive data sets, for example, overlapping rectangles in Franke 1985 and quadtrees in Mitasova et al., 1995. 

Table 1. Examples of bivariate spline functions, their corresponding smoothness seminorms and Euler-Lagrange equations.  denote first and second order partial derivative according to x, etc. Variable w denotes the weight of the given term.

Alternatively, the minimization in Equation (1) can be carried out numerically, for example, by using a finite difference multi-grid iteration method (Hutchinson 1989). Numerical solution enables incorporation of stream enforcement and other topographic features (Hutchinson 1989) with the computational cost proportional to the number of interpolated grid points. To interpolate large climate data sets, Hancock and Hutchinson (2006) discretized the bivariate thin plate smoothing spline equations using hierarchical biquadratic B-splines, and used a nested grid multigrid procedure to solve the system while optimizing the smoothness parameter using generalized cross-validation.

3. Spline Interpolation Properties

Spline functions minimize smoothness seminorms with weighted derivatives shown in Table 1. For example, thin plate spline (TPS) function (Duchon 1976) minimizes surface curvature (the respective I(F) includes second order derivatives, see Table 1) and simulates a thin steel sheet forced to pass through the data points. The equilibrium shape of the sheet minimizes the bending energy closely related to the surface curvature.

The stiffness of the simulated sheet in TPS can be controlled by adding the first order derivatives into the seminorm I(F), leading to a TPS with tension (Franke 1985; Hutchinson 1989; Mitas and Mitasova 1988, Table 1). Change in the tension tunes the surface from a stiff plate into an elastic membrane and in the limit of an infinite tension the surface resembles a rubber sheet with cusps at the data points (Figure 2). Adding higher order derivatives into I(F), leads to a function with regular second order derivatives (Mitas and Mitasova 1988) suitable for deriving slope, aspect and curvatures of the interpolated surface simultaneously with interpolation.  

Table 2. Expressions for the bivariate radial basis functions corresponding to models in Table 1.  is the modified Bessel function of zero order,  is the Euler constant, and is the exponential integral function. 

Regularized spline with tension (RST) was designed to synthesize the desired properties into a single function by including the sum of all derivatives up to infinity with rapidly decreasing weights in I(F) (Table 1). The resulting surface has a tunable tension (sheet stiffness) and regular derivatives of all orders and therefore is suitable for computing topographic parameters (gradient and curvatures) simultaneously with interpolation (Figure 3, Mitasova et al., 2005). The tension parameter can be further generalized to a tensor which enables modeling of anisotropy (Goovaerts 2019). Smoothing weights in Equation (1) can be spatially variable and are often based on data accuracy, a useful property for modeling surfaces from noisy data.

Effect of smoothing on the resulting surface can be mapped by computing a deviation map, derived from differences between the given and interpolated values in the given points (Figure 3a, b). Predictive accuracy of interpolation can be computed using cross-validation (Figure 4c). The tension and smoothing parameters can be selected empirically, based on the knowledge of the modeled phenomenon, or automatically by minimization of the predictive error estimated by a cross-validation procedure (Hofierka et al. 2002).

Several authors have demonstrated that splines are formally equivalent to universal kriging (Goovaerts 2019) with the choice of the covariance function determined by the smoothness seminorm I(F) (e.g., Wahba 1990). Therefore, many of the geostatistical concepts can be utilized  within the spline framework.

5. Multivariate and Partial Splines

function was applied for interpolation of chemical concentrations in water bodies or groundwater (Mitasova et al. 1995, see also animations showing concentration of nitrogen in Chesapeake Bay over a period of several months or concentration of pollutants in groundwater over several years).

The additional variable(s) can be used for interpolation of phenomena that are influenced by another effect or external condition, such as interpolation of precipitation  with influence of topography (Hofierka et al. 2002). Impact of elevation included as an additional variable in interpolation of annual precipitation in mountains of North Carolina is illustrated in Figure 5.  Incorporation of dependence on additional variables (Hutchinson 1998; Wahba 1990) has been widely used for modeling precipitation patterns including high resolution precipitation in Australia using ANUSPLIN (Nguyen et al, 2024).

The variational approach to spatial interpolation also offers numerous possibilities to incorporate additional conditions such as value constraints, prescribed derivatives at given points, and integral constraints (Wahba 1990).

5. Bivariate Multiquadrics

Multiquadric methods (e.g., Carlson and Foley 1991) are similar in formulation and comparable in performance with splines with radial basis functions as follows:

Equation 3

In three dimensions, in the limit of b→0, the basis functions are solutions of biharmonic and harmonic equations, respectively. Multiquadric methods offer high accuracy, differentiability, d-dimensional formulation, and, with segmentation, applicability to large datasets. Multiquadrics with anisotropy have been employed to model spatial distribution of groundwater seepage (Kennedy et al. 2008).

6. Summary

Interpolation by multi-variate splines combines convenient analytical properties based on  physical models such as low-bending energy elastic plates or higher-dimensional objects. Often, interpolated phenomena result from processes which indeed minimize some type of energy, with a typical example being a settled terrain with its slowly-evolving balance between gravitation, soil cohesion and impact of climate. For these cases, splines have proven to be rather effective. In addition, using auxiliary variables can help to describe phenomena where indirect impacts by another effect/phenomenon can be made more explicit.  Moreover, splines are very useful for calculations of local geometry and differential analysis which is often used as input to various process-based models. Variety of preprocessing tools such as data filtering, partitioning of large and complex data sets open further possibilities for applications.

• GRASS  v.surf.rst, v.vol.rst tools for bi-variate and tri-variate spline interpolation with simultaneous computation of gradients and curvatures
• Spline interpolation topic with assignments (Geospatial Modeling course at North Carolina State University) https://ncsu-geoforall-lab.github.io/geospatial-modeling-course/topics/spatial_interpolation.html
• ANUSPLIN 4: software package for analysis and interpolation of noisy multi-variate data using thin plate smoothing splines