目录
聚类可以用于识别、划分图像数据集，组织与导航。此外，我们还会对聚类后的图像进行相似性可视化。
6.1 K-means聚类
K-means 是一种将输入数据划分成k 个簇的简单的聚类算法。K-means 反复提炼初始评估的类中心，步骤如下：
- 以随机或猜测的方式初始化类中心\(u_i\)，i=1…k；
- 将每个数据点归并到离它距离最近的类中心所属的类\(c_i\)；
- 对所有属于该类的数据点求平均，将平均值作为新的类中心；
- 重复步骤（2）和步骤（3）直到收敛。
K-means 试图使类内总方差最小：
\(V=\sum_{i=1}^{k}\sum_{x_{j}\in c_{i}}(x_{j}-\mu_{i})^{2}\)
其中，\(x_{j}\)是输入数据，并且是矢量。该算法是启发式提炼算法，在很多情形下都适用，但是不能保证得到最优的结果。为了避免初始化类中心时没选取好类中心初值所造成的影响，该算法通常会初始化不同的类中心进行多次运算，然后选择方差V最小的结果。
K-means算法的缺陷是：必须预先设定聚类数k，如果选择不恰当则会导致聚类出来的结果很差。当样本集规模大时，收敛速度会变慢；对孤立点数据敏感，少量噪声就会对平均值造成较大影响；k的取值十分关键，对不同数据集，k选择没有参考性，需要大量实验。
K-means算法优点是：容易实现，可以并行计算，并且对于很多别的问题不需要任何调整就能够直接使用；聚类效果较优；算法的可解释度比较强；主要需要调参的参数仅仅是簇数k。
K-means算法的关键：在于初始中心的选择和距离公式。
6.1.1 SciPy聚类包
K-means算法很容易实现，可以使用Scipy矢量量化包scipy.clusterr.vq中有K-means的实现，下面是使用方法。
#导入scipy中K-means的相关工具
from scipy.cluster.vq import *
#randn是NumPy中的一个函数
from numpy import *
from pylab import *
#生成简单的二维数据:生成两类二维正态分布数据
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 ,我们可以检查是否有归类错误
code,distance = vq(features,centroids)
#可视化结果：画出这些数据点及最终的聚类中心：函数 where() 给出每个类的索引
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 图像聚类
from PCV.tools import imtools
import pickle
from scipy import *
from pylab import *
from PIL import Image
from scipy.cluster.vq import *
from PCV.tools import pca
from sympy import minimum
# Uses sparse pca codepath.
imlist = imtools.get_imlist('selectedfontimages/a_selected_thumbs')
# 获取图像列表和他们的尺寸
im = array(Image.open(imlist[0])) # open one image to get the size
m, n = im.shape[:2] # get the size of the images
imnbr = len(imlist) # get the number of images
print ("The number of images is %d" % imnbr)
# Create matrix to store all flattened images
immatrix = array([array(Image.open(imname)).flatten() for imname in imlist], 'f')
# PCA降维
V, S, immean = pca.pca(immatrix)
# 保存均值和主成分
#f = open('./a_pca_modes.pkl', 'wb')
f = open('./a_pca_modes.pkl', 'wb')
pickle.dump(immean,f)
pickle.dump(V,f)
f.close()
# get list of images
imlist = imtools.get_imlist('selectedfontimages/a_selected_thumbs/')
imnbr = len(imlist)
# load model file
with open('a_pca_modes.pkl','rb') as f:
immean = pickle.load(f)
V = pickle.load(f)
# create matrix to store all flattened images
immatrix = array([array(Image.open(im)).flatten() for im in imlist],'f')
# project on the 40 first PCs
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)
# plot clusters
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 在主成分上可视化图像
为了便于观察上面如何利用主成分进行聚类的，可以在一对主成分方向的坐标上可视化这些图像。一种方法是将图像投影到两个主成分上，改变投影为：
projected = array([dot(V[[0,2]],immatrix[i]-immean) for i in range(imnbr)])
以得到相应的坐标（在这里 V[[0,2]] 分别是第一个和第三个主成分）。当然，也可以将其投影到所有成分上，之后挑选出需要的列。
用PIL中的ImageDraw模块进行可视化。用下面的脚本可以生成如图所示的效果：
from PCV.tools import imtools, pca
from PIL import Image, ImageDraw
from pylab import *
from PCV.clustering import hcluster
from sympy import floor
imlist = imtools.get_imlist('selectedfontimages/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)])
# height and width
h, w = 1200, 1200
# create a new image with a white background
img = Image.new('RGB', (w, h), (255, 255, 255))
draw = ImageDraw.Draw(img)
# draw axis
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 coordinates to fit
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)
# paste thumbnail of each image
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('pca_font.png')
show()
6.1.4 像素聚类
将图像区域或像素合并成有意义的部分称为图像分割。单纯在像素水平上应用 K-means可以用于一些简单图像的图像分割，但是对于复杂图像得出的结果往往是毫无意义的。要产生有意义的结果，往往需要更复杂的类模型而非平均像素色彩或空间一致性。
下面在RGB三通道的像素值上运用K-means进行聚类：
from scipy.cluster.vq import *
from scipy import *
from pylab import *
from PIL import Image
# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)
steps = 50 # image is divided in steps*steps region
infile = 'pic/empire.jpg'
im = array(Image.open(infile))
dx = im.shape[0]//steps
dy = im.shape[1]//steps
# compute color features for each region
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') # make into array
# cluster
centroids, variance = kmeans(features, 3)
code, distance = vq(features, centroids)
# create image with cluster labels
codeim = code.reshape(steps, steps)
#codeim = imresize(codeim, im.shape[:2], 'nearest')
figure()
ax1 = subplot(121)
title(u'原图', fontproperties=font)
#ax1.set_title('Image')
axis('off')
imshow(im)
ax2 = subplot(122)
title(u'聚类后的图像', fontproperties=font)
#ax2.set_title('Image after clustering')
axis('off')
imshow(codeim)
show()
6.2 层次聚类
层次聚类（或凝聚式聚类）是另一种简单但有效的聚类算法，其思想是基于样本间成对距离建立一个简相似性树。 该算法首先将特征向量距离最近的两个样本归并为一组，并在树中建立一个”平均“节点，将这两个距离最近的样本作为该”平均“节点下的子节点；然后再剩下的包含任意平均节点的样本中寻找下一个最近的对，重复进行前面的操作。在每一个节点处保存了两个子节点之间的距离。遍历整个树，通过设定的阈值，遍历过程可以在比阈值大的节点位置终止，从而提取出聚类簇。
层次聚类技术是第二类重要的聚类方法。层次聚类方法对给定的数据集进行层次的分解，直到满足某种条件为止，传统的层次聚类算法主要分为两大类算法：
- 凝聚的层次聚类： AGNES算法(AGglomerative NESting) → \rightarrow→ 采用 自底向上 的策略。
- 最初将每个对象作为一个簇，然后这些簇根据某些准则被一步一步合并，两个簇间的距离可以由这两个不同簇中距离最近的数据点的相似度来确定；聚类的合并过程反复进行直到所有的对象满足簇数目。
- 分裂的层次聚类： DIANA算法(DIvisive ANALysis) → \rightarrow→采用 自顶向下 的策略。
- 首先将所有对象置于一个簇中，然后按照某种既定的规则逐渐细分为越来越小的簇(比如最大的欧式距离)，直到达到某个终结条件(簇数目或者簇距离达到阈值)。
创建文件 hcluster.py，将下面代码添加进去：
from itertools import combinations
from pylab import *
from sympy import sqrt
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()。对于每个子树，计算其所有节点特征向量的平均值，作为新的特征向量来表示该子树，并将每个子树视为一个对象。当然，还有其他将哪两个节点合并在一起的方案，比如在两个子树中使用对象间距离最小的单向锁，及在两个子树中用对象间距离最大的完全锁。选择不同的锁会生成不同类型的聚类树。
在一个简单的例子中观察该聚类过程：
from PCV.clustering import hcluster
from numpy import *
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())
在实际图像处理中的应用：
import os
from PCV.clustering import hcluster
from matplotlib.pyplot import *
from numpy import *
from PIL import Image
from sympy import minimum
# 创建图像列表
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, 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')
用100幅日落图像进行层次聚类，将RGB空间的512个小区间直方图作为每幅图像的特征向量。树中挨得近的图像具有相似的颜色分布。树状图的高和子部分由距离决定，这些都需要调整，以适应所选择的图像分辨率。随着坐标向下传递到下一级，会递归绘制出这些节点，上述代码用20×20像素绘制叶节点的缩略图，，使用 get_height() 和 get_depth() 这两个辅助函数可以获得树的高和宽。
6.3 谱聚类
谱： 方阵作为线性算子，它的所有特征值的全体统称为方阵的谱。方阵的谱半径为最大的特征值。矩阵A AA的谱半径是矩阵\(A^TA\)的最大特征值。
谱聚类：是一种基于图论的聚类方法，通过对样本数据的拉普拉斯矩阵的特征向量进行聚类，从而达到对样本数据聚类的谱。谱聚类可以理解为将高维空间的数据映射到低维，然后在低维空间用其它聚类算法（如KMeans）进行聚类。
谱聚类的过程：
给定一个n×n的相似性矩阵S SS，sij为相似性分数，可以创建一个矩阵，称为拉普拉斯矩阵：
\(L=I-D^{-1/2}SD^{-1/2}\)
其中，\(\text{I}\)是单位矩阵，\(\text{D}\)是对角矩阵，对角线上的元素是\(\text{S}\)对应行元素之和，\(D=diag(d_i), d_i=\sum_js_{ij}\)。拉普拉斯矩阵中的\( D^{-1/2}\)为：
\(\boldsymbol{D}^{-1/2}=\begin{bmatrix}\frac{1}{\sqrt{d_1}}&&&&\\&\frac{1}{\sqrt{d_2}}&&&\\&&\ddots&&\\&&&\frac{1}{\sqrt{d_n}}\end{bmatrix}\)
为了简介表示，使用较小的值并且要求\(s_{ij}\geqslant 0\)。
计算L的特征向量，并使用k个最大特征值对应的k个特征向量，构建出一个特征向量集，从而可以找到聚类簇。创建一个矩阵，该矩阵的各列是由之前求出的k个特征向量构成，每一行可以看作一个新的特征向量，长度为k。本质上，谱聚类算法是将原始空间中的数据转换成更容易聚类的新特征向量。在某些情况下，不会首先使用聚类算法。
谱聚类算法的代码:
# -*- 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('fontimages\\a_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)])
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)): # type: ignore
im = Image.open(imlist[ind[i]])
subplot(4, 10, i + 1)
imshow(array(im))
axis('equal')
axis('off')
show()
在上面的实验中，用两两间的欧式距离创建矩阵S，并对k个特征向量用常规的K-means进行聚类。注意，矩阵V包含的是对特征值进行排序后的特征向量。然后，绘制出这些聚类簇。观察到，上面分别显示出了五类，根据不同的特征向量，将相同的类聚集起来，形成这些聚类图像。
21万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
