图像白化处理，协方差求解，PCA白化

在人脸表情所用的特征中，LBP特征中是一个常用的特征。但是LBP特征提取后，维数太高，即使使用LBP直方图的特征仍然很高，因此需要降维。一个比较好的降维方法就是PCA。在OpenCV中，cvCalcPCA函数可以直接调用进行PCA降维。以前在人脸识别的程序中，我曾经调用过这个函数，这次再次使用，但是却出错了。主要有两个问题：

1 特征值和特征向量的行列数设置问题

如果按如下方式使用，也就是根据PCA的数学表达式进行参数设置：

AvgVector=cvCreateMat(1,m_pm.dim,CV_32FC1);

EigenValue_Row=cvCreateMat(1,m_pm.dim,),CV_32FC1);

EigenVector=cvCreateMat((m_pm.dim),m_pm.dim,CV_32FC1);

那么就会出现错误：

Incorrect size of input array。

根据错误提示，说明特征向量的列数和输入矩阵（特征矩阵）的列数一致，行数要和特征值的数目一致；而特征值的数目（特征值矩阵的行数或列数）不能大于输入矩阵的行数和列数。

因此，需要下列设定：

AvgVector=cvCreateMat(1,m_pm.dim,CV_32FC1);

EigenValue_Row=cvCreateMat(1,min(m_pm.dim,m_pm.numOfSamples),CV_32FC1);

或EigenValue_Row=cvCreateMat(min(m_pm.dim,m_pm.numOfSamples),1,CV_32FC1);

EigenVector=cvCreateMat(min(m_pm.dim,m_pm.numOfSamples),m_pm.dim,CV_32FC1);

// 计算特征值，对原始数据进行变换，得到主成分

cvCalcPCA(vector_feature,AvgVector,EigenValue_Row,EigenVector,CV_PCA_DATA_AS_ROW);

cvProjectPCA(vector_feature,AvgVector,EigenVector,vector_pca);

 

2 出现中断点

由于保存输入矩阵的特征是一维的，需要转换成二维的，所以要申请一段内存保存特征，在cvCalcPCA之前这段内存不能释放，需要在cvCalcPCA之后才能释放，所以正确 的代码为：

float* features=new float[lenOfFeatures];

...

CvMat* vector_feature=cvCreateMat(m_pm.numOfSamples,m_pm.dim,CV_32FC1);

cvSetData(vector_feature,features,vector_feature->step);

CvMat *AvgVector;

CvMat *EigenVector;

CvMat *EigenValue_Row;

CvMat* vector_pca=cvCreateMat(m_pm.numOfSamples,PCA_DIM,CV_32FC1);

 

AvgVector=cvCreateMat(1,m_pm.dim,CV_32FC1);

EigenValue_Row=cvCreateMat(1,min(m_pm.dim,m_pm.numOfSamples),CV_32FC1);

EigenVector=cvCreateMat(min(m_pm.dim,m_pm.numOfSamples),m_pm.dim,CV_32FC1);

// 计算特征值，对原始数据进行变换，得到主成分

cvCalcPCA(vector_feature,AvgVector,EigenValue_Row,EigenVector,CV_PCA_DATA_AS_ROW);

cvProjectPCA(vector_feature,AvgVector,EigenVector,vector_pca);

delete[] features;

在图形识别方面,主成分分析(Principal Comonents Analysis,PCA)算是比较快速而且又准确的方式之一,它可以对抗图形平移旋转的事件发生,并且藉由主要特征(主成分)投影过后的数据做数据的比对,在多个特征信息里面,取最主要的K个,做为它的特征依据,在这边拿前面共变量矩阵的数据来做沿用,主成分分析使用的方法为计算共变量矩阵,在加上计算共变量矩阵的特征值及特征向量,将特征值以及所对应的特征向量排序之后,取前面主要K个特征向量当做主要特征,而OpenCV也可以对高维度的向量进行主成分分析的计算


数据原始的分布情况,红点代表着它的平均数


将坐标系位移,以红点为主要的原点


计算共变量矩阵以及共变量的特征值以及特征向量,将特征向量排序后投影回原始数据的结果的结果,也就是说,对照上面的图片,EigenVector的作用是找到主轴后,将原本的坐标系做旋转了


