TrixiBottomTopography.jl API

TrixiBottomTopography

TrixiBottomTopography.jl is a supporting framework for Trixi.jl and TrixiShallowWater.jl, which can be used to approximate bottom topography functions using B-splines from real world data.

BicubicBSpline(x, y, Delta, Q, IP)

Two dimensional cubic B-spline structure that contains all important attributes to define a B-Spline interpolation function. These attributes 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. The value 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 in x and y-direction
• Q: Matrix which contains the control points
• IP: Coefficients matrix

BicubicBSpline(path::String; end_condition = "free", smoothing_factor = 0.0)

A function which reads in the x, y and z values for BilinearBSpline from a .txt file. The input values are:

• path: String of a path of the specific .txt file
• end_condition: String which can either be "free" or "not-a-knot" and defines which end condition should be considered. By default this is set to "free"
• smoothing_factor: Float64 $\geq$ 0.0 which specifies the degree of smoothing of the y values. By default this value is set to 0.0 which corresponds to no smoothing.

The .txt file has to have the following structure to be interpreted by this function:

• First line: comment # Number of x values
• Second line: integer which gives the number of x values
• Third line: comment # Number of y values
• Fourth line: integer which gives the number of y values
• Fifth line: comment # x values
• Following lines: the x values where each value has its own line
• Line after the x-values: comment # y values
• Following lines: y values where each value has its own line
• Line after the y-values: comment # z values
• Remaining lines: values for z where each value has its own line and is in th following order: z11, z12, ... z1n, z21, ... z2n, ..., zm1, ..., z_mn

An example can be found here

BicubicBSpline(x::Vector, y::Vector, z::Matrix; end_condition = "free", smoothing_factor = 0.0)

This function calculates the inputs for the structure BicubicBSpline. The input values are:

• x: Vector that contains equally spaced values in x-direction
• y: Vector that contains equally spaced values in y-direction where the spacing between the y-values has to be the same as the spacing between the x-values
• z: Matrix that contains the corresponding values in z-direction. Where the values are ordered in the following way:

\[\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}\]

• end_condition: a string which can either be "free" or "not-a-knot" and defines which end condition should be considered. By default this is set to "free".
• smoothing_factor: a Float64 $\geq$ 0.0 which specifies the degree of smoothing of the z values. By default this value is set to 0.0 which corresponds to no smoothing.

Bicubic B-spline interpolation is only possible if the dimensions of vectors x and y correspond with the dimensions of the matrix z.

First, the data is sorted via sort_data to guarantee that the x and y values are in ascending order with corresponding matrix z.

The patch size Delta is calculated by subtracting the second by the first x value. This can be done because we only consider equal space between consecutive x and y values. A patch is the area between two consecutive x and y values.

If a smoothing_factor $>$ 0.0 is set, the function calc_tps calculates new values for z which guarantee a resulting parametric B-spline surface with less curvature.

The coefficients matrix IP for bicubic B-splines is fixed to be

\[\begin{aligned} \begin{pmatrix} -1 & 3 & -3 & 1\\ 3 & -6 & 3 & 0\\ -3 & 0 & 3 & 0\\ 1 & 4 & 1 & 0 \end{pmatrix} \end{aligned}\]

To get the matrix of control points Q which is necessary to set up an interpolation function, we need to define a matrix Phi which maps the control points to a vector P. This can be done by solving the following system of linear equations for Q.

\[\underbrace{ \begin{bmatrix} z_{1,1} \\ z_{1,2} \\ \vdots \\ z_{1,n} \\ z_{2,1} \\ \vdots \\ z_{m,n} \\ 0 \\ \vdots \\ 0 \end{bmatrix} }_{\text{:=P} \in \mathbb{R}^{(m+2)(n+2)\times 1}} = \frac{1}{36} \Phi \cdot \underbrace{\begin{bmatrix} Q_{1,1} \\ Q_{1,2} \\ \vdots \\ Q_{1,n+2} \\ Q_{2,1} \\ \vdots \\ Q_{m+2,n+2} \end{bmatrix}}_{\text{:= Q} \in \mathbb{R}^{(m+2) \times (n+2)}}\]

