[笔记]Python计算机视觉编程《一》 基本的图像操作和处理（三）

1.3.5　图像平均

图像平均操作是减少图像噪声的一种简单方式，通常用于艺术特效。

我们可以简单地从图像列表中计算出一幅平均图像。假设所有的图像具有相同的大小，我们可以将这些图像简单地相加，然后除以图像的数目，来计算平均图像。下面的函数可以用于计算平均图像，将其添加到 imtool.py 文件里：

def compute_average(imlist): 
 """ 计算图像列表的平均图像 """
 # 打开第一幅图像，将其存储在浮点型数组中
 averageim = array(Image.open(imlist[0]), 'f')
 for imname in imlist[1:]:
    try:
        averageim += array(Image.open(imname))
    except:
        print imname + '...skipped'
        averageim /= len(imlist)
        # 返回 uint8 类型的平均图像
 return array(averageim, 'uint8')

该函数包括一些基本的异常处理技巧，可以自动跳过不能打开的图像。我们还可以

使用 mean() 函数计算平均图像。mean() 函数需要将所有的图像堆积到一个数组中；

也就是说，如果有很多图像，该处理方式需要占用很多内存。我们将会在下一节中

使用该函数

1.3.6 图像的主成分分析（PCA）

PCA（Principal Component Analysis，主成分分析） 是一个非常有用的降维技巧。

它可以在使用尽可能少维数的前提下，尽量多地保持训练数据的信息，在此意义上是一个最佳技巧。即使是一幅 100×100 像素的小灰度图像，也有 10 000 维，可以看成 10 000 维空间中的一个点。

一兆像素的图像具有百万维。由于图像具有很高的维数，在许多计算机视觉应用中，我们经常使用降维操作。

PCA 产生的投影矩阵可以被视为将原始坐标变换到现有的坐标系，坐标系中的各个坐标按照重要性递减排列。

为了对图像数据进行 PCA 变换，图像需要转换成一维向量表示。我们可以使用NumPy 类库中的flatten() 方法进行变换。

将变平的图像堆积起来，我们可以得到一个矩阵，矩阵的一行表示一幅图像。在计算主方向之前，所有的行图像按照平均图像进行了中心化。我们通常使用 SVD（Singular Value Decomposition，奇异值分解）方法来计算主成分；但当矩阵的维数很大时，SVD 的计算非常慢，所以此时通常不使用 SVD 分解。

下面就是 PCA 操作的代码：

from PIL import Image
from numpy import *
def pca(X):
   """ 主成分分析：
   输入：矩阵X，其中该矩阵中存储训练数据，每一行为一条训练数据
   返回：投影矩阵（按照维度的重要性排序）、方差和均值 """
   # 获取维数
   num_data,dim = X.shape
   # 数据中心化
   mean_X = X.mean(axis=0)
   X = X - mean_X
  if dim>num_data:
   # PCA- 使用紧致技巧
   M = dot(X,X.T) # 协方差矩阵
   e,EV = linalg.eigh(M) # 特征值和特征向量
   tmp = dot(X.T,EV).T # 这就是紧致技巧
   V = tmp[::-1] # 由于最后的特征向量是我们所需要的，所以需要将其逆转
   S = sqrt(e)[::-1] # 由于特征值是按照递增顺序排列的，所以需要将其逆转
   for i in range(V.shape[1]):
   V[:,i] /= S
  else:
   # PCA- 使用 SVD 方法
   U,S,V = linalg.svd(X)
   V = V[:num_data] # 仅仅返回前 nun_data 维的数据才合理
  # 返回投影矩阵、方差和均值
  return V,S,mean_X

该函数首先通过减去每一维的均值将数据中心化，然后计算协方差矩阵对应最大特征值的特征向量，此时可以使用简明的技巧或者 SVD 分解。这里我们使用了range() 函数，该函数的输入参数为一个整数 n，函数返回整数 0…(n-1) 的一个列表。你也可以使用 arange() 函数来返回一个数组，或者使用 xrange() 函数返回一个产生器（可能会提升速度）。我们在本书中贯穿使用 range() 函数。

如果数据个数小于向量的维数，我们不用 SVD 分解，而是计算维数更小的协方差矩阵 XXT 的特征向量。通过仅计算对应前 k（k 是降维后的维数）最大特征值的特征向量，可以使上面的 PCA 操作更快。由于篇幅所限，有兴趣的读者可以自行探索。矩阵 V 的每行向量都是正交的，并且包含了训练数据方差依次减少的坐标方向。

我们接下来对字体图像进行 PCA 变换。fontimages.zip 文件包含采用不同字体的字符 a 的缩略图。所有的 2359 种字体可以免费下载 1。

注 1：免费字体图像库由 Martin Solli 收集并上传（http://webstaff.itn.liu.se/~marso/）。

假定这些图像的名称保存在列表 imlist 中，跟之前的代码一起保存传在 pca.py 文件中，我们可以使用下面的脚本计算图像的主成分：

from PIL import Image
from numpy import *
from pylab import *
import pca
im = array(Image.open(imlist[0])) # 打开一幅图像，获取其大小
m,n = im.shape[0:2] # 获取图像的大小
imnbr = len(imlist) # 获取图像的数目
# 创建矩阵，保存所有压平后的图像数据
immatrix = array([array(Image.open(im)).flatten()for im in imlist],'f')
# 执行 PCA 操作
V,S,immean = pca.pca(immatrix)
# 显示一些图像（均值图像和前 7 个模式）
figure()
gray()
subplot(2,4,1) 
imshow(immean.reshape(m,n))
for i in range(7):
  subplot(2,4,i+2)
  imshow(V[i].reshape(m,n))
