保刚性图像变形算法-Siggraph 2004

图像处理（十三）保刚性图像变形算法-Siggraph 2004

保刚性图像变形算法

原文地址：http://blog.csdn.net/hjimce/article/details/45766321 

作者：hjimce

一、相关理论

图像变形可以说是很多图像、动画领域的一个非常常见的功能，就说ps、天天P图、美图秀秀、可牛等这些每个软件，有好多个功能都要用到图像变形，比如图像方向校正、图像全景、视频防抖等，在我的另外一篇博文全景矩形还原，就要用到图像变形算法。

变形算法实现图像姿态调整

变形算法实现物体方向位置调整

可以说ps中的一些图像扭曲都是通过变形方法实现的，比如这篇paper:《As-Rigid-As-Possible Image Registration for Hand-drawn Cartoon Animations》就有讲到它的变形算法与ps的效果分析，然而我今天要讲的算法无论是在视屏防抖《Content-Preserving Warps for 3D Video Stabilization》《Bundled Camera Paths for Video Stabilization》，还是在图像全景《Rectangling Panoramic Images via Warping》等都用到这个算法，这几篇都是siggraph上的文献，对于搞图形图像的人来说，siggraph上的paper是每年必看的文章。保相似的变形算法就相当于是一个经典又好用的图像变形算法，具有保相似性，避免图像变形过程中，发生严重错切变换。好了到底是哪篇文献这么好用，那就是2005 siggraph上的一篇文献《As-Rigid-As-Possible Shape Manipulation》，这篇文献在三维网格曲面图形变形方面也具有重要的作用。

好了废话不多说，这里开始讲解《As-Rigid-As-Possible Shape Manipulation》，因为这篇paper分为两个步骤，

第一步主要是实现尽量保证变形前后，每个三角形都于原来的对应三角形相似，又称之为网格模型的保相似变形。

第二步主要是为了尽量保证每个三角形与原三角形尽量全等。由于第一步中算法只是保相似，这样变形后三角形会有的被放大，有的被缩小，出现图像的局部放大与缩小，所以需要对每个三角形进行适合的缩放调整，这一步又称为拟合步，拟合变形后的每个三角形与原网格三角形全等。

一、理论讲解

(1)相对坐标的概念。开始这个算法之前，我需要先给讲解一个三角形中相对坐标的概念：给定一个三角形v0v1v2，如下图所示：

那么我们定义，v2相对于v0、v1的坐标为(x01,y01)，满足：

也就是说，我们可以通过已知的一个三角形，然后建立一个以三角形的一个顶点为原点，以其邻接边、及其对应的垂直旋转90度后的方向分别作为x轴、y轴，求取三角形另外一个顶点的相对于这个坐标系的坐标。OK，如果我给你一个三角形，我要你计算出(x01,y01)，你要怎么计算？这个是一个简单的数学问题，同时也是本文算法开始的第一步，具体x01、y01怎么计算，那是高中问题了，不解释。

(2)整体保相似。保相似图像变形算法的原理很简单，说白了就是尽量保证变形前后，相对坐标不变，也就是说保证每个三角形的相对坐标(x12,y12)、(x20,y20)、(x01,y01)不变。因为如果两个三角形的这个相对坐标一样，那么这两个三角形就是相似三角形。那么算法具体是怎么实现尽量保证，每个三角形变形前后的这三个相对坐标不变的呢？

假设变形后的三角形顶点坐标为v0’，v1’，v2’。因为我们是要尽量的保证相对坐标不变，我们可以根据下面的计算公式：

计算出v2desired，然后使得如下误差最小化：

上面的能量公式越小，就代表变形后的三角形与原三角形越相似。而对于有多个三角面片的网格模型，那么总误差能量公式为：

这个问题其实归结为求解一个最小二乘问题，大家还记得吗，最小二乘其实就是使得误差能量最小化，因此到了这里其实我们就可以猜到，这个算法最后归结为求解一个稀疏的线性系统。我们的目的是使得上面的总能量最小化，设相对坐标v’=,那么：

我们的问题是要使得E1{v’}最小化，而最小化说白了就是使得其偏导数为0，因此我们需要对E1{v’}求偏导。开始前先把公式5，换另外一种形式来表达：

其中，q表示控制顶点,u表示自由顶点。对上面的式子求偏导：

把上面的式子简写为：

需要注意的是：对于控制顶点不变的模型，G’、B是一个固定的稀疏矩阵，因此我们可以通过上式求解u，得到自由顶点的坐标。

(3)局部完全保相似,及局部缩放调整。每个三角形分开重新调整，使得其与原三角形更为相似，同时对其做缩放调整。

首先，在经过上面的整体相似变换后，三角形可能进行放大或者缩小，因为我们用了三角形的内在属性相对坐标，这样只能保证所有三角形的尽量相似，并不能保证全等。如果算法直到这里就结束了，那么变形后的图像将会出现局部放大与缩小。就拿上面的瘦脸图片来说，如果只经过上面相似变换步骤，我们只是要移动脸型轮廓的特征点，而最后变形出来的结果可能就会出现眼睛等局部位置的放大，或缩小，因此我们要进行大小调整，瘦脸后，眼睛的大小与原来图像的大小相等。

就像上面图像中，图(a)是原图，图(c)的结果就是只经过相似变换后的结果，可以看到(c)中的右手部变粗了许多，而左手部变细了许多，因此我们接着这一步要做的事就是调整大小，使得(a)变形后，使得图像变形后不会出现局部放大与缩小的情况，得到如(d)结果图。

