Converting images to greyscale using JuicyPixel

I would like to use Haskell to do some image processing so I have been writing small programs to test performance. I wrote a program to batch convert PNG images from color to greyscale. I compared the performance of the program below with a Python script using PIL. I found when compiled with -threaded -O2, the Haskell code takes four times longer than the python script (this comparison was done on a folder of 50 images).

Is there a way to speed this process up?

Haskell Code

import Codec.Picture
import System.Environment
import System.Directory
import System.FilePath
import Control.Monad

dynToImg :: DynamicImage -> Image PixelRGB8
dynToImg (ImageRGB8 img) = img

procAvg :: FilePath -> FilePath -> IO ()
procAvg inpath outpath = do
  Right inImage ) baseNames
      outfiles  = map (outfolder ) baseNames

  zipWithM_ procAvg infiles outfiles

avgImage :: Image PixelRGB8 -> Image Pixel8
avgImage = pixelMap pixelAvg

pixelAvg :: PixelRGB8 -> Pixel8
pixelAvg (PixelRGB8 r g b) =
  (r `div` 3 + g `div` 3 + b `div` 3 +
          (r `mod` 3 + g `mod` 3 + b `mod` 3) `div` 3)

Python Code

from PIL import Image
import sys
import os

if __name__ == "__main__":
    infolder = sys.argv[1]
    outfolder = sys.argv[2]

    fnames = os.listdir(infolder)

    for f in fnames:
        infile = os.path.join(infolder,f)
        outfile = os.path.join(outfolder,f)

        inimg = Image.open(infile)
        gray = inimg.convert('L')
        gray.save(outfile)
        inimg.close()

Python

What does your Python code look like? I'm sure that most of
the work is actually being done in C, not Python. JuicyPixels
allows you to write an arbitrary pixel mapping function in Haskell.
If the Python library called a user-supplied Python function to convert
every pixel I can assure you that it wouldn't be as fast as it is.

Other Libraries

That said, JuicyPixels may not be the fastest Haskell library to use.
For instance, Haskell can call routines in C libraries, and there are
Haskell bindings to the ImageMagick library via the imagemagick
package.

Another package you should have a look at is friday.
On assorted algorithms it compares favorable with ImageMagick,
It cannot produce SIMD instructions, so it's not as fast
as libraries such as OpenCV. On the other hand, it gives you the
tools to construct your own image transformations so it's a lot more
versitile than a library which gives you just a handful of canned
algorithms.

Improving pixelAvg

Obviously the function pixelAvg is the hot spot in the program, so
anything we can do to optimize it will have a significant impact on
the execution time. Here are some versions I tried:

pixelAvg (PixelRGB8 r g b) = 

-- (1) original version:
   (r `div` 3 + g `div` 3 + b `div` 3) + (r `mod` 3 + g `mod` 3 + b `mod` 3) `div` 3

-- (2) use quot and rem:
   (quot r 3 + quot g 3 + quot b 3) + (quot (rem r 3 + rem g 3 + rem b 3) 3)

-- (3) use quotRem:
   let (rq,rr) = quotRem r 3
       (gq,gr) = quotRem g 3
       (bq,br) = quotRem b 3
   in (rq+gq+bq + (rr+gr+br) `quot` 3)

-- (4) simplify the expression
   (div r 3 + div g 3 + div b 3)

-- (5) just use one division
   fromIntegral $ ((fromIntegral r + fromIntegral g + fromIntegral g) :: Int) `quot` 3

-- (6) use a lookup table
   V.unsafeIndex avgVector (fromIntegral r + fromIntegral g + fromIntegral b)

Here are some timings I obtained on a 4246 x 2856 image:

(1) 0m2.796s
(2) 0m2.768s
(3) 0m2.440s
(4) 0m2.297s
(5) 0m2.198s
(6) 0m2.117s

quot/rem avoids a sign check and therefore should be marginally faster than div/mod
while giving the same answers for positive arguments. However an even bigger speed-up
can be obtained by computing both the quotient and remainder with one function call
(either quotRem or divMod) as demonstrated by (3).

Note that the formula used by (4) only differs from the original formula by at most 2
and has a lot fewer operations. Formula (5) reduces the number of divisions to
just one by converting all three pixel values to Ints, adding and then dividing the
total by 3.

As to be expected, the lookup table method is is the overall fastest. Here's the rest
of the code that initializes the table:

import qualified Data.Vector.Unboxed as V
...
avgVector :: V.Vector Pixel8
avgVector = V.generate (3*256) (\x -> fromIntegral $ div x 3)

So there are significant gains to be made by optimizing pixelAvg, but there are
also limitations to what you can achieve. To get better results you'll likely
have to use a different library.

StackExchange Code Review Q#128335, answer score: 4