基于場(chǎng)向參數(shù)化的肌肉纖維紋理合成研究

馬健超，陳國(guó)棟，周 霆

(1.福州大學(xué) 物理與信息工程學(xué)院，福建 福州 350116；2.福州大學(xué) 計(jì)算機(jī)圖像圖形研究所，福建 福州 350116)

基于場(chǎng)向參數(shù)化的肌肉纖維紋理合成研究

馬健超1,2，陳國(guó)棟1,2，周 霆1,2

(1.福州大學(xué) 物理與信息工程學(xué)院，福建 福州 350116；2.福州大學(xué) 計(jì)算機(jī)圖像圖形研究所，福建 福州 350116)

虛擬手術(shù)系統(tǒng)研究是計(jì)算機(jī)圖形學(xué)的一個(gè)重要領(lǐng)域，為了提高肌肉纖維紋理合成的真實(shí)感與實(shí)時(shí)性，文章采用了基于場(chǎng)向參數(shù)化的方法，并結(jié)合最光滑方向場(chǎng)進(jìn)行了三維模型上的肌肉纖維紋理映射與合成，通過布置奇異點(diǎn)，模擬肌肉纖維損傷的狀態(tài)。該系統(tǒng)實(shí)現(xiàn)了場(chǎng)向參數(shù)化的可視化及肌肉纖維紋理合成的簡(jiǎn)單交互，實(shí)驗(yàn)結(jié)果表明，系統(tǒng)的實(shí)時(shí)性與真實(shí)感均得到提高。

虛擬手術(shù)；場(chǎng)向參數(shù)化；最光滑方向場(chǎng)；紋理合成

0 引言

隨著現(xiàn)代醫(yī)學(xué)技術(shù)和計(jì)算機(jī)技術(shù)的發(fā)展，借助計(jì)算機(jī)圖形技術(shù)分析處理醫(yī)學(xué)問題成為了醫(yī)學(xué)研究和計(jì)算機(jī)圖形技術(shù)的熱點(diǎn)，虛擬手術(shù)就是其中一個(gè)重要的學(xué)科發(fā)展方向，其具有沉浸性、交互性等特點(diǎn)[1]。

在計(jì)算機(jī)圖形學(xué)中，紋理的應(yīng)用是增強(qiáng)真實(shí)感的有效方法，對(duì)于三維表面紋理合成，實(shí)現(xiàn)的目標(biāo)是將二維紋理映射到三維表面，并且不產(chǎn)生走樣和扭曲。現(xiàn)有方法通常是把表面領(lǐng)域通過參數(shù)化方法展平，然后沿著表面局部坐標(biāo)系進(jìn)行采樣，將采樣得到的平面領(lǐng)域與樣本紋理的領(lǐng)域進(jìn)行比較[2]。Praun[3]等人提出用重疊紋理技術(shù)來合成三維表面紋理，該算法需要較多的人工干預(yù)；Wei[4]等人把表面紋理的顏色信息存儲(chǔ)在模型的頂點(diǎn)上，根據(jù)松弛算法生成方向場(chǎng)來定義三維表面局部坐標(biāo)系，但合成速度不理想；Turk[5]等人使用密集模型保存豐富紋理細(xì)節(jié)，對(duì)每個(gè)目標(biāo)領(lǐng)域進(jìn)行局部參數(shù)化，合成速度較慢。這些方法都從圖像特征的統(tǒng)計(jì)相似性來解決問題，但不能很好地繪制出肌肉纖維紋理的特點(diǎn)。骨骼肌由肌腹和肌腱構(gòu)成，肌腹由大量的橫紋纖維構(gòu)成，不同部位間纖維方向差異較大，而肌腱由腱腱纖維構(gòu)成[6]，肌肉損傷的表現(xiàn)為肌肉纖維沿筋膜面羽毛狀擴(kuò)展，或向鄰近肌肉擴(kuò)展[7]。

