Python图像处理库(2)

1.4　SciPy

SciPy（http://scipy.org/）是建立在 NumPy 基础上，用于数值运算的开源工具包。SciPy 提供很多高效的操作，可以实现数值积分、优化、统计、信号处理，以及对我们来说最重要的图像处理功能。接下来，本节会介绍 SciPy 中大量有用的模块。SciPy 是个开源工具包，可以从http://scipy.org/Download 下载。

1.4.1　图像模糊

图像的高斯模糊是非常经典的图像卷积例子。本质上，图像模糊就是将（灰度）图像 I 和一个高斯核进行卷积操作：

Iσ = I*Gσ

其中 * 表示卷积操作；Gσ 是标准差为 σ 的二维高斯核，定义为 :

$G_sigma=frac{1}{2pisigma}e^{-(x^2+y^2)/2sigma^2}$

高斯模糊通常是其他图像处理操作的一部分，比如图像插值操作、兴趣点计算以及很多其他应用。

SciPy 有用来做滤波操作的 scipy.ndimage.filters 模块。该模块使用快速一维分离的方式来计算卷积。你可以像下面这样来使用它：

from PIL import Image
from numpy import *
from scipy.ndimage import filters

im = array(Image.open('empire.jpg').convert('L'))
im2 = filters.gaussian_filter(im,5)

上面 guassian_filter() 函数的最后一个参数表示标准差。

图 1-9 显示了随着 σ 的增加，一幅图像被模糊的程度。σ 越大，处理后的图像细节丢失越多。如果打算模糊一幅彩色图像，只需简单地对每一个颜色通道进行高斯模糊：

im = array(Image.open('empire.jpg'))
im2 = zeros(im.shape)
for i in range(3):
  im2[:,:,i] = filters.gaussian_filter(im[:,:,i],5)
im2 = uint8(im2)

在上面的脚本中，最后并不总是需要将图像转换成 uint8 格式，这里只是将像素值用八位来表示。我们也可以使用：

im2 = array(im2,'uint8')

来完成转换。

关于该模块更多的内容以及不同参数的选择，请查看http://docs.scipy.org/doc/scipy/reference/ndimage.html 上 SciPy 文档中的 scipy.ndimage部分。

图 1-9：使用 scipy.ndimage.filters 模块进行高斯模糊：（a）原始灰度图像；（b）使用 σ=2 的高斯滤波器；（c）使用 σ=5 的高斯滤波器；（d）使用 σ=10 的高斯滤波器

1.4.2　图像导数

整本书中可以看到，在很多应用中图像强度的变化情况是非常重要的信息。强度的变化可以用灰度图像 I（对于彩色图像，通常对每个颜色通道分别计算导数）的 x 和 y 方向导数 Ix 和 Iy 进行描述。

图像的梯度向量为∇I = [Ix, Iy]T。梯度有两个重要的属性，一是梯度的大小：

$left|oldsymbol{ abla I} ight|=sqrt{{oldsymbol{I}_x}^2+{oldsymbol{I}_y}^2}$

它描述了图像强度变化的强弱，一是梯度的角度：

α=arctan2(Iy, Ix)

描述了图像中在每个点（像素）上强度变化最大的方向。NumPy 中的 arctan2() 函数返回弧度表示的有符号角度，角度的变化区间为 -π...π。

我们可以用离散近似的方式来计算图像的导数。图像导数大多数可以通过卷积简单地实现：

Ix=I*Dx 和 Iy=I*Dy

对于 Dx 和 Dy，通常选择 Prewitt 滤波器：

$D_x=egin{vmatrix}-1&0&1\-1&0&1\ -1&0&1end{vmatrix}$ 和 $D_y=egin{vmatrix}-1&-1&-1\0&0&0\ 1&1&1end{vmatrix}$

或者 Sobel 滤波器：

$D_x=egin{vmatrix}-1&0&1\-2&0&2\ -1&0&1end{vmatrix}$ 和 $D_y=egin{vmatrix}-1&-2&-1\0&0&0\ 1&2&1end{vmatrix}$

这些导数滤波器可以使用 scipy.ndimage.filters 模块的标准卷积操作来简单地实现，例如：

from PIL import Image
from numpy import *
from scipy.ndimage import filters

im = array(Image.open('empire.jpg').convert('L'))

# Sobel 导数滤波器
imx = zeros(im.shape)
filters.sobel(im,1,imx)

imy = zeros(im.shape)
filters.sobel(im,0,imy)

magnitude = sqrt(imx**2+imy**2)

上面的脚本使用 Sobel 滤波器来计算 x 和 y 的方向导数，以及梯度大小。sobel() 函数的第二个参数表示选择 x 或者 y 方向导数，第三个参数保存输出的变量。图 1-10 显示了用 Sobel 滤波器计算出的导数图像。在两个导数图像中，正导数显示为亮的像素，负导数显示为暗的像素。灰色区域表示导数的值接近于零。

