python计算机视觉学习———图像聚类

6.1 K-means聚类

K-means 是一种将输入数据划分成 k 个簇的简单的聚类算法。K-means 反复提炼初始评估的类中心。
k-means算法中的k代表类簇个数，means代表类簇内数据对象的均值（这种均值是一种对类簇中心的描述），因此，k-means算法又称为k-均值算法。k-means算法是一种基于划分的聚类算法，以距离作为数据对象间相似性度量的标准，即数据对象间的距离越小，则它们的相似性越高，则它们越有可能在同一个类簇。
距离的度量：聚类简单理解就是把相似的东西聚到一起，如何判断两个样本点是不是相似的呢，这个就要根据距离做判断，最常见的计算方式就是欧几里得距离（直接算两个点的欧式距离）和余弦相似度（先标准化），当使用欧式距离的时候，先要对数据进行标准化，
什么是标准化：比如现在有两个维度，x轴，y轴，x轴数据0.01,0.02,0.04，y轴数据100,200,300，当计算相似度的时候，x轴的差异无论怎么算都比较小，y轴的差异无论怎么算都比较大，那样我们潜意识里就认为相似度主要由y轴决定，实际算出来也是这样，所以说，在使用距离的度量的时候，基本情况下，都要对所有数据进行标准化，比如让x轴取值范围在0到1之间，y轴取值范围在0到1之间，让数据x和y基本在一个比较小的范围内浮动，比如说都是0到1，或者-1到1。

先把数据做标准化，然后再用距离的度量看一下什么样的两个样本点是相似的，再把相似的分到一簇。

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

K-means 试图使类内总方差最小：

xj 是输入数据，并且是矢量。该算法是启发式提炼算法，在很多情形下都适用，但
是并不能保证得到最优的结果。为了避免初始化类中心时没选取好类中心初值所造
成的影响，该算法通常会初始化不同的类中心进行多次运算，然后选择方差 V 最小
的结果。

6.1.1 SciPy 聚类包

SciPy 矢量量化包scipy.cluster.vq 中有 K-means 的实现。
首先，生成简单二维数据，再令K=2，对这些数据进行聚类。最后可视化

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

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

一个对二维数据用 K-means 进行聚类的示例。类中心标记为绿色大圆环，预测出的
类分别标记为蓝色星号和红色点。

6.1.2 图像聚类

用 K-means 对字体图像进行聚类。

# -*-coding:utf-8 -*-
from PCV.tools import imtools
import pickle
from scipy.cluster.vq import *
from numpy import *
from PIL import Image
from pylab import *
# 获取 selected-fontimages 文件下图像文件名，并保存在列表中
imlist = imtools.get_imlist('D:\\project\\a_selected_thumbs\\')
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 = 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()

6.1.3 在主成分上可视化图像

为了便于观察上面是如何利用主成分进行聚类的，我们可以在一对主成分方向的坐标上可视化这些图像。

# -*- coding: utf-8 -*-
from PCV.tools import imtools, pca
from PIL import Image, ImageDraw
from pylab import *
from PCV.clustering import hcluster

imlist = imtools.get_imlist('D:\\project\\a_selected_thumbs\\')
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)])
# projected = array([dot(V[[1, 2]], 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('D:\\project\\pca_font.png')
show()

这里用到了整数或floor向下取整除法运算//，通过移去小数点后面的部分，可以返回各个缩略图在白色背景中对应的整数坐标位置。这类图像说明这些字体图像在40维里的分布情况，对于选择一个好的描述子很有帮助。看到，二维投影后相似的字体图像距离较近。
这是生成聚类树的效果图：

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

6.2 层次聚类

层次聚类（或凝聚式聚类）是另一种简单但有效的聚类算法，其思想是基于样本间成对距离建立一个简相似性树。该算法首先将特征向量距离最近的两个样本归并为一组，并在树中创建一个“平均”节点，将这两个距离最近的样本作为该“平均”节点下的子节点；然后在剩下的包含任意平均节点的样本中寻找下一个最近的对，重复进行前面的操作。在每一个节点处保存了两个子节点之间的距离。遍历整个树，通过设定的阈值，遍历过程可以在比阈值大的节点位置终止，从而提取出聚类簇。
层次聚类可以分为凝聚(agglomerative)层次聚类和分裂(divsive)层次聚类。分裂层次聚类采用的就是"自顶而下"的思想，先将所有的样本都看作是同一个簇，然后通过迭代将簇划分为更小的簇，直到每个簇中只有一个样本为止。凝聚层次聚类采用的是"自底向上"的思想，先将每一个样本都看成是一个不同的簇，通过重复将最近的一对簇进行合并，直到最后所有的样本都属于同一个簇为止。
层次聚类有若干优点。例如，利用树结构可以可视化数据间的关系，并显示这些簇是如何关联的。在树中，一个好的特征向量可以给出一个很好的分离结果。另外一个优点是，对于给定的不同的阈值，可以直接利用原来的树，而不需要重新计算。不足之处是，对于实际需要的聚类簇，我们需要选择一个合适的阈值。

