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



2009年1月23日 星期五

OpenCV統計應用-共變數矩陣

共變數(Covariance),為兩個隨機變數的離均差除以母體個數,可以判斷兩個事件隨機變數的相依性,而單一變數的共變數,那就是變異數(Variance)了,共變數以及變異數簡單的定義如下

變異數


共變數


而當兩隨機變數為獨立事件的時候


再來提到的是共變數矩陣(Covariance matrix),代表著隨機變數所有共變數的種類,以兩組隨機變數來說所產生的共變數矩陣如下


而三組,四組以上隨機變數的共變數矩陣又是不同情形了,在這邊,隨機變數的個數代表著向量的維度,而他的數對則是同時發生的情形,下面就以座標的簡單例子來示範


這是個維度為2的座標數組,代表的是一個座標平面的點集合,而下面就是給它跑cvCalcCovarMatrix()共變數矩陣的計算方法

共變數矩陣1
#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 *Vector[10];
    CvMat *CovarMatrix;
    CvMat *AvgMatrix;
    IplImage *Image1=cvCreateImage(cvSize(450,450),IPL_DEPTH_8U,3);
    Image1->origin=1;
    for(int i=0;i<10;i++)
    {
        Vector[i]=cvCreateMat(1,2,CV_32FC1);
        cvSetReal1D(Vector[i],0,Coordinates[i*2]);
        cvSetReal1D(Vector[i],1,Coordinates[i*2+1]);
        cvCircle(Image1,cvPoint((int)(Coordinates[i*2]*100),(int)(Coordinates[i*2+1]*100)),0,CV_RGB(0,0,255),10,CV_AA,0);
    }

    CovarMatrix=cvCreateMat(2,2,CV_32FC1);
    AvgMatrix=cvCreateMat(1,2,CV_32FC1);
    cvCalcCovarMatrix((const CvArr **)Vector,10,CovarMatrix,AvgMatrix,CV_COVAR_SCALE+CV_COVAR_NORMAL);

    for(int i=0;i<2;i++)
    {
        for(int j=0;j<2;j++)
        {
            printf("%f ",cvGetReal2D(CovarMatrix,i,j));
        }
        printf("\n");
    }
    cvNamedWindow("Coordinates",1);
    cvShowImage("Coordinates",Image1);
    cvWaitKey(0);
}

執行結果:


計算方法:


上面的方法是用個別的向量來實作,向量的維度為2,而圖片顯示的結果為它們二維座標的分佈情況,cvCalcCovarMatrix()必須給它空的平均值向量矩陣來計算,而cvCalcCovarMatrix()它的參數被定義如下

#define CV_COVAR_SCRAMBLED 0
#define CV_COVAR_NORMAL 1
#define CV_COVAR_USE_AVG 2
#define CV_COVAR_SCALE 4
#define CV_COVAR_ROWS 8
#define CV_COVAR_COLS 16

由上面可以知道,它是由2的N次方所表達的,所以可以用合成參數的方式來表達cvCalcCovarMatrix()的共變數矩陣函式,也就是說,可以用CV_COVAR_SCRAMBLED+CV_COVAR_ROWS的組合,也可以用CV_COVAR_NORMAL+CV_COVAR_USE_AVG+CV_COVAR_ROWS之類的組合,而它所代表的含意分別是

CV_COVAR_SCRAMBLED
一種共變數矩陣的計算方式,不可以與CV_COVAR_NORMAL合用,表達方式如下

這計算方式在OpenCV說明文件提到為Eigenface的計算方式

CV_COVAR_NORMAL
共變數矩陣計算方式,不可與CV_COVAR_SCRAMBLED合用,表達方式如下

這個則是為一般共變數矩陣的計算方式

CV_COVAR_USE_AVG
不用共變數矩陣內建計算平均數的函式,而是自己給予平均數值,可與任何參數共用