For the first n $\cdot$ m lines, the matrix Phi is the same for the "free" end and the "not-a-knot" end condition. These lines have to address the following condition:

\[\begin{align*} z_{j,i} = \frac{1}{36} \Big( &Q_{j,i} + 4Q_{j+1,i} + Q_{j+2,i} + 4Q_{j,i+1} + 16Q_{j+1,i+1}\\ &+ 4Q_{j+2,i+1} + Q_{j,i+2} + 4Q_{j+1,i+2} + Q_{j+2,i+2} \Big) \end{align*}\]

for i = 1,...,n and j = 1,...,m.

The "free" end condition needs at least two values for the x and y vectors. The "free" end condition has the following additional requirements for the control points which have to be addressed by Phi:

• $Q_{j,1} - 2Q_{j,2} + Q_{j,3} = 0$ for j = 2,...,m+1
• $Q_{j,n} - 2Q_{j,n+1} + Q_{j,n+2} = 0$ for j = 2,...,m+1
• $Q_{1,i} - 2Q_{2,i} + Q_{3,i} = 0$ for i = 2,...,n+1
• $Q_{m,i} - 2Q_{m+1,i} + Q_{m+2,i} = 0$ for i = 2,...,n+1
• $Q_{1,1} - 2Q_{2,2} + Q_{3,3} = 0$
• $Q_{m+2,1} - 2Q_{m+1,2} + Q_{m,3} = 0$
• $Q_{1,n+2} - 2Q_{2,n+1} + Q_{3,n} = 0$
• $Q_{m,n} - 2Q_{m+1,n+1} + Q_{m+2,n+2} = 0$

The "not-a-knot" end condition needs at least four values for the x and y vectors.

• Continuity of the third x derivative between the leftmost and second leftmost patch
• Continuity of the third x derivative between the rightmost and second rightmost patch
• Continuity of the third y derivative between the patch at the top and the patch below
• Continuity of the third y derivative between the patch at the bottom and the patch above
• $Q_{1,1} - Q_{1,2} - Q_{2,1} + Q_{2,2} = 0$
• $Q_{m-1,1} + Q_{m,1} + Q_{m-1,2} - Q_{m,2} = 0$
• $Q_{1,n-1} + Q_{2,n} + Q_{1,n-1} - Q_{2,n} = 0$
• $Q_{m-1,n-1} - Q_{m,n-1} - Q_{m-1,n} + Q_{m,n} = 0$

A reference for the calculations in this script can be found in Chapter 2 of

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

BilinearBSpline(x, y, Delta, Q, IP)

Two dimensional bilinear B-spline structure that contains all important attributes to define a B-Spline interpolation function. These attributes 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. The value 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 in x and y-direction
• Q: Matrix which contains the control points
• IP: Coefficients matrix

BilinearBSpline(path::String)

A function which reads in the x, y and z values for BilinearBSpline from a .txt file. The input values are:

• path: String of a path of the specific .txt file

The .txt file has to have the following structure to be interpreted by this function:

• First line: comment # Number of x values
• Second line: integer which gives the number of x values
• Third line: comment # Number of y values
• Fourth line: integer which gives the number of y values
• Fifth line: comment # x values
• Following lines: the x values where each value has its own line
• Line after the x-values: comment # y values
• Following lines: y values where each value has its own line
• Line after the y-values: comment # z values
• Remaining lines: values for z where each value has its own line and is in the following order: z11, z12, ... z1n, z21, ... z2n, ..., zm1, ..., z_mn

An example can be found here

BilinearBSpline(x::Vector, y::Vector, z::Matrix)

This function calculates the inputs for the structure BilinearBSpline. The input values are:

• x: Vector that contains equally spaced values in x-direction
• y: Vector that contains equally spaced values in y-direction where the spacing between the y-values has to be the same as the spacing between the x-values
• z: Matrix that contains the corresponding values in z-direction. Where the values are ordered in the following way:

\[\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}\]