在凝聚层次聚类中，判定簇间距离的两个标准方法就是单连接(single linkage)和全连接(complete linkage)。单连接，是计算每一对簇中最相似两个样本的距离，并合并距离最近的两个样本所属簇。全连接，通过比较找到分布于两个簇中最不相似的样本(距离最远)，从而来完成簇的合并。
基于全连接的凝聚层次聚类主要包括下面几个步骤：

1、获取所有样本的距离矩阵

2、将每个数据点作为一个单独的簇

3、基于最不相似(距离最远)样本的距离，合并两个最接近的簇

4、更新样本的距离矩阵

5、重复2到4，直到所有样本都属于同一个簇为止。

聚类模式：
1)自底向上型（agglomerative）
2)自上向下型（divisive ）

层次聚类的合并算法通过计算两类数据点间的相似性，对所有数据点中最为相似的两个数据点进行组合，并反复迭代这一过程。简单的说层次聚类的合并算法是通过计算每一个类别的数据点与所有数据点之间的距离来确定它们之间的相似性，距离越小，相似度越高。并将距离最近的两个数据点或类别进行组合，生成聚类树。
欧几里德距离矩阵

（https://blog.csdn.net/sinat_29957455/article/details/80146093）博文有详细代码距离过程。
其中最后一步，这里详细阐述：

s1与s5合并之后得到新簇S6，S2与S3合并之后得到新簇S7，继续比较距离，继续循环，合并得到上图结果。

层次聚类算法在代码中体现：

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

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

# -*- coding: utf-8 -*-
import os
from PCV.clustering import hcluster
from matplotlib.pyplot import *
from numpy import *
from PIL import Image

# 创建图像列表
path = 'D:\\Python\\chapter6\\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, 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')

将 R、G、B 三个颜色通道作为特征向量，将其传递到 NumPy 的 histogramdd()中，该函数能够计算多维直方图（本例中是三维）。我们在每个颜色通道中使用 8 个小区间进行量化，将三个通道量化后的小区间拉成一行后便可用 512（8×8×8）维的特征向量描述每幅图像。为避免图像尺寸不一致，我们用“normed=True”归一化直方图，并将每个颜色通道范围设置为 0…255。将 reshape() 第一个参数设置为 -1会自动确定正确的尺寸，故可以创建一个输入数组来计算以 RGB 颜色值为行向量的直方图。

树状图的高和子部分由距离决定，随着坐标向下传递到下一级，会递归绘制出这些节点，上述代码用 20×20 像素绘制叶节点的缩略图，使用 get_height() 和 get_depth() 这两个辅助函数可以获得树的高和宽。

6.3 谱聚类

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

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

谱聚类原理阐述：
谱聚类是从图论中演化出来的算法，后来在聚类中得到了广泛的应用。它的主要思想是把所有的数据看做空间中的点，这些点之间可以用边连接起来。距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高，通过对所有数据点组成的图进行切图，让切图后不同的子图间边权重和尽可能的低，而子图内的边权重和尽可能的高，从而达到聚类的目的。
谱聚类基础之一：无向权重图
由于谱聚类是基于图论的，因此我们首先温习下图的概念。对于一个图G，我们一般用点的集合V和边的集合E来描述。即为G(V,E)。其中V即为我们数据集里面所有的点(v1,v2,…vn)。对于V中的任意两个点，可以有边连接，也可以没有边连接。我们定义权重wij为点vi和点vj之间的权重。由于我们是无向图，所以wij=wji。
对于有边连接的两个点vi和vj，wij>0,对于没有边连接的两个点vi和vj，wij=0。对于图中的任意一个点vi，它的度di定义为和它相连的所有边的权重之和，即

利用每个点度的定义，我们可以得到一个nxn的度矩阵D,它是一个对角矩阵，只有主对角线有值，对应第i行的第i个点的度数，定义如下：

利用所有点之间的权重值，我们可以得到图的邻接矩阵W，它也是一个nxn的矩阵，第i行的第j个值对应我们的权重wij。

除此之外，对于点集V的的一个子集A⊂V，我们定义：
|A|:=子集A中点的个数

谱聚类基础之二：相似矩阵
在上一节我们讲到了邻接矩阵W，它是由任意两点之间的权重值wij组成的矩阵。通常我们可以自己输入权重，但是在谱聚类中，我们只有数据点的定义，并没有直接给出这个邻接矩阵，那么怎么得到这个邻接矩阵呢？

