图像聚类
- 1 K-means聚类
- 1.1 Scipy聚类包
- 1.2 图像聚类
- 1.3 在主成分上可视化图像
- 1.4 像素聚类
- 2 层次聚类
- 3 谱聚类
本章将介绍几种聚类方法，并展示如何利用他们对图像进行聚类，从而寻找相似的图像组。聚类可以用于识别、划分图像数据集，组织和导航。此外，我们还会对聚类后的图像进行相似性可视化。
1 K-means聚类
K-means是一种将输入数据划分成k个簇的简单聚类算法。K-means反复提炼初始评估的类中心，步骤如下：
1）以随机或猜测的方式初始化类中心，i=1…k；
2）将每个数据点归并到离他距离最近的类中心所属的类；
3）对所有属于该类的数据点球平均，将平均值作为新的类中心；
4）重复步骤2）和步骤3）直到收敛。
K-means试图使类内总方差最小：xj是输入数据，并且是矢量。该算法是启发式提炼算法，在很多情况下都适用，但是并不能保证得到最优的结果。为了避免初始化类中心是没选取好类中心初值所造成的的影响，该算法通常会初始化不同的类中心进行多次运算，然后选择方差V最小的结果。
K-means算法最大的缺陷是必须余弦设定聚类数k，如果选择不恰当则会导致聚类出来的结果很差。其优点是很容易实现，可以并行计算，并且对于很多别的问题不需要任何调整就能够直接使用。
1.1 Scipy聚类包
尽管K-means算法很容易实现，但我们没有必要自己实现它，Scipy矢量量化包scipy.cluster.vq中有K-means的实现，下面是使用方法：
from pylab import *
from scipy.cluster.vq 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, 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()
1.2 图像聚类
我们用K-means对字体图像进行聚类。文件selectedfontimage.zip包含66幅来自该字体数据集fontimages的图像，用pickle模块载入模型文件，在主成分上对图像进行投影，然后用下面的方法聚类：
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
# 获取selected-fontimages 文件下图像文件名，并保存在列表中
imlist = imtools.get_imlist('D:\\123\图像处理\Image Processing\Image Processing\Chapter6\selectedfontimages\\a_selected_thumbs')
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)
# 保存均值和主成分
# f = open('D:\\123\图像处理\Image Processing\Image Processing\Chapter6\\font_pca_modes.pkl','wb')
# pickle.dump(immean, f)
# pickle.dump(V, f)
# f.close()
# 载入模型文件
with open('D:\\123\图像处理\Image Processing\Image Processing\Chapter6\\font_pca_modes.pkl','rb') as f:
immean = pickle.load(f)
V = pickle.load(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()
1.3 在主成分上可视化图像
为了便于观察上面是如何利用主成分进行聚类的，我们可以在一对主成分方向的坐标上可视化这些图像。一种方法是将图像投影到两个主成分上，改变投影为：
projected = array([dot(V[[0,2]],immatrix[i]-immean) for i in range(imnbr)]以得到相应的坐标（在这里V[[0,2]]分别是第一个和第三个主成分）。当然，你也可以将其投影到所有成分上，之后挑选出你需要的列：
from PCV.tools import imtools, pca
from PIL import Image, ImageDraw
from pylab import *
from PCV.clustering import hcluster
imlist = imtools.get_imlist('D:\\123\图像处理\Image Processing\Image Processing\Chapter6\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)])
# 高和宽
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:\\123\图像处理\Image Processing\Image Processing\Chapter6\\pca_font.png')
show()
1.4 像素聚类
在结束本节之前，我们来看一个对单幅图像中的像素而非全部图像进行聚类的例子。将图像区域或像素合并成有意义的部分称为图像分割。单纯在像素水平上应用 K-means可以用于一些简单图像的图像分割，但是对于复杂图像得出的结果往往是毫无意义的。要产生有意义的结果，往往需要更复杂的类模型而非平均像素色彩或空间一致性。
from scipy.cluster.vq import *
from scipy import *
from pylab import *
from PIL import Image
steps = 50 # 图像被划分成steps×steps的区域
im = array(Image.open('D:\\123\图像处理\Image Processing\Image Processing\Chapter6\empire.jpg'))
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()
从结果可以看出像素聚类还是具有一定的可行性的，算法将图像分为了天空，楼房的向阳面和向阴面。尽管利用窗口对图像进行了下采样，但结果仍然是由噪声的。如果图像某些区域没有空间一致性，则很难将其分开，如图中向阳面中的屋檐处。空间一致性和更好的分割方法以及其他的图像分割算法会在后面讨论，现在，让我们继续来看下一个基本的聚类算法。
2 层次聚类
层次聚类（或凝聚式聚类）是另一种简单但有效的聚类算法，其思想是基于样本间成对距离建立一个简相似性树。该算法首先将特征向量距离最近的两个样本归并为一组，并在树中创建一个“平均节点”，将这两个距离最近的样本作为该“平均”节点下的子节点；然后在剩下的包含任意平均节点的样本中寻找下一个最近的对，重复进行前面的操作。在每一个节点处保存了两个子节点之间的距离。遍历整个树，通过设定的阈值，遍历过程可以在比阈值大的节点位置终止，从而提取出聚类簇。
层次聚类由若干优点，例如，利用树结构可以可视化数据间的关系，并显示这些簇是如何关联的。在树中，一个好的特征向量可以给出一个很好的分离结果。另外一个优点是，对于给定的不同的阈值，可以直接利用原来的树，而不需要重新计算。不足之处是，对于实际需要的聚类簇，我们需要选择一个合适的阈值。
创建文件hcluster.py，将下面代码添加进去：
from itertools import combinations
from pylab import *
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()。对于每个子树，计算其所有节点特征向量的平均值，作为新的特征向量来表示该子树，并将每个子树视为一个对象。当然，还有其他将哪两个节点合并在一起的方案，比如在两个子树中使用对象间距离最小的单向锁，及在两个子树中用对象间距离最大的完全锁。选择不同的锁会生成不同类型的聚类树。
下面，我们在一个简单例子中观察该聚类过程：
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())
全连接的凝聚层次聚类的操作步骤：
1、获取所有样本的距离矩阵
2、将每个数据点作为一个单独的簇
3、基于最不相似(距离最远)样本的距离，合并两个最接近的簇
4、更新样本的距离矩阵
5、重复2到4，直到所有样本都属于同一个簇为止。
下面编写代码实现上面的操作：
1、获取样本
随机产生5个样本，每个样本包含3个特征(x，y，z)
import pandas as pd
import numpy as np
np.random.seed(1)
# 设置特征名称
variable = ["x","y","z"]
# 设置标号
labels = ["s1","s2","s3","s4","s5"]
# 产生一个（5，3）的数组
data = np.random.random_sample([5,3])*10
# 通过pandas将数组转化成一个DataFrame
df = pd.DataFrame(data,columns=variable,index=labels)
# 查看数据
print(df)
2、获取所有样本的距离矩阵
通过SciPy来计算距离矩阵，计算每个样本间两两的欧式距离，将矩阵矩阵用一个DataFrame进行保存，方便查看：
from scipy.spatial.distance import pdist,squareform
# 获得距离矩阵
"""
pdist:计算两两样本间的欧式距离，返回的是一个一维数组
squareform：将数组转换成一个对称矩阵
"""
dist_matrix = pd.DataFrame(squareform(pdist(df,metric="euclidean")),
columns=labels,index=labels)
print(dist_matrix)
3、获取全连接矩阵的关联矩阵
通过scipy的linkage函数，获取一个以全连接作为距离判定标准的关联矩阵(linkage matrix)
from scipy.cluster.hierarchy import linkage
# 以全连接作为距离判断标准，获取一个关联矩阵
row_clusters = linkage(dist_matrix.values,method="complete",metric="euclidean")
#将关联矩阵转换成为一个DataFrame
clusters = pd.DataFrame(row_clusters,columns=["label 1","label 2","distance","sample size"],
index=["cluster %d"%(i+1) for i in range(row_clusters.shape[0])])
print(clusters)
4、通过关联矩阵绘制树状图
使用scipy的dendrogram来绘制树状图
from scipy.cluster.hierarchy import dendrogram
import matplotlib.pyplot as plt
row_dendr = dendrogram(row_clusters,labels=labels)
plt.tight_layout()
plt.ylabel("欧式距离")
plt.show()
通过上面的树状图，可以直观的发现。首先是s1和s5合并，s2和s3合并，然后s2、s3、s4合并，最后再和s1、s5合并。
我们下面来看一个基于图像颜色信息对图像进行聚类的例子。文件sunsets.zip中包含100张图像。在这个例子中，我们用颜色直方图作为每幅图像的特征向量。虽然这样处理有些简单粗糙，但仍然能够很好地说明分层聚类地过程。在包含这些日落图像的文件夹中运行下面的代码：
import os
from PCV.clustering import hcluster
from matplotlib.pyplot import *
from numpy import *
from PIL import Image
# 创建图像列表
path = 'D:\\123\图像处理\Image Processing\Image Processing\Chapter6\\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.jpg')
3 谱聚类
谱聚类方法是一种有趣的聚类算法，与前面的K-means和层次聚类方法截然不同。
对于n个元素（如n幅图像），相似矩阵（或亲和矩阵，有时也称为距离矩阵）是一个n×n的矩阵，矩阵每个元素表示两两之间的相似性分数。谱聚类是由相似性矩阵构建谱矩阵而得名的。对该谱矩阵进行特征分解得到的特征向量可以用于降维，然后聚类。
谱聚类的优点之一是仅需输入相似性矩阵，并且可以采用你所想到的任意度量方式构建该相似性矩阵。像K-means和层次聚类需要计算那些特征向量求平均；为了计算平均值，会将特征或描述子限制为向量。而对于谱方法，特征向量就没有类别限制，只要有一个“距离”或“相似性”的概念即可。
下面说明谱聚类的过程。给定一个n×n的相似矩阵S，sij为相似性分数，我们可以创建一个矩阵，称为拉普拉斯矩阵：其中，I是单位矩阵，D是对角矩阵，对角线上的元素是S对应行元素之和，D=diag(di)，di=
sij。拉普拉斯矩阵的D-1/2为：
为了简洁表示，使用较小的值并且要求sij⩾0。
计算L的特征向量，并使用k个最大特征值对应的k个特征向量，构建出一个特征向量集，从而可以找到聚类簇。创建一个矩阵，该矩阵的各列是由之前求出的k个特征向量构成，每一行可以看作一个新的特征向量，长度为k。本质上，谱聚类算法是将原始空间中的数据转换成更容易聚类的新特征向量。在某些情况下，不会首先使用聚类算法。
我们来看看真实的例子中谱聚类算法的代码，使用K-means例子中的字体图像：
from PCV.tools import imtools, pca
from PIL import Image, ImageDraw
from pylab import *
from scipy.cluster.vq import *
# 获取selected-fontimages 文件下图像文件名，并保存在列表中
imlist = imtools.get_imlist('D:\\123\图像处理\Image Processing\Image Processing\Chapter6\selectedfontimages\\a_selected_thumbs')
imnbr = len(imlist)
# 创建矩阵，存储所有拉成一组形式后的图像
immatrix = array([array(Image.open(im)).flatten() for im in imlist],'f')
# PCA降维
V, S, immean = pca.pca(immatrix)
# 投影到前两个主成分上
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=1)
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()
实验中，用两两间的欧式距离创建矩阵S，并对k个特征向量用常规的K-means进行聚类。注意，矩阵V包含的是对特征值进行排序后的特征向量。然后，绘制出这些聚类簇。观察到，上面分别显示出了五类，根据不同的特征向量，将相同的类聚集起来，形成这些聚类图像。
