用repa计算图像直方图

Question

我正在用JuicyPixels-repa从磁盘加载一个RGB映像。不幸的是，图像的数组表示为Array F DIM3 Word8，其中内部维度是RGB像素。这与现有的repa图像处理算法有点不兼容，因为RGB图像是Array U DIM2 (Word8, Word8, Word8)。

我要计算图像的RGB直方图，我正在搜索一个带有签名的函数：

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

我如何折叠我的3d数组，以获得一个一维阵列为每个颜色通道？

编辑：

主要的问题并不是我不能为每个通道从DIM3转换为DIM2 (切片很容易)。问题是，我必须迭代源映像DIM2或DIM3，并且必须累积到一个不同的Shape、(Z:.256)、和范围的DIM1数组中。因此，我不能使用repa的foldS，因为它减少了一个维度，但具有相同的程度。

我也尝试过traverse，但是它迭代目标图像的范围，提供一个函数从源图像中获取像素，这将导致非常低效率的代码，对每个颜色值计算相同的像素。

一种很好的方法是在Vector上进行简单的折叠，将直方图类型作为累加器，但不幸的是，我没有基于U (未装箱)或基于V (向量)的数组，因此可以有效地获得Vector。我有一个Array F (外接指针)。

Answer 1

好吧，我找到了几分钟。下面，我将介绍四种解决方案，并使最糟糕的解决方案(中间两种，涉及O(n)数据转换)对您来说非常容易。

让我们确认Dumb解决方案

从显而易见的事情开始是合理的。您可以使用Data.List.foldl遍历行和列，从初始零数组构建直方图(未测试/部分代码如下)：

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

我不知道这会有多高的效率，但我怀疑它会做太多的边界检查。切换到unsafeIndex应该是合理的，假设延迟数组hist*做得很好，因为无论您选择实现incElem。

你可以构建你想要的数组

使用traverse，您实际上可以将JP样式数组转换为包含元素元组的DIM2数组：

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))

你能给我指一下你说过的使用这种格式的代码吗？这将是很简单的补丁JP包括这种类型的功能。

您可以构建您提到的无框向量

您提到了一个简单的解决方案，就是折叠未装箱的向量，但遗憾的是JP-repa没有提供无盒数组。幸运的是，转换很简单：

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

我们可以修补Repa

这只是一个问题，因为Repa没有我认为的普通traverse函数。Repa的遍历更多的是一种数组结构，它恰好为另一个数组提供索引函数。我们希望以如下形式进行遍历：

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

但粗的，这实际上只是一个畸形的褶皱。因此，让我们重命名它并重新排序参数：

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

这与(先前存在的) foldAllS操作形成了很好的对比：

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

注意我们的新折叠是如何具有两个关键特性的。结果类型不需要匹配元素类型，所以我们可以从直方图的元组开始。其次，我们的折叠版本传递索引，它允许您选择要更新的元组中的哪个直方图(如果有的话)。

你可以懒洋洋地使用最新的JuicyPixels-Repa

要获得首选的Repa数组格式，或获取未装箱的向量，只需使用新上传的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

现在，您可以折叠矢量或直接使用元组数组，这是您最初希望的。

我还丑陋地尝试了第一个，非常天真的解决方案：

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] ]

我对这段代码的性能太谨慎了(每个索引有3条遍历.我一定很累了)把它扔进JP-Repa，但如果你发现它运行良好，那就让我知道。