這里要跟大家分享的是2013年Siggraph上面的一篇paper，名為《Geodesics in Heat:A New Approach to Computing Distance Based on Heat Flow》，這篇paper沒有提供源代碼，但是因?yàn)樗惴ǖ乃枷胂喈?dāng)新穎，如果你之前有研究過其它的測(cè)地三角網(wǎng)格曲面上的測(cè)地距離算法，那么看到這篇paper后，你會(huì)非常的激動(dòng)，覺得這個(gè)算法相當(dāng)神奇，網(wǎng)格曲面上測(cè)地距離的計(jì)算方法又有了新的突破。因?yàn)榭吹竭@篇paper非常激動(dòng)，以至于我興奮得馬上去把代碼寫了一遍，雖然測(cè)地距離的計(jì)算方法，對(duì)于我來說沒什么用，但它的想法，它的思想值得我們學(xué)習(xí)。
《Geodesics in Heat:A New Approach to Computing Distance Based on Heat Flow》算法主要是通過向量場(chǎng)的方法，通過求解熱流，求解泊松方程。這篇paper我覺得應(yīng)該把它歸類為向量場(chǎng)類型的文獻(xiàn)，曲面上向量場(chǎng)的應(yīng)用非常廣泛，可以用于網(wǎng)格優(yōu)化重建、參數(shù)化紋理映射、網(wǎng)格變形……等，如今這篇paper又讓我明白向量場(chǎng)也可以求解測(cè)地距離。通過熱量傳播的方法，去求解測(cè)地距離，真是大牛啊
一、理論知識(shí)
在很早之前對(duì)于測(cè)地距離，大牛們就推導(dǎo)出了測(cè)地距離的求解歸結(jié)為求解eikonal equation( 程函方程 )：
且滿足邊界約束條件為：
上面符號(hào)Φ就是測(cè)地距離。
也就是不管是歐式空間還是曲面空間，其測(cè)地距離的梯度模長恒等于1，然后邊界條件：源點(diǎn)的測(cè)地距離值為0。
因此很多大牛，都是針對(duì)上面的程函方程的求解進(jìn)行研究，然而所提出的算法都是非線性的方法，因?yàn)樯厦娴姆匠叹褪欠蔷€性方程，在網(wǎng)格曲面上計(jì)算量非常大。而這篇paper作者的真正貢獻(xiàn)在于把非線性問題轉(zhuǎn)換為線性問題，到最后只需要求解一個(gè)泊松方程就可以了。開始這個(gè)算法之前，建議先看一下《Mesh Editing with Poisson-Based Gradient Field Manipulation》這篇利用泊松方程進(jìn)行三角網(wǎng)格模型編輯，其實(shí)這篇paper的思想應(yīng)該是借用了《Poisson Image Editing》，如果你已經(jīng)很熟悉《Mesh Editing with Poisson-Based Gradient Field Manipulation》那么看這篇paper將事半功倍。
1、時(shí)間離散化， 熱傳播方程時(shí)間離散化為：
2、空間離散化， 在網(wǎng)格曲面上，離散的拉普拉斯坐標(biāo)定義為：
Ai表示三角面片的面積的三份之一，j表示頂點(diǎn)i的鄰接頂點(diǎn)。對(duì)于有V個(gè)頂點(diǎn)的網(wǎng)格模型，我們可以列出上面的V個(gè)方程，寫出矩陣形式為：
其中，Lc為拉普拉矩陣，A為包含每個(gè)頂點(diǎn)面積的對(duì)角矩陣。
代入公式(3)可得：
在網(wǎng)格曲面上，對(duì)于一個(gè)給定的三角面片平面上，標(biāo)量場(chǎng)的梯度計(jì)算公式如下：
Af表示三角形的面積，N表示三角面片的法矢，ui就是標(biāo)量場(chǎng)，我們可以把網(wǎng)格曲面每個(gè)頂點(diǎn)的測(cè)地距離看成是一個(gè)標(biāo)量場(chǎng)。
頂點(diǎn)i的散度離散形式為：
懶得多說廢話了，總之到最后歸結(jié)為求解如下方程：
其中b我們需要先通過求解標(biāo)量場(chǎng)u，然后根據(jù)散度的計(jì)算公式，計(jì)算出頂點(diǎn)測(cè)地距離的散度。
二、算法實(shí)現(xiàn)
算法總流程：
示意圖：
其中，Φ就是我們要求的測(cè)地距離，t是一個(gè)無窮小的數(shù)，趨近于0
算法具體實(shí)現(xiàn)，先說明一下，我下面自己寫的代碼很亂，懶得整理，因?yàn)槲抑皇菫榱藢W(xué)習(xí)：
1、求解u
根據(jù)公式：
求解u，因此我們需要先算出A，t，Lc，還有δ。接著我將逐漸講解著四個(gè)參數(shù)的求解：
a、時(shí)間t計(jì)算 ：見文獻(xiàn)3.2.4，時(shí)間理論上來說是一個(gè)非常小的數(shù)，文中作者給出了其合適的值為：網(wǎng)格平均邊長的平方
代碼實(shí)現(xiàn)如下：
[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、面積對(duì)角矩陣A的計(jì)算：
[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計(jì)算：
[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()?);?? ????????}????? ????}?? ?? }??
鄰接頂點(diǎn)的余切cot權(quán)重計(jì)算：
[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、求解三角網(wǎng)格曲面的熱量場(chǎng)▽u (Heat flow)：
這一步直接根據(jù)公式：
求解就可以了。然后 對(duì) ▽u進(jìn)行歸一化，并取熱量場(chǎng)的反方向，即求測(cè)地距離的梯度場(chǎng)：
[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、對(duì)測(cè)地距離的梯度場(chǎng)X，求取散度。
根據(jù)公式：
求解測(cè)地距離標(biāo)量場(chǎng)梯度場(chǎng)的散度。
[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、最后求解通過求解泊松方程，求取網(wǎng)格曲面上每個(gè)頂點(diǎn)的測(cè)地距離。
[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);?? ????}??
可以說到了這里算法已經(jīng)結(jié)束了，當(dāng)然paper后面還有后續(xù)的處理，比如距離光順，還有邊界條件的等問題。但都屬于后處理，你如果要實(shí)現(xiàn)跟paper一模一樣的效果，還是得好好看一看后面的邊界條件問題。
邊界條件對(duì)于結(jié)果的影響還是蠻大的。
OK，做一下簡單的總結(jié)：算法整個(gè)過程就說白了就是求解一個(gè)泊松方程，說的更白一點(diǎn)，就是求解一個(gè)大型的方程組AX=b，因此我們的目標(biāo)就是先算出A、b，其中對(duì)于一個(gè)給定的網(wǎng)格曲面模型來說，A就是拉普拉斯矩陣，是固定的。b的求解就是整個(gè)算法成功的關(guān)鍵。本文地址：http://blog.csdn.net/hjimce/article/details/46415499? 更多資源請(qǐng)關(guān)注我的博客：http://blog.csdn.net/hjimce? ? ? ? ? ? ? ? 原創(chuàng)文章，轉(zhuǎn)載請(qǐng)保留本行信息。
參考文獻(xiàn)：
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》
總結(jié)
以上是生活随笔 為你收集整理的图形处理（七）基于热传播的测地距离计算-Siggraph 2013 的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
如果覺得生活随笔 網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔 推薦給好友。
