2009年1月31日 星期六

OpenCV統計應用-PCA主成分分析

在圖形識別方面,主成分分析(Principal Comonents Analysis,PCA)算是比較快速而且又準確的方式之一,它可以對抗圖形平移旋轉的事件發生,並且藉由主要特徵(主成分)投影過後的資料做資料的比對,在多個特徵資訊裡面,取最主要的K個,做為它的特徵依據,在這邊拿前面共變數矩陣的數據來做沿用,主成分分析使用的方法為計算共變數矩陣,在加上計算共變數矩陣的特徵值及特徵向量,將特徵值以及所對應的特徵向量排序之後,取前面主要K個特徵向量當做主要特徵,而OpenCV也可以對高維度的向量進行主成分分析的計算


數據原始的分佈情況,紅點代表著它的平均數


將座標系位移,以紅點為主要的原點


計算共變數矩陣以及共變數的特徵值以及特徵向量,將特徵向量排序後投影回原始資料的結果的結果,也就是說,對照上面的圖片,EigenVector的作用是找到主軸後,將原本的座標系做旋轉了


再來就是對它做投影,也就是降低維度的動作,將Y軸的數據全部歸零,投影在X軸上


投影完之後,在將它轉回原本的座標系


PCA主成分分析,與線性迴歸有異曲同工之妙,也就是說,這條投影過後的直線,可以稱做迴歸線,當它在做主軸旋轉的時候,所投影的結果為最小均方誤,在將它轉置回來的時候,就形成了一條迴歸直線了

OpenCV的PCA輸入必須要是單通道32位元浮點數格式或是單通道64位元浮點數格式的,參數為CV_32FC1或是CV_64FC1,程式寫法如下

PCA程式實作
#include <cv.h>
#include <highgui.h>
#include <stdio.h>
#include <stdlib.h>


float Coordinates[20]={1.5,2.3,
                                 3.0,1.7,
                                 1.2,2.9,
                                 2.1,2.2,
                                 3.1,3.1,
                                 1.3,2.7,
                                 2.0,1.7,
                                 1.0,2.0,
                                 0.5,0.6,
                                 1.0,0.9};

void PrintMatrix(CvMat *Matrix,int Rows,int Cols);

int main()
{
    CvMat *Vector1;
    CvMat *AvgVector;
    CvMat *EigenValue_Row;
    CvMat *EigenVector;

    Vector1=cvCreateMat(10,2,CV_32FC1);
    cvSetData(Vector1,Coordinates,Vector1->step);
    AvgVector=cvCreateMat(1,2,CV_32FC1);
    EigenValue_Row=cvCreateMat(2,1,CV_32FC1);
    EigenVector=cvCreateMat(2,2,CV_32FC1);

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

    printf("Original Data:\n");
    PrintMatrix(Vector1,10,2);

    printf("==========\n");
    PrintMatrix(AvgVector,1,2);

    printf("\nEigne Value:\n");
    PrintMatrix(EigenValue_Row,2,1);

    printf("\nEigne Vector:\n");
    PrintMatrix(EigenVector,2,2);

    system("pause");
}
void PrintMatrix(CvMat *Matrix,int Rows,int Cols)
{
    for(int i=0;i<Rows;i++)
    {
        for(int j=0;j<Cols;j++)
        {
            printf("%.2f ",cvGet2D(Matrix,i,j).val[0]);
        }
        printf("\n");
    }
}

執行結果:


這部份是把平均數,共變數矩陣,以及特徵值及特徵向量都計算出來了,全部都包在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特徵向量投影
#include <cv.h>
#include <highgui.h>
#include <stdio.h>
#include <stdlib.h>


float Coordinates[20]={1.5,2.3,
                                 3.0,1.7,
                                 1.2,2.9,
                                 2.1,2.2,
                                 3.1,3.1,
                                 1.3,2.7,
                                 2.0,1.7,
                                 1.0,2.0,
                                 0.5,0.6,
                                 1.0,0.9};

void PrintMatrix(CvMat *Matrix,int Rows,int Cols);

int main()
{
    CvMat *Vector1;
    CvMat *AvgVector;
    CvMat *EigenValue_Row;
    CvMat *EigenVector;

    Vector1=cvCreateMat(10,2,CV_32FC1);
    cvSetData(Vector1,Coordinates,Vector1->step);
    AvgVector=cvCreateMat(1,2,CV_32FC1);
    EigenValue_Row=cvCreateMat(2,1,CV_32FC1);
    EigenVector=cvCreateMat(2,2,CV_32FC1);

    cvCalcPCA(Vector1,AvgVector,EigenValue_Row,EigenVector,CV_PCA_DATA_AS_ROW);
    cvProjectPCA(Vector1,AvgVector,EigenVector,Vector1);

    printf("Project Original Data:\n");
    PrintMatrix(Vector1,10,2);

    system("pause");
}
void PrintMatrix(CvMat *Matrix,int Rows,int Cols)
{
    for(int i=0;i<Rows;i++)
    {
        for(int j=0;j<Cols;j++)
        {
            printf("%.2f ",cvGet2D(Matrix,i,j).val[0]);
        }
        printf("\n");
    }
}

