python 聚类效果图 python图像聚类

文章目录

• 第 6 章 图像聚类

• 引言
• 6.1 K-means聚类
• SciPy聚类包

• 图像聚类
• 在主成分上可视化图像
• 像素聚类

• 6.2　层次聚类

• 图像聚类

• 6.3　谱聚类
• 6.4 小结

第 6 章 图像聚类

引言

本章将介绍几种聚类方法，并展示如何利用它们对图像进行聚类，从而寻找相似的图像组。聚类可以用于识别、划分图像数据集，组织与导航。此外，我们还会对聚类后的图像进行相似性可视化。

6.1 K-means聚类

K-means 是一种将输入数据划分成 k 个簇的简单的聚类算法。 K-means 反复提炼初
始评估的类中心，步骤如下：

1. 以随机或猜测的方式初始化类中心 ui， i=1… k；
2. 将每个数据点归并到离它距离最近的类中心所属的类 ci；
3. 对所有属于该类的数据点求平均，将平均值作为新的类中心；
4. 重复步骤（ 2）和步骤（ 3）直到收敛。

K-means 试图使类内总方差最小：
$V=\sum_{i=1}^{k} \sum_{x_{j} \in c_{i}}\left(\boldsymbol{x}_{j}-\boldsymbol{\mu}_{i}\right)^{2}$

$x_j$

K-means 算法最大的缺陷是必须预先设定聚类数 k，如果选择不恰当则会导致聚类出来的结果很差。其优点是容易实现，可以并行计算，并且对于很多别的问题不需要任何调整就能够直接使用。

SciPy聚类包

二维数据 K-means 实现：

from pylab import *
from scipy.cluster.vq import *
import time
import matplotlib.font_manager as fm
#定义字体模板
myfont=fm.FontProperties(fname='C:/Windows/Fonts/simsun.ttc') 

def kMeans2Data():
    # 生成简单的二维数据
    class1 = 1.5 * randn(100,2)
    class2 = randn(100,2) + array([5,5])
    features = vstack((class1, class2))

    # 用k=2对这些数据进行聚类
    centroids, variance = kmeans(features, 2)

    # 由于Scipy中实现的K-means会默认计算20次，并为我们选择方法最下的结果，所有这里返回的方差并不是我们所需要的。
    # 现在，可以用Scipy包中的矢量量化函数对每个数据点进行归类
    code, distance = vq(features, centroids)

    # 通过上面得到的code，我们可以检查是否有归类错误。
    # 为了将其可视化，我们可以画出这些数据点及最终的聚类中心

    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()
    
    return 

def main():
    imagepath = []

    kMeans2Data()
    return

if __name__ =="__main__":

    startTime = time.time()
    main()
    endTime = time.time()
    print("the program running time is :",endTime - startTime)

输出如下：

上图两个绿色点为聚类的中心点，可以看到聚类基本准确。

图像聚类

让我们在 1.3.6 节的字体图像上，我们用 K-means 对这些字体图像进行聚类。文件selectedfontimages.zip 包含 66 幅来自该字体数据集 fontinages 的图像我们利用之前计算过的前 40 个主成分进行投影，用投影系数作为每幅图像的向量描述符。用 pickle 模块载入模型文件，在主成分上对图像进行投影，然后用下面的方法聚类：

def pictureCluster(imagepath):

    # 获取selected-fontimages 文件下图像文件名，并保存在列表中
    imlist = imtools.get_imlist(imagepath[0])
    imnbr = len(imlist)

    # 获取图像列表和他们的尺寸
    im = array(Image.open(imlist[0]))
    m, n = im.shape[:2]
    imnbr = len(imlist)

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

    # PCA降维
    V, S, immean = pca.pca(immatrix)

    # 载入模型文件
    with open(imagepath[1],'rb') as f:
        immean = pickle.load(f,encoding='bytes')
        #要加‘bytes’，否则会报错
        V = pickle.load(f,encoding='bytes')

    # 投影到前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 = where(code==k)[0]
        figure()
        gray()
        for i in range(minimum(len(ind), 40)):
            subplot(4,10,i+1)
            imshow(immatrix[ind[i]].reshape((25,25)))
            axis('off')
    show()

    return