為了提高虛擬手術(shù)中肌肉纖維仿真的真實(shí)感和實(shí)時(shí)性，本文從幾何學(xué)的角度來看待這個(gè)問題，采用最光滑方向場(chǎng)與場(chǎng)向參數(shù)化相結(jié)合，在Ray[8]等人的全局場(chǎng)面參數(shù)化方法的基礎(chǔ)上，允許加入新的奇異點(diǎn)，以此來仿真肌肉損傷，并省略了旋度校正和單位化約束，使用一個(gè)簡(jiǎn)單的凸二次能量來規(guī)劃問題，算法速度得到提升，通過簡(jiǎn)單的交互，可以實(shí)時(shí)改變最光滑方向場(chǎng)的方向，從而得到三維模型上不同方向的肌肉纖維紋理，可應(yīng)用到虛擬手術(shù)中。

本文實(shí)現(xiàn)的整體算法流程如圖1所示，通過讀取三維網(wǎng)格模型，如OBJ格式，使用重新定義的狄利克雷能量來獲取最光滑方向場(chǎng)，使用該方向場(chǎng)進(jìn)行參數(shù)化，得到相應(yīng)的紋理坐標(biāo)，最后使用模型加載庫Assimp及OpenGL等開源庫進(jìn)行紋理的合成。

圖1 整體算法流程圖

1 最光滑方向場(chǎng)

1.1 2-方向場(chǎng)

本文中的方向場(chǎng)指的是單位長(zhǎng)度的矢量場(chǎng)，方向場(chǎng)的度n指的是曲面上每一點(diǎn)有n等間隔分布的方向場(chǎng)。如圖2所示，左圖為n為1時(shí)，即傳統(tǒng)意義上的場(chǎng);右圖為n為2時(shí)，稱為2-方向場(chǎng)[9]，每一點(diǎn)由兩個(gè)方向相反的單位矢量組成。這些方向場(chǎng)通常存在奇異點(diǎn)，即方向場(chǎng)在該點(diǎn)無法光滑地變化，注意，本文產(chǎn)生的方向場(chǎng)將自動(dòng)地布置奇異點(diǎn)的個(gè)數(shù)與位置。

圖2 方向場(chǎng)的圖示

本文將曲面M上的一點(diǎn)的切空間視為復(fù)數(shù)域C的副本，其上的切矢量可用復(fù)數(shù)乘以一個(gè)基矢量來表示，若以實(shí)軸為基矢量，再通過平方，得到一個(gè)不需要考慮單位約束的方向場(chǎng)的表示方法，為：

u=z2=ei2kπ,k=0,1

(1)

1.2 光滑能量

通過使用狄利克雷能量來測(cè)量2-方向場(chǎng)ψ的光滑度[10]，▽表示協(xié)變導(dǎo)數(shù)，也就是曲面上的列維-奇維塔聯(lián)絡(luò)，表示為：

(2)

但方向場(chǎng)中，在奇異點(diǎn)處該能量變得無限大，導(dǎo)致該公式失去意義，而每一個(gè)2-方向場(chǎng)可寫成一個(gè)比例因子與單位方向場(chǎng)的乘積，即：

ψ=aφ

(3)

將所有的ψ中最小的狄利克雷能量定義為2-方向場(chǎng)的能量，即式(4)，此時(shí)奇異點(diǎn)處的能量將變?yōu)榱悖瑥亩鉀Q了上述的問題。

(4)

遍歷所有方向場(chǎng)，從而得到全局最小化的方向場(chǎng)，也就是最光滑方向場(chǎng)，由于經(jīng)過了縮放的2-方向場(chǎng)的集合與2-矢量場(chǎng)的集合是沒有區(qū)別的，所以可以替而求解。

(5)

這意味著一個(gè)最小特征值的問題，三維網(wǎng)格上所有方向場(chǎng)ψ都是分段線性的，將其表示為基底的復(fù)線性組合，V代表點(diǎn)集。

