圖像處理（一）圖像變形(1)矩形全景圖像還原-Siggraph 2014
作者：cvvision
最近發現，看過的文章，沒幾天就忘了，于是開始寫點東西記錄一下，所學習過的算法。廢話不多說，今天看了這篇文獻“Rectangling Panoramic Images via Warping”，是以做一下記錄。這篇文獻算法分為兩步：
1、第一步需要通過計算圖像最小能量線，對圖像進行預變性。
看懂這篇文章首先要看懂Seam Carving，這個算法作者有源碼，可以下載下傳下來好好解讀一下算法原理。通過這一步可以把圖像放大為矩形圖像，這一步又稱之為局部變形。在這一步每一次 疊代的過程中，源圖像上的每個像素點都會有可能移動位置，假設input圖像的源像素點I(x,y)經過Seam Carving算法放大後，位置移動到I'(x,y)，那麼我們可以計算得其運動場（光流場），u(x,y)=I‘(x,y)-I(x,y)。算法第一步 的目的便是要計算出input圖像變為local warping圖像後，每個像素點的u(x,y)運動向量。
2、第二步對local warping圖像矩形進行網格劃分
記錄劃分後，每個格子頂點的坐标v(x,y)，ok，通過第一步我們知道了每個像素點經過u(x,y)運動場，變形後才有了local warping 圖像，接着我們把圖像local warping通過
-u(x,y)運動場，變形回去，同時記錄下v(x,y)每個頂點的新位置v’(x,y)，其實這一步不需要對圖像的資料進行變形，隻需要定劃分的網格根據-u(x,y)進行直接變形，記錄下每個網格頂點的新位置。
3、接着我們以這些網格頂點新位置v’(x,y)為控制頂點，利用相似變換的變形方法，對最開始的input 圖像進行全局變形。變形的控制頂點為v’(x,y)，變形後的位置為v(x,y)，據此可以得到全局的圖像變形結果
學習這篇文獻，需要先好好學習：“Seam Carving for Content-Aware Image Resizing”和“As-rigid-as-possible_shape_manipulation”這兩篇文獻網上可以下載下傳到作者寫的源代碼，所 以比較容易學會，搞定這兩篇paper，這篇文獻就簡單多了。
http://www.cvvision.cn/5493.html
著作權歸作者所有。
商業轉載請聯系作者獲得授權，非商業轉載請注明出處。
連結：http://www.cvvision.cn/5493.html
來源：CV視覺網
一、相關理論圖像變形可以說是很多圖像、動畫領域的一個非常常見的功能，就說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的坐标。
二、算法流程
輸入：原三角網格模型，控制頂點
輸出：變形後的三角網格模型
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度後的向量，是以我們需要對投影結果做一個縮放處理
相關代碼實作如下：
// 周遊網格模型的每個三角面片
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)的相關矩陣。
void RigidMeshDeformer2D::PrecomputeOrientationMatrix()
{
std::vector<constraint,eigen::aligned_allocator > vConstraintsVec;
std::set::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矩陣:
//計算公式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)的相關矩陣：
void RigidMeshDeformer2D::PrecomputeFittingMatrices()
for ( unsigned int i = 0; i < nVerts; ++i )
{
Constraint c(i,vec2(0,0));
Eigen::VectorXd gUTestX( nVerts ), gUTestY(nVerts);
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 j = 0; j < nVerts; ++j )
mHX(i,j) = mHY(i,j) = 0.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(m_mHXPrime);
m_mLUDecompY=Eigen::FullPivLU(m_mHYPrime);
這篇paper的總思路就是分成三步，因為如果控制頂點、和自由頂點的沒有發生變化， 相關的矩陣隻要經過第一次求解就可以了，是以可以實作實時拖拽更新，這個我在三維圖形學網格曲面變形的時候，用到的拉普拉斯網格變形的方法思路大體相同， 需要對原網格模型做預處理操作，然後才能實作實時變形。具體這個算法要和圖像結合起來，就要對圖像進行三角剖分，然後計算每個三角形的放射變換矩陣，根據 仿射變換矩陣計算每個像素點值的變化，才能對圖像進行變形，具體就不貼代碼了，因為這一些代碼都是比較簡單的，沒有什麼技術含量在裡面。
參考文獻
1、http://www-ui.is.s.u-tokyo.ac.jp/~takeo/research/rigid/
2、《As-Rigid-As-Possible Shape Manipulation》
![Ubuntu20.04下配置OpenCV2.4.13[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)