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_<float>(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个模型的不同姿态：