本文中，我们通过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;

    //所有的形状向量(特征)对齐到原点，即每个向量的分量减去其平均值，每列是一个形状向量。
    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;
            }
        }

//优化缩放和旋转
    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) //当两个标准形状或者标准形状的误差小于某一值，这里是ftol则，停止迭代。如果达到最大迭代次数，也退出迭代。
            {

norm函数默认是矩阵各元素平方和的根范式
            if(norm(C,C_old) < ftol)
                break;
            }
        C_old = C.clone(); //记下当前的矩阵，和下一次进行比较
        for(int i = 0; i < N; i++)
            {

              Mat R = this->rot_scale_align(P.col(i),C); //求两个形状之间的误差满足最小二乘时的旋转矩阵。即相当于两个形状"最靠近"时，需要的旋转的仿射矩阵

rot_scale_align函数求每副图像的特征点向量和平均向量满足最小二乘法时候的旋转矩阵。即求得a,b值组成的旋转缩放矩阵(注意：这儿表示旋转和缩放矩阵的组合，便于线性表示)。

        rot_scale_align函数的代码如下，原理如下面的两个图，即当前特征向量中的每个特征点经过旋转缩放变化的值和标准形状(权威形状)对应的特征点值差的平方和最小，结果得到一个2*2的旋转缩放矩阵。

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);
    }

            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;
    }