Python计算机视觉编程第六章——图像聚类(K-means聚类，DBSCAN聚类，层次聚类，谱聚类，PCA主成分分析)

图像聚类

聚类概念

无监督学习：没有标签。（对于监督学习问题中，我们会被告知什么是正确答案，在无监督学习中，没有任何标签，或者都具有相同的标签，得到的数据如下图，图上有一系列点，但是它们没有标签，因此训练集可以写成$\left\{x^{(1)},x^{(1)},x^{(1)},...,x^{(m)}\right\}$，没有标签y，在无监督学习中，我们要将这系列无标签的数据输入到算法中，然后让算法找到一些隐含在数据中的结构，比如这个数据集中的点，可以分成两组分开的点集（簇），这种能够找到这些簇的算法称为聚类算法）

聚类：相似的东西分到一组。
   
可以看到这个图，在原始的数据集上是没有颜色的标记，没有告诉我们蓝色是哪个，绿色是哪个，红色是哪个，我们要把相似的东西归到一边，蓝色归到左下角，粉色归到右上角，也就是说根据数据分布的不同，把相似的往一边凑，凑完之后得到三个簇，也就是三个不同的类别。

难点：有监督问题，可以在预测值和真实值之间作比较，按某种方法评估出来一个指标，但是在无监督问题中，没有标签，评估就比较困难，还有就是调参，比如第一个参数得到一个数据结果，第二个参数得到一个数据结果，那第一个参数和第二个参数哪个好呢，很难比较，因为我们没有标签，不知道这两种情况分的怎么样，没有标准答案。

（一）K-means 聚类

定义：k-means算法中的k代表类簇个数，means代表类簇内数据对象的均值（这种均值是一种对类簇中心的描述），因此，k-means算法又称为k-均值算法。k-means算法是一种基于划分的聚类算法，以距离作为数据对象间相似性度量的标准，即数据对象间的距离越小，则它们的相似性越高，则它们越有可能在同一个类簇。

首先我们必须知道这几个概念：

1.要得到簇的个数，需要指定K值（k=2，输入的数据聚成2个堆；k=3，输入的数据聚成3个堆）

2.质心：数据的均值，即各个数据在各个维取平均值得到的值（例如现在有x轴，y轴，x轴上所有数据取均值，y轴所有数据取均值，这样我们就得到了一个质心。每个簇都可以得到一个质心，质心在后面迭代的时候要用到。）

3.距离的度量：聚类简单理解就是把相似的东西聚到一起，如何判断两个样本点是不是相似的呢，这个就要根据距离做判断，最常见的计算方式就是欧几里得距离（直接算两个点的欧式距离）和余弦相似度（先标准化），当使用欧式距离的时候，先要对数据进行标准化，

       什么是标准化：比如现在有两个维度，x轴，y轴，x轴数据0.01,0.02,0.04，y轴数据100,200,300，当计算相似度的时候，x轴的差异无论怎么算都比较小，y轴的差异无论怎么算都比较大，那样我们潜意识里就认为相似度主要由y轴决定，实际算出来也是这样，所以说，在使用距离的度量的时候，基本情况下，都要对所有数据进行标准化，比如让x轴取值范围在0到1之间，y轴取值范围在0到1之间，让数据x和y基本在一个比较小的范围内浮动，比如说都是0到1，或者-1到1。

先把数据做标准化，然后再用距离的度量看一下什么样的两个样本点是相似的，再把相似的分到一簇。

4.优化目标：

$$min\sum _{i=1}^k\sum _{x\in C_i}dist(c_i,x{)}^2$$

最外层i从1到k，表示一共有多少簇，k=3的时候，表示对3个簇分别进行优化，里面一层是针对对每一个簇来说，每一个样本点到质心的距离算出来，优化目标就是每一个样本点到质点的距离求和越小。为什么越小越好，因为越小越相似。我们的目的就是把相似的东西分到一组，假如现在图像有一个红色的点，这个点到红色簇质心的距离非常大，到蓝色簇质心的距离非常小，非常大是我们不希望看到的情况，所以在迭代优化的时候，把这个红色的点化成蓝色就好了。总的来说，就是对于每一个点，要划分到合适的地方，如何看合不合适，就看它到中心的距离，比如它到红色中心的距离比较大，到蓝色中心的距离比较小，所以说它是蓝色的可能性更大，就会把它规划到蓝色簇。

工作流程：
图解：
   

   
   

   

   
讲解：

