Procrustes 分析

Geometrical constraints (几何约束)

面部几何通常参数化为两个元素的组成：一个全局（刚性）转换和一般局部（非刚性）变形。全局转换解释为图像中人脸的整体布局，通常允许没有约束的改变（即，人脸可以在图像上的任何地方显示）。这包括每张图像人脸的x，y坐标，面内头部旋转，和图像中人脸的大小。另一方面，局部变形解释为面部形状交叉身份之间的不同和表情之间的不同。与全局转换相对照，这些局部变形通常有更多的限制，约束，主要是由于人脸特征的更高地结构化参数。全局转换是2D坐标 的一般函数，适用于任何类型的对象，然而局部变形是特定对象并且必须从训练集中学习。

补充：

1、“变换”（或者转换）和“变形”的区别。变换例如尺度，位置，大小。是控制一个人形状特征在图像中的大体位置。而变形就好比由一个人的面部表情，或者不同人脸的人脸形状内引起的。

2、如果我们知道V是一个空间，那么：

s=V.t()*p;//表示矢量p在V空间下的坐标投影。

p=V*s,表示相当于由坐标（权重）和与之对于的V的列向量，叠加，组合而成的一个形状矢量。

特别的，如果V.t()==V.inv(),即矩阵转置等于矩阵逆，则V为正交矩阵，也即V.t()*V.inv()=E，此时上式就是完全的可逆变换了。

在这一部分，我们描述面部结构几何模型的构造，因此称为形状模型。取决于应用，它可以捕捉单个个体的表情变化，和交叉不同的人脸形状或者两者的混合。这个模型用shape_model类实现，可以在shape_model.hpp和shape_mode.cpp文件中找到。下面的代码片段是shape_model类的头的一部分，用来说明它的基本的功能。

class shape_model{ //2d linear shape model//2维线性形状模型
public:
Mat p; //parameter vector (kx1) CV_32F   //参数矢量
Mat V; //linear subspace (2nxk) CV_32F  //线性子空间
Mat e; //parameter variance (kx1) CV_32F  //方差参数
Mat C; //connectivity (cx2) CV_32S
...
void calc_params(
const vector<Point2f> &pts, //points to compute parameters//用来计算参数的点集
const Mat &weight = Mat(), //weight/point (nx1) CV_32F  //点集的权重
const float c_factor = 3.0); //clamping factor   //钳位因子
...
vector<Point2f> //shape described by parameters //形状的参数化描述
calc_shape();
...
void train(
const vector<vector<Point2f> > &p, //N-example shapes //n个样本形状
const vector<Vec2i> &con = vector<Vec2i>(),//connectivity //连通性
const float frac = 0.95, //fraction of variation to retain  //用来保留的变化分数
const int kmax = 10); //maximum number of modes to retain //用来保留的模型的最大数量
...
}

说明:

1)p,参数矢量p存储着关于模型的形状编码

2)用来表示人脸形状变化的模型用子空间矩阵V和方差矢量e编码

3)连通性矩阵C也存储这个类中，因为仅属于人脸形状的可视化实例。

4)shape_model有三个基本函数：calc_param,calc_shape,和train

calc_param函数:将一个点集投影到可信的人脸形状空间，该函数对投影到该空间的每一个点提供一个单独的权重。

calc_shape函数:该函数使用人脸模型(由V和e编码)解码参数矢量p，产生一个点集。

train函数:该函数从脸形状的数据集中学习编码模型,这些人脸形状都是有相同数量的点组成。

5)参数frac和kmax是训练过程的参数，可是随时为数据指定该参数。

这个类的功能将在下面的部分做详细的阐述，这里我们开始描述一下Procrustes分析，严格地寄存点集的一种方法，接着描述表示局部表形的线性模型。程序在train_shape_model.cpp和visualize_shape_model.cpp文件下，训练和可视化模型性单独进行的，他们的使用将在本部分的结束部分概括。

Procrustes 分析

Procrustes分析的基本思想是最小化所有形状到平均形状的距离和，即最小化公式:

这里的X表示形状，由点集组成的一个矢量，如X=(x0,x1,x2,x3....,xn-1,y0,y1,y2....yn-1)或者交错排列X=(x0,y0,x1,y1,....xn-1,yn-1)这里的0,1,2...n-1为下标。

以下摘自原文：

为了建立人脸形状的变形模型，我们必须首先处理原始标记数据，去除属于全局刚性运动的成分。当用2维进行几何建模时，刚性运动通常表现为类似的变换。这种变换包括尺度，面内旋转和平移。下面的图像阐述了一套在类似变换下的运动类型。从点集中移除全局刚性运动的过程成为Procrustes分析。

从数学上，Procrustes分析的目的是同时找到一个权威的形状和相似变换，该变换用这个权威形状将每个数据实例对齐。这里的对齐是通过变换形状间的最小二乘距离来衡量 的。完成这个目的是个迭代的过程，通过下面的shape_model类来实现。