另一方面，在经过上面的整体相似变换后，因为我们用了整个三角网格模型进行整体调整，使得能量最小化，这个就像最小二乘一样，最后求解出来的结果并不能保证每个三角形与原三角形完全相似，也会发生错切变换，误差还是存在的。就像最小二乘，求解出来的结果并不能保证最后每个方程的残差为0。

在经过(2)中的相似变换后，因为我们是对网格模型的所有三角形一起调整的，保证总三角形错切变形能最小化，这个时候得到的每个三角形也不全是完全相似的一个三角形，且大小也发生了变化，如上图所示。我们得到了v0’v1’v2’，这个三角形我们称之为中介三角形，接着我们要把每个三角形分开拟合成跟原三角形v0v1v2完全全等的三角形，我们定义拟合阶段的三角形变形能为：

当然我们还需要保证变形前后的也与原三角形相似：

这样联立公式(9)、(10)，同时对其求偏导，使得变形能最小化：

记住，这一步是对每个三角形进行分开调整，这样才能得到与原网格模型的三角形完全相似。同时在满是上面式子后，对整个三角形以重心坐标为缩放原点，对每个三角形进行缩放调整

(4)整体拟合

上面那一部是对每个三角形分开调整，这个时候，每个三角形是完全独立的，我们需要对齐对整体连接调整，否则两个邻接三角形的共用边的两个顶点的位置各不一样，这样根本是不是三角网格模型。

整个三角网格模型的总能量为：

H是一个表示整个网格模型拓扑连接关系的矩阵，对上式换另外一种表达方式，并对其求偏导，求取最小值：

简化为：

然后重新求解，可得到自由顶点u的坐标。

二、算法流程

输入：原三角网格模型，控制顶点

输出：变形后的三角网格模型

Algorithm：

1、整体相似变换

(1)计算原网格模型的相对坐标。根据给定的原三角网格模型，遍历每一个三角面片，

并计算每个顶点，相对于另外两个顶点的坐标。以计算v2相对于边v0v1的坐标为例：

(x01,y01)表示的是一个相对坐标，如果三角形的三个顶点的坐标已知v0，v1，v2，其实计算这个相对坐标是非常容易的。

(1)建立局部坐标轴。先求出边向量v0v1，然后利用v0v1求出其单位向量即为局部坐标系的x轴，接着只需要把x轴旋转90度，就可以得到y轴了。

(2)计算相对坐标(x01,y01)，这一步首先需要把向量v0v2投影到x轴，y轴。这里需要注意的是我们上面的公式中x轴指的是v0v1，y轴指的是把v0v1旋转90度后的向量，因此我们需要对投影结果做一个缩放处理

相关代码实现如下：

[cpp] view plain copy

print ?

1. // 遍历网格模型的每个三角面片  
2. for ( unsigned int i = 0; i < nTris; ++i ) {  
3.     Triangle & t = m_vTriangles[i];//第i个三角面片  
4.     //计算当前三角形每个顶点相对于另外两个顶点的坐标  
5.     for ( int j = 0; j < 3; ++j ) {  
6.         unsigned int n0 = j;  
7.         unsigned int n1 = (j+1)%3;  
8.         unsigned int n2 = (j+2)%3;//三角面片三个顶点的索引  
9.   
10.         Eigen::Vector2d v0 = GetInitialVert( t.nVerts[n0] );  
11.         Eigen::Vector2d v1 = GetInitialVert( t.nVerts[n1] );  
12.         Eigen::Vector2d v2 = GetInitialVert( t.nVerts[n2] );//三角面片三个顶点的坐标  
13.   
14.         //建立局部坐标系的 x轴 与y 轴  
15.         Eigen::Vector2d v01( v1 - v0 );  
16.         Eigen::Vector2d v01N( v01 );    
17.         v01N.normalize();//边向量v0v1的方向向量   x轴  
18.         Eigen::Vector2d v01Rot90( v01(1), -v01(0) );  
19.         Eigen::Vector2d v01Rot90N( v01Rot90 );  v01Rot90N.normalize();//垂直于边v0v1的方向向量    y轴  
20.   
21.         //计算局部相对坐标  
22.         Eigen::Vector2d vLocal(v2 - v0);  
23.         float fX = vLocal.dot(v01) /SquaredLength(v01);//计算局部相对坐标x01  
24.         float fY = vLocal.dot(v01Rot90) /SquaredLength(v01Rot90);//计算局部相对坐标y01  
25.   
26.           
27.         Eigen::Vector2d v2test(v0 + fX * v01 + fY * v01Rot90);  
28.         float fLength = Length(v2test - v2);  
29.         t.vTriCoords[j] = Eigen::Vector2d(fX,fY);//局部相对坐标(x01,y01)  
30.     }  
31. }  

 // 遍历网格模型的每个三角面片
    for ( unsigned int i = 0; i < nTris; ++i ) {
        Triangle & t = m_vTriangles[i];//第i个三角面片
        //计算当前三角形每个顶点相对于另外两个顶点的坐标
        for ( int j = 0; j < 3; ++j ) {
            unsigned int n0 = j;
            unsigned int n1 = (j+1)%3;
            unsigned int n2 = (j+2)%3;//三角面片三个顶点的索引

            Eigen::Vector2d v0 = GetInitialVert( t.nVerts[n0] );
            Eigen::Vector2d v1 = GetInitialVert( t.nVerts[n1] );
            Eigen::Vector2d v2 = GetInitialVert( t.nVerts[n2] );//三角面片三个顶点的坐标

            //建立局部坐标系的 x轴 与y 轴
            Eigen::Vector2d v01( v1 - v0 );
            Eigen::Vector2d v01N( v01 );  
            v01N.normalize();//边向量v0v1的方向向量   x轴
            Eigen::Vector2d v01Rot90( v01(1), -v01(0) );
            Eigen::Vector2d v01Rot90N( v01Rot90 );  v01Rot90N.normalize();//垂直于边v0v1的方向向量    y轴

            //计算局部相对坐标
            Eigen::Vector2d vLocal(v2 - v0);
            float fX = vLocal.dot(v01) /SquaredLength(v01);//计算局部相对坐标x01
            float fY = vLocal.dot(v01Rot90) /SquaredLength(v01Rot90);//计算局部相对坐标y01


            Eigen::Vector2d v2test(v0 + fX * v01 + fY * v01Rot90);
            float fLength = Length(v2test - v2);
            t.vTriCoords[j] = Eigen::Vector2d(fX,fY);//局部相对坐标(x01,y01)
        }
    }

