OpenCV 例程300篇234. 特征提取之主成分分析（PCA）

Posted 2022-09-01 YouCans

『youcans 的 OpenCV 例程300篇 - 总目录』

【youcans 的 OpenCV 例程300篇】234. 特征提取之主成分分析（PCA）

5.1 特征提取的方法

初步获取的图像特征维数通常很大，而且往往包含一定的无关或冗余特征。特征提取是指从原始特征中通过数学变换得到一组新的特征，以降低特征维数，消除相关性，减少无用信息。

降维方法可以从事物之间错综复杂的关系中找出一些主要因素，从而能有效利用大量统计数据进行定量分析，解释变量之间的内在关系，得到对事物特征及其发展规律的一些深层次的启发。

众多原始变量之间往往具有一定的相关关系。这意味着相关变量所反映的信息有一定程度的重叠，因此可以用较少的综合指标聚合、反映众多原始变量所包含的全部信息或主要信息。

特征提取分为线性映射方法和非线性映射方法。

线性映射方法主要有：

• 主成分分析（PCA）：按均方误差损失最小化原则，将高维原始数据空间变换到低维特征向量空间。
• 线性判别函数（LDA）：将高维原始数据向线性判别超平面的法向量上投影，使样本数据在新投影空间具有最大的类间距离和最小的类内距离，实现最大区分度（高内聚，低耦合）。
• 多维标度分析 （MDS）：根据样本之间的距离关系或不相似度关系在低维空间里生成对样本的一种表示。度量型 MDS 把样本间的距离关系或不相似度关系看作一种定量的度量，尽可能的在低维空间里保持这种度量关系；非度量型 MDS 把样本间的距离关系或不相似度关系看作一种定性的关系，在低维空间里只需保持这种关系的顺序。

非线性映射方法主要有：

• 基于核的非线性降维：通过核函数对样本进行非线性变换后，再在变换空间进行线性映射，如核主成分分析（KPCA）、核线性判别函数（KLDA）。核方法通过选择不同的核函数类型，反映了对数据分布的不同假设，可以看作对数据引入了一种非线性距离度量。
• 二维化和张量化： 将数据映射到二维空间，如二维主成分分析（2DPCA）、二维线性判别分析（2DLDA）、二维典型相关分析（2DCCA）；
• 流形学习方法：从高维采样数据中恢复低维流形结构并求出相应的嵌入映射。基本思想是通过局部距离来定义非线性距离度量，在样本分布较密集的情况下可以实现各种复杂的非线性距离度量。 具体方法如等距特征映射 （ISOMap） ， 拉普拉斯特征映射 （LE）， 局部线性嵌入 （LPP）。

此外，还可以通过聚类分析、神经网络方法进行数据降维。本质上，非线性映射的思想和算法与神经网络是相通的。

5.2 主成分分析的数学方法

主成分分析（Principal Components Analysis，PCA）是一种基于统计的数据降维方法，又称主元素分析、主分量分析。主成分分析只需要特征值分解，就可以对数据进行压缩、去噪，应用非常广泛。

众多原始变量之间往往具有一定的相关关系。这意味着相关变量所反映的信息有一定程度的重叠，因此可以用较少的综合指标聚合、反映众多原始变量所包含的全部信息或主要信息。主成分分析方法研究特征变量之间的相关性、相似性，将一组相关性高的高维变量转换为一组彼此独立、互不相关的低维变量，从而降低数据的维数。

主成分分析方法的思想是，将高维特征（p维）映射到低维空间（k维）上，新的低维特征是在原有的高维特征基础上通过线性组合而重构的，并具有相互正交的特性，称为主成分特性。

通过正交变换构造彼此正交的新的特征向量，这些特征向量组成了新的特征空间。将特征向量按特征值排序后，样本数据集中所包含的全部方差，大部分就包含在前几个特征向量中，其后的特征向量所含的方差很小。因此，可以只保留前 k个特征向量，而忽略其它的特征向量，实现对数据特征的降维处理。

主成分分析的基本步骤是：对原始数据归一化处理后求协方差矩阵，再对协方差矩阵求特征向量和特征值；对特征向量按特征值大小排序后，依次选取特征向量，直到选择的特征向量的方差占比满足要求为止。

主成分分析方法得到的主成分变量具有几个特点：（1）每个主成分变量都是原始变量的线性组合；（2）主成分的数目大大少于原始变量的数目；（3）主成分保留了原始变量的绝大多数信息；（4）各主成分变量之间彼此相互独立。

