圖形處理（四）基於梯度場的網格編輯-Siggraph 2004_圖形處理

基於梯度場的網格編輯，對應的Paper為《Mesh Editing with Poisson-Based Gradient Field Manipulation》，是Siggraph 2004上的一篇paper，這篇paper與基於拉普拉斯的網格變形方法，統稱為基於微分域的網格變形演算法，這篇paper其實本質上最後的求解公式和基於拉普拉斯的網格變形方法一樣，之所以能夠siggraph，是因為它通過泊松梯度場的原理進行推導，演算法的巧妙之處在於它以頂點(x,y,z)中的每一維作為一個標量場。

這篇paper涉及到的概念：散度、梯度場、標量場、向量場等看起來很難的東西，說實話，對於這篇paper因為網上找不到原始碼，我把這篇paper看了好多遍，才把它的代碼寫出來。學這篇paper時是我第一次學習向量場的相關知識，向量場在三維演算法中非常重要，同時當時給我的感覺也真不是一般的難，我看了好多關於標量場、向量場的相關知識理論，才感覺慢慢理解。

一、相關理論

數學上的泊松方程：

其中f表示標量場，w表示梯度場。

三角網格曲面上的微分運算元離散化(引用自《勾畫式泊松網格編輯》)：

給定定義在網格曲面上的分段線性標量場f(v)=fi*φi(v),其中v為網格曲面上的任意一點；fi為標量場在網格曲面頂點vi處的函數值；φi(*)為分段線性基函數，它在頂點vi處取值為1 ，在其餘頂點處取值為0。我們有標量場f 對應的梯度運算元

其中▽φi(*)僅在頂點的鄰接三角形上有非零值，且由於φi(*)分段線性，▽φi(*)在各個鄰接三角形上為分段常值函數. 從幾何角度，可以容易地給
出▽φi(*)在三角形T =(vi，vj ，vk)上的定義:

其中，R90代表繞三角形法向量nT 旋轉90度，AT是三角形的面積。類似地，給定定義在三角網格曲面上的分段常值的向量場w，我們定義在頂點vi處w的散度為：

根據梯度運算元和散度運算元的定義，最後可以推匯出網格曲面上的標量場f在頂點vi處的拉普拉斯運算元為：

這篇paper是由浙大的牛人周坤提出來的，演算法最後跟拉普拉斯網格編輯的最後公式可以說是一樣的，然而它的標量場給我很大的啟示，這篇paper直接把(x,y,z)中的x，y，z分別當做一個標量場，然後對標量場求取梯度場，最後求取散度，然後通過泊松方程重建網格模型，實現網格變形。想要更深入的瞭解泊松重建，可以看看我的另外一篇博文《影像處理（十二）映像融合(1)Seamless cloning泊松複製-Siggraph 2004》

二、演算法實現

1、求取源網格曲面的梯度場，最後求取梯度場的散度。

