数字图像处理学习笔记 十一 表示与描述

前言综述：

将一幅图像分割成多个区域后，分割后的像素集经常以一种合适于计算机进一步处理的形式来表示和描述。

表示
表示一个区域的两种选择：

1. 根据区域外部特征（如边界）来表示区域
2. 根据内部特征（如组成该区域的像素）来表示区域

表示是直接具体地表示目标，好的表示方法应具有节省存储空间易于特征计算等优点。
表示的下一步工作是描述。表示方式的选择要使数据有利于描述工作的展开。
当关注的重点是形状特征时，可选择外部表示；当关注的重点是内部属性（如颜色和纹理）时，可选择内部表示。

描述
描述是较抽象地表示目标。好的描绘子都应尽可能对目标的大小、平移、旋转等不敏感，这样的描绘子比较通用。

表示和描述应该是个递进的关系，表示旨在以更（精确/方便/高效）的方式组织数据，而描述旨在从表示中总结某种模式以便于任务的完成。表示侧重于数据结构，而描述侧重于区域特性以及不同区域间的联系和差别。
常见的图像特征分为灰度、纹理和几何形状等。其中，灰度和纹理属于内部特征，几何形状属于外部特征。

（一）表示

分割技术会获得原始数据，其形式是沿着边界或包含在区域中的像素。标准做法是使用某种方案将分割后的数据精简为便于描绘子计算的表示。
以下介绍的算法要求一个区域的边界上的点以顺时针（或逆时针）方向排序。

1.1 边界追踪算法

该算法输出的是排序后的点序列，我们假设：

• 处理的是二值图像，其目标和背景点分别标位1和0

• 图像已使用值为0的边界填充，因而消除了目标与图像边界合并的可能性

给定一个二值区域R或其边界，追踪R的边界或给定边界的算法由如下步骤组成：

• 令起始点$b_0$​为图像中左上角标记为1的点。使用$c_0$表示$b_0$西侧的邻点。从$c_0$​开始按顺时针方向考察$b_0$​的8个邻点。令$b_1$表示所遇到的值为1的第一个邻点，并直接令$c_0$​是序列中$b_1$​之前的点。存储$b_0$和$b_1$​的位置，以便在步骤5中使用
• 令$b=b_1$​和$c=c_1$​
• 从c开始按顺时针方向行进，令b的8个邻点为n1,n2,…,n8 n1​,n2​,…,n8​找到标记为1的第一个nk
• 令$b=n_k$和$c=n_{k-1}$
• 重复步骤3和步骤4, 当算法停止时，所找到的b点的序列就构成了排列后的边界点的集合
  

1.2 链码

链码通过具有指定长度和方向的直线段的顺次连接来表示边界。该表示基于线段的4连接或8连接。每个线段的方向用数字编码方案编码。也称为佛雷曼链码。
链码的一次差分为，通过计算链码中分隔两个相邻像素的方向变化的数（按逆时针方向计算前一个数字变化到后一个数字所需的步数），其中第一元素是通过使用链码的最后一个元素和第一个元素间的转变来计算得到的。例如，4方向链码0231，得到的一次差分为3212，其中3为1逆时针转到0需3步；2为0逆时针转到2需2步；1为2逆时针转到3需1步；2为3逆时针转到1需2步。
边界的链码取决于起点。
起点归一化：将链码视为方向号码的一个循环序列，并重新定义起点，选择其中最小值整数序列作为最终编码（将编码作为自然数进行比较，前面0最多的就是最小的）。
旋转归一化：使用链码的一次差分。

• 若差分是通过计算链码中分割两个相邻像素的方向变换次数得到的。如4方向链码（按逆时针方向）10103322的一次差分为3133030。具体算法：逆时针从1->0需要3步，逆时针从0->1需要1步，逆时针从1->0需要3步，逆时针从0->3需要3步，以此类推。
• 若将链码作为循环序列，则差分的第一个元素是通过使用链码的最后一个元素和第一个元素间的转变得到的，此时4方向链码（按逆时针方向）10103322的一次差分为33133030。具体算法：逆时针从2->1需要3步，再接着依次计算从1->0，从0->1等等。

