Question

为了对图像进行操作，我们会使用一个叫做mahotas的程序包。它是一个开源程序包（经MIT许可，它可以用在任何项目里），是由本书的一个作者开发出来的。幸运的是，它是基于NumPy 的。所以到目前为止你所掌握的NumPy知识都可以用于图像处理。这里还有其他图像包，如Scikit-image （Skimage ），SciPy 中的ndimage （n 维图像）模块，以及OpenCV中的Python 绑定。所有这些都原生支持NumPy，所以你可以把来自不同程序包的函数混合在一起使用。

我们从引入mahotas 开始，在本章中用mh 这个缩写来表示它：

import mahotas as mh

现在用imread 读取一个图像文件：

image = mh.imread('imagefile.png')

如果imagefile.phg包含的是一个高为h 宽为w 的彩色图像，那么image 就是一个形为(h,w,3) 的数组。第一维是高度，第二维是宽度，第三维是红色/绿色/蓝色。其他系统可能会把宽度放在第一维，但这是一个数学上的惯例，所有基于NumPy的程序包都是这样使用的。数组元素的类型通常是np.unit8 （8位无符号整数）。这就是你的照相机所拍摄的或显示器全屏所能显示的图像。

然而，一些专业设备（主要是在科学领域）会拍出更高分辨率的图像。一般是12位或16位。mahotas可以处理所有这些图像，包括浮点数值 的图像。（并不是所有操作对于浮点数都有意义，但如果要这样做，mahotas也都可以支持。）在很多计算中，即使原始数据里包含的是无符号整数，把它们转化为浮点数也是有用处的，这样可以简化对舍入和溢出问题的处理。

提示 　mahotas可以使用各种不同的输入/输出后端。但不幸的是，它们不能读取所有现有的图像格式（有几百种格式，每一种又有一些变种）。但是，它们都支持对PNG和JPEG图像的读取。所以我们只聚焦在这些常见格式上。对于不常见格式的读取，你可以参考mahotas的文档。

mh.imread 的返回值是一个NumPy数组。这意味着你可以使用标准NumPy函数来处理图像。例如，从图像里减去像素均值通常是很有用处的操作。它有助于在不同光照条件下对图像进行归一化，这个操作可以用标准的mean 方法来实现：

image = image – image.mean()

我们可以用matplotlib 在屏幕上把图像展示出来。这个绘图工具库我们已经使用过几次了：

from matplotlib import pyplot as plt plt.imshow(image) plt.show()

这张图像遵循了是第一维高度、第二维宽度的惯例。它对彩色图像一样可以正确处理。在使用Python进行数值计算的时候，整个系统都能很好地协作，这使我们受益匪浅。

10.2.1　图像处理基础

我们从一个特意为本书而收集的小规模数据集开始。它有3个类别：建筑物、自然景色（风景），和文字图像。每一类里有30份图像，是用手机摄像头拍摄的。所以这些图像和用户上传到网站的那些图像是类似的。这个数据集可以从本书的网站上获得。在本章的后面，我们还将会看到一个更加困难的数据集，有更多的图像和更多的类别。

这个建筑物是该数据集中里一张图像。我们拿它作为一个例子。

也许你已经意识到了，图像处理是一个很大的领域。这里只会涉及一些非常基本的操作。有一些最基本的操作只用NumPy就可以实现，但其他的操作我们将使用mahotas。

1. 阈值

卡阈值是一种非常简单的操作：我们对像素值变换，大于一定阈值的是1 ，而其他小于阈值的是0 （或者用Booleans，把它转换为True 或False ）：

binarized = (image > threshold_value)

我们需要选择阈值的宽度值（代码中的threshold_value ）。如果图像都非常相似，那么我们可以基于统计选择一个，并把它应用于其他所有图像，否则必须基于像素值对每张图像都计算一个不同的阈值。

mahotas实现了一些选择阈值的方法。其中一个叫做Otsu，它是以它的发明者命名的。该方法的第一个必要步骤就是用rgb2gray 把图像转换为灰度图。

