计算机视觉编程 第六章

前言

本专栏按《python计算机视觉编程 ——Jan Erik Solem》一书为参考，第六章介绍用K-means方法进行图像聚类，进而寻找到相似的图像组合，此外还介绍了层次聚类和谱聚类两种聚类方法和示例

6.1 K-means聚类

在之前的博客中已经学习过K-means聚类，这里再简单提一下：它是一种无监督机器学习算法，用于将数据点划分成不同的群组（簇），使得同一簇内的数据点彼此相似，而不同簇之间的数据点差异较大。目标是将数据点划分到K个簇中，其中K是预先指定的聚类数量

K-means算法的步骤主要如下：

1. 初始化: 选择K个初始的聚类中心$u_i$，$i=\mathrm{1....}k$。可以是随机选择的数据点，或者其他方法来选择
2. 分配数据点: 对于每个数据点，计算其与各个聚类中心的距离，然后将该数据点分配到距离最近的聚类中心所属的簇中
3. 更新聚类中心: 对每个簇，计算其内部所有数据点的平均值（中心），然后将该平均值作为新的聚类中心。
4. 重复步骤2和步骤3: 迭代执行步骤2和步骤3，直到聚类中心不再发生明显变化，或者达到预定的迭代次数

K-means目的是使聚类的总方差最小$$V=\sum _{i=1}^k\sum _{x_j\in c_j}(x_j-u_i{)}^2$$其中$x_j$是输入的矢量数据。K-means的优点包括简单易实现、计算效率高等。然而，它也有一些局限性，比如对初始聚类中心的选择敏感，容易收敛到局部最优解，而不是全局最优解，在算法开始需要人为确定k的值。此外，K-means假定每个簇的形状是球状的，并且每个维度的权重相等，这在某些情况下可能不适用

6.1.1 SciPy聚类包K-means

使用SciPy聚类包可以很容易实现K-means算法，代码如下

from scipy.cluster.vq import *
from pylab import *
from numpy import *

class1 = 1.5 * random.randn(100, 2)
class2 = random.randn(100, 2) + array([5, 5])
features = vstack((class1, class2))

centroids, variance = kmeans(features, 2)

code, distance = vq(features, centroids)

figure()
ndx = where(code == 0)[0]
plot(features[ndx, 0], features[ndx, 1], '*')
ndx = where(code == 1)[0]
plot(features[ndx, 0], features[ndx, 1], 'r.')
plot(centroids[:, 0], centroids[:, 1], 'go')
axis('off')

show()

效果如上，类中心标记为绿色大圆环，预测出的类用蓝色星号和红色点标记出

6.1.2 图像聚类

这里使用的图像数据集是MNIST数据集，下载和解析请查看此网站，使用第一章的PCA主成分分析进行投影计算，并下述代码进行聚类

import pickle
from PIL import Image
from pylab import *
from numpy import *
from scipy.cluster.vq import *
from 视觉编程 import imtools


def pca(x):
    """

    :param x: 矩阵x，存储训练数据，每一行为一个训练数据
    :return: 投影矩阵、方差和均值
    """
    # 获取维数
    num_data, dim = x.shape
    # 数据中心化
    mean_x = x.mean(axis=0)
    x = x - mean_x
    if dim > num_data:
        # 协方差矩阵
        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] /= 5
    else:
        U, S, V = linalg.svd(x)
        V = V[:num_data]
    return V, S, mean_x


# 获取seleced-fontimages文件下图像文件名，并保存在列表中
imlist = imtools.get_imlist('./mnist_train/6/')
im = array(Image.open(imlist[0]))
m, n = im.shape[0:2]
imnbr = len(imlist)

# # 载入模型文件
# with open('a_pca_modes.pkl', 'rb') as f:
#     immean = pickle.load(f)
#     V = pickle.load(f)

# 创建矩阵，存储所有拉成一组形式后的图像
immartix = array([array(Image.open(im)).flatten() for im in imlist], 'f')

V, S, immean = pca(immartix)