实验：链码表示边界
预处理

img = cv2.imread('D:/Data/school.jpg')

gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

ret, binary =cv2.threshold(gray,127,255,cv2.THRESH_BINARY)

contours, hierarchy =cv2.findContours(binary,cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)

cv2.drawContours(img,contours,-1,(0,0,255),3)

差分归一化：

columns = [0]

for i in range(81):

     columns.append(contours[0][i]-contours[0][i - 1])

     print (len(columns))

     print (columns[1][0][0])

a = []

for i in range(81):

   if columns[i][0][0] == 0 and columns[i][0][1] ==  -1:

       a.append(6)

   elif columns[i][0][0] == 0 and columns[i][0][1] ==  1:

       a.append(2)

   elif columns[i][0][0] == 1 and columns[i][0][1] ==  1:

       a.append(1)

   elif columns[i][0][0] == 1 and columns[i][0][1] ==  0:

       a.append(0)

   elif columns[i][0][0] == 1 and columns[i][0][1] ==  -1:

       a.append(7)

   elif columns[i][0][0] == -1 and columns[i][0][1] ==  1:

       a.append(3)

   elif columns[i][0][0] == -1 and columns[i][0][1] ==  0:

       a.append(4)

   elif columns[i][0][0] == -1 and columns[i][0][1] ==  -1:

       a.append(5)

print(a)

处理效果：
结果分析：
对于不同情境下的边缘检测连接效果有区别，过亮和过暗的效果都不是特别精准。

1.3 使用最小周长多边形的多边形近似（MPP）

一副图像我们可以用多边形拟合它的形状
目的就是使用合适的最大可能单元大小，以最少的定点数来生产MPP，并且我们可以发现构成边界的多边形每个转向要么是一个凸定点，要么是一个凹顶点。
注意的是凹顶点（黑色）有一个对应的“镜像”顶点，位于凹顶点的对角处
因此，设计的算法就只需要关注这些顶点

MPP算法
令W（白）和B（黑）分别表示凸顶点镜像凹顶点
寻找MMP的算法使用两个“爬行”点：

• 白色的爬行垫（$W_c$​）：$W_c$​沿凸顶点(w)爬行
• 黑色的爬行垫（$B_c$​）：$B_c$​​沿镜像凹顶点（B）爬行。

算法首先令$W_c=B_c=V_0$（$V_0$​是一个MPP顶点），$V_L$​表示最后一个MPP顶点，$V_k$​表示正在考察的当前顶点。
$sgn(a,b,c)=det(A)$（$det(A)$是A的行列式）
在$V_L$和$V_k$两个爬行点之间存在如下三个条件之一：

1. $V_k$​位于通过$(V_L{V}_k)$的直线的正的一侧，即$sgn(W_L,W_c,V_k)>0$

2. $V_k$​​位于通过$(V_L{W}_c)$的直线的负的一侧，或$V_k$与$(V_L{W}_c)$共线，即$sgn(W_L,W_c,V_k)\le 0$。同时，$V_k$​​位于通过$(V_L{B}_c)$的直线的正的一侧，或者$V_k$​​与$(V_L{B}_c)$共线，即$sgn(W_L,W_c,V_k)\ge 0$

3. $V_k$​​位于通过$(V_L{B}_c)$的直线的负的一侧，即$sgn(W_L,W_c,V_k)<0$

如果条件1成立，则下一个MPP顶点是$W_c$​​并且我们令$V_L=W_c$​；然后我们令$W_c=B_c=V_L$​来重新初始化该算法，并在$V_L$​之后的下一个顶点继续执行算法

如果条件2成立，则$V_k$​​k​编程一个候选的MPP顶点。在这种情况下，如果$V_k$​​是凸顶点，则令$W_c=V_k$；否则，我们置
$B_c=V_k$​然后，我们使用列表中的下一顶点继续执行算法

如果条件3成立，则下一个MPP顶点是$B_c$，并且令$V_L=B_c$​；然后我们令$W_c=B_c=V_L$​来重新初始化该算法，并用$V_L$之后的下一个顶点继续执行算法

