Comonad在图像处理中的应用

分享 未结
1 1 1 390
小编 2020-04-21发布
收藏 点赞
来源: https://zhuanlan.zhihu.com/p/133460864

前几天我回答了一个关于comonad的问题Monad和Comonad到底是什么东西?。其中有讲到comonad的应用例子,但都还不够直观和实用。后来找到一个Comonad在图像处理中的应用的例子,觉得不错,在这里重新整理一下写出来。与大家分享一下这个例子。

我们已经知道,Comonad最通常的用法是作为一个无穷的stream,可以源源不断的通过扩展复制生成新的无穷的stream。只要提供从stream中计算得到值的cokleisli函数,我们就可以得到新的stream,或者新的有限的列表,或者一个累加的值。

例如下面就是生成无穷的自然数的stream,和对自然数的stream进行滑动窗口求和的Comonad的应用例子。

data Stream a = a :< Stream a deriving Show
instance Functor Stream where
  fmap f (x :< xs) = f x :< fmap f xs
instance Comonad Stream where
  extract (x :< _) = x
  duplicate s@(_:<xs) = s :< duplicate xs

generate :: (a -> a) -> a -> Stream a
generate f x = x :< generate f x =>> (f . extract)

nats :: Stream Integer
nats = generate (+1) 0

getWindowSum :: Integer -> Stream Integer -> Integer
getWindowSum n (x :< xs)
  | n == 1 = x
  | n >  1 = x + getWindowSum (n-1) xs
  | otherwise = undefined

tensWindows :: Stream Integer
tensWindows = nats =>> getWindowSum 10

这里nats调用generate生成了一个有无穷多个自然数的stream,其类型是Stream Integer。

nats生成自然数的stream的计算步骤效果如下所示,逐行累加起来就得到了自然数的stream:

0 0 0 0 0 0 0 0 ...
  1 1 1 1 1 1 1 1 ...
    1 1 1 1 1 1 1 1 ...
      1 1 1 1 1 1 1 1 ...
        1 1 1 1 1 1 1 1 ...
          1 1 1 1 1 1 1 1 ...
            1 1 1 1 1 1 1 1 ...

nats的值如下所示:

> nats
0 :< 1 :< 2 :< 3 :< 4 :< 5 :< 6 :< 7 :< ...

而getWindowSum函数则是取出自然数的stream的前10个自然数,求这10个自然数的和,得到一个新的无穷的stream。这和信号处理中的滑动滤波的窗口计算是一样的,其值如下所示:

> tensWindows
45 :< 55 :< 65 :< 75 :< 85 :< 95 :< 105 :< ...

我们现在来看图像处理,图像可以看成一个有限的二维阵列的数据结构,如下面的图所示,并不是一个无穷的stream。图像处理中用的最多的一种场景是图像滤波,根据当前图像的一个点的邻域像素的值计算出新图像上该点的值,这种场景是可以应用Comonad来计算的。

在图像滤波的计算中,我们关心的是图像上的一个点及其邻域的点,因此有一个焦点,如下图所示:

于是,我们可以把图像看成一个带焦点的二维阵列的数据结构。其中二维阵列可以用带宽和高参数的一维向量来表示,再加上焦点坐标就得到了我们要的图像数据结构定义。如下所示:

-- | 用一维向量表示的图像二维阵列
data BoxedImage a = BoxedImage
  { biWidth  :: !Int
  , biHeight :: !Int
  , biData   :: !(V.Vector a)
  }

instance Functor BoxedImage where
  fmap f (BoxedImage w h d) = BoxedImage w h (fmap f d)

-- | 带焦点的图像数据结构
data FocusedImage a = FocusedImage 
  { boxedImage :: !(BoxedImage a)
  , cx         :: !Int
  , cy         :: !Int
  }

instance Functor FocusedImage where
  fmap f (FocusedImage bi x y) = FocusedImage (fmap f bi) x y

-- | 带焦点的图像和一般图像的转换函数
focus :: BoxedImage a -> FocusedImage a
focus bi
  | biWidth bi > 0 && biHeight bi > 0 = FocusedImage bi 0 0
  | otherwise                         = error "Can not focus on empty image"

unfocus :: FocusedImage a -> BoxedImage a
unfocus (FocusedImage bi _ _) = bi

带焦点的图像数据结构FocusedImage a 我们可以将其看成一个Comonad。其extract函数就是直接把焦点处的值取出来。duplicate函数就是把这个图像的二维阵列中的每一个点的值都扩展为一个以该点为焦点的带焦点的图像数据结构,这个新的图像数据结构的二维阵列和原来的是一样的,保持了相似性。extend函数的参数f 是个函数,这个函数f 将新的带焦点的图像数据结构reduce为一个图像的像素值。这里我们回顾一下 extend f = fmap f . duplicate。

instance Comonad FocusedImage where
  extract (FocusedImage bi x y) = (biData bi) V.! (biWidth bi * y + x)
  duplicate (FocusedImage bi@(BoxedImage w h _) x y) = FocusedImage
      (BoxedImage w h $ V.generate (w * h) $ \i ->
          let (y', x') = i `quotRem` w
          in  FocusedImage bi x' y')
      x y

