這里要跟大家分享的是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 plaincopy
?? 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 plaincopy
?? 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 plaincopy
?? 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 plaincopy
?? void ?CHeatGeodesics::CotangentWeights(TriMesh*TMesh, int ?vIndex,vector< float >&vweight, float ?&WeightSum, bool ?bNormalize) ?? {????? ????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 plaincopy
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];?? ?????? }??
3、對測地距離的梯度場X，求取散度。
根據公式：
求解測地距離標量場梯度場的散度。
[cpp] ?view plaincopy
?? 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?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 plaincopy
gsi::SparseLinearSystem?*pSystemPos2?=? new ?gsi::SparseLinearSystem();?? ????gsi::Solver_TAUCS?*?pSolverPos2?=?new ?gsi::Solver_TAUCS(pSystemPos2);?? ????pSystemPos2->Resize(m_NumberV,?m_NumberV);?? ????pSystemPos2->ResizeRHS(1);?? ?????? ????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++)?? ????{?? ????????pRHSPos2[i]=float (-1.0)*m_vertices[i].VDivergence;???? ????}?? ?? ??????? ????pRHSPos2[0]=pRHSPos2[0]+10;?? ?????? ????pSystemPos2->SetRHS(0,?pRHSPos2);?? ????bool ?boksoveLs2=pSolverPos2->Solve();?? ????if ?(!boksoveLs2)assert(0);?? ????m_GeodesicsDistance.resize(m_NumberV);?? ????for ?(? int ?i=0;i<m_NumberV;++i)??? ????{????? ???????m_GeodesicsDistance[i]=(float )pSystemPos2->GetSolution(i,0);?? ????}??
可以說到了這里算法已經結束了，當然paper后面還有后續的處理，比如距離光順，還有邊界條件的等問題。但都屬于后處理，你如果要實現跟paper一模一樣的效果，還是得好好看一看后面的邊界條件問題。
邊界條件對于結果的影響還是蠻大的。
OK，做一下簡單的總結：算法整個過程就說白了就是求解一個泊松方程，說的更白一點，就是求解一個大型的方程組AX=b，因此我們的目標就是先算出A、b，其中對于一個給定的網格曲面模型來說，A就是拉普拉斯矩陣，是固定的。b的求解就是整個算法成功的關鍵。本文地址：http://blog.csdn.net/hjimce/article/details/46415499? 更多資源請關注我的博客：http://blog.csdn.net/hjimce? ? ? ? ? ? ? ? 原創文章，轉載請保留本行信息。
參考文獻：
1、《Geodesics in Heat:A New Approach to Computing Distance Based on Heat Flow》
2、《Mesh Editing with Poisson-Based Gradient Field Manipulation》
3、《Poisson Image Editing》
總結
以上是生活随笔 為你收集整理的图形处理（七）基于热传播的测地距离计算-Siggraph 2013 的全部內容，希望文章能夠幫你解決所遇到的問題。
如果覺得生活随笔 網站內容還不錯，歡迎將生活随笔 推薦給好友。