当算法再次到达第一个顶点时，算法结束，此时算法已经处理了多边形中的所有顶点。

1.4 边界标记

标记是边界的一维函数 表示。一种最简单的生成方法是以角度函数形式画出质心到边界的距离。如下图所示
无论如何生成标记，基本思想都是将边界表示简化为比原始二维边界更简单的一维函数。
标记实验：

    def MergeEdage(self,savepath):
        leftmatrix = np.array([[1,0,-1],[1,0,-1],[1,0,-1]]) 
        left = self.Filter(leftmatrix)
        rightmatrix = np.array([[-1,0,1],[-1,0,1],[-1,0,1]]) 
        right = self.Filter(rightmatrix)
        w_1,h_1 = left[0].shape# 
        full_edage = []
        for i in range(3):
            m = np.hstack((left[i][:,:int(w_1/2)],right[i][:,w_1-int(w_1/2):]))
            full_edage.append(m)

效果 ：

效果分析：
对于图像边界标记的效果来说，可以有效标记出文字这些有明显集体特征的图像，对于建筑物来说，棱角分明且线条较为粗壮的标记效果也很好，对于一些物体和背景颜色接近，边界效果不是很好的，则还需要改进。

1.5 边界线段

将边界分解为线段通常是很有用的。分解降低了边界的复杂性，从而简化了描述过程。当边界线包含一个或多个携带形状信息的明显凹度时，该方法很有用。使用由边界所围成区域的凸壳就成为边界鲁棒分解的有力工具。
集合S的凸壳H是包含S的最小凸集，集合的差H-S被称为S的凸缺D。
分段算法：追踪S的轮廓，并标记进入或离开一个凸缺的转变点，从打标记的位置进行分段。
基本方法： 将边界分成若干段，然后分别对每一段进行表示，从而降低了边界的复杂度，并简化表示过程，尤其是当边界具有多个凹点的时候这种方法更为有效。
基本步骤：

1. 构造边界的凸包
2. 跟踪区域凸包的边界，记录凸包边界进出区域的转变点即可实现对边界的分割

1.6 骨架

基本思想：表示一个平面区域结构形状的重要方法是把 它削减成图形。这种削减可以通过细化（也称为抽骨架）算法，获取区域的骨架来实现 。
Blum的中轴变换方法（MAT）：设:R是一个区域，B为R的边界点，对于R中 的点p，找p在B上“近”的邻居。如果p有多于一个的邻居，称它属于R的中轴（骨架）。

问题：计算量大 ，包括计算区域的每个内部点到其边界点的距离。

算法改进思想 ：

在保证产生正确骨架的同时，改进算法的效率。比较典型的是一类细化算法，它们不断删除区域边界点，但保证删除满足：

• 不删除端点

• 不破坏连通性

• 不造成对区域的过度腐蚀

函数binary：

ef to_binary(img):
   
    w, h = img.shape
    i_two = copy.deepcopy(img)
    for i in range(w):
        for j in range(h):
            if img[i, j] < 200:
                i_two[i, j] = 0
            else:
                i_two[i, j] = 255

    return i_two

该函数用来生成二值图像中所有区域的骨骼。
细化函数，根据算法，运算出中心点的对应值：

def h_thin(img):
    h, w = img.shape
    i_next = 1
    for j in range(w):
        for i in range(h):
            if i_next == 0:
                i_next = 1
            else:
                i_m = int(img[i -1, j]) + int(img[i, j]) + int(img[i + 1, j]) if 0 < i < h - 1 else 1
                if img[i, j] == 0 and i_m != 0:
                    a = [0] * 9
                    for k in range(3):
                        for l in range(3):
                            if -1 < (i - 1 + k) < h and -1 < (j - 1 + l) < w and img[i - 1 + k, j - 1 + l] == 255:
                                a[k * 3 + l] = 1
                    i_sum = a[0] * 1 + a[1] * 2 + a[2] * 4 + a[3] * 8 + a[5] * 16 + a[6] * 32 + a[7] * 64 + a[8] * 128
                    img[i, j] = l_array[i_sum] * 255
                    if l_array[i_sum] == 1:
                        i_next = 0

