The Haskell image processing benchmark game

Since by now we have several ongoing image processing library projects in Haskell I thought it’s time to spice things up a little. I therefore announce the Haskell image processing benchmark game. Right now the five contestants are: friday, unm-hip, yarr, repa and opencv but the contest is open and participation is just one pull request away! The disciplines are: reading an image, binary thresholding and mean filtering. Of course more disciplines will be added.
Performance is not the only goal one should strive for. Ease of use and number of predefined algorithms are two other important dimensions.
Let’s compare the different libraries with code examples for specific use cases.

The image type

Currently all benchmarks run on grey images so the first thing we need is a type for those.

-- friday
import Vision.Image.Grey (Grey)
type Image = Grey

-- unm-hip
import Data.Image.Boxed (GrayImage)
type Image = GrayImage

-- yarr
import Data.Yarr (UArray,D,L,Dim2)
import Data.Word (Word8)
type Image = UArray D L Dim2 Word8

-- repa
import Data.Array.Repa (Array,D,DIM2)
import Data.Word (Word8)
type Image = Array D DIM2 Word8

-- opencv
import Foreign.ForeignPtr (ForeignPtr)
data IPLImage
type Image = ForeignPtr IPLImage

For friday and unm-hip the preferred type for grey images was easy to find. There is only one. For yarr and repa not so much. In opencv there is only one image type for any kind of image regardless of the pixel type. From experience I can tell that this is not a good thing.

Reading images from disk

Next we want to read in an image. Let’s compare:

-- friday
import Vision.Image.Storage (load)
import Vision.Image.Type (convert)

readPng :: FilePath -> IO Image
readPng filepath = do
    Right image <- load Nothing filepath
    return (convert image)

-- unm-hip
import Data.Image.Boxed (readImage)

readPgm :: FilePath -> IO Image
readPgm = readImage

-- yarr
import Data.Yarr (delay)
import Data.Yarr.IO.Image (readImage)
import qualified Data.Yarr.IO.Image as Yarr (Image(Grey))

readPng :: FilePath -> IO Image
readPng filepath = do
    Yarr.Grey image <- readImage filepath
    return (delay image)

-- repa
import Data.Array.Repa.IO.DevIL (readImage,runIL)
import qualified Data.Array.Repa.IO.DevIL as Repa (Image(Grey))
import Data.Array.Repa (delay)

readPng :: FilePath -> IO Image
readPng filepath = runIL (do
   Repa.Grey image <- readImage filepath
   return (delay image))

-- opencv
import Foreign.C.String (CString,withCString)
import Foreign.Ptr (Ptr,FunPtr)
import Foreign.ForeignPtr (ForeignPtr,newForeignPtr)

foreign import ccall "&finalizeImage" opencv_finalizeImage ::
    FunPtr (Ptr IPLImage -> IO ())

foreign import ccall "readPng" opencv_readPng ::
    CString -> IO (Ptr IPLImage)

addImageFinalizer :: Ptr IPLImage -> IO (ForeignPtr IPLImage)
addImageFinalizer = newForeignPtr opencv_finalizeImage

readPng :: FilePath -> IO Image
readPng filepath = do
    imagePtr <- withCString filepath opencv_readPng
    addImageFinalizer imagePtr

-- opencv C code
void finalizeImage(IplImage* image){
    cvReleaseImage(&image);
}

IplImage* readPng(char* filepath){
    return cvLoadImage(filepath,CV_LOAD_IMAGE_GRAYSCALE);
}

Reading in an image in unm-hip is straight forward. Except that it does not know how to read PNG and can only read PGM. For yarr, repa and friday I had to look a little more stuff up than I’d have liked to. And opencv being a foreign library we need to manually make every image we create garbage collected by adding a finalizer.

Forcing an image

For benchmarking we also need a way to make sure the resulting images are fully evaluated. The friday and opencv image types already guarantee that. For the others we define an extra force function:

-- unm-hip
import Data.Image.Internal (maxIntensity)

force :: Image -> IO Image
force image = do
    maxIntensity image `seq` return image

-- yarr
import Data.Yarr (UArray,F,L,Dim2)
import Data.Yarr.Eval (dComputeP)

force :: Image -> IO (UArray F L Dim2 Word8)
force = dComputeP

-- repa
import Data.Array.Repa (Array,DIM2,computeP,deepSeqArray)
import Data.Array.Repa.Repr.Unboxed (U)
import Data.Word (Word8)

force :: Image -> IO (Array U DIM2 Word8)
force image = do
    forcedImage <- computeP image
    forcedImage `deepSeqArray` return forcedImage

Wait, what kind of hack is that forcing an unm-hip image? I couldn’t find a better way because the constructor of the image type is hidden. Of course this means that the benchmark is pointless, but I decided to include it anyway.

Binary thresholding

Now let’s look at the actual image processing algorithms. The threshold function should accept a grey value image and yield an image where every pixel that has a value above 127 in the original image has value 255 and all others have value 0.

-- friday
import qualified Vision.Image.Threshold as Friday (
    threshold)
import Vision.Image.Threshold (
    ThresholdType(BinaryThreshold))

threshold :: Image -> Image
threshold = Friday.threshold (>127) (BinaryThreshold 0 255)