Bilinear B-spline interpolation is only possible if we have at least two values in x and two values in y and the dimensions of vectors x and y correspond with the dimensions of the matrix z.

First of all the data is sorted which is done by sort_data to guarantee that the x and y values are in ascending order with corresponding matrix z.

The patch size Delta is calculated by subtracting the second by the first x value. This can be done because we only consider equal space between consecutive x and y values. A patch is the area between two consecutive x and y values.

For bilinear B-spline interpolation, the control points Q correspond with the z values.

The coefficients matrix IP for bilinear B-splines is fixed to be

\[\begin{aligned} \begin{pmatrix} -1 & 1\\ 1 & 0 \end{pmatrix} \end{aligned}\]

A reference for the calculations in this script can be found in Chapter 2 of

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

CubicBSpline(x, Delta, Q, IP)

One dimensional cubic B-spline structure that contains all important attributes to define a B-Spline interpolation function. Similar to LinearBSpline These attributes are:

• x: Vector of values in x-direction
• Delta: Length of a single patch in the given data set. A patch is the area between two consecutive x values. The value 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: Vector which contains the Control points
• IP: Coefficients matrix

CubicBSpline(path::String; end_condition = "free", smoothing_factor = 0.0)

A function that reads in the x and y values for CubicBSpline from a .txt file. The input values are:

• path: String of a path of the specific .txt file
• end_condition: String which can either be free or not-a-knot and defines which end condition should be considered. By default this is set to free
• smoothing_factor: Float64 $\geq$ 0.0 which specifies the degree of smoothing of the y values. By default this value is set to 0.0 which corresponds to no smoothing.

The .txt file has to have the following structure to be interpreted by this function:

• First line: comment # Number of x values
• Second line: integer which gives the number of x values
• Third line: comment # x values
• Following lines: the x values where each value has its own line
• Line after the x-values: comment # y values
• Remaining lines: y values where each value has its own line

Note that the number of x and y values have to be the same. An example can be found here

CubicBSpline(x::Vector, y::Vector; end_condition = "free", smoothing_factor = 0.0)

This function calculates the inputs for the structure CubicBSpline. The input values are:

• x: Vector that contains equally spaces values in x-direction
• y: Vector that contains values in y-direction
• end_condition: String that can either be free or not-a-knot and defines which end condition should be considered. By default this is set to "free".
• smoothing_factor: Float64 $\geq$ 0.0 which specifies the degree of smoothing of the y values. By default this value is set to 0.0 that corresponds to no smoothing.

First the data is sorted via sort_data to guarantee that the x values are in ascending order.

The patch size Delta is calculated by subtracting the second and first x value. This can be done because we only consider equally spaced x values. (A patch is the area between two consecutive x values)

If a smoothing_factor > 0.0 is set, the function spline_smoothing calculates new y values which guarantee a B-Spline with less curvature.

The coefficients matrix IP for linear B-splines is fixed to be

\[\begin{aligned} IP = \begin{pmatrix} -1 & 3 & -3 & 1\\ 3 & -6 & 3 & 0\\ -3 & 0 & 3 & 0\\ 1 & 4 & 1 & 0 \end{pmatrix} \end{aligned}\]

The "free" end condition requires the second and the second to last control points lie between the first and the third control point and that the second to last control points are between the third to last and the last control point. This procedure is only possible with at least two values in x data. The system of linear equations to determine the control points have the following form:

\[\begin{aligned} \underbrace{\begin{bmatrix} 0 \\ P_1 \\ P_2 \\ \vdots \\ P_{n-1} \\ P_n\\ 0 \end{bmatrix}}_{:= P^*_{\text{free}}} = \frac{1}{6} \underbrace{ \begin{bmatrix} 1 & -2 & 1 & 0 & ... & ... & 0 \\ 1 & 4 & 1 & 0 & ... & ... & 0\\ 0 & 1 & 4 & 1 & 0 & & \vdots\\ \vdots & 0 & \ddots & \ddots & \ddots & 0 & \vdots\\ \vdots & & 0 & 1 & 4 & 1 & 0\\ 0 & ... & ... & 0 & 1 & 4 & 1\\ 0 & ... & ... & 0 & 1 & -2 & 1 \end{bmatrix} }_{:= \Phi^*_{\text{free}}} \underbrace{\begin{bmatrix} Q_1 \\ Q_2 \\ Q_3 \\ \vdots \\ Q_n \\ Q_{n+1} \\ Q_{n+2} \end{bmatrix}}_{:= Q_{\text{free}}}, \end{aligned}\]