show()

这里的代码可能不可运行，因为我没找到合适的尺寸代码，resize可能会报错并且源图片也没有显示效果可能跟书上不同

注意，图像需要从一维表示重新转换成二维图像；

可以使用 reshape() 函数。

如图1-8 所示，运行该例子会在一个绘图窗口中显示 8 个图像。这里我们使用了 PyLab 库的 subplot() 函数在一个窗口中放置多个图像。

1.3.7　使用pickle模块

如果想要保存一些结果或者数据以方便后续使用，Python 中的 pickle 模块非常有用。

pickle 模块可以接受几乎所有的 Python 对象，并且将其转换成字符串表示，该过程叫做封装（pickling）。

从字符串表示中重构该对象，称为拆封（unpickling）。这些字符串表示可以方便地存储和传输。

我们来看一个例子。假设想要保存上一节字体图像的平均图像和主成分，可以这样

来完成：

#保存均值和主成分数据
f = open('font_pca_modes.pkl', 'wb')
pickle.dump(immean,f)
pickle.dump(V,f)
f.close()

在上述例子中，许多对象可以保存到同一个文件中。pickle 模块中有很多不同的协议可以生成 .pkl 文件；如果不确定的话，最好以二进制文件的形式读取和写入。

在其他 Python 会话中载入数据，只需要如下使用 load() 方法：

# 载入均值和主成分数据
f = open('font_pca_modes.pkl', 'rb')
immean = pickle.load(f)
V = pickle.load(f)
f.close()

注意，载入对象的顺序必须和先前保存的一样。Python 中有个用Ｃ语言写的优化版本，叫做 cpickle 模块，该模块和标准 pickle 模块完全兼容。关于 pickle 模块的更多内容，参见 pickle 模块文档页 http://docs.python.org/library/pickle.html。

在本书接下来的章节中，我们将使用 with 语句处理文件的读写操作。这是 Python 2.5 引入的思想，可以自动打开和关闭文件（即使在文件打开时发生错误）。

下面的例子使用 with() 来实现保存和载入操作：

# 打开文件并保存
with open('font_pca_modes.pkl', 'wb') as f:
 pickle.dump(immean,f) 
 pickle.dump(V,f)

和

# 打开文件并载入
with open('font_pca_modes.pkl', 'rb') as f:
 immean = pickle.load(f)
 V = pickle.load(f)

上面的例子乍看起来可能很奇怪，但 with() 确实是个很有用的思想。如果你不喜欢它，可以使用之前的 open 和 close 函数。

作为 pickle 的一种替代方式，NumPy 具有读写文本文件的简单函数。如果数据中不包含复杂的数据结构，比如在一幅图像上点击的点列表，NumPy 的读写函数会很有用。保存一个数组 x 到文件中，可以使用：

savetxt('test.txt',x,'%i')

最后一个参数表示应该使用整数格式。类似地，读取可以使用：

x = loadtxt('test.txt'

你可以从在线文档 http://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html 了解更多内容。

最后，NumPy 有专门用于保存和载入数组的函数。你可以在上面的在线文档里查看关于 save() 和 load() 的更多内容。

1.4 SciPy

SciPy（http://scipy.org/）是建立在 NumPy 基础上，用于数值运算的开源工具包。

SciPy 提供很多高效的操作，可以实现数值积分、优化、统计、信号处理，以及对我们来说最重要的图像处理功能。

接下来，本节会介绍 SciPy 中大量有用的模块。

SciPy 是个开源工具包，可以从 http://scipy.org/Download 下载。

1.4.1　图像模糊

图像的高斯模糊是非常经典的图像卷积例子。

图像卷积：我理解就是矩阵变换，使用卷积模板进行矩阵变换，使一个图像矩阵转为另一个图像矩阵。

本质上，图像模糊就是将（灰度）图像 I 和一个高斯核进行卷积操作：

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

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

很多其他应用。

SciPy 有用来做滤波操作的 scipy.ndimage.filters 模块。

一维分离

该模块使用快速一维分离的方式来计算卷积。你可以像下面这样来使用它：

from PIL import Image
from numpy import *
from scipy.ndimage import filters
from pylab import *
im = array(Image.open('data/empire.jpg').convert('L'))
im2 = filters.gaussian_filter(im, 5)
imshow(im2)
show()

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

图 1-9 显示了随着 σ 的增加，一幅图像被模糊的程度。

σ 越大，处理后的图像细节丢失越多。

每一个颜色通道分离

如果打算模糊一幅彩色图像，只需简单地对每一个颜色通道 即可RGB三色通道进行高斯模糊：

from PIL import Image
from numpy import *
from scipy.ndimage import filters
from pylab import *
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)
imshow(im2)
show()

在上面的脚本中，最后并不总是需要将图像转换成 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 进行描述。

描述了图像中在每个点（像素）上强度变化最大的方向。NumPy 中的 arctan2() 函数

返回弧度表示的有符号角度，角度的变化区间为 -π…π。

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

地实现：

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

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

from PIL import Image
from numpy import *
from scipy.ndimage import filters
im = array(Image.open('data/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)
imshow(imx)
subplot(1,2,1)
imshow(imy)
subplot(1,2,2)
show()

上面的脚本使用 Sobel 滤波器来计算 x 和 y 的方向导数，以及梯度大小。

sobel() 函数的第二个参数表示选择 x 或者 y 方向导数，第三个参数保存输出的变量。

图 1-10显示了用 Sobel 滤波器计算出的导数图像。

在两个导数图像中，正导数显示为亮的像素，负导数显示为暗的像素。灰色区域表示导数的值接近于零。