diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,3 @@ +## 0.1 + +- Initial revision. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2015, Dimitri Sabadie + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Dimitri Sabadie nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/smoothie.cabal b/smoothie.cabal new file mode 100644 --- /dev/null +++ b/smoothie.cabal @@ -0,0 +1,36 @@ +name: smoothie +version: 0.1.0.0 +synopsis: Smooth curves via several spline and polynomials. +description: This package exports several splines and curves you can use + to interpolate points in between. +homepage: https://github.com/phaazon/smoothie +license: BSD3 +license-file: LICENSE +author: Dimitri Sabadie +maintainer: Dimitri Sabadie +copyright: Dimitri Sabadie +category: Data +build-type: Simple + +cabal-version: >=1.10 + +data-files: CHANGELOG.md + +library + + ghc-options: -W -Wall -O2 -funbox-strict-fields + + default-extensions: DeriveFunctor + + exposed-modules: Data.Spline + , Data.Spline.CP + , Data.Spline.Polynomial + + + build-depends: base >= 4.7 && < 4.8 + , linear >= 1.16 && < 1.17 + , vector >= 0.10 && < 0.11 + + hs-source-dirs: src + + default-language: Haskell2010 diff --git a/src/Data/Spline.hs b/src/Data/Spline.hs new file mode 100644 --- /dev/null +++ b/src/Data/Spline.hs @@ -0,0 +1,54 @@ +----------------------------------------------------------------------------- +-- | +-- Copyright : (C) 2015 Dimitri Sabadie +-- License : BSD3 +-- +-- Maintainer : Dimitri Sabadie +-- Stability : experimental +-- Portability : portable +-- +-- A 'Spline s a' represents a curve in which 'a' is very likely to be +-- 'Additive' (see "linear") and 's' is the sampling type. +-- +-- The library exports two useful functions: 'spline' and 'smooth'. The former +-- enables you to create splines while the latter enables you to sample from +-- them using their control points. +---------------------------------------------------------------------------- + +module Data.Spline ( + -- * Spline + Spline + , spline + -- * Smoothing values along splines + , smooth + ) where + +import Data.List ( sortBy ) +import Data.Ord ( comparing ) +import Data.Spline.CP +import Data.Spline.Polynomial ( Polynomial(..), bsearchLower ) +import Data.Vector ( Vector, (!?), fromList ) + +-- |A 'Spline' is a collection of control points with associated polynomials. +-- Given two control points which indices are /i/ and /i+1/, interpolation on +-- the resulting curve is performed using the polynomial of indice /i/. Thus, +-- the latest control point is ignored and can be set to whatever the user wants +-- to, even 'undefined' – you should use hold, though. Yeah, don’t go filthy. +data Spline s a = Spline (Vector (CP s a)) (Vector (Polynomial s a)) + +-- |Create a spline using a list of control points and associated polynomials. +-- Since 'spline' sorts the list before creating the 'Spline', you don’t have to +-- ensure the list is sorted – even though you should, setting control points +-- with no order might be… chaotic. +spline :: (Ord a,Ord s) => [(CP s a,Polynomial s a)] -> Spline s a +spline = uncurry spline_ . unzip . dupLast . sortBy (comparing fst) + where + dupLast s = s ++ [last s] + spline_ cps polys = Spline (fromList cps) (fromList polys) + +-- |Smoothly interpolate a point on a spline. +smooth :: (Ord s) => Spline s a -> s -> Maybe a +smooth (Spline cps polys) s = do + i <- bsearchLower (\(CP s' _) -> compare s s') cps + p <- polys !? i + unPolynomial p s cps diff --git a/src/Data/Spline/CP.hs b/src/Data/Spline/CP.hs new file mode 100644 --- /dev/null +++ b/src/Data/Spline/CP.hs @@ -0,0 +1,23 @@ +----------------------------------------------------------------------------- +-- | +-- Copyright : (C) 2015 Dimitri Sabadie +-- License : BSD3 +-- +-- Maintainer : Dimitri Sabadie +-- Stability : experimental +-- Portability : portable +-- +---------------------------------------------------------------------------- + +module Data.Spline.CP ( + -- * Control points + CP(..) + ) where + +-- | A 'CP' is a **control point**. A curve passes through control points and +-- the shape of the curve is determined by the polynomials used to interpolate +-- values in between. +-- +-- 'CP s a' is a control point of sampling type 's' and carried type 'a'. In +-- most cases, 's' must be 'Ord' and 'a' must be 'Additive' and 'Fractional'. +data CP s a = CP !s !a deriving (Functor,Eq,Ord,Show) diff --git a/src/Data/Spline/Polynomial.hs b/src/Data/Spline/Polynomial.hs new file mode 100644 --- /dev/null +++ b/src/Data/Spline/Polynomial.hs @@ -0,0 +1,98 @@ +----------------------------------------------------------------------------- +-- | +-- Copyright : (C) 2015 Dimitri Sabadie +-- License : BSD3 +-- +-- Maintainer : Dimitri Sabadie +-- Stability : experimental +-- Portability : portable +-- +---------------------------------------------------------------------------- + +module Data.Spline.Polynomial ( + -- * Polynomial + Polynomial(unPolynomial) + -- * Polynomials for interpolation + , hold + , linear + , linearBy + , cosine + -- * Helpers + , bsearchLower + ) where + +import Control.Monad ( guard ) +import Data.Spline.CP +import Data.Vector as V ( Vector, (!?), length ) +import Linear ( Additive(lerp) ) + +-- |A 'Polynomial' is used to interpolate in between a spline’s control points. +newtype Polynomial s a = Polynomial { unPolynomial :: s -> Vector (CP s a) -> Maybe a} + +-- |Constant polynomial – a.k.a. /no interpolation/. +-- +-- Given two control points and a sample value in between, the 'hold' polynomial +-- won’t perform any interpolation but it just /holds/ the value carried by the +-- lower control point along the whole curve between the two control points. +hold :: (Ord s) => Polynomial s a +hold = Polynomial go + where + go s cps = do + li <- bsearchLower (\(CP s' _) -> compare s s') cps + CP _ r <- cps !? li + return r + +-- |Parametric linear polynomial. +-- +-- This form applies a pre-filter on the input before performing a linear +-- interpolation. Instead of: +-- +-- @ lerp x a b @ +-- +-- We have: +-- +-- @ lerp (pref x) a b @ +-- +-- This can be used to implement 1-degree splines if @pref = id@, basic cubic +-- non-hermitian splines if @pref = (^3)@, cosine splines if +-- @pref = \x -> (1 - cos (x*pi)) * 0.5@, and so on and so forth. +linearBy :: (Additive a,Fractional s,Ord s) => (s -> s) -> Polynomial s (a s) +linearBy pref = Polynomial go + where + go s cps = do + li <- bsearchLower (\(CP s' _) -> compare s s') cps + lower <- cps !? li + upper <- cps !? succ li + return $ lerp_ s lower upper + lerp_ x (CP s0 a) (CP s1 b) = lerp x' b a + where + x' = (pref x - s0) / (s1 - s0) + +-- |1-degree polynomial – a.k.a. /straight line interpolation/, or /linear +-- interpolation/. +-- +-- This polynomial connects control points with straight lines. +-- +-- Note: implemented with @linearBy id@. +linear :: (Additive a,Fractional s,Ord s) => Polynomial s (a s) +linear = linearBy id + +-- |Cosine polynomial. +cosine :: (Additive a,Floating s,Ord s) => Polynomial s (a s) +cosine = linearBy $ \x -> (1 - cos (x * pi)) * 0.5 + +-- |Helper binary search that search the ceiling index for the +-- value to be searched according to the predicate. +bsearchLower :: (a -> Ordering) -> Vector a -> Maybe Int +bsearchLower p v = go 0 (pred $ V.length v) + where + go start end = do + guard (start <= end) + ma <- v !? m + ma1 <- v !? succ m + case p ma of + LT -> go start (pred m) + EQ -> Just m + GT -> if p ma1 == LT then Just m else go (succ m) end + where + m = (end + start) `div` 2