再来就是对它做投影,也就是降低维度的动作,将Y轴的数据全部归零,投影在X轴上


投影完之后,在将它转回原本的坐标系


PCA主成分分析,与线性回归有异曲同工之妙,也就是说,这条投影过后的直线,可以称做回归线,当它在做主轴旋转的时候,所投影的结果为最小均方误,在将它转置回来的时候,就形成了一条回归直线了

OpenCV的PCA输入必须要是单信道32位浮点数格式或是单信道64位浮点数格式的,参数为CV_32FC1或是CV_64FC1,程序写法如下

PCA程序实作

[cpp] view plaincopy

1. #include <cv.h>  
2. #include <highgui.h>  
3. #include <stdio.h>  
4. #include <stdlib.h>  
5.   
6. float Coordinates[20]={1.5,2.3,  
7.                                  3.0,1.7,  
8.                                  1.2,2.9,  
9.                                  2.1,2.2,  
10.                                  3.1,3.1,  
11.                                  1.3,2.7,  
12.                                  2.0,1.7,  
13.                                  1.0,2.0,  
14.                                  0.5,0.6,  
15.                                  1.0,0.9};  
16.   
17. void PrintMatrix(CvMat *Matrix,int Rows,int Cols);  
18.   
19. int main()  
20. {  
21.     CvMat *Vector1;  
22.     CvMat *AvgVector;  
23.     CvMat *EigenValue_Row;  
24.     CvMat *EigenVector;  
25.   
26.     Vector1=cvCreateMat(10,2,CV_32FC1);  
27.     cvSetData(Vector1,Coordinates,Vector1->step);  
28.     AvgVector=cvCreateMat(1,2,CV_32FC1);  
29.     EigenValue_Row=cvCreateMat(2,1,CV_32FC1);  
30.     EigenVector=cvCreateMat(2,2,CV_32FC1);  
31.   
32.     cvCalcPCA(Vector1,AvgVector,EigenValue_Row,EigenVector,CV_PCA_DATA_AS_ROW);  
33.   
34.     printf("Original Data:\n");  
35.     PrintMatrix(Vector1,10,2);  
36.   
37.     printf("==========\n");  
38.     PrintMatrix(AvgVector,1,2);  
39.   
40.     printf("\nEigne Value:\n");  
41.     PrintMatrix(EigenValue_Row,2,1);  
42.   
43.     printf("\nEigne Vector:\n");  
44.     PrintMatrix(EigenVector,2,2);  
45.   
46.     system("pause");  
47. }  
48. void PrintMatrix(CvMat *Matrix,int Rows,int Cols)  
49. {  
50.     for(int i=0;i<Rows;i++)  
51.     {  
52.         for(int j=0;j<Cols;j++)  
53.         {  
54.             printf("%.2f ",cvGet2D(Matrix,i,j).val[0]);  
55.         }  
56.         printf("\n");  
57.     }  
58. }  

执行结果:


这部份是把平均数,共变量矩阵,以及特征值及特征向量都计算出来了,全部都包在cvCalcPCA()的函式里面,因此可以不必特地的用cvCalcCovarMatrix()求得共变量矩阵,也不需要再由共变量矩阵套用cvEigenVV()求得它的EigenValue以及EigenVector了,所以说,cvCalcPCA()=cvCalcCovarMatrix()+cvEigenVV(),不仅如此,cvCalcPCA()使用上更是灵活,当向量的维度数目比输入的数据那的时候(例如Eigenface),它的共变量矩阵就会自动转成CV_COVAR_SCRAMBLED,而当输入数据量比向量维度大的时候,它亦会自动转成CV_COVAR_NORMAL的形态,而OpenCV也提供了计算投影量cvProjectPCA(),以及反向投影的函式cvBackProjectPCA(),cvCalcPCA()的计算结果如下


详细计算方法可以参考"OpenCV统计应用-共变量矩阵"以及"OpenCV线性代数-cvEigenVV实作"这两篇,cvCalcPCA()第一个自变量为输入目标要计算的向量,整合在CvMat数据结构里,第二个自变量为空的平均数向量,第三个自变量为输出排序后的EigenValue,以列(Rows)为主的数值,第四个自变量为排序后的EigenVector,第五个自变量为cvCalcPCA()的参数,它的参数公式如下

