Add KernelInterpolation.jl extension (#73)

JoshuaLampert · andrewwinters5000 · web-flow · commit 65319da94cb7 · 2025-09-29T20:33:26.000+02:00
* add KernelInterpolation.jl extension

* add tests

* define kernel separately to be able to overwrite with trixi_include

* pass all kwargs in CubicBSpline

* format

* Apply suggestions from code review

Co-authored-by: Andrew Winters &lt;andrew.ross.winters@liu.se&gt;

* shuffle input data in tests

* actually run the tests

* add KernelInterpolation to test project

* add Random to test project

* fix test

* loosen tolerance

* loosen tolerance more

* add NEWS and README entries

* use evaluate_two_dimensional_interpolant

* fix link

* fix another link

* x- not y-direction

* try to include extension docstrings into references

* format

* need to use Base.get_extension

* load extensions

* try listing functions of extensions explicitly

* fix docstring of RBFInterpolation1D

* add name to AUTHORS

---------

Co-authored-by: Andrew Winters &lt;andrew.ross.winters@liu.se&gt;
diff --git a/.gitignore b/.gitignore
@@ -12,3 +12,5 @@ docs/src/changelog.md
 docs/src/code_of_conduct.md
 docs/src/contributing.md
 docs/src/index.md
+
+.vscode/
diff --git a/AUTHORS.md b/AUTHORS.md
@@ -18,6 +18,7 @@ The following people contributed major additions or modifications to TrixiBottom
 are listed in alphabetical order:
 
 * Maximilian D. Bertrand
+* Joshua Lampert
 * Hendrik Ranocha
 * Michael Schlottke-Lakemper
-* Andrew Winters
+* Andrew Winters
diff --git a/NEWS.md b/NEWS.md
@@ -9,6 +9,7 @@ for human readability.
 
 #### Added
 
+- Added scattered radial basis function interpolation using [KernelInterpolation.jl](https://github.com/JoshuaLampert/KernelInterpolation.jl) [#73]
 - Implementation of visualization routines `plot_topography` and `plot_topography_with_interpolation_knots`
 
 #### Changed
diff --git a/Project.toml b/Project.toml
@@ -10,13 +10,16 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
+KernelInterpolation = "95686e15-3c94-4443-8f08-76d06d509f7b"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 
 [extensions]
-TrixiBottomTopographyMakieExt = "Makie"
+KernelInterpolationExt = "KernelInterpolation"
+MakieExt = "Makie"
 
 [compat]
 Downloads = "1.6"
+KernelInterpolation = "0.3.5"
 LinearAlgebra = "1"
 Makie = "0.21, 0.22, 0.24"
 SparseArrays = "1"
diff --git a/README.md b/README.md
@@ -33,10 +33,11 @@ TrixiBottomTopography.jl does this by B-spline interpolation of the underlying d
 
 ## Functionalities
 
-This package contains the following three main functionalities:
+This package contains the following main functionalities:
 - Converting geographical data given in form of `.xyz` files, e.g., from the [DGM data set](https://www.opengeodata.nrw.de/produkte/geobasis/hm/) provided by [Geobasis NRW](https://www.bezreg-koeln.nrw.de/geobasis-nrw) to make it readable for TrixiBottomTopography.jl
 - Setting up a B-spline interpolation structure in one and two dimensions which contains all the relevant information to define a B-spline interpolation function with additional specifications
 - Using the B-spline structure to set up a B-spline interpolation function
+- Setting up a radial basis function interpolation based on scattered data leveraging [KernelInterpolation.jl](https://github.com/JoshuaLampert/KernelInterpolation.jl)
 
 A detailed description of the functionalities can be found in the [documentation](https://trixi-framework.github.io/TrixiBottomTopography.jl/stable/) to this package.
 
@@ -125,4 +126,4 @@ The full list of contributors can be found in [AUTHORS.md](AUTHORS.md).
 
 ## License and contributing
 TrixiBottomTopography.jl is published under the MIT license (see [LICENSE.md](LICENSE.md)).
-We are pleased to accept contributions from everyone, preferably in the form of a PR.
+We are pleased to accept contributions from everyone, preferably in the form of a PR.
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -2,6 +2,7 @@
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Changelog = "5217a498-cd5d-4ec6-b8c2-9b85a09b6e3e"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+KernelInterpolation = "95686e15-3c94-4443-8f08-76d06d509f7b"
 OrdinaryDiffEqLowStorageRK = "b0944070-b475-4768-8dec-fb6eb410534d"
 Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 TrixiBottomTopography = "86af9953-43df-404b-8eaa-d20d82623a82"
@@ -14,6 +15,7 @@ path = ".."
 CairoMakie = "0.13, 0.15"
 Changelog = "1.1"
 Documenter = "1.3"
+KernelInterpolation = "0.3.5"
 OrdinaryDiffEqLowStorageRK = "1.3"
 Trixi = "0.12, 0.13"
 TrixiShallowWater = "0.2"
diff --git a/docs/make.jl b/docs/make.jl
@@ -2,6 +2,9 @@ using Documenter
 import Pkg
 using TrixiBottomTopography
 using Changelog: Changelog
+# To load the extensions
+import CairoMakie
+import KernelInterpolation
 
 # Copy list of authors to not need to synchronize it manually.
 # Since the authors header exists twice we create a unique identifier for the docs section.
@@ -89,7 +92,9 @@ rm(joinpath(@__DIR__, "src", "changelog_tmp.md"))
 
 # Make documentation
 makedocs(;
-         modules = [TrixiBottomTopography],
+         modules = [TrixiBottomTopography,
+             Base.get_extension(TrixiBottomTopography, :MakieExt),
+             Base.get_extension(TrixiBottomTopography, :KernelInterpolationExt)],
          authors = "Andrew R. Winters <andrew.ross.winters@liu.se>, Michael Schlottke-Lakemper <michael@sloede.com>",
          sitename = "TrixiBottomTopography.jl",
          format = Documenter.HTML(;
@@ -110,12 +115,12 @@ makedocs(;
              "Advanced topics & developers" => ["Development" => "development.md",
                  "Style guide" => "styleguide.md",
                  "Testing" => "testing.md"],
+             "Reference" => "reference.md",
              "Changelog" => "changelog.md",
              "Authors" => "authors.md",
              "Contributing" => "contributing.md",
              "Code of Conduct" => "code_of_conduct.md",
-             "License" => "license.md",
-             "Reference" => "reference.md"])
+             "License" => "license.md"])
 
 # Documenter can also automatically deploy documentation to gh-pages.
 # See "Hosting Documentation" and deploydocs() in the Documenter manual
diff --git a/docs/src/reference.md b/docs/src/reference.md
@@ -7,3 +7,19 @@ CurrentModule = TrixiBottomTopography
 ```@autodocs
 Modules = [TrixiBottomTopography]
 ```
+
+## Makie.jl extension
+
+```@docs
+TrixiBottomTopography.evaluate_two_dimensional_interpolant
+TrixiBottomTopography.plot_topography
+TrixiBottomTopography.plot_topography_with_interpolation_knots
+```
+
+## KernelInterpolation.jl extension
+
+```@docs
+TrixiBottomTopography.RBFInterpolation
+TrixiBottomTopography.RBFInterpolation1D
+TrixiBottomTopography.RBFInterpolation2D
+```
diff --git a/examples/rhine_data_rbf_1d.jl b/examples/rhine_data_rbf_1d.jl
@@ -0,0 +1,38 @@
+###############################################################################
+# Script which uses the functionalities implemented in TrixiBottomTopography  #
+# to plot a one dimensional RBF interpolated section of the Rhine river with  #
+# data extracted along the x-direction.                                       #
+###############################################################################
+
+# Include packages
+using TrixiBottomTopography
+using KernelInterpolation
+
+# Define data path
+root_dir = pkgdir(TrixiBottomTopography)
+data = joinpath(root_dir, "examples", "data", "rhine_data_1d_20_x.txt")
+
+# Define RBF structure
+kernel = ThinPlateSplineKernel{1}()
+itp = RBFInterpolation1D(data, kernel)
+# Define interpolation function
+itp_func(x) = itp(x)
+
+# Evaluate the thin plate spline on a new set of nodes and plot
+if isdefined(Main, :Makie)
+    nodes = nodeset(itp)
+    x_knots = sort(unique(nodes[:, 1]))
+    y_knots = itp_func.(x_knots)
+
+    # Define interpolation points
+    n = 200
+    x_int_pts = Vector(LinRange(x_knots[1], x_knots[end], n))
+
+    # Get interpolated values
+    y_int_pts = itp_func.(x_int_pts)
+
+    plot_topography_with_interpolation_knots(x_int_pts, y_int_pts, x_knots, y_knots;
+                                             xlabel = "ETRS89 East",
+                                             ylabel = "DHHN2016 Height",
+                                             legend_position = :rb)
+end
diff --git a/examples/rhine_data_rbf_2d.jl b/examples/rhine_data_rbf_2d.jl
@@ -0,0 +1,39 @@
+##############################################################################
+# Script which uses the functionalities implemented in TrixiBottomTopography #
+# to plot a RBF interpolated section of the Rhine river.                     #
+##############################################################################
+
+# Include packages
+using TrixiBottomTopography
+using KernelInterpolation
+
+# Define data path
+root_dir = pkgdir(TrixiBottomTopography)
+data = joinpath(root_dir, "examples", "data", "rhine_data_2d_20.txt")
+
+# Define RBF structure
+kernel = ThinPlateSplineKernel{2}()
+itp = RBFInterpolation2D(data, kernel)
+# Define interpolation function
+itp_func(x, y) = itp([x; y])
+
+# Evaluate the thin plate spline on a new set of nodes and plot
+if isdefined(Main, :Makie)
+    nodes = nodeset(itp)
+    x_knots = sort(unique(nodes[:, 1]))
+    y_knots = sort(unique(nodes[:, 2]))
+    z_knots = evaluate_two_dimensional_interpolant(itp_func, x_knots, y_knots)
+
+    # Define interpolation points
+    n = 100
+    x_int_pts = Vector(LinRange(x_knots[1], x_knots[end], n))
+    y_int_pts = Vector(LinRange(y_knots[1], y_knots[end], n))
+
+    # Get interpolated matrix
+    z_int_pts = evaluate_two_dimensional_interpolant(itp_func, x_int_pts, y_int_pts)
+
+    plot_topography_with_interpolation_knots(x_int_pts, y_int_pts, z_int_pts,
+                                             x_knots, y_knots, z_knots;
+                                             xlabel = "ETRS89 East",
+                                             ylabel = "DHHN2016 Height", zlabel = "H")
+end
diff --git a/ext/KernelInterpolationExt.jl b/ext/KernelInterpolationExt.jl
@@ -0,0 +1,66 @@
+# Package extension for adding KernelInterpolation.jl-based features to TrixiBottomTopography.jl
+module KernelInterpolationExt
+
+using KernelInterpolation: NodeSet, interpolate
+using TrixiBottomTopography
+
+"""
+    RBFInterpolation(nodeset::NodeSet, z, args...; kwargs...)
+
+Create a radial basis function (RBF) interpolant using [KernelInterpolation.jl](https://github.com/JoshuaLampert/KernelInterpolation.jl)
+from a set of scattered nodes and corresponding elevation values `z`.
+The arguments `args` and keyword arguments `kwargs` are passed to the `interpolate` function of KernelInterpolation.jl.
+See the [KernelInterpolation.jl documentation](https://joshualampert.github.io/KernelInterpolation.jl/stable/) for more details.
+"""
+function TrixiBottomTopography.RBFInterpolation(nodeset::NodeSet, z, args...; kwargs...)
+    return interpolate(nodeset, z, args...; kwargs...)
+end
+
+"""
+    RBFInterpolation1D(x, y, args...; kwargs...)
+    RBFInterpolation1D(path, args...; kwargs...)
+
+Create a 1D radial basis function (RBF) interpolant using [KernelInterpolation.jl](https://github.com/JoshuaLampert/KernelInterpolation.jl).
+The arguments `args` and keyword arguments `kwargs` are passed to the `interpolate` function of KernelInterpolation.jl.
+See the [KernelInterpolation.jl documentation](https://joshualampert.github.io/KernelInterpolation.jl/stable/) for more details.
+
+The RBF interpolant can be created with [`RBFInterpolation`](@ref) from a set of scattered nodes and corresponding elevation values `y` using
+a `NodeSet` from KernelInterpolation.jl. Alternatively, the same interface as for [`LinearBSpline`](@ref) and
+[`CubicBSpline`](@ref) can be used with `RBFInterpolation1D`, i.e., either by providing vectors `x` of coordinates and `y`
+of elevation values or by providing a path to a text file containing the data `x` and `y`.
+"""
+function TrixiBottomTopography.RBFInterpolation1D(x, y, args...; kwargs...)
+    nodeset = NodeSet(vec([[xx] for xx in x]))
+    return TrixiBottomTopography.RBFInterpolation(nodeset, vec(y), args...; kwargs...)
+end
+
+function TrixiBottomTopography.RBFInterpolation1D(path::String, args...; kwargs...)
+    x, y = TrixiBottomTopography.parse_txt_1D(path)
+    return TrixiBottomTopography.RBFInterpolation1D(x, y, args...; kwargs...)
+end
+
+"""
+    RBFInterpolation2D(x, y, z, args...; kwargs...)
+    RBFInterpolation2D(path, args...; kwargs...)
+
+Create a 2D radial basis function (RBF) interpolant using [KernelInterpolation.jl](https://github.com/JoshuaLampert/KernelInterpolation.jl).
+The arguments `args` and keyword arguments `kwargs` are passed to the `interpolate` function of KernelInterpolation.jl.
+See the [KernelInterpolation.jl documentation](https://joshualampert.github.io/KernelInterpolation.jl/stable/) for more details.
+
+The RBF interpolant can be created with [`RBFInterpolation`](@ref) from a set of scattered nodes and corresponding elevation values `z` using
+a `NodeSet` from KernelInterpolation.jl. Alternatively, the same interface as for [`BilinearBSpline`](@ref) and
+[`BicubicBSpline`](@ref) can be used with `RBFInterpolation2D`, i.e., either by providing vectors `x` and `y` of coordinates and a matrix `z`
+of elevation values or by providing a path to a text file containing the data `x`, `y`, and `z`.
+"""
+function TrixiBottomTopography.RBFInterpolation2D(x, y, z, args...; kwargs...)
+    nodeset = NodeSet(vec([[xx; yy] for xx in x, yy in y]))
+    return TrixiBottomTopography.RBFInterpolation(nodeset, vec(Matrix(z')), args...;
+                                                  kwargs...)
+end
+
+function TrixiBottomTopography.RBFInterpolation2D(path::String, args...; kwargs...)
+    x, y, z = TrixiBottomTopography.parse_txt_2D(path)
+    return TrixiBottomTopography.RBFInterpolation2D(x, y, z, args...; kwargs...)
+end
+
+end
diff --git a/ext/MakieExt.jl b/ext/MakieExt.jl
@@ -1,5 +1,5 @@
 # Package extension for adding Makie-based features to TrixiBottomTopography.jl
-module TrixiBottomTopographyMakieExt
+module MakieExt
 
 # Required for visualization code
 using Makie: Makie
diff --git a/src/1D/spline_cache_1D.jl b/src/1D/spline_cache_1D.jl
@@ -103,13 +103,7 @@ Note that the number of `x` and `y` values have to be the same.
 An example can be found [here](https://gist.githubusercontent.com/maxbertrand1996/b05a90e66025ee1ebddf444a32c3fa01/raw/90d375c1ac11b26589aab1fe92bd0e6f6daf37b7/Rhine_data_1D_10.txt)
 """
 function LinearBSpline(path::String)
-    file = open(path)
-    lines = readlines(file)
-    close(file)
-
-    num_elements = parse(Int64, lines[2])
-    x = [parse(Float64, val) for val in lines[4:(3 + num_elements)]]
-    y = [parse(Float64, val) for val in lines[(5 + num_elements):end]]
+    x, y = parse_txt_1D(path)
 
     LinearBSpline(x, y)
 end
@@ -326,14 +320,8 @@ The .txt file has to have the following structure to be interpreted by this func
 Note that the number of `x` and `y` values have to be the same.
 An example can be found [here](https://gist.githubusercontent.com/maxbertrand1996/b05a90e66025ee1ebddf444a32c3fa01/raw/90d375c1ac11b26589aab1fe92bd0e6f6daf37b7/Rhine_data_1D_10.txt)
 """
-function CubicBSpline(path::String; end_condition = "free", smoothing_factor = 0.0)
-    file = open(path)
-    lines = readlines(file)
-    close(file)
-
-    num_elements = parse(Int64, lines[2])
-    x = [parse(Float64, val) for val in lines[4:(3 + num_elements)]]
-    y = [parse(Float64, val) for val in lines[(5 + num_elements):end]]
+function CubicBSpline(path::String; kwargs...)
+    x, y = parse_txt_1D(path)
 
-    CubicBSpline(x, y; end_condition = end_condition, smoothing_factor = smoothing_factor)
+    CubicBSpline(x, y; kwargs...)
 end
diff --git a/src/1D/spline_utils_1D.jl b/src/1D/spline_utils_1D.jl
@@ -91,3 +91,15 @@ function spline_smoothing(lambda::Number, Delta::Number, y::Vector)
 
     return inv(diagm(ones(n)) + lambda * K) * y
 end
+
+function parse_txt_1D(path)
+    file = open(path)
+    lines = readlines(file)
+    close(file)
+
+    num_elements = parse(Int64, lines[2])
+    x = [parse(Float64, val) for val in lines[4:(3 + num_elements)]]
+    y = [parse(Float64, val) for val in lines[(5 + num_elements):end]]
+
+    return x, y
+end
diff --git a/src/2D/spline_cache_2D.jl b/src/2D/spline_cache_2D.jl
@@ -130,20 +130,9 @@ The .txt file has to have the following structure to be interpreted by this func
 An example can be found [here](https://gist.githubusercontent.com/maxbertrand1996/7b1b943eac142d5bc836bb818fe83a5a/raw/74228e349e91fbfe1563479f99943b469f26ac62/Rhine_data_2D_10.txt)
 """
 function BilinearBSpline(path::String)
-    file = open(path)
-    lines = readlines(file)
-    close(file)
+    x, y, z = parse_txt_2D(path)
 
-    n = parse(Int64, lines[2])
-    m = parse(Int64, lines[4])
-
-    x = [parse(Float64, val) for val in lines[6:(5 + n)]]
-    y = [parse(Float64, val) for val in lines[(7 + n):(6 + n + m)]]
-    z_tmp = [parse(Float64, val) for val in lines[(8 + n + m):end]]
-
-    z = transpose(reshape(z_tmp, (n, m)))
-
-    BilinearBSpline(x, y, Matrix(z))
+    BilinearBSpline(x, y, z)
 end
 
 ######################################
@@ -566,20 +555,8 @@ The .txt file has to have the following structure to be interpreted by this func
 
 An example can be found [here](https://gist.githubusercontent.com/maxbertrand1996/7b1b943eac142d5bc836bb818fe83a5a/raw/74228e349e91fbfe1563479f99943b469f26ac62/Rhine_data_2D_10.txt)
 """
-function BicubicBSpline(path::String; end_condition = "free", smoothing_factor = 0.0)
-    file = open(path)
-    lines = readlines(file)
-    close(file)
-
-    n = parse(Int64, lines[2])
-    m = parse(Int64, lines[4])
-
-    x = [parse(Float64, val) for val in lines[6:(5 + n)]]
-    y = [parse(Float64, val) for val in lines[(7 + n):(6 + n + m)]]
-    z_tmp = [parse(Float64, val) for val in lines[(8 + n + m):end]]
-
-    z = transpose(reshape(z_tmp, (n, m)))
+function BicubicBSpline(path::String; kwargs...)
+    x, y, z = parse_txt_2D(path)
 
-    BicubicBSpline(x, y, Matrix(z); end_condition = end_condition,
-                   smoothing_factor = smoothing_factor)
+    BicubicBSpline(x, y, z; kwargs...)
 end
diff --git a/src/2D/spline_utils_2D.jl b/src/2D/spline_utils_2D.jl
diff --git a/src/TrixiBottomTopography.jl b/src/TrixiBottomTopography.jl
diff --git a/test/Project.toml b/test/Project.toml
diff --git a/test/runtests.jl b/test/runtests.jl
diff --git a/test/test_rbf_1d.jl b/test/test_rbf_1d.jl
diff --git a/test/test_rbf_2d.jl b/test/test_rbf_2d.jl
