AI学习笔记（四）相机模型、图像聚类算法

相机模型

相机与图像

坐标系

针孔相机模型存在四个坐标系：世界坐标系、摄像机坐标系、图像物理坐标系和图像像素坐标系。
假设：
$\cdot$世界坐标系的坐标为$Pw(Xw,Yw,Zw)$，其为客观三维世界的绝对坐标系，也称客观坐标系。就是物体在真实世界中的坐标。世界坐标系是随着物体的大小和位置变化的，单位是长度单位；
$\cdot$对应的摄像机坐标系的坐标为$Po(x,y,z)$，其以相机的光心为坐标系的原点，以平行于图像的x和y方向为x轴和y轴，z轴和光轴平行，x，y，z互相垂直，单位是长度单位；
$\cdot$对应的图像物理坐标系的坐标为$P^{\prime }(x^{\prime },y^{\prime })$，其以主光轴和图像平面焦点为坐标原点，x’和y’方向如图所示，单位是长度单位；
$\cdot$对应的图像像素坐标系的坐标为$p(u,v)$，其以图像的顶点为坐标原点，u和v方向平行于x’和y’方向，单位是以像素计。

相机成像

世界坐标系到摄像机坐标系

这两个坐标系处理旋转矩阵R，还存在平移矩阵t。其关系可表示为：
$$\left[\begin{array}{c}{X}_C\\ Y_C\\ Z_C\\ 1\end{array}\right]=\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right]=L_W\left[\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right]$$
欧式变换
欧式变换由两部分组成：旋转和平移。${\alpha }^{\prime }=R\alpha +t$

齐次坐标
多次连续的旋转和平移的情况下。
假设我们将向量a进行了两次欧式变换、旋转和平移分别为$R_1,t_1$和$R_2,t_2$，分别得到：
$b=R_1\ast a+t_{1}c=R_2\ast b+t_2\Rightarrow c=R_2\ast \left(R_1\ast a+t_1\right)+t_2$
$$\left[\begin{array}{c}{\alpha }^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{cc}R& t\\ 0& 1\end{array}\right]\left[\begin{array}{c}\alpha \\ 1\end{array}\right]=T\left[\begin{array}{c}\alpha \\ 1\end{array}\right]$$$$b=T_1\tilde{\alpha },\tilde{c}=T_2\tilde{b}\to \tilde{c}=T_2{T}_1\tilde{\alpha }$$

摄像机坐标系到图像物理坐标系

根据相似三角形：

$$\{\begin{array}{l}{X}^{\prime }=f\frac{X_C}{Z_C}\\ Y^{\prime }=f\frac{Y_C}{Z_C}\end{array}矩阵形式为Z_C\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{ccc}f& 0& 0\\ 0& f& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{X}_C\\ Y_C\\ Z_C\end{array}\right]$$

图像物理坐标系到图像像素坐标系

$\{\begin{array}{l}u=\frac{x}{dx}+u_0\\ v=\frac{y}{dy}+v_0\end{array}$
$dx$和$dy$表示：$x$方向和$y$方向的一个像素分别占多少个（可能是小数）长度单位。
齐次坐标下：
$\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{dx}& 0& u_0\\ 0& \frac{1}{dy}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$

摄像机坐标系到图像像素坐标系

$Z_C\ast \left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{ccc}f& 0& 0\\ 0& f& 0\\ 0& 0& 1\end{array}\right]\ast \left[\begin{array}{c}{X}_C\\ Y_C\\ Z_C\end{array}\right]\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{dx}& 0& u_0\\ 0& \frac{1}{dy}& v_0\\ 0& 0& 1\end{array}\right]\ast \left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$
$Z_C\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}\frac{f}{dx}& 0& u_0\\ 0& \frac{f}{dy}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{X}_C\\ Y_C\\ Z_C\end{array}\right]=\left[\begin{array}{ccc}{f}_x& 0& u_0\\ 0& f_y& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{X}_C\\ Y_C\\ Z_C\end{array}\right]=K\left[\begin{array}{c}{X}_C\\ Y_C\\ Z_C\end{array}\right]$

世界坐标系到图像像素坐标系

