Edit on GitHub

LaplaceInterpolation.jl Documentation

This package quickly interpolates data on a grid in one and higher dimensions.

Statement of Need

While numerous implementations of interpolation routines exist that fill missing data points on arbitrary grids, these are largely restricted to one and two dimensions and are slow to run. The implementation we propose is dimension-agnostic, based on a linear-time algorithm, and implements an approximate Matérn kernel interpolation (of which thin plate splines, polyharmonic splines, and radial basis functions are a special case).

Installation Instructions

This package is unregistered, so please install using

Pkg> add https://github.com/vishwas1984/LaplaceInterpolation.jl

Dependencies

Install the following dependencies

Pkg> add TestImages
Pkg> add Images
Pkg> add Plots

Getting started

Suppose we need to interpolate the following image 

using LaplaceInterpolation, TestImages, Images, Random, Plots

img = Float64.(Gray.(testimage("mandrill")))

For illustration purposes, we'll punch a few holes and randomize some data

rows, columns = (256, 512)
N = rows*columns

mat_original = convert(Array{Float64}, img)[1:rows,1:columns]

N2 = Int64(round(N/2))
No_of_nodes_discarded = Int64(round(0.9*N2))

discard1 = N2 .+ randperm(N2)[1:No_of_nodes_discarded]

cent = [(150, 150), (60, 100)]
rad = 30*ones(Int64, 2)
discard2 = punch_holes_2D(cent, rad, rows, columns);

discard = vcat(discard1, discard2)
mat_discard = copy(mat_original)
mat_discard[discard] .= 1

heatmap(mat_discard, title = "Image with Missing data", yflip = true, 
              c = :bone, clims = (0.0, 1.0))

Interpolating using the laplace and matern approximations, we get

restored_img_laplace = matern_2d_grid(mat_discard, discard, 1, 0.0)
restored_img_matern = matern_2d_grid(mat_discard, discard, 2, 0.0)

And plotting, we have

p1 = heatmap(mat_original, title = "Original Data", yflip = true, 
              c = :bone, clims = (0.0, 1.0))
p2 = heatmap(mat_discard, title = "Image with Missing data", yflip = true, 
              c = :bone, clims = (0.0, 1.0))
p3 = heatmap(restored_img_laplace, title = "Laplace Interpolated Image", yflip =
true, c = :bone, clims = (0.0, 1.0))
p4 = heatmap(restored_img_matern, title = "Matern, m = 2, eps = 0.0", yflip =
true, c = :bone, clims = (0.0, 1.0))
plot(p1, p2, p3, p4, layout = (2, 2), legend = false, size = (900, 500))

Mandrill_Random

The Notebooks folder contains this and other examples.

Mathematical Details

Radial basis functions and splines can be unified conceptually through the notion of Green's functions and eigenfunction expansions (Fasshauer, 2012). The general multivariate Matern kernels are of the form

\[K(\mathbf x ; \mathbf z) = K_{m-d/2}(\epsilon||\mathbf x -\mathbf z ||)(ϵ||\mathbf x - \mathbf z ||)^{m-d/2}\]

for $m > d/2$, where $K_ν$ is the modified Bessel function of the second kind, and can be obtained as Green’s kernels of

\[L = (ϵ^2I-Δ)^m \]

where $Δ$ denotes the Laplacian operator in $d$ dimensions. Polyharmonic splines, including thin plate splines, are a special case of the above, and this class includes the thin plate splines.

The discrete gridded interpolation seeks to find an interpolation $u (\mathbf x )$ that satisfies the differential operator in $d$ dimensions on the nodes $\mathbf x_i$ where there is no data and equals $y_i$ everywhere else. Discretely, one solves the matrix problem

\[\mathbf C (\mathbf u - \mathbf y ) - (1 - \mathbf C ) L \mathbf u = 0 \]

where $\mathbf{y}$ contains the $y_i$'s and placeholders where there is no data, $L$ denotes the discrete matrix operator and $C$ is a diagonal matrix that indicates whether node $\mathbf x_i$ is observed.

In $d-$ dimensions the matrix $A^{(d)}$ of size $M \times M$ expands the first order finite difference curvature and its $(i,j)$th entry is -1 when node j is in the set of neighbors of the node $\mathbf x_i$, and has the number of such neighbors on the diagonal. Note that if node $i$ is a boundary node, the $i$-th row of $A^{(d)}$ has $-1$s in the neighboring node spots and the number of such nodes on the diagonal. In general, the rows of $A^{(d)}$ sum to zero.

Denote by $L = A^{(d)}$ the discrete analog of the Laplacian operator. To use the Matern operator, one substitutes

\[L = B^{(d)}(m, ϵ) = ((A^{(d)})^m - ϵ^2 I).\]

