最近看完了 Willi Richert的《机器学习系统设计》。书虽然有点薄但也比较全，内容感觉有点偏文本处理，里面介绍了一些文本处理的方法和工具。综合起来看作为机器学习入门还是挺不错的，这里就简单记一下我做的笔记，方便回顾。书中的代码可以通过它说到的网站下载，这里是第8部分笔记。

第十章 图像的分类

机器学习有很多关于图像分析和计算机视觉方面的研究。这一章的重点其实是用mahotas包对图像进行预处理和特征提取，对于机器学习方法介绍的并不是太多，mahotas程序包是这本书的一个作者开发的。除了这个工具包外，还有scikit-image(Skimage)，Scipy中的ndimage(n维图像)模块，以及OpenCV。

图像处理

mahotas的的使用如下，首先是读取和显示图片：

1
2
3
4
5
6
7
8
9
>>> import mahotas as mh
>>> image = mh.imread('imagefile.png')

>>> image = image - image.mean() # 有助于在不同光照下进行图像的归一化。

>>> 图像的显示
>>> from matplotlib import pyplot as plt
>>> plt.imshow(image)
>>> plt.show()

阈值

卡阈值是一种非常简单的操作：我们对像素值变化，大于阈值的为1，小于阈值的为0。mahotas实现了一些选择阈值的方法。其中一个叫做Otsu，该方法的第一个必要步骤就是用rgb2gray把图像转换为灰度图。

1
2
3
>>> image = mh.colors.rgb2gray(image, dtype=np.uint8)
>>> plt.imshow(image) # 展示结果
>>> plt.gray() # 显示成灰度图

计算阈值：

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

高斯模糊

对图像进行模糊化经常被用于降噪，可以对后续的处理有所帮助，在mahotas中，只需要一个函数调用即可：

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

注意，这里并没有把灰度图转化成无符号整数，只是利用了浮点数结果。gaussian_filter函数的第二个参数是这个滤波器的大小（滤波器的标准差）。较大的值会导致结果更为模糊。通过模糊化再进行阈值分析，结果可能会含有更少的噪声。

加入椒盐噪声

加入噪声可能测试分类器对有噪声图像的处理能力，椒是黑色噪声值，盐是白色噪声值。

1
2
3
4
5
6
7
salt = np.random.random(image.shape) > .975
pepper = np.random.random(image.shape) > .975

image = mh.stretch(image)
image = np.maximum(salt*170, sep)
image = np.minimum(peper*30 + image * (~pepper), image)
## 数值定为170和30比纯白和纯黑要平滑一些。

聚焦中心

最后一个例子显示了如何用Numpy操作和一些滤波混合起来：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
# 首先把它切分成几个颜色通道
im = mh.imread('lenna.jpg')
r, g, b = im.transpose(2,0,1)

# 分别对三个通道进行滤波，并用mh.as_rgb()构建出一个合成图像。
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) # 这个函数接受3个二维数组，
# 进行对比度拉伸，并把每个二维数组转换成一个8位整形数组，然后把他们叠放起来

# 从中心到边界将两个图片混合在一起，首先我们需要构建一个权重数组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里增加一个虚设的第三维度

# 最后，我们把两幅图像组合起来，让中间成为焦点，让边上比较柔和。
ringed = mh.stretch(im * C + (1-C) * im12)

模式识别

对于图像的分类，如果直接把每个像素值都当作特征放到机器学习算法中，可能效果并不好，因为每个像素和最终结果之间的关系是非常间接的。传统的方法是从图像中计算特征，然后把这些特征用于分类。有一些方法可以直接对像素进行操作，他们有特征计算子模块。他们甚至试图自动从图像里面学出好的特征，这些都是当今学术界正在研究的课题。由于历史的原因，图像分类又叫做模式识别然而它就是将分类方法应用于图像。

特征设计和提取

使用mahotas可以很容易的提取图像的特征。有一个子模块叫做mahotas.features ,包含了一些特征计算函数。

一个很常用的特征集合叫做Haralick纹理特征。这些特征是基于纹理的：它们会对平滑的图像和带有模式的图像进行区分。用mahotas很容易计算这些特征：

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

mh.features.haralick函数返回了一个4*13的数组，第一维表示特征的四个方向（上，下，左，右）。如果对方向不感兴趣，也可以对整体计算平均值。