算法的基本流程如下：

（1）归一化处理，数据减去平均值；
（2）通过特征值分解，计算协方差矩阵；
（3）计算协方差矩阵的特征值和特征向量；
（4）将特征值从大到小排序；
（5）依次选取特征值最大的 k个特征向量作为主成分，直到其累计方差贡献率达到要求；
（6）将原始数据映射到选取的主成分空间，得到降维后的数据。

主成分分析方法的主要优点是：
（1）仅以方差衡量信息量，不受数据集以外的因素影响；　
（2）各主成分之间正交，可消除原始数据各变量之间的相互影响；
（3）方法简单，易于实现。

在图像处理中，把每幅二维图像拉伸为一维向量，即展平为一维数组。一组 m 幅图像就构造为一个 m 维向量，使用 Karhunen-Loève transform（KLT） 变换得到变换矩阵，选取特征值最大的 k个特征向量作为主成分，从而实现特征降维。

图像压缩过程是把一组原始图像变换成低维向量的过程，图像重建就是由低维向量变换重建图像组的过程。使用主成分分析进行图像压缩和重建会有少量信息损失，但可以把损失控制到很小。

对于一组 P 维向量 X，通过线性变换转换为另一组 K 维向量 Y。

向量 Xi 包含 m 个数据样本，记为：

$$\begin{array}{rl}X& =[X_1,X_2,...,X_p]\\ & =\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1p}\\ x_{21}& x_{22}& \cdots & x_{2p}\\ ⋮& ⋮& \ddots & ⋮\\ x_{m1}& x_{m2}& \cdots & x_{mp}\end{array}\right]\end{array}$$

计算平均值
$$u_j=E{X}_j=\frac{1}{m}\sum _{i=1}^m{X}_{ij},j=1,...p$$

计算协方差矩阵：
$$C_x=E(X-u)(X-u{)}^T$$

协方差矩阵 $C_x$ 是实对称矩阵，可以协方差矩阵的特征向量和对应的特征值。

通过变换矩阵 A 将 X 映射为由 Y 表示的向量：
$$Y=A\ast (X-u)$$
Y 的协方差矩阵用对角阵 A 和 Cx 表示，即：
$$C_y=A\ast C_x\ast A^T$$
则 Cy 是一个对角阵：
$$Cy=\left[\begin{array}{cccc}{\lambda }_1& & & 0\\ & {\lambda }_2& & \\ & & \ddots & \\ 0& & & {\lambda }_p\end{array}\right]$$
${\lambda }_i$ 是 $Cx$ 的特征值，将其按降序排列。

选择前 K 个特征向量，作为主成分特征，可以重建 p 维向量 X，并使均方误差最小。
$$\hat{X}=A_k^T\ast T+u$$

例程 14.15：特征描述之主分量

本例程的图像来自 R.C.Gonzalez 《数字图像处理（第四版）》P622 例11.16。

原始图像显示了对应于 6 个波段的 6 幅多光谱卫星图像：可见蓝光（450-520nm）、可见绿光（520-600nm）、可见红光（630-690nm）、近红外（760-900nm）、中红外（1550-1750nm）和热红外（10400-12500nm）。

本例的目的是说明如何使用主分量作为图像特征。