处理效果：
例1
例2

结果分析：
平面区域的结构形状被削减成图形提取出物体主要构成部分的骨架，对于各个组成部分也有一定的处理抽取出主要结构部分。

（二）边界描绘子

2.1 一些简单的描述子

边界的周长： 是简单的描述符之一。沿轮廓线计算像素的个数，给出了一个长度的近似估计。

边界的直径： $Diam(B)=ma{x}_{i，j}[D(pi,pj)]$
$p_i,p_j$是边界上的点，$D(p_i,p_j)$是定义$p_i$​和$p_j$之间的距离。直径的值和连接组成该直径两个端点的直线段（该直线称为边界的长轴）的方向是边界的有用描述子，边界的短轴定义为与长轴垂直的直线，且由边界与两个轴相交的4个外部点所组成的方框（该方框称为基本矩形），可以全完包围该边界。
边界上的最远点并不总是唯一，如圆和方形上的点。但通常会假设：直径是一个有用的描绘子，则最好应用到具有单个最远点对的边界。
连接边界最远点对的线段称为边界的长轴。边界短轴 定义为与长轴垂直的线段。长轴和短轴所形成的矩形完全包含了边界，该矩形称为基本矩形。长轴和短轴的比率是边界的偏心率。
偏心率也是一个有用的描绘子。
曲率： 定义为斜率的变化率。近似：用相邻边界线段的斜率差作为这两条线段交点处的曲率。

2.2 形态数

形态数定义： 最小循环首差链码

循环首差链码： 用相邻链码的差代替链码
边界的形状数一般是以4连接的佛雷曼链码为基础的，形状数被定义为最小值整数序列的一阶差分，形状数的阶n定义为表示形状数的数字的个数。 边界的形状数可以用函数fchcode中的c.diffmm给出，形状数的阶n使用length(c.diffmm)得到。对于闭合边界，n为偶数。

2.3 傅里叶描述子

图显示了xy平面内的一个K点数字边界，从任意点$（x_0,y_0)$开始，以逆时针方向在该边界上行进时，会遇到坐标对$（x_0,y_0）,(x_1,y_1),(x_2,y_2),\dots ,(x_{k-1},y_{k-1})$这些坐标可以表示为$x(k)=x_k,y(k)=y_k$​的形式，此外，每个坐标都可当做一个复数来处理:
$s(k)=x(k)+jy(k)$ ,式中有k=0,1,2,…,K−1

$s(k)$的傅里叶变换:$a(u)={\sum }_{k=0}^{k-1}s(k)e^{\frac{-j2\pi uk}{K}}$

式中u=0,1,2,…,K−1u=0,1,2,…,K−1复系数a(u)称为边界的傅里叶描述子。这些系数的傅里叶反变换可恢复$s(k)$:

$s(k)=\frac{1}{K}{\sum }_{k=0}^{k-1}s(k)e^{\frac{-j2\pi uk}{K}}$

然而，假设仅使用前P个傅里叶系数而不使用所有系数，等同于上式中令$a(u)=0,u>P-1$结果$s(k)$的如下近似：

$\widehat{s(k)}=\frac{1}{P}{\sum }_{k=0}^{P-1}s(k)e^{\frac{-j2\pi uk}{P}}$

尽管求$\widehat{s(k)}$的每个成分时仅使用了P项，但k的范围仍然是从0到K-1，也就是说，在近似边界中存在同样数量的点，但项数不像在每个点的重建中那么多。
由于高频成分说明精细细节，而低频成分决定全局形状，因此P越小，边界丢失的细节就越多。

2.4 统计矩

如图a，它显示了一段边界，图b显示了以任意变量r的一维函数$g(r)$描述的线段，该函数是这样获得的：先将该线段的两个端点连接，然后旋转该直线，直至其为水平线段，此时所有的点的坐标也旋转同样的角度

把$g(r)$归一化为单位面积，并把它当作直方图来处理，换句话说，$g(r_i)$现在作为值$r_i$出现的概率来处理，此时，r为一个随机变量，故n阶矩为