Importantly, $A$ is sparse, containing at most 5 nonzero entries per row when $d = 2$ and $7$ nonzero entries per row when $d = 3$ and so on. The Matern matrix $B^{(d)}(m, \epsilon)$ is also sparse, having $2(m+d)-1$ nonzero entries per row. The sparsity of the matrix allows for the interpolation to solve in linear time.

Function Index

One dimensional

LaplaceInterpolation.nablasq_1d_grid — Function
nablasq_1d_grid(n, h, bc)

Laplacian matrix on a 1D grid

Arguments:

- `n`: Number of points
- `h`: Aspect ratio
- `bc`: Boundary conditions (1 implies Neumann BC and 0 implies Do nothing BC)

Outputs:

- discrete Laplacian matrix
source
LaplaceInterpolation.matern_1d_grid — Function
matern_1d_grid(y, idx, m, eps, h, bc)

Matern Interpolation in one dimension

Arguments:

  • y: the vector of y's for which the values are known
  • idx: the vector of indices for which values are to be interpolated
  • m::Int64 = 1: Matern parameter m
  • eps = 0.0: Matern parameter eps
  • h = 1.0: aspect ratio
  • bc: Boundary conditions (1 implies Neumann BC and 0 implies Do nothing BC)

Outputs:

  • vector of interpolated data

Example:

x = 1:100
h = x[2] - x[1]
y = sin.(2 * pi * x * 0.2)
discard = randperm(100)[1:50]
# Laplace interpolation
y_lap = matern_1d_grid(y, discard, 1, 0.0, h)
# Matern interpolation
y_mat = matern_1d_grid(y, discard, 2, 0.1, h)
source

Two dimensional

LaplaceInterpolation.bdy_nodes — Function
bdy_nodes(Nx, Ny)

Find the nodes on the boundary of a 2D array of size (Nx, Ny)

Arguments

  • Nx::Int64: the number of points in the first dimension
  • Ny::Int64: the number of points in the second dimension

Outputs

  • vector containing the indices of coordinates on the boundary of the 2D rectangle
source

bdy_nodes(dims)

Boundary node computation, for arbitrary dimension

Arguments

  • dims::Tuple number of points in each direction, must be ints

Outputs

  • Vector{Int64}: vector containing the indices of coordinates

on the boundary of the hyperrectangle volume

source
LaplaceInterpolation.nablasq_2d_grid — Function
nablasq_2d_grid(Nx, Ny, h, k, bc)

Laplacian matrix on a 2D grid

Arguments:

  • Nx::Int64: Number of points in first dimension
  • Ny::Int64: Number of points in second dimension
  • h::Float64: Aspect ratio in first dimension
  • k::Float64: Aspect ratio in second dimension
  • bc: Boundary conditions (1 implies Neumann BC and 0 implies Do nothing BC)

Outputs:

  • discrete Laplacian matrix in 2D
source
LaplaceInterpolation.matern_2d_grid — Function
matern_2d_grid(mat, discard, m, eps, h, k, bc)

Arguments

  • mat: the matrix containing the image
  • idx: the linear indices of the nodes to be discarded
  • eps: Matern parameter eps
  • m: The Matern exponent (integer)
  • h: The aspect ratio in the first dimension
  • k: The aspect ratio in the second dimension
  • bc: Boundary conditions (1 implies Neumann BC and 0 implies Do nothing BC)

Outputs

  • matrix containing the interpolated image

Example:

x = y = 1:30
h = k = x[2] - x[1]
y = sin.(2 * pi * x * 0.2) * cos.(2 * pi * y * 0.3)
discard = randperm(900)[1:450]

For Laplace interpolation

y_lap = matern_2d_grid(y, discard, 1, 0.0, h, k)

For Matern interpolation

y_mat = matern_2d_grid(y, discard, 2, 0.1, h, k)
source

Three dimensional

LaplaceInterpolation.nablasq_3d_grid — Function
nablasq_3d_grid(Nx,Ny, Nz, h, k, l, bc)

Construct the 3D Laplace matrix

Arguments

  • Nx::Int64: The number of nodes in the first dimension
  • Ny::Int64: The number of nodes in the second dimension
  • Nz::Int64: The number of nodes in the third dimension
  • h::Float64: Grid spacing in the first dimension
  • k::Float64: Grid spacing in the second dimension
  • l::Float64: Grid spacing in the third dimension
  • bc: Boundary conditions (1 implies Neumann BC and 0 implies Do nothing BC)

Outputs

  • -nablasq (discrete Laplacian, real-symmetric positive-definite) on Nx×Ny×Nz grid
source
LaplaceInterpolation.matern_3d_grid — Function
matern_3d_grid(imgg, discard, m, eps, h, k, l, bc)

Interpolates a single punch

