陳國(guó)棟, 潘冠慈, 田 影
(福州大學(xué)物理與信息工程學(xué)院,福建 福州 350116)
虛擬手術(shù)又稱(chēng)手術(shù)仿真,是虛擬現(xiàn)實(shí)技術(shù)在醫(yī)學(xué)領(lǐng)域中的應(yīng)用。在虛擬手術(shù)場(chǎng)景中,要建立真實(shí)的視覺(jué)反饋系統(tǒng),除了要正確逼真地進(jìn)行器官模型建模外,還要通過(guò)人體器官的各種紋理細(xì)節(jié)來(lái)豐富其表現(xiàn)特征[1]。
基于樣圖的紋理合成技術(shù)在過(guò)去二十年一直是研究的熱點(diǎn)。傳統(tǒng)的算法主要以馬爾科夫隨機(jī)場(chǎng)(MRF)為模型通過(guò)像素點(diǎn)或紋理快的查找拼接來(lái)完成。2012年ImageNet比賽結(jié)果讓深度學(xué)習(xí)成為學(xué)術(shù)界和工業(yè)界研究的熱點(diǎn)。2015年Gatys[2]等人提出基于卷積神經(jīng)網(wǎng)絡(luò)(CNN)的全新紋理合成技術(shù)。該方法相對(duì)于普通的參數(shù)法最大的突破之一就是降低了對(duì)復(fù)雜濾波器的設(shè)計(jì)。圖像中的“深”層特征通過(guò)多層已經(jīng)訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)隱式地學(xué)習(xí)從而解決復(fù)雜的任務(wù)。文章提出的算法基于Gatys的研究理論并結(jié)合2017年Cohen[3]等人提出的紋理深度相關(guān)性理論進(jìn)行深入研究,在神經(jīng)網(wǎng)絡(luò)的高層引入馬爾科夫隨機(jī)場(chǎng)(MRF)以合成具有全局結(jié)構(gòu)特征以及局部和非局部結(jié)構(gòu)特征相結(jié)合的肺表面紋理圖像。
該算法是基于CNN來(lái)合成肺表面紋理圖像,使用已經(jīng)訓(xùn)好的VGG-19[4]網(wǎng)絡(luò)的16個(gè)卷積層和5個(gè)池化層的特征空間,并且去掉了最后3個(gè)全連接層。輸入圖像的大小為224×224,網(wǎng)絡(luò)架構(gòu)基于如下兩類(lèi)基本運(yùn)算:
(1)卷積濾波器大小為3×3×k,其中k是輸入特征圖的數(shù)量。卷積的步長(zhǎng)為1以確保輸出特征圖和輸入特征圖的空間維度一樣。
(2)不重疊2×2區(qū)域的池化層,采用平均池化來(lái)對(duì)特征圖進(jìn)行降采樣。這兩類(lèi)運(yùn)算交替進(jìn)行。
(1)
其中q和m是對(duì)應(yīng)特征圖的位置索引。則第l層的Gram矩陣可表示為:
(2)
Gram矩陣表示該層各個(gè)維度自身的特征(對(duì)角線(xiàn))以及各個(gè)維度之間的關(guān)系,最終生成的{G1,G2,…Gl}是整個(gè)CNN網(wǎng)絡(luò)產(chǎn)生的所有層關(guān)系矩陣的集合,該集合表示紋理圖的局部結(jié)構(gòu)特征。第l層對(duì)總損失函數(shù)的貢獻(xiàn)可表示為:
(3)
(4)
(5)
圖1 卷積操作得到特征圖
根據(jù)肺表面紋理的特征,從視覺(jué)上可以直觀的得出該紋理具有一定的相似性,即在全局部結(jié)構(gòu)上保持了紋理的規(guī)律性。Gram矩陣定義了特征間的相關(guān)性,即圖像的局部結(jié)構(gòu)信息,而圖像的非局部結(jié)構(gòu)信息則用特征內(nèi)的相關(guān)性來(lái)表示,以下給出不同層特征內(nèi)的深度相關(guān)矩陣的定義,其作用是引導(dǎo)紋理合成過(guò)程的優(yōu)化以保持整體結(jié)構(gòu)的規(guī)律性。
(6)
Rl,n∈Ql×Ml,公式(6)說(shuō)明了深度相關(guān)矩陣是將垂直移動(dòng)i個(gè)像素,水平移動(dòng)j個(gè)像素,并在重疊區(qū)域上應(yīng)用逐點(diǎn)乘法,最后再由重疊區(qū)域總量的倒數(shù)加權(quán)求和。其中wi,j的定義如下:
wi,j=[(Q-|i|)(M-|j|)]-1
(7)
基于上述深度相關(guān)矩陣,可以定義第l層特征內(nèi)的相關(guān)性損失函數(shù)Dcor:
(8)
計(jì)算出該損失函數(shù)對(duì)該層特征矩陣的偏導(dǎo)數(shù):
(9)
其中*表示卷積。與基于Gram矩陣定義的損失函數(shù)類(lèi)似,通過(guò)計(jì)算各層基于DCor的損失函數(shù),最后得出基于DCor的總損失函數(shù):
(10)
圖2 左邊為輸入樣圖,右邊為輸出結(jié)果
Gram矩陣本質(zhì)上屬于像素間的一一對(duì)應(yīng),計(jì)算量相當(dāng)大。為提高算法效率,僅在pool1,pool2層計(jì)算Gram矩陣。在網(wǎng)絡(luò)的更高層(relu3_1和relu4_1)引入馬爾科夫隨機(jī)場(chǎng)(MRF),將特征圖分成若干圖像塊,尋找圖像塊之間的匹配,由于不再是像素上的一一對(duì)應(yīng),計(jì)算效率顯著提高。MRF對(duì)于圖像高層語(yǔ)義的空間布局限制起到了關(guān)鍵作用,相當(dāng)于在DCor相關(guān)項(xiàng)的約束下進(jìn)行二次約束。
2016年Li等人[5]在基于MRF的圖像風(fēng)格轉(zhuǎn)移上給出了基于MRF的損失函數(shù),定義如下:
(11)
(12)
2.1.1 合成多樣性肺表面紋理圖
長(zhǎng)期抽煙或是生活在霧霾地區(qū)肺表面紋理會(huì)發(fā)生相應(yīng)變化,具體可表現(xiàn)為肺表面斑點(diǎn)狀紋理增多且大規(guī)模出現(xiàn)黑色斑塊紋理。當(dāng)肺出現(xiàn)病變時(shí)也會(huì)出現(xiàn)相應(yīng)的斑點(diǎn)或斑塊紋理,且局部特征更為突出。在2.1的基礎(chǔ)上加入兩個(gè)新的特征損失表達(dá)式,以合成多樣性的肺表面紋理圖。為了增加紋理的局部細(xì)節(jié)即表面斑點(diǎn)狀紋理圖定義如下l層的損失函數(shù):
(13)
(14)
總損失函數(shù)為:
(15)
隨著病變程度加深肺表面紋理黑色斑塊進(jìn)一步擴(kuò)大,隨機(jī)分布在肺表面,且紋理邊緣變得更明顯,常見(jiàn)的紋理狀態(tài)可呈白色糜爛狀或邊緣黑色加粗,因此根據(jù)紋理邊緣的特點(diǎn),可定義一種邊緣保留的損失函數(shù)。2016年Johnson等人[6]在圖像超分辨率重建中提出一種平滑度損失函數(shù),并將其定義為總變差項(xiàng),該項(xiàng)會(huì)使紋理過(guò)度平滑以至于邊緣被弱化。2017年Cohen等人[3]在其啟發(fā)下提出一種保留邊緣的損失函數(shù)。
(16)
其中σ為圖像標(biāo)準(zhǔn)差,δi,δj為特征圖上相鄰的索引(4鄰域像素)。根據(jù)紋理表面病變的特點(diǎn),為了突出局部特征,在特征圖的鄰域像素中引入局部標(biāo)準(zhǔn)差(LSD),LSD濾波器的大小為3×3,與卷積濾波器的大小保持一致。σδ定義為:
(17)
(18)
2.1.1.1 肺表面紋理合成綜述
(19)
算法對(duì)規(guī)定的網(wǎng)絡(luò)層使用指定的損失函數(shù),實(shí)驗(yàn)中具體的超參數(shù)如表1所示,算法運(yùn)行時(shí)間等數(shù)據(jù)如表2所示,總損失函數(shù)E在迭代過(guò)程中的誤差變化如圖5所示。BP算法中使用隨機(jī)梯度下降進(jìn)行參數(shù)優(yōu)化,并用有界約束的BFGS[7]算法進(jìn)行測(cè)試。圖像I被初始化為標(biāo)準(zhǔn)差為σN=1×10-3的高斯白噪聲。圖4給出了基于該算法得出的目標(biāo)紋理圖像,其中4.1是樣本紋理圖像,圖像來(lái)自網(wǎng)絡(luò),其中黑色斑塊為常年抽煙所致,病變主要體現(xiàn)在黑色斑塊不規(guī)則分布且顏色深淺不一,中間框圖部分為圖像殘缺部分,但在4.2-4.5的輸出圖像中不僅殘缺部分得到修復(fù),并且圖像的整體紋理結(jié)構(gòu)并沒(méi)被改變。4.2、4.3、4.4、4.5分別是γ=-0.75×10-4,γ=-0.75×10-3,γ=-0.75×10-5,γ=-0.75×10-2時(shí)的目標(biāo)紋理圖像,從4.4與4.5的對(duì)比可以看出γ取較大值時(shí)對(duì)紋理邊緣的影響較大,在整體結(jié)構(gòu)保持大致相同的情況下(EDCor和MRF約束),白色方框處可以看出紋理的局部細(xì)節(jié)得到增強(qiáng),帶有病灶的紋理邊緣更明顯(EDiv與ESmooth作用)。
圖3 算法框圖
損失函數(shù)對(duì)應(yīng)層層權(quán)重?fù)p失權(quán)重EGramPool10.50.5EGramPool20.50.5EDCorPool210.5×10-4EDivPool21-1×10-4ESmoothConv11-0.75×10-2
表2 算法運(yùn)行效率對(duì)比
圖4 γ取不同值的結(jié)果
圖5 總損失函數(shù)誤差
實(shí)驗(yàn)基于Matlab2015進(jìn)行編程,使用MatConvnet[8]來(lái)搭建CNN框架實(shí)現(xiàn)本課題的算法。硬件環(huán)境主要是NVIDIA GeForce GTX 750,8G RAM,Intel I7-4790。