id	summary	reporter	owner	description	type	status	priority	milestone	component	version	resolution	keywords	cc
15	Support for stencil-based traversal in the language front-end	blever	chak	"A ""stencil"" primitive would be a very useful in the language front-end of Accelerate.
This type of building-block is very useful in applications such as computer
vision which typically involves traversing an array where target elements are
updated by applying some function to neighbouring elements in a source array.

Below are propsed type definitions for stencil primitives that could be added to Accelerate:

{{{
-- |Traversal over a single array using stencils.
stencil
  :: (Ix dim, Ix dim', Ix dimS, Elem a, Elem b)
  => Acc (Array dim a)                                  -- Source array
  -> (Exp dim -> Exp dim')                              -- Function specifying extent of destination array
  -> (Exp dim -> Exp dimS)                              -- Function specifying extent of stencil array
  -> (Exp dimS -> Exp dim' -> Exp dim)                  -- Stencil specification function
  -> (Exp dim' -> Acc (Array dimS a) -> Acc (Scalar b)) -- Populate/update function
  -> Acc (Array dim' b)                                 -- Destination array

-- |Traversal over two arrays at once using stencils.
stencil2
  :: (Ix dim1, Ix dim2, Ix dim', Ix dimS1, Ix dimS2, Elem a, Elem b, Elem c)
  => Acc (Array dim1 a)                                                         -- First source array
  -> Acc (Array dim2 b)                                                         -- Second source array
  -> (Exp dim1 -> Exp dim2 -> Exp dim')                                         -- Function specifying extent of destination array
  -> (Exp dim1 -> Exp dimS1)                                                    -- Function specifying extent of the first stencil array
  -> (Exp dim2 -> Exp dimS2)                                                    -- Function specifying extent of the second stencil array
  -> (Exp dimS1 -> Exp dim' -> Exp dim1)                                        -- First stencil specification function
  -> (Exp dimS2 -> Exp dim' -> Exp dim2)                                        -- Second stencil specification function
  -> (Exp dim' -> Acc (Array dimS1 a) -> Acc (Array dimS2 b) -> Acc (Scalar c)) -- Populate/update function
  -> Acc (Array dim' c)                                                         -- Destination array
}}}

The semantics of the the above functions is best described through a Repa
implementation of ""stencil"" using the Repa traverse functions:

{{{
-- |Traversal over a single array using stencils (Repa version).
stencil
    :: forall sh sh' shS a b . (Shape sh, Shape sh', Shape shS, Elt a, Elt b)
    => Array sh a
    -> (sh -> sh')
    -> (sh -> shS)
    -> (shS -> sh' -> sh)
    -> (sh' -> Array shS a -> b)
    -> Array sh' b

stencil arr outExtent stencilExtent stencilFunc updateEle
    = traverse
        arr
        outExtent
        updateEle'
    where
        updateEle' _ ix = updateEle ix stencilArr
            where
                stencilArr = traverse arr stencilExtent (\get stencilIx -> get $ stencilFunc stencilIx ix) 


-- |Traversal over two arrays at once using stencils (Repa version).
stencil2
    :: forall sh1 sh2 sh' shS1 shS2 a b c . (Shape sh1, Shape sh2, Shape sh', Shape shS1, Shape shS2, Elt a, Elt b, Elt c)
    => Array sh1 a
    -> Array sh2 b
    -> (sh1 -> sh2 -> sh')
    -> (sh1 -> shS1)
    -> (sh2 -> shS2)
    -> (shS1 -> sh' -> sh1)
    -> (shS2 -> sh' -> sh2)
    -> (sh' -> Array shS1 a -> Array shS2 b -> c)
    -> Array sh' c

stencil2 arrA arrB outExtent stencilExtentA stencilExtentB stencilFuncA stencilFuncB updateEle
    = traverse2
        arrA
        arrB
        outExtent
        updateEle'
    where
        updateEle' _ _ ix = updateEle ix stencilArrA stencilArrB
            where
                stencilArrA = traverse arrA stencilExtentA (\get stencilIx -> get $ stencilFuncA stencilIx ix)
                stencilArrB = traverse arrB stencilExtentB (\get stencilIx -> get $ stencilFuncB stencilIx ix)
}}}

The ""stencil"" primitive has been proven useful in a number of computer vision computation steps. Below
are two supporting examples:

{{{
-- |The y gradient is computed by convolving columns with a [-1, 0, 1] kernel.
-- At image boundaries, i.e. the top and bottom rows, the top or bottom row pixels are
-- used, respectively. That is, pixels outside the image are clamped to the values of
-- the edge pixels.
gradientYCompute :: Repa.Array DIM2 Int -> Repa.Array DIM2 Int
gradientYCompute imgArr
    = Repa.stencil imgArr id (const $ extent kernel) stencilFunc updateEle
    where
        (Z :. h :. w) = extent imgArr

        kernel = force $ Repa.fromList (Z :. (2 :: Int)) ([1, -1] :: [Int])

        stencilFunc (Z :. i) (Z :. y :. x) | i == 0 = (Z :. clamp (y - 1) :. x)
                                           | i == 1 = (Z :. clamp (y + 1) :. x)

        updateEle (Z :. y :. x) stencilArr = Repa.sumAll $ Repa.zipWith (*) kernel stencilArr

        clamp v | v < 0     = 0
                | v > h - 1 = h - 1
                | otherwise = v


-- |Histogram image is computed by binning the orientation gradients (values are between 0
-- and 7) in a 4x4 neighbourhood patch for each pixel in the image. An orientation is only
-- counted if its associated magnitude is over a certain threshold. The resultant array is
-- 3-dimensional where the 3rd dimension's axis is orientation bins. The pixel of interest
-- is located in the lower right corner of the patch.
histogramImageCompute :: Repa.Array DIM2 Int -> Repa.Array DIM2 Int -> Repa.Array DIM3 Int
histogramImageCompute gradMagArr gradOriArr
    = Repa2.stencil2
        gradMagArr
        gradOriArr
        (\(Z :. y :. x) _ -> (Z :. y :. x :. 8))
        (const patchExtent)
        (const patchExtent)
        stencilPatch
        stencilPatch
        updateEle

    where
        patchExtent = (Z :. (4 :: Int) :. (4 :: Int))

        stencilPatch (Z :. j :. i) (Z :. y :. x :. bin) = (Z :. clamp (y + j - 3) :. clamp (x + i - 3))

        updateEle (Z :. y :. x :. bin) magPatch oriPatch
            = sumAll $ Repa.zipWith (\mag ori -> if mag > magnitudeThreshold && ori == bin then 1 else 0) magPatch oriPatch

        clamp x | x < 0     = 0
                | otherwise = x 

        noOrientation = 8       -- Only count orientation bin values betwen 0 and 7
        magnitudeThreshold = 8
}}}
"	feature request	closed	blocker	AutoMap, prototype	Accelerate language	0.7.1.0	fixed