mahotas还有一些其他的特征集合。线性二元模式（linear binary pattern）是另一个基于纹理的特征集合，它对光亮变化非常健壮。

设计一个新的特征并不是很难的事情，为了不同的图像分类目的，不同的特征表现可能差异很大。我们介绍一个边界寻找的例子，并把它加到特征中去，我们使用sobel滤波。从数学上说，我们用两个矩阵对图像滤波（求卷积）；竖向矩阵如下所示

$$
\begin{pmatrix}
1 & 0 & 1 \\
-2 & 0 & -2 \\
1 & 0 & 1 \\
\end{pmatrix}
$$

横向矩阵

$$
\begin{pmatrix}
1 & -2 & 1 \\
0 & 0 & 0 \\
1 & -2 & 1 \\
\end{pmatrix}
$$

然后我们把结果的平方相加，得到对每一点锐度的一个综合估计。直接用mahotas进行sobel滤波：

1
2
3
filtered = mh.sobel(image, just_filter = True)
# just_filter=True这个参数是必要的，否则他就会用阈值进行过滤，
# 并得到对边界位置的一个估计。

局部特征表示

在计算机视觉领域一个较新的进展就是基于局部特征的方法。*局部特征(local feature)是在图像的一小块区域内计算出来的特征。mahotas支持这类特征中的一种：加速稳健特征(Speeded Up Robust Feature)又叫做SURF；还有一些其他的特征，如尺度不变特征变换(Scale-Invariant Feature Trandform, SIFT)*。这些局部特征对于旋转和光照变化十分稳健。

在使用这些特征的时候，我们需要确定在哪里计算它们。这里有3个经常采用的选项：

  1. 随机计算；
  2. 在一个格子里计算；
  3. 检测图像中的兴趣区域（这种技术叫做关键点检测(keypoint detection)）

mahotas对这三种方式都提供了支持。如果你确定的兴趣点可以和图像里的重要区域相对应，那么兴趣点检测的效果将会是最好的。通常他们对人造图像的效果比自然风景图像要好。人造景观有较强的角度边界，以及高对比度的区域。这些通常会被自动检测器识别为兴趣区域。

1
2
3
4
5
6
7
8
9
10
11
from mathotas.features import surf

# 使用兴趣点方法
descriptors = surf.surf(image, descriptors_only = True)
# descriptors_only=True标志是说，我们只对他们的描述符感兴趣，对他们的像素位置、大小和其他信息并不感兴趣。

# 使用密度采样的方法
descriptors = surf.dense(image, spacing = 16)
# 它返回了描述符的值，这个值是从相互之间相距16像素的点集中计算出来的。
# 由于这些点的位置是固定的，所以我们对兴趣点的元信息并不是很感兴趣，它也不会默认返回。
# 不管在哪种情况下，它的结果（这些描述符）都是一个n*64的数组，（n是抽样点的个数）

我们不能直接把这些描述符传进机器学习分类器中，要使用图像中的描述符，可以采用词袋模型：把图像中看起来相似的区域聚成一组，把它们叫做视觉词语。词语个数一般对算法的最终效果并没有很大影响。不过，如果这个数字特别小，那么系统的效果不会很好，特别多的时候，整个系统表现也不会很好。在两个极端之间，可以试的区间会很大，经验法则，如果你有很多很多图片，采用诸如256，512，1024这样的数值应该可以得到不错的效果。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
alldescriptors = []
for im in images:
    im = mh.imread(im, as_grey=True)
    im = im.astype(np.uint8)
    alldescriptors.append(surf.surf(im,descriptors_only))
# 这样就可以得到超过100 000个局部描述符

# 使用k-means得到聚类中心
concatenated = np.concatenated(alldescriptors) # 把所有描述符放进一个数组中
concatenated = concatenated[::32] # 只使用第32个元素所组成的向量
from sklearn.cluster import Kmeans
k = 256
km = KMeans(k)
km.fit(concatenated)
# 此时km中包含所有关于聚类中心的信息。

# 构建特征向量
features = []
for d in alldescriptors:
    c = km.predict(d)
    features.append(
        np.array([np.sum(c == ci) for ci in range(k)])
    )
features = np.array(features)
# features[fi] 是一个跟图像位置fi相对应的直方图

上面得到的features就是最后算出来的特征空间，使用features进行分类可以获得一定的性能提升。