#define fl at<float>
Mat shape_model::procrustes(
const Mat &X, //interleaved raw shape data as columns //交错原始形状数据，按列排列
const int itol, //maximum number of iterations to try //尝试最大的迭代次数。
const float ftol) //convergence tolerance //收敛误差
{
int N = X.cols,n = X.rows/2; Mat Co,P = X.clone();//copy//复制,
/*----------形状对齐流程步骤1，将所有样本对齐到原点，即将每个样本的x,y坐标各自减去其平均值-----------------*/
for(int i = 0; i < N; i++){
Mat p = P.col(i); //i'th shape//第i个形状
float mx = 0,my = 0; //compute centre of mass...//计算块的中心
for(int j = 0; j < n; j++){ //for x and y separately//对x和y进行单独处理
mx += p.fl(2*j); my += p.fl(2*j+1);
}
mx /= n; my /= n;
for(int j = 0; j < n; j++){ //remove center of mass//减去块的中心
p.fl(2*j) -= mx; p.fl(2*j+1) -= my;
}
}
/*-----------end步骤1-------------------*/
/*-------形状对齐流程步骤3和4，计算变换后的形状的平均值X-,并进行标准化处理------------*/
for(int iter = 0; iter < itol; iter++){ 
Mat C = P*Mat::ones(N,1,CV_32F)/N; //compute normalized...//计算标准化(计算变换后形状的平均值)
normalize(C,C); //canonical shape //权威形状(对x-进行标准化处理)
/*-------------end步骤3和步骤4--------------*/
/*---判断是否收敛----*/
if(iter > 0){if(norm(C,Co) < ftol)break;} //converged?//收敛？当两个权威形状或者成为标准形状的误差小于某一值这里是ftol则，停止迭代。
Co = C.clone(); //remember current estimate//记下当前的误差，进行和下一次进行比较
for(int i = 0; i < N; i++){//对于每一个形状
Mat R = this->rot_scale_align(P.col(i),C);//求两个形状之间的误差满足最小二乘时的旋转矩阵。即相当于两个形状"最靠近"时，需要的旋转的仿射矩阵
for(int j = 0; j < n; j++){ //apply similarity transform//应用相似变换，这对形状中的每一个点，应用仿射矩阵。
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; //returned procrustes aligned shapes//返回procurstes分析对齐的形状
}

其中：

我们要分析该函数的功能，我们首先要知道传入该函数的数据的格式，这里我们主要看一下const Mat  &X参数，我们进行Procrustes分析，是为了满足后面的train训练，我们看一下train函数(shape_model.cpp中，为shape_model类的成员函数)中的procruste函数的调用:

 //vectorize points//矢量化点集
  Mat X = this->pts2mat(points);//这里的points为传入的参数类型vector<vector<Point2f> >
  int N = X.cols,n = X.rows/2;

  //align shapes//对齐形状
  Mat Y = this->procrustes(X);

因此我们知道，它是通过pts2mat函数实现的点击到Mat矩阵的转换,该函数实现如下：

Mat
shape_model::
pts2mat(const vector<vector<Point2f> > &points)
{
  int N = points.size(); assert(N > 0);
  int n = points[0].size();
  for(int i = 1; i < N; i++)assert(int(points[i].size()) == n); //检测是否所有的特征点个数均为76
  Mat X(2*n,N,CV_32F);
  for(int i = 0; i < N; i++){
    Mat x = X.col(i),y = Mat(points[i]).reshape(1,2*n); y.copyTo(x);
  }return X;
}

说明：

1、X为2*n×N的矩阵，这里的n为特征点的个数76，N为图像的个数，或者说样本，实例的个数。

2、vector<Point2f> 类型的对象作为参数传入Mat构造函数形成一个临时Mat对象，该对象调用reshape函数。故实现为(x0,y0,x1,y1,....xn-1,yn-1)T 的交错的列矢量。

procrustes函数

算法开始减去每个形状块的中心，接着在计算权威形状，所有形状均值归一化和旋转、尺度化到最佳匹配权威形状中迭代处理。估计权威形状的标准化步骤有必要固定尺度的问题，防止将所有形状的尺度缩小为0。这个固定尺度的选择是任意的，这里我们标准化权威矢量的为1.0，这是OpenCV的normalize函数的默认处理。计算面内选择和尺度化，将每一个形状对齐到当前估计的权威形状，通过函数rot_scale_align函数来完成。

Mat shape_model::rot_scale_align(
const Mat &src, //[x1;y1;...;xn;yn] vector of source shape //原形状矢量
const Mat &dst) //destination shape//目的形状
{
//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);
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;//solve linear system//解决线性系统
return (Mat_<float>(2,2) << a,-b,b,a);
}

这个函数最小化下面的待旋转形状和权威形状的最小二乘不同。数学上表示如下：

这里最小二乘问题的解决采用闭合形式的求解，展示在下面图像等式的右手边。既然不是求解尺度化和面内旋转，这两个操作在尺度化2维旋转矩阵是非线性关系，我们求解变量(a,b)。这两个变量与尺度旋转矩阵有关：

Procureste对原始注释形状数据分析的可视化效果在下面的图像中阐述。每一个面部特征用唯一的颜色显示。转换标准化之后，人脸的结构开始显现出来，这里人脸特征簇的位置围绕他们均值的位置。迭代尺度和旋转标准化过程之后，特征簇变的更加紧凑并且他们的分布变的由面部变形导致的变化更具有代表性。最后一点是重要的，因为这些变形我们将在接下来的部分尝试模型化。因此，Procruster分析的角色可以被认为是对原始数据的操作，这个操作将得到更好的人脸局部变形模型用来学习。

（以此为未对齐形状，平移对齐，Procrustes对齐)