输出结果如下：

在主成分上可视化图像

为了便于观察上面是如何利用主成分进行聚类的，我们可以在一对主成分方向的坐标上可视化这些图像。一种方法是将图像投影到两个主成分上，改变投影为:$\text { projected }=\operatorname{array}([\operatorname{dot}(V[[0,2]] \text {, immatrix[i]-immean }) \text { for } i \text { in range(imnbr) }])$
以得到相应的坐标（在这里 V[[0,2]] 分别是第一个和第三个主成分）。当然，你也可以将其投影到所有成分上，之后挑选出你需要的列。

def PCAcluster(imagepath,savepath):
    
    imlist = imtools.get_imlist(imagepath[0])
    imnbr = len(imlist)

    # Load images, run PCA.
    immatrix = array([array(Image.open(im)).flatten() for im in imlist], 'f')
    V, S, immean = pca.pca(immatrix)

    # Project on 2 PCs.
    projected = array([dot(V[[0, 1]], immatrix[i] - immean) for i in range(imnbr)])

    # 高和宽
    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 = 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):
        nodeim = Image.open(imlist[i])
        nodeim.thumbnail((25, 25))
        ns = nodeim.size
        box = (scaled[i][0] - ns[0] // 2, scaled[i][1] - ns[1] // 2,
            scaled[i][0] + ns[0] // 2 + 1, scaled[i][1] + ns[1] // 2 + 1)
        img.paste(nodeim, box)

    tree = hcluster.hcluster(projected)
    hcluster.draw_dendrogram(tree, imlist, filename='fonts.png')

    figure()
    imshow(img)
    axis('off')
    img.save(savepath[0]+'pca_font.png')
    show()

    return

输出图像如下：

可以看到形状相似的大致被分到了一块。

tree = hcluster.hcluster(projected)
hcluster.draw_dendrogram(tree, imlist, filename=‘fonts.png’)

上图显示了对字体图像进行层次聚类后的树状图。

像素聚类

对单幅图像中的像素而非全部图像进行聚类的例子。

将图像区域或像素合并成有意义的部分称为图像分割，它是第 9 章的主题。除了在一些简单的图像上，单纯在像素水平上应用 K-means 得出的结果往往是毫无意义的。要产生有意义的结果，往往需要更复杂的类模型而非平均像素色彩或空间一致性。现在，我们仅会在 RGB 三通道的像素值上运用 K-means 进行聚类，分割问题的处理方法会在之后谈到（ 9.2 节）给予关注及给出细节部分。

下面的代码示例载入一幅图像，用一个步长为 steps 的方形网格在图像中滑动，每滑一次对网格中图像区域像素求平均值，将其作为新生成的低分辨率图像对应位置处的像素值，并用 K-means 进行聚类：

def pixelCluster(imagepath):

    steps = 50   # 图像被划分成steps×steps的区域
    im = array(Image.open(imagepath[2]))
    dx = im.shape[0] // steps
    dy = im.shape[1] // steps

    # 计算每个区域的颜色特征
    features = []
    for x in range(steps):
        for y in range(steps):
            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, steps)
    #codeim = imresize(codeim, im.shape[:2], 'nearest')

    figure()
    ax1 = subplot(121)
    title('Original Image'), axis('off')
    imshow(im)

    ax2 = subplot(122)
    title('After Clustering'), axis('off')
    imshow(codeim)

    show()

    return

输出图像如下：

对比原始图像与聚类后的图像，可以看到主要的颜色大致就是三类。蓝色的天、影阴、和红色的墙。

6.2　层次聚类

层次聚类（或凝聚式聚类）是另一种简单但有效的聚类算法，其思想是基于样本间成对距离建立一个简相似性树。该算法首先将特征向量距离最近的两个样本归并为一组，并在树中创建一个“平均”节点，将这两个距离最近的样本作为该“平均”节点下的子节点；然后在剩下的包含任意平均节点的样本中寻找下一个最近的对，重复进行前面的操作。在每一个节点处保存了两个子节点之间的距离。遍历整个树，通过设定的阈值，遍历过程可以在比阈值大的节点位置终止，从而提取出聚类簇。

