Blurring Image with convolution matrices - artifacts (Haskell)

It is great that you are learning image processing in Haskell, but there are too many problems in your implementation and I don't want others to be learning from it. I'll first go through the problems and then provide a better solution for doing Gaussian blur in Haskell.

Issue with integer overflow that was already identified in the other answer by @Daniel Wagner is definitely the problem. However I would have solved it by changing the type of pxPlus instead:

pxPlus :: PixelRGBF -> PixelRGBF -> PixelRGBF

but that's besides the point.

Much bigger problem in the original implementation is the performance, it is simply atrocious. Just to give you an idea, implementation in the question takes ~16 sec on my laptop, while the solution I provide towards the bottom takes 0.007 sec

Here are just some of the main reasons:

• !! is O(n) operation for a list, thus making indexing into the kernel ridiculously slow. And because convolution is O(m*n*k^2) complexity already, where m and n are image dimensions and k is the kernel side the actual implementation is of O(m*n*k^4) complexity. That is forgoing other performance limitations of lists like boxed elements and cache unfriendliness.
• blur function is applied to an image that has an alpha channel, which is rightfully ignored, but for some reason it is just discarded and set to 255. It is best to either drop it completely before supplying the image to blur function, or use the transparency value from the original image.
• Gaussian convolution is separable, so instead of applying a 5x5 kernel it is faster to apply 1x5 kernel onto the rows and then 5x1 to the columns. This approach reduces complexity to O(m*n*k)
• Convolution can be implemented much faster in general, with help of various optimizations such as fast border resolution with unsafe internal indexing, parallelization etc. All of which are pretty hard to get right, that's why it is better to use a library that provides such functionality.
• Gaussian blur kernel values that were used are designed for 8bit channels, but you are using them as floating point values. The point of 8bit in the old days was to make it run as fast as possible on old hardware because floating point numbers were slow. That's why all those tutorials online use such crappy kernels. Using a mix of Word8 and Double incorrectly like it's done in the question results in a much worse quality for the filter than it really has to be. In other words you get an error due to bad values in the filter and then that error is further amplified during summation, because rounding should happen after summing the floating point values, not before.
• Another very common mistake is that the filter is applied to the non-linear sRGB color space, which is most likely the color space the source image is encoded in. It requires an inverse gamma correction first
• Padding with zeros results in an image that has this white frame and of course working with slightly bigger image adds a bit of overhead.

JuicyPixels is a great package for encoding/decoding images, but it is definitely not suitable for more complex image processing tasks. What you need is either an actual image processing library or an array manipulation library that has already a convolution algorithm implemented. I have a library HIP that is designed precisely for this sort of stuff, but it is undergoing a major rewrite, so I will provide an answer with an array library massiv and massiv-io which actually use JuicyPixels underneath for reading/writing images.

Here is the fast and accurate way to do gaussian blur:

import Data.Massiv.Array as A
import Data.Massiv.Array.Unsafe (makeUnsafeConvolutionStencil)
import Data.Massiv.Array.IO as A

blurImageF :: (ColorModel cs Float) => Image S cs Float -> Image S cs Float
blurImageF = computeAs S . applyStencil padding blurStencil5x1f
           . computeAs S . applyStencil padding blurStencil1x5f
  where
    padding = noPadding -- decides what happens at the border
{-# INLINE blurImageF #-}

blurStencil1x5f :: (Floating e, ColorModel cs e) => Stencil Ix2 (Pixel cs e) (Pixel cs e)
blurStencil1x5f = makeUnsafeConvolutionStencil (Sz2 1 5) (0 :. 2) stencil
  where
    stencil f = f (0 :. -2) 0.03467403390152031 .
                f (0 :. -1) 0.23896796340399287 .
                f (0 :.  0) 0.45271600538897480 .
                f (0 :.  1) 0.23896796340399287 .
                f (0 :.  2) 0.03467403390152031
    {-# INLINE stencil #-}
{-# INLINE blurStencil1x5f #-}

blurStencil5x1f :: (Floating e, ColorModel cs e) => Stencil Ix2 (Pixel cs e) (Pixel cs e)
blurStencil5x1f = makeUnsafeConvolutionStencil (Sz2 5 1) (2 :. 0) stencil
  where
    stencil f = f (-2 :. 0) 0.03467403390152031 .
                f (-1 :. 0) 0.23896796340399287 .
                f ( 0 :. 0) 0.45271600538897480 .
                f ( 1 :. 0) 0.23896796340399287 .
                f ( 2 :. 0) 0.03467403390152031
    {-# INLINE stencil #-}
{-# INLINE blurStencil5x1f #-}

Values in the kernel were computed using very accurate integral approximation and optimal standard deviation for the kernel size of σ=5/6

In order to see this in action we first read the image, which is automatically converted to linear sRGB colorspace with Float precision (FYI the original image that was supplied in the question was in Y'CbCr color space, not sRGB with alpha channel) Then we apply the blur to it and concat together with cropped original (applying convolution without padding reduces the size). After which we convert the constructed image to Y'CbCr color space and write it to JPEG file:

λ> img <- readImageAuto "4ZYKa.jpg" :: IO (Image S (SRGB 'Linear) Float)
λ> let imgBlurred = blurImageF img
λ> displayImage imgBlurred -- this will open the blurred image with default viewer
λ> imgCropped <- extractM (2 :. 2) (size imgBlurred) img
λ> imgBoth <- appendM 1 imgCropped imgBlurred
λ> let out = convertPixel <$> imgBoth :: Image DL (Y'CbCr SRGB) Word8
λ> writeImage "out.jpg" out

Note that this will be quite a bit slower in GHCi. It needs to be compiled with -O2 -threaded -with-rtsopts=-N flags to get the optimal performance.

You write

| p1 + p2 > 255 = 255
| p1 + p2 < 0 = 0

to protect against overflow, but these conditions are never true. Word8 values are always between 0 and 255. A correct way to check for overflow looks like this:

| p1 + p2 < p1 = 255

But better would be to avoid overflow in the first place. I note that your blur matrix's coefficients add up to 256/255. Perhaps you could start by fixing that.