python
    # # 14.15 特征描述之主成分分析
    # 读取光谱图像组
    img = cv2.imread("../images/Fig1138a.tif", flags=0)
    height, width = img.shape[:2]  # (564, 564)
    nBands = 6  # 光谱波段种类
    snBands = ['a','b','c','d','e','f']  # Fig1138a~f
    imgMulti = np.zeros((height, width, nBands))  # (564, 564, 6)
    mbMatrix = np.zeros((img.size, nBands))  # (318096, 6)
    print(imgMulti.shape, mbMatrix.shape)

    # 显示光谱图像组
    fig1 = plt.figure(figsize=(9, 6))  # 原始图像，6 个不同波段
    fig1.suptitle("Spectral image of multi bands by NASA")
    for i in range(nBands):
        path = "../images/Fig1138.tif".format(snBands[i])
        imgMulti[:,:,i] = cv2.imread(path, flags=0)  # 灰度图像
        ax1 = fig1.add_subplot(2,3,i+1)
        ax1.set_xticks([]), ax1.set_yticks([])
        ax1.imshow(imgMulti[:,:,i], 'gray')  # 绘制光谱图像 snBands[i]
    plt.tight_layout()

    # 主成分分析 (principal component analysis)
    mbMean = np.zeros((nBands,))
    for i in range(nBands):
        mbArray = imgMulti[:,:,i].flatten()  # 转为一维数组
        mbMean[i] = mbArray.mean()
        mbMatrix[:,i] = mbArray - mbMean[i]  # 数据标准化 (318096, 6)

    # cov = np.cov(mbMatrix.transpose(),rowvar=0)  # 以列为变量计算协方差矩阵 (6,6)
    cov = np.cov(mbMatrix, rowvar=0)  # 以列为变量计算协方差矩阵 (6,6)
    # eigenvalues: 特征值，一维数组 (M,)，eigenvectors: 特征向量，二维矩阵 (M,M)
    eigenValues, eigenVectors = np.linalg.eig(cov)  # 计算特征值和特征向量
    principal = np.matmul(mbMatrix, eigenVectors)  # 主元素变换，投影到特征向量方向 (318096, 6)
    print(eigenValues.shape, eigenVectors.shape, principal.shape)
    print("Eigenvalues:\\n", eigenValues.round(4))  # 特征值，从大到小
    # print("Eigenvectors:\\n", eigenVectors.round(4))

    # 显示主成分变换图像
    fig2 = plt.figure(figsize=(9, 6))  # 主元素图像，q=6
    fig2.suptitle("Principal component images")
    imgPCA = np.zeros((height, width, nBands))  # (564, 564, 6)
    for i in range(nBands):
        pca = principal[:, i].reshape(-1, img.shape[1])  # 主元素图像 (564, 564)
        imgPCA[:,:,i] = np.uint8(cv2.normalize(-pca, (height, width), 0, 255,  cv2.NORM_MINMAX))
        ax2 = fig2.add_subplot(2,3,i+1)
        ax2.set_xticks([]), ax2.set_yticks([])
        ax2.imshow(imgPCA[:,:,i], 'gray')  # 绘制主成分图像
    plt.tight_layout()

    # 保留的主成分数量
    remainRatio = 0.95  # 设定的主成分累计方差贡献率
    remainValue = remainRatio * sum(eigenValues)
    IndexDes = np.argsort(-eigenValues)  # 特征值降序排列的索引
    topKValue = 0.0
    for i in range(len(eigenValues)):
        topKValue += eigenValues[IndexDes[i]]  # 降序排列的第 i 个特征值
        print("k=, topKValue = :.2f, topKRatio = :.4f"
              .format(i, topKValue, topKValue/sum(eigenValues)))
        if topKValue > remainValue:
            K = i + 1  # （0，1，2）-> K=3
            break
    print("number of PCA features: K=", K)  # 主成分方差贡献率 95% 时的特征维数 K=3
    indexPCA = IndexDes[:K]  # 选择特征值最大的 K 个特征向量的索引
    eigenVectPCA = eigenVectors[:, indexPCA]  # 选择 K 个主要特征向量组成降维特征矩阵 (P=6, K=3)
    print("PCA eigenvalues:\\n", eigenValues[indexPCA].round(4))  # 特征值，从大到小
    print("PCA eigenvectors:\\n", eigenVectPCA.round(4))

运行结果：

(564, 564, 6) (318096, 6)
(6,) (6, 6) (318096, 6)
Eigenvalues:
[10344.3048 2965.8977 1400.635 203.4555 94.2774 31.0373]
k=0, topKValue = 10344.30, topKRatio = 0.6878
k=1, topKValue = 13310.20, topKRatio = 0.8850
k=2, topKValue = 14710.84, topKRatio = 0.9781
number of PCA features: K= 3
PCA eigenvalues:
[10344.3048 2965.8977 1400.635 ]
PCA eigenvectors:
[[-0.489 -0.0124 -0.2301]
[-0.4777 0.0394 -0.3012]
[-0.4899 -0.022 -0.315 ]
[ 0.1375 0.7986 0.0431]
[-0.2188 0.5981 0.0165]
[-0.4753 -0.0486 0.8689]]

【本节完】

版权声明：
本例程的图像来自 R.C.Gonzalez 《数字图像处理（第四版）》P622 例11.16。
youcans@xupt 原创作品，转载必须标注原文链接：(https://blog.csdn.net/youcans/article/details/125761655)
Copyright 2022 youcans, XUPT
Crated：2022-7-15

234. 特征提取之主成分分析（PCA）
235. 特征提取之主成分分析（sklearn）
2