Thin plate spline interpolation and fitting

ALGLIB, a free and commercial open source numerical library, provides the best open-source implementation of thin plate splines (TPSs) and related algorithms (biharmonic splines and multiquadrics). ALGLIB is available in multiple programming languages, including C++, C#, Java, and Python.

Our multi-threaded, SIMD-optimized implementation of the TPS algorithm can interpolate and smooth scattered spatial data, compute model values and derivatives, and convert scattered data into a gridded form. Unlike other open-source implementations, it is scalable to big datasets with hundreds of thousands of points.

Using the DDM-RBF solver

Programming languages supported

ALGLIB supports many programming languages, including C++, C#, Java, Python, and others:

• C++ version. Highly optimized, portable, self-contained library (x86/x64/ARM; works on Windows, Linux, and POSIX systems)
• C# version. Unified C# API for two backends: a fully managed C# implementation and a high-performance native computational core written in C
• Java, Python, and other languages. A wrapper for a high-performance native computational core written in C. Seamless integration with the target programming language

A distinctive feature of ALGLIB is that it provides the same API in all programming languages. This is achieved through our exclusive automatic code translation and wrapper generation technology.

RBF types supported

The DDM-RBF solver is a fast large-scale model construction algorithm employing state-of-the-art mathematical methods. It can handle hundreds of thousands of points and supports the following basis functions:

• Thin plate splines. The TPS kernel f=r²·log(r) is the best default choice, as there are no parameters to tune. It is C1 continuous and offers good numerical stability. This option is activated by the rbfsetalgothinplatespline function.
• Biharmonic splines. The biharmonic kernel f=r also provides good numerical stability and has no tunable parameters. Unlike thin plate splines, however, it is merely C0 continuous. This option is activated by the rbfsetalgobiharmonic function.
• Multiquadrics. The multiquadric kernel f=sqrt(r²+a²) is infinitely smooth when a is nonzero, but it requires well-separated data points. This option is activated by the rbfsetalgomultiquadricauto function (when a is proportional to the mean interpoint distance) or the rbfsetalgomultiquadricmanual function (requiring the manual choice of a).
• Gaussian RBFs. The rbf sub-package also supports Gaussian RBFs, but these RBFs are handled by the hierarchical RBF (HRBF) algorithm, which takes a completely different approach to accelerate computations (sparse linear algebra). This algorithm is discussed in a separate article.

API overview

RBF interpolation functionality is provided by the rbf subpackage (the link points to the ALGLIB Reference Manual entry for the subpackage).

The first step is model construction, which involves the following steps:

• the model object is initialized using rbfcreate
• the dataset is added using rbfsetpoints or rbfsetpointsandscales
• the solver is configured with one of the solver configuration functions (see the previous section)
• the model is built by rbfbuildmodel

After the model is built, the following operations can be performed:

• one-off RBF evaluation with rbfcalc and its specialized versions (convenience wrappers for 1..3-dimensional problems or thread-safe rbftscalcbuf)
• gridded RBF evaluation with rbfgridcalc2v or rbfgridcalc3v
• model serialization with rbfserialize/rbfunserialize
• model unpacking with rbfunpack (which extracts model coefficients)

Interpolation vs. smoothing

All RBF algorithms in ALGLIB accept the smoothing parameter LambdaV. A zero value means that we solve an exact RBF interpolation problem. Such problems arise in many areas, including geospatial applications, surface reconstruction, and others.

When the smoothing parameter LambdaV is nonzero, the solver fits a smoothing RBF spline to the data. Smoothing may be used to actually smooth the data, e.g., to suppress noise. Another reason for employing smoothing is its regularizing effect on the model.

Generally, RBF interpolation requires nodes to be distinct and well separated, but industrial data often contain badly separated nodes. Some RBF kernels (Gaussian and multiquadric) are extremely sensitive to such nodes. Other kernels are also sensitive, albeit to a lesser extent. In such cases, specifying a small, nonzero LambdaV will help the solver deal with the data.

Mathematical background: Fast algorithms for polyharmonic splines and RBFs

The most time- and memory-consuming element of RBF interpolation is model construction. Our DDM-RBF solver accelerates its computations by using the following:

• approximate cardinal basis function (ACBF) preconditioning, which transforms the original basis functions into a better-conditioned internal form
• a domain decomposition method that separates a data set into smaller pieces and handles them separately
• parallelism and SIMD support (in both the C++ and C# versions of ALGLIB)

These improvements result in a parallel large-scale interpolation/fitting method with O(N) memory requirements and O(N²) running time.

Note #1
For comparison, a straightforward solution with a dense linear solver (as in SciPy and other open-source implementations) requires O(N³) time and O(N²) memory.