馬健超1,2，陳國棟1,2，周霆1,2
?。?.福州大學 物理與信息工程學院，福建 福州 350116；2.福州大學 計算機圖像圖形研究所，福建 福州 350116）
摘要：虛擬手術系統研究是計算機圖形學的一個重要領域，為了提高肌肉纖維紋理合成的真實感與實時性，文章采用了基于場向參數化的方法，并結合最光滑方向場進行了三維模型上的肌肉纖維紋理映射與合成，通過布置奇異點，模擬肌肉纖維損傷的狀態。該系統實現了場向參數化的可視化及肌肉纖維紋理合成的簡單交互，實驗結果表明，系統的實時性與真實感均得到提高。
關鍵詞：虛擬手術；場向參數化；最光滑方向場；紋理合成
中圖分類號：TP391.41文獻標識碼：ADOI： 10.19358/j.issn.1674-7720.2017.09.006
引用格式：馬健超，陳國棟，周霆.基于場向參數化的肌肉纖維紋理合成研究［J］.微型機與應用，2017,36（9）：18-21.
0引言
隨著現代醫學技術和計算機技術的發展，借助計算機圖形技術分析處理醫學問題成為了醫學研究和計算機圖形技術的熱點，虛擬手術就是其中一個重要的學科發展方向，其具有沉浸性、交互性等特點［1］。
在計算機圖形學中，紋理的應用是增強真實感的有效方法，對于三維表面紋理合成，實現的目標是將二維紋理映射到三維表面，并且不產生走樣和扭曲?，F有方法通常是把表面領域通過參數化方法展平，然后沿著表面局部坐標系進行采樣，將采樣得到的平面領域與樣本紋理的領域進行比較［2］。Praun［3］等人提出用重疊紋理技術來合成三維表面紋理，該算法需要較多的人工干預；Wei［4］等人把表面紋理的顏色信息存儲在模型的頂點上，根據松弛算法生成方向場來定義三維表面局部坐標系，但合成速度不理想；Turk［5］等人使用密集模型保存豐富紋理細節，對每個目標領域進行局部參數化，合成速度較慢。這些方法都從圖像特征的統計相似性來解決問題，但不能很好地繪制出肌肉纖維紋理的特點。骨骼肌由肌腹和肌腱構成，肌腹由大量的橫紋纖維構成，不同部位間纖維方向差異較大，而肌腱由腱腱纖維構成［6］，肌肉損傷的表現為肌肉纖維沿筋膜面羽毛狀擴展，或向鄰近肌肉擴展［7］。
為了提高虛擬手術中肌肉纖維仿真的真實感和實時性，本文從幾何學的角度來看待這個問題，采用最光滑方向場與場向參數化相結合，在Ray［8］等人的全局場面參數化方法的基礎上，允許加入新的奇異點，以此來仿真肌肉損傷，并省略了旋度校正和單位化約束，使用一個簡單的凸二次能量來規劃問題，算法速度得到提升，通過簡單的交互，可以實時改變最光滑方向場的方向，從而得到三維模型上不同方向的肌肉纖維紋理，可應用到虛擬手術中。
本文實現的整體算法流程如圖1所示，通過讀取三維網格模型，如OBJ格式，使用重新定義的狄利克雷能量來獲取最光滑方向場，使用該方向場進行參數化，得到相應的紋理坐標，最后使用模型加載庫Assimp及OpenGL等開源庫進行紋理的合成。
1最光滑方向場
1.12-方向場
本文中的方向場指的是單位長度的矢量場，方向場的度n指的是曲面上每一點有n等間隔分布的方向場。如圖2所示，左圖為n為1時，即傳統意義上的場;右圖為n為2時，稱為2方向場［9］，每一點由兩個方向相反的單位矢量組成。這些方向場通常存在奇異點，即方向場在該點無法光滑地變化，注意，本文產生的方向場將自動地布置奇異點的個數與位置。
本文將曲面M上的一點的切空間視為復數域C的副本，其上的切矢量可用復數乘以一個基矢量來表示，若以實軸為基矢量，再通過平方，得到一個不需要考慮單位約束的方向場的表示方法，為：
u=z2=ei2kπ,k=0,1(1)
1.2光滑能量
通過使用狄利克雷能量來測量2方向場ψ的光滑度［10］，表示協變導數，也就是曲面上的列維奇維塔聯絡，表示為：
但方向場中，在奇異點處該能量變得無限大，導致該公式失去意義，而每一個2方向場可寫成一個比例因子與單位方向場的乘積，即：
ψ=aφ(3)
將所有的ψ中最小的狄利克雷能量定義為2方向場的能量，即式(4)，此時奇異點處的能量將變為零，從而解決了上述的問題。
遍歷所有方向場，從而得到全局最小化的方向場，也就是最光滑方向場，由于經過了縮放的2方向場的集合與2矢量場的集合是沒有區別的，所以可以替而求解。
![TRD(N0DELW}J8QOZ{6KE3]E.png](http://files.chinaaet.com/images/2017/06/06/6363238612712658213226982.png "TRD(N0DELW}J8QOZ{6KE3]E.png")
這意味著一個最小特征值的問題，三維網格上所有方向場ψ都是分段線性的，將其表示為基底的復線性組合，V代表點集。
式(5)的問題轉化為求解特征向量的問題。等式(7)中，A是一個相對于分段線性基底Ψi的埃爾米特矩陣，M是埃爾米特質量矩陣，使用逆冪迭代法來求解，得出最光滑方向場。
Au=λMu(7)
2場向參數化
2.1坐標函數
首先需要選擇參數化后的坐標函數的表達形式，Ray［8］等人和Zhang［11］等人提出使用復數域的矢量值函數，其角度分量提供最終的坐標，但每一點有一個四階的懲罰項，導致了非凸規劃問題，所以本算法忽略了這一項，同樣使用一個復數域的矢量值函數ψ，用其幅角α作為參數化后的坐標，為：
α=argψ(8)
設定最光滑方向場為X，每點的單位化函數為：
φ=eiα(9)
角α只沿著方向X以速率ν來變化，通過微分，得到
dφ=iωφ(10)
這意味著通過求解該等式可以得到角速度ω，但需要方向場總是可積的，因此考慮它的L2殘差：
![3F~RM[1]RJ((LLIW[I}M}~T.png](http://files.chinaaet.com/images/2017/06/06/6363238613734639995106661.png "3F~RM[1]RJ((LLIW[I}M}~T.png")
其中算子c=d－iω恰好是曲面M上的一個聯絡，ε可以被認為是在一個復變函數上的狄利克雷能量，所以定義單位化函數φ能量為比例函數a下，所有可能的最小的狄利克雷能量：
在單位化函數φ上的能量極小化是等價于在所有復變函數ψ上的狄利克雷能量極小化：
![]QHQ98(AFL509K}5_C%@GWE.png](http://files.chinaaet.com/images/2017/06/06/6363238614647235424604889.png "]QHQ98(AFL509K}5_C%@GWE.png")
因此對應為一個標準的求解特征值的問題，ψ為特征函數，Δs為薛定諤算子。
Δsψ=λψ(14)
2.2不可定向特征
方向場中存在不可定向的點，即在奇異點有不可定向的特征，引進Kalberer［12］等人提出的二次覆蓋的概念。除了奇異點之外，二次覆蓋看起來就像原曲面的副本，但算法中并沒有實際地去構造二次覆蓋，只是用這一思想來解決奇異點的問題。圖3為二次覆蓋于原曲面的示意圖，MD為二次覆蓋，M為原曲面，τ稱為片交換函數，ρ將二次覆蓋映射回原曲面，曲面M上的線場，在MD上任意選擇其中一個方向的場。
2.3離散化
三角網格上，每一個頂點i∈V，用單位切矢量Xi和目標線頻率νi∈R+共同定義了矢量Zi=νiXi，離散化算法的具體步驟如下：
輸入：二維單純復形K（三角網格）。
過程：
（1）將Z投影到每一條邊上，得到角位移ω；
?。?）構建類拉普拉斯矩陣A，元素由ω決定；
（3）求矩陣A最小特征值對應特征向量ψ 。
輸出：每一個面上，通過ψ和ω得到計算紋理坐標α。
首先求得角位移ωij，即邊矢量與輸入矢量場的內積
的平均值。
ωij=∫ijω=12(〈eij,Zi〉+〈eij,Zj〉)(15)
然后每個三角形構建類余切拉普拉斯矩陣，即薛定諤算子的離散化。該矩陣A是正定對稱，與余切拉普拉斯矩陣有相同結構，對于每個正則邊ij∈E，非對角塊為：
Aij=－wij［eiωij］,sij=－1
－wij［eiωij］,sij=+1
Aji=ATij(16)
矩陣A的對角塊為：
Aii=∑ij∈E［ωij］(17)
塊對角集總質量矩陣為B，用Aijk表示三角形ijk的面積，對角元素為關聯三角形總面積的三分之一。
Bii=13∑ijk∈FAijk(18)
存在ψi=ai+bii的某個值來使ε最小，用ai,bi∈R值交替組成的一個矢量x，通過求解式(19)的最小特征值對應的特征向量x，使用喬里斯基分解矩陣A，然后應用逆冪法來求解。
Ax=λBx(19)
求得ψ后，如果僅通過式(8)來求最終的坐標函數α，將產生畸變與扭曲，因此采用旋轉形式進行調整，最終坐標為：
αjki:=arg(ψi)
αkij:=arg(ψi)+σij
αijk:=arg(ψi)+σik(20)
3結果與分析
3.1編程環境與開源庫
本課題使用系統為Window 7，編程環境為VS2012，處理器為2.40 GHz，內存8 GB，顯卡為AMD Radeon 6550M。采用了開源庫SuiteSparse來實現對稀疏矩陣的運算，選用了模型加載庫Assimp，選用的三維模型格式是OBJ格式，還有OpenGL及其擴展庫GLUT和GLEW處理窗口的創建、基本的交互，使用著色語言GLSL編寫了頂點著色器和片段著色器。
3.2紋理合成結果與分析
圖4展示了整個紋理合成的過程，使用模型的面片數為6 142片，（a）為導入模型后顯示的三角網格模型；（b）為場向參數化的可視化，可看到奇異點（分叉點）；（c）為模型表面紋理合成后的效果；（d）為紋理合成的細節圖，可以看到紋理的方向與方向場一致；（e） 為所使用的紋理圖片。圖5為圖4的2方向場經過旋轉90°后的可視化結果，從右圖可以看到肌肉纖維紋理的方向已經發生了改變，但是同樣與方向場方向保持一致。
圖6展示的是參數化中奇異點與紋理合成后細節的對比，使用的模型的面片數為15 418片，可以看到紋理在奇異點處發生了分叉，但總體還是保持連續的，以此來模擬肌肉的輕微損傷，左圖為場向參數化的可視化，右圖為奇異點處的紋理細節。
通過對不同面片數的模型進行實驗，得到了以下的對比數據，表1表明了面片數從1 K~16 K的幾種情況下，紋理合成所需要的時間基本控制在0.5 s以內，這樣的性能是符合進行實時操作和交互的。表2通過與其他算法的對比，列出了部分數據，表明了本文所使用的算法性能得到提升，計算所需時間平均縮短為算法1的35.2%。
4結論
本文基于場向參數化的方法，結合了最光滑方向場，并使用了一系列開源庫，實現了對肌肉纖維紋理的合成，最終的系統能夠實現實時操作與交互，在保證真實感的同時，也證明了所使用算法的性能有較大的提高，為虛擬手術中肌肉手術的場景的進一步研究提供了條件。
參考文獻
?。?］ 邢英杰, 張少華, 劉曉冰. 虛擬手術系統技術現狀［J］. 計算機工程與應用, 2004, 40(7):88-90.
?。?］ 韓建偉. 基于樣本的三維表面紋理快速合成技術［D］. 杭州：浙江大學, 2009.
［3］ PRAUN E,FINKELSTEIN A,HOPPE H. Lapped textures［C］. Proceeding of ACM SIGGRAPH 2000.New York,USA:ACM Press/AddisonWelsey Publishing Co.2000:465-470.
?。?］ WEI L,LEVOY M. Texture synthesis over arbitrary manifold surfaces［C］.Proceedings of ACM SIGGRAPH 2001.New York,USA:ACM,2001:355-360.
［5］ TURK G. Texture synthesis on surfaces［C］. Proceedings of ACM SIGGRAPH 2001.NY, ACM,2001:347-354.
?。?］ 席占國, 喬亞亞, 沈素紅. 超聲診斷肌肉損傷的臨床價值［J］. 醫藥前沿, 2014(18):206207.
?。?］ 劉亞娟, 冉艮龍, 葉倫,等. 大腿肌肉損傷的MRI診斷［J］. 西南國防醫藥, 2015, 25(11):1222-1224.
?。?］ RAY N. Periodic global parameterization［J］. Acm Transactions on Graphics, 2006, 25(4):1460-1485.
［9］ VAXMAN A, CAMPEN M, DIAMANTI O, et al. Directional field synthesis, design, and processing［C］． SIGGRAPH Asia， 2016:1-30.
［10］ LIU B B, WENG Y L, WANG J N, et al. Orientation field guided texture synthesis［J］. Journal of Computer Science and Technology, 2013, 28(5):827-835.
?。?1］ ZHANG M, HUANG J, LIU X, et al. A wavebased anisotropic quadrangulation method［J］. Acm Transactions on Graphics, 2010, 29(4):157-166.
?。?2］ KLBERER F, NIESER M, POLTHIER K. Stripe Parameterization of Tubular Surfaces［M］. Topological Methods in Data Analysis and Visualization, 2010.