(2)构造公式(7)的相关矩阵。

[cpp] view plain copy

print ?

1. void RigidMeshDeformer2D::PrecomputeOrientationMatrix()  
2. {  
3.   
4.     std::vector<Constraint,Eigen::aligned_allocator<Constraint> > vConstraintsVec;  
5.     std::set<Constraint>::iterator cur(m_vConstraints.begin()), end(m_vConstraints.end());  
6.     while ( cur != end )  
7.         vConstraintsVec.push_back( *cur++ );  
8.   
9.   
10.     unsigned int nVerts = (unsigned int)m_vDeformedVerts.size();  
11.     m_mFirstMatrix.resize( 2*nVerts, 2*nVerts);  
12.     for ( unsigned int i = 0; i < 2*nVerts; ++i )  
13.         for ( unsigned int j = 0; j < 2*nVerts; ++j )  
14.             m_mFirstMatrix(i,j) = 0.0;  
15.   
16.     size_t nConstraints = vConstraintsVec.size();  
17.     unsigned int nFreeVerts = nVerts - nConstraints;  
18.   
19.   
20.     unsigned int nRow = 0;  
21.     m_vVertexMap.resize(nVerts);  
22.     for ( unsigned int i = 0; i < nVerts; ++i ) {  
23.         Constraint c(i, vec2(0,0));  
24.         if ( m_vConstraints.find(c) != m_vConstraints.end() )  
25.             continue;  
26.         m_vVertexMap[i] = nRow++;  
27.     }  
28.     for ( unsigned int i = 0 ; i < nConstraints; ++i )  
29.         m_vVertexMap[vConstraintsVec[i].nVertex ] = nRow++;  
30.   
31.   
32.   
33.   
34.   
35.   
36.     size_t nTriangles = m_vTriangles.size();  
37.     for ( unsigned int i = 0; i < nTriangles; ++i ) {  
38.         Triangle & t = m_vTriangles[i];  
39.           
40.         double fTriSumErr = 0;  
41.         for ( int j = 0; j < 3; ++j ) {  
42.             double fTriErr = 0;  
43.   
44.             int n0x = 2 * m_vVertexMap[ t.nVerts[j] ];  
45.             int n0y = n0x + 1;  
46.             int n1x = 2 * m_vVertexMap[ t.nVerts[(j+1)%3] ];  
47.             int n1y = n1x + 1;  
48.             int n2x = 2 * m_vVertexMap[ t.nVerts[(j+2)%3] ];  
49.             int n2y = n2x + 1;  
50.             float x = t.vTriCoords[j](0);  
51.             float y = t.vTriCoords[j](1);  
52.   
53.   
54.   
55.   
56.             m_mFirstMatrix(n0x,n0x)+= 1 - 2*x + x*x + y*y;  
57.             m_mFirstMatrix(n0x,n1x)+= 2*x - 2*x*x - 2*y*y;    
58.             m_mFirstMatrix(n0x,n1y)+= 2*y;        
59.             m_mFirstMatrix(n0x,n2x)+= -2 + 2*x;   
60.             m_mFirstMatrix(n0x,n2y)+= -2 * y;     
61.   
62.   
63.   
64.           
65.             m_mFirstMatrix(n0y,n0y)+= 1 - 2*x + x*x + y*y;  
66.             m_mFirstMatrix(n0y,n1x)+= -2*y;           
67.             m_mFirstMatrix(n0y,n1y)+= 2*x - 2*x*x - 2*y*y;  
68.             m_mFirstMatrix(n0y,n2x)+= 2*y;        
69.             m_mFirstMatrix(n0y,n2y)+= -2 + 2*x;   
70.   
71.   
72.               
73.             m_mFirstMatrix(n1x,n1x)+= x*x + y*y;  
74.             m_mFirstMatrix(n1x,n2x)+= -2*x;       
75.             m_mFirstMatrix(n1x,n2y)+= 2*y;        
76.   
77.   
78.   
79.             m_mFirstMatrix(n1y,n1y)+= x*x + y*y;  
80.             m_mFirstMatrix(n1y,n2x)+= -2*y;               
81.             m_mFirstMatrix(n1y,n2y)+= -2*x;               
82.   
83.   
84.   
85.       
86.             m_mFirstMatrix(n2x,n2x)+= 1;  
87.             m_mFirstMatrix(n2y,n2y)+= 1;  
88.   
89.         }  
90.     }  
91.   
92.   
93.   
94.   
95.     Eigen::MatrixXd mG00( 2*nFreeVerts, 2*nFreeVerts );  
96.     ExtractSubMatrix( m_mFirstMatrix, 0, 0, mG00 );  
97.   
98.       
99.     Eigen::MatrixXd mG01( 2 * (int)nFreeVerts, 2 * (int)nConstraints );  
100.     ExtractSubMatrix( m_mFirstMatrix, 0, 2*nFreeVerts, mG01 );  
101.     Eigen::MatrixXd mG10( 2 * (int)nConstraints, 2 * (int)nFreeVerts );  
102.     ExtractSubMatrix( m_mFirstMatrix, 2*nFreeVerts, 0, mG10 );  
103.   
104.     //   
105.     Eigen::MatrixXd mGPrime = mG00 + mG00.transpose();  
106.     Eigen::MatrixXd mB = mG01 + mG10.transpose();  
107.   
108.       
109.     Eigen::MatrixXd mGPrimeInverse=mGPrime.inverse();  
110.   
111.   
112.     Eigen::MatrixXd mFinal = mGPrimeInverse * mB;  
113.     mFinal *= -1;  
114.   
115.     m_mFirstMatrix = mFinal;      
116. }  

