PCA人脸识别

“PCA人脸识别”

        做人脸识别研究的，PCA人脸识别应该必须了解的，这是基础，因为真正的计算机自动人脸识别的研究就是从PCA人脸识别开始的，20世纪80~90年代——经典PCA(Eigenface)和LDA(Fisherface)相继诞生。我们所研究的人脸识别，一般都基于标准人脸库进行训练、测试识别。

1.PCA人脸识别方法

将PCA方法用于人脸识别，其实是假设所有的人脸都处于一个低维线性空间，而且不同的人脸在这个空间中具有可分性。其具体做法是由高维 图像空间经PCA变换后得到一组新的正交基，对这些正交基做一定的取舍，保留其中的一部分生成低维的人脸空间，也即是人脸的特征子空间。PCA人脸识别算法步骤包括：

a.人脸图像预处理

b.读入人脸库，训练形成特征子空间 【特征值、特征向量的求法】

c.把训练图像和测试图像投影到上一步骤中的特征子空间上 【矩阵相乘】

d.选择一定的距离函数进行判别  【欧氏距离，挑最小的匹配】

2.PCA人脸识别流程

a.读入人脸库，读入每一个二维的人脸图像并转化为一维的向量，每个人选定一定数量的人脸照片构成训练集【共20张】，则训练集是一个36000*20的矩阵。测试集共10张图像，每次选一张，则测试集是一个36000*1的矩阵。

样本集：

测试集：

代码：

[cpp]  view plain  copy

 

1. void load_data(double *T,IplImage *src,int k)  
2. {  
3.     int i,j;  
4.   
5.     //一副图像压缩成一维的，存在T的一列里  
6.     for (i=0;i<IMG_HEIGHT;i++)  
7.     {  
8.         for (j=0;j<IMG_WIDTH;j++)  
9.         {  
10.             T[(i*IMG_WIDTH+j)*TRAIN_NUM+k-1]= (double)(unsigned char)src->imageData[i*IMG_WIDTH+j];  
11.         }  
12.     }  
13. }  

b.计算 PCA变换的生成矩阵Q。首先计算训练集的协方差矩阵X，其中x1,x2,...,xn为第i副图像的描述，即xi为一个36000*1的列向量。

，

由于这个矩阵太大36000*36000，求特征值和特征向量比较坑，所以改为求 P=X^TX 的特征向量和特征值，且有如下性质：

设e是矩阵P的特征值λ对应的特征向量，则有：

这里，X*e也是矩阵Q的特征值λ对应的特征向量，可以如此变换。

代码：

[cpp]  view plain  copy

 

1. void calc_mean(double *T,double *m)  
2. {  
3.     int i,j;  
4.     double temp;  
5.   
6.     for (i=0;i<IMG_WIDTH*IMG_HEIGHT;i++)  
7.     {  
8.         temp=0;  
9.         for (j=0;j<TRAIN_NUM;j++)  
10.         {  
11.             temp = temp + T[i*TRAIN_NUM+j];  
12.         }  
13.         m[i] = temp/TRAIN_NUM;  
14.     }  
15. }  
16.   
17. void calc_covariance_matrix(double *T,double *L,double *m)  
18. {  
19.     int i,j,k;  
20.     double *T1;  
21.   
22.     //T = T -m  
23.     for (i=0;i<IMG_WIDTH*IMG_HEIGHT;i++)  
24.     {  
25.         for (j=0;j<TRAIN_NUM;j++)  
26.         {  
27.             T[i*TRAIN_NUM+j] = T[i*TRAIN_NUM+j] - m[i];  
28.         }  
29.     }  
30.   
31.     T1 = (double *)malloc(sizeof(double)*IMG_HEIGHT*IMG_WIDTH*TRAIN_NUM);  
32.   
33.     //L = T' * T  
34.     matrix_reverse(T,T1,IMG_WIDTH*IMG_HEIGHT,TRAIN_NUM);  
35.     matrix_mutil(L,T1,T,TRAIN_NUM,IMG_HEIGHT*IMG_WIDTH,TRAIN_NUM);  
36.   
37.     free(T1);  
38. }  

c.计算生成矩阵P的特征值和特征向量，并挑选合适的特征值和特征向量，构造特征子空间变化矩阵。这里P是实对称矩阵，可以先进行Household变换将P变成对角矩阵，然后使用QR迭代算法求解特征值和特征向量，迭代次数60，误差eps=0.000001，代码：

[cpp]  view plain  copy

 