除了rgb2gray ，我们还可以通过调用image.mean(2) 得到红、绿、蓝通道的均值。然而，它的结果和rgb2gray 并不相同，这是因为rgb2gray 对不同的颜色使用了不同的权重，给出的是一个在主观上更让人感到愉悦的结果。我们的眼睛对3种基本色的敏感程度并不相同。

image = mh.colors.rgb2gray(image, dtype=np.uint8) plt.imshow(image) # 展示图像

matplotlib 会默认把这个单通道图像显示为假彩色图像，较高的值用红色，较低的值用蓝色。对于自然图像，灰度图则更为适合。你可以用下述方法调用：

plt.gray()

现在图像显示成灰度图了。注意只有像素值的解释和显示方式变化了，这幅图的其他地方并没有改变。我们可以继续进行处理，并计算阈值：

thresh = mh.thresholding.otsu(image) print(thresh) imshow(image > thresh)

把这个方法应用到前面那张图的时候，该方法发现了阈值164，它可以把建筑物、停放的车辆以及上面的天空分离出来。

这个结果本身可能就很有用（如果想评估这幅阈值图的一些属性），它也可以用于后续处理。最后的结果是一个二值图像，可以用于选择感兴趣的区域。

这个结果仍然不是非常好。我们可以在这幅图上进行一些操作来进一步调优。例如，我们可以运行close 操作器，去除上边角落里面的一些噪声：

otsubin = (image <= thresh)
otsubin = mh.close(otsubin, np.ones((15,15)))

在这里，我们是在把低于阈值的区域封闭起来，我们也可以把这个阈值操作反过来做，在反向图中进行open 操作：

otsubin = (image > thresh) otsubin = mh.open(otsubin, np.ones((15,15)))

不管哪种情况，这个操作都会对一个结构元素进行处理，该元素定义了我们要关闭的区域的类型。在这里，我们使用的结构元素是一个15×15的方块。

这仍然不够完美，因为在停车区域里有一些光点没有去掉。我们会在本章的后面进一步提升效果。

Otsu阈值可以把天空区域识别得比建筑物更加明亮。另一种阈值方法叫做Ridley-Calvard 方法（也是以它的发明者命名的）：

thresh = mh.thresholding.rc(image) print(thresh)

该方法返回了一个更小的阈值137.7，并且把建筑物的细节区分了出来。

这样是好是坏，取决于你想区分出什么来。

2. 高斯模糊

对图像进行模糊化似乎有点奇怪，但它经常被用于降噪，可以对后续的处理有所帮助。在mahotas里，只需要一个函数调用即可：

image = mh.colors.rgb2gray(image) im8 = mh.gaussian_filter(image,8)

注意，我们并没有把灰度图转化为无符号整数，只是利用了浮点数结果。gaussian_filter 函数的第二个参数是这个滤波器的大小（滤波器的标准差）。较大的值会导致结果更为模糊，如我们在下图中所看到的那样（显示出了大小为8、16和32的滤波效果）：

我们可以对左边那个图用Otsu阈值卡一下（代码和之前一样）。现在这个结果就是建筑区域和天空区域的一个完美划分。在一些细节被平滑掉的同时，停车区域内的光点也被平滑掉了。这个结果是天空的近似轮廓，不存在任何伪造。通过模糊化，我们把和总体布局无关的细节都去掉了。看看下面这张图：

3. 不同效果的滤波

通过图像处理技术获得令人愉悦的图像效果，可以追溯到数字图像刚刚出现的时候。但在最近，它已经变成了很多有趣应用的基础，其中最有名的可能要属Instagram。

我们将要使用图像处理中的一幅传统图像，Lenna图像。它可以从本书网站（或其他很多图像处理网站）下载：

im = mh.imread('lenna.jpg', as_grey=True)

10.2.2　加入椒盐噪声

我们可以在这个结果上进行很多后续处理，如果我们愿意的话。例如，我们可以在图像里加入一点椒盐噪声，来模拟扫描噪点。我们生成一个跟原始图像相同宽度高度的随机数组。其中只有1%的值是True 。

salt = np.random.random(lenna.shape) > .975 pepper = np.random.random(lenna.shape) > .975