$u_n（r)={\sum }_{i=0}^{k-1}(r_i-m{)}^{n}g(r_i)$

式中，
$m={\sum }_{i=0}^{k-1}{r}_{i}g(r_i)$

K是边界上的点数，$u_n（r)$直接与 $g(r)$的形状相关。例如，二阶矩$u_2(r)$度量曲线关于r的均值扩展程度，而三阶矩$u_3(r)$度量曲线关于均值的对称性。
在文献中n阶矩通常用符号${\mu }_n$表示，直接使用变量计算的矩被称为原始矩（raw moment），移除均值后计算的矩被称为中心矩（central moment）。变量的一阶原始矩等价于数学期望（expectation）、二至四阶中心矩被定义为方差（variance，度量曲线关于r的均值的扩展程度）、偏度（skewness，度量曲线关于均值的对称性）和峰度（kurtosis）。

（三）区域描绘子

3.1 简单的区域描述子

• 区域的面积： 区域的面积定义为该区域中像素的数量

• 区域的周长： 区域的周长是其边界的长度

• 致密性: 致密性定义为$\frac{周{长}^2}{面积}$

• 圆周率: 圆周率即一个区域的面积与具有相同周长的一个圆（最致密形状）的面积之比。
  周长为P的一个圆的面积$\frac{P^2}{4\pi }$。因此，圆周率$R_c$​由下式给出：
  $$R_c=\frac{4\pi A}{P^2}$$

• 其他:用作区域描述子的其他简单测度包括灰度值的均值和中值，最小灰度值和最大灰度值，以及其高于和低于均值的像素数

3.2 拓扑描绘子

拓扑特性对于图像平面区域的整体描述很有用。拓扑学研究未受任何变形影响的图形的性质，前提是该图形未被撕裂或粘连。
如果一个拓扑描绘子由该区域内的孔洞数量来定义，那么这种性质明显不受拉伸或旋转变换的影响。拓扑特性与距离或基于距离测度概念的任何特性无关。
另一个对区域描述有用的拓扑特性是连通分量的数量。
例如，图(a)显示了一个带有两个孔洞的区域。如果一个拓扑描述子由该区域内的孔洞数量来定义，那么这种性质明显不受拉伸或旋转变换的影响。
另一个对区域描述有用的拓扑特性是连通分量的数量。如图4(b)显示了一个具有3个连通分量的区域。
图形中孔洞的数量 H 和连通分量的数量 C ，可用于定义欧拉数 E ：E = C - H

欧拉数也是一种拓扑特性。例如，图所示的区域有分别等于 0 和 −1 的欧拉数。

3.3 纹理

描绘区域的一种重要方法是量化该区域的纹理内容。图像处理中用于描绘区域纹理的三种主要方法是统计法、结构法和频谱法。

1. 统计法：获得诸如平滑、粗糙、粒状等纹理特征。
2. 结构法：处理图像像元的排列，如基于规则间距平行线的纹理描述。
3. 频谱法：基于傅里叶频谱的特性，主要用于检测图像中的全局周期性，方法是识别频谱中的高能量的窄波峰。

统计法： 该方法是基于灰度直方图的统计特性。令z表示灰度的一个随机变量，令p(zi)，i = 0,1,2,…，L-1为相应的直方图，L为不同灰度级的数量。关于均值的z的第n阶矩为：

${\mu }_n（z）={\sum }_{i=0}^{L-1}(Z_i-m{)}^{n}P(Z_i)$

其中，m是z的均值（平均灰度）：

$m={\sum }_{i=0}^{L-1}{Z}_{i}P(Z_i)$
其中， ${\mu }_0=1$,${\mu }_1=0$。二阶矩${\mu }_2(z)$为方差${\sigma }^2$在纹理描述中特别重要。它是灰度对比度的度量，可用于建立相对平灰度的描述子。
对于恒定灰度区域为0（该区域方差为0），而对于较大的$sigm{a}^2(z)$值，其接近于1，因为对灰度级图像方差值增大而增大的。