1. void cstrq(double a[],int n,double q[],double b[],double c[])  
2. {  
3.     int i,j,k,u,v;  
4.     double h,f,g,h2;  
5.     for (i=0; i<=n-1; i++)  
6.         for (j=0; j<=n-1; j++)  
7.         { u=i*n+j; q[u]=a[u];}  
8.         for (i=n-1; i>=1; i--)  
9.         { h=0.0;  
10.         if (i>1)  
11.             for (k=0; k<=i-1; k++)  
12.             { u=i*n+k; h=h+q[u]*q[u];}  
13.             if (h+1.0==1.0)  
14.             { c[i]=0.0;  
15.             if (i==1) c[i]=q[i*n+i-1];  
16.             b[i]=0.0;  
17.             }  
18.             else  
19.             { c[i]=sqrt(h);  
20.             u=i*n+i-1;  
21.             if (q[u]>0.0) c[i]=-c[i];  
22.             h=h-q[u]*c[i];  
23.             q[u]=q[u]-c[i];  
24.             f=0.0;  
25.             for (j=0; j<=i-1; j++)  
26.             { q[j*n+i]=q[i*n+j]/h;  
27.             g=0.0;  
28.             for (k=0; k<=j; k++)  
29.                 g=g+q[j*n+k]*q[i*n+k];  
30.             if (j+1<=i-1)  
31.                 for (k=j+1; k<=i-1; k++)  
32.                     g=g+q[k*n+j]*q[i*n+k];  
33.             c[j]=g/h;  
34.             f=f+g*q[j*n+i];  
35.             }  
36.             h2=f/(h+h);  
37.             for (j=0; j<=i-1; j++)  
38.             { f=q[i*n+j];  
39.             g=c[j]-h2*f;  
40.             c[j]=g;  
41.             for (k=0; k<=j; k++)  
42.             { u=j*n+k;  
43.             q[u]=q[u]-f*c[k]-g*q[i*n+k];  
44.             }  
45.             }  
46.             b[i]=h;  
47.             }  
48.         }  
49.         for (i=0; i<=n-2; i++) c[i]=c[i+1];  
50.         c[n-1]=0.0;  
51.         b[0]=0.0;  
52.         for (i=0; i<=n-1; i++)  
53.         { if ((b[i]!=0.0)&&(i-1>=0))  
54.         for (j=0; j<=i-1; j++)  
55.         { g=0.0;  
56.         for (k=0; k<=i-1; k++)  
57.             g=g+q[i*n+k]*q[k*n+j];  
58.         for (k=0; k<=i-1; k++)  
59.         { u=k*n+j;  
60.         q[u]=q[u]-g*q[k*n+i];  
61.         }  
62.         }  
63.         u=i*n+i;  
64.         b[i]=q[u]; q[u]=1.0;  
65.         if (i-1>=0)  
66.             for (j=0; j<=i-1; j++)  
67.             { q[i*n+j]=0.0; q[j*n+i]=0.0;}  
68.         }  
69.         return;  
70. }  
71.   
72. //q:特征向量，b：特征值  
73. int csstq(int n,double b[],double c[],double q[],double eps,int l)  
74. {  
75.     int i,j,k,m,it,u,v;  
76.     double d,f,h,g,p,r,e,s;  
77.     c[n-1]=0.0; d=0.0; f=0.0;  
78.     for (j=0; j<=n-1; j++)  
79.     { it=0;  
80.     h=eps*(fabs(b[j])+fabs(c[j]));  
81.     if (h>d) d=h;  
82.     m=j;  
83.     while ((m<=n-1)&&(fabs(c[m])>d)) m=m+1;  
84.     if (m!=j)  
85.     { do  
86.     { if (it==l)  
87.     { printf("fail\n");  
88.     return(-1);  
89.     }  
90.     it=it+1;  
91.     g=b[j];  
92.     p=(b[j+1]-g)/(2.0*c[j]);  
93.     r=sqrt(p*p+1.0);  
94.     if (p>=0.0) b[j]=c[j]/(p+r);  
95.     else b[j]=c[j]/(p-r);  
96.     h=g-b[j];  
97.     for (i=j+1; i<=n-1; i++)  
98.         b[i]=b[i]-h;  
99.     f=f+h; p=b[m]; e=1.0; s=0.0;  
100.     for (i=m-1; i>=j; i--)  
101.     { g=e*c[i]; h=e*p;  
102.     if (fabs(p)>=fabs(c[i]))  
103.     { e=c[i]/p; r=sqrt(e*e+1.0);  
104.     c[i+1]=s*p*r; s=e/r; e=1.0/r;  
105.     }  
106.     else  
107.     { e=p/c[i]; r=sqrt(e*e+1.0);  
108.     c[i+1]=s*c[i]*r;  
109.     s=1.0/r; e=e/r;  
110.     }  
111.     p=e*b[i]-s*g;  
112.     b[i+1]=h+s*(e*g+s*b[i]);  
113.     for (k=0; k<=n-1; k++)  
114.     { u=k*n+i+1; v=u-1;  
115.     h=q[u]; q[u]=s*q[v]+e*h;  
116.     q[v]=e*q[v]-s*h;  
117.     }  
118.     }  
119.     c[j]=s*p; b[j]=e*p;  
120.     }  
121.     while (fabs(c[j])>d);  
122.     }  
123.     b[j]=b[j]+f;  
124.     }  
125.     for (i=0; i<=n-1; i++)  
126.     { k=i; p=b[i];  
127.     if (i+1<=n-1)  
128.     { j=i+1;  
129.     while ((j<=n-1)&&(b[j]<=p))  
130.     { k=j; p=b[j]; j=j+1;}  
131.     }  
132.     if (k!=i)  
133.     { b[k]=b[i]; b[i]=p;  
134.     for (j=0; j<=n-1; j++)  
135.     { u=j*n+i; v=j*n+k;  
136.     p=q[u]; q[u]=q[v]; q[v]=p;  
137.     }  
138.     }  
139.     }  
140.     return(1);  
141. }  
142.   
143. void matrix_reverse(double *src,double *dest,int row,int col)   //转置  
144. {  
145.     int i,j;  
146.   
147.     for(i = 0;i < col;i++)  
148.     {  
149.            for(j = 0;j < row;j++)  
150.            {  
151.          dest[i * row + j] = src[j * col + i];  
152.            }  
153.          }  
154. }  
155.   
156. void matrix_mutil(double *c,double *a,double *b,int x,int y,int z)  //矩阵乘法  
157. {  
158.     int i,j,k;  
159.     for (i=0;i<x;i++)  
160.     {  
161.         for (k=0;k<z;k++)  
162.         {  
163.             for (j=0;j<y;j++)  
164.             {  
165.                 c[i*z+k] +=a[i*y+j]*b[j*z+k];  
166.             }  
167.         }  
168.     }  
169. }  

