Procrustes Analysis普氏分析法

 

 

     選取N幅同類目標物體的二維圖像，並用上一篇博文的方法標註輪廓點，這樣就獲得訓練樣本集：

     因爲圖像中目標物體的形狀和位置存在較大誤差，所以所獲得的數據並不具備仿射不變性，須要對其進行歸一化處理。這裏採用Procrustes分析方法對樣本集中的全部形狀集合進行歸一化。形狀和位置的載體仍是樣本點的空間座標。

     普氏分析法是一種用來分析形狀分佈的方法。數學上來說，就是不斷迭代，尋找標準形狀(canonical shape)，並利用最小二乘法尋找每一個樣本形狀到這個標準形狀的仿射變化方式。（可參照維基百科的GPA算法）

 

本書中，兩個形狀的歸一化過程（一個形狀爲canonical shape，另外一個爲樣本形狀）：

(1)    求每一個樣本點i(i=1,2..,n)在N幅圖像中的均值

 

(2)    對全部形狀的大小進行歸一化，即將每一個樣本點減去其對應均值

 

(3)    根據去中心化數據，計算每幅圖像中形狀的重心，對於第i幅圖像，其重心爲：

 

(4)    根據重心和角度，將標準和樣本形狀對齊在一塊兒，使得兩個形狀的普氏距離最小，下式爲普氏距離定義：

 

 這個第(4)步的具體作法，不斷迭代如下過程：

(a)經過計算每幅圖像中全部歸一化樣本點的平均值獲得每一個圖像的標準形狀canonical shape。

(b)利用最小二乘法求每一個圖像中樣本形狀到標準形狀的旋轉角度。根據普氏距離的定義，也就是求：

 

其中的a和b表示仿射變換裏旋轉變化的參數：

 

對上式求偏導數，能夠獲得所求的a和b：

 

 (c)根據旋轉參數，對樣本形狀作旋轉變化，獲得和標準形狀對齊的新的形狀

 

 (d)重複以上步驟，直到達到指定循環次數或者先後兩次迭代之間canonical shape的絕對範數知足必定閾值

 

   Procrustesanalysis的做用能夠看做是一種對原始數據的預處理，目的是爲了獲取更好的局部變化模型做爲後續模型學習的基礎。以下圖所示，每個人臉特徵點能夠用一種單獨的顏色表示；通過歸一化變化，人臉的結構愈來愈明顯，即臉部特徵簇的位置愈來愈接近他們的平均位置；通過一系列迭代，尺度和旋轉的歸一化操做，這些特徵簇變得更加緊湊，它們的分佈愈來愈能表達人臉表情的變化。【剔除剛性部分、保留柔性部分】

 

 

      上圖爲不一樣大小、不一樣長寬比的矩形，通過歸一化過程後，各個樣本點分佈服從必定機率分佈趨勢。

 

具體代碼實現：

函數：Mat shape_model::procrustes(const Mat &X,const int itol,constfloat ftol)

矩陣X：2n*N維矩陣，即有N副圖像，每幅圖像有n個特徵點（2n是由於x和y座標）

Itol：迭代次數

Ftol：求絕對範數的閾值

返回值：矩陣P，存放對齊後的形狀座標

 1 Mat shape_model::procrustes(const Mat &X,const int itol,const float ftol)
 2 {
 3   /* X.cols:特徵點個數，X.rows: 圖像數量*2 */
 4   int N = X.cols,n = X.rows/2; 
 5 
 6   //remove centre of mass
 7   Mat P = X.clone();
 8   for(int i = 0; i < N; i++)
 9   {
10     /*取X第i個列向量*/
11     Mat p = P.col(i);
12     float mx = 0,my = 0;
13     for(int j = 0; j < n; j++)
14     {
15         mx += p.fl(2*j);
16         my += p.fl(2*j+1);
17     }
18     /*分別求圖像集，2維空間座標x和y的平均值*/
19     mx /= n;
20     my /= n;
21     /*對x,y座標去中心化*/
22     for(int j = 0; j < n; j++)
23     {
24         p.fl(2*j) -= mx; 
25         p.fl(2*j+1) -= my;
26     }
27   }
28   //optimise scale and rotation
29   Mat C_old;
30   for(int iter = 0; iter < itol; iter++)
31   {   
32     /*計算（含n個形狀）的重心*/  
33     Mat C = P*Mat::ones(N,1,CV_32F)/N; 
34     /*C爲2n*1維矩陣，含n個重心，對n個重心歸一化處理*/
35     normalize(C,C);
36     if(iter > 0)
37     {
38       if(norm(C,C_old) < ftol)//norm:求絕對範數，小於閾值，則退出循環
39         break;
40     }
41     C_old = C.clone();
42     for(int i = 0; i < N; i++)
43     {
44       //求當前形狀與歸一化重心之間的旋轉角度，即上式a和b
45       Mat R = this->rot_scale_align(P.col(i),C); 
46       for(int j = 0; j < n; j++)
47       {
48         float x = P.fl(2*j,i),y = P.fl(2*j+1,i);
49         /*仿射變化*/
50         P.fl(2*j  ,i) = R.fl(0,0)*x + R.fl(0,1)*y;
51         P.fl(2*j+1,i) = R.fl(1,0)*x + R.fl(1,1)*y;
52       }
53     }
54   }
55 return P;
56 }http://blog.csdn.net/jinshengtao/article/details/43376049