Ticket #15 (closed feature request: fixed)

Opened 4 years ago

Last modified 4 years ago

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

Changed 4 years ago by chak

  • milestone changed from milestone1 to Automap, prototype

Changed 4 years ago by chak

Those stencil functions certainly make sense. However, I am wondering whether they should be slightly more specialised. In both of the two use cases that you give, the extent of the stencil is constant (and not dependent on the extent of the source array(s)). Do you have any examples, where this is not the case? Otherwise, I would propose to use the following signature instead (here just for the unary version):

-- |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 dimS                                           -- Extent of the 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

This slightly specialised version would enable much more efficient code. For example, on a GPU, you would like to keep the stencil array entirely in shared memory. That is not hard if the stencil size is constant (and not too big), but would be virtually impossible if the stencil size is truly variable.

Changed 4 years ago by chak

  • status changed from new to assigned

Changed 4 years ago by chak

  • status changed from assigned to closed
  • resolution set to fixed

I implemented the following:

-- |Map a stencil over an array.  In contrast to 'map', the domain of a stencil function is an
--  entire /neighbourhood/ of each array element.  Neighbourhoods are sub-arrays centred around a
--  focal point.  They are not necessarily rectangular, but they are symmetric in each dimension
--  and have an extent of at least three in each dimensions — due to the symmetry requirement, the
--  extent is necessarily odd.  The focal point is the array position that is determined by the
--  stencil.
--
--  For those array positions where the neighbourhood extends past the boundaries of the source
--  array, a boundary condition determines the contents of the out-of-bounds neighbourhood
--  positions.
--
stencil :: (Ix dim, Elem a, Elem b, Stencil dim a stencil)
        => (stencil -> Exp b)          -- ^stencil function
        -> Boundary a                  -- ^boundary condition
        -> Acc (Array dim a)           -- ^source array
        -> Acc (Array dim b)           -- ^destination array
Note: See TracTickets for help on using tickets.