层次聚类有若干优点。例如，利用树结构可以可视化数据间的关系，并显示这些簇是如何关联的。在树中，一个好的特征向量可以给出一个很好的分离结果。另外一个优点是，对于给定的不同的阈值，可以直接利用原来的树，而不需要重新计算。

层次聚类算法代码实现：

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:
            closet = 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 < closet:
                    closet = 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=closet)
                node.remove(ni)
                node.remove(nj)
                node.append(new_node)


        return node[0]

我们为树节点创建了两个类，即ClusterNode和ClusterLeafNode，这两个类将用于创建聚类树，其中函数hcluster()用于创建树。首先创建一个包含叶节点的列表，然后根据选择的距离度量方式将距离最近的对归并到一起，返回的终节点即为树的根。对于一个行为特征向量的矩阵，运行hcluster()会创建和返回聚类树。

距离度量的选择依赖于实际的特征向量，利用欧式距离L2（同时提供了L1距离度量函数），可以创建任意距离度量函数，并将它作为参数传递给hcluster()。对于每个子树，计算其所有节点特征向量的平均值，作为新的特征向量来表示该子树，并将每个子树视为一个对象。当然，还有其他将哪两个节点合并在一起的方案，比如在两个子树中使用对象间距离最小的单向锁，及在两个子树中用对象间距离最大的完全锁。选择不同的锁会生成不同类型的聚类树。

为了从树中提取聚类簇，需要从顶部遍历树直至一个距离小于设定阈值的节点终止，这通过递归很容易做到。 ClusterNode 的 extract_clusters() 方法用于处理该过程，如果节点间距离小于阈值，则用一个列表返回节点，否则调用子节点（叶节点通常返回它们自身）。调用该函数会返回一个包含聚类簇的子树列表。对于每一个子聚类簇，为了得到包含对象 id 的叶节点，需要遍历每个子树，并用方法 get_cluster_elements() 返回一个包含叶节点的列表。

二维数据点全连接的凝聚层次聚类：

def cluster2D():
    class1 = 1.5 * randn(100, 2)
    class2 = 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())

    return

输出如下：

在这个对二维数据聚类的简单例子中，一个类中的值小于 100，另外一个大于等于 100。

图像聚类

def drawPictureCluster(imagepath):
    """ 绘制出图聚类 """

    imlist = [os.path.join(imagepath[3], f) for f in os.listdir(imagepath[3]) 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, normed=True, range=[(0, 255), (0, 255), (0, 255)])
        features[i] = h.flatten()
    tree = hcluster.hcluster(features)

    # 设置一些（任意的）阈值以可视化聚类簇
    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()

    hcluster.draw_dendrogram(tree, imlist, filename='sunset.pdf')

    return

输出图像如下：

上图为用100 幅日落图像进行层次聚类的示例聚类簇：

上图为用 100 幅日落图像进行层次聚类，将 RGB 空间的 512 个小区间直方图作为每幅图像的特征向量。树中挨的相近的图像具有相似的颜色分布。

6.3　谱聚类

谱聚类方法是一种有趣的聚类算法，与前面 K-means 和层次聚类方法截然不同。对于 n 个元素（如 n 幅图像）， 相似矩阵（或亲和矩阵，有时也称距离矩阵）是一个n× n 的矩阵，矩阵每个元素表示两两之间的相似性分数。谱聚类是由相似性矩阵构建谱矩阵而得名的。对该谱矩阵进行特征分解得到的特征向量可以用于降维，然后聚类。

谱聚类的优点之一是仅需输入相似性矩阵，并且可以采用你所想到的任意度量方式构建该相似性矩阵。像 K-means 和层次聚类需要计算那些特征向量求平均；为了计算平均值，会将特征或描述子限制为向量。而对于谱方法，特征向量就没有类别限制，只要有一个“距离”或“相似性”的概念即可。

