This is an alternative site for discovering Elm packages. You may be looking for the official Elm package site instead.

# Interpolate.Bicubic

This module uses bicubic splines to interpolate surfaces. Bicubic interpolation is the 2D equivalent of cubic interpolation.

For example, if we have the following data:

``````y=|
9 | 1.9  2.1  5.5
6 | 2.2  3.1  5.4
3 | 3.8  4.0  4.3
0 | 1.1  2.3  3.8
-----------------
x=  2    4    6
``````

We could construct a 2D spline this way:

``````data = rows [ [ 1.1, 2.3, 3.8 ]
, [ 3.8, 4.0, 4.3 ]
, [ 2.2, 3.1, 5.4 ]
, [ 1.9, 2.1, 5.5 ]
]
|> Maybe.withDefault emptyData

start = { x = 2, y = 0 }
end = { x = 6, y = 9 }
delta = { x = 2, y = 3 }

-- These two splines are equivalent
splineOne = withRange start end data
splineTwo = withDelta start delta data
``````

# Creating data sets

rows : List (List Float) -> Maybe Data

Construct a two-dimensional data set. The input is given in x-major order, so a grid of values like this:

``````f(0,0) = 1
f(1,0) = 2
f(0,1) = 3
f(1,1) = 4
``````

Should be passed to the function in this format:

``````rows [ [1,2], [3,4] ] -- returns Just Data
``````

If every row has the same length, the function returns a `Data` object that can be used to build a spline. If the rows are uneven (or empty), it returns `Nothing`.

``````rows [ [1,1,1], [1,1,1] ] -- returns Just Data

rows [ [1], [1], [1] ] -- returns Just Data

rows [ [1,1,1], [1] ] -- returns Nothing

rows [ [], [] ] -- returns Nothing

rows [ ] -- returns Nothing
``````
emptyData : Data

An empty data set. Useful as a default value when working with the `rows` function. If you create a spline from it, the spline will be zero everywhere.

type Data = Data (Matrix Float)

# Creating splines

withRange : Vector -> Vector -> Data -> Spline

Construct a spline, given the positions of the lower left (min-x, min-y) and upper right (max-x, max-y) data samples, and a data set.

withDelta : Vector -> Vector -> Data -> Spline

Construct a spline, given the position of the lower left data sample and the dimensions of a grid cell.

type alias Vector = { x : Float , y : Float }

Stores data with an x and a y component.

type Spline = Spline { coefficients : Matrix Coefficients , start : Vector , delta : Vector }

# Interpolating

valueAt : Vector -> Spline -> Float

Evaluate the spline at the given point

gradientAt : Vector -> Spline -> Vector

Compute the gradient of the spline at the given point. The gradient is the x and y partial derivatives of the spline.

laplacianAt : Vector -> Spline -> Float

Compute the Laplacian of the spline at the given point. The Laplacian is the divergence of the gradient. It is computed by adding the x and y second partial derivatives.

surfaceAt : Vector -> Spline -> LocalSurface

Returns `valueAt`, `gradientAt`, and `laplacianAt` results in a single record.

type alias LocalSurface = { value : Float , gradient : Vector , laplacian : Float }
``````module Interpolate.Bicubic
exposing
( Data
, Vector
, Spline
, LocalSurface
, rows
, emptyData
, withRange
, withDelta
, valueAt
, surfaceAt
, laplacianAt
)

{-| This module uses [bicubic splines](https://en.wikipedia.org/wiki/Bicubic_interpolation)
to interpolate surfaces. Bicubic interpolation is the 2D equivalent
of cubic interpolation.

For example, if we have the following data:

y=|
9 | 1.9  2.1  5.5
6 | 2.2  3.1  5.4
3 | 3.8  4.0  4.3
0 | 1.1  2.3  3.8
-----------------
x=  2    4    6

We could construct a 2D spline this way:

data = rows [ [ 1.1, 2.3, 3.8 ]
, [ 3.8, 4.0, 4.3 ]
, [ 2.2, 3.1, 5.4 ]
, [ 1.9, 2.1, 5.5 ]
]
|> Maybe.withDefault emptyData

start = { x = 2, y = 0 }
end = { x = 6, y = 9 }
delta = { x = 2, y = 3 }

-- These two splines are equivalent
splineOne = withRange start end data
splineTwo = withDelta start delta data

# Creating data sets
@docs rows, emptyData, Data

# Creating splines
@docs withRange, withDelta, Vector, Spline

# Interpolating
@docs valueAt, gradientAt, laplacianAt, surfaceAt, LocalSurface
-}

import Array exposing (Array)
import Interpolate.Cubic as Cubic
import Matrix exposing (Matrix, Quad)

{-| Construct a two-dimensional data set. The input is given in x-major order,
so a grid of values like this:

f(0,0) = 1
f(1,0) = 2
f(0,1) = 3
f(1,1) = 4

Should be passed to the function in this format:

rows [ [1,2], [3,4] ] -- returns Just Data

If every row has the same length, the function returns a `Data` object that
can be used to build a spline. If the rows are uneven (or empty), it returns
`Nothing`.

rows [ [1,1,1], [1,1,1] ] -- returns Just Data

rows [ [1], [1], [1] ] -- returns Just Data

rows [ [1,1,1], [1] ] -- returns Nothing

rows [ [], [] ] -- returns Nothing

rows [ ] -- returns Nothing
-}
rows : List (List Float) -> Maybe Data
rows =
Matrix.fromRows >> Maybe.map Data

{-| An empty data set. Useful as a default value when working
with the `rows` function. If you create a spline from it, the
spline will be zero everywhere.
-}
emptyData : Data
emptyData =
Data defaultMatrix

{-| Evaluate the spline at the given point
-}
valueAt : Vector -> Spline -> Float
valueAt =
evaluate cubic

of the spline at the given point. The gradient is the x and y
partial derivatives of the spline.
-}
gradientAt : Vector -> Spline -> Vector
evaluate del

