圖形處理（十）測地極座標參數化_圖形處理

一、相關理論

測地極座標參數化的paper是我研究生學習階段，寫的第二篇沒有源碼的外文文獻，年少無知，看文獻的也是沒有完全看懂，就開始寫代碼，然後通過代碼不斷的根據文獻演算法，進行調整，終於把paper的代碼實現了。在此把paper的演算法詳細解讀一下，以供正在苦逼的搞三維圖形演算法的同道中人學習。

三維的演算法比較難，網上資料也比較少，國內除了到知網、萬方等資料庫可以有文獻看，代碼方面可以說沒有，國外要找到提高源碼的paper也很不容易，所以要把三維CAGD的一些演算法學懂實屬不易，在此祭奠曾經走過的三維苦逼演算法學習之路。本篇博文主要講解文獻：《Geodesic polar coordinates on polygonal》，這裡中文我又把它翻譯為：測地極座標參數化，離散冪映射演算法中每個頂點參數值取決於該頂點的法矢，因此網格曲率變化較大時經常使得參數化結果發生較大扭曲。相比於離散指數映射演算法，Eivind等提出的測地極座標參數化方法具有更好的穩定性和魯棒性。

給定具有n個頂點的三角網格模型M=(V,E,T),V為頂點集，E為邊集，T為三角面片集合。設Tijk=[vi,vj,vk],設s為網格曲面上的任意點，切向量X為基方向，如圖。

二、演算法講解

在網格曲面上，對於給定的參數化種子點S，我們把網格曲面上的頂點分為兩類：1-ring、k-ring(k>=2)，其中1-ring表示與s直接相連的頂點，k-ring表示不與s相鄰的頂點。對於測地極座標參數化方法來說，我們需要根據這兩種頂點，使用不同的計算公式。

1、1-ring參數化方法

對於網格曲面上的任意點s，根據s點位於網格三角面片的位置可分為三角面片內點、邊點、頂點三種情況：

①s為三角面片T0內的點，如圖 s1點，此時只需以T0作為參數化平面，以s點為原點，即可求取T0三個頂點的極座標

②s位於網格邊上E0，如圖s2點，此時只需以E0作為鉸鏈，固定其中的一個平面，以E0作為鉸鏈，把另一個鄰接三角面片旋轉到同一個平面上，

③s為網格曲面的頂點，使用:《Free-form shape design using triangulated surfaces》提出的極座標映射：即保證s點鄰接邊的長度不變，鄰接角度為的縮放比例為，其中θi為頂點的鄰接角，即對於s點鄰接頂點的局部參數化滿足：

其中v和w為s點的鄰接頂點，根據公式可對s點的鄰接三角面片作局部展平參數化。

 

 2、k-ring參數化

如圖所示，設網格曲面上有4點，分別為s、vi、vj、vk，s點為源點，vi、vj、vk為同一個三角面片的三個頂點，且vi測地極座標未知，vj、vk為極座標已知的點，分別為(Rj,θj)、(Rk,θk)，目標是求取vi的極座標，即vi在參數空間Tp上的映射點。

測地極座標參數化方法的基本思想是：以vi'(0,0)為原點，對vi的一環鄰域點作局部參數化，接著分別以vk'、vj'為圓點，以Rk、Rj為半徑作圓，則可以獲得兩圓的交點s'。以s'作為偽源點，通過該偽源點計算vi的極座標。

設ek=vk'-vi'，ej=vj'-vi'，ekj=vk'-vj'，則偽源點s'可以表示為：

s'=xjej+xkek

式中，

，

根據(8)式可計算s'在vi局部參數空間的座標。由此可得vi映射至Tp的極徑：

極角計算公式：

θi=(1-α)θj+αθk

式中，α=ϕij/ϕkj，ϕkj是向量vk's'與向量vj's'之間的夾角，ϕij是向量vi's'與向量vj's'之間的夾角。

測地極座標方法：首先對源點的一環鄰域點作局部參數化，接著用極徑作為Dijkstra演算法往外擴散的依據。Dijkstra演算法擴散過程中，用公式(9)(10)不斷更新活動頂點的極座標。

三、演算法實現

下面是paper中的虛擬碼： [cpp]  view plain  copy <span style="font-size:18px;">Algorithm 1: Pseudo code for computing DGPC on a mesh   1: for i = 1, · · · ,n do   2: U[i] =∞   3: end for   4: initializeNeighbourhood(s)   5: candidates.push( neighbourhood( s ) )   6: while candidates.notEmpty() do   7: j = candidates.getSmallestNode()   8: for i ∈ neighbours(j) do   9: newUi = computeDistance(i)   10: if U[i] / newUi > 1 + eps then   11: U[i] = newUi   12: theta[i] = computeAngle(i)   13: if newUi < Umax then   14: candidates.push(i)   15: end if   16: end if   17: end for   18: end while</span>  

參考文獻：

1、《Differential coordinates for interactive mesh editing》

2、《Interactive decal compositing with discrete exponential maps》

3、《Part-based representation and editing of 3D surface models》

4、《Geodesic polar coordinates on polygonal》 **********************作者：hjimce   時間：2015.5.10  聯絡QQ：1393852684   地址：http://blog.csdn.net/hjimce 原創文章，著作權，轉載請保留本行資訊********************