Procrustes Analysis普氏分析法

     选取N幅同类目标物体的二维图像，并用上一篇博文的方法标注轮廓点，这样就得到训练样本集：

     由于图像中目标物体的形状和位置存在较大偏差，因此所得到的数据并不具有仿射不变性，需要对其进行归一化处理。这里采用Procrustes分析方法对样本集中的所有形状集合进行归一化。形状和位置的载体还是样本点的空间坐标。

     普氏分析法是一种用来分析形状分布的方法。数学上来讲，就是不断迭代，寻找标准形状(canonical shape)，并利用最小二乘法寻找每个样本形状到这个标准形状的仿射变化方式。（可参照维基百科的GPA算法）

本书中，两个形状的归一化过程（一个形状为canonical shape，另一个为样本形状）：

(1)    求每个样本点i(i=1,2..,n)在N幅图像中的均值

(2)    对所有形状的大小进行归一化，即将每个样本点减去其对应均值

(3)    根据去中心化数据，计算每幅图像中形状的重心，对于第i幅图像，其重心为：

(4)    根据重心和角度，将标准和样本形状对齐在一起，使得两个形状的普氏距离最小，下式为普氏距离定义：

 这个第(4)步的具体做法，不断迭代以下过程：

(a)通过计算每幅图像中所有归一化样本点的平均值得到每个图像的标准形状canonical shape。

(b)利用最小二乘法求每个图像中样本形状到标准形状的旋转角度。根据普氏距离的定义，也就是求：

其中的a和b表示仿射变换里旋转变化的参数：

对上式求偏导数，可以得到所求的a和b：

 (c)根据旋转参数，对样本形状做旋转变化，得到和标准形状对齐的新的形状

 (d)重复以上步骤，直到达到指定循环次数或者前后两次迭代之间canonical shape的绝对范数满足一定阈值

   Procrustesanalysis的作用可以看作是一种对原始数据的预处理，目的是为了获取更好的局部变化模型作为后续模型学习的基础。如下图所示，每一个人脸特征点可以用一种单独的颜色表示；经过归一化变化，人脸的结构越来越明显，即脸部特征簇的位置越来越接近他们的平均位置；经过一系列迭代，尺度和旋转的归一化操作，这些特征簇变得更加紧凑，它们的分布越来越能表达人脸表情的变化。【剔除刚性部分、保留柔性部分】

 

      上图为不同大小、不同长宽比的矩形，经过归一化过程后，各个样本点分布服从一定概率分布趋势。

 

具体代码实现：

函数：Mat shape_model::procrustes(const Mat &X,const int itol,constfloat ftol)

矩阵X：2n*N维矩阵，即有N副图像，每幅图像有n个特征点（2n是因为x和y坐标）

Itol：迭代次数

Ftol：求绝对范数的阈值

返回值：矩阵P，存放对齐后的形状坐标

Mat shape_model::procrustes(const Mat &X,const int itol,const float ftol)
{
  /* X.cols:特征点个数，X.rows: 图像数量*2 */
  int N = X.cols,n = X.rows/2; 

  //remove centre of mass
  Mat P = X.clone();
  for(int i = 0; i < N; i++)
  {
    /*取X第i个列向量*/
    Mat p = P.col(i);
    float mx = 0,my = 0;
    for(int j = 0; j < n; j++)
    {
        mx += p.fl(2*j);
        my += p.fl(2*j+1);
    }
    /*分别求图像集，2维空间坐标x和y的平均值*/
    mx /= n;
    my /= n;
    /*对x,y坐标去中心化*/
    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++)
  {   
    /*计算（含n个形状）的重心*/  
    Mat C = P*Mat::ones(N,1,CV_32F)/N; 
    /*C为2n*1维矩阵，含n个重心，对n个重心归一化处理*/
    normalize(C,C);
    if(iter > 0)
    {
      if(norm(C,C_old) < ftol)//norm:求绝对范数，小于阈值，则退出循环
        break;
    }
    C_old = C.clone();
    for(int i = 0; i < N; i++)
    {
      //求当前形状与归一化重心之间的旋转角度，即上式a和b
      Mat R = this->rot_scale_align(P.col(i),C); 
      for(int j = 0; j < n; j++)
      {
        float x = P.fl(2*j,i),y = P.fl(2*j+1,i);
        /*仿射变化*/
        P.fl(2*j  ,i) = R.fl(0,0)*x + R.fl(0,1)*y;
        P.fl(2*j+1,i) = R.fl(1,0)*x + R.fl(1,1)*y;
      }
    }
  }
return P;
}

http://blog.csdn.net/jinshengtao/article/details/43376049