下面说明谱聚类的过程。给定一个 n× n 的相似矩阵 S，$s_{ij}$ 为相似性分数，我们可以创建一个矩阵，称为拉普拉斯矩阵 ：$\boldsymbol{L}=\boldsymbol{I}-\boldsymbol{D}^{-1 / 2} \boldsymbol{S} \boldsymbol{D}^{-1 / 2}$
其中，$I$是单位矩阵，D是对角矩阵，对角线上的元素是S对应行元素之和，$D=\operatorname{diag}\left(d_{i}\right)$， $d_{i}=\sum_{j} s_{i j}$。拉普拉斯矩阵中的 $D^{-1/2}$ 为：
$\boldsymbol{D}^{-1 / 2}=\left[\begin{array}{llll} \frac{1}{\sqrt{d_{1}}} & & & \\ & \frac{1}{\sqrt{d_{2}}} & & \\ & & \ddots & \\ & & & \frac{1}{\sqrt{d_{n}}} \end{array}\right]$

为了使表述更简洁，对于相似性矩阵中的元素 sij，我们使用较小的值并且要求$s_{ij}$

计算 L 的特征向量，并使用 k 个最大特征值对应的 k 个特征向量，构建出一个特征向量集（记住，我们可能并没有以任何东西来开始），从而可以找到聚类簇。创建一个矩阵，该矩阵的各列是由之前求出的 k 个特征向量构成，每一行可以看做一个新的特征向量，长度为 k。这些新的特征向量可以用诸如 K-means 方法进行聚类，生成最终的聚类簇。本质上，谱聚类算法是将原始空间中的数据转换成更容易聚类的新特征向量。在某些情况下，不会首先使用聚类算法。

我们来看看真实的例子中谱聚类算法的代码。

def spectralCluster(imagepath):

    # 获取selected-fontimages 文件下图像文件名，并保存在列表中
    imlist = imtools.get_imlist(imagepath[4])
    imnbr = len(imlist)

    # Load images, run PCA.
    immatrix = array([array(Image.open(im)).flatten() for im in imlist], 'f')
    V, S, immean = pca.pca(immatrix)

    # Project on 2 PCs.
    projected = array([dot(V[[0, 1]], immatrix[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个特征向量（eigenvector）中创建特征向量（feature vector） # 叠加特征向量作为数组的列
    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()
        gray()
        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()

    return

输出结果如下：

第一类	第二类
第三类	第四类
第五类

在本例中，我们用两两间的欧式距离创建矩阵 S，并对 k 个特征向量（ eignvector）用常规的 K-means 进行聚类（在该例中， k=5）。注意，矩阵 V 包含的是对特征值进行排序后的特征向量。最后，绘制出这些聚类簇。

最后总结一下谱聚类的优缺点：

谱聚类的优点：
    1、仅仅需要输入相似性矩阵，并且可以采用所想到的任何度量方式构建该相似性矩阵。正如，K-means和层次聚类需要计算特征向量求平均；而对于谱聚类，特征向量没类别限制，只要有一个“距离”或者“相似性”的概念即可。
    2、当聚类的类别个数较小的时候，谱聚类的效果会很好，但是当聚类的类别个数较大的时候，则不建议使用谱聚类；
    3、谱聚类算法使用了降维的技术，所以更加适用于高维数据的聚类；
    4、谱聚类算法建立在谱图理论基础上，与传统的聚类算法相比，它具有能在任意形状的样本空间上聚类且收敛于全局最优解。

谱聚类的缺点：
    1、谱聚类对相似度图的改变和聚类参数的选择非常的敏感；
    2、谱聚类适用于均衡分类问题，即各簇之间点的个数相差不大，对于簇之间点个数相差悬殊的聚类问题，谱聚类则不适用。

6.4 小结

本章是引入了机器学习中的无监督学习的一个方法——K聚类。并将K聚类应用到了图像聚类上。通过SVD提取出图像的主特征，利用该特征对图像进行聚类。这是本章的核心，之后提出了像素聚类、图像聚类、还有谱聚类、层次聚类。