Python计算机视觉 第6章-图像聚类
6.1 K-means聚类
K-means 聚类 是一种常用的无监督学习算法，用于将数据集划分为 ( K ) 个簇。
算法步骤
- 选择 ( K ) 个初始簇中心（可以是随机选择或其他启发式方法）。
- 将每个数据点分配到距离其最近的簇中心所属的簇。
- 计算每个簇的新簇中心，通常是簇中所有点的均值。
- 重复步骤 2 和 3，直到簇中心不再改变或达到预设的迭代次数。
目标函数
K-means 试图最小化以下目标函数，即每个点到其簇中心的距离平方和：
J = ∑ i = 1 K ∑ x ∈ C i ∥ x − μ i ∥ 2 J = \sum_{i=1}^{K} \sum_{x \in C_i} \| x - \mu_i \|^2 J=i=1∑Kx∈Ci∑∥x−μi∥2
其中，( C_i ) 是第 ( i ) 个簇，( \mu_i ) 是第 ( i ) 个簇的簇中心。
优点
- 算法简单且易于实现。
- 计算效率较高，适合大规模数据集。
缺点
- 需要预先指定簇的数量 ( K )。
- 对初始簇中心的选择敏感，可能导致局部最优。
- 不适用于形状不规则或不同密度的簇。
6.1.1 SciPy聚类包
SciPy矢量量化包scipy.cluster.vq 中有 K-means 的实现，下面是使用方法。
实验代码如下：
from matplotlib import pyplot as plt
import numpy as np
from scipy.cluster.vq import kmeans, vq
# 生成数据
class1 = 1.5 * np.random.randn(100, 2)
class2 = np.random.randn(100, 2) + np.array([5, 5])
features = np.vstack((class1, class2))
# 聚类
centroids, variance = kmeans(features, 2)
code, distance = vq(features, centroids)
# 可视化
plt.figure()
# 绘制属于第 0 类的数据点
ndx = np.where(code == 0)[0]
plt.plot(features[ndx, 0], features[ndx, 1], '*')
# 绘制属于第 1 类的数据点
ndx = np.where(code == 1)[0]
plt.plot(features[ndx, 0], features[ndx, 1], 'r.')
# 绘制簇中心
plt.plot(centroids[:, 0], centroids[:, 1], 'go')
plt.axis('off')
plt.show()
6.1.2 图像聚类
用下面的方法聚类：
import imtools
import pickle
from scipy.cluster.vq import kmeans, vq, whiten
from numpy import array, dot
from PIL import Image
# 获取 selected_fontimages 文件下图像文件名，并保存在列表中
imlist = imtools.get_imlist('selected_fontimages/')
imnbr = len(imlist)
# 载入模型文件
with open('a_pca_modes.pkl', 'rb') as f:
immean = pickle.load(f)
V = pickle.load(f)
# 创建矩阵，存储所有拉成一组形式后的图像
immatrix = array([array(Image.open(im)).flatten() for im in imlist], 'f')
# 投影到前40个主成分上
immean = immean.flatten()
projected = array([dot(V[:40], immatrix[i] - immean) for i in range(imnbr)])
# 进行 k-means 聚类
projected = whiten(projected)
centroids, distortion = kmeans(projected, 4)
code, distance = vq(projected, centroids)
利用下面的代码可以可视化聚类后的结果：
# 绘制聚类簇
for k in range(4):
ind = np.where(code == k)[0]
plt.figure()
plt.gray()
for i in range(min(len(ind), 40)):
plt.subplot(4, 10, i + 1)
plt.imshow(immatrix[ind[i]].reshape((25, 25)))
plt.axis('off')
plt.show()
6.1.3 在主成分上可视化图像
为了便于观察上面是如何利用主成分进行聚类的，我们可以在一对主成分方向的坐标上可视化这些图像。一种方法是将图像投影到两个主成分上，改变投影为：
projected = np.array([np.dot(V[[0, 2]], immatrix[i] - immean) for i in range(imnbr)])
以得到相应的坐标（在这里V[[0,2]]分别是第一个和第三个主成分）。当然，你也可以将其投影到所有成分上，之后挑选出你需要的列。
结果如下：
我们用PIL中的ImageDraw模块进行可视化。假设你有如上所示投影后的图像，及保存有图像文件名的列表，利用下面简短的脚本可以生成如图6-3所示的结果：
from PIL import Image, ImageDraw
import numpy as np
# 高和宽
h, w = 1200, 1200
# 创建一幅白色背景图
img = Image.new('RGB', (w, h), (255, 255, 255))
draw = ImageDraw.Draw(img)
# 绘制坐标轴
draw.line((0, h / 2, w, h / 2), fill=(255, 0, 0))
draw.line((w / 2, 0, w / 2, h), fill=(255, 0, 0))
# 缩放以适应坐标系
scale = np.abs(projected).max(0)
scaled = np.floor(np.array([(p / scale) * (w / 2 - 20, h / 2 - 20) + (w / 2, h / 2) for p in projected]))
# 粘贴每幅图像的缩略图到白色背景图片
for i in range(imnbr):
nodeim = Image.open(imlist[i])
nodeim.thumbnail((25, 25))
ns = nodeim.size
img.paste(nodeim, (int(scaled[i][0] - ns[0] // 2), int(scaled[i][1] - ns[1] // 2),
int(scaled[i][0] + ns[0] // 2 + 1), int(scaled[i][1] + ns[1] // 2 + 1)))
img.save('pca_font.jpg')
这里，我们用到了整数或floor向下取整除法运算符//，通过移去小数点后面的部分，可以返回各个缩略图在白色背景中对应的整数坐标位置。
这类图像说明这些字体图像在40维里的分布情况，对于选择一个好的描述子很有帮助。可以很清楚地看到，二维投影后相似的字体图像距离较近。
6.1.4 像素聚类
将图像区域或像素合并成有意义的部分称为图像分割。
除了在一些简单的图像上，单纯在像素水平上应用K-means得出的结果往往是毫无意义的。要产生有意义的结果，往往需要更复杂的类模型而非平均像素色彩或空间一致性。现在，我们仅会在RGB三通道的像素值上运用K-means进行聚类。
下面的代码示例载入一幅图像，用一个步长为steps的方形网格在图像中滑动，每滑一次对网格中图像区域像素求平均值，将其作为新生成的低分辨率图像对应位置处的像素值，并用K-means进行聚类：
from scipy.cluster.vq import kmeans, vq
from scipy.misc import imresize
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
steps = 50 # 图像被划分成 steps×steps 的区域
im = np.array(Image.open('empire.jpg'))
dx = im.shape[0] / steps
dy = im.shape[1] / steps
# 计算每个区域的颜色特征
features = []
for x in range(steps):
for y in range(steps):
R = np.mean(im[int(x*dx):(int(x+1)*dx), int(y*dy):(int(y+1)*dy), 0])
G = np.mean(im[int(x*dx):(int(x+1)*dx), int(y*dy):(int(y+1)*dy), 1])
B = np.mean(im[int(x*dx):(int(x+1)*dx), int(y*dy):(int(y+1)*dy), 2])
features.append([R, G, B])
features = np.array(features, 'f') # 变为数组
# 聚类
centroids, variance = kmeans(features, 3)
code, distance = vq(features, centroids)
# 用聚类标记创建图像
codeim = code.reshape(steps, steps)
codeim = imresize(codeim, im.shape[:2], interp='nearest')
plt.figure()
plt.imshow(codeim)
plt.show()
K-means 的输入是一个有steps×steps行的数组，数组的每一行有3列，各列分别为区域块R、G、B三个通道的像素平均值。为可视化最后的结果,我们用SciPy的imresize() 函数在原图像坐标中显示这幅steps×steps的图像。参数interp指定插值方法；我们在这里采用最近邻插值法，以便在类间进行变换时不需要引入新的像素值。
图6-4显示了用50×50和100×100窗口对两幅相对比较简单的示例图像进行像素聚类后的结果。注意，K-means标签的次序是任意的（在这里的标签指最终结果中图像的颜色）。正如你所看到的，尽管利用窗口对它进行了下采样，但结果仍然是有噪声的。如果图像某些区域没有空间一致性，则很难将它们分开，如下方图中小男孩和草坪的图。空间一致性和更好的分割方法以及其他的图像分割算法会在后面讨论。现在，让我们继续来看下一个基本的聚类算法。
6.2 层次聚类
层次聚类（或凝聚式聚类）是另一种简单但有效的聚类算法，其思想是基于样本间成对距离建立一个简相似性树。该算法首先将特征向量距离最近的两个样本归并为一组，并在树中创建一个“平均”节点，将这两个距离最近的样本作为该“平均”节点下的子节点；然后在剩下的包含任意平均节点的样本中寻找下一个最近的对，重复进行前面的操作。在每一个节点处保存了两个子节点之间的距离。遍历整个树，通过设定的阈值，遍历过程可以在比阈值大的节点位置终止，从而提取出聚类簇。
层次聚类有若干优点。例如，利用树结构可以可视化数据间的关系，并显示这些簇是如何关联的。在树中，一个好的特征向量可以给出一个很好的分离结果。另外一个优点是，对于给定的不同的阈值，可以直接利用原来的树，而不需要重新计算。
不足之处是，对于实际需要的聚类簇，我们需要选择一个合适的阈值。关键代码如下：
from itertools import combinations
from numpy import array, sqrt, sum, abs
class ClusterNode(object):
def __init__(self, vec, left, right, distance=0.0, count=1):
self.left = left
self.right = right
self.vec = vec
self.distance = distance
self.count = count # 只用于加权平均
def extract_clusters(self, dist):
""" 从层次聚类树中提取距离小于 dist 的子树簇群列表 """
if self.distance < dist:
return [self]
return self.left.extract_clusters(dist) + self.right.extract_clusters(dist)
def get_cluster_elements(self):
""" 在聚类子树中返回元素的 id """
return self.left.get_cluster_elements() + self.right.get_cluster_elements()
def get_height(self):
""" 返回节点的高度，高度是各分支的和 """
return self.left.get_height() + self.right.get_height()
def get_depth(self):
""" 返回节点的深度，深度是每个子节点取最大再加上它的自身距离 """
return max(self.left.get_depth(), self.right.get_depth()) + self.distance
class ClusterLeafNode(object):
def __init__(self, vec, id):
self.vec = vec
self.id = id
def extract_clusters(self, dist):
return [self]
def get_cluster_elements(self):
return [self.id]
def get_height(self):
return 1
def get_depth(self):
return 0
def L2dist(v1, v2):
return sqrt(sum((v1 - v2) ** 2))
def L1dist(v1, v2):
return sum(abs(v1 - v2))
def hcluster(features, distfcn=L2dist):
""" 用层次聚类对行特征进行聚类 """
# 用于保存计算出的距离
distances = {}
# 每行初始化为一个簇
node = [ClusterLeafNode(array(f), id=i) for i, f in enumerate(features)]
while len(node) > 1:
closest = float('Inf')
# 遍历每对，寻找最小距离
for ni, nj in combinations(node, 2):
if (ni, nj) not in distances:
distances[ni, nj] = distfcn(ni.vec, nj.vec)
d = distances[ni, nj]
if d < closest:
closest = d
lowestpair = (ni, nj)
ni, nj = lowestpair
# 对两个簇求平均
new_vec = (ni.vec + nj.vec) / 2.0
# 创建新的节点
new_node = ClusterNode(new_vec, left=ni, right=nj, distance=closest)
node.remove(ni)
node.remove(nj)
node.append(new_node)
return node[0]
我们为树节点创建了两个类，即ClusterNode和ClusterLeafNode，这两个类将用于创建聚类树，其中函数hcluster()用于创建树。首先创建一个包含叶节点的列表，然后根据选择的距离度量方式将距离最近的对归并到一起，返回的终节点即为树的根。对于一个行为特征向量的矩阵，运行hcluster()会创建和返回聚类树。
我们为树节点创建了两个类，即ClusterNode和ClusterLeafNode，这两个类将用于创建聚类树，其中函数hcluster()用于创建树。首先创建一个包含叶节点的列表，然后根据选择的距离度量方式将距离最近的对归并到一起，返回的终节点即为树的根。对于一个行为特征向量的矩阵，运行hcluster()会创建和返回聚类树。
为了从树中提取聚类簇，需要从顶部遍历树直至一个距离小于设定阈值的节点终止，这通过递归很容易做到。ClusterNode的extract_clusters()方法用于处理该过程，如果节点间距离小于阈值，则用一个列表返回节点，否则调用子节点（叶节点通常返回它们自身）。调用该函数会返回一个包含聚类簇的子树列表。对于每一个子聚类簇，为了得到包含对象id的叶节点，需要遍历每个子树，并用方法get_cluster_ elements() 返回一个包含叶节点的列表。
下面，我们在一个简单的例子中观察该聚类过程。首先创建一些二维数据点（和之前K-means一样）：
class1 = 1.5 * randn(100, 2)
class2 = randn(100, 2) + array([5, 5])
features = vstack((class1, class2))
对这些数据点进行聚类，设定阈值（这里的阈值设定为5），从列表中提取这些聚类簇，并在控制台打印出来：
import hcluster
tree = hcluster.hcluster(features)
clusters = tree.extract_clusters(5)
print('number of clusters', len(clusters))
for c in clusters:
print(c.get_cluster_elements())
打印结果应该与下面类似：
理想情况下，你应该得到两一个聚类簇，但是在实际数据中，你可能会得到三类或更多这主要依赖于实际生成的二维数据。在这个对二维数据聚类的简单例子中，一个类中的值应该小于100，另外一个应该大于等于100。
6.3 谱聚类
谱聚类是一种常用的聚类技术，基于图论中的谱图理论。以下是谱聚类的基本过程：
-
构建相似矩阵：
给定一个 ( n \times n ) 的相似矩阵 ( S )，其中 ( S_{ij} ) 表示样本 ( i ) 和样本 ( j ) 之间的相似性分数。 -
构建拉普拉斯矩阵：
- 度矩阵 ( D )：构建一个对角矩阵 ( D )，其中 ( D_{ii} ) 是第 ( i ) 行（或列）在相似矩阵 ( S ) 中的和，即
D i i = ∑ j S i j D_{ii} = \sum_{j} S_{ij} Dii=j∑Sij - 拉普拉斯矩阵 ( L )：根据 ( D ) 和 ( S ) 构建拉普拉斯矩阵 ( L )。常用的拉普拉斯矩阵有三种形式：
- 标准拉普拉斯矩阵：
L = D − S L = D - S L=D−S - 归一化拉普拉斯矩阵（对称归一化）：
L s y m = I − D − 1 / 2 S D − 1 / 2 L_{sym} = I - D^{-1/2} S D^{-1/2} Lsym=I−D−1/2SD−1/2 - 归一化拉普拉斯矩阵（随机游走归一化）：
L r w = I − D − 1 S L_{rw} = I - D^{-1} S Lrw=I−D−1S
- 标准拉普拉斯矩阵：
- 度矩阵 ( D )：构建一个对角矩阵 ( D )，其中 ( D_{ii} ) 是第 ( i ) 行（或列）在相似矩阵 ( S ) 中的和，即
-
计算特征值和特征向量：
对拉普拉斯矩阵 ( L ) 进行特征分解，得到特征值和特征向量。选择前 ( k ) 个最小的特征值对应的特征向量（通常是前 ( k ) 个特征向量）。 -
聚类：
使用这些特征向量构成的矩阵来进行聚类。通常，特征向量会被用来进行 ( k )-均值聚类，得到最终的簇分配。 -
结果分析：
分析聚类结果，检查每个簇的样本分布，并评估聚类的效果。
我们来看看真实的例子中谱聚类算法的代码：
from scipy.cluster.vq import *
import numpy as np
from numpy import linalg
from PIL import Image
n = len(projected)
# 计算距离矩阵
S = np.array([[np.sqrt(np.sum((projected[i] - projected[j]) ** 2))
for i in range(n)] for j in range(n)], 'f')
# 创建拉普拉斯矩阵
rowsum = np.sum(S, axis=0)
D = np.diag(1 / np.sqrt(rowsum))
I = np.identity(n)
L = I - np.dot(D, np.dot(S, D))
# 计算矩阵L的特征向量
U, sigma, V = linalg.svd(L)
k = 5
# 从矩阵L的前k个特征向量（eigenvector）中创建特征向量（feature vector）
# 叠加特征向量作为数组的列
features = np.array(V[:k]).T
# k-means 聚类
features = whiten(features)
centroids, distortion = kmeans(features, k)
code, distance = vq(features, centroids)
# 绘制聚类簇
for c in range(k):
ind = np.where(code == c)[0]
plt.figure()
for i in range(min(len(ind), 39)):
im = Image.open(path + imlist[ind[i]])
plt.subplot(4, 10, i + 1)
plt.imshow(np.array(im))
plt.axis('equal')
plt.axis('off')
plt.show()
在本例中，我们用两两间的欧式距离创建矩阵S，并对k个特征向量（eignvector）用常规的K-means进行聚类（在该例中，k=5）。 注 意 ， 矩 阵 V包含的是对特征值进行排序后的特征向量。最后，绘制出这些聚类簇。图6-8显示了运行后的聚类簇；需要记住的是，在K-means阶段，每次运行的结果可能不同。
在该例中选择k有些技巧。很多人会认为这里只有两类（即白宫的两侧），以及其他一些垃圾图像。用k=2可以得到类似图6-9的结果，其中一个聚类簇是包含很多白宫一侧的图像，另一个聚类簇是白宫另一侧的图像和其他所有垃圾图像。将k设定为一个较大的值，比如k=10，则有些聚类簇可能只包含一幅图像（很可能是垃圾图像），另一些是真实的聚类簇。图6-10给出了上面示例代码运行的结果。在该例中，仅有两个真实的聚类簇，每个聚类簇包含白宫一个侧面的图像。
这里展示的谱聚类算法有很多不同的版本供选择，它们对如何构造矩阵L和如何处理特征向量有各自不同的思想。