執行結果:


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


因此,它所呈現的結果就是將座標旋轉後的特徵向量投影,cvProjectPCA()的地一個引數為輸入原始向量資料,第二個引數為輸入空的(或是以設定好的)平均數向量,第三個引數為輸入以排序好的特徵向量EigenVector,第四個引數為輸出目標投影向量,而投影之外,OpenCV還可以給他做反向投影回原始的資料

PCA反向投影
#include <cv.h>
#include <highgui.h>
#include <stdio.h>
#include <stdlib.h>


float Coordinates[20]={1.5,2.3,
                                 3.0,1.7,
                                 1.2,2.9,
                                 2.1,2.2,
                                 3.1,3.1,
                                 1.3,2.7,
                                 2.0,1.7,
                                 1.0,2.0,
                                 0.5,0.6,
                                 1.0,0.9};

int main()
{
    CvMat *Vector1;
    CvMat *AvgVector;
    IplImage *Image1=cvCreateImage(cvSize(450,450),IPL_DEPTH_8U,3);
    Image1->origin=1;

    Vector1=cvCreateMat(10,2,CV_32FC1);
    cvSetData(Vector1,Coordinates,Vector1->step);

    AvgVector=cvCreateMat(1,2,CV_32FC1);

    CvMat *EigenValue_Row=cvCreateMat(2,1,CV_32FC1);
    CvMat *EigenVector=cvCreateMat(2,2,CV_32FC1);

    cvCalcPCA(Vector1,AvgVector,EigenValue_Row,EigenVector,CV_PCA_DATA_AS_ROW);
    cvProjectPCA(Vector1,AvgVector,EigenVector,Vector1);
    cvBackProjectPCA(Vector1,AvgVector,EigenVector,Vector1);

    printf("Back Project Original Data:\n");
    for(int i=0;i<10;i++)
    {
        printf("%.2f ",cvGetReal2D(Vector1,i,0));
       printf("%.2f ",cvGetReal2D(Vector1,i,1));

        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);

        printf("\n");
    }
    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);

    printf("==========\n");
    printf("%.2f ",cvGetReal2D(AvgVector,0,0));
    printf("%.2f ",cvGetReal2D(AvgVector,0,1));

    cvNamedWindow("Coordinates",1);
    cvShowImage("Coordinates",Image1);
    cvWaitKey(0);
}

執行結果:


而反向投影,只不過是將投影向量還原成原始資料罷了,可以用來做為新進資料反向投影後用來比對的步驟,以下是它的計算公式推導


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反向投影資料結構)



10 意見:

匿名 提到...

你好,作者计算的pca都是数据,维数都比较小的。但是在实际当中我们处理的维数都比较大,我处理了两张60*60的图像,可pca运行了20多分钟还没结束,想问下这是什么原因?谢谢!

yester 提到...

如果是只有兩張圖像
cvCalcPCA()會自動轉換成
CV_COVAR_SCRAMBLED
而且共變數矩陣會是2*2的大小
我看過Source code轉換的程式應該是沒有問題
如果我有時間在試看看大維度的...

匿名 提到...

你好~我是個新手我想請問一下...小弟最近再用opencv做PCA...實在搞不懂createMat跟createImage差在哪裡?!使得我沒辦法對一張圖做PCA的train...大大的輸入是點
...可以教教我怎麼弄成圖片嗎??

匿名 提到...

最近使用openCV來做我的FYP。
我的FYP是Face Recognition!請問一下,如果用opencv做eigenfaces呢?另外,如果我要做dimension reduction(PCA的強項),是否選K個eigenvalue,然后,如何應用在程式里?

匿名 提到...

y大你好:

我目前是統計所的學生,有空時也會寫一些

使用C語言撰寫統計程式的文章。

我認為y大你寫的很仔細,不只統計方面解說

程式方面也寫的很簡潔。

請你繼續更新下去喔!

匿名 提到...

问个问题:
用PCA做图像融合,求协方差矩阵的特征向量时,用OpenCV和其他工具(matlab)求出的特征向量符号相反,比如:u=[0.33,0.233,-0.5]和 u=[-0.33,-0.233,0.5];这对于投影(Projection)后得到的图像的显示结果差别很大,请问:特征向量的正负符号有什么选取标准,在做图像相关的PCA变换时正负号的物理意义是什么呢?
hk-teaching@163.com

匿名 提到...

想請問一下
你那些程式 執行之後
都沒有 最上面的
那些圖片
只有下面的程式圖片而已

想請問一下 PCA 做出來之後
應該是類似 上面那些分群的圖案吧
可是下面的都沒有
是沒有提供嗎

匿名 提到...

你好,我最近在使用cvCalcPCA,
可是每當data的維度大於data的數目時(即data的cols>rows),
就會出現runtime的exception,
我是用CV_PCA_DATA_AS_ROW,
想請問y大有出現過這個問題嗎?
謝謝

Unknown 提到...

想請問大大, PCA 可以輸出投影的角度嗎?

匿名 提到...

大大你好!我想請問一下,計算完PCA或特徵值以後,正負號所代表的意思,謝謝
jay82217@yahoo.com.tw

Copyright 2008-2009,yester