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

Interpolate.Cubic

This library interpolates cubic splines for one-dimensional sets of data.

It computes "natural splines", which means the second derivative at the endpoints is set to zero.

For more information on the mathematics used, check out this paper.

Creating splines

withRange : Float -> Float -> List Float -> Spline

Compute a spline, given the minimum and maximum values of x and a list of data for f(x). The data should be evenly spaced and in order of increasing x.

For example, if we have the data

f(2) = 1
f(3) = 5.2
f(4) = 3.2
f(5) = 0.8

Then we would generate a spline by calling

fSpline = withRange 2 6 [ 1, 5.2, 3.2, 0.8 ]

If there is only one data point, then the spline will be a horizontal line passing through that point. If the data is empty, the spline will be zero everywhere.

withDelta : Float -> Float -> List Float -> Spline

Same as withRange, except instead of passing the endpoint as the second argument, you pass the x-distance between data points.

fSpline = withDelta 2 1 [ 1, 5.2, 3.2, 0.8 ]
-- equivalent to withRange 2 6 [ 1, 5.2, 3.2, 0.8 ]
type Spline = Spline { coefficients : Array Coefficients , start : Float , dx : Float }

Interpolating

valueAt : Float -> Spline -> Float

Given x1 and a spline f(x), return f(x1)

slopeAt : Float -> Spline -> Float

Return the first derivative of the curve at the given point

concavityAt : Float -> Spline -> Float

Return the second derivative of the curve at the given point

curveAt : Float -> Spline -> LocalCurve

Return a record containing the value, slope, and concavity of the curve at the given point

type alias LocalCurve = { value : Float , slope : Float , concavity : Float }
module Interpolate.Cubic
    exposing
        ( Spline
        , LocalCurve
        , withRange
        , withDelta
        , valueAt
        , slopeAt
        , concavityAt
        , curveAt
        )

{-| This library interpolates cubic splines for one-dimensional sets of data.

It computes "natural splines", which means the second derivative at the endpoints
is set to zero.

For more information on the mathematics used, check out
[this paper.](http://web.archive.org/web/20090408054627/http://online.redwoods.cc.ca.us/instruct/darnold/laproj/Fall98/SkyMeg/Proj.PDF)

# Creating splines
@docs withRange, withDelta, Spline

# Interpolating
@docs valueAt, slopeAt, concavityAt, curveAt, LocalCurve
-}

import Array exposing (Array)


{-| Given `x1` and a spline `f(x)`, return `f(x1)`
-}
valueAt : Float -> Spline -> Float
valueAt =
    evaluate cubic


{-| Return the first derivative of the curve at the given point
-}
slopeAt : Float -> Spline -> Float
slopeAt =
    evaluate firstDerivative


{-| Return the second derivative of the curve at the given point
-}
concavityAt : Float -> Spline -> Float
concavityAt =
    evaluate secondDerivative


{-| Return a record containing the value, slope, and concavity of the curve
at the given point
-}
curveAt : Float -> Spline -> LocalCurve
curveAt =
    let
        curve x coeff =
            { value = cubic x coeff
            , slope = firstDerivative x coeff
            , concavity = secondDerivative x coeff
            }
    in
        evaluate curve


cubic : Float -> Coefficients -> Float
cubic x { a, b, c, d } =
    a * x ^ 3 + b * x ^ 2 + c * x + d


firstDerivative : Float -> Coefficients -> Float
firstDerivative x { a, b, c, d } =
    3 * a * x ^ 2 + 2 * b * x + c


secondDerivative : Float -> Coefficients -> Float
secondDerivative x { a, b, c, d } =
    6 * a * x + 2 * b


evaluate : (Float -> Coefficients -> a) -> Float -> Spline -> a
evaluate f x (Spline spline) =
    let
        maxIndex =
            Array.length spline.coefficients - 1

        index =
            (x - spline.start)
                / spline.dx
                |> floor
                |> clamp 0 maxIndex

        offset =
            x - spline.dx * (toFloat index) - spline.start

        coeff =
            Array.get index spline.coefficients
                |> Maybe.withDefault { a = 0, b = 0, c = 0, d = 0 }
    in
        f offset coeff


{-| -}
type alias LocalCurve =
    { value : Float
    , slope : Float
    , concavity : Float
    }


{-| Compute a spline, given the minimum and maximum values of `x` and a
list of data for `f(x)`. The data should be evenly spaced and in order of
increasing `x`.

For example, if we have the data

    f(2) = 1
    f(3) = 5.2
    f(4) = 3.2
    f(5) = 0.8

Then we would generate a spline by calling

    fSpline = withRange 2 6 [ 1, 5.2, 3.2, 0.8 ]

If there is only one data point, then the spline will be a horizontal line
passing through that point. If the data is empty, the spline will be zero
everywhere.
-}
withRange : Float -> Float -> List Float -> Spline
withRange start end heights =
    let
        n =
            List.length heights |> toFloat

        dx =
            (end - start) / (n - 1)
    in
        if 1 < n then
            withDelta start dx heights
        else
            { coefficients = degenerateCoefficients heights
            , start = 0
            , dx = 0
            }
                |> Spline


{-| Same as `withRange`, except instead of passing the endpoint as the
second argument, you pass the x-distance between data points.

    fSpline = withDelta 2 1 [ 1, 5.2, 3.2, 0.8 ]
    -- equivalent to withRange 2 6 [ 1, 5.2, 3.2, 0.8 ]
-}
withDelta : Float -> Float -> List Float -> Spline
withDelta start dx heights =
    { coefficients = findCoefficients dx heights
    , start = start
    , dx = dx
    }
        |> Spline


degenerateCoefficients : List Float -> Array Coefficients
degenerateCoefficients heights =
    { a = 0
    , b = 0
    , c = 0
    , d = List.head heights |> Maybe.withDefault 0
    }
        |> Array.repeat 1


findCoefficients : Float -> List Float -> Array Coefficients
findCoefficients dx heights =
    let
        concavity y0 y1 y2 =
            6 * (y0 + y2 - 2 * y1) / dx ^ 2

        sweep y ( uPrev, vPrev ) =
            ( 1 / (4 - uPrev), (y - vPrev) / (4 - uPrev) )

        backSub ( u, v ) mNext =
            v - u * mNext

        piece ( y0, m0 ) ( y1, m1 ) =
            { a = (m1 - m0) / 6 / dx
            , b = m0 / 2
            , c = (y1 - y0) / dx - (m1 + 2 * m0) * dx / 6
            , d = y0
            }
    in
        mapTriple concavity heights
            |> List.scanl sweep ( 0, 0 )
            |> scanr backSub 0
            |> List.map2 (,) heights
            |> mapPair piece
            |> Array.fromList


mapTriple : (a -> a -> a -> b) -> List a -> List b
mapTriple f x0 =
    let
        tail =
            List.tail x0

        doubleTail =
            tail |> Maybe.andThen List.tail
    in
        case ( tail, doubleTail ) of
            ( Just x1, Just x2 ) ->
                List.map3 f x0 x1 x2

            otherwise ->
                []


mapPair : (a -> a -> b) -> List a -> List b
mapPair f x0 =
    case (List.tail x0) of
        Just x1 ->
            List.map2 f x0 x1

        otherwise ->
            []


scanr : (a -> b -> b) -> b -> List a -> List b
scanr f init =
    List.reverse >> List.scanl f init >> List.reverse


{-| -}
type Spline
    = Spline
        { coefficients : Array Coefficients
        , start : Float
        , dx : Float
        }


type alias Coefficients =
    { a : Float
    , b : Float
    , c : Float
    , d : Float
    }