三阶矩是直方图偏斜度的度量，而四阶矩是直方图相对平坦度的度量，五阶矩和更高阶矩不容易与直方图形状联系起来，但它们的确提供了纹理内容的进一步量化辨别。
使用直方图计算得到的纹理度量不懈怠像素彼此之间的相对位置的信息，但是这个信息很重要，所以我们使用一种方法叫共生矩阵。

令Q是定义两个像素彼此相对位置的一个算子，并考虑一副具有L个可能灰度级的图像f，令G为一个矩阵，其元素$g_{ij}$​是灰度为$Z_i$和$Z_j$​的像素对出现在f中有Q所指定的位置处的次数，按照这种方法形成的矩阵成为灰度共生矩阵。
通俗来讲就是：找到图像f中可能的灰度级数L，生成一个L*L 的矩阵G，矩阵中的元素是看其行号和列号在图像f中有没有Q指定的相对位置，出现几次，该位置就是几。如矩阵G的（6,2）位置，对应的图像像素出现了三次（像图中圈出的那样），则该位置的值为3，其余的类似。

$P_{ij}=g_{ij}/n$

满足Q的像素对的总数n，等于G的元素之和。故在这里插入图片描述
是满足Q的一个值为(zi,zj)的点对的概率估计。这些概率的值域为[0,1]，且它们的和为1。
结构方法

结构法的基本思想是，一个简单的“纹理基元”可借助一些规则用于形成更复杂的纹理模式，这些规则限制基元（或这些基元）的可能排列的数量。

频谱方法

对纹理描述有用的傅里叶频谱的三个特征：

1. 频谱中突出的尖峰给出纹理模式的主要方向；
2. 频率平面中尖峰的位置给出模式的基本空间周期；
3. 采用滤波方法消除任何周期成分而留下非周期性图像元素，然后采用统计技术来描述。

3.4 不变矩

几何矩是由Hu(Visual pattern recognition by moment invariants)在1962年提出的，具有平移、旋转和尺度不变性。
大小为M*N的数字图像f(x.y)的二维(p + q)阶矩定义：

$$m_{pq}=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}{x}^p{y}^{q}f(x,y)$$

p和q均为0,1,2,…的整数。相应的(p + q)阶中心矩定义:

$${\mu }_{pq}=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}(x-\bar{x}{)}^p(y-\bar{y}{)}^{q}f(x,y)$$
其中：
$$\bar{x}=\frac{m_{10}}{m_{00}},\bar{y}=\frac{m_{01}}{m_{00}}$$

归一化的中心矩定义:
$${\eta }_{pq}=\frac{{\mu }_{pq}}{{\mu }_{00}^r}$$
其中：
$$\gamma =\frac{p+q}{2}+1,p+q=2,3...$$
利用二阶和三阶归一化中心矩构造了7个不变矩:

这7个不变矩构成一组特征量，Hu.M.K在1962年证明了他们具有旋转，缩放和平移不变性。

实际上，在对图片中物体的识别过程中，只有M1 和 M2不变性保持的比较好，其他的几个不变矩带来的误差比较大，有学者认为只有基于二阶矩的不变矩对二维物体的描述才是真正的具有旋转、缩放和平移不变性（ 和刚好都是由二阶矩组成的）。
由Hu矩组成的特征量对图片进行识别，优点就是速度很快，缺点是识别率比较低，这一部分原因是由于Hu不变矩只用到低阶矩（最多也就用到三阶矩），对于图像的细节未能很好的描述出来，导致对图像的描述不够完整。
Hu不变矩一般用来识别图像中大的物体，对于物体的形状描述得比较好，图像的纹理特征不能太复杂，像识别水果的形状，或者对于车牌中的简单字符的识别效果会相对好一些。

(四) 使用主成分进行描绘

4.1 概念介绍

本节内容适用于边界和区域。也是描绘一组空间上已配准图像的基础。
主成分变换，是指由原始图像数据协方差矩阵的特征值和特征向量建立起来的变换核，将光谱特征空间原始数据向量投影到平行于地物集群椭球体各结构轴的主成分方向，突出和保留主要的物类别信息，用来进行图像增强、特征选择和图像压缩的处理方法。
假设有n幅已配准的图像，他们堆叠在一起。如下图所示：
对任意给定的坐标(i,j)，都有n个像素，每一幅图像在该位置上都有一个像素。这些像素以列向量的形式排列：

