OpenCV 脸部跟踪(2)
前面一篇文章中提到，我们在一副脸部图像上选取76个特征点，以及这些特征点的连通性信息来描述脸部形状特征，本文中我们会把这些特征点映射到一个标准形状模型。
通常，脸部形状特征点能够参数化分解为两个变量，一个是全局的刚体变化，一个是局部的变形。全局的刚体变化主要是指脸部能够在图像中移动，旋转，缩放，局部的变形则是指脸部的表情变化，不同人脸的特征等等。
下面我们通过train函数，一步步了解下如何把标记点特征数据转化为标准形状模型：
本文参考了下面两篇文章：
http://blog.csdn.net/raby_gyl/article/details/13148193
http://blog.csdn.net/raby_gyl/article/details/13024867
train函数的输入为n个样本图像的采样特征点，该点集会被首先转化为行152，列为样本数量的矩阵表示，另外还有连通性点集索引，以及方差的置信区间以及保留模型的最大数量，后两个参数主要在模型空间降维时候使用。
void train(const vector<vector<Point2f> > &p, //N-example shapes
const vector<Vec2i> &con = vector<Vec2i>(), //point-connectivity
const float frac = 0.95, //fraction of variation to retain
const int kmax = 10) //maximum number of modes to retain
{
Mat X = this->pts2mat(points);
二维向量组成的特征点集points经过函数pts2mat转化为矩阵X，X矩阵为152行，N列，N表示样本图像的数量，每一列为一副图像的76个特征点的x，y坐标。
//N = 2914*2 n = 76，共2914副图像特征点，但因为我们使用镜像，所以采样图像特征数据翻了一倍。
int N = X.cols,n = X.rows/2;
Mat Y = this->procrustes(X);
上面的函数中，我们通过Procrustes analysis来处理特征点，Procrustes analysis算法可以参考：http://en.wikipedia.org/wiki/Procrustes_analysis
在数学上，Procruster analysis就是寻找一个标准形状，然后把所有其它特征点数据都和标准形状对齐，对齐的时候采用最小平方距离，用迭代的方法不断逼近。
下面通过代码来了解如何实现Procrustes analysis， http://www.cnblogs.com/mikewolf2002/p/3659964.html
经过Procrustes analysis变化后得到矩阵Y，矩阵Y仍是152*5828的矩阵，但它的每列表示的特征向量都接近于标准形状。下面我们要做的就是根据Y计算刚体变化矩阵R。R是一个152*4的矩阵。
//计算得到刚体变化矩阵R
Mat R = this->calc_rigid_basis(Y);
计算R的原理如下：
通过Procrustes analysis对齐的特征向量，我们要用一个统一的矩阵把平移和旋转统一起来表示（成为线性表示），然后把该矩阵追加到局部变形空间，注意对该矩阵表示，我们最后进行了史密斯正交处理。这儿x1,y1,....xn,yn是Y矩阵的平均列向量，也就是标准形状(权威形状)。
该函数的代码如下：
Mat shape_model::calc_rigid_basis(const Mat &X)
{
//计算均值形状
int N = X.cols,n = X.rows/2;
Mat mean = X*Mat::ones(N,1,CV_32F)/N;
注意这儿的x1，y1,..., xn, yn是Procrustes结果的平均值
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正交化处理
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函数剩余的代码，前面计算了刚体变化矩阵R，下面计算非刚体变化。
//计算非刚体变化
Y是152*5828的矩阵，它是procrustes分析的结果，R是刚体变化矩阵152*4，它的转置就是4*152 。
我们先了解一个矩阵相乘的概念，假设V是一个空间，那么：
s=V.t()*p; //表示矢量p在V空间下的坐标投影。
p=V*s,表示由坐标（权重）和与之对应的V的列向量，叠加，组合而成的一个形状矢量。
特别的，如果V.t()==V.inv(),即矩阵转置等于矩阵逆，则V为正交矩阵，也即V.t()*V.inv()=E，此时上式就是完全的可逆变换了。
则，R*P为152*5828
Mat dY = Y - R*P; //dy变量的每一列表示减去均值的Procrustes对齐形状，投影刚体运动
这儿Y-R*P后的dY表示Y减去刚性变化，则dY表示非刚性变化。
非刚性空间dY是152*5828的高维矩阵，我们需要对其进行降维操作，这儿采用矩阵奇异分解法来进行降维操作：
奇异值分解SVD有效的应用到形状数据的协方差矩阵（即，dY.t()*dY），OpenCV的SVD类的w成员存储着数据变化性的主要方向的变量，从最大到最小排序。一个选择子空间维数的普通方法是选择保存数据总能量分数frac的方向最小集（即占总能量的比例为frac），这是通过svd.w记录表示的，因为这些记录是从最大的到最小的排序的，它充分地用来评估子空间，通过用变化性方向的最大值k来评估能量。他们自己的方向存储在SVD类的u成员内。svd.w和svd.u成分一般分别被成为特征值和特征矢量。 http://www.cnblogs.com/mikewolf2002/category/351010.html
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个特征向量，这儿k=1，表示第一个特征向量已经占到了95%的主成分量，m=14，为最大保留14个特征向量（如果使用了—mirror参数k=1,如果没有，则为14)。
Mat D = svd.u(Rect(0,0,k,2*n));
D矩阵为152行，1列。
把全局刚体运动和局部变形运动结合起来，组成一个统一矩阵V。
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; //把数据投影到子空间，Q为5*5828矩阵
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为参数方差，注意方差是在子空间标准化后的坐标上计算的，标准化是关于第一维的坐标（即尺度)。这样可以阻止相对大尺度的样本数据主导估计。也注意到一个负值被分配到刚性子空间坐标的方差（即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<int>(i,0) = con[i][0];
C.at<int>(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<int>(i,0) = i; C.at<int>(i,1) = i+1;
}
C.at<int>(n-1,0) = n-1; C.at<int>(n-1,1) = 0;
}
}
形状模型类主要成员如下：
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<Point2f> &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<Point2f> calc_shape(); //shape described by parameters @p
工程文件：FirstOpenCV40,
程序的运行参数为：annotations.yaml shapemodle.yaml
程序执行后，可以看到我们只保留了14个模型。
我们也可以使用下面的运行参数：annotations.yaml shapemodle.yaml –mirror
这时候，每副图像的特征点，会生成一个y轴对称的镜像特征点集，这时训练的采样数目翻倍，为5828。此时竟然值保留了一个模型。
在工程文件FirstOpenCV41中，我们可视化了生成的模型，会连续显示14个模型的不同姿态：