挑选合适的特征值和特征向量，其实就是挑特征值大于1的【关于挑选，可以排序选前k个，也可以设阈值】：

[cpp]  view plain  copy

 

1. void pick_eignevalue(double *b,double *q,double *p_q,int num_q)  
2. {  
3.     int i,j,k;  
4.   
5.     k=0;//p_q的列  
6.     for (i=0;i<TRAIN_NUM;i++)//col  
7.     {  
8.         if (b[i]>1)  
9.         {  
10.             for (j=0;j<TRAIN_NUM;j++)//row  
11.             {  
12.                 p_q[j*num_q+k] = q[j*TRAIN_NUM+i];//按列访问q,按列存储到p_q  
13.   
14.             }  
15.             k++;  
16.         }  
17.     }  
18. }  

d.把训练图像和测试图像投影到特征空间中。每一幅人脸图像投影到子空间以后，就对应与子空间的一个点。同样，子空间中的任一点也对应于一副图像。这些子空间的点在重构以后的图像很像人脸，所以他们被成为特征脸Eigenface。有了这样一个由特征脸组成的降维子空间，任何一副人脸图像都可以向其做投影并获得一组坐标系数，这组系数表明了该图像在子空间中的位置，这样原来的人脸图像识别问题就转化为依据子空间的训练样本点进行分类的问题。

【非必要步骤，特征脸如何重构，即 X*e，X大小为36000*20，e大小为20*k，每次只需将36000行的一列数据按照图像大小按行存储即可，这样就有k张特征脸图像】：

 

[cpp]  view plain  copy

 

1. double  *temp;  
2.     IplImage *projected;  
3.     char res[20]={0};   //file name  
4.     temp = (double *)malloc(sizeof(double)*IMG_HEIGHT*IMG_WIDTH*num_q);//按列存取  
5.     projected = cvCreateImage(cvSize(IMG_WIDTH,IMG_HEIGHT),IPL_DEPTH_8U,1);  
6.     //求特征脸  
7.     matrix_mutil(temp,T,p_q,IMG_WIDTH*IMG_HEIGHT,TRAIN_NUM,num_q);  
8.       
9.     for (i=0;i<num_q;i++)  
10.     {  
11.         sprintf(res,"%d.jpg",i);  
12.         for (j=0;j<IMG_HEIGHT;j++)  
13.         {  
14.             for (k=0;k<IMG_WIDTH;k++)  
15.             {  
16.                 projected->imageData[j*IMG_WIDTH+k] = (unsigned char)abs(temp[(j*IMG_WIDTH+k)*num_q+i]);  
17.             }  
18.         }  
19.         cvSaveImage(res,projected);  
20.     }  

1. PCA人脸识别操作流程

    在平时的研究中，我总结的PCA人脸识别的主要流程如下图所示：

                                

                                                    图 1 PCA人脸识别流程图

      通过上图的PCA人脸识别流程可以看出，PCA方法可以总结为以下几个阶段：训练样本、特征提取、构造特征空间、投影计算。

