movarec corner detection in accelerate

movarec.hs

module Moravec where
import Data.Array.Accelerate as A
import Data.Array.Accelerate.Interpreter as A
import Data.Array.Accelerate.IO as IO

moravec :: Exp DIM2               -- ^ neighbor region size;
        -> Exp DIM2               -- ^ patch sizes;
        -> Acc (Array DIM2 Float) -- ^ a grayscale image;
        -> Acc (Array DIM2 Float) -- ^ corner intensities.
moravec n p i = fold2 min (1/0) $ A.map unzero
              $ fold2 (+) 0 $ A.map sqr $ A.zipWith (-) nRepl iRepl
  where
    sqr x = x * x
    unzero x = cond (x ==* 0) (1/0) x -- a dirty hack to exclude original patch
    fold2 f z = A.fold f z . A.fold f z

    nRepl = A.backpermute intermediateSh fn i
    iRepl = A.backpermute intermediateSh fi i

    fn :: Exp DIM6 -> Exp DIM2 -- neighborhood index permutation
    fn ix = lift (Z :. (j + u + y) :. (i + v + x))
      where
        j, i, u, v, y, x :: Exp Int
        Z :. j :. i :. u :. v :. y :. x = unlift ix

    fi :: Exp DIM6 -> Exp DIM2 -- image index permutation
    fi ix = lift (Z :. (j + y) :. (i + x))
      where
        j, i, u, v, y, x :: Exp Int
        Z :. j :. i :. u :. v :. y :. x = unlift ix

    intermediateSh = lift (Z :. rh :. rw :. nh :. nw :. ph :. pw)

    rh, rw :: Exp Int
    rh = ih - (nh - 1) - (ph - 1) -- result image height (crop mode)
    rw = iw - (nw - 1) - (pw - 1) -- result image width  (crop mode)

    nh, nw, ph, pw, ih, iw :: Exp Int
    Z :. nh :. nw = unlift n         -- neighborhood height and width
    Z :. ph :. pw = unlift p         -- patch height and width
    Z :. ih :. iw = unlift (shape i) -- image  height and width

main :: IO ()
main = do
  let path suf = "image" Prelude.++ suf Prelude.++ ".bmp"

  putStrLn "Loading original image..."
  Right orig <- readImageFromBMP (path "orig")

  putStrLn "Color to grayscale convertion..."
  let gray = A.run $ grayscale (use orig)
  writeImageToBMP (path "gray") $ A.run $ A.map rgba32OfLuminance (use gray)

  putStrLn "Finding points of interests"
  let corners = A.run $ moravec (index2 3 3) (index2 3 3) (use gray)
  writeImageToBMP (path "corners") $ A.run $ A.map rgba32OfLuminance (use corners)