void RigidMeshDeformer2D::PrecomputeOrientationMatrix()
{

    std::vector<Constraint,Eigen::aligned_allocator<Constraint> > vConstraintsVec;
    std::set<Constraint>::iterator cur(m_vConstraints.begin()), end(m_vConstraints.end());
    while ( cur != end )
        vConstraintsVec.push_back( *cur++ );


    unsigned int nVerts = (unsigned int)m_vDeformedVerts.size();
    m_mFirstMatrix.resize( 2*nVerts, 2*nVerts);
    for ( unsigned int i = 0; i < 2*nVerts; ++i )
        for ( unsigned int j = 0; j < 2*nVerts; ++j )
            m_mFirstMatrix(i,j) = 0.0;

    size_t nConstraints = vConstraintsVec.size();
    unsigned int nFreeVerts = nVerts - nConstraints;


    unsigned int nRow = 0;
    m_vVertexMap.resize(nVerts);
    for ( unsigned int i = 0; i < nVerts; ++i ) {
        Constraint c(i, vec2(0,0));
        if ( m_vConstraints.find(c) != m_vConstraints.end() )
            continue;
        m_vVertexMap[i] = nRow++;
    }
    for ( unsigned int i = 0 ; i < nConstraints; ++i )
        m_vVertexMap[vConstraintsVec[i].nVertex ] = nRow++;






    size_t nTriangles = m_vTriangles.size();
    for ( unsigned int i = 0; i < nTriangles; ++i ) {
        Triangle & t = m_vTriangles[i];

        double fTriSumErr = 0;
        for ( int j = 0; j < 3; ++j ) {
            double fTriErr = 0;

            int n0x = 2 * m_vVertexMap[ t.nVerts[j] ];
            int n0y = n0x + 1;
            int n1x = 2 * m_vVertexMap[ t.nVerts[(j+1)%3] ];
            int n1y = n1x + 1;
            int n2x = 2 * m_vVertexMap[ t.nVerts[(j+2)%3] ];
            int n2y = n2x + 1;
            float x = t.vTriCoords[j](0);
            float y = t.vTriCoords[j](1);




            m_mFirstMatrix(n0x,n0x)+= 1 - 2*x + x*x + y*y;
            m_mFirstMatrix(n0x,n1x)+= 2*x - 2*x*x - 2*y*y;  
            m_mFirstMatrix(n0x,n1y)+= 2*y;      
            m_mFirstMatrix(n0x,n2x)+= -2 + 2*x; 
            m_mFirstMatrix(n0x,n2y)+= -2 * y;   




            m_mFirstMatrix(n0y,n0y)+= 1 - 2*x + x*x + y*y;
            m_mFirstMatrix(n0y,n1x)+= -2*y;         
            m_mFirstMatrix(n0y,n1y)+= 2*x - 2*x*x - 2*y*y;
            m_mFirstMatrix(n0y,n2x)+= 2*y;      
            m_mFirstMatrix(n0y,n2y)+= -2 + 2*x; 



            m_mFirstMatrix(n1x,n1x)+= x*x + y*y;
            m_mFirstMatrix(n1x,n2x)+= -2*x;     
            m_mFirstMatrix(n1x,n2y)+= 2*y;      



            m_mFirstMatrix(n1y,n1y)+= x*x + y*y;
            m_mFirstMatrix(n1y,n2x)+= -2*y;             
            m_mFirstMatrix(n1y,n2y)+= -2*x;             




            m_mFirstMatrix(n2x,n2x)+= 1;
            m_mFirstMatrix(n2y,n2y)+= 1;

        }
    }




    Eigen::MatrixXd mG00( 2*nFreeVerts, 2*nFreeVerts );
    ExtractSubMatrix( m_mFirstMatrix, 0, 0, mG00 );


    Eigen::MatrixXd mG01( 2 * (int)nFreeVerts, 2 * (int)nConstraints );
    ExtractSubMatrix( m_mFirstMatrix, 0, 2*nFreeVerts, mG01 );
    Eigen::MatrixXd mG10( 2 * (int)nConstraints, 2 * (int)nFreeVerts );
    ExtractSubMatrix( m_mFirstMatrix, 2*nFreeVerts, 0, mG10 );

    // 
    Eigen::MatrixXd mGPrime = mG00 + mG00.transpose();
    Eigen::MatrixXd mB = mG01 + mG10.transpose();


    Eigen::MatrixXd mGPrimeInverse=mGPrime.inverse();


    Eigen::MatrixXd mFinal = mGPrimeInverse * mB;
    mFinal *= -1;

    m_mFirstMatrix = mFinal;    
}