#define CV_PCA_DATA_AS_ROW 0
#define CV_PCA_DATA_AS_COL 1
#define CV_PCA_USE_AVG 2

分别代表以列为主,以栏为主的参数设定,以及使用自己定义的平均数,CV_PCA_DATA_AS_ROW与CV_PCA_DATA_AS_COL的参数不可以同时使用,而对于主成分分析的EigenVector对原始数据投影的程序范例如下

PCA特征向量投影

[cpp] view plaincopy

1. #include <cv.h>  
2. #include <highgui.h>  
3. #include <stdio.h>  
4. #include <stdlib.h>  
5.   
6. float Coordinates[20]={1.5,2.3,  
7.                                  3.0,1.7,  
8.                                  1.2,2.9,  
9.                                  2.1,2.2,  
10.                                  3.1,3.1,  
11.                                  1.3,2.7,  
12.                                  2.0,1.7,  
13.                                  1.0,2.0,  
14.                                  0.5,0.6,  
15.                                  1.0,0.9};  
16.   
17. void PrintMatrix(CvMat *Matrix,int Rows,int Cols);  
18.   
19. int main()  
20. {  
21.     CvMat *Vector1;  
22.     CvMat *AvgVector;  
23.     CvMat *EigenValue_Row;  
24.     CvMat *EigenVector;  
25.   
26.     Vector1=cvCreateMat(10,2,CV_32FC1);  
27.     cvSetData(Vector1,Coordinates,Vector1->step);  
28.     AvgVector=cvCreateMat(1,2,CV_32FC1);  
29.     EigenValue_Row=cvCreateMat(2,1,CV_32FC1);  
30.     EigenVector=cvCreateMat(2,2,CV_32FC1);  
31.   
32.     cvCalcPCA(Vector1,AvgVector,EigenValue_Row,EigenVector,CV_PCA_DATA_AS_ROW);  
33.     cvProjectPCA(Vector1,AvgVector,EigenVector,Vector1);  
34.   
35.     printf("Project Original Data:\n");  
36.     PrintMatrix(Vector1,10,2);  
37.   
38.     system("pause");  
39. }  
40. void PrintMatrix(CvMat *Matrix,int Rows,int Cols)  
41. {  
42.     for(int i=0;i<Rows;i++)  
43.     {  
44.         for(int j=0;j<Cols;j++)  
45.         {  
46.             printf("%.2f ",cvGet2D(Matrix,i,j).val[0]);  
47.         }  
48.         printf("\n");  
49.     }  
50. }  

执行结果:


cvProjectPCA(),它的公式定义如下


因此,它所呈现的结果就是将坐标旋转后的特征向量投影,cvProjectPCA()的地一个自变量为输入原始向量数据,第二个自变量为输入空的(或是以设定好的)平均数向量,第三个自变量为输入以排序好的特征向量EigenVector,第四个自变量为输出目标投影向量,而投影之外,OpenCV还可以给他做反向投影回原始的数据

PCA反向投影

[cpp] view plaincopy