which is solved for $Q_{\text{free}}$.

The "not-a-knot" end condition requires the continuity of the third derivative in the second and second to last fit knot. This end condition is only possible with at least four values in x data. The system of linear equations to determine the control points has the following form:

\[\begin{aligned} \underbrace{\begin{bmatrix} 0 \\ P_1 \\ P_2 \\ \vdots \\ P_{n-1} \\ P_n\\ 0 \end{bmatrix}}_{:= P^*_{\text{not-a-knot}}} = \frac{1}{6} \underbrace{ \begin{bmatrix} -1 & 4 & -6 & 4 & -1 & 0 &... & 0 \\ 1 & 4 & 1 & 0 & ... & ... & ... & 0\\ 0 & 1 & 4 & 1 & 0 & & &\vdots\\ \vdots & 0 & \ddots & \ddots & \ddots & 0 & &\vdots\\ \vdots & & 0 & \ddots & \ddots & \ddots & 0 &\vdots\\ \vdots & & & 0 & 1 & 4 & 1 & 0\\ 0 & ... & ... & ... & 0 & 1 & 4 & 1\\ 0 & ... & 0 & -1 & 4 & -6 & 4 & -1 \end{bmatrix} }_{:= \Phi^*_{\text{not-a-knot}}} \underbrace{\begin{bmatrix} Q_1 \\ Q_2 \\ Q_3 \\ \vdots \\ \vdots \\ Q_n \\ Q_{n+1} \\ Q_{n+2} \end{bmatrix}}_{:= Q_{\text{not-a-knot}}}. \end{aligned}\]

which is solved for $Q_{\text{not-a-knot}}$.

For both cases $P_1,...,P_n = y_1,...,y_n$.

A reference for the calculations in this script can be found in Chapter 1 of

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

LinearBSpline(x, Delta, Q, IP)

One dimensional linear B-spline structure which contains all important attributes to define a B-Spline interpolation function. These attributes are:

• x: Vector of values in x-direction
• Delta: Length of a single patch in the given data set. A patch is the area between two consecutive x values. The value 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: Vector which contains the control points
• IP: Coefficients matrix

LinearBSpline(path::String)

A function that reads in the x and y values for LinearBSpline from a .txt file. The input values are:

• path: String of a path of the specific .txt file

The .txt file has to have the following structure to be interpreted by this function:

• First line: comment # Number of x values
• Second line: integer which gives the number of x values
• Third line: comment # x values
• Following lines: the x values where each value has its own line
• Line after the x-values: comment # y values
• Remaining lines: y values where each value has its own line

Note that the number of x and y values have to be the same. An example can be found here

LinearBSpline(x::Vector, y::Vector)

This function calculates the inputs for the structure LinearBSpline. The input values are:

• x: A vector that contains equally spaced values in x-direction
• y: A vector that contains values in y-direction

Linear B-spline interpolation is only possible if the data set has at least two values in x.

First the data is sorted via sort_data to guarantee that the x values are in ascending order.

The patch size Delta is calculated by subtracting the second and first x values. This can be done because we only consider equally spaced x values. A patch is the area between two consecutive x values.

For linear B-spline interpolation, the control points Q correspond with the values in y.

The coefficients matrix IP for linear B-splines is fixed to be

\[\begin{aligned} \begin{pmatrix} -1 & 1\\ 1 & 0 \end{pmatrix} \end{aligned}\]

A reference for the calculations in this script can be found in Chapter 1 of

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

calc_tps(lambda::Number, x::Vector, y::Vector, z::Matrix)

The inputs to this function are:

• lambda: Smoothing factor which specifies the degree of the smoothing that should take place
• x: Vector of x values
• y: Vector of x values
• z: Matrix with the z values to be smoothed where the values of z correspond to the indexing (y,x)

This function uses the thin plate spline approach to perform the smoothing. To do so, the following linear equations system has to be solved for coeff:

\[\begin{aligned} \underbrace{ \begin{bmatrix} K & P \\ P^T & O \end{bmatrix} }_{:= L} \underbrace{\begin{bmatrix} w \\ a \end{bmatrix}}_{\text{:= coeff}} = \underbrace{\begin{bmatrix} z\\o \end{bmatrix}}_{\text{:= rhs}} \end{aligned}\]

First, the inputs are restructured using the function restructure_data and saved in the variables x_hat, y_hat and z_hat.

Then the matrix L can be filled by setting K = tps_base_func(||(x_hat[i], y_hat[i]) - (x_hat[j], y_hat[j])||) where || || is the Euclidian norm, P = [1 x y] and O = $3\times 3$ zeros matrix.

Afterwards the vector rhs is filled by setting z = z_hat and o = a vector with three zeros.

Now the system is solved to redeem the vector coeff. This vector is then used to calculate the smoothed values for z and save them in H_f by the following function:

\[\begin{align*} z\_smth[i] = &a[1] + a[2]x\_hat[i] + a[3]y\_hat[i] \\ + &\sum_{j = 0}^p w[j] tps\_base\_func(\|(x\_hat[i], y\_hat[i]) - (x\_hat[j], y\_hat[j]) \|) \end{align*}\]

here p is the number of entries in z_hat.

A reference to the calculations can be found in the lecture notes of

• Gianluca Donato and Serge Belongie (2001), Approximate Thin Plate Spline Mappings DOI: 10.1007/3-540-47977-5_2

convert_dgm_1d(path_read::String, path_write::String; excerpt = 1, direction = "x", section = 1)

Function to convert DGM data files into one dimensional readable files.

Inputs:

• path_read: String of the path of the DGM data which should be converted
• path_write: String of the path where the new file should be saved. (Needs to also include the name of the file)
• excerpt: Optional integer that specifies a stride through of the data that will be extracted. E.g. if excerpt is set to 10, only every 10th x and y value are considered with their corresponding z values. The default value is 1 which means that every value is taken.
• direction: Optional String that specifies if the one dimensional data should be read from the x or y direction. By default this is set to the x direction.
• section: Optional integer which can be between 1 and 1000 and specifies which section of the other dimension should be chosen. By default this values is set to 1 which means that if direction is set to x, the corresponding z values are taken with respect to the first y value.

convert_dgm_2d(path_read::String, path_write::String; excerpt = 1)

Function to convert DGM data files into two dimensional readable files.

Inputs:

• path_read: String of the path of the DGM data which should be converted
• path_write: String of the path where the new file should be saved. (Needs to also include the name of the file)
• excerpt: Optional integer that specifies a stride through of the data that will be extracted. E.g. if excerpt is set to 10, only every 10th x and y value are considered with their corresponding z values. The default value is 1 which means that every value is taken.

TrixiBottomTopography.default_example()

Function which calls the example file "rhinedatabicubic-nak.jl" to check if including the package has worked. The example script creates a bicubic B-spline with the "not-a-knot" boundary condition. If a Makie.jl backend such as GLMakie.jl is included a plot will be generated as well.

  restructure_data(x::Vector, y::Vector, z::Matrix)

This function restructures the input values

• x: a vector with n values in x-direction
• y: a vector with m values in y-direction
• z: a m $\times$ n matrix with values in z-direction where the values of z correspond to the indexing (y,x)

The output is of the following form:

\[\begin{aligned} \begin{bmatrix} x_1 & y_1 & z_{1,1}\\ x_2 & y_1 & z_{1,2}\\ & \vdots & \\ x_n & y_1 & z_{1,n}\\ x_1 & y_2 & z_{2,1}\\ & \vdots & \\ x_n & y_m & z_{m,n} \end{bmatrix} \end{aligned}\]