基本思想是，距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高，不过这仅仅是定性，我们需要定量的权重值。一般来说，我们可以通过样本点距离度量的相似矩阵S来获得邻接矩阵W。

构建邻接矩阵W的方法有三类。ϵ-邻近法，K邻近法和全连接法。

对于ϵ-邻近法，它设置了一个距离阈值ϵ，然后用欧式距离sij度量任意两点xi和xj的距离。即相似矩阵的sij=||xi−xj||2, 然后根据sij和ϵ的大小关系，来定义邻接矩阵W如下：

从上式可见，两点间的权重要不就是ϵ,要不就是0，没有其他的信息了。距离远近度量很不精确，因此在实际应用中，我们很少使用ϵ-邻近法。

第二种定义邻接矩阵W的方法是K邻近法，利用KNN算法遍历所有的样本点，取每个样本最近的k个点作为近邻，只有和样本距离最近的k个点之间的wij>0。但是这种方法会造成重构之后的邻接矩阵W非对称，我们后面的算法需要对称邻接矩阵。为了解决这种问题，一般采取下面两种方法之一：

第三种定义邻接矩阵W的方法是全连接法，相比前两种方法，第三种方法所有的点之间的权重值都大于0，因此称之为全连接法。可以选择不同的核函数来定义边权重，常用的有多项式核函数，高斯核函数和Sigmoid核函数。最常用的是高斯核函数RBF，此时相似矩阵和邻接矩阵相同：

谱聚类基础之三：拉普拉斯矩阵
单独把拉普拉斯矩阵(Graph Laplacians)拿出来介绍是因为后面的算法和这个矩阵的性质息息相关。它的定义很简单，拉普拉斯矩阵L=D−W。D即为我们第二节讲的度矩阵，它是一个对角矩阵。而W即为我们第二节讲的邻接矩阵，它可以由我们第三节的方法构建出。

拉普拉斯矩阵有一些很好的性质如下：

1）拉普拉斯矩阵是对称矩阵，这可以由D和W都是对称矩阵而得。

2）由于拉普拉斯矩阵是对称矩阵，则它的所有的特征值都是实数。

3）对于任意的向量f,我们有
　　　　
　　　　这个利用拉普拉斯矩阵的定义很容易得到如下：
　　　　
　　　　4) 拉普拉斯矩阵是半正定的，且对应的n个实数特征值都大于等于0，即0=λ1≤λ2≤…≤λn， 且最小的特征值为0，这个由性质3很容易得出。
谱聚类基础之四：无向图切图
对于无向图G的切图，我们的目标是将图G(V,E)切成相互没有连接的k个子图，每个子图点的集合为：A1,A2,…Ak，它们满足Ai∩Aj=∅,且A1∪A2∪…∪Ak=V.

对于任意两个子图点的集合A,B⊂V, A∩B=∅, 我们定义A和B之间的切图权重为：

那么对于我们k个子图点的集合：A1,A2,…Ak，我们定义切图cut为：

谱聚类过程主要有两步，第一步是构图，将采样点数据构造成一张网图，表示为G(V,E)，V表示图中的点，E表示点与点之间的边，如下图：

第二步是切图，即将第一步构造出来的按照一定的切边准则，切分成不同的图，而不同的子图，即我们对应的聚类结果，举例如下：

下面以Ncut总结谱聚类算法流程。

输入：样本集D=(x1,x2,…,xn)，相似矩阵的生成方式, 降维后的维度k1, 聚类方法，聚类后的维度k2
　　　　输出： 簇划分C(c1,c2,…ck2).

1) 根据输入的相似矩阵的生成方式构建样本的相似矩阵S

2）根据相似矩阵S构建邻接矩阵W，构建度矩阵D

3）计算出拉普拉斯矩阵L

4）构建标准化后的拉普拉斯矩阵
　　　　这里写图片描述
　　　　5）计算
　　　　这里写图片描述
　　　　最小的k1个特征值所各自对应的特征向量f
　　　　6) 将各自对应的特征向量f组成的矩阵按行标准化，最终组成n×k1维的特征矩阵F

7）对F中的每一行作为一个k1维的样本，共n个样本，用输入的聚类方法进行聚类，聚类维数为k2。

8）得到簇划分C(c1,c2,…ck2).　　　
谱聚类算法实现：

# -*- coding: utf-8 -*-
from PCV.tools import imtools, pca
from PIL import Image, ImageDraw
from pylab import *
from scipy.cluster.vq import *

imlist = imtools.get_imlist('D:\\project\\a_selected_thumbs\\')
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)])  # P131 Fig6-3左图
# projected = array([dot(V[[1, 2]], immatrix[i] - immean) for i in range(imnbr)])  # P131 Fig6-3右图

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