[cpp]  view plain  copy //計算各個頂點的梯度   void CScaleDeformBrush::Get_Faces_Gradient()   {       int fn=m_BaseMesh->faces.size();       m_BaseMesh->need_adjacentfaces();          #pragma omp parallel for       for (int i=0;i<fn;i++)       {           TriMesh::Face &f=m_BaseMesh->faces[i];           vec vij=m_BaseMesh->vertices[f[1]]-m_BaseMesh->vertices[f[0]];           vec vik=m_BaseMesh->vertices[f[2]]-m_BaseMesh->vertices[f[0]];           vec normalf=vij CROSS vik;           float areaf=0.5f*len(normalf);           normalize(normalf);           for (int k=0;k<3;k++)           {              m_Face_Gradient[i][k]=vec(0,0,0);              for (int j=0;j<3;j++)              {               vec ei=m_BaseMesh->vertices[f[(j+2)%3]]-m_BaseMesh->vertices[f[(j+1)%3]];               vec gradient=float(m_BaseMesh->vertices[f[j]][k]*0.5f/areaf)*(normalf CROSS ei);               m_Face_Gradient[i][k]=m_Face_Gradient[i][k]+gradient;               }              }       }      }   void CScaleDeformBrush::Compute_Divergence()   {           //計算頂點的散度       m_BaseMesh->need_adjacentfaces();       int vn=m_BaseMesh->vertices.size();       #pragma omp parallel for       for (int i=0;i<vn;i++)       {              for (int j=0;j<3;j++)           {               m_vertices[i].VDivergence[j]=0.0f;           }           vector<int>&adjacentface=m_BaseMesh->adjacentfaces[i];           for (int j=0;j<adjacentface.size();j++)           {               TriMesh::Face &f=m_BaseMesh->faces[adjacentface[j]];               for (int k=0;k<3;k++)               {                   if (f[k]==i)                   {                       vec ei=m_BaseMesh->vertices[f[(k+2)%3]]-m_BaseMesh->vertices[f[(k+1)%3]];                       vec e1=m_BaseMesh->vertices[f[(k+1)%3]]-m_BaseMesh->vertices[f[k]];                       vec e2=m_BaseMesh->vertices[f[(k+2)%3]]-m_BaseMesh->vertices[f[k]];                       double cot_angle1=Cot_angle(e2,ei);                       double cot_angle2=Cot_angle(-1.0f*e1,ei);                       for (int xyz=0;xyz<3;xyz++)                       {                           m_vertices[i].VDivergence[xyz]+=0.5*(cot_angle1*(e1 DOT m_Face_Gradient[adjacentface[j]][xyz])+cot_angle2*(e2 DOT m_Face_Gradient[adjacentface[j]][xyz]));                       }                       break;                   }               }           }          }      }   //計算v1 v2 之間夾角的餘切值   double CScaleDeformBrush::Cot_angle(vec v1,vec v2)   {       vec vivo=v1;       vec vjvo=v2;       double dotvector=vivo DOT vjvo;       dotvector=dotvector/sqrt(len2(vivo)*len2(vjvo)-dotvector*dotvector);       return dotvector;   }  

2、構建泊松方程的係數，矩陣A，也就是計算拉普拉斯矩陣

[cpp]  view plain  copy //鄰接頂點的餘切權重計算   void CScaleDeformBrush::CotangentWeights(TriMesh*TMesh,int vIndex,vector<double>&vweight,double &WeightSum,bool bNormalize)//計算一階鄰近點的各自cottan權重   {          int NeighborNumber=TMesh->neighbors[vIndex].size();       vweight.resize(NeighborNumber);       WeightSum=0;       vector<int>&NeiV=TMesh->neighbors[vIndex];       for (int i=0;i<NeighborNumber;i++)       {           int j_nei=NeiV[i];           vector<int>tempnei;           Co_neighbor(TMesh,vIndex,j_nei,tempnei);           double cotsum=0.0;           for (int j=0;j<tempnei.size();j++)           {               vec vivo=TMesh->vertices[vIndex]-TMesh->vertices[tempnei[j]];               vec vjvo=TMesh->vertices[j_nei]-TMesh->vertices[tempnei[j]];               double dotvector=vivo DOT vjvo;               dotvector=dotvector/sqrt(len2(vivo)*len2(vjvo)-dotvector*dotvector);               cotsum+=dotvector;           }           vweight[i]=cotsum/2.0;           WeightSum+=vweight[i];       }          if ( bNormalize )        {           for (int k=0;k<NeighborNumber;++k)           {               vweight[k]/=WeightSum;           }           WeightSum=1.0;       }   }      //擷取兩頂點的共同鄰接頂點   void CScaleDeformBrush::Co_neighbor(TriMesh *Tmesh,int u_id,int v_id,vector<int>&co_neiv)   {       Tmesh->need_adjacentedges();       vector<int>&u_id_ae=Tmesh->adjancetedge[u_id];        int en=u_id_ae.size();       Tedge Co_Edge;       for (int i=0;i<en;i++)       {           Tedge &ae=Tmesh->m_edges[u_id_ae[i]];           int opsi=ae.opposite_vertex(u_id);           if (opsi==v_id)           {               Co_Edge=ae;               break;           }       }       for (int i=0;i<Co_Edge.m_adjacent_faces.size();i++)       {           TriMesh::Face af=Tmesh->faces[Co_Edge.m_adjacent_faces[i]];           for (int j=0;j<3;j++)           {               if((af[j]!=u_id)&&(af[j]!=v_id))               {                   co_neiv.push_back(af[j]);               }           }       }   }   //計算拉普拉斯矩陣   void CScaleDeformBrush::Get_Laplace_Matrix()   {       int vn=m_BaseMesh->vertices.size();       int count0=0;       vector<int>begin_N(vn);       for (int i=0;i<vn;i++)       {              begin_N[i]=count0;           count0+=m_BaseMesh->neighbors[i].size()+1;       }       typedef Eigen::Triplet<double> T;       std::vector<T> tripletList(count0);       for(int i=0;i<vn;i++)       {           VProperty & vi = m_vertices[i];           tripletList[begin_N[i]]=T(i,i,-vi.VSumWeight);           int nNbrs = vi.VNeighbors.size();           for (int k = 0;k<nNbrs;++k)            {               tripletList[begin_N[i]+k+1]=T(vi.VNeighbors[k],i,vi.VNeiWeight[k]);           }       }       m_Laplace_Matrix.resize(vn,vn);        m_Laplace_Matrix.setFromTriplets(tripletList.begin(), tripletList.end());      }  
3、添加邊界約束條件，並求解泊松方程，更新變形結果。

