1.4.3 形态学：对象计数
形态学（或数学形态学）是度量和分析基本形状的图像处理方法的基本框架与集合。
形态学通常用于处理二值图像，但是也能够用于灰度图像。
二值图像是指图像的每个像素只能取两个值，通常是 0 和 1。
二值图像通常是，在计算物体的数目，或者度量其大小时，对一幅图像进行阈值化后的结果。
你可以从 http://en.wikipedia.org/wiki/Mathematical_morphology 大体了解形态学及其处理图像的方式。
scipy.ndimage 中 的 morphology 模 块 可 以 实 现 形 态 学 操 作。
你 可 以 使 用 scipy.ndimage 中的 measurements 模块来实现二值图像的计数和度量功能。
下面通过一个简单的例子介绍如何使用它们。
考虑在图 1-12a1 里的二值图像，计算该图像中的对象个数可以通过下面的脚本实现：
from scipy.ndimage import measurements,morphology # 载入图像，然后使用阈值化操作，以保证处理的图像为二值图像 im = array(Image.open('houses.png').convert('L')) im = 1*(im<128) labels, nbr_objects = measurements.label(im) print "Number of objects:", nbr_objects
上面的脚本首先载入该图像，通过阈值化方式来确保该图像是二值图像。通过和 1相乘，脚本将布尔数组转换成二进制表示。然后，我们使用 label() 函数寻找单个的物体，并且按照它们属于哪个对象将整数标签给像素赋值。
图 1-12b 是 labels 数组的图像。
图像的灰度值表示对象的标签。可以看到，在一些对象之间有一些小的连接。进行二进制开（binary open）操作，我们可以将其移除：
# 形态学开操作更好地分离各个对象 im_open = morphology.binary_opening(im,ones((9,5)),iterations=2) labels_open, nbr_objects_open = measurements.label(im_open) print("Number of objects:", nbr_objects_open)
binary_opening() 函数的第二个参数指定一个数组结构元素。该数组表示以一个像素为中心时，使用哪些相邻像素。在这种情况下，我们在 y 方向上使用 9 个像素（上面 4 个像素、像素本身、下面 4 个像素），在 x 方向上使用 5 个像素。你可以指定任意数组为结构元素，数组中的非零元素决定使用哪些相邻像素。参数iterations 决定执行该操作的次数。你可以尝试使用不同的迭代次数 iterations 值，看一下对象的数目如何变化。你可以在图 1-12c 与图 1-12d 中查看经过开操作后的图像，以及相应的标签图像。正如你想象的一样，binary_closing() 函数实现相反的操作。我们将该函数和在 morphology 和 measurements 模块中的其他函数的用法留作练习。
你可以从 scipy.ndimage 模块文档http://docs.scipy.org/doc/scipy/reference/ndimage.html 中了解关于这些函数的更多知识。
图 1-12：形态学示例。使用二值开操作将对象分开，然后计算物体的数目：（a）为原始二值图像；（b）为对应原始图像的标签图像，其中灰度值表示物体的标签；（c）为使用开操作后的二值图像；（d）为开操作后图像的标签图像
1.4.4 一些有用的SciPy模块
SciPy 中包含一些用于输入和输出的实用模块。下面介绍其中两个模块：
- io
- misc
1. 读写.mat文件
如果你有一些数据，或者在网上下载到一些有趣的数据集，这些数据以 Matlab的 .mat 文件格式存储，那么可以使用 scipy.io 模块进行读取。
data = scipy.io.loadmat('test.mat')
上面代码中，data 对象包含一个字典，字典中的键对应于保存在原始 .mat 文件中的变量名。
由于这些变量是数组格式的，因此可以很方便地保存到 .mat 文件中。你仅需创建一个字典（其中要包含你想要保存的所有变量）,然后使用 savemat() 函数：
data = {} data['x'] = x scipy.io.savemat('test.mat',data)
因为上面的脚本保存的是数组 x，所以当读入到 Matlab 中时，变量的名字仍为 x。
关 于 scipy.io 模 块 的 更 多 内 容， 请 参 见 在 线 文 档 http://docs.scipy.org/doc/scipy/reference/io.html
2. 以图像形式保存数组
因为我们需要对图像进行操作，并且需要使用数组对象来做运算，所以将数组直接保存为图像文件 1 非常有用。本书中的很多图像都是这样的创建的。
imsave() 函数 可以从 scipy.misc 模块中载入。
要将数组 im 保存到文件中，可以使用下面的命令：
from scipy.misc import imsave imsave('test.jpg',im)
scipy.misc 模块同样包含了著名的 Lena 测试图像：
lena = scipy.misc.lena()
该脚本返回一个 512×512 的灰度图像数组。
1.5 高级示例：图像去噪
我们通过一个非常实用的例子——图像的去噪——来结束本章。图像去噪是在去除图像噪声的同时，尽可能地保留图像细节和结构的处理技术。我们这里使用 ROF（Rudin-Osher-Fatemi）去噪模型。该模型最早出现在文献 [28] 中。图像去噪对于很多应用来说都非常重要；这些应用范围很广，小到让你的假期照片看起来更漂亮，大到提高卫星图像的质量。ROF 模型具有很好的性质：使处理后的图像更平滑，同时保持图像边缘和结构信息。
注 1：所有 Pylab 图均可保存为多种图像格式，方法是点击图像窗口中的“保存”按钮。
ROF 模型的数学基础和处理技巧非常高深，不在本书讲述范围之内。在讲述如何基于 Chambolle 提出的算法 [5] 实现 ROF 求解器之前，本书首先简要介绍一下 ROF模型。
一幅（灰度）图像 I 的全变差（Total Variation，TV）定义为梯度范数之和。
其中范数 || 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 部分。纯化论者和经验丰富的程序员们知道这些区别，他们能够选择自己喜欢的方式。但是，为了使本书的内容和例子更容易被读者接受，我们不打算这样做。请读者注意。