如果extend函数的参数f 这个函数只关心图像中的焦点及其邻域的点的值(一般邻域是3 x 3,5 x 5),则extend f就是这个带焦点的图像数据结构的滤波处理函数,得到一个新的带焦点的图像数据结构,这个新的图像就是滤波处理后的图像。

例如一个常见的gauss模糊滤波的处理,其计算核大小是3 x 3的,如下所示:

我们就可以定义如下的计算核函数,可以看到其只关心3 x 3 大小的邻域的点的值:

blur :: Integral a => V.Vector a -> a
blur v = fromIntegral
       $ (`quot` 16)
       $ ( nbi 0     + nbi 1 * 2 + nbi 2
         + nbi 3 * 2 + nbi 4 * 4 + nbi 5 * 2
         + nbi 6     + nbi 7 * 2 + nbi 8
         )
  where {-# INLINE nbi #-}
        nbi :: Int -> Int
        nbi i = fromIntegral $ v V.! i

然后定义如下的获取邻域点的值的函数,注意这个neighbour函数对超出图像边界的邻域点做了特殊处理,使用了回绕的处理方式,这是一种效果较好的处理方式:

neighbour :: Int -> Int -> FocusedImage a -> a
neighbour dx dy (FocusedImage bi x y) = (biData bi) V.! (biWidth bi * y' + x')
  where x' = wrap (x + dx) 0 (biWidth bi - 1)
        y' = wrap (y + dy) 0 (biHeight bi - 1)
        {-# INLINE wrap #-}
        wrap i lo hi
          | i < lo = lo - i            -- 坐标回绕处理
          | i > hi = hi - (i - hi)     -- 坐标回绕处理
          | otherwise = i

再定义一个取遍焦点的邻域的所有点的值的函数,把这几个函数组合起来,就得到了带焦点的图像数据结构的高斯滤波处理函数blurImage:

filterImage :: Integral a => (V.Vector a -> a) -> Int -> FocusedImage a -> a
filterImage filter kernelW fimg@(FocusedImage bi x y) = filter $
    V.generate kernelSize $ \i -> neighbour (nbx i) (nby i) fimg
  where nbx i = (-(kernelW `quot` 2)) + (i `rem`  kernelW);
        nby i = (-(kernelW `quot` 2)) + (i `quot` kernelW);
        {-# INLINE kernelSize #-}
        kernelSize = kernelW * kernelW

{-# INLINE blurImage #-}
blurImage :: Integral a => FocusedImage a -> a
blurImage = filterImage blur 3

然后在main函数中调用blurImage函数,就可以对图像进行高斯模糊处理了,这里的readImage读取jpg 文件并将其转换为灰度图的带焦点的图像数据结构FocusedImage a,writePng则将带焦点的图像数据结构FocusedImage a 写入到png 文件中。这两个函数用了JuicyPixels包中提供的readImage和writePng两个函数,JuicyPixels是非常优秀的加载、保存、转换图像的的Haskell实现。使用它,减少了自己写加载和保存图像的函数的麻烦。

main :: IO ()
main = do
  putStrLn "Hello, Haskell!"
  img <- readImage "data/yingwuhua.jpg"
  -- writePng "output.png" . unfocus . (=>> ebossImage) . focus $ img
  writePng "output.png" . unfocus . (=>> blurImage) . focus $ img

文章的题图经过这个高斯模糊滤波处理后的灰度图效果如下:

我们还可以定义其他滤波处理函数,比如浮雕效果的滤波处理。我们先定义计算核函数,计算核大小同样是3 x 3的:

eboss :: Integral a => V.Vector a -> a
eboss v = fromIntegral
       $ (+ 127)
       $ ( nbi 0 * 3 `quot` 2 + 0  + 0
         + 0                  + 0  + 0
         + 0                  + 0  + nbi 8 * (-3) `div` 2
         )
  where {-# INLINE nbi #-}
        nbi :: Int -> Int
        nbi i = fromIntegral $ v V.! i

重复利用前面定义的取遍焦点的邻域的所有点的函数,就得到了浮雕效果的滤波处理函数:

{-# INLINE ebossImage #-}
ebossImage :: Integral a => FocusedImage a -> a
ebossImage = filterImage eboss 3

这里由于Haskell的惰性计算特性,eboss计算核函数中的0 项不参与计算,计算速度比高斯模糊滤波处理快了60%左右,其灰度图的效果图如下:

Comonad的cokleisli函数是可以组合的,因此我们也可以把这些带焦点的图像的滤波处理函数组合起来,得到一个多种滤波组合的功能更丰富的滤波处理函数:

reduceNoise :: Integral a => FocusedImage a -> a
reduceNoise fimg =
  let !original   = extract fimg
      !blurred    = blurImage fimg
      !edged       = fromIntegral original - fromIntegral blurred :: Int
      !threshold  = if edged < 7 && edged < (-7) then 0 else edged
  in fromIntegral $ fromIntegral blurred + threshold

由于上面的滤波处理函数的计算核是一样大小的,其组合起来所耗费的滤波处理时间比单个的滤波处理时间差别很小。

这是Comonad的一个很有意思的应用例子,大家可以在此基础上做些进一步的尝试。

参考和引用:

Image Processing with Comonads


jaspervdj.be



回帖