(6)

式(5)的問題轉(zhuǎn)化為求解特征向量的問題。等式(7)中，A是一個(gè)相對(duì)于分段線性基底Ψi的埃爾米特矩陣，M是埃爾米特質(zhì)量矩陣，使用逆冪迭代法來求解，得出最光滑方向場(chǎng)。

Au=λMu

(7)

2 場(chǎng)向參數(shù)化

2.1 坐標(biāo)函數(shù)

首先需要選擇參數(shù)化后的坐標(biāo)函數(shù)的表達(dá)形式，Ray[8]等人和Zhang[11]等人提出使用復(fù)數(shù)域的矢量值函數(shù)，其角度分量提供最終的坐標(biāo)，但每一點(diǎn)有一個(gè)四階的懲罰項(xiàng)，導(dǎo)致了非凸規(guī)劃問題，所以本算法忽略了這一項(xiàng)，同樣使用一個(gè)復(fù)數(shù)域的矢量值函數(shù)ψ，用其幅角α作為參數(shù)化后的坐標(biāo)，為：

α=argψ

(8)

設(shè)定最光滑方向場(chǎng)為X，每點(diǎn)的單位化函數(shù)為：

φ=eiα

(9)

角α只沿著方向X以速率ν來變化，通過微分，得到

dφ=iωφ

(10)

這意味著通過求解該等式可以得到角速度ω，但需要方向場(chǎng)總是可積的，因此考慮它的L2-殘差：

ε(φ)=∫M▽c=∫M|(d-iω)φ|2

(11)

其中算子▽c=d-iω恰好是曲面M上的一個(gè)聯(lián)絡(luò)，ε可以被認(rèn)為是在一個(gè)復(fù)變函數(shù)上的狄利克雷能量，所以定義單位化函數(shù)φ能量為比例函數(shù)a下，所有可能的最小的狄利克雷能量：

(12)

在單位化函數(shù)φ上的能量極小化是等價(jià)于在所有復(fù)變函數(shù)ψ上的狄利克雷能量極小化：

(13)

因此對(duì)應(yīng)為一個(gè)標(biāo)準(zhǔn)的求解特征值的問題，ψ為特征函數(shù)，Δs為薛定諤算子。

Δsψ=λψ

(14)

2.2 不可定向特征

方向場(chǎng)中存在不可定向的點(diǎn)，即在奇異點(diǎn)有不可定向的特征，引進(jìn)Kalberer[12]等人提出的二次覆蓋的概念。除了奇異點(diǎn)之外，二次覆蓋看起來就像原曲面的副本，但算法中并沒有實(shí)際地去構(gòu)造二次覆蓋，只是用這一思想來解決奇異點(diǎn)的問題。圖3為二次覆蓋于原曲面的示意圖，MD為二次覆蓋，M為原曲面，τ稱為片交換函數(shù)，ρ將二次覆蓋映射回原曲面，曲面M上的線場(chǎng)，在MD上任意選擇其中一個(gè)方向的場(chǎng)。

圖3 二次覆蓋與原曲面

2.3 離散化

三角網(wǎng)格上，每一個(gè)頂點(diǎn)i∈V，用單位切矢量Xi和目標(biāo)線頻率νi∈R+共同定義了矢量Zi=νiXi，離散化算法的具體步驟如下：

輸入：二維單純復(fù)形K(三角網(wǎng)格)。

過程：

(1)將Z投影到每一條邊上，得到角位移ω；

(2)構(gòu)建類拉普拉斯矩陣A，元素由ω決定；

(3)求矩陣A最小特征值對(duì)應(yīng)特征向量ψ。

輸出：每一個(gè)面上，通過ψ和ω得到計(jì)算紋理坐標(biāo)α。

首先求得角位移ωij，即邊矢量與輸入矢量場(chǎng)的內(nèi)積的平均值。

