圖形處理（七）基於熱傳播的測地距離計算-Siggraph 2013_圖形處理

這裡要跟大家分享的是2013年Siggraph上面的一篇paper，名為《Geodesics in Heat:A New Approach to Computing Distance Based on Heat Flow》，這篇paper沒有提供原始碼，但是因為演算法的思想相當新穎，如果你之前有研究過其它的測地三角網格曲面上的測地距離演算法，那麼看到這篇paper後，你會非常的激動，覺得這個演算法相當神奇，網格曲面上測地距離的計算方法又有了新的突破。因為看到這篇paper非常激動，以至於我興奮得馬上去把代碼寫了一遍，雖然測地距離的計算方法，對於我來說沒什麼用，但它的想法，它的思想值得我們學習。

《Geodesics in Heat:A New Approach to Computing Distance Based on Heat Flow》演算法主要是通過向量場的方法，通過求解熱流，求解泊松方程。這篇paper我覺得應該把它歸類為向量場類型的文獻，曲面上向量場的應用非常廣泛，可以用於網格最佳化重建、參數化紋理映射、網格變形……等，如今這篇paper又讓我明白向量場也可以求解測地距離。通過熱量傳播的方法，去求解測地距離，真是大牛啊

一、理論知識

在很早之前對於測地距離，大牛們就推匯出了測地距離的求解歸結為求解eikonal equation(程函方程)：

且滿足邊界約束條件為：

上面符號Φ就是測地距離。

也就是不管是歐式空間還是曲面空間，其測地距離的梯度模長恒等於1，然後邊界條件：源點的測地距離值為0。

因此很多大牛，都是針對上面的程函方程的求解進行研究，然而所提出的演算法都是非線性方法，因為上面的方程就是非線性方程，在網格曲面上計算量非常大。而這篇paper作者的真正貢獻在於把非線性問題轉換為線性問題，到最後只需要求解一個泊松方程就可以了。開始這個演算法之前，建議先看一下《Mesh Editing with Poisson-Based Gradient Field Manipulation》這篇利用泊松方程進行三角網格模型編輯，其實這篇paper的思想應該是借用了《Poisson Image Editing》，如果你已經很熟悉《Mesh Editing with Poisson-Based Gradient Field Manipulation》那麼看這篇paper將事半功倍。

1、時間離散化，熱傳播方程時間離散化為：

2、空間離散化，在網格曲面上，離散的拉普拉斯座標定義為：

Ai表示三角面片的面積的三份之一，j表示頂點i的鄰接頂點。對於有V個頂點的網格模型，我們可以列出上面的V個方程，寫出矩陣形式為：

其中，Lc為拉普拉矩陣，A為包含每個頂點面積的對角矩陣。

代入公式(3)可得：

在網格曲面上，對於一個給定的三角面片平面上，標量場的梯度計算公式如下：

Af表示三角形的面積，N表示三角面片的法矢，ui就是標量場，我們可以把網格曲面每個頂點的測地距離看成是一個標量場。

頂點i的散度離散形式為：

懶得多說廢話了，總之到最後歸結為求解如下方程：

其中b我們需要先通過求解標量場u，然後根據散度的計算公式，計算出頂點測地距離的散度。

二、演算法實現

演算法總流程：

示意圖：

其中，Φ就是我們要求的測地距離，t是一個無窮小的數，趨近於0

演算法具體實現，先說明一下，我下面自己寫的代碼很亂，懶得整理，因為我只是為了學習：

1、求解u

根據公式：

求解u，因此我們需要先算出A，t，Lc，還有δ。接著我將逐漸講解著四個參數的求解：

a、時間t計算：見文獻3.2.4，時間理論上來說是一個非常小的數，文中作者給出了其合適的值為：網格平均邊長的平方

代碼實現如下：

[cpp]  view plain  copy //計算時間步長，文獻中時間步長為：網格模型的邊長平均，然後平方   void CHeatGeodesics::Set_Time()   {       m_BaseMesh->need_edge();       int en=m_BaseMesh->m_edges.size();       float sumlength=0;       for (int i=0;i<en;i++)       {           sumlength+=m_BaseMesh->m_edges[i].length();       }       sumlength=sumlength/en;//邊長的平均長度       sumlength=sumlength*sumlength;       m_Time_Step=sumlength;   }   b、面積對角矩陣A的計算：

[cpp]  view plain  copy //計算Geodesics in Heat 文獻中的包含頂點Vi的面積矩陣   void CHeatGeodesics::Get_Matrix_A_areas()   {       m_BaseMesh->need_adjacentfaces();       gsi::SparseMatrix &A=m_A_areas;       A.Resize(m_NumberV,m_NumberV);       //計算每個三角面片的面積       int fn=m_BaseMesh->faces.size();       vector<float>&face_area=m_Faces_Area;       face_area.resize(fn);       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]];           float areaf=0.5*len(vij CROSS vik);           face_area[i]=areaf;       }       //計算拉普拉斯運算元中每個頂點所佔據的面積，即鄰接三角面片面積總和的三分之一       for (int i=0;i<m_NumberV;i++)       {           vector<int>&af=m_BaseMesh->adjacentfaces[i];           int n_af=af.size();           float sumarea=0.0;           for (int j=0;j<n_af;j++)           {               sumarea+=face_area[af[j]];           }           //包含頂點的面積為：鄰接三角面片面積總和的三分之一           sumarea=sumarea/3.0;           A(i,i) =sumarea;       }   }   C、拉普拉斯矩陣Lc計算：