若每个图像的大小都为MN，则在n幅图像中，包含所有图像的n维向量共有MN个，即图像的个数决定向量的维数，图像的大小决定向量的个数
K（K = MN）个向量的均值向量mx可通过样本的平均值来近似：
$$m_x=\frac{1}{k}\sum _{k=1}^k{x}_k$$
为了从样本中获得Cx的无偏估计，用K-1代替K。因为Cx是一个实对称矩阵，故总可以找到n个正交特征向量。
主分量变换（霍特林变换）如下：
$$y=A(x-m_k)$$

A是由Cx的特征向量组成，特征向量的排序方式为A的第一行对应于最大特征值的特征向量，最后一行对应于最小特征值的特征向量。
y的均值my为零。y的元素是不相关的，故y的协方差矩阵Cy是对角阵。主对角线元素就是Cx的特征值。故Cx和Cy有相同的特征值。
因为矩阵A的行向量是正交的，故$A^{-1}=A^T$。任何向量x都能通过下式由其对应的y来恢复：
$$x=A^{T}y+m_x$$

然而往往不使用Cx的所有特征向量，而是有对应的前k个最大特征值的k个特征x向量来形成矩阵$A_k$，得到一个k*n的变换矩阵，向量y成为k维矩阵。这就是主成分变换的魅力所在。
由$A_k$重建的向量x（是样本x是近似）由下式表示
$$\hat{x}=A_k^{T}y+m_x$$

4.2 图像的主成分分析（PCA）和pickle模块的使用

由于图像具有很高的维数，在许多计算机视觉应用中，需要经常使用降维操作。PCA（Principal Component Analysis，主成分分析）可以在使用尽可能少维数的前提下，尽量多地保持训练数据的信息。思想是把数据投影到方向向量使数据集的特征向量到方向向量的垂线长度最短。
PCA 产生的投影矩阵可以被视为将原始坐标变换到现有的坐标系，坐标系中的各个坐标按照重要性递减排列。

为了对图像数据进行 PCA 变换，变换前，图像需要转换成一维向量表示，使用 NumPy 类库中的 flatten() 方法进行变换。
PCA计算步骤：

1. 去平均 （即为数据中心化）

2. 计算协方差矩阵

3. 计算协方差矩阵的特征向量和特征值

4. 将特征值从小到大排列

5. 保留最上面的n个特征向量

6. 将数据转换到上述n个特征向量构建新的空间

其中，降维过的方向向量就是新的数据集、特征向量就是数据集维度个数，通过特征向量个数限制降维的维数。
通常使用 SVD（Singular Value Decomposition，奇异值分解）方法来计算主成分；但当矩阵的维数很大时， SVD 的计算非常慢，所以此时通常不使用 SVD 分解。如果数据个数小于向量的维数，不用 SVD 分解，而是计算维数更小的协方差矩阵X X^{T} 的特征向量。通过仅计算对应前 k（k 是降维后的维数）最大特征值的特征向量，可以使上面的 PCA 操作更快。

代码如下：

def pca(X):
    """ 主成分分析：
输入：矩阵X，其中该矩阵中存储训练数据，每一行为一条训练数据
返回：投影矩阵（按照维度的重要性排序）、方差和均值 """
    # 获取维数
    num_data,dim = X.shape
    # 数据中心化
    mean_X = X.mean(axis=0)
    X = X - mean_X
 
    if dim > num_data:
        # PCA- 使用紧致技巧
        M = dot(X, X.T)  # 协方差矩阵
        e, EV = linalg.eigh(M)  # 特征值和特征向量
        tmp = dot(X.T, EV).T  # 这就是紧致技巧
        V = tmp[::-1]  # 逆转
        S = sqrt(e)[::-1]  # 逆转
        for i in range(V.shape[1]):
            V[:, i] /= S
    else:
        # PCA- 使用 SVD 方法
        U, S, V = linalg.svd(X)
        V = V[:num_data]  # 仅仅返回前 nun_data 维的数据才合理
    return V, S, mean_X

