《Master Opencv...读书笔记》非刚性人脸跟踪 III

   上篇文章中，我们获得了人脸的各种表情模式，也就是一堆标注点的形变参数。这次我们需要训练一中人脸特征（团块模型），它能够对人脸的不同部位(即“标注点”)分别进行描述，作为后面人脸跟踪、表情识别的区分依据。本次博文的主要内容：

a.      介绍下人脸特征检测器大概有哪些类别

b.      详细介绍随机梯度法，并介绍在人脸团块特征提取时的应用

c.      为了提高训练/跟踪的健壮性，利用上一讲对输入的图像进行大小、角度的约束

• 人脸特征检测器综述

      人脸特征检测与普通的物体检测非常相似，opencv工具也为人脸检测提供了一系列复杂功能接口。但它们的不同之处，大致可以从以下三个维度进行考虑：

1.      精确性与健壮性程度（Precisionversus robustness）

对于普通的目标检测，只是从图像中找到目标对象的大致位置。而对于人脸特征的检测需要对特征的位置进行高度精确的估计。少量像素位置的错误估计也许对于普通目标跟踪无关紧要，但是对于人脸特征跟踪这意味着后面无法区分人脸的表情（微笑与皱眉）。

2.      数据支撑程度（Ambiguityfrom limited spatial support）

大家普遍认为，对于普通的目标检测，它的ROI区域包含了足够多的图像信息，我们可以通过它非常稳定的从图像中剔除非目标区域。但是人脸特征的情况则不容乐观，图像中包含人脸特征的区域非常有限，少量的数据不足以支撑人脸特征的结构表达。比如人脸轮廓的特征，非常容易和其他包含强烈边缘的图像块混淆。

3.计算复杂度（Computationalcomplexity）

  普通的目标检测侧重寻找图像中所有符合要求的目标对象实体的位置，而人脸跟踪则需要知道每个脸部特征的位置（通常包含20-100个特征点才能描述人脸，可认为是几何特征点），有效评估每个人脸特征点检测的计算能力对于实时的人脸跟踪非常重要。

   由于以上差别，经常有普通的目标检测算法被用于人脸跟踪，但是这种做法并不被业内人事认可，毕竟它们无法精确表达人脸结构。

这次我们将用线性的图像团块来表达人脸特征，虽然它的构造非常简单，通过精心设计它的学习过程，我们了解到这种表达方式在人脸跟踪时能合理的估计人脸特征的位置。此外，由于设计简单，它的计算速度非常快，使得实时人脸跟踪成为了可能。

     而在特征检测器训练算法主要有两类：侧重于描述对象外观结构的算法（generative）和侧重于从多个物体中区分出目标对象的算法（discriminative）【我认为一个尽可能复原整个对象特征，一个仅仅是根据某几个特征来挑选目标】。Generative算法的优点是结果模型通过对特定对象的属性进行编码使得其图像的细节在视觉上可以被观察到，比如特征脸Eigenfaces：http://blog.csdn.net/jinshengtao/article/details/18599165。discriminative算法的优点是模型集合的全集包含所有对象实例，直接面向求解的问题，比较代表性的算法SVM：http://blog.csdn.net/jinshengtao/article/details/40900865。

（注意：特征脸和支持向量机算法设计的初衷是分类，而不是检测与图像配准。但是这两个算法所运用的数学技巧证明它们可以被运用于人脸跟踪领域）

    虽然很多场合上述两种算法都能取得不错的效果，但是本次以图像块对人脸特征进行建模，利用discriminative类型的算法得到的效果更好。接下来就详细介绍这种团块模型(Patch Model)的训练过程.

•  discriminative patch models

