Below is an example of creating a Stencil, which, when mapped over a 2-dimensional array, will compute an average of all elements in a 3x3 square for each element in that array.

Note - Make sure to add an INLINE pragma, otherwise performance will be terrible.

average3x3Stencil :: (Default a, Fractional a) => Stencil Ix2 a a
average3x3Stencil = makeStencil (Sz (3 :. 3)) (1 :. 1) $ \ get ->
  (  get (-1 :. -1) + get (-1 :. 0) + get (-1 :. 1) +
     get ( 0 :. -1) + get ( 0 :. 0) + get ( 0 :. 1) +
     get ( 1 :. -1) + get ( 1 :. 0) + get ( 1 :. 1)   ) / 9
{-# INLINE average3x3Stencil #-}

In order to see the affect of padding we can simply apply an identity stencil to an array:

>>> import Data.Massiv.Array as A
>>> a = computeAs P $ resize' (Sz2 2 3) (Ix1 1 ... 6)
>>> applyStencil noPadding idStencil a
Array DW Seq (Sz (2 :. 3))
  [ [ 1, 2, 3 ]
  , [ 4, 5, 6 ]
  ]
>>> applyStencil (Padding (Sz2 1 2) (Sz2 3 4) (Fill 0)) idStencil a
Array DW Seq (Sz (6 :. 9))
  [ [ 0, 0, 0, 0, 0, 0, 0, 0, 0 ]
  , [ 0, 0, 1, 2, 3, 0, 0, 0, 0 ]
  , [ 0, 0, 4, 5, 6, 0, 0, 0, 0 ]
  , [ 0, 0, 0, 0, 0, 0, 0, 0, 0 ]
  , [ 0, 0, 0, 0, 0, 0, 0, 0, 0 ]
  , [ 0, 0, 0, 0, 0, 0, 0, 0, 0 ]
  ]

It is also a nice technique to see the border resolution strategy in action:

>>> applyStencil (Padding (Sz2 2 3) (Sz2 2 3) Wrap) idStencil a
Array DW Seq (Sz (6 :. 9))
  [ [ 1, 2, 3, 1, 2, 3, 1, 2, 3 ]
  , [ 4, 5, 6, 4, 5, 6, 4, 5, 6 ]
  , [ 1, 2, 3, 1, 2, 3, 1, 2, 3 ]
  , [ 4, 5, 6, 4, 5, 6, 4, 5, 6 ]
  , [ 1, 2, 3, 1, 2, 3, 1, 2, 3 ]
  , [ 4, 5, 6, 4, 5, 6, 4, 5, 6 ]
  ]
>>> applyStencil (Padding (Sz2 2 3) (Sz2 2 3) Edge) idStencil a
Array DW Seq (Sz (6 :. 9))
  [ [ 1, 1, 1, 1, 2, 3, 3, 3, 3 ]
  , [ 1, 1, 1, 1, 2, 3, 3, 3, 3 ]
  , [ 1, 1, 1, 1, 2, 3, 3, 3, 3 ]
  , [ 4, 4, 4, 4, 5, 6, 6, 6, 6 ]
  , [ 4, 4, 4, 4, 5, 6, 6, 6, 6 ]
  , [ 4, 4, 4, 4, 5, 6, 6, 6, 6 ]
  ]
>>> applyStencil (Padding (Sz2 2 3) (Sz2 2 3) Reflect) idStencil a
Array DW Seq (Sz (6 :. 9))
  [ [ 6, 5, 4, 4, 5, 6, 6, 5, 4 ]
  , [ 3, 2, 1, 1, 2, 3, 3, 2, 1 ]
  , [ 3, 2, 1, 1, 2, 3, 3, 2, 1 ]
  , [ 6, 5, 4, 4, 5, 6, 6, 5, 4 ]
  , [ 6, 5, 4, 4, 5, 6, 6, 5, 4 ]
  , [ 3, 2, 1, 1, 2, 3, 3, 2, 1 ]
  ]
>>> applyStencil (Padding (Sz2 2 3) (Sz2 2 3) Continue) idStencil a
Array DW Seq (Sz (6 :. 9))
  [ [ 1, 3, 2, 1, 2, 3, 2, 1, 3 ]
  , [ 4, 6, 5, 4, 5, 6, 5, 4, 6 ]
  , [ 1, 3, 2, 1, 2, 3, 2, 1, 3 ]
  , [ 4, 6, 5, 4, 5, 6, 5, 4, 6 ]
  , [ 1, 3, 2, 1, 2, 3, 2, 1, 3 ]
  , [ 4, 6, 5, 4, 5, 6, 5, 4, 6 ]
  ]

>>> import Data.Massiv.Array as A
>>> a = computeAs P $ iterateN (Sz2 2 5) (* 2) (1 :: Int)
>>> a
Array P Seq (Sz (2 :. 5))
  [ [ 2, 4, 8, 16, 32 ]
  , [ 64, 128, 256, 512, 1024 ]
  ]
>>> applyStencil noPadding (sumStencil (Sz2 1 2)) a
Array DW Seq (Sz (2 :. 4))
  [ [ 6, 12, 24, 48 ]
  , [ 192, 384, 768, 1536 ]
  ]
>>> [2 + 4, 4 + 8, 8 + 16, 16 + 32] :: [Int]
[6,12,24,48]

>>> import Data.Massiv.Array as A
>>> a = computeAs P $ iterateN (Sz2 2 2) (+1) (0 :: Int)
>>> a
Array P Seq (Sz (2 :. 2))
  [ [ 1, 2 ]
  , [ 3, 4 ]
  ]
>>> applyStencil (Padding 0 2 (Fill 0)) (productStencil 2) a
Array DW Seq (Sz (3 :. 3))
  [ [ 24, 0, 0 ]
  , [ 0, 0, 0 ]
  , [ 0, 0, 0 ]
  ]
>>> applyStencil (Padding 0 2 Reflect) (productStencil 2) a
Array DW Seq (Sz (3 :. 3))
  [ [ 24, 64, 24 ]
  , [ 144, 256, 144 ]
  , [ 24, 64, 24 ]
  ]

>>> import Data.Massiv.Array as A
>>> a = computeAs P $ iterateN (Sz2 3 4) (+1) (10 :: Double)
>>> a
Array P Seq (Sz (3 :. 4))
  [ [ 11.0, 12.0, 13.0, 14.0 ]
  , [ 15.0, 16.0, 17.0, 18.0 ]
  , [ 19.0, 20.0, 21.0, 22.0 ]
  ]
>>> applyStencil noPadding (avgStencil (Sz2 2 3)) a
Array DW Seq (Sz (2 :. 2))
  [ [ 14.0, 15.0 ]
  , [ 18.0, 19.0 ]
  ]
>>> Prelude.sum [11.0, 12.0, 13.0, 15.0, 16.0, 17.0] / 6 :: Double
14.0

Here is a sample implementation of max pooling.

>>> import Data.Massiv.Array as A
>>> a <- computeAs P <$> resizeM (Sz2 9 9) (Ix1 10 ..: 91)
>>> a
Array P Seq (Sz (9 :. 9))
  [ [ 10, 11, 12, 13, 14, 15, 16, 17, 18 ]
  , [ 19, 20, 21, 22, 23, 24, 25, 26, 27 ]
  , [ 28, 29, 30, 31, 32, 33, 34, 35, 36 ]
  , [ 37, 38, 39, 40, 41, 42, 43, 44, 45 ]
  , [ 46, 47, 48, 49, 50, 51, 52, 53, 54 ]
  , [ 55, 56, 57, 58, 59, 60, 61, 62, 63 ]
  , [ 64, 65, 66, 67, 68, 69, 70, 71, 72 ]
  , [ 73, 74, 75, 76, 77, 78, 79, 80, 81 ]
  , [ 82, 83, 84, 85, 86, 87, 88, 89, 90 ]
  ]
>>> computeWithStrideAs P (Stride 3) $ mapStencil Edge (maxStencil (Sz 3)) a
Array P Seq (Sz (3 :. 3))
  [ [ 30, 33, 36 ]
  , [ 57, 60, 63 ]
  , [ 84, 87, 90 ]
  ]

Here is a sample implementation of min pooling.

>>> import Data.Massiv.Array as A
>>> a <- computeAs P <$> resizeM (Sz2 9 9) (Ix1 10 ..: 91)
>>> a
Array P Seq (Sz (9 :. 9))
  [ [ 10, 11, 12, 13, 14, 15, 16, 17, 18 ]
  , [ 19, 20, 21, 22, 23, 24, 25, 26, 27 ]
  , [ 28, 29, 30, 31, 32, 33, 34, 35, 36 ]
  , [ 37, 38, 39, 40, 41, 42, 43, 44, 45 ]
  , [ 46, 47, 48, 49, 50, 51, 52, 53, 54 ]
  , [ 55, 56, 57, 58, 59, 60, 61, 62, 63 ]
  , [ 64, 65, 66, 67, 68, 69, 70, 71, 72 ]
  , [ 73, 74, 75, 76, 77, 78, 79, 80, 81 ]
  , [ 82, 83, 84, 85, 86, 87, 88, 89, 90 ]
  ]
>>> computeWithStrideAs P (Stride 3) $ mapStencil Edge (minStencil (Sz 3)) a
Array P Seq (Sz (3 :. 3))
  [ [ 10, 13, 16 ]
  , [ 37, 40, 43 ]
  , [ 64, 67, 70 ]
  ]

>>> import Data.Massiv.Array as A
>>> a = computeAs P $ iterateN (Sz2 3 4) (+1) (10 :: Int)
>>> a
Array P Seq (Sz (3 :. 4))
  [ [ 11, 12, 13, 14 ]
  , [ 15, 16, 17, 18 ]
  , [ 19, 20, 21, 22 ]
  ]
>>> applyStencil noPadding (foldlStencil (flip (:)) [] (Sz2 3 2)) a
Array DW Seq (Sz (1 :. 3))
  [ [ [20,19,16,15,12,11], [21,20,17,16,13,12], [22,21,18,17,14,13] ]
  ]
>>> applyStencil (Padding (Sz2 1 0) 0 (Fill 10)) (foldlStencil (flip (:)) [] (Sz2 3 2)) a
Array DW Seq (Sz (2 :. 3))
  [ [ [16,15,12,11,10,10], [17,16,13,12,10,10], [18,17,14,13,10,10] ]
  , [ [20,19,16,15,12,11], [21,20,17,16,13,12], [22,21,18,17,14,13] ]
  ]

>>> import Data.Massiv.Array as A
>>> a = computeAs P $ iterateN (Sz2 3 4) (+1) (10 :: Int)
>>> a
Array P Seq (Sz (3 :. 4))
  [ [ 11, 12, 13, 14 ]
  , [ 15, 16, 17, 18 ]
  , [ 19, 20, 21, 22 ]
  ]
>>> applyStencil noPadding (foldrStencil (:) [] (Sz2 2 3)) a
Array DW Seq (Sz (2 :. 2))
  [ [ [11,12,13,15,16,17], [12,13,14,16,17,18] ]
  , [ [15,16,17,19,20,21], [16,17,18,20,21,22] ]
  ]

Here is how to create a 2D horizontal Sobel Stencil:

sobelX :: Num e => Stencil Ix2 e e
sobelX = makeConvolutionStencil (Sz2 3 3) (1 :. 1) $
           \f -> f (-1 :. -1) (-1) . f (-1 :. 1) 1 .
                 f ( 0 :. -1) (-2) . f ( 0 :. 1) 2 .
                 f ( 1 :. -1) (-1) . f ( 1 :. 1) 1
{-# INLINE sobelX #-}