2、局部完全相似，及大小缩放调整。

计算公式(11)的F、C矩阵:

[cpp] view plain copy

print ?

1. //计算公式11 F、C矩阵  
2. void RigidMeshDeformer2D::PrecomputeScalingMatrices( unsigned int nTriangle )  
3. {  
4.      //遍历每个三角形  
5.     Triangle & t = m_vTriangles[nTriangle];  
6.   
7.       
8.     t.mF =Eigen::MatrixXd(4,4);  
9.     t.mC = Eigen::MatrixXd(4,6);  
10.   
11.       
12.     double x01 = t.vTriCoords[0](0);  
13.     double y01 = t.vTriCoords[0](1);  
14.     double x12 = t.vTriCoords[1](0);  
15.     double y12 = t.vTriCoords[1](1);  
16.     double x20 = t.vTriCoords[2](0);  
17.     double y20 = t.vTriCoords[2](1);  
18.   
19.     double k1 = x12*y01 + (-1 + x01)*y12;  
20.     double k2 = -x12 + x01*x12 - y01*y12;  
21.     double k3 = -y01 + x20*y01 + x01*y20;  
22.     double k4 = -y01 + x01*y01 + x01*y20;  
23.     double k5 = -x01 + x01*x20 - y01*y20 ;  
24.   
25.     double a = -1 + x01;  
26.     double a1 = pow(-1 + x01,2) + pow(y01,2);  
27.     double a2 = pow(x01,2) + pow(y01,2);  
28.     double b =  -1 + x20;  
29.     double b1 = pow(-1 + x20,2) + pow(y20,2);  
30.     double c2 = pow(x12,2) + pow(y12,2);  
31.   
32.     double r1 = 1 + 2*a*x12 + a1*pow(x12,2) - 2*y01*y12 + a1*pow(y12,2);  
33.     double r2 = -(b*x01) - b1*pow(x01,2) + y01*(-(b1*y01) + y20);  
34.     double r3 = -(a*x12) - a1*pow(x12,2) + y12*(y01 - a1*y12);  
35.     double r5 = a*x01 + pow(y01,2);  
36.     double r6 = -(b*y01) - x01*y20;  
37.     double r7 = 1 + 2*b*x01 + b1*pow(x01,2) + b1*pow(y01,2) - 2*y01*y20;  
38.   
39.     //计算F矩阵  
40.   
41.     // F矩阵的第一行  
42.     t.mF(0,0)= 2*a1 + 2*a1*c2 + 2*r7;  
43.     t.mF(0,1)= 0;  
44.     t.mF(0,2)= 2*r2 + 2*r3 - 2*r5;  
45.     t.mF(0,3)= 2*k1 + 2*r6 + 2*y01;  
46.   
47.     // F矩阵的第二行  
48.     t.mF(1,0)= 0;  
49.     t.mF(1,1)= 2*a1 + 2*a1*c2 + 2*r7;  
50.     t.mF(1,2)= -2*k1 + 2*k3 - 2*y01;  
51.     t.mF(1,3)= 2*r2 + 2*r3 - 2*r5;  
52.   
53.     // F矩阵的第三行  
54.     t.mF(2,0)= 2*r2 + 2*r3 - 2*r5;  
55.     t.mF(2,1)= -2*k1 + 2*k3 - 2*y01;  
56.     t.mF(2,2)= 2*a2 + 2*a2*b1 + 2*r1;  
57.     t.mF(2,3)= 0;  
58.   
59.     //F矩阵的第四行  
60.     t.mF(3,0)= 2*k1 - 2*k3 + 2*y01;  
61.     t.mF(3,1)= 2*r2 + 2*r3 - 2*r5;  
62.     t.mF(3,2)= 0;  
63.     t.mF(3,3)= 2*a2 + 2*a2*b1 + 2*r1;  
64.   
65.     // 求逆  
66.     Eigen::MatrixXd mFInverse=t.mF.inverse();  
67.     mFInverse *= -1.0;  
68.     t.mF =  mFInverse;  
69.   
70.     //计算C矩阵  
71.   
72.     // C矩阵的第一行  
73.     t.mC(0,0)= 2*k2;  
74.     t.mC(0,1)= -2*k1;  
75.     t.mC(0,2)= 2*(-1-k5);  
76.     t.mC(0,3)= 2*k3;  
77.     t.mC(0,4)= 2*a;  
78.     t.mC(0,5)= -2*y01;  
79.   
80.     // C矩阵的第二行  
81.     t.mC(1,0)= 2*k1;  
82.     t.mC(1,1)= 2*k2;  
83.     t.mC(1,2) = -2*k3;  
84.     t.mC(1,3)= 2*(-1-k5);  
85.     t.mC(1,4)= 2*y01;  
86.     t.mC(1,5)= 2*a;  
87.   
88.     // C矩阵的第三行  
89.     t.mC(2,0)= 2*(-1-k2);  
90.     t.mC(2,1)= 2*k1;  
91.     t.mC(2,2)= 2*k5;  
92.     t.mC(2,3)= 2*r6;  
93.     t.mC(2,4)= -2*x01;  
94.     t.mC(2,5)= 2*y01;  
95.   
96.     //C矩阵的第四行  
97.     t.mC(3,0)= 2*k1;  
98.     t.mC(3,1)= 2*(-1-k2);  
99.     t.mC(3,2)= -2*k3;  
100.     t.mC(3,3)= 2*k5;  
101.     t.mC(3,4)= -2*y01;  
102.     t.mC(3,5)= -2*x01;  
103.   
104.   
105. }  