1.      是什么patch models？

   它是一种特征模版，当用它覆盖在原始图像上进行搜索时，含有人脸特征的区域在该模版上会有强烈的反馈，而不含人脸特征的区域反馈则较弱，可以用如下数学公式表达：

  上式中矩阵I_i表示第i幅样本图像，I_i(.)表示手工标注的包含人脸特征的样本区域（TrainingPatch）的范围（之前是手工标准的几何特征点，现在由点化面啦~）；矩阵P就是团块模型(Patch Model)，其长和宽分别为w和h；矩阵R表示理想的反馈结果。

  我们本节的目标是像来训练团块特征模版（为了完整描述人脸所有细节特征，通常需要多个这样的团块特征模版，因此称为correlation-basedpatch models）。这里也就是利用最小二乘法从上式中求解矩阵P（矩阵R和I是不变量，P是自变量）。上述公式表达了在所有经特征模版P扫描过的区域而产生的响应图像，平均来说最接近理想响应图像(二者的差最小)。

  那么理想反馈图像R该如何选取呢？最直接的方式，构造一个周围全0而中间非零的矩阵。并且这种做法的前提是假设手工标注的训练样本所包含的人脸特征位于其窗口中央。但是人脸特征在图像中的位置都是手工标注的，因此总是有错误的标注（或者说是偏差）。为了克服这个问题，可以将R设计成一个由中间向外逐渐衰减的函数。这里选取二维高斯分布来表达R，即等价于认为手工标注的错误服从高斯分布。下图为检测人脸左眼眼角的特征模版在手工标注的不同训练样本上产生的不同响应图像，右侧为理想响应图像：

2.      如何训练和获取团块模型（stochasticgradient descent）

(1)    为什么要用随机梯度下降算法

    对于上述公式，如果用最小二乘法求解矩阵P，那它的计算量实在太大。因为问题解的数量和团块模型中像素数量的个数一样多，比如若团块模型的大小为40*40，则方程的解（自变量）的数量将有1600个（维），这么大的计算量对于一个实时跟踪系统是不能忍受的。

    比较有效的替代方案是采用随机梯度下降法，该算法将团块模型的解集看成一幅地势图，通过不断迭代获得地势图梯度方向的近似估计，并每次向梯度的反方向前进，直到走到目标函数的极小值（达到阈值或迭代次数上限）。另外，随机梯度法每次随机选取样本，只需要很少的样本就能达到最优解，非常适合实时性要求较高的系统。

上面的三维图仅仅是为了演示，真实的维数为团块模型中像素的数量

(2)    梯度下降法理论及算法步骤

     在讲解随机梯度法如何应用于本节系统前，先看看教科书是如何介绍梯度下降法（最速下降法）。梯度下降法是求解无约束最优化问题的一种最常用方法，它也是一种迭代方法，每一步需要求解目标函数的梯度向量。

    假设f(x)在Rⁿ上是具有一阶连续偏导数的函数，要求解的无约束最优化问题是：

     x^*表示目标函数f(x)的极小点。梯度下降法是一种迭代算法，选取适当的初值x⁽⁰⁾，不断迭代，更新x的值，进行目标函数的极aaa小化，直到收敛。由于负梯度方向是使函数值下降最快的方向，在迭代的每一步，以负方向更新x的值，从而达到减少函数值的目的。

   这里，为f(x)在x^(k)的梯度。

   求出第k+1次迭代值x^(k+1)：

   其中，p_k是搜索方向，取负方向，λ_{k是步长，由一维搜索确定，即λk使得：}

梯度下降法的求解步骤：

实例：利用随机梯度法拟合直线[求方程系数]

     如下图所示，给定m组点的集合R，找到一条曲线或者曲面，对其进行拟合。

    如果我们用分量a,b描述横坐标x_i，则可用以下模型来表达对函数的估计：

   我们建立误差函数，以此评价上面建立的模型对数据的描述是否准确：

  _{为了尽可能的拟合上面曲线，我们希望所选取的参数θ能使得误差函数最小：}

   这里将参数θ的集合看成一个场（地势图），假设随机站在该曲面的一点，要以最快的速度到达最低点，我们当然会沿着坡度最大的方向往下走（梯度的反方向）。用数学描述就是一个求偏导数的过程：

   根据梯度下降算法，参数θ的更新过程如下，a为学习速率：

   下面给出matlab仿真代码：

clc;
clear;

matrix_A = [ 1 4; 2 5; 5 1; 4 2];
matrix_y = [19 26 19 20];
weights = ones(1,2)';
[r,c] = size(matrix_A);