CV_COVAR_SCALE
對共變數矩陣的數據做向量個數總和的純量積,用的是除法計算如下的共變數矩陣

可與任何參數共用,而如果沒這參數則是無除以N的計算

CV_COVAR_ROWS
將共變數矩陣的輸入值用矩陣的方式表達,而不是用向量的方式,矩陣表達方式以列(Rows)為主,並且參數不可與CV_COVAR_COLS合用

CV_COVAR_COLS
將共變數矩陣的輸入值用矩陣的方式表達,而不是用向量的方式,矩陣表達方式以欄(Columns)為主,並且參數不可與CV_COVAR_ROWS共用

共變數矩陣有幾個規則,也就是,輸入一定要是方陣,平均數的長度要是向量的維度,而平均數的長度也一定要是向量的大小,然後一定要是用單通道CV_32FC1或是CV_64FC1做為輸入,cvCalcCovarMatrix()函式,第一個引數則必須要用(const CvArr **)強制型別轉換,第二個引數為輸入向量的數目,第三個引數為空的或是非空平均數向量,第四個引數為cvCalcCovarMatrix()這函式要輸入的參數,而下面,則是使用矩陣方式表達共變數矩陣的範例

共變數矩陣2
#include <cv.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 *Vector[1];
    CvMat *Vector1;
    CvMat *CovarMatrix;
    CvMat *avg;

    Vector1=cvCreateMat(10,2,CV_32FC1);
    cvSetData(Vector1,Coordinates,Vector1->step);
    Vector[0]=Vector1;
    CovarMatrix=cvCreateMat(2,2,CV_32FC1);
    avg=cvCreateMat(1,2,CV_32FC1);

    cvCalcCovarMatrix((const CvArr **)Vector,10,CovarMatrix,avg,CV_COVAR_SCALE+CV_COVAR_NORMAL+CV_COVAR_ROWS);

    for(int i=0;i<2;i++)
    {
        for(int j=0;j<2;j++)
        {
            printf("%f ",cvGetReal2D(CovarMatrix,i,j));
        }
        printf("\n");
    }
    system("pause");
}

執行結果:


上面的程式碼,如果想用矩陣來表達向量的方式計算共變數矩陣,就必須要將第一個引數設為二維陣列當做輸入,輸入的方式就如上面程式的寫法,而其他地方則是沒什麼差異.

共變數矩陣在OpenCV內,不但可以做到計算主成分分析(Principal Cmponents Analysis,PCA),以及另一個Mahalanobis距離的計算

cvCalcCovarMatrix()
計算共變數矩陣,輸入可為多個IplImage或CvMat資料結構,可輸入高維度資料的向量,有多種輸入方式,在資料輸入方面,可以用單一CvMat資料結構,以列(Rows)為主或以行(Columns)為主,使用CV_COVAR_ROWS以及CV_COVAR_COLS的參數輸入,而要做多個IplImage或CvMat資料結構輸入則不需提供CV_COVAR_ROWS,CV_COVAR_COLS的參數,在計算方面,又分為以維度為主的共變數矩陣參數輸入CV_COVAR_NORMAL以及以個數為主的共變數矩陣CV_COVAR_SCRAMBLED,而他也可一自行定義平均值CV_COVAR_USE_AVG,以及除以是否除以共變數矩陣的共變數個數CV_COVAR_SCALE,而這些參數可以自行組合運算,第一個引數為輸入目標向量,第二個引數為輸出目標共變數矩陣,第三個引數為輸入或輸出平均數向量,第四個引數為cvCalcCovarMatrix()共變數計算函式的參數輸入
cvCalcCovarMatrix(輸入多個IplImage或CvMat資料結構,輸出目標共變數矩陣,輸入/出平均數向量,共變數矩陣參數或代號)



2009年1月2日 星期五

OpenCV統計應用-直方圖比較

cvCompareHist(),是比較兩個統計直方圖的分布,總共有四個方法,被定義如下:

#define CV_COMP_CORREL 0
#define CV_COMP_CHISQR 1
#define CV_COMP_INTERSECT 2
#define CV_COMP_BHATTACHARYYA 3

而這些方法分別為相關係數,卡方,交集法以及在做常態分布比對的Bhattacharyya距離,這些方法都是用來做統計直方圖的相似度比較的方法,而且,都是根據統計學的概念,這邊就簡單的拿來用灰階統計直方圖來比較,而這部份的比較方式,是由圖形的色彩結構來著手,下面就簡單的用三種情況來分析它們距離比較的方式

直方圖比較實作
#include <cv.h>
#include <highgui.h>
#include <stdio.h>
#include <stdlib.h>

int HistogramBins = 256;
float HistogramRange1[2]={0,255};
float *HistogramRange[1]={&HistogramRange1[0]};

int main()
{
    IplImage *Image1=cvLoadImage("RiverBank.jpg",0);
    IplImage *Image2=cvLoadImage("DarkClouds.jpg",0);

    CvHistogram *Histogram1=cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,HistogramRange);
    CvHistogram *Histogram2=cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,HistogramRange);

    cvCalcHist(&Image1,Histogram1);
    cvCalcHist(&Image2,Histogram2);

    cvNormalizeHist(Histogram1,1);
    cvNormalizeHist(Histogram2,1);

    printf("CV_COMP_CORREL : %.4f\n",cvCompareHist(Histogram1,Histogram2,CV_COMP_CORREL));
    printf("CV_COMP_CHISQR : %.4f\n",cvCompareHist(Histogram1,Histogram2,CV_COMP_CHISQR));
    printf("CV_COMP_INTERSECT : %.4f\n",cvCompareHist(Histogram1,Histogram2,CV_COMP_INTERSECT));
    printf("CV_COMP_BHATTACHARYYA : %.4f\n",cvCompareHist(Histogram1,Histogram2,CV_COMP_BHATTACHARYYA));

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

原始圖片:





 

執行結果:

(1)RiverBank.jpg & DarkClouds.jpg

Output:
CV_COMP_CORREL : -0.1407
CV_COMP_CHISQR : 0.6690
CV_COMP_INTERSECT : 0.4757
CV_COMP_BHATTACHARYYA : 0.4490

(2)RiverBank.jpg & RiverBank.jpg

Output:
CV_COMP_CORREL : 1
CV_COMP_CHISQR : 0
CV_COMP_INTERSECT : 1
CV_COMP_BHATTACHARYYA : 0

(3)Black.jpg & White.jpg

Output:
CV_COMP_CORREL : 1
CV_COMP_CHISQR : 1
CV_COMP_INTERSECT : 0
CV_COMP_BHATTACHARYYA : 1

這邊的直方圖比較,則是將它們的統計直方圖用cvNormalizeHist()正規化成1,在由正規化的統計分佈來做直方圖的比較,從上面的輸出結果可以推測出,卡方法以及Bhattacharyya是數值越小圖形越相似,而相關係數則是看圖形的分佈程度,因此第三個Black.jpg&White.jpg所顯示相關係數的結果才會是1,這邊的圖形比較用的是統計學的方法,而一般的比較方式還有歐幾里德距離的方式,在這個cvCompareHist()的函式,也可以實作出多通道的多維度直方圖比較,而且支援CV_HIST_ARRAY及CV_HIST_SPARSE這兩種直方圖資料結構的格式,而這些比對方法的相關公式如下


相關係數法

而它的公式推導如下



再來是卡方的方式

這邊的卡方法跟一般的卡方檢定不太一樣,下面是一般適合度檢定(Goodness of fit test)的公式

o為觀察者次數,e為期望值次數


交集的方式就比較簡單了,兩個直方圖取最小的做累加



再來就是常態分配比對的Bhattacharyya距離






Copyright 2008-2009,yester