OpenCV統計應用-Mahalanobis距離

Mahalanobis距離是一個可以準確找出資料分布上面極端值(Outliers)的統計方法,使用線性迴歸的概念,也就是說他使用的是共變數矩陣以及該資料分布的平均數來找尋極端值的產生,而可以讓一群資料系統具有穩健性(Robust),去除不必要的雜訊訊息,這邊拿前面共變數矩陣的資料為例,並且新增了兩個點座標向量來做Mahalanobis距離的比較



加入兩個座標點找尋極端值


Mahalanobis距離實作
#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};

float Coordinates2[2]={1.3,1.5};
float Coordinates3[2]={200,100};

int main()
{
    CvMat *Vector[1];
    CvMat *Vector1;
    CvMat *CovarMatrix;
    CvMat *InvertCovarMatrix;
    CvMat *AvgVector;
    CvMat *Vector2;
    CvMat *Vector3;

    Vector1=cvCreateMat(10,2,CV_32FC1);
    cvSetData(Vector1,Coordinates1,Vector1->step);
    Vector[0]=Vector1;
    CovarMatrix=cvCreateMat(2,2,CV_32FC1);
    InvertCovarMatrix=cvCreateMat(2,2,CV_32FC1);

    AvgVector=cvCreateMat(1,2,CV_32FC1);
    Vector2=cvCreateMat(1,2,CV_32FC1);
    Vector3=cvCreateMat(1,2,CV_32FC1);
    cvSetData(Vector2,Coordinates2,Vector2->step);
    cvSetData(Vector3,Coordinates3,Vector3->step);

    cvCalcCovarMatrix((const CvArr **)Vector,10,CovarMatrix,AvgVector,CV_COVAR_SCALE+CV_COVAR_NORMAL+CV_COVAR_ROWS);
    cvInvert(CovarMatrix,InvertCovarMatrix,CV_SVD_SYM);

    printf("\nVector2 Mahalanobis Distance\n");
    printf("%f\n",cvMahalanobis(Vector2,AvgVector,InvertCovarMatrix));

    printf("\nVector3 Mahalanobis Distance\n");
    printf("%f\n",cvMahalanobis(Vector3,AvgVector,InvertCovarMatrix));

    system("pause");
}

由上面可以看得出來第一筆座標向量Vector2,他的座標位在(1.3,1.5)的位置,第二筆座標向量Vector2他的座標在(200,100),對於原始資料分布來說,第二筆座標向量它已經遠遠離開了這些座標分佈,因此可以從執行結果看的出來Vector3的Mahalanobis距離遠大於Vector2,而在做Mahalanobis距離的同時,也必須要將共變數矩陣做反矩陣運算,cvMahalanobis()第一個引數為目標要判斷是否是極端值的向量CvMat資料結構,第二個引數為該資料分佈的平均數向量或是該資料分佈的某筆資料的CvMat資料結構,第三個引數為該資料分佈共變數矩陣的反矩陣

而Mahalanobis距離的原理如下


當資料分佈使的共變數矩陣為單位向量I的話,就會退化成歐機理得距離(Euclidean Distance)


而他的計算方式如下



由上面可知,在一個簡單線性迴歸的資料模型下,第二筆資料Vector3為這個線性迴歸的極端值

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