# Implementing recursive filters with Haskell/Repa

I recently learned Haskell, and I am trying to apply it to the code I use in order to get a feeling for the language.

I really like the Repa library since I manipulate a lot of multi-dimensional data. Yet it is, in my sense, missing a lot of signal processing tools. For instance, I use quite a lot of Gaussian filter combinations, and approximate them with the Canny-Deriche recursive filter, so that the complexity remains linear and independent of the deviation.

So, this is the code that I have started to write, only for Gaussian convolutions, but in any dimensions.

import Control.Monad
import Data.List (tails)

import Data.Array.Repa as R hiding ((++), map, zipWith)
import Data.Array.Repa.Eval
import Data.Array.Repa.Eval.Gang
import Data.Array.Repa.Repr.Unboxed
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as UM

-- |Apply a recursive Canny-Deriche filter on one dimension of a Repa array.
deriche :: Shape sh =>
Double -> Int -> Int -> Bool ->
Array U sh Double ->
IO (Array U sh Double)
deriche s o d bc a
-- Convolution with a dirac is a nop
| s == 0 = return a
-- Convolution only make sense for positive sigmas
| s < 0 = error "Sigma for the filter must be positive"
-- Make sure the selected dimension exists
| d >= rank (extent a) =
error "Invalid array dimension selected for the filter"
| otherwise =
case o of
0 -> let k  = (1-ena) * (1-ena) / (1 + 2*alpha*ena - ens)
a0 = k
a1 = k * ena * (alpha - 1)
a2 = k * ena * (alpha + 1)
a3 = -k * ens
par = 1
in do
b <- UM.new $size$ extent a :: IO (UM.IOVector Double)
innerLoop a0 a1 a2 a3 0 [] True b
liftM (fromUnboxed $extent a) (U.unsafeFreeze b) otherwise -> error "Unimplemented filter order" where alpha = 1.695 / s expAlpha = exp alpha ena = exp (-alpha) ens = ena * ena b1 = -2 * ena b2 = ens ndim = rank sha dims = reverse$ listOfShape sha
ddim = dims!!d
stride = product $tails dims !! (d+1) sha = extent a innerLoop :: Double -> Double -> Double -> Double -> Int -> [Int] -> Bool -> UM.IOVector Double -> IO () innerLoop a0 a1 a2 a3 l path pl b | l < ndim && l == d = innerLoop a0 a1 a2 a3 (l+1) (0:path) pl b | l < ndim && pl = let len = dims!!l threads = gangSize theGang chunkLen = len quot threads chunkOver = len rem threads splitIdx thread | thread < chunkOver = thread * (chunkLen + 1) | otherwise = thread * chunkLen + chunkOver fill :: Int -> Int -> IO () fill ix end = mapM_ (\i -> innerLoop a0 a1 a2 a3 (l+1) (i:path) False b) [ix..end-1] in gangIO theGang$ \thread ->
end = splitIdx (thread + 1)
in fill start end
| l < ndim =
mapM_ (\i -> innerLoop a0 a1 a2 a3 (l+1) (i:path) pl b) [0..dims!!l-1]
| otherwise = do
-- First pass
if bc
then undefined --pass1 0 (a linearIndex idx)
else pass1 0 idx 0 0 0
-- Second pass
if bc
then undefined -- pass2 ...
else pass2 (ddim-1) (idx+(ddim-1)*stride) 0 0 0 0
where idx = toIndex sha (shapeOfList \$ cop path)
cop = id
pass1 !i !j !xp !yp !yb
| i == ddim = return ()
| otherwise = do
UM.unsafeWrite b j yc
pass1 (i+1) (j+stride) xc yc yp
where yc = a0*xc + a1*xp - b1*yp - b2*yb
xc = a unsafeLinearIndex j
{-# INLINE pass1 #-}
pass2 !i !j !xn !xa !yn !ya
| i < 0 = return ()
| otherwise = do
UM.unsafeWrite b j (bj + yc)
pass2 (i-1) (j-stride) xc xn yc yn
where yc = a2*xn + a3*xa - b1*yn - b2*ya
xc = a unsafeLinearIndex j
{-# INLINE pass2 #-}
{-# INLINE innerLoop #-}

-- |Sequentially apply a Canny-Deriche filter on two dimensions.
deriche2D :: Double -> Int -> Bool ->
Array U DIM2 Double ->
IO (Array U DIM2 Double)
deriche2D s o bc = deriche s o 0 bc >=> deriche s o 1 bc


Since I am a complete novice in Haskell, I would appreciate any comment you might give. In particular the code I developed is quite un-repa-ish due to the recursive filter. Also I have no clue as to when I should use bangs, since adding them to any array variable has no positive effect. Any idea on how to share the arrays for deriche2D would be greatly appreciated.

Also, I am compiling with -rtsopts -O2 -threaded -fllvm and running with +RTS -N4.