Python计算机视觉编程——1.5 　高级示例：图像去噪

我们通过一个非常实用的例子——图像的去噪——来结束本章。图像去噪是在去除图像噪声的同时，尽可能地保留图像细节和结构的处理技术。我们这里使用ROF （Rudin-Osher-Fatemi）去噪模型。该模型最早出现在文献[28] 中。图像去噪对于很多应用来说都非常重要；这些应用范围很广，小到让你的假期照片看起来更漂亮， 大到提高卫星图像的质量。ROF 模型具有很好的性质：使处理后的图像更平滑，同时保持图像边缘和结构信息。 ROF 模型的数学基础和处理技巧非常高深，不在本书讲述范围之内。在讲述如何基于Chambolle 提出的算法[5] 实现ROF 求解器之前，本书首先简要介绍一下ROF 模型。 一幅（灰度）图像I 的全变差（Total Variation，TV）定义为梯度范数之和。在连续表示的情况下，全变差表示为： () dJII x d= # （1.1） 在离散表示的情况下，全变差表示为： () JII x d =/ 其中，上面的式子是在所有图像坐标x=[x, y] 上取和。 在Chambolle 提出的ROF 模型里，目标函数为寻找降噪后的图像U，使下式最小： || || 2 (), min JI U UU 2 m - + 其中范数|| ||I U - 是去噪后图像U 和原始图像I 差异的度量。也就是说，本质上该模型使去噪后的图像像素值“平坦”变化，但是在图像区域的边缘上，允许去噪后的图像像素值“跳跃”变化。 按照论文[5] 中的算法，我们可以按照下面的代码实现ROF 模型去噪： from numpy import * def denoise(im,U_init,tolerance=0.1,tau=0.125,tv_weight=100): """ 使用A. Chambolle（2005）在公式（11）中的计算步骤实现Rudin-Osher-Fatemi（ROF）去噪模型 输入：含有噪声的输入图像（灰度图像）、U 的初始值、TV 正则项权值、步长、停业条件 输出：去噪和去除纹理后的图像、纹理残留""" m,n = im.shape # 噪声图像的大小 # 初始化 U = U_init Px = im # 对偶域的x 分量 Py = im # 对偶域的y 分量 error = 1 while (error > tolerance): Uold = U # 原始变量的梯度 GradUx = roll(U,-1,axis=1)-U # 变量U 梯度的x 分量 GradUy = roll(U,-1,axis=0)-U # 变量U 梯度的y 分量 # 更新对偶变量 PxNew = Px + (tau/tv_weight)*GradUx PyNew = Py + (tau/tv_weight)*GradUy NormNew = maximum(1,sqrt(PxNew**2+PyNew**2)) Px = PxNew/NormNew # 更新x 分量（对偶） Py = PyNew/NormNew # 更新y 分量（对偶） # 更新原始变量 RxPx = roll(Px,1,axis=1) # 对x 分量进行向右x 轴平移 RyPy = roll(Py,1,axis=0) # 对y 分量进行向右y 轴平移 DivP = (Px-RxPx)+(Py-RyPy) # 对偶域的散度 U = im + tv_weight*DivP # 更新原始变量 # 更新误差 error = linalg.norm(U-Uold)/sqrt(n*m); return U,im-U # 去噪后的图像和纹理残余 在这个例子中，我们使用了roll() 函数。顾名思义，在一个坐标轴上，它循环“滚动”数组中的元素值。该函数可以非常方便地计算邻域元素的差异，比如这里的导数。我们还使用了linalg.norm() 函数，该函数可以衡量两个数组间（这个例子中是指图像矩阵U 和Uold）的差异。我们将这个denoise() 函数保存到rof.py 文 件中。 下面使用一个合成的噪声图像示例来说明如何使用该函数： from numpy import * from numpy import random from scipy.ndimage import filters import rof # 使用噪声创建合成图像 im = zeros((500,500)) im[100:400,100:400] = 128 im[200:300,200:300] = 255 im = im + 30*random.standard_normal((500,500)) U,T = rof.denoise(im,im) G = filters.gaussian_filter(im,10) # 保存生成结果 from scipy.misc import imsave imsave('synth_rof.pdf',U) imsave('synth_gaussian.pdf',G) 原始图像和图像的去噪结果如图1-13 所示。正如你所看到的，ROF 算法去噪后的图像很好地保留了图像的边缘信息。 图1-13：使用ROF 模型对合成图像去噪：（a）为原始噪声图像；（b）为经过高斯模糊的图像 （σ=10）；（c）为经过ROF 模型去噪后的图像 下面看一下在实际图像中使用ROF 模型去噪的效果： from PIL import Image from pylab import * import rof im = array(Image.open('empire.jpg').convert('L')) U,T = rof.denoise(im,im) figure() gray() imshow(U) axis('equal') axis('off') show() 经过ROF 去噪后的图像如图1-14c 所示。为了方便比较，该图中同样显示了模糊后的图像。可以看到，ROF 去噪后的图像保留了边缘和图像的结构信息，同时模糊了“噪声”。 图1-14：使用ROF 模型对灰度图像去噪：（a）为原始噪声图像；（b）为经过高斯模糊的图像（σ=5）；（c）为经过ROF 模型去噪后的图像 练习 (1) 如图1-9 所示，将一幅图像进行高斯模糊处理。随着σ 的增加，绘制出图像轮廓。在你绘制出的图中，图像的轮廓有何变化？你能解释为什么会发生这些变化吗？ (2) 通过将图像模糊化，然后从原始图像中减去模糊图像，来实现反锐化图像掩模操作（http://en.wikipedia.org/wiki/Unsharp_masking）。反锐化图像掩模操作可以实现图像锐化效果。试在彩色和灰度图像上使用反锐化图像掩模操作，观察该操作的效果。 (3) 除了直方图均衡化，商图像是另一种图像归一化的方法。商图像可以通过除以模糊后的图像I/(I *Gσ) 获得。尝试使用该方法，并使用一些样本图像进行验证。 (4) 使用图像梯度，编写一个在图像中获得简单物体（例如，白色背景中的正方形） 轮廓的函数。 (5) 使用梯度方向和大小检测图像中的线段。估计线段的长度以及线段的参数，并在原始图像中重新绘制该线段。 (6) 使用label() 函数处理二值化图像，并使用直方图和标签图像绘制图像中物体的大小分布。 (7) 使用形态学操作处理阈值化图像。在发现一些参数能够产生好的结果后，使用morphology 模块里面的center_of_mass() 函数寻找每个物体的中心坐标，将其在图像中绘制出来。 代码示例约定 从第2 章起，我们假定PIL、NumPy 和Matplotlib 都包括在你所创建的每个文件和每个代码例子的开头： from PIL import Image from numpy import * from pylab import * 这种约定使得示例代码更清晰，同时也便于读者理解。除此之外，我们使用SciPy 模块时，将会在代码示例中显式声明。 一些纯化论者会反对这种将全体模块导入的方式，坚持如下使用方式： import numpy as np import matplotlib.pyplot as plt 这种方式能够保持命名空间（知道每个函数从哪儿来）。因为我们不需要PyLab 中的NumPy 部分，所以该例子只从Matplotlib 中导入pyplot 部分。纯化论者和经验丰富的程序员们知道这些区别，他们能够选择自己喜欢的方式。但是，为了使本书的内容和例子更容易被读者接受，我们不打算这样做。 请读者注意。