1. #include <cv.h>  
2. #include <highgui.h>  
3. #include <stdio.h>  
4. #include <stdlib.h>  
5.   
6. float Coordinates[20]={1.5,2.3,  
7.                                  3.0,1.7,  
8.                                  1.2,2.9,  
9.                                  2.1,2.2,  
10.                                  3.1,3.1,  
11.                                  1.3,2.7,  
12.                                  2.0,1.7,  
13.                                  1.0,2.0,  
14.                                  0.5,0.6,  
15.                                  1.0,0.9};  
16.   
17. int main()  
18. {  
19.     CvMat *Vector1;  
20.     CvMat *AvgVector;  
21.     IplImage *Image1=cvCreateImage(cvSize(450,450),IPL_DEPTH_8U,3);  
22.     Image1->origin=1;  
23.   
24.     Vector1=cvCreateMat(10,2,CV_32FC1);  
25.     cvSetData(Vector1,Coordinates,Vector1->step);  
26.   
27.     AvgVector=cvCreateMat(1,2,CV_32FC1);  
28.   
29.     CvMat *EigenValue_Row=cvCreateMat(2,1,CV_32FC1);  
30.     CvMat *EigenVector=cvCreateMat(2,2,CV_32FC1);  
31.   
32.     cvCalcPCA(Vector1,AvgVector,EigenValue_Row,EigenVector,CV_PCA_DATA_AS_ROW);  
33.     cvProjectPCA(Vector1,AvgVector,EigenVector,Vector1);  
34.     cvBackProjectPCA(Vector1,AvgVector,EigenVector,Vector1);  
35.   
36.     printf("Back Project Original Data:\n");  
37.     for(int i=0;i<10;i++)  
38.     {  
39.         printf("%.2f ",cvGetReal2D(Vector1,i,0));  
40.        printf("%.2f ",cvGetReal2D(Vector1,i,1));  
41.   
42.         cvCircle(Image1,cvPoint((int)(cvGetReal2D(Vector1,i,0)*100),(int)(cvGetReal2D(Vector1,i,1)*100)),0,CV_RGB(0,0,255),10,CV_AA,0);  
43.   
44.         printf("\n");  
45.     }  
46.     cvCircle(Image1,cvPoint((int)((cvGetReal2D(AvgVector,0,0))*100),(int)((cvGetReal2D(AvgVector,0,1))*100)),0,CV_RGB(255,0,0),10,CV_AA,0);  
47.   
48.     printf("==========\n");  
49.     printf("%.2f ",cvGetReal2D(AvgVector,0,0));  
50.     printf("%.2f ",cvGetReal2D(AvgVector,0,1));  
51.   
52.     cvNamedWindow("Coordinates",1);  
53.     cvShowImage("Coordinates",Image1);  
54.     cvWaitKey(0);  
55. }  

执行结果:


而反向投影,只不过是将投影向量还原成原始数据罢了,可以用来做为新进数据反向投影后用来比对的步骤,以下是它的计算公式推导


cvBackProjectPCA()的自变量输入与cvProjectPCA()是一样的,只不过是里面的计算公式不同,其实也只是cvProjectPCA()的反运算

cvCalcPCA()
计算多笔高维度数据的主要特征值及特征向量,先自动判断维度大于数据量或维度小于数据量来选择共变量矩阵的样式,在由共变量矩阵求得已排序后的特征值及特征向量,cvCalcPCA()函式的参数为CV_PCA_DATA_AS_ROW以列为主的资料排列,CV_PCA_DATA_AS_COL以行为主的数据排列,CV_PCA_USE_AVG自行定义平均数数据,CV_PCA_DATA_AS_ROW以及CV_PCA_DATA_AS_COL不可以同时使用,cvCalcPCA()的第一个自变量为输入CvMat维度向量数据结构,第二个自变量为输入空的CvMat平均数向量数据结构(或是自行定义平均数向量),第三个自变量为输入以列为主的特征值CvMat数据结构,第四个自变量为输入特征向量CvMat数据结构,第四个自变量为输入cvCalcPCA()函式的参数
cvCalcPCA(输入目标向量数据CvMat数据结构,输入或输出向量平均数CvMat数据结构,输入以列为主EigenValue的CvMat数据结构,输入EigenVector的CvMat数据结构,目标参数或代号)

cvProjectPCA()
将计算主成分分析的结果做投影的运算,主要作用是将最具意义的k个EgienVector与位移后的原始资料做矩阵乘法,投影过后的结果将会降低维度,cvProjectPCA()第一个自变量为输入CvMat数据结构原始向量数据数据,第二个自变量为输入CvMat数据结构平均数向量,第三个自变量为输入要降维的特征向量,第四个自变量为输出CvMat数据结构原始数据投影的结果
cvProjectPCA(输入CvMat原始向量数据数据结构,输入CvMat原始向量平均数数据结构,输入CvMat降低维度的特征向量数据结构,输出CvMat投影向量数据结构)

cvBackProjectPCA()
将投影向量转回原始向量维度坐标系,cvProjectPCA()第一个自变量为输入CvMat数据结构投影向量数据,第二个自变量为输入CvMat数据结构平均数向量,第三个自变量为输入特征向量,第四个自变量为输出CvMat数据结构反向投影的结果
cvProjectPCA(输入CvMat投影向量数据结构,输入CvMat原始向量平均数数据结构,输入CvMat特征向量数据结构,输出CvMat反向投影数据结构)

另附：http://blog.csdn.net/yang_xian521/article/details/7445536