$\left[\begin{array}{c}{X}_C\\ Y_C\\ Z_C\end{array}\right]=R\left[\begin{array}{c}{X}_W\\ Y_W\\ Z_W\end{array}\right]+t$
$Z_C\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=K\left(R{P}_W+t\right)=K\left(R\left[\begin{array}{c}{X}_W\\ Y_W\\ Z_W\end{array}\right]+t\right)$
$Z_C\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{dx}& 0& u_0\\ 0& \frac{1}{dy}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right]$

镜头畸变

畸变的起因于分类

• 镜头由于制造精度以及组装工艺的偏差会引入畸变，导致原始图像的失真。
• 镜头的畸变分为径向畸变和切向畸变两类。
  径向畸变（泰勒公式前几项）
  $$\{\begin{array}{l}{x}_{distorted}=x(1+k_1{r}^2+k_2{r}^4+k_3{r}^6)\\ y_{distorted}=y(1+k_1{r}^2+k_2{r}^4+k_3{r}^6)\end{array}$$切向畸变
  $$\{\begin{array}{l}{x}_{distorted}=x+\left[2p_{1}xy+p_2(r^2+2x^2)\right]\\ y_{distorted}=y+\left[p_1(r^2+2y^2)+2p_{2}xy\right]\end{array}$$
  

径向畸变

由透镜的形状引起的畸变称为径向畸变，透镜径向畸变后点位偏移示意图

切向畸变

切向畸变是由于透镜本身与相机传感器平面（成像平面）或图像平面不平行而产生的。这种情况多是由于透镜被粘贴到镜头模组上的安装偏差导致。
切向畸变示意图：

畸变矫正

• 径向畸变和切向畸变模型中一共有5个畸变参数，在opencv中他们被排列成一个人5*1的矩阵，依次包含$k_1,k_2,p_1,p_2,k_3$，经常被定义为Mat矩阵的形式，如Mat distCoeffs=Mat(1.5, CV_32FC1, Scalar::all(0));
• 求这5个参数就是相机标定中需要确定的5个畸变系数；
• 求得这5个参数后，就可以矫正由于镜头畸变引起的图像的变形失真。
  

透视变换

透视变换使将图片投影到一个新的平面(Viewing Plane)，也称作投影映射（Projective Mapping）。我们常说的仿射变换是透射变换的一个特例。透射变换的目的就是把现实中为直线的物体，在图片中可能呈现为斜线，通过透视变换转换成直线的变换。
仿射变换(Affine Transformation或Affine Map)，又称为仿射映射，是指在几何中，图像进行从一个向量空间进行一次现行变换和一次平移，变换到另一个向量空间的过程。

通用的变换公式为：
$$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$$
x、y是原始图片坐标，对应得到变换后的图片坐标(X’;Y’;Z’)其中：
$$\{\begin{array}{l}{X}^{\prime }=\frac{X}{Z}\\ Y^{\prime }=\frac{Y}{Z}\\ Z^{\prime }=\frac{Z}{Z}\end{array}\{\begin{array}{l}{X}^{\prime }=\frac{a_{11}x+a_{12}y+a_{13}}{a_{31}x+a_{32}x+a_{33}}\\ Y^{\prime }=\frac{a_{21}x+a_{22}y+a_{23}}{a_{31}x+a_{32}x+a_{33}}\\ Z^{\prime }=1\end{array}$$
当$a_{33}=1$，展开上面的公式，得到一个点的情况：
$$\{\begin{array}{l}{a}_{11}x+a_{12}y+a_{13}-a_{31}x{X}^{\prime }+a_{32}{X}^{\prime }y=X^{\prime }\\ a_{21}x+a_{22}y+a_{23}-a_{31}x{Y}^{\prime }+a_{32}y{Y}^{\prime }=Y^{\prime }\end{array}$$
源点四个坐标分别为$A:(x_0,y_0),(x_1,y_1),(x_2,y_2),(x_3,y_3)$
目标点四个坐标分别为$A:(X_0^{\prime },Y_0^{\prime }),(X_1^{\prime },Y_1^{\prime }),(X_2^{\prime },Y_2^{\prime }),(X_3^{\prime },Y_3^{\prime })$
$$\left[\begin{array}{cccccccc}{x}_0& y_0& 1& 0& 0& 0& -x_0{X}_0^{\prime }& -y_0{X}_0^{\prime }\\ 0& 0& 0& x_0& y_0& 1& -x_0{X}_0^{\prime }& -y_0{X}_0^{\prime }\\ x_1& y_1& 1& 0& 0& 0& -x_1{X}_1^{\prime }& -y_1{X}_1^{\prime }\\ 0& 0& 0& x_1& y_1& 1& -x_1{X}_1^{\prime }& -y_1{X}_1^{\prime }\\ x_2& y_2& 1& 0& 0& 0& -x_2{X}_2^{\prime }& -y_2{X}_2^{\prime }\\ 0& 0& 0& x_2& y_2& 1& -x_2{X}_2^{\prime }& -y_2{X}_2^{\prime }\\ x_3& y_3& 1& 0& 0& 0& -x_3{X}_3^{\prime }& -y_3{X}_3^{\prime }\\ 0& 0& 0& x_3& y_3& 1& -x_3{X}_3^{\prime }& -y_3{X}_3^{\prime }\end{array}\right]\left[\begin{array}{c}{a}_{11}\\ a_{12}\\ a_{13}\\ a21\\ a_{22}\\ a_{23}\\ a_{31}\\ a_{32}\end{array}\right]=\left[\begin{array}{c}{x}_0^{\prime }\\ y_0^{\prime }\\ x_1^{\prime }\\ y_1^{\prime }\\ x_2^{\prime }\\ y_2^{\prime }\\ x_3^{\prime }\\ y_3^{\prime }\end{array}\right]$$
透视变换代码如下：

import numpy as np

def WarpPerspectiveMatrix(src, dst):
    assert src.shape[0] == dst.shape[0] and src.shape[0] >= 4

    nums = src.shape[0]
    A = np.zeros((2 * nums, 8))  # A*warpMatrix=B
    B = np.zeros((2 * nums, 1))
    for i in range(0, nums):
        A_i = src[i, :]
        B_i = dst[i, :]
        A[2 * i, :] = [A_i[0], A_i[1], 1, 0, 0, 0,
                       -A_i[0] * B_i[0], -A_i[1] * B_i[0]]
        B[2 * i] = B_i[0]

        A[2 * i + 1, :] = [0, 0, 0, A_i[0], A_i[1], 1,
                           -A_i[0] * B_i[1], -A_i[1] * B_i[1]]
        B[2 * i + 1] = B_i[1]

    A = np.mat(A)
    # 用A.I求出A的逆矩阵，然后与B相乘，求出warpMatrix
    warpMatrix = A.I * B  # 求出a_11, a_12, a_13, a_21, a_22, a_23, a_31, a_32

    # 之后为结果的后处理
    warpMatrix = np.array(warpMatrix).T[0]
    warpMatrix = np.insert(
        warpMatrix,
        warpMatrix.shape[0],
        values=1.0,
        axis=0)  # 插入a_33 = 1
    warpMatrix = warpMatrix.reshape((3, 3))
    return warpMatrix


if __name__ == '__main__':
    print('warpMatrix')
    src = [[10.0, 457.0], [395.0, 291.0], [624.0, 291.0], [1000.0, 457.0]]
    src = np.array(src)

    dst = [[46.0, 920.0], [46.0, 100.0], [600.0, 100.0], [600.0, 920.0]]
    dst = np.array(dst)

    warpMatrix = WarpPerspectiveMatrix(src, dst)
    print(warpMatrix)

运行j结果如下：

warpMatrix
[[-5.01338334e-01 -1.35357643e+00  5.82386716e+02]
 [ 1.06858966e-15 -4.84035391e+00  1.38781980e+03]
 [ 4.33680869e-19 -4.14856327e-03  1.00000000e+00]]

opencv中的透视变换如下：

import cv2
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread('photo1.jpg')
result3 = img.copy()
# 注意这里src和dst的输入并不是图像，而是图像对应的顶点坐标。
src = np.float32([[207, 151], [517, 285], [17, 601], [343, 731]])
dst = np.float32([[0, 0], [337, 0], [0, 488], [337, 488]])
# 生成透视变换矩阵；进行透视变换
m = cv2.getPerspectiveTransform(src, dst)
result = cv2.warpPerspective(result3, m, (337, 488))
plt.subplot(121), plt.imshow(img), plt.title('src'), plt.axis("off")
plt.subplot(122), plt.imshow(result), plt.title('result'), plt.axis("off")
plt.show()

运行结果如下：

图像聚类算法

分类与聚类

分类

分类其实是从特定的数据中挖掘模式，做出判断的过程。
分类的主要过程：
(1). 训练数据集存在一个类标记号，判断它是正向数据集（起积极作用，不垃圾邮件），还是负向数据集（起抑制作用，垃圾邮件）；
(2).然后需要对数据集进行学习训练，并构建一个训练的模型；
(3).通过该模型对预测数据集进行预测，并计算其结果的性能。

聚类

从广义上来说，聚类就是将数据集中在某些方面相似的成员放在一起。一个聚类就是一些数据实例的集合，其中处于相同聚类中的数据元素彼此相似，但是处于不同聚类中的元素彼此不同。
由于在聚类中那些表示数据类别的分类或分组信息是没有的，即这些数据是没有标签的，所以聚类通常被称为无监督学习（Unsupervised Learning）。
聚类的目的也是数据分类，但是事先是不知道如何去分的，完全是算法子集来判断各条数据之间的相似性，相似的就放到一起。
在聚类的结论出来之前，完全不知道每一类有什么特点，一定要根据聚类的结果通过人的经验来分析，看看聚成的这一类大概有什么特点。
总之，聚类主要是“物以类聚”，通过相似性把相似元素聚集在一起，它没有标签；而分类通过标签来训练得到一个模型，对新数据进行预测的过程，其数据存在标签。
聚类样本间的属性包括有序属性和无序属性。

聚类常见的算法

• 原型聚类（包含3个子类算法：K均值聚类算法、学习向量化、高斯混合聚类）
• 层次聚类
• 密度聚类

K-Means聚类

K-Means聚类是最常用的聚类算法，最初起源与信号处理，其目标是将数据点划分为K个类簇，找到每个簇的中心并使其度量最小化。
该算法最大的优点在于简单，便于理解，运算速度较快，缺点是只能应用于连续型数据，并且要在聚类前指定聚集的类簇数，是一种原型聚类算法。

K-Means聚类算法的分析流程

第一步：确定K值，即数据集聚集成K个类簇或小组；
第二步：从数据集中速记选择K个数据点作为质心（Centroid）或数据中心；
第三步：分别计算每个点到每个质心之间的距离，并将每个点划分到离最近质心的小组；
第四步：当每个质心都聚集了一些点后，重新定义算法选出新的质心（对于每个簇，计算其均值，记得到新的K个质心点）；
第五步：迭代执行第三步到第四步，直到迭代终止条件满足位置（分类结果不再变化）

K-Means聚类算法举例

第一步：选取K=2为例；
第二步：假设选取P1、P2作为初始的质心
第三步：分别计算每个点到质心的距离：
经比较可得，第一次分组结果是：

• 组A：P1
• 组B：P2、P3、P4、P5、P6
  第四步：当每个质心都聚集了一些点后，重新定义算法选出新的质心；组A只有P1点，则Pnew1_1=P1，组B为5个点的均值有Pnew1_2=(6.2, 5.6)，P2-P6称为新的离散点；
  再次计算到质心的距离：
  
  经比较可得，第二次分组结果是：
• 组A：P1、P2、P3
• 组B：P4、P5、P6（虚拟质心这时候消失）
  按照上一次的方法选出两个新的虚拟质心：
  ----Pnew2_1(1.33,1) ，Pnew2_2(9,8.33)
  再次计算到质心的距离：
  
  经比较可得，第三次分组结果是：
• 组A：P1、P2、P3
• 组B：P4、P5、P6
  此时，这次分组的结果和上次没有任何变化了，说明已经收敛，聚类结束。
  以运动员每分钟助攻数与得分数为例，在用numpy实现K-Means聚类如下：

# coding=utf-8
import matplotlib.pyplot as plt
import numpy as np
from sklearn.cluster import KMeans

"""
第一部分：数据集
X表示二维矩阵数据，篮球运动员比赛数据
总共20行，每行两列数据
第一列表示球员每分钟助攻数：assists_per_minute
第二列表示球员每分钟得分数：points_per_minute
"""
X = [[0.0888, 0.5885],
     [0.1399, 0.8291],
     [0.0747, 0.4974],
     [0.0983, 0.5772],
     [0.1276, 0.5703],
     [0.1671, 0.5835],
     [0.1306, 0.5276],
     [0.1061, 0.5523],
     [0.2446, 0.4007],
     [0.1670, 0.4770],
     [0.2485, 0.4313],
     [0.1227, 0.4909],
     [0.1240, 0.5668],
     [0.1461, 0.5113],
     [0.2315, 0.3788],
     [0.0494, 0.5590],
     [0.1107, 0.4799],
     [0.1121, 0.5735],
     [0.1007, 0.6318],
     [0.2567, 0.4326],
     [0.1956, 0.4280]
     ]

"""
第二部分：KMeans聚类
clf = KMeans(n_clusters=3) 表示类簇数为3，聚成3类数据，clf即赋值为KMeans
y_pred = clf.fit_predict(X) 载入数据集X，并且将聚类的结果赋值给y_pred
"""
clf = KMeans(n_clusters=3)
y_pred = clf.fit_predict(X)
# 输出完整Kmeans函数，包括很多省略参数
print(clf)
# 输出聚类预测结果
print("y_pred = ", y_pred)

"""
第三部分：可视化绘图
"""
# 获取数据集的第一列和第二列数据 使用for循环获取 n[0]表示X第一列
x = [n[0] for n in X]
# print(x)
y = [n[1] for n in X]
# print(y)

'''
绘制散点图
参数：x横轴; y纵轴; c=y_pred聚类预测结果; marker类型:o表示圆点,*表示星型,x表示点;
'''
plt.scatter(x, y, c=y_pred, marker='x')
plt.title("Kmeans-Basketball Data")  # 绘制标题
plt.xlabel("assists_per_minute")  # 绘制x轴和y轴坐标
plt.ylabel("points_per_minute")
plt.show()  # 显示图形

运行结果如下：

KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
       n_clusters=3, n_init=10, n_jobs=None, precompute_distances='auto',
       random_state=None, tol=0.0001, verbose=0)
y_pred =  [1 2 1 1 1 1 1 1 0 0 0 1 1 1 0 1 1 1 1 0 0]

K-Means聚类与图像处理

在图像中，通过K-Means聚类算法可以实现图像分割、图像聚类、图像识别等操作。我们通过K-Means可以将这些像素点聚类成K个簇，然后使用每个簇内的质心来替换簇内所有的像素点，这样就能实现在不改变分辨率情况下量化压缩图像颜色，实现图像颜色层级分割。
用opencv中用Kmeans函数对图像进行聚类代码如下：
例子1：

# coding: utf-8

"""
在OpenCV中，Kmeans()函数原型如下所示：
retval, bestLabels, centers = kmeans(data, K, bestLabels, criteria, attempts, flags[, centers])
    data表示聚类数据，最好是np.flloat32类型的N维点集
    K表示聚类类簇数
    bestLabels表示输出的整数数组，用于存储每个样本的聚类标签索引
    criteria表示算法终止条件，即最大迭代次数或所需精度。在某些迭代中，一旦每个簇中心的移动小于criteria.epsilon，算法就会停止
    attempts表示重复试验kmeans算法的次数，算法返回产生最佳紧凑性的标签
    flags表示初始中心的选择，两种方法是cv2.KMEANS_PP_CENTERS ;和cv2.KMEANS_RANDOM_CENTERS
    centers表示集群中心的输出矩阵，每个集群中心为一行数据
"""

import cv2
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread('lenna.png', 0)  # 读取原始图像灰度颜色
print(img.shape)
rows, cols = img.shape[:]  # 获取图像高度、宽度
data = img.reshape((rows * cols, 1))  # 图像二维像素转换为一维
data = np.float32(data)
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)  # 定义中心 (type,max_iter,epsilon)
flags = cv2.KMEANS_RANDOM_CENTERS  # 设置标签
compactness, labels, centers = cv2.kmeans(data, 4, None, criteria, 10, flags)  # K-Means聚类 聚集成4类
dst = labels.reshape((img.shape[0], img.shape[1]))  # 生成最终图像
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签

# 显示图像
titles = [u'原始图像', u'聚类图像']
images = [img, dst]
for i in range(2):
    plt.subplot(1, 2, i + 1), plt.imshow(images[i], 'gray'),
    plt.title(titles[i])
    plt.xticks([]), plt.yticks([])
plt.show()

运行结果如下：

例子2：

# coding: utf-8
import cv2
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread('lenna.png')
print(img.shape)

data = img.reshape((-1, 3))  # 图像二维像素转换为一维
data = np.float32(data)
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)  # 定义中心 (type,max_iter,epsilon)
flags = cv2.KMEANS_RANDOM_CENTERS  # 设置标签
compactness, labels2, centers2 = cv2.kmeans(data, 2, None, criteria, 10, flags)  # K-Means聚类 聚集成2类
compactness, labels4, centers4 = cv2.kmeans(data, 4, None, criteria, 10, flags)  # K-Means聚类 聚集成4类
compactness, labels8, centers8 = cv2.kmeans(data, 8, None, criteria, 10, flags)  # K-Means聚类 聚集成8类
compactness, labels16, centers16 = cv2.kmeans(data, 16, None, criteria, 10, flags)  # K-Means聚类 聚集成16类
compactness, labels64, centers64 = cv2.kmeans(data, 64, None, criteria, 10, flags)  # K-Means聚类 聚集成64类

# 图像转换回uint8二维类型
centers2 = np.uint8(centers2)
res = centers2[labels2.flatten()]
dst2 = res.reshape((img.shape))

centers4 = np.uint8(centers4)
res = centers4[labels4.flatten()]
dst4 = res.reshape((img.shape))

centers8 = np.uint8(centers8)
res = centers8[labels8.flatten()]
dst8 = res.reshape((img.shape))

centers16 = np.uint8(centers16)
res = centers16[labels16.flatten()]
dst16 = res.reshape((img.shape))

centers64 = np.uint8(centers64)
res = centers64[labels64.flatten()]
dst64 = res.reshape((img.shape))

# 图像转换为RGB显示
img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
dst2 = cv2.cvtColor(dst2, cv2.COLOR_BGR2RGB)
dst4 = cv2.cvtColor(dst4, cv2.COLOR_BGR2RGB)
dst8 = cv2.cvtColor(dst8, cv2.COLOR_BGR2RGB)
dst16 = cv2.cvtColor(dst16, cv2.COLOR_BGR2RGB)
dst64 = cv2.cvtColor(dst64, cv2.COLOR_BGR2RGB)

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签

# 显示图像
titles = [u'原始图像', u'聚类图像 K=2', u'聚类图像 K=4',
          u'聚类图像 K=8', u'聚类图像 K=16', u'聚类图像 K=64']
images = [img, dst2, dst4, dst8, dst16, dst64]
for i in range(6):
    plt.subplot(2, 3, i + 1), plt.imshow(images[i], 'gray'),
    plt.title(titles[i])
    plt.xticks([]), plt.yticks([])
plt.show()

运行结果如下：

K-Means聚类算法优缺点

优点
1、是解决聚类问题的一种经典算法：简单、快速
2、对处理大数据集，该算法保持可伸缩性和高效率；
3、当结果簇是密集的，它的效果较好
缺点
1、在簇的平均值可被定义的情况下才能被使用，可能不适用于某些应用；
2、必须事先给出k(要生成簇的数目)，而且对初值敏感，对于不同的初始值，可能导致不同结果。
3、不适合于法线非凸形状的簇或者大小差别很大的簇
4、对噪声和孤立点数据敏感

层次聚类

层次聚类是一种很直观的算法，顾名思义就是要一层一层地进行聚类。
层次法(Hierarchical methods) 先计算样本之间的距离。每次距离最近的点合并到同一个类。然后，在计算类与类之间的距离，将距离最近的类合并为一个大类。不停地合并，知道合成一个类。其中类与类距离的计算方法有：最短距离法，最长距离法，中间距离法，类平均法等。比如最短距离法，将类与类的距离定义为类与类之间样本的最短距离。
层次聚类算法根据层次分解的顺序分为：自下而上和自上而下，即凝聚的层次聚类算法和分裂的层次聚类算法(agglomerative和divise)，也可以理解为自下而上法(bottom-up)和自上而下法(top-down)。

层次聚类算法的流程

凝聚型层次聚类的策略是先j将每个对象作为一个簇，然后合并这些原子粗为越来越大的簇，知道所有对象都在一个簇中，或者某个中间条件满足。绝大多数层次聚类属于凝聚型层次剧烈，它们只是在簇间相似度定义上有所不同。这里给出采用最小距离的邻居层次聚类算法流程：
1、将每个对象看做一类，计算两两之间的最小距离；
2、将距离最小的两个类合并成一个新类；
3、重新计算新类与所有类之间的距离；
4、重复2、3知道所有类最后合并成一类。

层次聚类算法的特点

特点

• 凝聚的层次聚类并没有类似K均值的全局目标函数，没有局部极小问题或是很难选择初始点的问题；
• 合并的操作往往是最终的，一旦合并两个簇之后就不会撤销
• 计算的代价也比较昂贵。
  优点
• 1、距离和规则相似度容易定义，限制少；
• 2、不需要预先限定聚类数；
• 3、可以发现类的层次关系；
• 4、可以聚类成其他形状。
  缺点
• 1、j计算复杂度太高；
• 2、奇异值也能产生较大影响；
• 3、算法很可能聚类成链状。
  
  
  Z的第一行：[0,3]意思是类别0和类别3距离最近，首先聚成一类，并且自动定义类别为5=(len(X)-1+1)
  
  Z的第二行：[4,5]意思是类别4和类别5距离最近，4、5聚成一类，类别为6=(len(X)-1+2)
  
  第三行，第四行以此类推：
  因为类别5有两个样本，加上类别4形成类别6，有3个样本；
  类别7是类别1、2聚类形成，有两个样本；
  类别6、7聚成一类后，类别8有5个样本，这样X全部样本参与聚类，聚类完成。
  Z第4列有样本的个数，当最下面一行中的样本数达到样本总数时，聚类就完成了。
  树状图分类判断
  想分两类时，就从上往下数有两个竖线是进行切割，那么所对应的竖线下面所连接的为一类
  
  想分三类时，就从上往下数有三根竖线时进行切割，那么所对应的竖线下面所连接的为一类
  
  K_Means与层次聚类，每一种聚类方法都有特定的数据结构，对于服从高斯分布的数据用K-Means来进行聚类效果比较好。而对于类别之间存在层结构的数据，用层次聚类会比较好。
  层次聚类代码实现如下：

# 导入相应的包
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
from matplotlib import pyplot as plt
'''
linkage(y, method=’single’, metric=’euclidean’) 共包含3个参数:
1. y是距离矩阵,可以是1维压缩向量（距离向量），也可以是2维观测向量（坐标矩阵）。
若y是1维压缩向量，则y必须是n个初始观测值的组合，n是坐标矩阵中成对的观测值。
2. method是指计算类间距离的方法。
'''
'''
fcluster(Z, t, criterion=’inconsistent’, depth=2, R=None, monocrit=None)
1.第一个参数Z是linkage得到的矩阵,记录了层次聚类的层次信息;
2.t是一个聚类的阈值-“The threshold to apply when forming flat clusters”。
'''
X = [[1, 2], [3, 2], [4, 4], [1, 2], [1, 3]]
Z = linkage(X, 'ward')
f = fcluster(Z, 4, 'distance')
fig = plt.figure(figsize=(5, 3))
dn = dendrogram(Z)
print(Z)
plt.show()

运行结果如下：

密度聚类DBSCAN

具体算法

需要两个参数$\epsilon$(eps)和形成高密度区域所需要的的最少点数(minPts)

• 它由一个任意未被访问的点开始，然后探索这个点的$\epsilon$-邻域里有足够的点，则建立一个新的聚类，否则这个点被标注为杂音。
• 足以，这个点之后可能被发现了在其它点的$\epsilon$-邻域里，而该$\epsilon$-邻域可能有足够的点，届时这个点会被加入到该聚类中。
  优点
  1、对噪声不敏感
  2、能发现任意形状的聚类
  缺点
  1、聚类的结果与参数有很大的关系；
  2、用固定参数识别聚类，但当聚类的稀疏成都不同是，相同的判定准则可能会破坏聚类的自然结构，即较稀的聚类会被划分为多个类和密度较大且离得较近的类会被合并成一个聚类。
  密度聚类代码实现如下：

mport matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets
from sklearn.cluster import DBSCAN

iris = datasets.load_iris()
X = iris.data[:, :4]  # #表示我们只取特征空间中的4个维度
print(X.shape)
# 绘制数据分布图
'''
plt.scatter(X[:, 0], X[:, 1], c="red", marker='o', label='see')
plt.xlabel('sepal length')
plt.ylabel('sepal width')
plt.legend(loc=2)
plt.show()
'''

dbscan = DBSCAN(eps=0.4, min_samples=9)
dbscan.fit(X)
label_pred = dbscan.labels_

# 绘制结果
x0 = X[label_pred == 0]
x1 = X[label_pred == 1]
x2 = X[label_pred == 2]
plt.scatter(x0[:, 0], x0[:, 1], c="red", marker='o', label='label0')
plt.scatter(x1[:, 0], x1[:, 1], c="green", marker='*', label='label1')
plt.scatter(x2[:, 0], x2[:, 1], c="blue", marker='+', label='label2')
plt.xlabel('sepal length')
plt.ylabel('sepal width')
plt.legend(loc=2)
plt.show()

运行结果如下：

谱聚类

步骤
1、根据数据构造一个图结构(Graph)，Graph的每一个节点对应一个数据点，将相似的点连接起来，并且边的权重用于表示数据之间的相似度。把这个Graph用邻接矩阵的形式表示起来，记为W.
2、把W的每一类元素加起来得到N个数，把他们放在对角线上（其他地方都是零），组成一个N*N的矩阵，记为D，并零L=D-W.
3、求出L的前K个特征值，以及对应的特征向量。
4、把这K个特征（列）向量排列在一起组成一个$N\ast k$的矩阵，将其中每一行看做k维空间中的一个向量，并使用K-means算法进行聚类。聚类的结果中每一行所属的类别就是原来Graph中节点亦即最初的N个数据点分别所属的类别。
简单抽闲谱聚类过程，主要分为两步：
1、构图，将采样点数据构造成一张网图。
2、切图，即将第一步构造出来的按照一定的切边准则，切分成不同的图，而不同的子图，即我们对应的聚类结果。

谱聚类代码实现如下：

# encoding=utf-8
import matplotlib.pyplot as plt
import numpy as np
from numpy import linalg as LA
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.metrics.pairwise import rbf_kernel
from sklearn.preprocessing import normalize
 
 
def similarity_function(points):
    """
    相似性函数，利用径向基核函数计算相似性矩阵，对角线元素置为０
    对角线元素为什么要置为０我也不清楚，但是论文里是这么说的
    :param points:
    :return:
    """
    res = rbf_kernel(points)
    for i in range(len(res)):
        res[i, i] = 0
    return res
 
 
def spectral_clustering(points, k):
    """
    谱聚类
    :param points: 样本点
    :param k: 聚类个数
    :return: 聚类结果
    """
    W = similarity_function(points)
    # 度矩阵D可以从相似度矩阵W得到，这里计算的是D^(-1/2)
    # D = np.diag(np.sum(W, axis=1))
    # Dn = np.sqrt(LA.inv(D))
    # 本来应该像上面那样写，我做了点数学变换，写成了下面一行
    Dn = np.diag(np.power(np.sum(W, axis=1), -0.5))
    # 拉普拉斯矩阵：L=Dn*(D-W)*Dn=I-Dn*W*Dn
    # 也是做了数学变换的，简写为下面一行
    L = np.eye(len(points)) - np.dot(np.dot(Dn, W), Dn)
    eigvals, eigvecs = LA.eig(L)
    # 前k小的特征值对应的索引，argsort函数
    indices = np.argsort(eigvals)[:k]
    # 取出前k小的特征值对应的特征向量，并进行正则化
    k_smallest_eigenvectors = normalize(eigvecs[:, indices])
    # 利用KMeans进行聚类
    return KMeans(n_clusters=k).fit_predict(k_smallest_eigenvectors)
 
'''
sklearn中的make_blobs()函数---为聚类产生数据集
  n_samples是待生成的样本的总数。
  n_features是每个样本的特征数。
  centers表示类别数。
  cluster_std表示每个类别的方差，例如我们希望生成2类数据，其中一类比另一类具有更大的方差，可以将cluster_std设置为[1.0,3.0]。
'''

X, y = make_blobs(n_samples=100, n_features=2, centers=2)
labels = spectral_clustering(X, 3)
 
# 画图
plt.style.use('ggplot')
# 原数据
fig, (ax0, ax1) = plt.subplots(ncols=2)
ax0.scatter(X[:, 0], X[:, 1], c=y)
ax0.set_title('raw data')
# 谱聚类结果
ax1.scatter(X[:, 0], X[:, 1], c=labels)
ax1.set_title('Spectral Clustering')
plt.show()

运行结果如下：