需要保存一些结果或者数据以方便后续使用，Python 中的 pickle 模块非常有用。pickle 模块可以接受几乎所有的 Python 对象，并且将其转换成字符串表示，该过程叫做封装（pickling）。从字符串表示中重构该对象，称为拆封（unpickling）。 这些字符串表示可以方便地存储和传输。

举例说明：

import pickle
 
dataList = [[1, 1, 'yes'],
            [1, 1, 'yes'],
            [1, 0, 'no'],
            [0, 1, 'no'],
            [0, 1, 'no']]
dataDic = {0: [1, 2, 3, 4],
           1: ('a', 'b'),
           2: {'c': 'yes', 'd': 'no'}}
 
# 使用dump()将数据序列化到文件中
fw = open('dataFile.pkl', 'wb')
# Pickle the list using the highest protocol available.
pickle.dump(dataList, fw, -1)
# Pickle dictionary using protocol 0.
pickle.dump(dataDic, fw)
fw.close()
 
# 使用load()将数据从文件中序列化读出
print('.............  open() close()  ................')
fr = open('dataFile.pkl', 'rb')
data1 = pickle.load(fr)
print(data1)
data2 = pickle.load(fr)
print(data2, '\n')
fr.close()
 
# 或者使用with()思想,打开文件并保存
with open('dataFile_01.pkl', 'wb') as f:
    pickle.dump(dataList, f, -1)
    pickle.dump(dataDic, f)
# 或者使用with()思想,打开文件并读取
with open('dataFile_01.pkl', 'rb') as f:
    dt_01 = pickle.load(f)
    dt_02 = pickle.load(f)
    print('.............  with()  ................')
    print(dt_01)
    print(dt_02, '\n')
 
# 使用dumps()和loads()举例
print('.............  使用dumps()和loads()举例  ................')
p = pickle.dumps(dataList)
print(pickle.loads(p))
p = pickle.dumps(dataDic)
print(pickle.loads(p))

修改部分数据再测试：

if __name__ == '__main__':
    # 产生三类数据
    x, y = make_blobs(n_samples=600, n_features=3, centers=3)
    # 绘制三维图像
    fig = plt.figure(figsize=(8, 6))
    ax = Axes3D(fig)
    ax.scatter(x[:, 0], x[:, 1], x[:, 2], c='r')
    plt.savefig('3-D distribution of Datapoints')
    plt.show()
    #  # 绘制二维图像
    fig = plt.figure(figsize=(8, 6))
    ax1 = fig.add_subplot(2, 2, 1)
    ax1.scatter(x[:, 0], x[:, 1], c=y)
    ax2 = fig.add_subplot(2, 2, 2)
    ax2.scatter(x[:, 1], x[:, 2], c=y)
    ax3 = fig.add_subplot(2, 2, 3)
    ax3.scatter(x[:, 0], x[:, 2], c=y)
    plt.savefig('2-D distribution of Datapoints')
    plt.show()

    #  step1 中心化
    x_mean = np.mean(x, axis=0)
    x_center = x - x_mean
    #  step2 计算协方差矩阵
    covX = np.cov(x_center.T)
    #  step3 对协方差矩阵进行特征分解
    featureValue, featureVector = np.linalg.eig(covX)
    #  step4 取最大的几个特征值
    index = np.argsort(-featureValue)
    selectVec = featureVector.T[index[:2]]
    print(selectVec.shape)
    final_matrix = np.dot(x_center, selectVec.T)
    print(final_matrix.shape)

    # 绘制最终二维图像
    fig = plt.figure(figsize=(8, 6))
    plt.scatter(final_matrix[:, 0], final_matrix[:, 1], c=y)
    plt.title('2-D distribution after PCA')
    plt.savefig('2-D distribution of Datapoints after PCA')
    plt.show()

PCA模块可以降维对象，并且将其转换成字符串表示，封装和拆封很顺利，字符串也更方便存储和传输。