Calculate an image histogram with repa

I'm loading an RGB image from disk with JuicyPixels-repa. Unfortunately the Array representation of the image is Array F DIM3 Word8 where the inner dimension is the RGB pixels. That's a bit incompatibe with existing repa imageprocessing algorithms where an RGB image is Array U DIM2 (Word8, Word8, Word8).

I want to calculate the RGB histograms of the image, I'm searching a function with the Signature:

type Hist = Array U DIM1 Int
histogram:: Array F DIM3  Word8 -> (Hist, Hist, Hist)

how can I fold my 3d array to get a 1d array for each colorchannel?

Edit:

The main problem is not that I'm not able to convert from DIM3 to DIM2 for each channel (easy done with slicing). The problem is that I have to iterate the source image DIM2 or DIM3 and have to accumulate to an DIM1 array of a different Shape (Z:.256) and extent. So I can't use repa's foldS as it reduces the dimension by one, but with the same extent.

I also experimented with traverse but it iterates over the extent of the destination image, providing a function to get pixels from the source image, that would lead to very inefficient code, counting the same pixels for each colorvalue.

A good way would be a simple folding over a Vector with the histogram type as accumulator, but unfortunately I have no U (unboxed) or V (vector) based array, from which I can efficiently get a Vector. I have an Array F (foreign pointer).

Ok, I found a few minutes. Below, I cover four solutions and have made the worst solutions (the middle two, involving O(n) data conversion) really easy for you.

Lets Acknowledge the Dumb Solution

It's reasonable to start with the obvious. You could use Data.List.foldl to traverse the rows and columns, building up your histograms from initial zero arrays (untested/partial code follows):

foldl (\(histR, histG, histB) (row,col) ->
            let r = arr ! (Z:.row:.col:.0)
                g = arr ! (Z:.row:.col:.1)
                b = arr ! (Z:.row:.col:.2)
            in (incElem r histR, incElem g histG, incElem b histB)
         (zero,zero,zero)
         [ (row,col) | row <- [0..nrRow-1], col <- [0..nrCol-1] ]
 ...
 where (Z:.nrRow:.nrCol:._) = extent arr

I'm not sure how efficient this will be, but suspect that it will do too much bounds checking. Switching to unsafeIndex should do reasonably, assuming the delayed arrays, hist*, do well due to however you'd pick to implement incElem.

You Can Build the Array You Want

Using traverse you can actually convert JP-Repa style arrays into DIM2 arrays with tuples for elements:

main = do
    let arr = R.fromFunction (Z:.a:.b:.c) (\(Z:.i:.j:.k) -> i+j-k)
        a =4 :: Int
        b = 4 :: Int
        c= 4 :: Int
        new = R.traverse arr
                       (\(Z:.r:.c:._) -> Z:.r:.c) -- the extent
                       (\l idx -> (l (idx:.0)
                                  ,l (idx:.1)
                                  ,l (idx :. 2)))
    print (R.computeS new :: R.Array R.U DIM2 (Int,Int,Int))

Could you point me to the body of code you talked about that uses this format? It would be simple to patch JP-Repa to include a function of this type.

You can build the Unboxed Vector You Mentioned

You mentioned an easy solution is to fold over unboxed vectors, but lamented that JP-repa doesn't provide an unboxed array. Luckily, conversion is simple:

toUnboxed :: Img a -> VU.Vector Word8
toUnboxed = R.toUnboxed . R.computeUnboxedS . R.delay . imgData

We Could Patch Repa

This is really only a problem because Repa doesn't have what I consider a normal traverse function. Repa's traverse is more of an array construction that happens to provide an indexing function into another array. We want traverse in the form:

newTraverse :: Array r sh e -> a -> (a -> sh -> e -> a) -> a

but of coarse this is actually just a malformed fold. So lets rename it and reorder the arguments:

foldAllIdxS :: (sh -> a - > e -> a) -> a -> Array r sh e -> a

which contrasts nicely with the (preexisting) foldAllS operation:

foldAllS :: (a -> a -> a) -> a -> Array r sh a -> a

Notice how our new fold has two critical characteristics. The result type is not required to match the element type, so we could start with a tuple of Histograms. Second, our version of fold passes the index, which allows you to select which histogram in the tuple to update (if any).

You can lazily use the latest JuicyPixels-Repa

To acquire your preferred Repa array format, or to acquire an unboxed vector, you can just use the newly uploaded JuicyPixels-Repa-0.6.

someImg <- readImage path :: IO (Either String (Img RGBA))
let img = either (error "Blah") id someImg
    uvec = toUnboxed img
    tupleArr = collapseColorChannel img

Now you can fold over the vector or use the tuple array directly, as you originally desired.

I also took an ugly stab at fleshing out the first, horribly naive, solution:

histograms :: Img a -> (Histogram, Histogram, Histogram, Histogram)
histograms (Img arr) =
    let (Z:.nrRow:.nrCol:._) = R.extent arr
        zero = R.fromFunction (Z:.256) (\_ -> 0 :: Word8)
        incElem idx x = RU.unsafeTraverse x id (\l i -> l i + if i==(Z:.fromIntegral idx) then 1 else 0)
    in Prelude.foldl (\(hR, hG, hB, hA) (row,col) ->
             let r = R.unsafeIndex arr (Z:.row:.col:.0)
                 g = R.unsafeIndex arr (Z:.row:.col:.1)
                 b = R.unsafeIndex arr (Z:.row:.col:.2)
                 a = R.unsafeIndex arr (Z:.row:.col:.3)
             in (incElem r hR, incElem g hG, incElem b hB, incElem a hA))
          (zero,zero,zero,zero)
          [ (row,col) | row <- [0..nrRow-1], col <- [0..nrCol-1] ]

I'm too wary of the performance of this code (3 traversals per index... I must be tired) to throw it into JP-Repa, but if you find it works well then let me know.

Recommended topics

Hot tags