GitHub

B-spline interpolation structure

Once the underlying data is in the correct format, we can start defining B-spline structures which we will use later to define the interpolation functions.

One dimensional structures

For the one dimensional case, the structures LinearBSpline and CubicBSpline are available. They contain all relevant values to define linear and cubic B-spline interpolation functions corresponding to linear and cubic B-spline interpolation. These are:

  • x: A vector of values in x-direction.
  • Delta: The length of a single patch in the given data set. A patch is the area between two consecutive x values. Delta corresponds to the distance between two consecutive values in x-direction. As we are only considering Cartesian grids, Delta is equal for all patches.
  • Q: A vector that contains the control points.
  • IP: The coefficients matrix.

To populate the structure, the outer constructor functions LinearBSpline(data_path) and CubicBSpline(data_path) are implemented. These constructors use the files in data_path to obtain the values which will be stored in the corresponding structure, as well as LinearBSpline(x,y) and CubicBSpline(x,y) which use given vectors x and y.

To get a better idea of the constructor functions, we consider the example rhine_data_cubic-nak.jl from the examples folder of this repo. This particular example reads one dimensional bottom topography data from a .txt file and constructs a cubic B-spline interpolation with not-a-knot end condition and smoothing of the data.

# Include packages
using TrixiBottomTopography

# Define data path
root_dir = pkgdir(TrixiBottomTopography)
data = joinpath(root_dir, "examples", "data", "rhine_data_1d_20_x.txt")

# Define B-spline structure
spline_struct = CubicBSpline(data; end_condition = "not-a-knot", smoothing_factor = 999)
CubicBSpline{Vector{Float64}, Float64, Vector{Float64}, Matrix{Int64}}([357000.0, 357020.0, 357040.0, 357060.0, 357080.0, 357100.0, 357120.0, 357140.0, 357160.0, 357180.0  …  357800.0, 357820.0, 357840.0, 357860.0, 357880.0, 357900.0, 357920.0, 357940.0, 357960.0, 357980.0], 20.0, [46.91426684256018, 47.397904854124576, 46.99689666444544, 46.50889411070535, 46.73154903008692, 46.862924837363984, 46.74903591362994, 46.40702226591983, 46.31904390610859, 46.654018220049146  …  39.07284424693978, 39.1231081223708, 39.13738000648763, 39.136760208084766, 39.10325245367299, 39.08822484618267, 39.09410545889054, 39.1030246596446, 39.08484073202447, 39.009411959609665], [-1 3 -3 1; 3 -6 3 0; -3 0 3 0; 1 4 1 0])

For the cubic case, we can also set the optional parameters end_condition, which defines (as the name suggests) the end condition of the spline. Available end conditions are the not-a-knot and the free end condition. By default, end_condition is set to free. If you are not familiar with the differences and influence of these end conditions, see Chapter 1 of

  • Quentin Agrapart & Alain Batailly (2020), Cubic and bicubic spline interpolation in Python. hal-03017566v2.

Besides the end condition, we can also specify a smoothing_factor for the cubic B-spline. This smoothing parameter defines a trade-off for the resulting cubic B-spline interpolation between how well the B-spline models (or fits) the original data and minimizing the curvature by defining new y values. This procedure is called spline smoothing. There is no general approach to determine which smoothing_factor is best suited for a given problem. It must be determined by the user via trial and error. However, as a general rule, the larger the smoothing factor, the less curvature will be present in the resulting B-spline. To understand the underlying maths of the smoothing spline procedure, please see:

Two dimensional structures

For the two dimensional case, the structures BilinearBSpline and BicubicBSpline are implemented, which contain all relevant values to define bilinear and bicubic B-spline interpolation functions. These are:

  • x: Vector of values in x-direction.
  • y: Vector of values in y-direction.
  • Delta: Length of one side of a single patch in the given data set. A patch is the area between two consecutive x and y values. Delta corresponds to the distance between two consecutive values in x-direction. The implementation only considers Cartesian grids, so Delta is equal for all patches in x and y-direction.
  • Q: Matrix which contains the control points.
  • IP: Coefficients matrix.

To populate the structure, the outer constructor functions BilinearBSpline(data_path) and BicubicBSpline(data_path) are implemented. They use the files in data_path to obtain the values which will be stored in the corresponding structure. Further, BilinearBSpline(x,y,z) and BicubicBSpline(x,y,z) which use given vectors x and y and matrix z. The x, y and z data are organized in the following form:

\[\begin{aligned} \begin{matrix} & & x_1 & x_2 & ... & x_n\\ & & & & &\\ y_1 & & z_{11} & z_{12} & ... & z_{1n}\\ y_1 & & z_{21} & z_{22} & ... & z_{2n}\\ \vdots & & \vdots & \vdots & \ddots & \vdots\\ y_m & & z_{m1} & z_{m2} & ... & z_{mn} \end{matrix} \end{aligned}\]

To better understand the constructor functions, we consider the example rhine_data_bicubic-nak.jl of this repo. This example reads in two dimensional bottom topography data from a .txt file and creates a bicubic B-spline interpolation with not-a-knot end condition and smoothing of the data.

# Include packages
using TrixiBottomTopography

# Define data path
root_dir = pkgdir(TrixiBottomTopography)
data = joinpath(root_dir, "examples", "data", "rhine_data_2d_20.txt")

# Define B-spline structure
spline_struct = BicubicBSpline(data; end_condition = "not-a-knot", smoothing_factor = 9999)
BicubicBSpline{Vector{Float64}, Vector{Float64}, Float64, Matrix{Float64}, Matrix{Int64}}([357000.0, 357020.0, 357040.0, 357060.0, 357080.0, 357100.0, 357120.0, 357140.0, 357160.0, 357180.0  …  357800.0, 357820.0, 357840.0, 357860.0, 357880.0, 357900.0, 357920.0, 357940.0, 357960.0, 357980.0], [5.646019e6, 5.646039e6, 5.646059e6, 5.646079e6, 5.646099e6, 5.646119e6, 5.646139e6, 5.646159e6, 5.646179e6, 5.646199e6  …  5.646819e6, 5.646839e6, 5.646859e6, 5.646879e6, 5.646899e6, 5.646919e6, 5.646939e6, 5.646959e6, 5.646979e6, 5.646999e6], 20.0, [46.600451320002854 47.282399695552705 … 47.09971997322787 47.10735082392629; 45.82016087061047 46.502109246160636 … 46.963741530201055 46.97137238089921; … ; 45.99186284178017 45.686484464847375 … 39.098457400873656 39.23988800954714; 46.243985904572256 45.938607527639306 … 39.06537978820518 39.206810396878865], [-1 3 -3 1; 3 -6 3 0; -3 0 3 0; 1 4 1 0])

For the bicubic spline, we can set the optional parameters end_condition, which defines (as in the one dimensional case) the end condition of the spline. Again, the available end conditions are the not-a-knot and the free end condition. By default, end_condition is set to free. If you are not familiar with the differences and influence of these end conditions, see Chapter 2 of

  • Quentin Agrapart & Alain Batailly (2020), Cubic and bicubic spline interpolation in Python. hal-03017566v2

Besides the end condition, we can also specify a smoothing_factor for the bicubic spline. Just as in the one dimensional constructor, this smoothing provides a trade-off of the bicubic B-spline interpolation between how well it models (or fits) the original topography data and minimizing the curvature of the spline by defining new z values. This procedure is called thin plate spline. There is no general approach to which smoothing_factor is best suited for the problem, and it must be determined via trial and error by the user. To understand the underlying maths of the thin plate spline, please see: