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

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資料結構,輸出目標共變數矩陣,輸入/出平均數向量,共變數矩陣參數或代號)

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距離

OpenCV統計應用-直方圖反向投影

影像處理的統計直方圖,可以知道一張圖片在該色彩空間的數據分布狀況,而這邊,就要介紹到直方圖反向投影的函式,直方圖反向投影,也就是將數據分布的狀況依照Look-up table的方式對應回去,實際上,這個函式是跟前面介紹到的cvLUT()是一樣的,只不過,差別是差異在cvLUT()的第三個引數改變成CvHistogram資料結構的輸入,直方圖反向投影,cvCalcBackProject()的第一個引數是輸入單通道IplImage資料結構,第二個引數是輸出單通道IplImage反向投影圖形資料結構,第三個引數是選定要被反向投影的CvHistogram直方圖資料結構,而cvCalcBackProject()把前面提到的Look-up table的計算方式包在cvCalcBackProject()函式的底層,因此,它可以整合CvHistogram這個資料結構做更多的應用,下面這個就是修改前面的範例"OpenCV統計應用-CvHistogram直方圖資料結構",來做直方圖反向投影的程式

灰階直方圖反向投影
#include <cv.h>
#include <highgui.h>
#include <stdio.h>

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

CvPoint Point1;
CvPoint Point2;

int main()
{
    IplImage *GrayImage1;
    IplImage *Image1;
    IplImage *Image2;
    IplImage *BackProjectImage;
    CvHistogram *Histogram1;
    IplImage *HistogramImage1;

    Image1=cvLoadImage("Riverbank.jpg",1);
    Image2=cvCreateImage(cvGetSize(Image1),IPL_DEPTH_8U,3);
    GrayImage1=cvLoadImage("Riverbank.jpg",0);
    BackProjectImage=cvCreateImage(cvGetSize(Image1),IPL_DEPTH_8U,1);
    Histogram1 = cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,HistogramRange);
    HistogramImage1 = cvCreateImage(cvSize(256,300),8,3);

    cvSetZero(HistogramImage1);
    HistogramImage1->origin=1;
    HistogramBinWidth=256/HistogramBins;

    cvCalcHist(&GrayImage1,Histogram1);
    cvNormalizeHist(Histogram1,5000);

    cvThreshHist(Histogram1,50);
    cvCalcBackProject(&GrayImage1,BackProjectImage,Histogram1);
    cvCopy(Image1,Image2,BackProjectImage);

    for(int i=0;i<HistogramBins;i++)
    {
        printf("%f\n",cvQueryHistValue_1D(Histogram1,i));
        Point1=cvPoint(i*HistogramBinWidth,0);
        Point2=cvPoint((i+1)*HistogramBinWidth,(int)        cvQueryHistValue_1D(Histogram1,i));

        cvRectangle(HistogramImage1,Point1,Point2,CV_RGB(127,127,127));
    }

    cvNamedWindow("Histogram1",1);
    cvNamedWindow("Riverbank",1);
    cvNamedWindow("Back Project RiverBank",1);
    cvShowImage("Riverbank",Image1);
    cvShowImage("Back Project RiverBank",Image2);
    cvShowImage("Histogram1",HistogramImage1);
    cvWaitKey(0);
}

執行結果:


這邊就是拿前面灰階去除較小直方圖區塊的程式碼做修改,然後將前面去除最小區塊的部份,對應到彩色的圖片去了,顯示的結果會是,只要是影像裡面灰階值分佈數量比較少的數據,全部都被對應成黑色的像素,也就是全部都變成0,而cvCalcBackProject(),就跟cvLUT()一樣直接拿直方圖的數據去對應,假設一個直方圖從頭開始的數據為254,129,80,70....那麼只要是圖片內像素值為1的數據就會對應到254,1的圖片裡面像素值的數據就會直接變成254,像素值為2的就會對應到129,2的數據就會直接變成129,以此類推,因此,用cvThreshHist()去除小於50的直方圖區段,讓小於50的全部歸0,再來,對他做一個圖片的反向投影,所對應出來的結果雖然不是0或255,可是它卻可以直接拿來當做是遮罩,也就是說,直接拿來給cvCopy()來做對應,因此,反向投影的結果就出來啦,在遮罩的部份就要參考"資料結構操作與運算-圖形的Mask遮罩實作",而這段程式碼則是用到"OpenCV統計應用-CvHistogram資料結構操作"裡面的部份.

在OpenCV Documentation的部份有提到cvCalcBackProject()可以對HSV色彩空間的Hue做反向投影,到底是怎麼實作出來呢？在OpenCV的Sample Code裡面有一個camshift.c的程式,就是用到這個反向投影的函式,而它的反向投影的簡單範例,原理就如下所示

HSV色彩空間反向投影
#include <cv.h>
#include <highgui.h>
#include <stdio.h>

IplImage *Image1,*Image2;
IplImage *HSVImage;
IplImage *HueImage;
IplImage *BackProjectHueImage,*BackProjectImage;
CvHistogram *Histogram1;
IplImage *HistogramImage1;
CvPoint Point1,Point2;

int HueValue=0;

int HistogramBins = 180;
int HistogramBinWidth;
float HistogramRange1[2]={0,180};
float *HistogramRange[1]={&HistogramRange1[0]};
void onTrackbar(int position);

int main()
{

    Image1 =cvLoadImage("Riverbank.jpg",1);
    HSVImage = cvCreateImage( cvGetSize(Image1),8,3);
    HueImage = cvCreateImage( cvGetSize(Image1),8,1);

    BackProjectHueImage = cvCreateImage( cvGetSize(Image1),8,1);
    BackProjectImage = cvCreateImage( cvGetSize(Image1),8,3);

    Histogram1 = cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,HistogramRange);
    HistogramImage1 = cvCreateImage(cvSize(180,300),8,3);
    HistogramImage1->origin=1;

    cvCvtColor( Image1, HSVImage, CV_BGR2HSV );
    cvSplit(HSVImage,HueImage,0,0,0);
    cvCalcHist( &HueImage, Histogram1);
    cvNormalizeHist(Histogram1,5000);
    cvZero( HistogramImage1 );
    cvNot(HistogramImage1,HistogramImage1);
    HistogramBinWidth = HistogramImage1->width/HistogramBins;
    for(int i=0;i<HistogramBins;i++)
    {

        Point1=cvPoint(i,0);
        Point2=cvPoint(i,(int)cvQueryHistValue_1D(Histogram1,i));
        printf("%d\n",(int)cvQueryHistValue_1D(Histogram1,i));
        cvLine(HistogramImage1,Point1,Point2,CV_RGB(127,127,127));
    }
    cvNamedWindow("Riverbank",1 );
    cvNamedWindow("Hue Histogram",1);
    cvCreateTrackbar("Hue Thresh","Riverbank",&HueValue,250,onTrackbar);
    cvShowImage("Riverbank",Image1);
    cvShowImage("Hue Histogram",HistogramImage1);
    cvWaitKey(0);

}

void onTrackbar(int position)
{
    IplImage *Image2=cvCreateImage( cvGetSize(Image1),8,3);
    CvHistogram *Histogram2= cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,HistogramRange);
    cvCopyHist(Histogram1,&Histogram2);

    cvThreshHist(Histogram2,position);
    cvCalcBackProject(&HueImage, BackProjectHueImage, Histogram2);
    cvCopy(Image1,Image2,BackProjectHueImage);

    cvZero( HistogramImage1 );
    cvNot(HistogramImage1,HistogramImage1);
    HistogramBinWidth = HistogramImage1->width/HistogramBins;
    for(int i=0;i<HistogramBins;i++)
    {

        Point1=cvPoint(i,0);
        Point2=cvPoint(i,(int)cvQueryHistValue_1D(Histogram2,i));
        printf("%d\n",(int)cvQueryHistValue_1D(Histogram2,i));
        cvLine(HistogramImage1,Point1,Point2,CV_RGB(127,127,127));
    }
    cvShowImage("Hue Histogram",HistogramImage1);
    cvShowImage("Riverbank",Image2);
}

執行結果:


在OpenCV裡面,HSV色彩空間,色調(Hue)值的範圍在0~180,飽和度(Saturation)的範圍在0~255,亮度(Value)的範圍在0~255,而這邊就只取色調(Hue)值在做反向投影,開啟了一個Track bar的功能,並且利用cvCvtColor()將BGR的色彩空間轉換成HSV,並且用cvSplit()通道分割取Hue通道的圖片,計算Hue值的直方圖,在onTrackbar()的部份,則是用Trackbar來調整去除cvThreshHist()的最小區塊的臨界值,去除之後在反向投影到原始的Hue的圖片,在由反向投影的結果當做遮罩,直接跟彩色圖片做對應.

而cvCalcBackProject()不單單只有這樣的功能,它可以對多維度空間的色彩直方圖做對應,cvCalcBackProject()提供了一個當CvHistogram資料結構維度為3的時候的一個反向投影,下面的這個例子就以HSV的色彩空間為例,建構一個三維的CvHistogram資料結構

HSV三維直方圖反向投影
#include <cv.h>
#include <highgui.h>
#include <stdio.h>

IplImage *Image1,*Image2;
IplImage *HSV;
IplImage *HueImage,*SaturationImage,*ValueImage;
IplImage *ImageArray[3];
IplImage *BackProjectImage;
CvHistogram *Histogram1;
IplImage *HistogramImage1;
CvPoint Point1,Point2;

int HueValue=0;

int HistogramBins[3] ={180,256,256};
int HistogramBinWidth;
float HistogramRange1[6]={0,180,0,255,0,255};
float *HistogramRange[3]={&HistogramRange1[0],&HistogramRange1[2],&HistogramRange1[4]};
void onTrackbar(int position);

int main()
{

    Image1 =cvLoadImage("Riverbank.jpg",1);
    HSV = cvCreateImage( cvGetSize(Image1),8,3 );
    HueImage = cvCreateImage( cvGetSize(Image1),8,1 );
    SaturationImage = cvCreateImage( cvGetSize(Image1),8,1 );
    ValueImage = cvCreateImage( cvGetSize(Image1),8,1);
    ImageArray[0]=HueImage;
    ImageArray[1]=SaturationImage;
    ImageArray[2]=ValueImage;

    BackProjectImage = cvCreateImage( cvGetSize(Image1),8,3 );

    Histogram1 = cvCreateHist(3,HistogramBins,CV_HIST_SPARSE,HistogramRange);


    cvCvtColor( Image1, HSV, CV_BGR2HSV );
    cvSplit(HSV,HueImage,SaturationImage,ValueImage,0);

    cvCalcHist( ImageArray, Histogram1);


    cvNamedWindow("Riverbank",1 );
    cvCreateTrackbar("Hue Thresh","Riverbank",&HueValue,200,onTrackbar);
    cvShowImage("Riverbank",Image1);
    cvWaitKey(0);

}

void onTrackbar(int position)
{
    CvHistogram *Histogram2= cvCreateHist(3,HistogramBins,CV_HIST_SPARSE,HistogramRange);
    IplImage *Image2=cvCreateImage( cvGetSize(Image1),8,3 );
    IplImage *BackProjectImage = cvCreateImage( cvGetSize(Image1),8,1 );

    cvCopyHist(Histogram1,&Histogram2);

    cvThreshHist(Histogram2,position);
    cvCalcBackProject( ImageArray, BackProjectImage, Histogram2);
    cvCopy(Image1,Image2,BackProjectImage);

    cvShowImage("Riverbank",Image2);
}

執行結果:


在三維空間的作法上面,就要參考到前面"OpenCV統計應用-CvHistogram直方圖資料結構"關於三維空間製作的部份,除了用cvCvtColor()將色彩空間轉換,用cvSplit()將通道做分割,還要做個圖形陣列(ImageArray)來讓cvCalcHist()這個函式做運算,計算出來的結果為一個CvHistogram的三維空間稀疏矩陣直方圖,而在onTrackbar()的部份,cvCalcBackProject()直方圖反向投影亦是同樣要用ImageArray做輸入,而輸出則是一個單通道的圖形,在稀疏矩陣裡面,由於維度為三維,所以他所形成的統計直方圖數值都是極小,所以門檻值只要一點點就快要全部都分佈了,而這個三維空間的反向投影可以如此建構,是基於Look-up table的功能來實現,只不過他的缺點是,每一個維度的Look-up table範圍是0~255,因此如果是像Hue值的範圍0~180,它的數值就會被模糊化,也就是數據會被些許位移,但是仍不會影響它出來結果的精確度

在這個直方圖反向投影的部份,也可以結合連通成分來做去除某一門檻值的連通分量

cvCalcBackProject()
將統計直方圖的分布數據根據Look-up table對應回去,也就是說,當今天CvHistogram資料結構內的數據分佈是243,110,0,60...則使用cvCalcBackProject()函式單通道的圖片像素值會是,當遇到像素值為1的時候變243,像素值為2的時候變110,依此類推,cvCalcBackProject()直方圖反向投影可以根據多維度設計,而cvCalcBackProject()第一個引數為輸入單通道IplImage或CvMat資料結構,第二個引數為輸入單通道反向投影IplImage或CvMat資料結構,第三個引數為輸入CvHistogram資料結構
cvCalcBackProject(輸入單通道IplImage或CvMat資料結構,輸入單通道反向投影IplImage或CvMat資料結構,輸入CvHistogram資料結構)

OpenCV統計應用-影像增強,亮度/對比實作

在一般顯示螢幕以及圖形處理的應用軟體上,都會有一個亮度/對比的色彩(Brightness/Contrast)調整,它是屬於影像增強的部份,在OpenCV裡面的Sample Code裡面就有這樣的灰階程式的實作,在這邊就修改了OpenCV的Sample Code,來做色彩增強的亮度/對比的程式,而在一般的亮度/對比來講亮度(Brightness)的範圍為0~200而對比(Contrast)亦是0~200,它們由一條線性函數的公式所定義,對比所代表的是斜率,亮度則是偏移量,這條線性公式代表的是Look-up table的對應,它的數學式定義如下

原始的亮度對比數值範圍為-100~100之間,C代表對比,B代表亮度


對於對比率(Contrast ratio)來講,delta範圍應該落在0~255,這邊將對比率的公式做重新的調整



對比率代表著斜率的α值,而亮度則是決定線性公式位移的情況,也就是β值,而Y=αX+β這個線性公式它所表達的情況如下


α值的範圍落在0~255之間,而它的情況如下


再來下面是用虛擬碼的方式表達亮度/對比的演算法


下面就是亮度/對比的程式了

亮度/對比實作
#include <cv.h>
#include <highgui.h>
#include <stdio.h>

int BrightnessPosition = 100;
int ContrastPosition = 100;

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

IplImage *Image1,*Image2;
CvHistogram *Histogram1;
IplImage *HistogramImage;

uchar LookupTableData[256];
CvMat *LookupTableMatrix;
IplImage *LookupTableImage;
CvPoint Point1,Point2;


void OnTrackbar(int Position)
{
    int Brightness=BrightnessPosition-100;
    int Contrast=ContrastPosition -100;
    double Delta;
    double a,b;
    int y;

    //Brightness/Contrast Formula
    if(Contrast>0)
    {
        Delta=127*Contrast/100;
        a=255/(255-Delta*2);
        b=a*(Brightness-Delta);

        for(int x=0;x<256;x++)
        {
            y=(int)(a*x+b);

            if(y<0)
                y=0;
            if(y>255)
                y=255;

            LookupTableData[x]=(uchar)y;
        }
    }
    else
    {
        Delta=-128*Contrast/100;
        a=(256-Delta*2)/255;
        b=a*Brightness+Delta;

        for(int x=0;x<256;x++)
        {
            y=(int)(a*x+b);

            if(y<0)
                y=0;
            if(y>255)
                y=255;

            LookupTableData[x]=(uchar)y;
        }
    }
    //End

    //Look up table sketch
    cvSetZero(LookupTableImage);
    cvNot(LookupTableImage,LookupTableImage);
    Point2=cvPoint(0,LookupTableData[0]);
    for(int i=0;i<256;i++)
    {
        Point1=cvPoint(i,LookupTableData[i]);
        cvLine(LookupTableImage,Point1,Point2,CV_RGB(0,0,0),3);
        Point2=Point1;
    }
    cvLUT(Image1,Image2,LookupTableMatrix);
    //End

    //Gray Level Histogram
    cvCalcHist(&Image2,Histogram1);
    cvNormalizeHist(Histogram1,3000);

    cvSetZero(HistogramImage);
    cvNot(HistogramImage,HistogramImage);
    HistogramBinWidth = HistogramImage->width/HistogramBins;
    for(int i=0;i<HistogramBins;i++)
    {
        Point1=cvPoint(i*HistogramBinWidth,0);
        Point2=cvPoint((i+1)*HistogramBinWidth,(int)cvQueryHistValue_1D(Histogram1,i));
        cvRectangle(HistogramImage,Point1,Point2,CV_RGB(0,0,0),CV_FILLED);
    }
    //End

    cvShowImage("Gray Level Histogram",HistogramImage);
    cvShowImage("Brightness/Contrast",Image2);
    cvShowImage("Image Enhance",LookupTableImage);
    cvZero(Image2);
}

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

    Histogram1=cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,HistogramRange);
    HistogramImage = cvCreateImage(cvSize(320,200),8,1);

    LookupTableImage=cvCreateImage(cvSize(256,256),8,3);
    LookupTableMatrix=cvCreateMatHeader(1,256,CV_8UC1);
    cvSetData(LookupTableMatrix,LookupTableData,0);

    LookupTableImage->origin=1;
    HistogramImage->origin=1;

    cvNamedWindow("Brightness/Contrast",1);
    cvNamedWindow("Gray Level Histogram",1);
    cvNamedWindow("Image Enhance",1);

    cvCreateTrackbar("brightness","Brightness/Contrast",&BrightnessPosition,200,OnTrackbar);
    cvCreateTrackbar("contrast","Brightness/Contrast",&ContrastPosition,200,OnTrackbar);

    OnTrackbar(0);

    cvWaitKey(0);
}

執行結果:


這隻程式同樣也是用到CvHistogram資料結構,使用到兩個拉軸(Trackbar),以及Look-up table的應用,在//Brightness/Contrast Formula的註解內所包的就是亮度/對比演算法虛擬碼的實作,再來就是把它的線性系統化出來,也就是Y=αX+β的函數方程式,這個方程式,當然同等於Look up table,而之後,在把他們灰階直方圖的分布畫出來,在main()裡面,當然是先讀取目標圖片轉成灰階,初始化繪製直方圖與線性系統圖片的空間,創立三個視窗介面,設立兩個拉軸,並且將拉軸的事件函式設定成同一個的副程式的名稱.而對於影像增強(Image Enhance)這個視窗介面,它所代表的含意如下



X軸代表為是原始灰階的輸入值,而Y軸代表的是灰階值所對應的結果,而X軸跟Y軸的範圍都是0~255,而這條直線公式也會受到斜率(α)以及平移(β)的結果改面灰階值輸入以及輸出的對應,它是將一張原始灰階圖片的每一個像素值做線性函式的對應,使得每個灰階值對應出來的結果產生了變化,由下面可以知道它(LUT)對應的關係

(a)亮度條為0因此小於100的灰階值都為0而灰階值方圖也像左偏移


(b)亮度條為100,因此大於156的灰階值都為255,而灰階值方圖也都向右偏移


(c)對比為0,這個時候斜率α為0,因此輸入的0~255的灰階值輸出都固定為128,因此整張圖片都是灰階值128的影像,而灰階直方圖則是所有數據都集中在128


(d)對比為100,這個時候斜率為255,而這樣的圖片又可以叫做二值化圖片,因為輸出結果非黑即白,而移動亮度則是在平移二值化的門檻值,由灰階值方圖可以得知,所有數據都被分開到0跟255兩邊


上面所表達的其實就是一種Look-up table的表達方式,藉由一個輸入灰階值的矩陣,對應岀另一個不同的灰階值數據,因此改變了原始灰階值的數據,而整張圖片也因此產生了變化