sort_data(x::Vector, y::Vector, z::Matrix)

This function sorts the inputs vectors x and y in a ascending order and also reorders the input matrix z accordingly.

Therefore, first the x values are sorted with the matrix z accordingly and afterwards the y values are sorted with the matrix z accordingly.

The sorted x, y and z values are returned.

sort_data(x::Vector,y::Vector)

Sorts the input vectors x and y so that x is in ascending order and the y values accordingly to still correspond to the x values.

spline_interpolation(b_spline::BicubicBSpline, x::Number, y::Number)

The inputs are the BicubicBSpline object and the variable x and y at which the spline will be evaluated.

The parameters i and j indicate the patch in which (x,y) is located. This information is also used to get the correct control points from Q. A patch is the area between two consecutive b_spline.x and b_spline.y values.

my is an interim variable which maps x to the interval $[0,1]$ for further calculations. ny does the same for y.

To evaluate the spline at (x,y), we have to calculate the following:

\[\begin{aligned} c_{i,j,3}(\mu_i(x),\nu_j(y)) = \frac{1}{36} \begin{bmatrix} \nu_j^3(y) \\ \nu_j^2(y) \\ \nu_j(y) \\ 1 \end{bmatrix}^T \underbrace{\begin{bmatrix} -1 & 3 & -3 & 1\\ 3 & -6 & 3 & 0\\ -3 & 0 & 3 & 0\\ 1 & 4 & 1 & 0 \end{bmatrix}}_{\text{IP}} \begin{bmatrix} Q_{i,j} & Q_{i+1,j} & Q_{i+2,j} & Q_{i+3,j}\\ Q_{i,j+1} & Q_{i+1,j+1} & Q_{i+2,j+1} & Q_{i+3,j+1}\\ Q_{i,j+2} & Q_{i+1,j+2} & Q_{i+2,j+2} & Q_{i+3,j+2}\\ Q_{i,j+3} & Q_{i+1,j+3} & Q_{i+2,j+3} & Q_{i+3,j+3} \end{bmatrix} \underbrace{\begin{bmatrix} -1 & 3 & -3 & 1\\ 3 & -6 & 0 & 4\\ -3 & 3 & 3 & 1\\ 1 & 0 & 0 & 0 \end{bmatrix}}_{\text{IP}} \begin{bmatrix} \mu_i^3(x) \\ \mu_i^2(x) \\ \mu_i(x) \\ 1 \end{bmatrix} \end{aligned}\]

A reference for the calculations in this script can be found in Chapter 2 of

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

spline_interpolation(b_spline::BilinearBSpline, x::Number, y::Number)

The inputs are the BilinearBSpline object and the variable x and y at which the spline will be evaluated.

The parameters i and j indicate the patch in which (x,y) is located. This information is also used to get the correct control points from Q. A patch is the area between two consecutive b_spline.x and b_spline.y values.

my is an interim variable that maps x to the interval $[0,1]$ for further calculations. ny does the same for y.

To evaluate the spline at (x,y), we have to calculate the following:

\[\begin{aligned} c_{i,j,1}(\mu_i(x),\nu_j(y)) = \begin{bmatrix} \nu_j(y)\\ 1 \end{bmatrix}^T \underbrace{\begin{bmatrix} -1 & 1\\ 1 & 0 \end{bmatrix}}_{\text{IP}} \begin{bmatrix} Q_{i,j} & Q_{i,j+1}\\ Q_{i+1,j} & Q_{i+1,j+1} \end{bmatrix} \underbrace{\begin{bmatrix} -1 & 1\\ 1 & 0 \end{bmatrix}}_{\text{IP}^T} \begin{bmatrix} \mu_i(x) \\ 1\end{bmatrix} \end{aligned}\]

A reference for the calculations in this script can be found in Chapter 2 of

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

spline_interpolation(b_spline::CubicBSpline, x::Number)

The inputs are the CubicBSpline object and a variable x at which the spline will be evaluated.