我们现在添加椒（意思是几乎为黑色的值）盐（意思是几乎为白色的值）噪声：

lenna = mh.stretch(lenna) lenna = np.maximum(salt*170, sep) lenna = np.minimum(pepper*30 + lenna*(~pepper), lenna)

我们把数值170 当做白色，30 当做黑色。这比极端值255 和0 要平滑一些。然而，这些都是选项，是由人的主观选择和做事方式决定的。

聚焦中心

最后这个例子显示了如何把NumPy操作和一些滤波混合起来，得到一个有趣的结果。我们从Lenna图像开始，并把它切分成几个颜色通道：

im = mh.imread('lenna.jpg') r,g,b = im.transpose(2,0,1)

现在我们分别对这3个通道进行滤波，并用mh.as_rgb 构建出一个合成图像。这个函数接收3个二维数组，进行对比度拉伸，并把每个二维数组转换为一个8位整型数组，然后把它们叠放起来：

r12 = mh.gaussian_filter(r, 12.) g12 = mh.gaussian_filter(g, 12.) b12 = mh.gaussian_filter(b, 12.) im12 = mh.as_rgb(r12,g12,b12)

然后我们从中心到边界将两个图片混合在一起。首先我们需要构建一个权重数组W ，它包含每个像素的归一化结果，这是每个像素到中心的距离：

h,w = r.shape # 高度和宽度 Y,X = np.mgrid[:h,:w]

我们使用了np.mgrid 对象，它返回一个大小为（h,w） 的数组，里面的值分别对应于y 和x 坐标：

Y = Y-h/2. # 中心在h/2 Y = Y / Y.max() # 归一化到-1 .. +1 X = X-w/2. X = X / X.max()

现在用一个高斯函数赋予中心区域一个高权重：

W = np.exp(-2.*(X**2+ Y**2)) # 再次归一化到0...1 W = W - C.min() W = W / C.ptp() W = C[:,:,None] # 这会在W里增加一个虚设的第三维度

注意这些操作是如何通过NumPy数组实现的，它们并不是mahotas特有的方法。这就是Python NumPy生态系统的一个优点：你在学习纯机器学习方法时所学到的所有操作，在一个完全不同的背景下，依然是有用的。

最后，我们把两幅图像组合起来，让中间成为焦点，让边上较为柔和。

ringed = mh.stretch(im*C + (1-C)*im12)

既然你已经了解图像滤波的一些基本技术，就可以基于它生成新的滤波器了。这一点与其说是一种技术，还不如说是一种艺术。

10.2.3　模式识别

在对图像进行分类的时候，我们从一个较大规模的数组（包含像素值）开始。如今，几百万像素的图像已经很常见了。

我们可以把所有这些数字都当做特征放进学习算法。但这并不是一个很好的主意。因为每个像素（甚至每个像素组）和最终结果之间的关系是非常间接的。相反，传统方法是从图像中计算特征，然后把这些特征用于分类。

有一些方法可以直接对像素进行操作。它们有特征计算子模块。它们甚至试图自动从图像里学出好特征。这些都是当今学术界正在研究的课题。

我们前面使用了一个建筑类别的例子。这里还有一些文本和风景类别的例子。

提示　模式识别就是图像分类

由于历史的原因，图像分类又叫做模式识别。然而它就是将分类方法应用于图像。自然，图像也有它自己的特定问题，我们将会在本章里处理这些问题。

10.2.4　计算图像特征

采用mahotas，计算图像特征就变得非常容易。有一个子模块叫做mahotas.features ，包含了一些特征计算函数。

一个很常用的特征集合叫做Haralick纹理特征 。和特征处理里面的众多方法一样，这个方法也是以它的发明者命名的。这些特征是基于纹理的：它们会对平滑的图像和带有模式的图像进行区分。用mahotas很容易计算这些特征：

haralick_features = np.mean(mh.features.haralick(image),0)

mh.features.haralick 函数返回了一个4x13数组。第一维代表计算特征的4个可能方向（上、下、左、右）。如果我们对方向不感兴趣，那可以采用整体方向的平均值。基于这个函数，我们很容易构建一个分类系统。