1. 由于是无监督问题，不知道每一个点该属于哪一个簇，假设k=2（ k=2,分两堆;k=3,分三堆，此处以k=2举例 ），首先会在(a)图中初始化两个点。
2. 如(b)图所示，一个黄色的，一个粉色的，这两个点是随机初始化的。然后基于初始化的这两个点，算每一个其他的点是属于粉色还是黄色，如何判断，就是算这个点到黄色点的距离（如下图），假设为$d_1$，这个点到红色点的距离，假设为$d_2$，$d_1$<$d_2$，认为当前点属于黄色，因为距离越小越相似。以此类推。把每一个样本点都算一遍。一些离黄色近的化成黄色，一些离粉色近的化成粉色。得到了( c )图。
   
3. 可以看到c图离预期效果还有一段距离。所以要进行更新，更新衡量的依据，也就是两个质心，因为这两个质心是随机选择的，肯定是不太准的。如何更新，就是把所有黄色的点拿出来，重新算质心，所有粉色的点同样算质心。得到了(d)图。
4. 质心更新完之后，重新遍历样本中所有的点，比如原来黄色的点到粉色质心的距离$d_2$更小（如下图），现在就更新到粉色。以此类推。得到了(e)图。
   
5. 再按照之前的套路，更新质心，得到(f)图。
6. 再进行遍历，找每个点的归属，然后不断的更新下去，直到更新到某一步，所有的样本点不再发生变化。不再发生变化，更新基本结束，因为再更新还是原来状态，黄色的还是属于黄色，粉色的还是属于粉色。

优势：

1. 简单（原理很简单，理解起来也很简单）
2. 快速（直接考虑k值就可以了）
3. 适合常规数据集（比如上图，两堆看起来就能分开）

劣势：

1. k值难确定（当拿到一个数据，可以做可视化展示，但是由于标签不知道是什么，所以不知道究竟分成几个簇，普遍情况下需要设置多组，然后看效果）。
         
         
         
   可以看到，好像分成几个簇都能划分出来，所以到底分成多少个簇合适，这个比较难确定。（小插曲：可能会有疑惑，为什么第八张图和第九张图与前面的分布有点不一样，是因为博主有点强迫症哈，本来想用第七张图与第六章图对比一下就可以了，结果发现七张图不好排版，但是这时候博主已经关掉了这个页面，所以第二次运行，又会重新设置初值，不同的初值得到的结果是不一样的，算是误打误撞发现新大陆吧）
2. 复杂度与样本呈现线性相关（复杂度与样本个数相关，所有的样本都要和质心计算，如果与样本过多，计算越多，复杂度就相当高了）
3. 很难发现任意形状的簇（并不能保证一个数据是一个非常规整的簇，举例，如下图，可以认为里面是一个簇，外面是一个簇，但是用K-means算法没办法把里面和外面分开，而是会分为左边一个簇，右边一个簇）
   
   下面来验证一下，这是一个笑脸，理想情况下应该分成两个眼睛，一个嘴巴，一个圆圈脸，四个簇，但是看下运行效果：
      
      
      
   一直迭代，最后得到结果：
   

下面介绍K-means算法怎么实际运用在数据处理的。详细讲解见代码注释

# beer dataset
import pandas as pd
beer = pd.read_csv('data.txt', sep=' ')   #导入数据，这是20个啤酒数据，包括四个属性
beer

X = beer[["calories","sodium","alcohol","cost"]]   #聚类算法的输入是x，x是所有的特征，聚类算法不需要构建标签

from sklearn.cluster import KMeans
km = KMeans(n_clusters=3).fit(X)     #用三个堆做聚类
km2 = KMeans(n_clusters=2).fit(X)    #用两个堆做聚类

km.labels_          #显示每个数据属于哪个类别，三个堆，所以有0,1,2

array([0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 2, 0, 0, 2, 1])

km2.labels_         #两个堆，所以有0,1

array([1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0])

beer['cluster'] = km.labels_   
beer['cluster2'] = km2.labels_
beer.sort_values('cluster')

from pandas.tools.plotting import scatter_matrix
%matplotlib inline
cluster_centers = km.cluster_centers_
cluster_centers_2 = km2.cluster_centers_
beer.groupby("cluster").mean()  #分别算三个cluster的均值，观察不同类别上的差异性

beer.groupby("cluster2").mean()   #分别算三个cluster的均值，观察不同类别上的差异性

