# Ticket #15 (closed feature request: fixed)

## Support for stencil-based traversal in the language front-end

Reported by: | blever | Owned by: | chak |
---|---|---|---|

Priority: | blocker | Milestone: | AutoMap, prototype |

Component: | Accelerate language | Version: | 0.7.1.0 |

Keywords: | Cc: |

### Description

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

## Change History

**Note:**See TracTickets for help on using tickets.