2. PCA人脸识别方法原理介绍

       Karhunen-Loeve(K-L)变换或主成分分析(Principal Component Analysis，PCA)是普遍使用的一种技术,这种技术的主要作用是将有续信号转换为一组不相关的表示系数。主成分分析(PCA)形成了K-L变换的基础，主要用于数据的紧凑表示。在数据挖掘的应用中，它主要应用于简化大维数的数据集合，减少特征空间维数，可以用较小的存储代价和计算复杂度获得较高的准确性。

    PCA方法通过消除数据的相关性，找到一个空间，使得各个类别的数据在该空间上能够很好地分离。如下图所示：

                                              

                                                  图 2 PCA降维和分类示意图

在图2中，有一些离散的二维分布点，其中五角星表示一类集合，小圆圈表示另一类集合，假设这两个类别可以用特征X和特征Y进行描述，由图可知，在X轴和Y轴上这两个类别的投影是重叠的，表明这些点的两个特征X和Y没有表现出突出的识别性。但是两个类的投影在Z轴上区分度较大，显示出很好的识别性。PCA就是这样的一个工具，它可以产生非常好的降维效果，这种方法也可以用在图像处理的其他的一些研究中，如图像压缩、分类以及特征选择等。

    接下来详细介绍一下PCA的理论原理：

    通过对主成分分析方法做了简单的说明以后，接下来，我们来具体分析一下K-L变换的原理。

    假设在图像集合f(m,n)中，每张图像可以用堆叠的方式表示成一个mn维列向量：

                            　　　　　　　  (1)

其表示方法为：从图像矩阵的第一列开始取，依次取到最后一列，每一列首尾相连，构成一个mn维列向量。将每一幅人脸图像表示

列向量以后，然后依次将每一个图像列向量转置成行向量，构成一个人脸样本矩阵，该样本矩阵即图像集合：

                                                     (2)

由式(2)可知样本矩阵有L行，其中每一行的数据就代表一张人脸样本图像，L表示训练样本的总个数。  

训练样本的协方差矩阵为：

                              　　　　　　　　　　　　                                           (3)

式中mf是所有训练样本的平均值向量，也即是所有样本的平均脸。(3)式中的[Cf]阵为L阶实对称方阵，则其一定存在L个相互正交的属于各个特征值的特征向量，即有：

                          　　　                                             　 (4)

将上述特征值进行降序排列，并取每个特征值对应的特征向量构成一个正交矩阵，也即是一个L维的正交空间。按照文献中的说明，此处的矩阵[Cf]的维数很大，求解其特征值和特征向量比较复杂，这个时候需要对（3）式进行变形，简化求解，其变形结果如下：

                     　　　　　　　　　　　　　　                                      (5)

此时，对（5）式中的矩阵进行特征值和特征向量的求解，将求解出的特征向量和特征值经过SVD奇异值分解，得到原训练样本的特征向量，这样就可以构造出最终的人脸投影空间。

                    　　　　　　　　　　　                              （6）

其中vi就是（5）中协方差矩阵的特征向量，ｐ是特征向量的个数。将特征向量转化为矩阵，矩阵就可以表示图像，也即是所谓的特征脸，如下：

                        

                                              图 3 PCA特征脸

到这里为止，我们就找到了PCA人脸识别需要的投影特征空间了：

　　　　　　　　　　　　　　　　　　（７）

接下来就可以通过Ｋ－Ｌ变换式进行投影计算了，得出各个样本在空间Wpca上的投影特征了：

                                           　　　　　                                 　　　　　　　　　　       （8）

其中［A］就是空间Wpca，其实（8）式中不用（f-mf）也行，直接换成f就可以，即表示将原人脸样本在空间Wpca上进行投影，将投影后的特征系数存在矩阵g中，如果没有记错的话，g中每一行表示一个人脸样本的特征系数（或者每一列表示一个人脸样本的特征稀疏，这个与个人的计算方法有关，主要是看矩阵的转置或者不转置）。

    到这里PCA人脸识别的主要步骤已经介绍完了，剩下的就是识别过程了，这个比较容易理解，首先将训练样本在空间[A]上进行投影，得到投影样本的特征系数，然后将测试样本也在空间[A]上进行投影，得到每个测试杨样本的投影特征系数，此时，只需将测试某个样本的特征系数与训练样本投影特征系数进行欧式距离度量，看要测试的那个样本与训练集中哪个样本的欧式距离最近，就可以将该测试样本归为与之距离最近那个样本的类别。比如：训练样本a1属于S1类（S1类中包含很多样本，a1只是其中的一个样本），如果测试样本b1和a1的距离最近，那么就将b1归为S1类。如此下来，对所有的样本进行分类，就可以得到PCA人脸识别率了。