The parameter i indicates the patch in which the variable x is located. This parameter is also used to get the correct control points from Q. A patch is the area between two consecutive b_spline.x values.

kappa is an interim variable which maps t to the interval $[0,1]$ for further calculations.

To evaluate the spline at x, we have to calculate the following:

\[\begin{aligned} c_{i,3}\left(\kappa_i(x) \right) = \frac{1}{6} \begin{bmatrix} \kappa_i(x)^3\\ \kappa_i(x)^2\\ \kappa_i(x) \\1 \end{bmatrix}^T \underbrace{\begin{bmatrix} -1 & 3 & -3 & 1\\ 3 & -6 & 3 & 0\\ -3 & 0 & 3 & 0\\ 1 & 4 & 1 & 0 \end{bmatrix}}_{\text{IP}} \begin{bmatrix} Q_{i,\text{free}}\\ Q_{i+1,\text{free}}\\ Q_{i+2,\text{free}}\\ Q_{i+3,\text{free}} \end{bmatrix} \end{aligned}\]

A reference for the calculations in this script can be found in Chapter 1 of

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

spline_interpolation(b_spline::LinearBSpline, x::Number)

The inputs are the LinearBSpline object and a variable x at which the spline will be evaluated.

The parameter i indicates the patch in which the variable x is located. This parameter is also used to get the correct control points from Q. A patch is the area between two consecutive b_spline.x values.

kappa is an interim variable which maps x to the interval $[0,1]$ for further calculations.

To evaluate the spline at x, we have to calculate the following:

\[\begin{aligned} c_{i,1}(\kappa_i(x)) = \begin{bmatrix} \kappa_i(x)\\ 1 \end{bmatrix}^T \begin{bmatrix} -1 & 1\\1 & 0 \end{bmatrix} \begin{bmatrix} Q_i\\Q_{i+1} \end{bmatrix} \end{aligned}\]

A reference for the calculations in this script can be found in Chapter 1 of

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

spline_smoothing(lambda::Number, Delta::Number, y::Vector)

The inputs to this function are:

• lambda: Smoothing factor which specifies the degree of the smoothing that should take place
• Delta: Step size of a patch (A patch is the area between two consecutive x values)
• y: Data values to be smoothed

The goal is to find a new interpolation values $\hat{y}$ for $y$, so that for given $\lambda$, the following equation is minimized:

\[\begin{aligned} \text{PSS} = \sum_{i = 1}^{n} \left( y_i - \underbrace{S(t_i)}_{=\hat{y}_i} \right)^2 + \lambda \int_{x_1}^{x_n} (S''(t))^2 dt, \end{aligned}\]

where $S(t)$ is a cubic spline function. $\hat{y}$ is determined as follows:

\[\begin{aligned} \hat{y} = (I+\lambda K)^{-1} y \end{aligned}\]

where $I$ is the $n \times n$ identity matrix and $K = \Delta_2^T W^{-1} \Delta_2$ with

\[\begin{aligned} \Delta_2 = \begin{pmatrix} 1/\Delta & -2/\Delta & 1/\Delta & ... & 0\\ 0 & \ddots & \ddots & \ddots & 0\\ 0 & ... & 1/\Delta & -2/\Delta & 1/\Delta \end{pmatrix} \in \mathbb{R}^{(n-2) \times n} \end{aligned}\]

and

\[\begin{aligned} W = \begin{pmatrix} 2/3 \Delta & 1/6 \Delta & 0 & ... & 0\\ 1/6 \Delta & 2/3 \Delta & 1/6 \Delta & ... & 0\\ 0 & \ddots & \ddots & \ddots & 0\\ 0 & ... & 0 & 2/3 \Delta & 1/6 \Delta \end{pmatrix} \in \mathbb{R}^{n \times n} \end{aligned}\]

• Germán Rodríguez (2001), Smoothing and non-parametric regression

tps_base_func(r::Number)

Thin plate spline basis function.

• Gianluca Donato and Serge Belongie (2001) Approximate Thin Plate Spline Mappings DOI: 10.1007/3-540-47977-5_2