-- unm-hip
threshold :: Image -> Image
threshold = toBinaryImage (<127)

-- yarr
import Data.Yarr.Flow (dmap)

threshold :: Image -> Image
threshold = dmap (\value -> if value > 127 then 255 else 0)

-- repa
import qualified Data.Array.Repa as Repa (map)

threshold :: Image -> Image
threshold = Repa.map (\value -> if value > 127 then 255 else 0)

-- opencv
foreign import ccall "threshold" opencv_threshold :: Ptr IPLImage -> IO (Ptr IPLImage)

threshold :: Image -> IO Image
threshold = withImage opencv_threshold

withImage :: (Ptr IPLImage -> IO (Ptr IPLImage)) -> Image -> IO Image
withImage f image = withForeignPtr image (\imagePtr -> do
    f imagePtr >>= addImageFinalizer)

-- opencv C code
IplImage* threshold(IplImage* image){
IplImage* output = cvCreateImage(cvGetSize(image),IPL_DEPTH_8U,1);
cvThreshold(image,output,127,255,CV_THRESH_BINARY);
return output;
}

Binary thresholding is a simple map that is straight forward to implement in yarr and repa. The other three friday, unm-hipand opencv provide their own specialized functions for binary thresholding. unm-hips toBinaryImage chooses reasonable default values for you. friday provides its own mini language for creating filtering operations that I had to learn. opencv takes as its fifth parameter an integer that represents the thresholding algorithm to be used and depending on that gives different meaning to the other parameters. I personally find that very confusing.

Mean filter

Now the task is to apply a five by five mean filter:

-- friday
import Vision.Image.Filter (apply,blur,SeparableFilter)

mean :: Image -> Image
mean = apply (blur 2 :: SeparableFilter GreyPixel GreyPixel GreyPixel)

-- unm-hip
import Data.Image.Convolution (convolveRows,convolveCols)

mean :: Image -> Image
mean = fmap (/25) . convolveCols [1,1,1,1,1] . convolveRows [1,1,1,1,1]

-- yarr
import Data.Yarr.Convolution (
    dConvolveLinearDim2WithStaticStencil,
    dim2St,Dim2Stencil)
import Data.Yarr.Utils.FixedVector (N1,N5)

mean :: Image -> IO Image
mean image = do

    rowConvolvedImage <- dComputeP (
        dConvolveLinearDim2WithStaticStencil
            rowMeanStencil (dmap fromIntegral image))

    convolvedImage <- dComputeP (
        dConvolveLinearDim2WithStaticStencil
            colMeanStencil (rowConvolvedImage :: UArray F L Dim2 Int))

    return (dmap (fromIntegral . (`div` 25)) (
        convolvedImage :: UArray F L Dim2 Int))

rowMeanStencil :: Dim2Stencil N5 N1 Int (Int -> Int -> IO Int) Int
rowMeanStencil = [dim2St| 1 1 1 1 1 |]

colMeanStencil :: Dim2Stencil N1 N5 Int (Int -> Int -> IO Int) Int
colMeanStencil = [dim2St| 1
                          1
                          1
                          1
                          1 |]

-- repa
import Data.Array.Repa.Algorithms.Convolve (
    convolveOutP,outClamp)
import Data.Array.Repa (
    Array,DIM2,delay,computeP,Z(Z),(:.)((:.)))
import Data.Array.Repa.Repr.Unboxed (
    U,fromListUnboxed)

mean :: Image -> IO Image
mean image = do
    intImage <- computeP (Repa.map fromIntegral image)
    rowConvolvedImage <- convolveOutP outClamp rowMeanStencil intImage
    convolvedImage <- convolveOutP outClamp colMeanStencil rowConvolvedImage
    grayImage <- computeP (Repa.map (fromIntegral . (`div` 25)) convolvedImage)
    return (delay (grayImage :: Array U DIM2 Word8))

rowMeanStencil :: Array U DIM2 Int
rowMeanStencil = fromListUnboxed (Z:.1:.5) [1,1,1,1,1]

colMeanStencil :: Array U DIM2 Int
colMeanStencil = fromListUnboxed (Z:.5:.1) [1,1,1,1,1]

-- opencv
foreign import ccall "mean" opencv_mean ::
    Ptr IPLImage -> IO (Ptr IPLImage)

mean :: Image -> IO Image
mean = withImage opencv_mean

-- opencv C code
IplImage* mean(IplImage* image){
    IplImage* output = cvCreateImage(cvGetSize(image),IPL_DEPTH_8U,1);
    cvSmooth(image,output,CV_BLUR,5,5,0.0,0.0);
    return output;
}

friday provides it’s own set of functions for filtering. Using unm-hip feels like functional image processing should feel like. Using repa and yarr to implement a five by five mean filter required quite a bit of reading documentation and examples. Again the opencv function cvSmooth accepts as its third parameter an integer representing the smoothing algorithm to be used and changes the meaning of the other parameters or completely ignores them depending on that. This is why I think we need one or more Haskell image processing library, so let’s work together to make it happen.

Oh, right, there are benchmark results as well! Disclaimer: While I have tried my best, I might have used the libraries in the wrong way. Please do correct my mistakes.

EDIT: After using deepseq to force unm-hip images, rewriting the mean function in friday to pointful style and adjusting the yarr and repa functions these are the new results.