{-| Compute the [Laplacian](https://en.wikipedia.org/wiki/Laplace_operator)
of the spline at the given point. The Laplacian is the divergence
of the gradient. It is computed by adding the x and y second partial
derivatives.
-}
laplacianAt : Vector -> Spline -> Float
laplacianAt =
evaluate delDotDel

{-| Returns `valueAt`, `gradientAt`, and `laplacianAt` results in
a single record.
-}
surfaceAt : Vector -> Spline -> LocalSurface
surfaceAt =
evaluate
(\pt coeff ->
{ value = cubic pt coeff
, gradient = del pt coeff
, laplacian = delDotDel pt coeff
}
)

evaluate : (Vector -> Coefficients -> a) -> Vector -> Spline -> a
evaluate f point (Spline spline) =
let
maxIndex =
{ x = Matrix.width spline.coefficients - 1
, y = Matrix.height spline.coefficients - 1
}

recentered =
pointMap (-) point spline.start

normalized =
pointMap (/) recentered spline.delta

index =
pointMap (floor >> flip (clamp 0)) normalized maxIndex

offset =
pointMap (toFloat >> (*)) index spline.delta
|> pointMap (-) recentered
in
Matrix.get index.x index.y spline.coefficients
|> Maybe.withDefault defaultMatrix
|> f offset

cubic : Vector -> Coefficients -> Float
cubic { x, y } =
(\i j factor ->
factor * x ^ i * y ^ j
)

del : Vector -> Coefficients -> Vector
del { x, y } coeff =
let
dz_dx i j factor =
if (0 < i) then
factor * i * x ^ (i - 1) * y ^ j
else
0

dz_dy i j factor =
if (0 < j) then
factor * j * x ^ i * y ^ (j - 1)
else
0
in
{ x = addMonomials dz_dx coeff
, y = addMonomials dz_dy coeff
}

delDotDel : Vector -> Coefficients -> Float
delDotDel { x, y } coeff =
let
d2z_dx2 i j factor =
if (1 < i) then
factor * i * (i - 1) * x ^ (i - 2) * y ^ j
else
0

d2z_dy2 i j factor =
if (1 < j) then
factor * j * (j - 1) * x ^ i * y ^ (j - 2)
else
0

derivTotal i j factor =
(d2z_dx2 i j factor) + (d2z_dy2 i j factor)
in

addMonomials : (Float -> Float -> Float -> Float) -> Coefficients -> Float
addMonomials monomial coeff =
Matrix.indexedMap (\i j -> monomial (toFloat i) (toFloat j)) coeff
|> Matrix.sum

{-| Construct a spline, given the positions of the lower left (min-x, min-y)
and upper right (max-x, max-y) data samples, and a data set.
-}
withRange : Vector -> Vector -> Data -> Spline
withRange start end (Data data) =
let
deltaFor dim elem =
((elem end) - (elem start)) / (toFloat (dim data) - 1)

delta =
{ x = deltaFor Matrix.width .x
, y = deltaFor Matrix.height .y
}

bigEnough =
(Matrix.width data > 1) && (Matrix.height data > 1)
in
if bigEnough then
withDelta start delta (Data data)
else
{ coefficients = degenerateCoefficients data
, start = { x = 0, y = 0 }
, delta = { x = 0, y = 0 }
}
|> Spline

{-| Construct a spline, given the position of the lower left data sample
and the dimensions of a grid cell.
-}
withDelta : Vector -> Vector -> Data -> Spline
withDelta start delta (Data data) =
{ coefficients = findCoefficients delta data
, start = start
, delta = delta
}
|> Spline

degenerateCoefficients : Matrix Float -> Matrix Coefficients
degenerateCoefficients data =
let
height =
Matrix.get 0 0 data
|> Maybe.withDefault 0
in
Matrix.repeat 4 4 0
|> Matrix.set 0 0 height
|> Matrix.repeat 1 1

findCoefficients : Vector -> Matrix Float -> Matrix Coefficients
findCoefficients delta data =
findDerivatives data
|> Maybe.map (Matrix.map (fromDerivatives delta))
|> Maybe.withDefault (degenerateCoefficients data)

findDerivatives : Matrix Float -> Maybe (Matrix Derivatives)
findDerivatives data =
let
xDerivs =
Matrix.rows data
|> List.map (Cubic.withDelta 0 1)
|> List.map (slopes (Matrix.width data))
|> Matrix.fromRows
|> withDefaultMatrix

yDerivs =
Matrix.columns data
|> List.map (Cubic.withDelta 0 1)
|> List.map (slopes (Matrix.height data))
|> Matrix.fromColumns
|> withDefaultMatrix

xyDerivs =
Matrix.columns xDerivs
|> List.map (Cubic.withDelta 0 1)
|> List.map (slopes (Matrix.height data))
|> Matrix.fromColumns
|> withDefaultMatrix

pack z dz_dx dz_dy dz_dx_dy =
{ z = z
, dz_dx = dz_dx
, dz_dy = dz_dy
, dz_dx_dy = dz_dx_dy
}
in
Matrix.map4 pack data xDerivs yDerivs xyDerivs

slopes : Int -> Cubic.Spline -> List Float
slopes n spline =
(toFloat >> flip Cubic.slopeAt spline)
|> Array.initialize n
|> Array.toList

fromDerivatives : Vector -> Quad Derivatives -> Coefficients
fromDerivatives delta d =
let
derivsMatrix =
[ [ d.n.w.z
, d.n.e.z
, d.s.w.z
, d.s.e.z
, d.n.w.dz_dx
, d.n.e.dz_dx
, d.s.w.dz_dx
, d.s.e.dz_dx
, d.n.w.dz_dy
, d.n.e.dz_dy
, d.s.w.dz_dy
, d.s.e.dz_dy
, d.n.w.dz_dx_dy
, d.n.e.dz_dx_dy
, d.s.w.dz_dx_dy
, d.s.e.dz_dx_dy
]
]
|> Matrix.fromColumns
|> withDefaultMatrix

factors =
[ [ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ]
, [ 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ]
, [ -3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ]
, [ 2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ]
, [ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 ]
, [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 ]
, [ 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0 ]
, [ 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0 ]
, [ -3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0, 0, 0, 0, 0 ]
, [ 0, 0, 0, 0, -3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0 ]
, [ 9, -9, -9, 9, 6, 3, -6, -3, 6, -6, 3, -3, 4, 2, 2, 1 ]
, [ -6, 6, 6, -6, -3, -3, 3, 3, -4, 4, -2, 2, -2, -2, -1, -1 ]
, [ 2, 0, -2, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0 ]
, [ 0, 0, 0, 0, 2, 0, -2, 0, 0, 0, 0, 0, 1, 0, 1, 0 ]
, [ -6, 6, 6, -6, -4, -2, 4, 2, -3, 3, -3, 3, -2, -1, -2, -1 ]
, [ 4, -4, -4, 4, 2, 2, -2, -2, 2, -2, 2, -2, 1, 1, 1, 1 ]
]
|> Matrix.fromRows
|> withDefaultMatrix

scale i j factor =
factor / delta.x ^ (toFloat i) / delta.y ^ (toFloat j)
in
Matrix.times factors derivsMatrix
|> Maybe.map squarify
|> Maybe.map (Matrix.indexedMap scale)
|> withDefaultMatrix

squarify : Matrix Float -> Matrix Float
squarify matrix =
let
slices flat =
[ Array.slice 0 4 flat
, Array.slice 4 8 flat
, Array.slice 8 12 flat
, Array.slice 12 16 flat
]
|> List.map Array.toList
in
Matrix.columns matrix
|> List.concat
|> Array.fromList
|> slices
|> Matrix.fromRows
|> withDefaultMatrix

withDefaultMatrix : Maybe (Matrix Float) -> Matrix Float
withDefaultMatrix =
Maybe.withDefault defaultMatrix

defaultMatrix : Matrix Float
defaultMatrix =
Matrix.repeat 1 1 0

pointMap : (a -> b -> c) -> { x : a, y : a } -> { x : b, y : b } -> { x : c, y : c }
pointMap op a b =
{ x = op a.x b.x
, y = op a.y b.y
}

{-| -}
type alias LocalSurface =
{ value : Float
, gradient : Vector
, laplacian : Float
}

{-| -}
type Data
= Data (Matrix Float)

{-| Stores data with an x and a y component.
-}
type alias Vector =
{ x : Float
, y : Float
}

{-| -}
type Spline
= Spline
{ coefficients : Matrix Coefficients
, start : Vector
, delta : Vector
}

type alias Coefficients =
Matrix Float

type alias Derivatives =
{ z : Float
, dz_dx : Float
, dz_dy : Float
, dz_dx_dy : Float
}
```
```