//计算公式11 F、C矩阵
void RigidMeshDeformer2D::PrecomputeScalingMatrices( unsigned int nTriangle )
{
     //遍历每个三角形
    Triangle & t = m_vTriangles[nTriangle];


    t.mF =Eigen::MatrixXd(4,4);
    t.mC = Eigen::MatrixXd(4,6);


    double x01 = t.vTriCoords[0](0);
    double y01 = t.vTriCoords[0](1);
    double x12 = t.vTriCoords[1](0);
    double y12 = t.vTriCoords[1](1);
    double x20 = t.vTriCoords[2](0);
    double y20 = t.vTriCoords[2](1);

    double k1 = x12*y01 + (-1 + x01)*y12;
    double k2 = -x12 + x01*x12 - y01*y12;
    double k3 = -y01 + x20*y01 + x01*y20;
    double k4 = -y01 + x01*y01 + x01*y20;
    double k5 = -x01 + x01*x20 - y01*y20 ;

    double a = -1 + x01;
    double a1 = pow(-1 + x01,2) + pow(y01,2);
    double a2 = pow(x01,2) + pow(y01,2);
    double b =  -1 + x20;
    double b1 = pow(-1 + x20,2) + pow(y20,2);
    double c2 = pow(x12,2) + pow(y12,2);

    double r1 = 1 + 2*a*x12 + a1*pow(x12,2) - 2*y01*y12 + a1*pow(y12,2);
    double r2 = -(b*x01) - b1*pow(x01,2) + y01*(-(b1*y01) + y20);
    double r3 = -(a*x12) - a1*pow(x12,2) + y12*(y01 - a1*y12);
    double r5 = a*x01 + pow(y01,2);
    double r6 = -(b*y01) - x01*y20;
    double r7 = 1 + 2*b*x01 + b1*pow(x01,2) + b1*pow(y01,2) - 2*y01*y20;

    //计算F矩阵

    // F矩阵的第一行
    t.mF(0,0)= 2*a1 + 2*a1*c2 + 2*r7;
    t.mF(0,1)= 0;
    t.mF(0,2)= 2*r2 + 2*r3 - 2*r5;
    t.mF(0,3)= 2*k1 + 2*r6 + 2*y01;

    // F矩阵的第二行
    t.mF(1,0)= 0;
    t.mF(1,1)= 2*a1 + 2*a1*c2 + 2*r7;
    t.mF(1,2)= -2*k1 + 2*k3 - 2*y01;
    t.mF(1,3)= 2*r2 + 2*r3 - 2*r5;

    // F矩阵的第三行
    t.mF(2,0)= 2*r2 + 2*r3 - 2*r5;
    t.mF(2,1)= -2*k1 + 2*k3 - 2*y01;
    t.mF(2,2)= 2*a2 + 2*a2*b1 + 2*r1;
    t.mF(2,3)= 0;

    //F矩阵的第四行
    t.mF(3,0)= 2*k1 - 2*k3 + 2*y01;
    t.mF(3,1)= 2*r2 + 2*r3 - 2*r5;
    t.mF(3,2)= 0;
    t.mF(3,3)= 2*a2 + 2*a2*b1 + 2*r1;

    // 求逆
    Eigen::MatrixXd mFInverse=t.mF.inverse();
    mFInverse *= -1.0;
    t.mF =  mFInverse;

    //计算C矩阵

    // C矩阵的第一行
    t.mC(0,0)= 2*k2;
    t.mC(0,1)= -2*k1;
    t.mC(0,2)= 2*(-1-k5);
    t.mC(0,3)= 2*k3;
    t.mC(0,4)= 2*a;
    t.mC(0,5)= -2*y01;

    // C矩阵的第二行
    t.mC(1,0)= 2*k1;
    t.mC(1,1)= 2*k2;
    t.mC(1,2) = -2*k3;
    t.mC(1,3)= 2*(-1-k5);
    t.mC(1,4)= 2*y01;
    t.mC(1,5)= 2*a;

    // C矩阵的第三行
    t.mC(2,0)= 2*(-1-k2);
    t.mC(2,1)= 2*k1;
    t.mC(2,2)= 2*k5;
    t.mC(2,3)= 2*r6;
    t.mC(2,4)= -2*x01;
    t.mC(2,5)= 2*y01;

    //C矩阵的第四行
    t.mC(3,0)= 2*k1;
    t.mC(3,1)= 2*(-1-k2);
    t.mC(3,2)= -2*k3;
    t.mC(3,3)= 2*k5;
    t.mC(3,4)= -2*y01;
    t.mC(3,5)= -2*x01;


}

并求解公式11，得每个三角面片新顶点的位置

3、全局连接拟合。

构造公式(16)的相关矩阵：

[cpp] view plain copy

