Implementing recursive filters with Haskell/Repa

Question

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 ->
                                     let start = splitIdx 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
                            bj <- UM.unsafeRead b j
                            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.