bundles/org.simantics.scl.runtime/scl/Interpolation.scl

   1 import "Prelude"
   2 include "Vector"
   3
   4
   5 @private
   6 @inline
   7 chooseInterval :: Vector Double -> Double -> (<e> a) -> (Integer -> <e> a) -> (<e> a) -> <e> a
   8 chooseInterval table x down inRange up
   9     | x < table!0        = down
  10     | x >= table!(len-1) = up
  11     | otherwise          = search 0 (len-1)
  12   where
  13     len = length table
  14     search low high | high <= low+1 = inRange low
  15                     | otherwise = do
  16        middle = (low + high) `div` 2
  17        if x < table!middle
  18        then search low middle
  19        else search middle high
  20
  21
  22 interpolateNearest :: Vector Double -> Vector Double -> Double -> Double
  23 interpolateNearest xs ys x = chooseInterval xs x
  24     (ys!0)
  25     (\i -> if x - xs!i < xs!(i+1) - x then ys!i else ys!(i+1))
  26     (ys!(length xs-1))
  27
  28 interpolateLinear :: Vector Double -> Vector Double -> Double -> Double
  29 interpolateLinear xs ys x = chooseInterval xs x
  30     (ys!0)
  31     (\i -> lin (xs!i) (xs!(i+1)) (ys!i) (ys!(i+1)))
  32     (ys!(length ys-1))
  33   where
  34     lin x0 x1 y0 y1 = y0 + (y1-y0) * (x-x0) / (x1-x0)
  35
  36 interpolateCubic :: Vector Double -> Vector Double -> Double -> Double
  37 interpolateCubic xs ys x = chooseInterval xs x
  38     (ys!0)
  39     (\i -> let a = max 0 (i-1) ; b = min (i+2) (length xs-1) in
  40         cubic (xs!a) (xs!i) (xs!(i+1)) (xs!b)
  41               (ys!a) (ys!i) (ys!(i+1)) (ys!b))
  42     (ys!(length ys-1))
  43   where
  44     cubic x0 x1 x2 x3 y0 y1 y2 y3 = let
  45         l = x2-x1
  46         dx1 = (y2-y0) * l / (x2-x0)
  47         dx2 = (y3-y1) * l / (x3-x1)
  48         t = (x-x1) / l
  49       in
  50         y1 + t * (dx1 + t * ((-3*y1 + 3*y2 - 2*dx1 - dx2) + t * (2*y1 - 2*y2 + dx1 + dx2)))