print ?

1. void RigidMeshDeformer2D::PrecomputeFittingMatrices()  
2. {  
3.   
4.     std::vector<Constraint,Eigen::aligned_allocator<Constraint> > vConstraintsVec;  
5.     std::set<Constraint>::iterator cur(m_vConstraints.begin()), end(m_vConstraints.end());  
6.     while ( cur != end )  
7.         vConstraintsVec.push_back( *cur++ );  
8.   
9.       
10.     unsigned int nVerts = (unsigned int)m_vDeformedVerts.size();  
11.     size_t nConstraints = vConstraintsVec.size();  
12.     unsigned int nFreeVerts = nVerts - nConstraints;  
13.   
14.       
15.     unsigned int nRow = 0;  
16.     m_vVertexMap.resize(nVerts);  
17.     for ( unsigned int i = 0; i < nVerts; ++i )   
18.     {  
19.         Constraint c(i,vec2(0,0));  
20.         if ( m_vConstraints.find(c) != m_vConstraints.end() )  
21.             continue;  
22.         m_vVertexMap[i] = nRow++;  
23.     }  
24.     for ( unsigned int i = 0 ; i < nConstraints; ++i )  
25.         m_vVertexMap[vConstraintsVec[i].nVertex ] = nRow++;  
26.       
27.   
28.   
29.     Eigen::VectorXd gUTestX( nVerts ), gUTestY(nVerts);  
30.     for ( unsigned int i = 0; i < nVerts; ++i )   
31.     {  
32.         Constraint c(i,vec2(0,0));  
33.         if ( m_vConstraints.find(c) != m_vConstraints.end() )  
34.             continue;  
35.         int nRow = m_vVertexMap[i];  
36.         gUTestX[nRow] = m_vInitialVerts[i].vPosition(0);  
37.         gUTestY[nRow] = m_vInitialVerts[i].vPosition(1);  
38.     }  
39.     for ( unsigned int i = 0; i < nConstraints; ++i ) {  
40.         int nRow = m_vVertexMap[ vConstraintsVec[i].nVertex ];  
41.         gUTestX[nRow] = vConstraintsVec[i].vConstrainedPos[0];  
42.         gUTestY[nRow] = vConstraintsVec[i].vConstrainedPos[1];  
43.     }  
44.   
45.   
46.   
47.     Eigen::MatrixXd mHX( nVerts, nVerts );  
48.     Eigen::MatrixXd mHY( nVerts, nVerts );  
49.     for ( unsigned int i = 0; i < nVerts; ++i )  
50.         for ( unsigned int j = 0; j < nVerts; ++j )  
51.             mHX(i,j) = mHY(i,j) = 0.0;  
52.   
53.       
54.     size_t nTriangles = m_vTriangles.size();  
55.     for ( unsigned int i = 0; i < nTriangles; ++i ) {  
56.         Triangle & t = m_vTriangles[i];  
57.         double fTriSumErr = 0;  
58.         for ( int j = 0; j < 3; ++j ) {  
59.             double fTriErr = 0;  
60.   
61.             int nA = m_vVertexMap[ t.nVerts[j] ];  
62.             int nB = m_vVertexMap[ t.nVerts[(j+1)%3] ];  
63.   
64.       
65.             mHX(nA,nA)+= 2;  
66.             mHX(nA,nB)+= -2;  
67.             mHX(nB,nA)+= -2;  
68.             mHX(nB,nB)+= 2;  
69.   
70.       
71.             mHY(nA,nA)+= 2;  
72.             mHY(nA,nB)+= -2;  
73.             mHY(nB,nA)+= -2;  
74.             mHY(nB,nB)+= 2;           
75.         }  
76.     }  
77.   
78.   
79.     Eigen::MatrixXd mHX00( (int)nFreeVerts, (int)nFreeVerts );  
80.     Eigen::MatrixXd mHY00( (int)nFreeVerts, (int)nFreeVerts );  
81.     ExtractSubMatrix( mHX, 0, 0, mHX00 );  
82.     ExtractSubMatrix( mHY, 0, 0, mHY00 );  
83.   
84.   
85.     Eigen::MatrixXd mHX01( (int)nFreeVerts, (int)nConstraints );  
86.     Eigen::MatrixXd mHX10( (int)nConstraints, (int)nFreeVerts );  
87.     ExtractSubMatrix( mHX, 0, nFreeVerts, mHX01 );  
88.     ExtractSubMatrix( mHX, nFreeVerts, 0, mHX10 );  
89.   
90.       
91.     Eigen::MatrixXd mHY01( (int)nFreeVerts, (int)nConstraints );  
92.     Eigen::MatrixXd mHY10( (int)nConstraints, (int)nFreeVerts );  
93.     ExtractSubMatrix( mHY, 0, nFreeVerts, mHY01 );  
94.     ExtractSubMatrix( mHY, nFreeVerts, 0, mHY10 );  
95.   
96.   
97.     m_mHXPrime = mHX00;  
98.     m_mHYPrime = mHY00;  
99.   
100.     m_mDX = mHX01;  
101.     m_mDY = mHY01;  
102.   
103.   
104.     m_mLUDecompX=Eigen::FullPivLU<Eigen::MatrixXd>(m_mHXPrime);  
105.     m_mLUDecompY=Eigen::FullPivLU<Eigen::MatrixXd>(m_mHYPrime);  
106. }  