[cpp]  view plain  copy //計算Geodesics in Heat 文獻中的Lc矩陣，即拉普拉斯運算元   void CHeatGeodesics::Get_Matrix_Lc()   {       Ls.resize(m_NumberV,m_NumberV);       int m_nEdges=10000 ;       Ls.reserve(m_nEdges+m_NumberV);       for (int i = 0;i<m_NumberV;++i)        {           VProperty & vi = m_vertices[i];           int nNbrs = vi.VNeighbors.size();           for (int k = 0;k<nNbrs;++k)            {               Ls.insert(i, vi.VNeighbors[k]) = vi.VNeiWeight[k];           }           Ls.insert(i, i) = -vi.VSumWeight;       }       Ls.finalize();       gsi::SparseMatrix &A=m_Lc;       A.Resize(m_NumberV,m_NumberV);       for(int k=0;k<m_NumberV;++k)        {           for ( SparseMatrixType::InnerIterator it(Ls,k); it; ++it)           {               A.Set( it.row(), it.col(), it.value() );           }          }      }  

鄰接頂點的餘切cot權重計算：

[cpp]  view plain  copy //鄰接頂點的餘切權重計算   void CHeatGeodesics::CotangentWeights(TriMesh*TMesh,int vIndex,vector<float>&vweight,float &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);           float 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]];               float 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 CHeatGeodesics:: UniformWeights(TriMesh*TMesh,int vIndex,vector<float>&vweight,float &WeightSum,bool bNormalize)   {       int NeighborNumber=TMesh->neighbors[vIndex].size();       vweight.resize(NeighborNumber);       WeightSum = 0;       for (int j = 0; j <NeighborNumber; ++j )       {           vweight[j] = 1;           WeightSum += vweight[j];       }          if ( bNormalize )       {           for ( int k = 0; k < NeighborNumber; ++k )               vweight[k] /= WeightSum;           WeightSum=1.0;       }   }   //擷取兩頂點的共同鄰接頂點   void CHeatGeodesics::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=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]);               }           }       }   }  

最後求解方程組，就可以把u求出來了。

2、求解三角網格曲面的熱量場▽u (Heat flow)：

這一步直接根據公式：

求解就可以了。然後對▽u進行歸一化，並取熱量場的反方向，即求測地距離的梯度場：

[cpp]  view plain  copy for (int i=0;i<fn;i++)   {       TriMesh::Face &f=m_BaseMesh->faces[i];       for (int j=0;j<3;j++)       {           vec ei=m_BaseMesh->vertices[f[(j+2)%3]]-m_BaseMesh->vertices[f[(j+1)%3]];           vec FNormal=normalize(m_BaseMesh->FaceNormal[i]);           vec gradient=float(m_vertices[f[j]].VU*0.5/m_Faces_Area[i])*(FNormal CROSS ei);           FGradient[i]=FGradient[i]+gradient;       }       normalize(FGradient[i]);//文獻的重要兩步 即梯度單位化和梯度反向       FGradient[i]=float(-1.0)*FGradient[i];       //length0[i]=len(FGradient[i]);   }  
3、對測地距離的梯度場X，求取散度。

根據公式：

求解測地距離標量場梯度場的散度。

[cpp]  view plain  copy //計算頂點的散度   for (int i=0;i<vn;i++)   {          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 gradient=float(0.5/m_Faces_Area[adjacentface[j]])*(m_BaseMesh->FaceNormal[adjacentface[j]] CROSS ei);                  //計算散度                  m_vertices[i].VDivergence+=(FGradient[adjacentface[j]] DOT gradient)*m_Faces_Area[adjacentface[j]];*/                      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]];                   float cot_angle1=Cot_angle(e2,ei);                   float cot_angle2=Cot_angle(float(-1.0)*e1,ei);               m_vertices[i].VDivergence+=0.5*(cot_angle1*(e1 DOT FGradient[adjacentface[j]])+cot_angle2*(e2 DOT FGradient[adjacentface[j]]));                   break;               }           }       }      }   4、最後求解通過求解泊松方程，求取網格曲面上每個頂點的測地距離。

[cpp]  view plain  copy gsi::SparseLinearSystem *pSystemPos2 = new gsi::SparseLinearSystem();       gsi::Solver_TAUCS * pSolverPos2 = new gsi::Solver_TAUCS(pSystemPos2);       pSystemPos2->Resize(m_NumberV, m_NumberV);       pSystemPos2->ResizeRHS(1);       //由於以下解方程組使用Solver_TAUCS::TAUCS_LLT ，即半正定矩陣模式       m_Lc.Multiply(-1.0);          m_Lc(0,0) =m_Lc(0,0) +10;          pSystemPos2->SetMatrix(m_Lc);       if ( ! pSystemPos2->Matrix().IsSymmetric() )assert(0);       pSolverPos2->OnMatrixChanged();       pSolverPos2->SetStoreFactorization(true);       pSolverPos2->SetSolverMode( gsi::Solver_TAUCS::TAUCS_LLT );       pSolverPos2->SetOrderingMode( gsi::Solver_TAUCS::TAUCS_METIS );       gsi::Vector pRHSPos2 ;       pRHSPos2.Resize(m_NumberV);       // 右邊項約束添加       for (int i=0;i<m_NumberV;i++)       { &