OpenCV臉部跟蹤(2)

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

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

前面一篇文章中提到，我們在一副臉部圖像上選取76個(gè)特征點(diǎn)來描述臉部形狀特征，本文中我們會把這些特征點(diǎn)映射到一個(gè)標(biāo)準(zhǔn)形狀模型。

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

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

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

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

//該函數(shù)用人臉模型V和e，把向量p轉(zhuǎn)化為點(diǎn)集
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來處理特征點(diǎn)，移去全局剛性變化，Procrustes analysis算法可以參考：http://en.wikipedia.org/wiki/Procrustes_analysis

在數(shù)學(xué)上，Procruster analysis就是尋找一個(gè)標(biāo)準(zhǔn)形狀，然后把所有其它特征點(diǎn)數(shù)據(jù)都和標(biāo)準(zhǔn)形狀對齊，對齊的時(shí)候采用最小平方距離，用迭代的方法不斷逼近。下面通過代碼來了解如何實(shí)現(xiàn)Procrustes analysis。

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

)
{

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

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

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

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

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

該函數(shù)的代碼：

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對齊的特征向量，我們要用一個(gè)統(tǒng)一的矩陣把平移和旋轉(zhuǎn)統(tǒng)一起來表示（成為線性表示），然后把該矩陣追加到局部變形空間，注意對該矩陣表示，我們最后進(jìn)行了史密斯正交處理。

我們通過函數(shù) 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函數(shù)的實(shí)現(xiàn)：

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

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

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

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是樣本的數(shù)目，n是76，表示76個(gè)特征點(diǎn)。
int N = X.cols,n = X.rows/2;

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

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

臉部局部變形我們用一個(gè)線性模型表示，主要的思想如下圖所示：一個(gè)有N個(gè)面部特征組成面部形狀，被建模成一個(gè)2N維空間的點(diǎn)。我們要盡量找到一個(gè)低維的超平面，在這個(gè)平面內(nèi)，所有的面部形狀都在里面，由于這個(gè)超平面只是2N維空間的子集，占用剛少的空間，處理起來更快?？梢韵氲玫剑绻@個(gè)子空間來自于一個(gè)人，則子空間的點(diǎn)，表現(xiàn)這個(gè)人的各種表情變化。前面的教程中，我們知道PCA算法能夠找到低維子空間，但PCA算法需要指定子空間的維數(shù)，在啟發(fā)式算法中有時(shí)候這個(gè)值很難選擇。在本程序中，我們通過SVD算法來模擬PCA算法。

//compute non-rigid transformation

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

奇異值分解SVD有效的應(yīng)用到形狀數(shù)據(jù)的協(xié)方差矩陣（即，dY.t()*dY），OpenCV的SVD類的w成員存儲著數(shù)據(jù)變化性的主要方向的變量，從最大到最小排序。一個(gè)選擇子空間維數(shù)的普通方法是選擇保存數(shù)據(jù)總能量分?jǐn)?shù)frac的方向最小集（即占總能量的比例為frac），這是通過svd.w記錄表示的，因?yàn)檫@些記錄是從最大的到最小的排序的，它充分地用來評估子空間，通過用變化性方向的最大值k來評估能量。他們自己的方向存儲在SVD類的u成員內(nèi)。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;

達(dá)到了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個(gè)特征向量
Mat D = svd.u(Rect(0,0,k,2*n));

把全局剛體運(yùn)動和局部變形運(yùn)動結(jié)合起來，注意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);

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

//compute variance (normalized wrt scale)

Mat Q = V.t()*X; //把數(shù)據(jù)投影到子空間
for(int i = 0; i < N; i++) //normalize coordinates w.r.t scale
{ //用第一個(gè)坐標(biāo)縮放，防止太大的縮放值影響臉部識別
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,

程序的運(yùn)行參數(shù)為：annotations.yaml shapemodle.yaml

程序執(zhí)行后，可以看到我們只保留了14個(gè)模型。

我們也可以使用下面的運(yùn)行參數(shù)：annotations.yaml shapemodle.yaml –mirror

這時(shí)候，每副圖像的特征點(diǎn)，會生成一個(gè)y軸對稱的鏡像特征點(diǎn)集，這時(shí)訓(xùn)練的采樣數(shù)目翻倍，為5828。

在工程文件FirstOpenCV41中，我們可視化了生成的模型，會連續(xù)顯示14個(gè)模型的不同姿態(tài)：

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