圖4 參數(shù)化可視化與紋理合成結(jié)果

(15)

然后每個(gè)三角形構(gòu)建類余切-拉普拉斯矩陣，即薛定諤算子的離散化。該矩陣A是正定對(duì)稱，與余切-拉普拉斯矩陣有相同結(jié)構(gòu)，對(duì)于每個(gè)正則邊ij∈E，非對(duì)角塊為：

(16)

矩陣A的對(duì)角塊為：

(17)

塊對(duì)角集總質(zhì)量矩陣為B，用Aijk表示三角形ijk的面積，對(duì)角元素為關(guān)聯(lián)三角形總面積的三分之一。

(18)

存在ψi=ai+bii的某個(gè)值來使ε最小，用ai,bi∈R值交替組成的一個(gè)矢量x，通過求解式(19)的最小特征值對(duì)應(yīng)的特征向量x，使用喬里斯基分解矩陣A，然后應(yīng)用逆冪法來求解。

Ax=λBx

(19)

求得ψ后，如果僅通過式(8)來求最終的坐標(biāo)函數(shù)α，將產(chǎn)生畸變與扭曲，因此采用旋轉(zhuǎn)形式進(jìn)行調(diào)整，最終坐標(biāo)為：

(20)

3 結(jié)果與分析

3.1 編程環(huán)境與開源庫

本課題使用系統(tǒng)為Window 7，編程環(huán)境為VS2012，處理器為2.40 GHz，內(nèi)存8 GB，顯卡為AMD Radeon 6550M。采用了開源庫SuiteSparse來實(shí)現(xiàn)對(duì)稀疏矩陣的運(yùn)算，選用了模型加載庫Assimp，選用的三維模型格式是OBJ格式，還有OpenGL及其擴(kuò)展庫GLUT和GLEW處理窗口的創(chuàng)建、基本的交互，使用著色語言GLSL編寫了頂點(diǎn)著色器和片段著色器。

3.2 紋理合成結(jié)果與分析

圖4展示了整個(gè)紋理合成的過程，使用模型的面片數(shù)為6 142片，(a)為導(dǎo)入模型后顯示的三角網(wǎng)格模型；(b)為場(chǎng)向參數(shù)化的可視化，可看到奇異點(diǎn)(分叉點(diǎn))；(c)為模型表面紋理合成后的效果；(d)為紋理合成的細(xì)節(jié)圖，可以看到紋理的方向與方向場(chǎng)一致；(e) 為所使用的紋理圖片。圖5為圖4的2-方向場(chǎng)經(jīng)過旋轉(zhuǎn)90°后的可視化結(jié)果，從右圖可以看到肌肉纖維紋理的方向已經(jīng)發(fā)生了改變，但是同樣與方向場(chǎng)方向保持一致。

圖5 旋轉(zhuǎn)方向場(chǎng)后的可視化與對(duì)應(yīng)的紋理合成結(jié)果

圖6展示的是參數(shù)化中奇異點(diǎn)與紋理合成后細(xì)節(jié)的對(duì)比，使用的模型的面片數(shù)為15 418片，可以看到紋理在奇異點(diǎn)處發(fā)生了分叉，但總體還是保持連續(xù)的，以此來模擬肌肉的輕微損傷，左圖為場(chǎng)向參數(shù)化的可視化，右圖為奇異點(diǎn)處的紋理細(xì)節(jié)。

表1 不同模型性能對(duì)比

圖6 奇異點(diǎn)與紋理細(xì)節(jié)

通過對(duì)不同面片數(shù)的模型進(jìn)行實(shí)驗(yàn)，得到了以下的對(duì)比數(shù)據(jù)，表1表明了面片數(shù)從1 K～16 K的幾種情況下，紋理合成所需要的時(shí)間基本控制在0.5 s以內(nèi)，這樣的性能是符合進(jìn)行實(shí)時(shí)操作和交互的。表2通過與其他算法的對(duì)比，列出了部分?jǐn)?shù)據(jù)，表明了本文所使用的算法性能得到提升，計(jì)算所需時(shí)間平均縮短為算法1的35.2%。