图 1-10：使用 Sobel 导数滤波器计算导数图像：（a）原始灰度图像；（b）x 导数图像；（c）y导数图像；（d）梯度大小图像

上述计算图像导数的方法有一些缺陷：在该方法中，滤波器的尺度需要随着图像分辨率的变化而变化。为了在图像噪声方面更稳健，以及在任意尺度上计算导数，我们可以使用高斯导数滤波器：

Ix=I*Gσx 和 Iy=I*Gσy

其中，Gσx和 Gσy 表示 Gσ 在 x 和 y 方向上的导数，Gσ 为标准差为 σ 的高斯函数。

我们之前用于模糊的 filters.gaussian_filter() 函数可以接受额外的参数，用来计算高斯导数。可以简单地按照下面的方式来处理：

sigma = 5 # 标准差

imx = zeros(im.shape)
filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)

imy = zeros(im.shape)
filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)

该函数的第三个参数指定对每个方向计算哪种类型的导数，第二个参数为使用的标准差。你可以查看相应文档了解详情。图 1-11 显示了不同尺度下的导数图像和梯度大小。你可以和图 1-9 中做相同尺度模糊的图像做比较。

图 1-11：使用高斯导数计算图像导数：x 导数图像（上），y 导数图像（中），以及梯度大小图像（下）；（a）为原始灰度图像，（b）为使用 σ=2 的高斯导数滤波器处理后的图像，（c）为使用 σ=5 的高斯导数滤波器处理后的图像，（d）为使用 σ=10 的高斯导数滤波器处理后的图像

1.4.3　形态学：对象计数

形态学（或数学形态学）是度量和分析基本形状的图像处理方法的基本框架与集合。形态学通常用于处理二值图像，但是也能够用于灰度图像。二值图像是指图像的每个像素只能取两个值，通常是 0 和 1。二值图像通常是，在计算物体的数目，或者度量其大小时，对一幅图像进行阈值化后的结果。你可以从 http://en.wikipedia.org/wiki/Mathematical_morphology 大体了解形态学及其处理图像的方式。

scipy.ndimage 中的 morphology 模块可以实现形态学操作。你可以使用 scipy.ndimage 中的measurements 模块来实现二值图像的计数和度量功能。下面通过一个简单的例子介绍如何使用它们。

考虑在图 1-12a3 里的二值图像，计算该图像中的对象个数可以通过下面的脚本实现：

3这个图像实际上是图像“分割”后的结果。如果你想知道该图像是如何创建的，可以查看 9.3 节。

from scipy.ndimage import measurements,morphology

# 载入图像，然后使用阈值化操作，以保证处理的图像为二值图像
im = array(Image.open('houses.png').convert('L'))
im = 1*(im<</span>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。

读写.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。
以图像形式保存数组

因为我们需要对图像进行操作，并且需要使用数组对象来做运算，所以将数组直接保存为图像文件 4 非常有用。本书中的很多图像都是这样的创建的。

imsave() 函数可以从 scipy.misc 模块中载入。要将数组 im 保存到文件中，可以使用下面的命令：
```
from scipy.misc import imsave
imsave('test.jpg',im)
```
scipy.misc 模块同样包含了著名的 Lena 测试图像：
```
lena = scipy.misc.lena()
```
该脚本返回一个 512×512 的灰度图像数组。

4所有 Pylab 图均可保存为多种图像格式，方法是点击图像窗口中的“保存”按钮。

1.5　高级示例：图像去噪

我们通过一个非常实用的例子——图像的去噪——来结束本章。图像去噪是在去除图像噪声的同时，尽可能地保留图像细节和结构的处理技术。我们这里使用 ROF（Rudin-Osher-Fatemi）去噪模型。该模型最早出现在文献 [28] 中。图像去噪对于很多应用来说都非常重要；这些应用范围很广，小到让你的假期照片看起来更漂亮，大到提高卫星图像的质量。ROF 模型具有很好的性质：使处理后的图像更平滑，同时保持图像边缘和结构信息。

ROF 模型的数学基础和处理技巧非常高深，不在本书讲述范围之内。在讲述如何基于 Chambolle 提出的算法 [5] 实现 ROF 求解器之前，本书首先简要介绍一下 ROF 模型。

一幅（灰度）图像 I 的全变差（Total Variation，TV）定义为梯度范数之和。在连续表示的情况下，全变差表示为：

J(oldsymbol{I})=intleft|
ablaoldsymbol{I}
ight| ext{dx} 　　　　　　　　　　　　（1.1）

在离散表示的情况下，全变差表示为：

$J(oldsymbol{I})=sum_{ ext{x}}left| ablaoldsymbol{I} ight|$

其中，上面的式子是在所有图像坐标 x=[x, y] 上取和。

在 Chambolle 提出的 ROF 模型里，目标函数为寻找降噪后的图像 U，使下式最小：

$min_Uleft|left|oldsymbol{I}-oldsymbol{U} ight| ight|^2+2lambda J(oldsymbol{U}),$

其中范数 ||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 模型去噪后的图像