for j = 1:1:2000000  
    weights_old = weights;
    for k = 1:1:r     
        alpha = 4/(j+k+1) + 0.01;                                      %学习速率
        index = floor(random('unif',1,r,1,1));                         %随机选取
        deta_J = matrix_A(index,:) * weights - matrix_y(:,index)';     %偏导数
        weights = weights - alpha * deta_J * matrix_A(index,:)';
    end
    temp = sum((weights-weights_old)*(weights-weights_old)');
    if(sqrt(temp)<0.000001)                                            % 这里是判断收敛的条件，当然可以有其他方法来做  
        break;  
    end  
end

%画图
%原始坐标
plot3(matrix_A(:,1),matrix_A(:,2),matrix_y','r*');
hold on
%函数图像
x=1:1:5;
y=1:1:5;
z=weights(1)*x+weights(2)*y;
plot3(x,y,z);
xlabel('a')
ylabel('b')
zlabel('y')
grid on
axis square

三维图像，视角自己调整，否则看着非常不像是拟合的直线。。。

(3)    梯度上升法与梯度下降法的关系

    其实，二者可以互通，完全是一回事情，修改下目标函数。梯度上升是求最大值的，你把一个本来用梯度下降的函数取负，就变成梯度上升了。下面公式中a为步长。

梯度上升法迭代公式：

梯度上升法介绍参见：http://blog.csdn.net/jinshengtao/article/details/40021139

(4)    利用随机梯度下降法求解团块模型

    回到patch model的定义，我们对下面的函数求偏导数：

   这样P的更新过程与上面例子类似，p=p-a*D，然后就是不断更新P和求D的迭代，直到D满足一定条件退出。

3.      代码实现

（1 ）      输入输出

  团块模型的学习由train函数完成：

void 
patch_model::
train(const vector<Mat> &images,
      const Size psize,
      const float var,
      const float lambda,
      const float mu_init,
      const int nsamples,
      const bool visi)

输入参数含义：

images：包含多个样本图像的矩阵向量（原始含有人像的图像）

psize：团块模型窗口的大小

var：手工标注错误的方差（生成理想图像时使用）

lambda：调整的参数（调整上一次得到的团块模型的大小，以便于当前目标函数偏导数作差）

mu_init：初始步长（构造梯度下降法求团块模型时的更新速率）

nsamples：随机选取的样本数量（梯度下降算法迭代的次数）

visi：训练过程是否可观察标志

输出：

P：得到训练后的团块模型（针对某一个特征的团块模型，并不能描述完整的人脸）

（2 ）      具体实现步骤

。 设置idealresponse map大小：

。设置idealresponse map生成函数：

。利用最大最小值法，对理想反馈图像进行归一化处理：normalize(F,F,0,1,NORM_MINMAX)

b.  随机梯度下降法求最优的团块模型

。 给定初始更新速率：

。随机选取原始样本图像，并将其转换成灰度图，并对灰度图求对数，提高图像在不同光照和对比度情况下的健壮性：

I = this->convert_image(images[i]);

。计算目标函数的偏导数D： 公式同上不罗列了，只是对I有归一化处理。

。更新团块模型P：

c.      代码实现：

求灰度图：

Mat
patch_model::
convert_image(const Mat &im)
{
  Mat I; 
  if(im.channels() == 1)
  {
    if(im.type() != CV_32F)
	   im.convertTo(I,CV_32F); 
    else 
	   I = im;
  }else
  {
    if(im.channels() == 3)
	{
      Mat img; 
	  cvtColor(im,img,CV_RGB2GRAY);
      if(img.type() != CV_32F)
	     img.convertTo(I,CV_32F); 
      else 
	     I = img;
    }
	else
	{
	    cout << "Unsupported image type!" << endl; abort();
    }
  }
  I += 1.0; 
  log(I,I); 
  return I;
}

求团块模型：

void 
patch_model::
train(const vector<Mat> &images,
      const Size psize,
      const float var,
      const float lambda,
      const float mu_init,
      const int nsamples,
      const bool visi)
{
  int N = images.size(),n = psize.width*psize.height;
  int dx = wsize.width-psize.width,dy = wsize.height-psize.height;
  Mat F(dy,dx,CV_32F);
  
  //生成服从高斯分布的理想反馈图像F
  Size wsize = images[0].size();
  if((wsize.width < psize.width) || (wsize.height < psize.height))
  {
    cerr << "Invalid image size < patch size!" << endl; 
	throw std::exception();
  }
 
  for(int y = 0; y < dy; y++)
  {   
    float vy = (dy-1)/2 - y;
    for(int x = 0; x < dx; x++)
	{ 
	  float vx = (dx-1)/2 - x;
      F.fl(y,x) = exp(-0.5*(vx*vx+vy*vy)/var);
    }
  }
  normalize(F,F,0,1,NORM_MINMAX);

  //allocate memory
  Mat I(wsize.height,wsize.width,CV_32F);                  //被选中的样本灰度图像
  Mat dP(psize.height,psize.width,CV_32F);                 //目标函数的偏导数，大小同团块模型
  Mat O = Mat::ones(psize.height,psize.width,CV_32F)/n;    //生成团块模型的归一化模版
  P = Mat::zeros(psize.height,psize.width,CV_32F);         //团块模型

  //利用随机梯度下降法求最优团块模型
  RNG rn(getTickCount());
  double mu=mu_init,step=pow(1e-8/mu_init,1.0/nsamples);
  for(int sample = 0; sample < nsamples; sample++)
  { 
    int i = rn.uniform(0,N);                               // i为随机选中的样本图像标记
    I = this->convert_image(images[i]); 
	dP = 0.0;
    for(int y = 0; y < dy; y++)
	{
      for(int x = 0; x < dx; x++)
	  {
          Mat Wi = I(Rect(x,y,psize.width,psize.height)).clone();
          Wi -= Wi.dot(O); 
		  normalize(Wi,Wi);
          dP += (F.fl(y,x) - P.dot(Wi))*Wi;
      }
    }
    
    P += mu*(dP - lambda*P); mu *= step;
	
    if(visi)                                             //求解过程是否可视化
	{
      Mat R; 
	  matchTemplate(I,P,R,CV_TM_CCOEFF_NORMED);          //利用归一化相关系数匹配法在样本图像上寻找与团块模型匹配的区域
      Mat PP; 
	  normalize(P,PP,0,1,NORM_MINMAX);
      normalize(dP,dP,0,1,NORM_MINMAX);
      normalize(R,R,0,1,NORM_MINMAX);
      imshow("P",PP);                                    //归一化的团块模型
	  imshow("dP",dP);                                   //归一化的目标函数偏导数
	  imshow("R",R);                                     //与团块模型匹配的区域
      if(waitKey(10) == 27)
	     break;
    }
  }
  return;
}

•  考虑全局几何变换的情况

1.      为什么要考虑全局几何变换，它与联合团块模型有什么关系

    目前，我们上面所使用的样本图像，都是基于假设所有的人脸都位于样本图像的中央，并且全局的大小及旋转的角度都已经归一化（人脸中正，没有偏离）。但是，实际人脸跟踪时，人脸可以以任意角度、尺度出现在任何位置。解决训练与测试条件差异导致团块提取失败的方法有两个：一种是在训练阶段手工修改样本图像的缩放尺度、旋转角度到我们所期望的程度，但是由于联合分布团块模型的组成形式比较简单，对于手工设置的样本数据无法产生有效的响应图像response map；另一种方法，由于联合分布团块模型对图像尺度和旋转角度的小范围扰动具有不变性，即前后两帧图像之间人脸的表情变化相对较小，所以我们可以利用上一帧中预估的全局几何变化来约束当前帧人脸的大小尺度和角度，具体做法从产生联合团块模型的图像帧中挑选一个参照帧作为第一个全局几何变化。

   上面是直接翻译的，我认为第一种方法就是在训练阶段提前设置人脸器官尺和角度变化的范围，然后再去训练团块模型，这样做的结果肯定不会很好（因为实际跟踪时情况很复杂，无法假设）；第二种的做法，每次训练团块模型前，利用上一帧的全局几何约束对当前帧做处理，然后再训练联合团块模型，这么做的前提是必须知道第一帧的全局几何约束如何获取。

（第一帧的全局几何约束怎么做的呢？很简单，就是只生成包含尺度因素的人脸模型参考点，并且此参考点矩阵一直是保持不变的，具体解释见下面的reference矩阵）

2.      如何实现

    patch_models类中存储着每个脸部特征的联合团块模型，包括从参考帧中训练的团块模型。人脸跟踪的代码将会直接调用patch_models类中的接口获取人脸特征。（注意是patch_models比上一节多了一个s！）

    如何计算上一帧的全局几何变化？都有以下函数实现：

void 
patch_models::
train(ft_data &data,
      const vector<Point2f> &ref,
      const Size psize,
      const Size ssize,
      const bool mirror,
      const float var,
      const float lambda,
      const float mu_init,
      const int nsamples,
      const bool visi)

入参含义：

data：存放手工标注的数据，包括坐标点集、样本图像名、点的连接关系等等

ref：指定大小和旋转角度的人脸特征的参考点集，即上面的reference

psize：团块模型的大小

sszie：搜索窗口的大小，即在样本图像上可以搜索特征模版（团块模型）的范围，后面有个wsize = psize + ssize，表示标注点附近的图像区域（特征模版是在该图像区域内搜索的）

mirror：是否使用镜像样本数据

再后面的几个参数略，见patch_model的train函数入参说明。

    其中人脸参考点集ref的产生，首先人工指定人脸的宽度width，并设置参数矩阵smodel.p为全0；然后根据联合变化矩阵V中的尺度向量基smodel.V.col(0)的范围min~max，通过缩放计算在该尺度范围内人脸的大小，存入参数矩阵p的第一列smodel.p.fl(0)；最后，通过联合变化矩阵V和参数矩阵，生产人脸参考坐标点集ref。

    这里的V和p的概念可以见上一篇博文，简单概括一下V是2n*(k+4)矩阵，包含n个坐标点的(k+4)种变化，前4列为刚性变化，后k列为非刚性变化（就是k种表情模式），p是(k+4)*1矩阵，将原始坐标投影到展示的空间里去的投影矩阵，包括尺度、旋转、多个表情模式。现在p矩阵除了尺度项，其余全零（可认为旋转角度均为0），那么当V与p做乘法时，则对n个手工标注的特征点仅做尺度和旋转的约束，生成包含n个特征点的人脸参考模型ref。这里一个特征点在后面将对应一个团块模型。

vector<Point2f>  shape_model::calc_shape()
{
  Mat s = V*p; int n = s.rows/2; vector<Point2f> pts;
  for(int i = 0; i < n; i++)pts.push_back(Point2f(s.fl(2*i),s.fl(2*i+1)));
  return pts;
}

float                                      //scaling factor
calc_scale(const Mat &X,                   //scaling basis vector
       const float width)              //width of desired shape
{
  int n = X.rows/2; float xmin = X.at<float>(0),xmax = X.at<float>(0);
  for(int i = 0; i < n; i++){
    xmin = min(xmin,X.at<float>(2*i));
    xmax = max(xmax,X.at<float>(2*i));
  }return width/(xmax-xmin);
}


int width  = parse_face_width(argc,argv);
//generate reference shape ,p = Mat::zeros(e.rows,1,CV_32F);
smodel.p = Scalar::all(0.0); 
smodel.p.fl(0) = calc_scale(smodel.V.col(0),width); 
vector<Point2f> r = smodel.calc_shape();

  本节的目的是增强训练后团块特征的鲁棒性，而上面讲了如何生成人脸模型的参考点集，那么我们如何将参考点集应用到团块模型训练之前的样本图像中？说白了也就是如何通过合理的几何变换让样本图像训练出最好的团块模型。为此，我们首先要求出这么个变换矩阵，然后再作仿射变换。我们可以参考上一篇博文中Procrustes analysis算法计算每幅样本图像中特征点集的几何结构到参考人脸模型的变换矩阵，再利用这个矩阵对样本图像进行仿射变换（缩放、旋转、平移）。具体代码实现如下：

  //reference为人脸特征的参考模型，对样本图像做尺度和角度的约束
  int n = ref.size(); 
  reference = Mat(ref).reshape(1,2*n);
  Size wsize = psize + ssize;                             //wsize归一化的样本图像大小

  //n个特征点将对应n个团块
  patches.resize(n);
  for(int i = 0; i < n; i++)                                //遍历n个特征点
  {
    if(visi)cout << "training patch " << i << "..." << endl;
    vector<Mat> images(0);
    for(int j = 0; j < data.n_images(); j++)               //遍历所有样本图像
	{
      Mat im = data.get_image(j,0); 
      vector<Point2f> p = data.get_points(j,false);        //获取手工标注的样本点
      Mat pt = Mat(p).reshape(1,2*n);
      Mat S = this->calc_simil(pt);                  //计算样本点到参考模型的变化矩阵
Mat A(2,3,CV_32F)；
      A.fl(0,0) = S.fl(0,0);    //构造仿射变换矩阵，A前两列为缩放旋转、最后一列为平移 
	  A.fl(0,1) = S.fl(0,1);
      A.fl(1,0) = S.fl(1,0); 
	  A.fl(1,1) = S.fl(1,1);
      A.fl(0,2) = pt.fl(2*i  ) - (A.fl(0,0) * (wsize.width-1)/2 + A.fl(0,1)*(wsize.height-1)/2);
      A.fl(1,2) = pt.fl(2*i+1) - (A.fl(1,0) * (wsize.width-1)/2 + A.fl(1,1)*(wsize.height-1)/2);
      Mat I; 
	  warpAffine(im,I,A,wsize,INTER_LINEAR+WARP_INVERSE_MAP); //对样本图像进行仿射变化
      images.push_back(I);
    }
    patches[i].train(images,psize,var,lambda,mu_init,nsamples,visi);  //从样本图像中训练团块模型
  }

对单幅样本图像仿射变换矩阵的计算并没有迭代，采用更加简单的方式：

(1)   计算手工标注点的重心

(2)   对手工标注点去中心化

(3)   利用最小二乘法计算样本点到参考坐标的旋转角度（可参考上一篇文章）

Mat patch_models::calc_simil(const Mat &pts)
{
  //compute translation
  int n = pts.rows/2; float mx = 0,my = 0;
  for(int i = 0; i < n; i++){
    mx += pts.fl(2*i); my += pts.fl(2*i+1);
  }  
  Mat p(2*n,1,CV_32F); mx /= n; my /= n;
  for(int i = 0; i < n; i++){
    p.fl(2*i) = pts.fl(2*i) - mx; p.fl(2*i+1) = pts.fl(2*i+1) - my;
  }
  //compute rotation and scale
  float a=0,b=0,c=0;
  for(int i = 0; i < n; i++){
    a += reference.fl(2*i) * reference.fl(2*i  ) + 
      reference.fl(2*i+1) * reference.fl(2*i+1);
    b += reference.fl(2*i) * p.fl(2*i  ) + reference.fl(2*i+1) * p.fl(2*i+1);
    c += reference.fl(2*i) * p.fl(2*i+1) - reference.fl(2*i+1) * p.fl(2*i  );
  }
  b /= a; c /= a;
  float scale = sqrt(b*b+c*c),theta = atan2(c,b); 
  float sc = scale*cos(theta),ss = scale*sin(theta);
  return (Mat_<float>(2,3) << sc,-ss,mx,ss,sc,my);
}

    经过上述变化，我们就有了每个标注点对应的经过仿射变换的图像区域，我们在该图像区域搜索、训练团块模型。下面分别给出每个标注点对应的样本图像区域和训练后的归一化团块模型：

   上图中是7个标注点在对应8幅原始图像上的训练区域（由于设定了大小、又进行了平移旋转缩放等操作，所以效果类似于裁剪出一个团块模型的训练区域）

     上图是第一个标注点对应的团块模型训练过程，需要经过1000次迭代（随机梯度下降法），我们看到经过不到100次迭代其实就已经收敛了。

    下面是讲76个训练得到的团块模型“贴”在手动标注点后的结果，可以直观地看到这些团块模型与标注点的关系。

• 总结

概括一下本章团块特征提取算法都经历了哪些操作：

1.      获取手工标注的样本点、样本图像名称、形状模型(v,e,c)

2.      指定大小与旋转角度，通过形状模型的联合变化矩阵v，生成人脸特征点的参考模型ref

3.      计算每幅样本图像的标注点到参考特征点的旋转角度S

4.      利用旋转构造仿射变化矩阵，对样本图像进行仿射变换A

5.      利用随机梯度下降法对新生成的样本图像求解每个特征点对应的团块模型patch_model

 终于还剩最后一节了，老外的东西还真是麻烦。完整实现代码见留言板链接。

《Master Opencv...读书笔记》非刚性人脸跟踪 III

原文：http://blog.csdn.net/jinshengtao/article/details/43974311