OpenCV臉部跟蹤(2)

【原文：http://www.cnblogs.com/mikewolf2002/p/3474188.html】 前面一篇文章中提到，我們在一副臉部圖像上選取76個特征點來描述臉部形狀特征，本文中我們會把這些特征點映射到一個標準形狀模型。 通常，臉部形狀特征點能夠參數化分解為兩個變量，一個是全

【原文：http://www.cnblogs.com/mikewolf2002/p/3474188.html】

前面一篇文章中提到，我們在一副臉部圖像上選取76個特征點來描述臉部形狀特征，本文中我們會把這些特征點映射到一個標準形狀模型。

通常，臉部形狀特征點能夠參數化分解為兩個變量，一個是全局的剛體變化，一個是局部的變形。全局的剛體變化主要是指臉部能夠在圖像中移動，旋轉，縮放，局部的變形則是指臉部的表情變化，不同人臉的特征等等。

形狀模型類主要成員如下：

class shape_model
{ //2d linear shape model
public:
Mat p; //parameter vector (kx1) CV_32F，參數向量
Mat V; //shape basis (2nxk) CV_32F, line subspace,線性子空間
Mat e; //parameter variance (kx1) CV_32F 參數方差
Mat C; //connectivity (cx2) CV_32S 連通性

//把一個點集投影到一個可信的臉部形狀空間
void calc_params(const vector &pts, //points to compute parameters from
const Mat weight = Mat(), //weight of each point (nx1) CV_32F 點集的權重
const float c_factor = 3.0); //clamping factor

//該函數用人臉模型V和e，把向量p轉化為點集
vector calc_shape(); //shape described by parameters @p

...

void train(const vector > &p, //N-example shapes
const vector &con = vector(), //point-connectivity
const float frac = 0.95, //fraction of variation to retain
const int kmax = 10); //maximum number of modes to retain
...

}

本文中，我們通過Procrustes analysis來處理特征點，移去全局剛性變化，Procrustes analysis算法可以參考：http://en.wikipedia.org/wiki/Procrustes_analysis

在數學上，Procruster analysis就是尋找一個標準形狀，然后把所有其它特征點數據都和標準形狀對齊，對齊的時候采用最小平方距離，用迭代的方法不斷逼近。下面通過代碼來了解如何實現Procrustes analysis。

