【python计算机视觉编程——6.图像聚类】

6.图像聚类

6.1 K-means聚类

from scipy.cluster.vq import *
import numpy as np
from pylab import *
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['font.sans-serif'] = ['SimHei']  # 黑体字体
matplotlib.rcParams['axes.unicode_minus'] =False  # 显示负号

class1 = 1.5 * np.random.randn(100,2)  # 生成一个100x2的数组，每个元素都是从均值为0、标准差为1.5的正态分布中抽取的。
class2 = np.random.randn(100,2)+np.array([5,5])  # 生成一个100x2的数组，每个元素都是从均值为5、标准差为1的正态分布中抽取的，并且整体上移[5, 5]。
features=np.vstack((class1,class2))  # 将这两个生成的类别数据合并成一个200x2的数组。

# 对features数据集进行k-均值聚类，指定为2个簇。返回簇的质心和每个簇的方差。
# 是一个形状为(2, 2)的数组，其中每一行表示一个簇的质心(中心点)。由于指定了2个簇，centroids会有2行，每行包含2个值(对应每个簇的二维坐标)。
# variance:表示数据点到其对应簇的质心的总方差。方差越小，说明数据点与质心的距离越近，聚类效果越好。
centroids,variance=kmeans(features,2)

# 将每个数据点分配给离其最近的质心，并计算每个点到其分配的质心的距离。code包含簇的分配信息，distance包含距离。
code,distance=vq(features,centroids)

figure()
ndx=where(code==0)[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')
show()

6.1.2 图像聚类

import pickle
from scipy.cluster.vq import *
import os
from PIL import Image
from pylab import *

def get_imlist(path):        
    return [os.path.join(path,f) for f in os.listdir(path) if f.endswith('.jpg')]

imlist=get_imlist(r'D:\pyFile\Python计算机视觉编程\data\selectedfontimages\a_selected_thumbs')
imnbr=len(imlist)

with open('a_pca_modes.pkl','rb') as f:
    immean=pickle.load(f,encoding='iso-8859-1')  # 从文件中加载图像的均值向量
    V=pickle.load(f,encoding='iso-8859-1')  # 从文件中加载PCA变换的特征向量（主成分）

    
# 将每个图像打开并展平为一维向量，然后将这些向量堆叠成一个矩阵。每行代表一张图像的特征向量
immatrix=np.array([array(Image.open(im)).flatten()
                  for im in imlist],'f')


# 将图像均值向量展平为一维向量
immean=immean.flatten()
# 将每张图像从原始空间投影到PCA空间。V[:40](取前40行)表示使用前40个主成分来进行投影
# dot(V[:40],immatrix[i]-immean).shape  :  (40,)
# imnbr:66
projected=np.array([dot(V[:40],immatrix[i]-immean) for i in range(imnbr)])
# projected.shape:(66, 40)


# 白化处理是一种数据预处理技术,旨在将数据的特征转换为具有相同方差的标准正态分布(均值为0,方差为1)
projected=whiten(projected)  # 对投影后的数据进行白化处理，使得每个特征的方差为1
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)): # 对簇中的前 40 张图像进行迭代（如果簇中的图像少于 40 张，则显示所有图像）
        subplot(4,10,i+1)
        imshow(immatrix[ind[i]].reshape((25,25)))
        axis('off')
show()

类1

类2

类3

类4

6.1.3 在主成分上可视化图像

上述代码中的一句修改如下，将图像投影到两个主成分上，使用 PCA 的前两个主成分（V[[0,2]]）将每张图像从原始空间投影到降维后的空间。每个图像的投影是减去均值后的特征向量与主成分的点积。

projected=np.array([dot(V[:2],immatrix[i]-immean) for i in range(imnbr)])

类1

类2

类3

类4

from PIL import Image,ImageDraw
h,w=1200,1200  # 定义背景图像的高度和宽度为 1200 像素

img=Image.new('RGB',(w,h),(255,255,255))  # 创建一个新的 RGB 图像，背景颜色为白色
draw=ImageDraw.Draw(img)  #创建一个可以在图像上绘制的 ImageDraw 对象

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)  # 计算 PCA 投影数据的绝对值的最大值，以用于数据的缩放
# 将 PCA 投影数据按比例缩放，并将其平移到图像的中心。floor 函数将坐标值向下取整，以便图像能正确地绘制在整数像素上
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
    
#    将图像粘贴到背景图的指定位置。位置根据图像的缩放和 PCA 投影位置计算得出
    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)))

imshow(img)
img.save('pca_font.jpg')

6.1.4 像素聚类

from scipy.cluster.vq import *

steps=50
im=array(Image.open('sun.jpg'))
figure()
subplot(121)
imshow(im)
dx=im.shape[0]/steps
dy=im.shape[1]/steps

features=[]
for x in range(steps):
    for y in range(steps):
        R=mean(im[int(x*dx):int((x+1)*dx),int(y*dy):int((y+1)*dy),0])
        G=mean(im[int(x*dx):int((x+1)*dx),int(y*dy):int((y+1)*dy),1])
        B=mean(im[int(x*dx):int((x+1)*dx),int(y*dy):int((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=np.array(Image.fromarray(codeim).resize(im.shape[:2]))

subplot(122)
imshow(codeim)
show()

6.2 层次聚类

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):
        """ 从层次聚类树中提取距离小于 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 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]

class1=1.5*np.random.randn(100,2)
class2=np.random.randn(100,2)+np.array([5,5])
features=np.vstack((class1,class2))

tree=ClusterLeafNode.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 PIL import Image,ImageDraw
from pylab import *

path=r'D:\pyFile\Python计算机视觉编程\data\sunsets\flickr-sunsets-small'
imlist=[os.path.join(path,f) for f in os.listdir(path) if f.endswith('.jpg')]

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,density=True,range=[(0,255),(0,255),(0,255)])
    features[i]=h.flatten()
tree=ClusterLeafNode.hcluster(features)

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函数到ClusterNode类中

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

        l1=self.distance*s
        draw.line((x,top+h1,x+l1,top+h1),fill=(0,0,0))
        draw.line((x,bottom-h2,x+l1,bottom-h2),fill=(0,0,0))

        self.left.draw(draw,x+l1,top+h1,s,imlist,im)
        self.right.draw(draw,x+l1,bottom-h2,s,imlist,im)

添加draw函数到ClusterLeafNode中

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

draw_dendrogram(tree,imlist,filename='sunset.pdf')

clusters=tree.extract_clusters(0.23*tree.distance)

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

6.3 谱聚类

from scipy.cluster.vq import *
import pickle
import os
from PIL import Image
from pylab import *

imlist=get_imlist(r'D:\pyFile\Python计算机视觉编程\data\selectedfontimages\a_selected_thumbs')
imnbr=len(imlist)

with open('a_pca_modes.pkl','rb') as f:
    immean=pickle.load(f,encoding='iso-8859-1')
    V=pickle.load(f,encoding='iso-8859-1')

immatrix=np.array([array(Image.open(im)).flatten()
                  for im in imlist],'f')

immean=immean.flatten()
projected=np.array([dot(V[:40],immatrix[i]-immean) for i in range(imnbr)])

projected=whiten(projected)
centroids,distortion=kmeans(projected,4)

code,distance=vq(projected,centroids)

n=len(projected)
S=array([[np.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))

U,sigma,V=linalg.svd(L)

k=5

features=array(V[:k]).T
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)):
#         print(imlist[ind[i]])
        im=Image.open(imlist[ind[i]])
        subplot(4,10,i+1)
        imshow(array(im))
        axis('equal')
        axis('off')
show()

类1

类2

类3

类4

类5