即時更新函數：

[cpp]  view plain  copy void CScaleDeformBrush::Update_V_Position()   {       Get_Faces_Gradient();//求梯度       int fn=m_BaseMesh->faces.size();       if(!m_ScaleFace.empty())       for (int i=0;i<fn;i++)       {           if(m_ScaleFace[i])           {               for (int j=0;j<3;j++)               {                   m_Face_Gradient[i][j]=1.1f*m_Face_Gradient[i][j];               }               m_BaseMesh->faces[i].beSelect=false;           }       }       Compute_Divergence();//求散度       if(!m_MatricesCholesky)//構建拉普拉斯矩陣       {           double a=m_Laplace_Matrix.coeff(0,0) +1;           m_Laplace_Matrix.coeffRef(0,0)=a;           m_MatricesCholesky=new Eigen::SimplicialCholesky<SparseMatrixType>(m_Laplace_Matrix);//矩陣分解       }       int vn=m_BaseMesh->vertices.size();       for (int i=0;i<3;i++)       {           Eigen::VectorXd rhs_xyz(vn);           for (int j=0;j<vn;j++)           {               rhs_xyz[j]=m_vertices[j].VDivergence[i];//方程組右邊           }           rhs_xyz[0]=rhs_xyz[0]+1.0f*m_BaseMesh->vertices[0][i];           Eigen::VectorXd xyz=m_MatricesCholesky->solve(rhs_xyz);//求解方程           for (int j=0;j<vn;j++)           {               m_BaseMesh->vertices[j][i]=xyz[j];//更新結果           }       }       m_ScaleFace.clear();       m_ScaleFace.resize(fn,false);       m_BaseMesh->normals.clear();       m_BaseMesh->FaceNormal.clear();      }  

接著我們來看一下用這個演算法實現的簡單局部編輯結果：

上面的即時局部縮放演算法我是通過另外一篇paper《Differential-Based Geometry and Texture Editing with Brushes》的思想實現的，這篇paper基本上就是拷貝《Mesh Editing with Poisson-Based Gradient Field Manipulation》的思想，唯一的創新點在於它的即時互動設計方面，因為我是為了實現即時縮放刷，所以縮放的思想就是根據《Differential-Based Geometry and Texture Editing with Brushes》進行寫代碼的。

上面是用了上面的演算法進行簡單的即時編輯。

這篇paper後面還有後續的演算法調整，比如梯度方向調整、還有實現網格融合、幾何紋理Transfer。其中梯度方向調整是實現保特徵變形的必備條件，因此如果你想要實現完整的演算法，就要對梯度方向進行調整，這個可以參考我的以一篇博文《基於旋轉不變數的網格變形》。在這裡，我就不詳細講方向調整了，方向調整有專門的演算法，paper很多。

須知：基於梯度域的變形方法和拉普拉斯網格變形演算法一樣，微分座標不具有旋轉不變的特點，在變形的時候，會發生曲面細節扭曲，需要對微分座標，或者梯度方向進行調整，才能實現保特徵變形，要實現旋轉不變的變形，可以參考我的另外一篇博文《基於旋轉不變數的網格變形》以此實現旋轉不變的特點。

本文地址：http://blog.csdn.net/hjimce/article/details/46415291    作者：hjimce     聯絡qq：1393852684   更多資源請關注我的部落格：http://blog.csdn.net/hjimce                  原創文章，轉載請保留本行資訊

參考文獻：

1、《Mesh Editing with Poisson-Based Gradient Field Manipulation》

2、微分網格處理技術

3、勾畫式泊松網格編輯