//Procrustes分析的基本思想是最小化所有形狀到平均形狀的距離和
Mat shape_model::procrustes(const Mat &X,
const int itol, //最大迭代次數
const float ftol //精度

)
{

X矩陣就是多副樣本圖像76個特征點組成的矩陣，共152行，列數為圖像的個數，每列表示一個樣本圖像特征點的x，y坐標。
int N = X.cols,n = X.rows/2;

//remove centre of mass
//所有的形狀向量(特征)對齊到原點，即每個向量的分量減去其平均值，每列是一個形狀向量。
Mat P = X.clone();
for(int i = 0; i < N; i++)
{
Mat p = P.col(i); //第i個向量
float mx = 0,my = 0;
for(int j = 0; j < n; j++) //x，y分別計算得到平均值。
{
mx += p.fl(2*j);
my += p.fl(2*j+1);
}
mx /= n; my /= n;
for(int j = 0; j < n; j++)
{
p.fl(2*j) -= mx;
p.fl(2*j+1) -= my;
}
}

//optimise scale and rotation
Mat C_old;
for(int iter = 0; iter < itol; iter++)
{

注意下邊的一行代碼，會把每個圖像對齊到原點特征點x，y分別加起來，求平均值，得到一個152*1的矩陣，然后對該矩陣進行歸一化處理。
Mat C = P*Mat::ones(N,1,CV_32F)/N; //計算形狀變換后的平均值
normalize(C,C); //canonical shape (對x-進行標準化處理)
if(iter > 0) //converged?//收斂？當兩個標準形狀或者標準形狀的誤差小于某一值這里是ftol則，停止迭代。
{

norm函數默認是矩陣各元素平方和的根范式
if(norm(C,C_old) < ftol)
break;
}
C_old = C.clone(); //remember current estimate//記下當前的矩陣，和下一次進行比較
for(int i = 0; i < N; i++)
{

rot_scale_align函數求每副圖像的特征點向量和平均向量滿足最小乘法時候的旋轉矩陣。即求得a,b值組成的旋轉矩陣。

該函數的代碼：

Mat shape_model::rot_scale_align(const Mat &src, const Mat &dst)
 {
 //construct linear system
 int n = src.rows/2;
 float a=0,b=0,d=0;
 for(int i = 0; i < n; i++)
 {
 d += src.fl(2*i) * src.fl(2*i ) + src.fl(2*i+1) * src.fl(2*i+1); //x*x+y*y
 a += src.fl(2*i) * dst.fl(2*i ) + src.fl(2*i+1) * dst.fl(2*i+1);
 b += src.fl(2*i) * dst.fl(2*i+1) - src.fl(2*i+1) * dst.fl(2*i );
 }
 a /= d; 
 b /= d;
 //solved linear system
 return (Mat_(2,2) << a,-b,b,a);
 }

通過Procrustes analysis對齊的特征向量，我們要用一個統一的矩陣把平移和旋轉統一起來表示（成為線性表示），然后把該矩陣追加到局部變形空間，注意對該矩陣表示，我們最后進行了史密斯正交處理。

我們通過函數 calc_rigid_basis得到該矩陣表示：

Mat shape_model::calc_rigid_basis(const Mat &X)
{
//compute mean shape
int N = X.cols,n = X.rows/2;
Mat mean = X*Mat::ones(N,1,CV_32F)/N;

//construct basis for similarity transform
Mat R(2*n,4,CV_32F);
for(int i = 0; i < n; i++)
{
R.fl(2*i,0) = mean.fl(2*i );
R.fl(2*i+1,0) = mean.fl(2*i+1);
R.fl(2*i,1) = -mean.fl(2*i+1);
R.fl(2*i+1,1) = mean.fl(2*i );
R.fl(2*i,2) = 1.0;
R.fl(2*i+1,2) = 0.0;
R.fl(2*i,3) = 0.0;
R.fl(2*i+1,3) = 1.0;
}
//Gram-Schmidt orthonormalization
for(int i = 0; i < 4; i++)
{
Mat r = R.col(i);
for(int j = 0; j < i; j++)
{
Mat b = R.col(j);
r -= b*(b.t()*r);
}
normalize(r,r);
}
return R;
}

下面我們看看train函數的實現：

兩篇參考的翻譯：http://blog.csdn.net/raby_gyl/article/details/13148193

http://blog.csdn.net/raby_gyl/article/details/13024867

該函數的輸入為n個樣本圖像的采樣特征點，該點集會被首先轉化為行152，列為樣本數量的矩陣表示，另外還有連通性點集索引，以及方差的置信區間以及保留模型的最大數量。

void train(const vector > &p, //N-example shapes
const vector &con = vector(), //point-connectivity
const float frac = 0.95, //fraction of variation to retain
const int kmax = 10) //maximum number of modes to retain
{
//vectorize points
Mat X = this->pts2mat(points);

N是樣本的數目，n是76，表示76個特征點。
int N = X.cols,n = X.rows/2;

//align shapes
Mat Y = this->procrustes(X);

//compute rigid transformation 計算得到剛體變化矩陣R
Mat R = this->calc_rigid_basis(Y);

臉部局部變形我們用一個線性模型表示，主要的思想如下圖所示：一個有N個面部特征組成面部形狀，被建模成一個2N維空間的點。我們要盡量找到一個低維的超平面，在這個平面內，所有的面部形狀都在里面，由于這個超平面只是2N維空間的子集，占用剛少的空間，處理起來更快。可以想得到，如果這個子空間來自于一個人，則子空間的點，表現這個人的各種表情變化。前面的教程中，我們知道PCA算法能夠找到低維子空間，但PCA算法需要指定子空間的維數，在啟發式算法中有時候這個值很難選擇。在本程序中，我們通過SVD算法來模擬PCA算法。

//compute non-rigid transformation

Y是152*1294的矩陣，它是procrustes分析的結果，R是剛體變化矩陣152*4，它的轉置就是4*152
Mat P = R.t()*Y; //原始的位置
Mat dY = Y - R*P; //dy變量的每一列表示減去均值的Procrustes對齊形狀，投影剛體運動

奇異值分解SVD有效的應用到形狀數據的協方差矩陣（即，dY.t()*dY），OpenCV的SVD類的w成員存儲著數據變化性的主要方向的變量，從最大到最小排序。一個選擇子空間維數的普通方法是選擇保存數據總能量分數frac的方向最小集（即占總能量的比例為frac），這是通過svd.w記錄表示的，因為這些記錄是從最大的到最小的排序的，它充分地用來評估子空間，通過用變化性方向的最大值k來評估能量。他們自己的方向存儲在SVD類的u成員內。svd.w和svd.u成分一般分別被成為特征值和特征矢量。
SVD svd(dY*dY.t());
int m = min(min(kmax,N-1),n-1);
float vsum = 0;
for(int i = 0; i < m; i++)
vsum += svd.w.fl(i);
float v = 0;
int k = 0;

達到了95%的主成分量，退出，frac=0.95
for(k = 0; k < m; k++)
{
v += svd.w.fl(k);
if(v/vsum >= frac){k++; break;}
}
if(k > m) k = m;

取前k個特征向量
Mat D = svd.u(Rect(0,0,k,2*n));

把全局剛體運動和局部變形運動結合起來，注意V的第一列是縮放，第三、四列分別是x，y偏移。

//combine bases
V.create(2*n,4+k,CV_32F);
Mat Vr = V(Rect(0,0,4,2*n)); //剛體子空間
R.copyTo(Vr); //非剛體子空間
Mat Vd = V(Rect(4,0,k,2*n));
D.copyTo(Vd);

最后我們要注意的一點是如何約束子空間坐標，以使得子空間內的面部形狀都是有效的。在下面的圖中，我們可以看到，對于子空間內的圖像，如果在某個方向改變坐標值，當坐標值小時候，它仍是一個臉的形狀，但是變化值大時候，就不知道是什么玩意了。防止出現這種情況的最簡單方法，就是把變化的值clamp在一個范圍內，通常是現在± 3 的范圍，這樣可以cover到99.7%的臉部變化。clamping的值通過下面的代碼計算：

//compute variance (normalized wrt scale)

Mat Q = V.t()*X; //把數據投影到子空間
for(int i = 0; i < N; i++) //normalize coordinates w.r.t scale
{ //用第一個坐標縮放，防止太大的縮放值影響臉部識別
float v = Q.fl(0,i);
Mat q = Q.col(i);
q /= v;
}
e.create(4+k,1,CV_32F);
pow(Q,2,Q);
for(int i = 0; i < 4+k; i++)
{
if(i < 4)
e.fl(i) = -1; //no clamping for rigid coefficients
else
e.fl(i) = Q.row(i).dot(Mat::ones(1,N,CV_32F))/(N-1);
}
//store connectivity
if(con.size() > 0)
{ //default connectivity
int m = con.size();
C.create(m,2,CV_32F);
for(int i = 0; i < m; i++)
{
C.at(i,0) = con[i][0];
C.at(i,1) = con[i][1];
}
}
else
{ //user-specified connectivity
C.create(n,2,CV_32S);
for(int i = 0; i < n-1; i++)
{
C.at(i,0) = i; C.at(i,1) = i+1;
}
C.at(n-1,0) = n-1; C.at(n-1,1) = 0;
}
}

工程文件：FirstOpenCV40,

程序的運行參數為：annotations.yaml shapemodle.yaml

程序執行后，可以看到我們只保留了14個模型。

我們也可以使用下面的運行參數：annotations.yaml shapemodle.yaml –mirror

這時候，每副圖像的特征點，會生成一個y軸對稱的鏡像特征點集，這時訓練的采樣數目翻倍，為5828。

在工程文件FirstOpenCV41中，我們可視化了生成的模型，會連續顯示14個模型的不同姿態：

聲明：本網頁內容旨在傳播知識，若有侵權等問題請及時與本網聯系，我們將在第一時間刪除處理。TEL:177 7030 7066 E-MAIL:11247931@qq.com