void RigidMeshDeformer2D::PrecomputeFittingMatrices()
{

    std::vector<Constraint,Eigen::aligned_allocator<Constraint> > vConstraintsVec;
    std::set<Constraint>::iterator cur(m_vConstraints.begin()), end(m_vConstraints.end());
    while ( cur != end )
        vConstraintsVec.push_back( *cur++ );


    unsigned int nVerts = (unsigned int)m_vDeformedVerts.size();
    size_t nConstraints = vConstraintsVec.size();
    unsigned int nFreeVerts = nVerts - nConstraints;


    unsigned int nRow = 0;
    m_vVertexMap.resize(nVerts);
    for ( unsigned int i = 0; i < nVerts; ++i ) 
    {
        Constraint c(i,vec2(0,0));
        if ( m_vConstraints.find(c) != m_vConstraints.end() )
            continue;
        m_vVertexMap[i] = nRow++;
    }
    for ( unsigned int i = 0 ; i < nConstraints; ++i )
        m_vVertexMap[vConstraintsVec[i].nVertex ] = nRow++;



    Eigen::VectorXd gUTestX( nVerts ), gUTestY(nVerts);
    for ( unsigned int i = 0; i < nVerts; ++i ) 
    {
        Constraint c(i,vec2(0,0));
        if ( m_vConstraints.find(c) != m_vConstraints.end() )
            continue;
        int nRow = m_vVertexMap[i];
        gUTestX[nRow] = m_vInitialVerts[i].vPosition(0);
        gUTestY[nRow] = m_vInitialVerts[i].vPosition(1);
    }
    for ( unsigned int i = 0; i < nConstraints; ++i ) {
        int nRow = m_vVertexMap[ vConstraintsVec[i].nVertex ];
        gUTestX[nRow] = vConstraintsVec[i].vConstrainedPos[0];
        gUTestY[nRow] = vConstraintsVec[i].vConstrainedPos[1];
    }



    Eigen::MatrixXd mHX( nVerts, nVerts );
    Eigen::MatrixXd mHY( nVerts, nVerts );
    for ( unsigned int i = 0; i < nVerts; ++i )
        for ( unsigned int j = 0; j < nVerts; ++j )
            mHX(i,j) = mHY(i,j) = 0.0;


    size_t nTriangles = m_vTriangles.size();
    for ( unsigned int i = 0; i < nTriangles; ++i ) {
        Triangle & t = m_vTriangles[i];
        double fTriSumErr = 0;
        for ( int j = 0; j < 3; ++j ) {
            double fTriErr = 0;

            int nA = m_vVertexMap[ t.nVerts[j] ];
            int nB = m_vVertexMap[ t.nVerts[(j+1)%3] ];


            mHX(nA,nA)+= 2;
            mHX(nA,nB)+= -2;
            mHX(nB,nA)+= -2;
            mHX(nB,nB)+= 2;


            mHY(nA,nA)+= 2;
            mHY(nA,nB)+= -2;
            mHY(nB,nA)+= -2;
            mHY(nB,nB)+= 2;         
        }
    }


    Eigen::MatrixXd mHX00( (int)nFreeVerts, (int)nFreeVerts );
    Eigen::MatrixXd mHY00( (int)nFreeVerts, (int)nFreeVerts );
    ExtractSubMatrix( mHX, 0, 0, mHX00 );
    ExtractSubMatrix( mHY, 0, 0, mHY00 );


    Eigen::MatrixXd mHX01( (int)nFreeVerts, (int)nConstraints );
    Eigen::MatrixXd mHX10( (int)nConstraints, (int)nFreeVerts );
    ExtractSubMatrix( mHX, 0, nFreeVerts, mHX01 );
    ExtractSubMatrix( mHX, nFreeVerts, 0, mHX10 );


    Eigen::MatrixXd mHY01( (int)nFreeVerts, (int)nConstraints );
    Eigen::MatrixXd mHY10( (int)nConstraints, (int)nFreeVerts );
    ExtractSubMatrix( mHY, 0, nFreeVerts, mHY01 );
    ExtractSubMatrix( mHY, nFreeVerts, 0, mHY10 );


    m_mHXPrime = mHX00;
    m_mHYPrime = mHY00;

    m_mDX = mHX01;
    m_mDY = mHY01;


    m_mLUDecompX=Eigen::FullPivLU<Eigen::MatrixXd>(m_mHXPrime);
    m_mLUDecompY=Eigen::FullPivLU<Eigen::MatrixXd>(m_mHYPrime);
}

这篇paper的总思路就是分成三步，因为如果控制顶点、和自由顶点的没有发生变化，相关的矩阵只要经过第一次求解就可以了，因此可以实现实时拖拽更新，这个我在三维图形学网格曲面变形的时候，用到的拉普拉斯网格变形的方法思路大体相同，需要对原网格模型做预处理操作，然后才能实现实时变形。具体这个算法要和图像结合起来，就要对图像进行三角剖分，然后计算每个三角形的放射变换矩阵，根据仿射变换矩阵计算每个像素点值的变化，才能对图像进行变形，具体就不贴代码了，因为这一些代码都是比较简单的，没有什么技术含量在里面。

参考文献

1、http://www-ui.is.s.u-tokyo.ac.jp/~takeo/research/rigid/

2、《As-Rigid-As-Possible Shape Manipulation》

************* 作者：hjimce     联系qq：1393852684   更多资源请关注我的博客：http://blog.csdn.net/hjimce                  原创文章，转载请保留本行信息*****************

</article>