一、K-means聚类
K-means是一种将输入数据划分成k个簇的简单的聚类算法。K-means反复提炼初始评估的类中心,步骤如下:
- 以随机或猜测的方式初始化类中心u,,i=1…k;
- 将每个数据点归并到离它距离最近的类中心所属的类c;
- 对所有属于该类的数据点求平均,将平均值作为新的类中心;
- 重复步骤(2)和步骤(3)直到收敛。
K-means 试图使类内总方差最小:
1 SciPy聚类包
为便于说明,我们先生成简单的二维数据:
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()
运行结果:类中心标记为绿色大圆环,预测出的 类分别标记为蓝色星号和红色点。
2 图像聚类
我们利用之前计算过的前 40 个主 成分进行投影,用投影系数作为每幅图像的向量描述符。用 pickle 模块载入模型文 件,在主成分上对图像进行投影,然后用下面的方法聚类:
import imtools
import pickle
from scipy.cluster.vq import *
from PIL import Image
from pylab import *
# 获取 selected-fontimages 文件下图像文件名,并保存在列表中
imlist = imtools.get_imlist('E:/data/a_selected_thumbs')
imnbr = len(imlist)
# 载入模型文件
with open('../font_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()
运行结果:
3在主成分上可视图像
改变投影为:
projected = array([dot(V[[0,2]],immatrix[i]-immean) for i in range(imnbr)])
用
PIL
中的
ImageDraw
模块进行可视化:
# 高和宽
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]))
# 粘贴每幅图像的缩略图到白色背景图片
for i in range(imnbr):
nodeim = Image.open(imlist[i])
nodeim.thumbnail((25,25))
ns = nodeim.size
img.paste(nodeim,(int(scaled[i][0]-ns[0]//2),int(scaled[i][1]-ns[1]//2),int(scaled[i][0]+ns[0]//2+1),int(scaled[i][1]+ns[1]//2+1)))
img.save('pca_font.jpg')
在成对主成分上投影的字体图像。左图用的是第一个和第二个主成分,右图用的是
第二个和第三个主成分:
4 像素聚类
下面的代码示例载入一幅图像,用一个步长为
steps
的方形网格在图像中滑动,每
滑一次对网格中图像区域像素求平均值,将其作为新生成的低分辨率图像对应位置
处的像素值,并用
K-means
进行聚类:
from PIL import Image, ImageDraw
from pylab import *
from scipy.cluster.vq import *
# from scipy.misc import imresize
steps = 50 # 图像被划分成 steps×steps 的区域
im = array(Image.open('E:/data/empire.jpg'))
dx = int(im.shape[0] / steps)
dy = int(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],interp='nearest')
# im = imresize(im, (h, int(w * aspect_ratio)), interp='bicubic')
codeim = np.array(Image.fromarray(codeim).resize(im.shape[:2]))
figure()
imshow(codeim)
show()
二、层次聚类
层次聚类(或凝聚式聚类)是另一种简单但有效的聚类算法,其思想是基于样本间 成对距离建立一个简相似性树。该算法首先将特征向量距离最近的两个样本归并为 一组,并在树中创建一个“平均”节点,将这两个距离最近的样本作为该“平均” 节点下的子节点;然后在剩下的包含任意平均节点的样本中寻找下一个最近的对, 重复进行前面的操作。在每一个节点处保存了两个子节点之间的距离。遍历整个树, 通过设定的阈值,遍历过程可以在比阈值大的节点位置终止,从而提取出聚类簇。
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:
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]
运行结果:
number of clusters 2
[168, 192, 128, 145, 143, 149, 116, 144, 130, 199, 118, 123, 103, 181, 165, 193, 117, 159, 190, 135, 174, 109, 108, 184, 137, 171, 115, 177, 189, 197, 146, 166, 161, 183, 141, 152, 179, 107, 195, 125, 188, 111, 127, 140, 191, 113, 100, 121, 136, 185, 194, 178, 101, 163, 169, 150, 126, 182, 142, 158, 186, 106, 151, 147, 112, 167, 132, 164, 134, 162, 120, 157, 170, 104, 129, 155, 156, 176, 133, 122, 153, 160, 198, 131, 124, 148, 105, 180, 139, 154, 102, 196, 175, 119, 114, 173, 138, 110, 187]
[52, 94, 33, 39, 90, 23, 28, 49, 69, 91, 26, 30, 71, 89, 66, 3, 27, 76, 87, 93, 96, 172, 61, 59, 84, 72, 46, 43, 50, 15, 42, 55, 78, 16, 14, 79, 88, 63, 29, 56, 95, 8, 44, 57, 10, 45, 2, 7, 1, 19, 62, 0, 99, 17, 81, 18, 80, 20, 74, 6, 98, 12, 13, 75, 21, 67, 9, 25, 54, 85, 41, 51, 64, 31, 73, 22, 32, 92, 58, 83, 35, 77, 86, 4, 70, 40, 65, 37, 82, 97, 48, 5, 47, 11, 24, 36, 38, 68, 34, 53, 60]
图像聚类
def draw_dendrogram(node, imlist, filename='clusters.jpg'):
""" 绘制聚类树状图,并保存到文件中 """
# 高和宽
rows = node.get_height() * 20
cols = 1200
# 距离缩放因子,以便适应图像宽度
s = float(cols - 150) / node.get_depth()
# 创建图像,并绘制对象
im = Image.new('RGB', (cols, rows), (255, 255, 255))
draw = ImageDraw.Draw(im)
# 初始化树开始的线条
draw.line((0, rows / 2, 20, rows / 2), fill=(0, 0, 0))
# 递归地画出节点
node.draw(draw, 20, (rows / 2), s, imlist, im)
im.save(filename)
im.show()
#画树状图时对于每个节点用了 draw() 方法
def draw(self, draw, x, y, s, imlist, im):
""" 用图像缩略图递归地画出叶节点 """
h1 = int(self.left.get_height() * 20 / 2)
h2 = int(self.right.get_height() * 20 / 2)
top = y - (h1 + h2)
bottom = y + (h1 + h2)
# 子节点垂直线
draw.line((x, top + h1, x, bottom - h2), fill=(0, 0, 0))
# 水平线
ll = self.distance * s
draw.line((x, top + h1, x + ll, top + h1), fill=(0, 0, 0))
draw.line((x, bottom - h2, x + ll, bottom - h2), fill=(0, 0, 0))
# 递归地画左边和右边的子节点
self.left.draw(draw, x + ll, top + h1, s, imlist, im)
self.right.draw(draw, x + ll, bottom - h2, s, imlist, im)
#在画实际图像缩略图时,叶节点有自己的方法
def draw(self, draw, x, y, s, imlist, im):
nodeim = Image.open(imlist[self.id])
nodeim.thumbnail([20, 20])
ns = nodeim.size
im.paste(nodeim, [int(x), int(y - ns[1] // 2), int(x + ns[0]), int(y + ns[1] - ns[1] // 2)])
用
100
幅日落图像进行层次聚类的示例聚类簇,阈值集合设定为树中最大节点距离的 23%
三、谱聚类
谱聚类方法是一种有趣的聚类算法。
对于n个元素(如n幅图像),相似矩阵(或亲和矩阵,有时也称距离矩阵)是一个n×n的矩阵,矩阵每个元素表示两两之间的相似性分数。谱聚类是由相似性矩阵构建谱矩阵而得名的。对该谱矩阵进行特征分解得到的特征向量可以用于降维,然后聚类。
下面说明谱聚类的过程。给定一个n×n的相似矩阵S,为相似性分数,我们可以创建一个矩阵,称为拉普拉斯矩阵:
拉普拉斯矩阵中的为:
import pickle
from scipy.cluster.vq import *
from PIL import Image
from pylab import *
# 获取 selected-fontimages 文件下图像文件名,并保存在列表中
from py_cvp import imtools
imlist = imtools.get_imlist('E:/data/a_selected_thumbs')
imnbr = len(imlist)
# 载入模型文件
with open('../font_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)])
projected = array([dot(V[[1,2]],immatrix[i]-immean) for i in range(imnbr)])
# 进行 k-means 聚类
projected = whiten(projected)
centroids,distortion = kmeans(projected,4)
code,distance = vq(projected,centroids)
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()
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()