表2 不同算法性能對(duì)比

4 結(jié)論

本文基于場(chǎng)向參數(shù)化的方法，結(jié)合了最光滑方向場(chǎng)，并使用了一系列開源庫，實(shí)現(xiàn)了對(duì)肌肉纖維紋理的合成，最終的系統(tǒng)能夠?qū)崿F(xiàn)實(shí)時(shí)操作與交互，在保證真實(shí)感的同時(shí)，也證明了所使用算法的性能有較大的提高，為虛擬手術(shù)中肌肉手術(shù)的場(chǎng)景的進(jìn)一步研究提供了條件。

[1] 邢英杰, 張少華, 劉曉冰. 虛擬手術(shù)系統(tǒng)技術(shù)現(xiàn)狀[J]. 計(jì)算機(jī)工程與應(yīng)用, 2004, 40(7):88-90.

[2] 韓建偉. 基于樣本的三維表面紋理快速合成技術(shù)[D]. 杭州：浙江大學(xué), 2009.

[3] PRAUN E,FINKELSTEIN A,HOPPE H. Lapped textures[C]. Proceeding of ACM SIGGRAPH 2000.New York,USA:ACM Press/Addison-Welsey Publishing Co.2000:465-470.

[4] WEI L,LEVOY M. Texture synthesis over a-rbitrary 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.

[6] 席占國(guó), 喬亞亞, 沈素紅. 超聲診斷肌肉損傷的臨床價(jià)值[J]. 醫(yī)藥前沿, 2014(18):206-207.

[7] 劉亞娟, 冉艮龍, 葉倫,等. 大腿肌肉損傷的MRI診斷[J]. 西南國(guó)防醫(yī)藥, 2015, 25(11):1222-1224.

[8] 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.

[11] ZHANG M, HUANG J, LIU X, et al. A wave-based anisotropic quadrangulation method[J]. Acm Transactions on Graphics, 2010, 29(4):157-166.

The research of muscle fiber texture synthesis based on field-aligned parameterization

Ma Jianchao1,2，Chen Guodong1,2，Zhou Ting1,2

(1.College of Physics and Information Engineering,Fuzhou University,Fuzhou 350116,China;2.Institute of Computer Image and Graphics,Fuzhou University,Fuzhou 350116,China)

The research of virtual surgery system is an important research filed of computer graphics. In order to improve the realism and real-time of muscle fiber texture synthesis, this paper adopted a method which based on field-aligned parameterization, combined with the smoothest direction field, the algorithm achieved the muscle fiber texture mapping and synthesis on the 3D models. It can mimic the state of muscle fiber injury by arranging singular points. This system implements the visualization of field-aligned parameterization and simple interaction of muscle fiber texture synthesis. The results show that the performance and realism of the system has been great improved.

virtual surgery; field-aligned parameterization; smoothest direction field; texture synthesis

TP391.41

A

10.19358/j.issn.1674- 7720.2017.09.006

馬健超，陳國(guó)棟，周霆.基于場(chǎng)向參數(shù)化的肌肉纖維紋理合成研究[J].微型機(jī)與應(yīng)用，2017,36(9)：18-21.

2016-09-15)

馬健超(1991-)，男，碩士研究生，主要研究方向：計(jì)算機(jī)圖形學(xué)、數(shù)字幾何處理。

陳國(guó)棟(1979-)，男，博士，副研究員，主要研究方向：計(jì)算機(jī)圖形學(xué)、虛擬現(xiàn)實(shí)技術(shù)。

周霆(1973-)，女，博士，副教授，主要研究方向：信號(hào)處理、圖像處理和傳輸。