Arguments

  • imgg: the matrix containing the image
  • discard::Union{Vector{CartesianIndex{3}}}, Vector{Int64}}: the linear or Cartesian indices of the values to be filled
  • m::Int64 = 1 : Matern parameter
  • eps::Float64 = 0.0: Matern parameter eps
  • h = 1.0: Aspect ratio in the first dimension
  • k = 1.0: Aspect ratio in the second dimension
  • l = 1.0: Aspect ratio in the third dimension
  • bc: Boundary conditions (1 implies Neumann BC and 0 implies Do nothing BC)

Outputs

  • array containing the restored image
source
LaplaceInterpolation.matern_w_punch — Function
matern_w_punch(imgg, Qh_min, Qh_max, Qk_min, Qk_max, Ql_min, Ql_max, m, eps, h, k, l, symm)

Interpolate, in serial, multiple punches

Arguments

  • imgg: the matrix containing the image
  • Qh_min, Qh_max, Qk_min, Qk_max, Ql_min, Ql_max::Int64: the extents in h,k,l resp
  • m::Int64 = 1 : Matern parameter
  • eps::Float64 = 0.0: Matern parameter eps
  • h = 1.0: Aspect ratio in the first dimension
  • k = 1.0: Aspect ratio in the second dimension
  • l = 1.0: Aspect ratio in the third dimension

Outputs

  • array containing the interpolated image

Example

h = k = l = 1.0
symm = 'A'
Nx = Ny = Nz = 61
radius = (0.5, 0.5, 0.5 )
Qh_min = Qk_min = Ql_min = -3.0
Qh_max = Qk_max = Ql_max = 3.0
xpoints = ypoints = zpoints = LinRange(Qh_min, Qh_max, Nx)
imgg = rand(Nx, Ny, Nz)
m = 1
eps = 0.0
interp = matern_w_punch(imgg, #Qh_min, Qh_max, Qk_min, Qk_max, Ql_min, Ql_max, radius,
                      radius, xpoints, ypoints, zpoints, m, eps, 
                      h, k, l, symm)
source

Arbitrary dimensions

LaplaceInterpolation.nablasq_arb — Function
nablasq_arb(dims, aspect_ratios)

Compute the interpolating (Laplace) matrix for a volume with given dimensions

Arguments

  • dims::Tuple number of points in each dimension, must be ints
  • aspect_ratios::Tuple the aspect ratio in each dimension, in units

Outputs

  • A::Matrix{Float64}: matrix of size prod(dims) containing the Laplacian

Note that the boundary condition is average, and Neumann has not been implemented.

source
LaplaceInterpolation.interp — Function
interp(data, ind, m, eps, aspect_ratios)

Compute a Matern or Laplace interpolation on a multidimensional dataset

Arguments

  • data::Array: Array containing the data
  • ind::Array: Array of indices for which the data are to be interpolated
  • m::Int64: Matern parameter m (m = 1 for Laplace interpolation)
  • eps::Float64: Matern parameter epsilon (eps = 0.0 for Laplace)
  • aspect_ratios::Tuple: Tuple containing the aspect ratios in each of the dimensions

Outputs

  • u::Array: Array (same size as data) containing the interpolated result.
source

Spherical Punching

LaplaceInterpolation.punch_holes_3D — Function
punch_holes_3D(centers, radius, xpoints, ypoints, zpoints)

Punch holes in a 3D volume dataset

Arguments

  • centers::Vector{T}: the vector containing the centers of the punches
  • radius::Float64: the radius of the punch
  • Nx::Int64: the number of points in the x-direction, this code is hard-coded to start from one.
  • Ny::Int64: the number of points in the y-direction
  • Nz::Int64: the number of points in the z-direction

Outputs

  • absolute_indices::Vector{Int64}: vector containing the indices of coordinates

inside the punch

source
LaplaceInterpolation.punch_holes_2D — Function
punch_holes_2D(centers, radius, xpoints, ypoints)

Punch holes in a 2D dataset

Arguments

  • centers::Union{Vector}: the vector containing the punch centers
  • radius::Vector: the tuple containing the punch radii
  • Nx::Int64: the number of points in the x direction
  • Ny::Int64: the number of points in the y direction

Outputs

  • absolute_indices::Vector{Int64}: vector containing the indices of coordinates

inside the punch

source
LaplaceInterpolation.punch_3D_cart — Function
punch_3D_cart(center, radius, xpoints, ypoints, zpoints; <kwargs>)

This punch gives outputs either as linear or Cartesian indices

Arguments

  • center::Tuple{T}: the tuple containing the center of a round punch
  • radius::Union{Tuple{Float64},Float64}: the radii/radius of the punch
  • x::Vector{T} where T<:Real: the vector containing the x coordinate
  • y::Vector{T} where T<:Real: the vector containing the y coordinate
  • z::Vector{T} where T<:Real: the vector containing the z coordinate

Optional

  • linear::Bool = false can return a linear index

Outputs

  • inds::Vector{Int64}: vector containing the indices of coordinates

inside the punch

source