mahotas里还实现了一些其他的特征集合。线性二元模式（linear binary pattern）是另一个基于纹理的特征集合，它对光亮变化非常健壮。另外还有一些其他类型的特征，包括局部特征，在本章后面将会进行讨论。

提示　特征并不只是为分类而存在

基于特征对百万像素图像降维的方法，还可以用在其他机器学习情景中，包括聚类、回归，或维度归约。通过计算几百个特征，然后在结果中运行降维算法，你可以把一张百万像素的图像变成少数几个维度的特征，你构建的可视化工具甚至可以是二维的。

有了这些特征，我们就可以使用标准的分类方法进行分类，例如支持向量机：

images = glob('simple-dataset/*.jpg') features = [] labels = [] for im in images: features.append(mh.features.haralick(im).mean(0)) labels.append(im[:-len('00.jpg')]) features = np.array(features) labels = np.array(labels)

这三种图像类别具有非常不一样的纹理。建筑物图像有清晰的边沿，以及颜色相似的大区块（像素值很少恰好一样，但差异很微小）。文本图像是由很多鲜明的明暗过度所组成的，在白色中有一些较小的黑色区域。自然风景图像有一些类分形的平滑变化。因此，基于纹理的分类器应该可以做得不错。不过由于数据集很小，我们用逻辑回归只得到了79%的正确率。

10.2.5　设计你自己的特征

图像特征并不是什么神奇的事情。它就是从图像里计算出来的数值。本文里已经给出了一些特征集合。它们通常有一个额外优点，就是对很多不重要的因素都具有不变性。例如线性二元模式就对“像素值乘以一个数或与常量相加”这个操作完全不变。这使得它对于图像光照变化十分健壮。

然而，你的特定问题也可能会受益于一些特别设计的特征。例如，要从自然图像中区分出文字，定义轮廓分明的文字特征就很重要。我们并不是要知道文字本身是什么（它们可能是轮廓分明的或方块的），而是说文本图像会有很多边界。因此，我们希望引入一种“锐度特征”。有一些方法可以实现这个功能（无限多）。机器学习系统的一个优点就是，我们只需要写出一些想法，然后就可以让系统找出哪些是好的，哪些不太好。

我们开始介绍另一个传统图像处理操作：边界寻找。在这里，我们将使用sobel滤波 。从数学上说，我们用两个矩阵对图像滤波（求卷积）；竖向矩阵如下图所示：

横向矩阵：

然后我们把结果的平方相加，得到对每一点锐度的一个综合估计（在其他情况下，你可能需要把竖边和横边区分开，并以另外一种方式使用；当然，这取决于具体应用）。mahotas支持sobel滤波，如下所示：

filtered = mh.sobel(image, just_filter=True)

just_filter=True 这个参数是必要的，否则它就会用阈值进行过滤，并得到对边界位置的一个估计。下图显示了应用这个滤波器的结果（左图），以及卡阈值方法得到的结果（右图）：

基于这个操作元，我们可以定义一个全局特征，作为它的综合锐度：

def edginess_sobel(image): edges = mh.sobel(image, just_filter=True) edges = edges.ravel() return np.sqrt(np.dot(edges, edges))

在最后一行里，我们使用了一个技巧来计算均方根——用内积函数np.dot 和用np.sum(edges ** 2) 是等价的，但运算速度要快很多（我们仅需要确保先把数组分解开了）。自然，我们还可以想出很多不同方式来得到相似的结果。比如一个明显的例子是，使用卡阈值的方式计算出大于阈值的像素比例。

我们很容易把这个特征添加到之前的管道里：

features = [] for im in images: image = mh.imread(im) features.append(np.concatenate( mh.features.haralick(im).mean(0), # 用我们的特征构建一个单元素列表，与np.concatenate相匹配 [edginess_sobel(im)], ))

用这个结构很容易把特征集合融合进来。在使用了所有这些特征之后，我们得到了84%的正确率。

这个例子完美解释了这样一个原则：好算法只是比较容易的那个部分。你总可以找到一个前沿的分类方法来实现。但真正的秘密和附加价值通常是在特征设计和特征工程里面。这就是数据本身知识的价值所在。