使用卷积矩阵模糊图像 - 伪影 (Haskell)
Blurring Image with convolution matrices - artifacts (Haskell)
我正在尝试使用 Haskell(使用 JuicyPixels)进行一些图像处理。我已经完成了将高斯模糊应用于图像的这个功能,并且在处理后的图像中有一些奇怪的伪像:
blur :: Image PixelRGBA8 -> Image PixelRGBA8
blur img@Image {..} = generateImage blurrer imageWidth imageHeight
where blurrer x y | x >= (imageWidth - offset) || x < offset
|| y >= (imageHeight - offset) || y < offset = whitePx
| otherwise = do
let applyKernel i j p | j >= matrixLength = applyKernel (i + 1) 0 p
| i >= matrixLength = p
| otherwise = do
let outPixel = pxMultNum
(pixelAt img (x + j - offset) (y + i - offset))
(gblurMatrix !! i !! j)
applyKernel i (j+1) (outPixel `pxPlus` p)
applyKernel 0 0 zeroPx
gblurMatrix = [[1 / 255, 4 / 255, 6 / 255, 4 / 255, 1 / 255],
[4 / 255, 16 / 255, 24 / 255, 16 / 255, 4 / 255],
[6 / 255, 24 / 255, 36 / 255, 24 / 255, 6 / 255],
[4 / 255, 16 / 255, 24 / 255, 16 / 255, 4 / 255],
[1 / 255, 4 / 255, 6 / 255, 4 / 255, 1 / 255]]
matrixLength = length gblurMatrix
offset = matrixLength `div` 2
whitePx = PixelRGBA8 255 255 255 255
zeroPx = PixelRGBA8 0 0 0 255
pxPlus :: PixelRGBA8 -> PixelRGBA8 -> PixelRGBA8
pxPlus (PixelRGBA8 r1 g1 b1 a1) (PixelRGBA8 r2 g2 b2 a2) = PixelRGBA8 (s r1 r2) (s g1 g2) (s b1 b2) 255
where s p1 p2 | p1 + p2 > 255 = 255
| p1 + p2 < 0 = 0
| otherwise = p1 + p2
pxMultNum :: PixelRGBA8 -> Double -> PixelRGBA8
pxMultNum (PixelRGBA8 r g b a) q = PixelRGBA8 (m r) (m g) (m b) (m a)
where m px = pxify $ fromIntegral px * q -- pxify is just synonym to function rounding a = (floor (a + 0.5))
这里是输入图片
and output image
高斯模糊矩阵中 1/256
系数的蓝色区域看起来有点少,但如何完全去除它们?
修复:我猜这是因为 Pixel8
类型 0-255 和 Int,但我只是通过将内核中心的 36
更改为 35
来修复整个问题,所以总和矩阵系数是(我相信)255/255 = 1.
你写
| p1 + p2 > 255 = 255
| p1 + p2 < 0 = 0
防止溢出,但这些条件永远不会成立。 Word8
值 总是 在 0
和 255
之间。检查溢出的正确方法如下所示:
| p1 + p2 < p1 = 255
但最好首先避免溢出。我注意到您的模糊矩阵的系数加起来为 256/255
。也许你可以先解决这个问题。
很高兴你在Haskell学习图像处理,但是你的实现问题太多,我不希望其他人学习它。我将首先解决问题,然后在 Haskell.
中提供更好的高斯模糊解决方案
@Daniel Wagner 在另一个答案中已经确定的整数溢出问题绝对是问题所在。但是我会通过更改 pxPlus
的类型来解决它:
pxPlus :: PixelRGBF -> PixelRGBF -> PixelRGBF
但这不是重点。
原来的实现中更大的问题是性能,简直太糟糕了。只是为了给你一个想法,这个问题的实现在我的笔记本电脑上需要大约 16 秒,而我在底部提供的解决方案需要 0.007 秒
以下是一些主要原因:
!!
是列表的 O(n)
操作,因此使索引到内核的速度慢得离谱。因为卷积已经是 O(m*n*k^2)
复杂度,其中 m
和 n
是图像维度,k
是内核端,实际实现是 O(m*n*k^4)
复杂度。那就是放弃列表的其他性能限制,例如盒装元素和缓存不友好。
- blur 函数应用于具有 alpha 通道的图像,该通道理所当然地被忽略,但由于某种原因它被丢弃并设置为
255
。最好在将图像提供给模糊函数之前将其完全删除,或者使用原始图像的透明度值。
- 高斯卷积是可分离的,因此与其应用 5x5 内核,不如将
1x5
内核应用于行,然后将 5x1 应用于列更快。这种方法将复杂性降低到 O(m*n*k)
- 在各种优化的帮助下,卷积通常可以更快地实现,例如使用不安全的内部索引的快速边界解析、并行化等。所有这些都很难做到正确,这就是为什么最好使用提供此类功能的库。
- 使用的高斯模糊内核值是为 8 位通道设计的,但您将它们用作浮点值。过去 8 位的要点是使它 运行 在旧硬件上尽可能快,因为浮点数很慢。这就是为什么所有这些在线教程都使用如此糟糕的内核。像在问题中那样错误地混合使用
Word8
和 Double
会导致过滤器的质量比实际需要的质量差得多。换句话说,由于过滤器中的错误值,您会得到一个错误,然后在求和期间该错误会进一步放大,因为舍入应该发生在对浮点值求和之后,而不是之前。
- 另一个很常见的错误是 the filter is applied to the non-linear sRGB color space, which is most likely the color space the source image is encoded in. It requires an inverse gamma correction first
- 用零填充会产生具有此白框的图像,当然,使用稍大的图像会增加一些开销。
JuicyPixels
是 encoding/decoding 图像的一个很好的包,但它绝对不适合更复杂的图像处理任务。您需要的是一个实际的图像处理库或一个已经实现了卷积算法的数组处理库。我有一个库 HIP
that is designed precisely for this sort of stuff, but it is undergoing a major rewrite, so I will provide an answer with an array library massiv
and massiv-io
,它实际上在下面使用 JuicyPixels
用于 reading/writing 图像。
这是进行高斯模糊的快速准确的方法:
import Data.Massiv.Array as A
import Data.Massiv.Array.Unsafe (makeUnsafeConvolutionStencil)
import Data.Massiv.Array.IO as A
blurImageF :: (ColorModel cs Float) => Image S cs Float -> Image S cs Float
blurImageF = computeAs S . applyStencil padding blurStencil5x1f
. computeAs S . applyStencil padding blurStencil1x5f
where
padding = noPadding -- decides what happens at the border
{-# INLINE blurImageF #-}
blurStencil1x5f :: (Floating e, ColorModel cs e) => Stencil Ix2 (Pixel cs e) (Pixel cs e)
blurStencil1x5f = makeUnsafeConvolutionStencil (Sz2 1 5) (0 :. 2) stencil
where
stencil f = f (0 :. -2) 0.03467403390152031 .
f (0 :. -1) 0.23896796340399287 .
f (0 :. 0) 0.45271600538897480 .
f (0 :. 1) 0.23896796340399287 .
f (0 :. 2) 0.03467403390152031
{-# INLINE stencil #-}
{-# INLINE blurStencil1x5f #-}
blurStencil5x1f :: (Floating e, ColorModel cs e) => Stencil Ix2 (Pixel cs e) (Pixel cs e)
blurStencil5x1f = makeUnsafeConvolutionStencil (Sz2 5 1) (2 :. 0) stencil
where
stencil f = f (-2 :. 0) 0.03467403390152031 .
f (-1 :. 0) 0.23896796340399287 .
f ( 0 :. 0) 0.45271600538897480 .
f ( 1 :. 0) 0.23896796340399287 .
f ( 2 :. 0) 0.03467403390152031
{-# INLINE stencil #-}
{-# INLINE blurStencil5x1f #-}
内核中的值是使用非常精确的积分近似值和 σ=5/6
内核大小的最佳标准偏差计算得出的
为了在实际操作中看到这一点,我们首先读取图像,该图像会自动转换为 线性 sRGB 颜色space,精度为 Float
(仅供参考问题中提供的原始图像是 Y'CbCr
颜色 space,而不是带 alpha 通道的 sRGB
)然后我们对其应用模糊并与裁剪后的原始图像连接在一起(应用卷积而不填充减小了尺寸)。之后我们将构建的图像转换为 Y'CbCr
颜色 space 并将其写入 JPEG 文件:
λ> img <- readImageAuto "4ZYKa.jpg" :: IO (Image S (SRGB 'Linear) Float)
λ> let imgBlurred = blurImageF img
λ> displayImage imgBlurred -- this will open the blurred image with default viewer
λ> imgCropped <- extractM (2 :. 2) (size imgBlurred) img
λ> imgBoth <- appendM 1 imgCropped imgBlurred
λ> let out = convertPixel <$> imgBoth :: Image DL (Y'CbCr SRGB) Word8
λ> writeImage "out.jpg" out
请注意,这在 GHCi 中会慢很多。它需要使用 -O2 -threaded -with-rtsopts=-N
标志进行编译以获得最佳性能。
我正在尝试使用 Haskell(使用 JuicyPixels)进行一些图像处理。我已经完成了将高斯模糊应用于图像的这个功能,并且在处理后的图像中有一些奇怪的伪像:
blur :: Image PixelRGBA8 -> Image PixelRGBA8
blur img@Image {..} = generateImage blurrer imageWidth imageHeight
where blurrer x y | x >= (imageWidth - offset) || x < offset
|| y >= (imageHeight - offset) || y < offset = whitePx
| otherwise = do
let applyKernel i j p | j >= matrixLength = applyKernel (i + 1) 0 p
| i >= matrixLength = p
| otherwise = do
let outPixel = pxMultNum
(pixelAt img (x + j - offset) (y + i - offset))
(gblurMatrix !! i !! j)
applyKernel i (j+1) (outPixel `pxPlus` p)
applyKernel 0 0 zeroPx
gblurMatrix = [[1 / 255, 4 / 255, 6 / 255, 4 / 255, 1 / 255],
[4 / 255, 16 / 255, 24 / 255, 16 / 255, 4 / 255],
[6 / 255, 24 / 255, 36 / 255, 24 / 255, 6 / 255],
[4 / 255, 16 / 255, 24 / 255, 16 / 255, 4 / 255],
[1 / 255, 4 / 255, 6 / 255, 4 / 255, 1 / 255]]
matrixLength = length gblurMatrix
offset = matrixLength `div` 2
whitePx = PixelRGBA8 255 255 255 255
zeroPx = PixelRGBA8 0 0 0 255
pxPlus :: PixelRGBA8 -> PixelRGBA8 -> PixelRGBA8
pxPlus (PixelRGBA8 r1 g1 b1 a1) (PixelRGBA8 r2 g2 b2 a2) = PixelRGBA8 (s r1 r2) (s g1 g2) (s b1 b2) 255
where s p1 p2 | p1 + p2 > 255 = 255
| p1 + p2 < 0 = 0
| otherwise = p1 + p2
pxMultNum :: PixelRGBA8 -> Double -> PixelRGBA8
pxMultNum (PixelRGBA8 r g b a) q = PixelRGBA8 (m r) (m g) (m b) (m a)
where m px = pxify $ fromIntegral px * q -- pxify is just synonym to function rounding a = (floor (a + 0.5))
这里是输入图片
and output image
高斯模糊矩阵中 1/256
系数的蓝色区域看起来有点少,但如何完全去除它们?
修复:我猜这是因为 Pixel8
类型 0-255 和 Int,但我只是通过将内核中心的 36
更改为 35
来修复整个问题,所以总和矩阵系数是(我相信)255/255 = 1.
你写
| p1 + p2 > 255 = 255
| p1 + p2 < 0 = 0
防止溢出,但这些条件永远不会成立。 Word8
值 总是 在 0
和 255
之间。检查溢出的正确方法如下所示:
| p1 + p2 < p1 = 255
但最好首先避免溢出。我注意到您的模糊矩阵的系数加起来为 256/255
。也许你可以先解决这个问题。
很高兴你在Haskell学习图像处理,但是你的实现问题太多,我不希望其他人学习它。我将首先解决问题,然后在 Haskell.
中提供更好的高斯模糊解决方案@Daniel Wagner 在另一个答案中已经确定的整数溢出问题绝对是问题所在。但是我会通过更改 pxPlus
的类型来解决它:
pxPlus :: PixelRGBF -> PixelRGBF -> PixelRGBF
但这不是重点。
原来的实现中更大的问题是性能,简直太糟糕了。只是为了给你一个想法,这个问题的实现在我的笔记本电脑上需要大约 16 秒,而我在底部提供的解决方案需要 0.007 秒
以下是一些主要原因:
!!
是列表的O(n)
操作,因此使索引到内核的速度慢得离谱。因为卷积已经是O(m*n*k^2)
复杂度,其中m
和n
是图像维度,k
是内核端,实际实现是O(m*n*k^4)
复杂度。那就是放弃列表的其他性能限制,例如盒装元素和缓存不友好。- blur 函数应用于具有 alpha 通道的图像,该通道理所当然地被忽略,但由于某种原因它被丢弃并设置为
255
。最好在将图像提供给模糊函数之前将其完全删除,或者使用原始图像的透明度值。 - 高斯卷积是可分离的,因此与其应用 5x5 内核,不如将
1x5
内核应用于行,然后将 5x1 应用于列更快。这种方法将复杂性降低到O(m*n*k)
- 在各种优化的帮助下,卷积通常可以更快地实现,例如使用不安全的内部索引的快速边界解析、并行化等。所有这些都很难做到正确,这就是为什么最好使用提供此类功能的库。
- 使用的高斯模糊内核值是为 8 位通道设计的,但您将它们用作浮点值。过去 8 位的要点是使它 运行 在旧硬件上尽可能快,因为浮点数很慢。这就是为什么所有这些在线教程都使用如此糟糕的内核。像在问题中那样错误地混合使用
Word8
和Double
会导致过滤器的质量比实际需要的质量差得多。换句话说,由于过滤器中的错误值,您会得到一个错误,然后在求和期间该错误会进一步放大,因为舍入应该发生在对浮点值求和之后,而不是之前。 - 另一个很常见的错误是 the filter is applied to the non-linear sRGB color space, which is most likely the color space the source image is encoded in. It requires an inverse gamma correction first
- 用零填充会产生具有此白框的图像,当然,使用稍大的图像会增加一些开销。
JuicyPixels
是 encoding/decoding 图像的一个很好的包,但它绝对不适合更复杂的图像处理任务。您需要的是一个实际的图像处理库或一个已经实现了卷积算法的数组处理库。我有一个库 HIP
that is designed precisely for this sort of stuff, but it is undergoing a major rewrite, so I will provide an answer with an array library massiv
and massiv-io
,它实际上在下面使用 JuicyPixels
用于 reading/writing 图像。
这是进行高斯模糊的快速准确的方法:
import Data.Massiv.Array as A
import Data.Massiv.Array.Unsafe (makeUnsafeConvolutionStencil)
import Data.Massiv.Array.IO as A
blurImageF :: (ColorModel cs Float) => Image S cs Float -> Image S cs Float
blurImageF = computeAs S . applyStencil padding blurStencil5x1f
. computeAs S . applyStencil padding blurStencil1x5f
where
padding = noPadding -- decides what happens at the border
{-# INLINE blurImageF #-}
blurStencil1x5f :: (Floating e, ColorModel cs e) => Stencil Ix2 (Pixel cs e) (Pixel cs e)
blurStencil1x5f = makeUnsafeConvolutionStencil (Sz2 1 5) (0 :. 2) stencil
where
stencil f = f (0 :. -2) 0.03467403390152031 .
f (0 :. -1) 0.23896796340399287 .
f (0 :. 0) 0.45271600538897480 .
f (0 :. 1) 0.23896796340399287 .
f (0 :. 2) 0.03467403390152031
{-# INLINE stencil #-}
{-# INLINE blurStencil1x5f #-}
blurStencil5x1f :: (Floating e, ColorModel cs e) => Stencil Ix2 (Pixel cs e) (Pixel cs e)
blurStencil5x1f = makeUnsafeConvolutionStencil (Sz2 5 1) (2 :. 0) stencil
where
stencil f = f (-2 :. 0) 0.03467403390152031 .
f (-1 :. 0) 0.23896796340399287 .
f ( 0 :. 0) 0.45271600538897480 .
f ( 1 :. 0) 0.23896796340399287 .
f ( 2 :. 0) 0.03467403390152031
{-# INLINE stencil #-}
{-# INLINE blurStencil5x1f #-}
内核中的值是使用非常精确的积分近似值和 σ=5/6
为了在实际操作中看到这一点,我们首先读取图像,该图像会自动转换为 线性 sRGB 颜色space,精度为 Float
(仅供参考问题中提供的原始图像是 Y'CbCr
颜色 space,而不是带 alpha 通道的 sRGB
)然后我们对其应用模糊并与裁剪后的原始图像连接在一起(应用卷积而不填充减小了尺寸)。之后我们将构建的图像转换为 Y'CbCr
颜色 space 并将其写入 JPEG 文件:
λ> img <- readImageAuto "4ZYKa.jpg" :: IO (Image S (SRGB 'Linear) Float)
λ> let imgBlurred = blurImageF img
λ> displayImage imgBlurred -- this will open the blurred image with default viewer
λ> imgCropped <- extractM (2 :. 2) (size imgBlurred) img
λ> imgBoth <- appendM 1 imgCropped imgBlurred
λ> let out = convertPixel <$> imgBoth :: Image DL (Y'CbCr SRGB) Word8
λ> writeImage "out.jpg" out
请注意,这在 GHCi 中会慢很多。它需要使用 -O2 -threaded -with-rtsopts=-N
标志进行编译以获得最佳性能。