圖像變形可以說是很多圖像、動畫領(lǐng)域的一個非常常見的功能，就說ps、天天P圖、美圖秀秀、可牛等這些每個軟件，有好多個功能都要用到圖像變形，比如圖像方向校正、圖像全景、視頻防抖等，在我的另外一篇博文全景矩形還原，就要用到圖像變形算法。
變形算法實現(xiàn)圖像姿態(tài)調(diào)整
變形算法實現(xiàn)物體方向位置調(diào)整
可以說ps中的一些圖像扭曲都是通過變形方法實現(xiàn)的，比如這篇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是每年必看的文章。保相似的變形算法就相當于是一個經(jīng)典又好用的圖像變形算法，具有保相似性，避免圖像變形過程中，發(fā)生嚴重錯切變換。好了到底是哪篇文獻這么好用，那就是2005 siggraph上的一篇文獻《As-Rigid-As-Possible Shape Manipulation》，這篇文獻在三維網(wǎng)格曲面圖形變形方面也具有重要的作用。
好了廢話不多說，這里開始講解《As-Rigid-As-Possible Shape Manipulation》，因為這篇paper分為兩個步驟，
第一步主要是實現(xiàn)盡量保證變形前后，每個三角形都于原來的對應(yīng)三角形相似，又稱之為網(wǎng)格模型的保相似變形。
第二步主要是為了盡量保證每個三角形與原三角形盡量全等。由于第一步中算法只是保相似，這樣變形后三角形會有的被放大，有的被縮小，出現(xiàn)圖像的局部放大與縮小，所以需要對每個三角形進行適合的縮放調(diào)整，這一步又稱為擬合步，擬合變形后的每個三角形與原網(wǎng)格三角形全等。
一、理論講解
(1)相對坐標的概念。 開始這個算法之前，我需要先給講解一個三角形中相對坐標的概念：給定一個三角形v0v1v2，如下圖所示：
那么我們定義，v2相對于v0、v1的坐標為(x01,y01)，滿足：
也就是說，我們可以通過已知的一個三角形，然后建立一個以三角形的一個頂點為原點，以其鄰接邊、及其對應(yīng)的垂直旋轉(zhuǎn)90度后的方向分別作為x軸、y軸，求取三角形另外一個頂點的相對于這個坐標系的坐標。OK，如果我給你一個三角形，我要你計算出(x01,y01)，你要怎么計算？這個是一個簡單的數(shù)學(xué)問題，同時也是本文算法開始的第一步，具體x01、y01怎么計算，那是高中問題了，不解釋。
(2)整體保相似。 保相似圖像變形算法的原理很簡單，說白了就是盡量保證變形前后，相對坐標不變，也就是說保證每個三角形的相對坐標(x12,y12)、(x20,y20)、(x01,y01)不變。因為如果兩個三角形的這個相對坐標一樣，那么這兩個三角形就是相似三角形。那么算法具體是怎么實現(xiàn)盡量保證，每個三角形變形前后的這三個相對坐標不變的呢？
假設(shè)變形后的三角形頂點坐標為v0'，v1'，v2'。因為我們是要盡量的保證相對坐標不變，我們可以根據(jù)下面的計算公式：
計算出v2desired，然后使得如下誤差最小化：
上面的能量公式越小，就代表變形后的三角形與原三角形越相似。而對于有多個三角面片的網(wǎng)格模型，那么總誤差能量公式為：
這個問題其實歸結(jié)為求解一個最小二乘問題，大家還記得嗎，最小二乘其實就是使得誤差能量最小化，因此到了這里其實我們就可以猜到，這個算法最后歸結(jié)為求解一個稀疏的線性系統(tǒng)。我們的目的是使得上面的總能量最小化，設(shè)相對坐標v'=,那么：
我們的問題是要使得E 1{v’} 最小化，而最小化說白了就是使得其偏導(dǎo)數(shù)為0，因此我們需要對 E 1{v’} 求偏導(dǎo)。開始前先把公式5，換另外一種形式來表達：
其中，q表示控制頂點,u表示自由頂點。對上面的式子求偏導(dǎo)：
把上面的式子簡寫為：
需要注意的是：對于控制頂點不變的模型，G'、B是一個固定的稀疏矩陣，因此我們可以通過上式求解u，得到自由頂點的坐標。
(3)局部完全保相似,及局部縮放調(diào)整。 每個三角形分開重新調(diào)整，使得其與原三角形更為相似，同時對其做縮放調(diào)整。。如果算法直到這里就結(jié)束了，那么變形后的圖像將會出現(xiàn)局部放大與縮小。就拿上面的瘦臉圖片來說，如果只經(jīng)過上面相似變換步驟，我們只是要移動臉型輪廓的特征點，而最后變形出來的結(jié)果可能就會出現(xiàn)眼睛等局部位置的放大，或縮小，因此我們要進行大小調(diào)整，瘦臉后，眼睛的大小與原來圖像的大小相等。
就像上面圖像中，圖 (a) 是原圖，圖(c) 的結(jié)果就是只經(jīng)過相似變換后的結(jié)果，可以看到(c) 中的右手部變粗了許多，而左手部變細了許多，因此我們接著這一步要做的事就是調(diào)整大小，使得(a) 變形后，使得圖像變形后不會出現(xiàn)局部放大與縮小的情況，得到如(d)結(jié)果圖。
另一方面， 在經(jīng)過上面的整體相似變換后，因為我們用了整個三角網(wǎng)格模型進行整體調(diào)整，使得能量最小化，這個就像最小二乘一樣，最后求解出來的結(jié)果并不能保證每個三角形與原三角形完全相似，也會發(fā)生錯切變換，誤差還是存在的。 就像最小二乘，求解出來的結(jié)果并不能保證最后每個方程的殘差為0。
在經(jīng)過(2)中的相似變換后，因為我們是對網(wǎng)格模型的所有三角形一起調(diào)整的，保證總?cè)切五e切變形能最小化，這個時候得到的每個三角形也不全是完全相似的一個三角形，且大小也發(fā)生了變化，如上圖所示。我們得到了v0'v1'v2'，這個三角形我們稱之為中介三角形，接著我們要把每個三角形分開擬合成跟原三角形v0v1v2完全全等的三角形，我們定義擬合階段的三角形變形能為：
當然我們還需要保證變形前后的也與原三角形相似：
這樣聯(lián)立公式(9)、(10)，同時對其求偏導(dǎo)，使得變形能最小化：
記住，這一步是對每個三角形進行分開調(diào)整，這樣才能得到與原網(wǎng)格模型的三角形完全相似。同時在滿是上面式子后，對整個三角形以重心坐標為縮放原點，對每個三角形進行縮放調(diào)整
(4)整體擬合
上面那一部是對每個三角形分開調(diào)整，這個時候，每個三角形是完全獨立的，我們需要對齊對整體連接調(diào)整，否則兩個鄰接三角形的共用邊的兩個頂點的位置各不一樣，這樣根本是不是三角網(wǎng)格模型。
整個三角網(wǎng)格模型的總能量為：
H是一個表示整個網(wǎng)格模型拓撲連接關(guān)系的矩陣，對上式換另外一種表達方式，并對其求偏導(dǎo)，求取最小值：
簡化為：
然后重新求解，可得到自由頂點u的坐標。
二、算法流程
輸入 ：原三角網(wǎng)格模型，控制頂點
輸出 ：變形后的三角網(wǎng)格模型
Algorithm：
1、整體相似變換
(1)計算原網(wǎng)格模型的相對坐標 。根據(jù)給定的原三角網(wǎng)格模型，遍歷每一個三角面片，
并計算每個頂點，相對于另外兩個頂點的坐標。以計算v2相對于邊v0v1的坐標為例：
(x01,y01)表示的是一個相對坐標，如果三角形的三個頂點的坐標已知v0，v1，v2，其實計算這個相對坐標是非常容易的。
(1)建立局部坐標軸。先求出邊向量v0v1，然后利用v0v1求出其單位向量即為局部坐標系的x軸，接著只需要把x軸旋轉(zhuǎn)90度，就可以得到y(tǒng)軸了。
(2)計算相對坐標(x01,y01)，這一步首先需要把向量v0v2投影到x軸，y軸。這里需要注意的是我們上面的公式中x軸指的是v0v1，y軸指的是把v0v1旋轉(zhuǎn)90度后的向量，因此我們需要對投影結(jié)果做一個縮放處理
相關(guān)代碼實現(xiàn)如下：
[cpp] ?view plaincopy
?? for ?(?unsigned? int ?i?=?0;?i?<?nTris;?++i?)?{?? ????Triangle?&?t?=?m_vTriangles[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]?);?? ?? ?????????? ????????Eigen::Vector2d?v01(?v1?-?v0?);?? ????????Eigen::Vector2d?v01N(?v01?);???? ????????v01N.normalize();?? ????????Eigen::Vector2d?v01Rot90(?v01(1),?-v01(0)?);?? ????????Eigen::Vector2d?v01Rot90N(?v01Rot90?);??v01Rot90N.normalize();?? ?? ?????????? ????????Eigen::Vector2d?vLocal(v2?-?v0);?? ????????float ?fX?=?vLocal.dot(v01)?/SquaredLength(v01); ?? ????????float ?fY?=?vLocal.dot(v01Rot90)?/SquaredLength(v01Rot90); ?? ?? ?????????? ????????Eigen::Vector2d?v2test(v0?+?fX?*?v01?+?fY?*?v01Rot90);?? ????????float ?fLength?=?Length(v2test?-?v2);?? ????????t.vTriCoords[j]?=?Eigen::Vector2d(fX,fY);?? ????}?? }??
(2)構(gòu)造公式(7)的相關(guān)矩陣。
[cpp] ?view plaincopy
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、局部完全相似，及大小縮放調(diào)整。
計算公式(11)的F、C矩陣:
[cpp] ?view plaincopy
?? 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;?? ?? ?????? ?? ?????? ????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;?? ?? ?????? ????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;?? ?? ?????? ????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;?? ?? ?????? ????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;?? ?? ?????? ?? ?????? ????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;?? ?? ?????? ????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;?? ?? ?????? ????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;?? ?? ?????? ????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、全局連接擬合。
構(gòu)造公式(16)的相關(guān)矩陣：
[cpp] ?view plaincopy
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的總思路就是分成三步，因為如果控制頂點、和自由頂點的沒有發(fā)生變化，相關(guān)的矩陣只要經(jīng)過第一次求解就可以了，因此可以實現(xiàn)實時拖拽更新，這個我在三維圖形學(xué)網(wǎng)格曲面變形的時候，用到的拉普拉斯網(wǎng)格變形的方法思路大體相同，需要對原網(wǎng)格模型做預(yù)處理操作，然后才能實現(xiàn)實時變形。具體這個算法要和圖像結(jié)合起來，就要對圖像進行三角剖分，然后計算每個三角形的放射變換矩陣，根據(jù)仿射變換矩陣計算每個像素點值的變化，才能對圖像進行變形，具體就不貼代碼了，因為這一些代碼都是比較簡單的，沒有什么技術(shù)含量在里面。本文地址：http://blog.csdn.net/hjimce/article/details/45766321? ? 作者：hjimce ? ? 聯(lián)系qq：1393852684 ?? 更多資源請關(guān)注我的博客：http://blog.csdn.net/hjimce ??? ? ? ? ? ? ? ?原創(chuàng)文章，轉(zhuǎn)載請保留本行信息
總結(jié)
以上是生活随笔 為你收集整理的图像处理（十三）保刚性图像变形算法-Siggraph 2004 的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
如果覺得生活随笔 網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔 推薦給好友。