centers = beer.groupby("cluster").mean().reset_index()    #计算中心点
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 14
import numpy as np
colors = np.array(['red', 'green', 'blue', 'yellow'])   #指定绘图颜色
plt.scatter(beer["calories"], beer["alcohol"],c=colors[beer["cluster"]])

plt.scatter(centers.calories, centers.alcohol, linewidths=3, marker='+', s=300, c='black')#中心点

plt.xlabel("Calories")
plt.ylabel("Alcohol") #展示卡路里和酒精两个维度上

scatter_matrix(beer[["calories","sodium","alcohol","cost"]],s=100, alpha=1, c=colors[beer["cluster"]], figsize=(10,10))
plt.suptitle("With 3 centroids initialized")   ##展示四个维度

from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
X_scaled          #标准化，使数据更规整

km = KMeans(n_clusters=3).fit(X_scaled)
beer["scaled_cluster"] = km.labels_
beer.sort_values("scaled_cluster")   #标准化后得到的结果

beer.groupby("scaled_cluster").mean()   #分别算归一化后三个cluster的均值，观察不同类别上的差异性

pd.scatter_matrix(X, c=colors[beer.scaled_cluster], alpha=1, figsize=(10,10), s=100)  #显示四个维度

聚类评估：轮廓系数
s ( i ) = b ( i ) − a ( i ) m a x { a ( i ) , b ( i ) } s(i)=\frac{b(i)-a(i)}{max\left \{ a(i),b(i) \right \}} s(i)=max{ a(i),b(i)}b(i)−a(i)​