# 投影到前40个主成分上
immean = immean.flatten()
projected = array([dot(V[:40], immartix[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 = where(code == k)[0]
    figure()
    gray()
    for i in range(minimum(len(ind), 40)):
        subplot(4, 10, i + 1)
        imshow(immartix[ind[i]].reshape((m, n)), cmap='gray')  # 将形状修改为原始图像的大小
        axis('off')
show()

效果如上，只能说一言难尽，反正我是没有看出来聚类的明显区别，但是把k值变为3，主成分数目改成30时，产生了比较明显的聚类效果，如下

6.1.3 在主成分上可视化图像

由于效果不是很理想，尝试使用可视化的方法进行主成分聚类，一种方法是将图像投影到两个主成分上，改变投影为

projected = array([dot(V[:30], immartix[i] - immean) for i in range(imnbr)])

在之前代码的基础上去除K-means的部分加上以下部分，这里选取了前四十个图像作为样本

# 绘制坐标轴
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 = abs(projected).max(0)
scaled = floor(array([(p/scale)*(w/2-20,h/2-20)+(w/2,h/2) for p in projected])).astype(int)

for i in range(imnbr)[:40]:
    nodeim = Image.open(imlist[i])
    # 图像灰度转换
    nodeim_array = array(nodeim)
    nodeim_array = 255 - nodeim_array
    nodeim = Image.fromarray(nodeim_array)

    nodeim.thumbnail((m, n))
    ns = nodeim.size

    # 计算矩形框的左上角和右下角坐标
    box_left = scaled[i][0] - ns[0] // 2
    box_upper = scaled[i][1] - ns[1] // 2
    box_right = box_left + ns[0]
    box_lower = box_upper + ns[1]
    box = (box_left, box_upper, box_right, box_lower)

    img.paste(nodeim, box)

img.save('pca_font.jpg')

还是能看到相似的字体距离都比较近的，综合来看使用第一个和第二个主成分的效果好

主成分1，2	主成分2，3	主成分1，3｜

6.1.4像素聚类

接下来尝试对一整副图像进行聚类，就能得到非常nice的像素风格图像，但是只有黑白的。下面就是用一个步长为steps的方形网格在图像中滑动，每滑一次对网格中图像区域像素求平均值，将其作为新生的低分辨率图像对应位置处的像素值，并用K-means进行聚类:

from scipy.cluster.vq import *
from PIL import Image
from pylab import *

steps = [50,100,150]  # 图像被划分为steps*steps的区域
im = array(Image.open('./filelist/tp3.jpg'))
subplot(2,2,1)
imshow(im), axis('off')
for i in range(3):
    dx = im.shape[0] // steps[i]
    dy = im.shape[1] // steps[i]
    # 计算每个区域的颜色特征
    features = []
    for x in range(steps[i]):
        for y in range(steps[i]):
            R = mean(im[x * dx:(x + 1) * dx, y * dy:(y + 1) * dy, 0])
            G = mean(im[x * dx:(x + 1) * dx, y * dy:(y + 1) * dy, 1])
            B = mean(im[x * dx:(x + 1) * dx, y * dy:(y + 1) * dy, 2])
            features.append([R, G, B])
    features = array(features, 'f')  # 变成数组
    # 聚类
    centroids, variance = kmeans(features, 3)
    code, distance = vq(features, centroids)
    # 用聚类标记创建图像
    codeim = code.reshape(steps[i], steps[i])

    subplot(2,2,i+2)
    imshow(codeim), axis('off')

show()

上面这几张图分别对应于5050、100100和150*150窗口大小对图像的聚类效果，窗口越大和原图像素越接近，尽管用窗口进行下采样，但是可以看到还是有噪声的

6.2 层次聚类

层次聚类是一种基于树状结构的聚类方法，它将数据集中的对象逐步合并为越来越大的簇，直到所有对象最终合并成一个大簇，形成一个聚类树或者称为“谱系图”，这种做法有点类似数据结构中的霍夫曼编码。层次聚类可以分为两种主要类型：凝聚型和分裂型。这两种类型在聚类过程中的策略不同，就像他们的名字一样前者是从树枝到树根，后者则反之，这里主要介绍前。层次聚类的步骤通常包括以下几个步骤：

1. 距离/相似度计算： 计算数据点之间的距离或相似度，这用于确定合并或分裂哪些簇
2. 簇合并/分裂： 根据距离/相似度计算，决定哪些簇应该合并或分裂。对于凝聚型，选择距离最小的两个簇合并，对于分裂型，选择合适的簇进行分裂
3. 更新距离矩阵： 在每次合并/分裂后，需要更新距离矩阵以反映新的簇间距离或相似度。
   重复步骤2和3： 反复执行合并/分裂和更新距离矩阵，直到所有数据点都合并成一个大簇（凝聚型）或每个数据点都成为一个独立的簇（分裂型）
4. 构建聚类树： 在此过程中，每次的合并或分裂就合并起来就形成一个层次聚类树

利用本书附的hcluster.py中的ClusterNode和ClusterLeafNode两个类创建聚类树:每个 ClusterNode 对象代表一个聚类子树，其由两个子节点（左子节点和右子节点）组成，以及它们之间的距离信息。它的主要作用是支持层次聚类过程中的节点合并操作，以及提供从树中提取子树簇的方法。它还可以绘制聚类树的结构，将子节点连接起来并显示相关信息.每个 ClusterLeafNode 对象代表一个单独的数据点，具有其自己的特征向量和标识。它的主要作用是支持层次聚类过程中的数据点初始化，并提供从树中提取单一元素簇的方法。它还可以绘制图像缩略图，将数据点的视觉表示绘制在聚类树的结构上

下面用一个示例来观察聚类的过程：

from 视觉编程 import hcluster
class1 = 1.5 * random.randn(100,2)
class2 = random.randn(100,2) + array([5,5])
features = vstack((class1,class2))

tree = hcluster.hcluster(features)
clusters = tree.extract_clusters(5)

print('number of clusters', len(clusters))
for c in clusters:
    print(c.get_cluster_elements())

上面是两次的实验结果，一次3个类一次4个类，且四个类的效果稍好点。在四个聚类中第一个类都满足大于100，后三个则反之，与书上的理想案例算是比较接近的

图像聚类

这里是用基于图像颜色信息的聚类案例，使用此网盘中的sunset文件进行操作

# 创建图像列表
path = './sunsets/flickr-sunsets-small'
imlist = [os.path.join(path, f) for f in os.listdir(path) if f.endswith('.jpg')]
# 提取特征向量，每个颜色通道量化成8个小区间
features = zeros([len(imlist), 512])
for i, f in enumerate(imlist):
    im = array(Image.open(f))
    # 多维直方图
    h, edges = histogramdd(im.reshape(-1, 3), 8, density=True,
                           range=[(0, 255), (0, 255), (0, 255)])
    features[i] = h.flatten()
tree = hcluster.hcluster(features)
# 保存树状图
hcluster.draw_dendrogram(tree, imlist, filename='sunset.pdf')
# 设置一些阈值，得到可视化聚类簇
clusters = tree.extract_clusters(0.23 * tree.distance)
# 绘制聚类簇中元素超过3个的图像
for c in clusters:
    elements = c.get_cluster_elements()
    nbr_elements = len(elements)
    if nbr_elements > 3:
        figure()
        for p in range(minimum(nbr_elements, 20)):
            subplot(4, 5, p + 1)
            im = array(Image.open(imlist[elements[p]]))
            imshow(im)
            axis('off')
show()

聚类簇	树状图

和书上的结果一样，聚类效果标准

6.3 谱聚类

谱聚类的原理是利用数据的图谱信息，有效地处理数据中的非线性和复杂结构。它的核心思想是将数据转换为图形表示，然后在图上进行聚类操作。因此不论是图像还是数据都能够在谱聚类下进行聚类，谱聚类的步骤如下：

1. 数据表示和相似度矩阵S构建：首先将原始数据转换为图形表示。通常是使用欧几里德距离、余弦相似度等计算数据点之间的相似性分数$s_{ij}$。这些相似度值构成了相似度矩阵，矩阵中的每个元素表示对应数据点之间的相似度
2. 构建拉普拉斯矩阵L：从相似度矩阵中构建拉普拉斯矩阵，以捕捉了数据点之间的连接关系，其表示为$L=I-D^{-1/2}S{D}^{-1/2}$，$I$为单位矩阵，$D$为对角矩阵
3. 特征值分解：对构建的拉普拉斯矩阵进行特征值分解。特征值分解将拉普拉斯矩阵分解为特征值和对应的特征向量。这些特征向量反映了数据在图结构上的分布。通常情况下，选择前k个最小的非零特征值对应的特征向量，其中k是聚类的数量
4. 降维和聚类：特征值分解得到的特征向量构成了一个新的低维表示空间。可以将这些特征向量作为数据的新特征，在这个新的低维空间中进行聚类操作
5. 聚类结果映射：聚类完成后，可以将结果映射回原始数据空间，得到每个数据点所属的聚类标签

谱聚类的优点有：能够处理复杂的数据结构，且不受数据维度的限制，适用于高维数据。然而，谱聚类的计算复杂度较高，特别是在特征值分解步骤。参数选择对结果影响较大，在选择相似度度量、拉普拉斯矩阵类型、聚类数量等参数问题上都需要一定的技术

下面是一个实际谱聚类算法代码

imlist = imtools.get_imlist('./mnist_train/7/')
imnbr = len(imlist)
im = array(Image.open(imlist[0]))
m, n = im.shape[0:2]

# 创建矩阵，存储所有拉成一组形式后的图像
immartix = array([array(Image.open(im)).flatten() for im in imlist], 'f')

# 灰度逆转
immartix = 255 - immartix

V, S, immean = pca(immartix)
#
# 投影到前40个主成分上
immean = immean.flatten()
projected = array([dot(V[[1,2]], immartix[i] - immean) for i in range(imnbr)])

n = len(projected)
# 计算距离矩阵
S = array([[sqrt(sum((projected[i] - projected[j]) ** 2)) for i in range(n)] for j in range(n)], 'f')

# 创建拉普拉斯矩阵
rowsum = sum(S, axis=0)
D = diag(1 / sqrt(rowsum))
I = identity(n)
L = I - dot(D, dot(S, D))

# 计算矩阵L的特征向量
U, sigma, V = linalg.svd(L)
k = 5
# 从矩阵L的前k个特征向量中创建特征向量
# 叠加特征向量作为数组的列
features = 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 = where(code == c)[0]
    figure()
    for i in range(minimum(len(ind), 39)):
        im = Image.open(imlist[ind[i]])
        subplot(4, 10, i + 1)
        imshow(array(im))
        axis('equal')
        axis('off')
show()

由于图像大小有限，这里只放了四张图，但是还是能看到聚类的效果是不错的，但是有极个别的数字独树一帜比如第二类的第二行右数第四个和最下角的两个明显与其他数字有差别，但是整体来看效果还是可以的