目录
一、K-means聚类
K-means 是一种将输入数据划分成 k 个簇的简单的聚类算法。K-means 会反复提炼初始评估的类中心,其步骤为:
- 以随机或猜测的方式初始化类中心
,i=1…k
- 将每个数据点归并到离它距离最近的类中心所属的类
- 对所有属于该类的数据点求平均,将平均值作为新的类中心
- 重复步骤(2)和步骤(3)直到收敛
K-means 试图使类内总方差最小:
该算法是启发式提炼算法,在很多情形下都适用,但并不能保证得到最优的结果。K-means 算法的缺点是必须预先设定聚类数 k,如果选择不恰当则会导致聚出来的结果很差。其优点是容易实现,可以并行计算,并且对于很多别的问题不需要任何调整就能够直接使用。
1.SciPy聚类包
SciPy 矢量量化包 scipy.cluster.vq 中有 K-means 的实现方法,其使用方法为:
from matplotlib import pyplot as plt
import numpy as np
from scipy import randn
from scipy.cluster.vq import *
class1 = 1.5 * randn(100,2)
class2 = randn(100,2) + np.array([5,5])
features = np.vstack((class1,class2))
centroids,variance = kmeans(features,2) #聚类
code,distance = vq(features,centroids) #归类
#可视化
plt.figure()
ndx = np.where(code==0)[0]
plt.plot(features[ndx,0],features[ndx,1],'*')
ndx = np.where(code==1)[0]
plt.plot(features[ndx,0],features[ndx,1],'r.')
plt.plot(centroids[:,0],centroids[:,1],'go')
plt.axis('off')
plt.show()
得出的结果为:
该实验对其类中心标记为绿色大圆环,预测出的类分别标记为蓝色星号和红色点。
2.图像聚类
import os
from PIL import Image
import imtools
import pickle
from matplotlib import pyplot as plt
import numpy as np
from scipy.cluster.vq import *
def get_imlist(path):
return [os.path.join(path,f) for f in os.listdir(path) if f.endswith('.jpg')]
imglist = get_imlist('D:/BaiduNetdiskDownload/a_selected_thumbs/')
imnbr = len(imglist)
# 载入模型文件
with open('font_pca_modes.pkl','rb') as f:
immean = pickle.load(f)
V = pickle.load(f)
immatrix = np.array([np.array(Image.open(im)).flatten() for im in imglist],'f')
immean = immean.flatten()
projected = np.array([np.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 = np.where(code==k)[0]
plt.figure()
plt.gray()
for i in range(np.minimum(len(ind),40)):
plt.subplot(4,10,i+1)
plt.imshow(immatrix[ind[i]].reshape((25,25)))
plt.axis('off')
plt.show()
结果为:
第一张图是使用主成分数目为20,聚类数k为4所得得结果;第一张图是使用主成分数目为20,聚类数k为6所得得结果;第一张图是使用主成分数目为30,聚类数k为4所得得结果。
3.在主成分上可视化图像
可以将上述的结果可视化显示到一个坐标轴上面,其实现代码为:
projected = np.array([np.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 = np.floor(np.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(imglist[i])
nodeim.thumbnail((25,25))
ns = nodeim.size
# import pdb;pdb.set_trace()
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)))
plt.imshow(img)
plt.axis('off')
plt.show()
可视化得结果为:
第一张图为使用第一个和第二个主成分得到的结果;第二张图为使用第二个和第三个主成分得到的结果。
4.像素聚类
其次还可以对单幅图像中的像素而非全部图像进行聚类。将图像区域或像素合并成有意义的部分称为图像分割。下面是对一幅图像进行像素聚类得例子:
from PIL import Image
from matplotlib import pyplot as plt
import numpy as np
from scipy.cluster.vq import *
# from scipy.misc import imresize
steps = 50 # 图像被划分成 steps×steps 的区域
im = np.array(Image.open(r'D:\test\pil.png'))
dx = im.shape[0] / steps
dy = im.shape[1] / steps
# 计算每个区域的颜色特征
features = []
for x in range(steps):
for y in range(steps):
# import pdb; pdb.set_trace()
R = np.mean(im[int(x)*int(dx):(int(x)+1)*int(dx),int(y)*int(dy):(int(y)+1)*int(dy),0])
G = np.mean(im[int(x)*int(dx):(int(x)+1)*int(dx),int(y)*int(dy):(int(y)+1)*int(dy),1])
B = np.mean(im[int(x)*int(dx):(int(x)+1)*int(dx),int(y)*int(dy):(int(y)+1)*int(dy),2])
features.append([R,G,B])
features = np.array(features,'f') # 变为数组
# 聚类
centroids,variance = kmeans(features,3)
code,distance = vq(features,centroids)
# 用聚类标记创建图像
codeim = code.reshape(steps,steps)
codeim = np.array(Image.fromarray(codeim).resize(im.shape[:2]))
# codeim = imresize(codeim,im.shape[:2],interp='nearest')
plt.figure()
plt.subplot(121)
plt.imshow(im)
plt.axis('off')
plt.subplot(122)
plt.imshow(codeim)
plt.axis('off')
plt.show()
遇到了如下的错误:
原因在于imresize已经被新版的scipy遗弃了,可以使用numpy来进行替代。
结果为:
上面的第一张图是使用k=3,窗口大小为50*50进行聚类得到的结果;下图是将窗口大小修改为100*100得到的结果。
二、层次聚类
1.基本概念
层次聚类是另一种简单但有效的聚类算法,其思想是基于样本间成对距离建立一个简相似性树。层次聚类可以利用树结构可以可视化数据间的关系,并显示这些簇是如何关联的。其次对于给定的不同的阈值,可以直接利用原来的树,而不需要重新计算。缺点在于需要选择一个合适的阈值来应对实际需要的聚类簇。其实现方法为:
from itertools import combinations
import numpy as np
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):
if self.distance < dist:
return [self]
return self.left.extract_clusters(dist) + self.right.extract_clusters(dist)
def get_cluster_elements(self):
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 np.sqrt(sum((v1-v2)**2))
def L1dist(v1,v2):
return sum(abs(v1-v2))
def hcluster(features,distfcn=L2dist):
distances = {}
# 每行初始化为一个簇
node = [ClusterLeafNode(np.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]
from scipy import randn
import hcluster
import numpy as np
class1 = 1.5 * randn(100,2)
class2 = randn(100,2) + np.array([5,5])
features = np.vstack((class1,class2))
tree = hcluster.ClusterLeafNode.hcluster(features)
clusters = tree.extract_clusters(5)
print ('number of clusters', len(clusters))
for c in clusters:
print (c.get_cluster_elements())
结果为:
该代码首先创建一个包含叶节点的列表,然后根据选择的距离度量方式将距离最近的对归并到一起,返回的终节点即为树的根。距离度量的选择依赖于实际的特征向量。为了从树中提取聚类簇,需要从顶部遍历树直至一个距离小于设定阈值的节点终止,可以通过递归来实现。
2.图像聚类
接着就是将其进阶为对图像进行聚类操作
import os
from PIL import Image
from matplotlib import pyplot as plt
import numpy as np
import hcluster
# 创建图像列表
path = 'D:/BaiduNetdiskDownload/flickr-sunsets-small/'
imlist = [os.path.join(path,f) for f in os.listdir(path) if f.endswith('.jpg')]
# 提取特征向量,每个颜色通道量化成 8 个小区间
features = np.zeros([len(imlist), 512])
for i,f in enumerate(imlist):
im = np.array(Image.open(f))
# 多维直方图
h,edges = np.histogramdd(im.reshape(-1,3),8,normed=True,
range=[(0,255),(0,255),(0,255)])
features[i] = h.flatten()
tree = hcluster.ClusterLeafNode.hcluster(features)
其使用了RGB三个颜色通道作为特征向量。使用了树状图来可视化聚类树,其实现代码为:
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()
最终可以得到结果为:
三、谱聚类
谱聚类方法是一种有趣的聚类算法,谱聚类是由相似性矩阵构建谱矩阵而得名的。对该谱矩阵进行特征分解得到的特征向量可以用于降维,然后聚类。相似矩阵是一个 n×n 的矩阵,矩阵每个元素表示两两之间的相似性分数。
其优点在于仅需输入相似性矩阵。其过程为给定一个 n×n 的相似矩阵 S,为相似性分数,那么可以得到一个拉普拉斯矩阵为:
其中的为:
计算 L 的特征向量,并使用 k 个最大特征值对应的 k 个特征向量,构建出一个特征向量集,从而可以找到聚类簇。在本质上,谱聚类算法是将原始空间中的数据转换成更容易聚类的新特征向量。其实现方法为:
import pickle
import os
from PIL import Image
from matplotlib import pyplot as plt
import numpy as np
from scipy.cluster.vq import *
def get_imlist(path):
return [os.path.join(path,f) for f in os.listdir(path) if f.endswith('.jpg')]
imglist = get_imlist('D:/BaiduNetdiskDownload/a_selected_thumbs/')
imnbr = len(imglist)
# 载入模型文件
with open('font_pca_modes-34.pkl','rb') as f:
immean = pickle.load(f)
V = pickle.load(f)
immatrix = np.array([np.array(Image.open(im)).flatten() for im in imglist],'f')
immean = immean.flatten()
projected = np.array([np.dot(V[:30],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 = np.array([[ np.sqrt(sum((projected[i]-projected[j])**2)) for i in range(n) ] for j in range(n)], 'f')
# import pdb; pdb.set_trace()
# 创建拉普拉斯矩阵
rowsum = np.sum(S,axis=0)
D = np.diag(1 / np.sqrt(rowsum))
I = np.identity(n)
L = I - np.dot(D,np.dot(S,D))
# 计算矩阵 L 的特征向量
U,sigma,V = np.linalg.svd(L)
k = 5
features = np.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 = np.where(code==c)[0]
plt.figure()
for i in range(np.minimum(len(ind),39)):
im = Image.open(imglist[ind[i]])
plt.subplot(4,10,i+1)
plt.imshow(np.array(im))
plt.axis('equal')
plt.axis('off')
plt.show()
结果为:
其中用两两间的欧式距离创建矩阵 S,并对 k 个特征向量用常规的 K-means 进行聚类,K值为5。其次也可以在没有任何特征向量或没有严格定义相似性的例子使用该算法。