$$s(i)=\{\begin{array}{cc}1-\frac{a(i)}{b(i)}& a(i)<b(i)\\ 0& a(i)=b(i)\\ \frac{b(i)}{a(i)}-1& a(i)>b(i)\end{array}$$

1. 计算样本i到同簇其他样本的平均距离a(i)。a(i)越小，说明样本i越应该被聚类到该簇。将a(i)称为样本i的簇内不相似度。
2. 计算样本i到其他某簇Cj 的所有样本的平均距离bij，称为样本i与簇Cj 的不相似度。定义为样本i的簇间不相似度：b(i) =min{bi1, bi2, …, bik}，b(i)越大越好，越大说明当前样本越不可能被分到其他簇。
3. si接近1，则说明样本i聚类合理。（ b(i)-a(i)理想情况下等于b(i)，max{a(i),b(i)}理想情况下等于b(i) ）
4. si接近-1，则说明样本i更应该分类到另外的簇
5. 若si 近似为0，则说明样本i在两个簇的边界上。

from sklearn import metrics
score = metrics.silhouette_score(X,beer.cluster)#不做归一化得到的结果 si值
score_scaled = metrics.silhouette_score(X,beer.scaled_cluster)#归一化之后得到的结果 si值
print(score_scaled, score)#发现做完归一化结果反而低了，说明归一化有时候并不能使结果变好，比如某些特征确实很重要，那做完归一化就消除了差异

0.673177504646 0.179780680894

scores = []
for k in range(2,20):    #遍历K值，看效果
    labels = KMeans(n_clusters=k).fit(X).labels_
    score = metrics.silhouette_score(X, labels)
    scores.append(score)

scores   #选择轮廓系数比较高的比较合适

plt.plot(list(range(2,20)), scores) #画图更直观展示 
plt.xlabel("Number of Clusters Initialized")
plt.ylabel("Sihouette Score")

所以当k=2的时候合适。

下面介绍K-means算法怎么实际运用在图像处理的。

# -*- coding: utf-8 -*-
# 导包
from skimage import io
from sklearn.cluster import KMeans
import numpy as np

# 读取图片
image = io.imread("D:\\Python\\chapter6\\empire.jpg")

# 显示图片
io.imshow(image)

# 获取图片压缩前的信息
print ('image的类型：', type(image))  # numpy.ndarray类型
print ('image的尺寸：', image.shape)  # 显示尺寸
print ('image的宽度：', image.shape[0])  # 图片宽度
print ('image的高度：', image.shape[1])  # 图片高度
print ('image的通道数：', image.shape[2])  # 图片通道数
print ('image的总像素个数：', image.size)  # 显示总像素个数
print ('image的最大像素值：', image.max())  # 最大像素值
print ('image的最小像素值：', image.min())  # 最小像素值
print ('image的像素平均值：', image.mean())  # 像素平均值

# rows*cols*channel = 482*500*3
rows = image.shape[0]
cols = image.shape[1]
channel = image.shape[2]

# 样本数*channel = 241000*3
# 每个样本在不同通道都有存在1个点，即此时每个样本对应于3个点
image = image.reshape(image.shape[0] * image.shape[1], channel)

# 样本数*channel = 241000*1
# 使用Kmeans算法将3通道变为1通道（将原来很多的颜色用少量的颜色来表示），即此时每个样本只存在1个点
# n_clusters：K值，把集合分成K个簇；n_init：指定CPU个数；max_iter：最大的迭代次数
kmeans = KMeans(n_clusters=128, n_init=10, max_iter=200)
kmeans.fit(image)

clusters = np.asarray(kmeans.cluster_centers_, dtype=np.uint8)
# labels_：每个点的标签
labels = np.asarray(kmeans.labels_, dtype=np.uint8)
# rows*cols*channel = 482*500*1
labels = labels.reshape(rows, cols)

np.save('D:\\Python\\chapter6\\codebook_test.npy', clusters)
# 保存压缩后的图片
io.imsave('D:\\Python\\chapter6\\compressed_test.jpg', labels)
# 读取压缩后的图片
newimage = io.imread("D:\\Python\\chapter6\\compressed_test.jpg")

# 获取图片压缩后的信息
print ('newimage的类型：', type(newimage))  # numpy.ndarray类型
print ('newimage的尺寸：', newimage.shape)  # 显示尺寸
print ('newimage的宽度：', newimage.shape[0])  # 图片宽度
print ('newimage的高度：', newimage.shape[1])  # 图片高度
# print ('newimage的通道数：', newimage.shape[2])  #图片通道数
print ('newimage的总像素个数：', newimage.size)  # 显示总像素个数
print ('newimage的最大像素值：', newimage.max())  # 最大像素值
print ('newimage的最小像素值：', newimage.min())  # 最小像素值
print ('newimage的像素平均值：', newimage.mean())  # 像素平均值

# 显示压缩后的图片
io.imshow(newimage)

控制台输出：

(‘image\xe7\x9a\x84\xe7\xb1\xbb\xe5\x9e\x8b\xef\xbc\x9a’, <type ‘numpy.ndarray’>)
(‘image\xe7\x9a\x84\xe5\xb0\xba\xe5\xaf\xb8\xef\xbc\x9a’, (800L, 569L, 3L))
(‘image\xe7\x9a\x84\xe5\xae\xbd\xe5\xba\xa6\xef\xbc\x9a’, 800L)
(‘image\xe7\x9a\x84\xe9\xab\x98\xe5\xba\xa6\xef\xbc\x9a’, 569L)
(‘image\xe7\x9a\x84\xe9\x80\x9a\xe9\x81\x93\xe6\x95\xb0\xef\xbc\x9a’, 3L)
(‘image\xe7\x9a\x84\xe6\x80\xbb\xe5\x83\x8f\xe7\xb4\xa0\xe4\xb8\xaa\xe6\x95\xb0\xef\xbc\x9a’, 1365600)
(‘image\xe7\x9a\x84\xe6\x9c\x80\xe5\xa4\xa7\xe5\x83\x8f\xe7\xb4\xa0\xe5\x80\xbc\xef\xbc\x9a’, 255)
(‘image\xe7\x9a\x84\xe6\x9c\x80\xe5\xb0\x8f\xe5\x83\x8f\xe7\xb4\xa0\xe5\x80\xbc\xef\xbc\x9a’, 0)
(‘image\xe7\x9a\x84\xe5\x83\x8f\xe7\xb4\xa0\xe5\xb9\xb3\xe5\x9d\x87\xe5\x80\xbc\xef\xbc\x9a’, 132.4720262155829)
(‘newimage\xe7\x9a\x84\xe7\xb1\xbb\xe5\x9e\x8b\xef\xbc\x9a’, <type ‘numpy.ndarray’>)
(‘newimage\xe7\x9a\x84\xe5\xb0\xba\xe5\xaf\xb8\xef\xbc\x9a’, (800L, 569L))
(‘newimage\xe7\x9a\x84\xe5\xae\xbd\xe5\xba\xa6\xef\xbc\x9a’, 800L)
(‘newimage\xe7\x9a\x84\xe9\xab\x98\xe5\xba\xa6\xef\xbc\x9a’, 569L)
(‘newimage\xe7\x9a\x84\xe6\x80\xbb\xe5\x83\x8f\xe7\xb4\xa0\xe4\xb8\xaa\xe6\x95\xb0\xef\xbc\x9a’, 455200)
(‘newimage\xe7\x9a\x84\xe6\x9c\x80\xe5\xa4\xa7\xe5\x83\x8f\xe7\xb4\xa0\xe5\x80\xbc\xef\xbc\x9a’, 167)
(‘newimage\xe7\x9a\x84\xe6\x9c\x80\xe5\xb0\x8f\xe5\x83\x8f\xe7\xb4\xa0\xe5\x80\xbc\xef\xbc\x9a’, 0)
(‘newimage\xe7\x9a\x84\xe5\x83\x8f\xe7\xb4\xa0\xe5\xb9\xb3\xe5\x9d\x87\xe5\x80\xbc\xef\xbc\x9a’, 60.489773725834795)

观察到Kmeans算法将3通道变为1通道，原来每个像素点取值是0到256，现在是0到167，从文件夹中查看图像属性发现压缩后图片的大小从215kb变成了150kb，Kmeans算法会把类似的颜色分别放在K个簇中——也就是说，每个簇的颜色都变成了一种。因此，我们只需要保留每个像素的标签（表明该像素在哪个簇中），以及每个簇的颜色编码即可完成图像的压缩。图像压缩效果比较好。

扩展：
DBSCAN聚类

基本概念：是一种基于密度的空间聚类算法。
核心对象：若某个点的密度达到算法设定的阈值则其为核心点。(即画一个圈，圈中点的个数多于阈值 ，也就是r 邻域内点的数量不小于 minPts(阈值)）
ϵ-邻域的距离阈值：设定的半径r（DBSCAN算法需要指定两个参数，一个是半径，一个是阈值，不需要设置k值，究竟聚成多少堆，由算法决定）
直接密度可达：若某点p在点q的 r 邻域内，且q是核心点则p-q直接密度可达。

密度可达：若有一个点的序列q0、q1、…qk，对任意qi-qi-1是直接密度可达的 ，则称从q0到qk密度可达（如图，q0与q1直接密度可达，q1与q2直接密度可达，q0与q2密度可达），这实际上是直接密度可达的“传播”。

密度相连：若从某核心点p出发，点q和点k都是密度可达的 ,则称点q和点k是密度相连的。
边界点: 属于某一个类的非核心点,不能发展下线了。（比如：q0作为核心点发展了q1，q1作为核心点发展了q2，但是q2里面没有点了，则q2是边界点。）

噪声点：不属于任何一个类簇的点，从任何一个核心点出发都是密度不可达的。(例如：从任何一个核心点出发都不能把q3圈进来，）

工作流程：
讲解见注释

 标记所有对象为unvisited
 Do
 随机选择一个 unvited 对象 p；
 标记 p 为visited；
 if p 的 ϵ-（半径范围内） 领域内至少有 MinPts 个对象    #如果p是核心对象
      创建一个新簇 C，并把p添加到C；
      令 N 为 p 的 ϵ- 领域中的对象集合
      for N 中每个点 p                                #假设p里面有A,B,C,D四个点，则遍历这四个点，以A举例说明
           if p 是 unvisited                          #如果p中的A没有被标记
               标记 p 为visited                        #把P中的A标记为被访问过
               if p 的ϵ-领域至少有 MinPts 个对象，把这些对象添加到N；   #一开始p发展了A,B,C,D,现在A,B,C,D继续发展下线，都添加到N,直到所有发展的下线不是核心对象了
               如果 p 还不是任何簇的成员，把 p 添加到 C；     #把p,A,B,C,D以及它们的下线都放到C中，这样第一个簇就好了
           End for；
           输出 C；
 Else 标记 p 为噪声；
 Until 没有标记为unvisited 的对象。

参数选择：

1. 半径ϵ:可以根据K距离（给定数据集P={p(i); i=0,1,…n}，计算点P(i)到集合D的子集S中所有点之间的距离，距离按照从小到大的顺序排序，d(k)就被称为k-距离。）来设定：找突变点（ P(i)到p(1)的距离d1，P(i)到p(2)的距离d2，P(i)到p(3)的距离d3，假设d1=0.1，d2=0.11，d3=0.12，而d4=0.3，d5=0.32，则d4就是突变点，根据这个突变点可以认为前面的这个距离（d=0.12）比较合适 )
   (半径大，簇的个数可能就小了，半径小，簇的个数可能就多了，对结果有影响)
2. MinPts：点的个数，也就是密度。一般取的小一些，多次尝试。

优势：

1. 不需要指定簇个数（基于密度，算法自己找）

2. 可以发现任意形状的簇 （如图，它可以发现K-means发现不出来的笑脸）
   

3. 擅长找到离群点（还有一些没有颜色的点就是离群点，不满足epsilon = 1.00，minPoints = 4，这是上图设置的两个参数。
   ）

4. 两个参数就够了

劣势：

1. 高维数据处理有些困难（理论是这样，但是博主的实力现在还做不了这种实验，但是可以从另一个角度想，高维数据，点之间稀疏，密度就很难定义）
2. 参数难以选择（参数对结果影响很大，如下图所示）
   当epsilon = 1.00，minPoints = 4时
   
   当epsilon = 0.92，minPoints = 2时
   

1.1 SciPy 聚类包

scipy.cluster是scipy下的一个做聚类的package, 共包含了两类聚类方法:

1. 矢量量化(scipy.cluster.vq):支持vector quantization 和 k-means 聚类方法
2. 层次聚类(scipy.cluster.hierarchy):支持hierarchical clustering 和 agglomerative clustering(凝聚聚类)

# coding=utf-8
from pylab import *
from scipy.cluster.vq import *

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)

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

title(u'2维数据点聚类', fontproperties=font)
axis('off')
show()

讲解：首先生成简单的二维数据，用 k=2 对这些数据进行聚类，由于 SciPy 中实现的 K-means 会计算若干次（默认为 20 次），并为我们选择方差最小的结果，所以这里返回的方差并不是我们真正需要的。用 SciPy 包中的矢量量化函数code,distance = vq(features,centroids)对每个数据点进行归类，为了将其可视化，可以画出这些数据点及最终的聚类中心，函数 where() 给出每个类的索引，绘制出的结果如上图所示。类中心标记为绿色大圆环，预测出的类分别标记为蓝色星号和红色点。

1.2 图像聚类

用 K-means 对这些字体图像进行聚类。

# -*- coding: utf-8 -*-
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

# Uses sparse pca codepath.
imlist = imtools.get_imlist('D:\\Python\\chapter6\\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('D:\\Python\\chapter6\\fontimages\\font_pca_modes.pkl', 'wb')
pickle.dump(immean, f)
pickle.dump(V, f)
f.close()

# 获取 selected-fontimages 文件下图像文件名，并保存在列表中 
imlist = imtools.get_imlist('D:\\Python\\chapter6\\selectedfontimages\\a_selected_thumbs\\')
imnbr = len(imlist)

# 载入模型文件 
with open('D:\\Python\\chapter6\\fontimages\\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()

字体图像聚成 4 类(k=4)。

1.1 在主成分上可视化图像

为了便于观察上面是如何利用主成分进行聚类的，我们可以在一对主成分方向的坐标上可视化这些图像。一种方法是将图像投影到两个主成分上，改变投影为：

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

以得到相应的坐标（在这里 V[[0,2]] 分别是第一个和第三个主成分）。当然，也可以将其投影到所有成分上，之后挑选出需要的列。

# -*- coding: utf-8 -*-
from PCV.tools import imtools, pca
from PIL import Image, ImageDraw
from pylab import *
from PCV.clustering import hcluster

imlist = imtools.get_imlist('D:\\Python\\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)])  
# projected = array([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 = 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:\\Python\\chapter6\\pca_font.png')
show()

在上图成对主成分上投影的字体图像中。左图用的是第一个和第二个主成分，右图用的是第二个和第三个主成分。这里，我们用到了整数或 floor 向下取整除法运算符 //，通过移去小数点后面的部 分，可以返回各个缩略图在白色背景中对应的整数坐标位置。可以很清楚地看到，二维投影后相似的字体图像距离较近。

其实学到这，博主（小白）主成分是什么都不是很懂，所以在这个补一下主成分的知识。

主成分分析（PCA）

用途：降维中最常用的一种手段
目标：提取最有价值的信息（基于方差,找最大的方差方向，因为越大的方差方向，就会使数据降维后越分的开，数据分的开，就能更好的进行分类任务。）
问题：降维后的数据的意义不知道（只有机器知道，因为在数据变换中经过一系列指针变换，物理意义会变换没了，但是并不影响最后结果，但是我们也是想要最后结果，而不是阶段性目标）

向量的表示及基变换
向量的表示
内积：${\left(a_1,a_2,...,a_n\right)}^T.{\left(b_1,b_2,...,b_n\right)}^T=a_1$