基于變分光流估計(jì)的肺部4D-CT圖像超分辨率重建

耿鳳歡 劉 慧 郭 強(qiáng) 尹義龍

1(山東財(cái)經(jīng)大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院 濟(jì)南 250014) 2(山東省數(shù)字媒體技術(shù)重點(diǎn)實(shí)驗(yàn)室 濟(jì)南 250014) 3 (山東大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院 濟(jì)南 250013) (gfh0928@outlook.com)

基于變分光流估計(jì)的肺部4D-CT圖像超分辨率重建

耿鳳歡1,2劉 慧1,2郭 強(qiáng)1,2尹義龍1,3

1(山東財(cái)經(jīng)大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院 濟(jì)南 250014)2(山東省數(shù)字媒體技術(shù)重點(diǎn)實(shí)驗(yàn)室 濟(jì)南 250014)3(山東大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院 濟(jì)南 250013) (gfh0928@outlook.com)

由于受到掃描時(shí)間和照射劑量的限制，肺部4D-CT數(shù)據(jù)中縱向采樣率遠(yuǎn)小于面內(nèi)采樣率.為了得到更高質(zhì)量的肺部圖像，從醫(yī)學(xué)圖像固有的自相似性出發(fā)，提出了一種基于局部和全局相結(jié)合的變分光流估計(jì)的圖像序列超分辨率重建技術(shù)，用于提高4D-CT圖像重建質(zhì)量.首先，構(gòu)建了一個(gè)用于求解肺部4D-CT不同相位圖像之間的光流場(chǎng)的變分光流模型；然后，利用快速交替方向乘子法求解該模型，得到不同相位圖像之間的光流場(chǎng)；最后，基于光流場(chǎng)，并利用非局部迭代反投影超分辨率重建算法，實(shí)現(xiàn)了高分辨率肺部圖像的重建.實(shí)驗(yàn)結(jié)果表明：與已有算法相比，本方法在增強(qiáng)圖像紋理結(jié)構(gòu)的同時(shí)更好地保留了圖像的輪廓.

4D-CT圖像；超分辨率重建；光流估計(jì)；交替方向乘子法；迭代反投影

放射治療是醫(yī)學(xué)上治療肺癌最常用的手段.肺部4D-CT圖像因其能夠提供必要的呼吸運(yùn)動(dòng)信息來引導(dǎo)醫(yī)生進(jìn)行精確地放射治療，從而得到越來越多的關(guān)注.肺部4D-CT圖像是由呼吸運(yùn)動(dòng)周期中的多個(gè)相位各自對(duì)應(yīng)的3D-CT圖像按照獲取時(shí)間的先后順序排序得到的，能夠包含呼吸周期內(nèi)肺部器官以及腫瘤區(qū)域的空間結(jié)構(gòu)(可通過圖像分割算法[1]得到)以及運(yùn)動(dòng)信息.這些信息可以輔助醫(yī)生精確定位靶區(qū)，有助于病人的放射治療.但是，眾所周知CT圖像的獲取伴隨著高劑量照射[2-3]，為了降低對(duì)病人的輻射量，往往只能降低沿縱向的采樣，從而導(dǎo)致肺部4D-CT圖像的層間分辨率遠(yuǎn)低于層內(nèi)分辨率.在觀察3D肺部圖像的冠狀圖時(shí)，為了還原真實(shí)肺部形狀，通常會(huì)將3D圖像沿Z軸方向進(jìn)行插值放大.傳統(tǒng)插值算法會(huì)使圖像變模糊.因此，如何提高4D-CT圖像超分辨率重建圖像質(zhì)量是亟待解決的問題.

目前，研究者們常利用3類算法實(shí)現(xiàn)4D-CT圖像的超分辨率重建：傳統(tǒng)插值算法、投影重建以及基于字典學(xué)習(xí)的方法.

1) 傳統(tǒng)插值算法.如雙三次插值(Bicubic)算法[4]，是根據(jù)插值公式，利用鄰域內(nèi)已知像素值得到待插值點(diǎn)像素值，其實(shí)現(xiàn)簡(jiǎn)單且運(yùn)行速度快，但由于插值過程中未考慮圖像的各向異性特性，導(dǎo)致放大后圖像的輪廓和紋理比較模糊，降低了圖像質(zhì)量.

2) 投影重建算法.大部分是利用4D-CT不同相位間的互補(bǔ)信息來提高當(dāng)前相位圖像的分辨率，如基于全搜索(full search)的運(yùn)動(dòng)估計(jì)超分辨率重建方法[5]利用迭代反投影思想及全搜索得到的不同相位圖像間的仿射變換信息重建出高分辨率圖像，但是該方法的缺點(diǎn)是對(duì)不同相位圖像間運(yùn)動(dòng)估計(jì)的精度不高；此外，還有基于配準(zhǔn)的凸集投影算法[6]，這類算法的精度很大程度取決于光流估計(jì)誤差的大小.

3) 基于字典學(xué)習(xí)的超分辨重建方法[7-9].它們需要大量的訓(xùn)練數(shù)據(jù)集和較長(zhǎng)運(yùn)算時(shí)間，很難滿足實(shí)時(shí)性的需求.

綜合比較上述3類方法，其中投影重建類算法能夠在較短的計(jì)算時(shí)間內(nèi)得到可靠的超分辨率圖像，較好地迎合了4D-CT超分辨率重建的需求.常用的投影包括迭代反投影算法[10]、凸集投影算法[11]以及最大后驗(yàn)估計(jì)算法[12].其中迭代反投影算法因收斂速度快得到了廣泛的應(yīng)用.

因此，為了提高重建圖像的質(zhì)量和速度，我們采用迭代反投影方法，充分利用4D-CT不同相位間的互補(bǔ)信息實(shí)現(xiàn)快速的CT圖像超分辨率重建，并在2個(gè)方面進(jìn)行了創(chuàng)新：1)鑒于影響超分辨率重建算法效果的主要因素是圖像配準(zhǔn)的精度，本文采用了局部和全局相結(jié)合(combined local and global, CLG)的變分光流估計(jì)模型獲得更加精準(zhǔn)稠密的光流場(chǎng)，并采用快速交替方向乘子法[13-14]實(shí)現(xiàn)光流模型的快速求解；2)采用基于邊緣的非局部迭代反投影算法，有效增強(qiáng)超分辨率圖像中的邊緣和紋理細(xì)節(jié).

本文算法主要包含2個(gè)步驟：1)將待重建相位圖像作為參考相位，采用CLG變分光流估計(jì)算法獲取其他相位與參考相位之間的光流場(chǎng)；2)利用其他相位圖像的互補(bǔ)信息，采用改進(jìn)的非局部迭代反投影算法迭代求解重建圖像.

1 光流模型及求解

CLG變分光流估計(jì)模型是用于計(jì)算視頻中相鄰幀間光流場(chǎng)的模型，4D-CT數(shù)據(jù)中不同相位間光流場(chǎng)亦可用此模型求解.本節(jié)主要介紹針對(duì)肺部4D-CT圖像的CLG變分光流估計(jì)模型，以及利用快速交替方向乘子算法求解該光流模型的過程.

1.1CLG變分光流模型

光流估計(jì)[15]是利用圖像序列中參考幀與當(dāng)前幀之間像素值數(shù)據(jù)的時(shí)域相關(guān)性來確定當(dāng)前幀中各像素位置相對(duì)于參考幀的“運(yùn)動(dòng)”.光流計(jì)算模型都是基于光流基本等式，通過引入光流約束條件來實(shí)現(xiàn)光流求解的.CLG變分光流估計(jì)[16]是一種對(duì)噪聲具有魯棒性且更加準(zhǔn)確的光流計(jì)算模型.首先，該模型中引入了對(duì)光流的l1范數(shù)約束，從而保證了模型對(duì)噪聲的魯棒性；其次，模型引入了對(duì)數(shù)據(jù)項(xiàng)的雙邊濾波操作，以減小相鄰像素之間的相互影響；最后，該模型基于圖像局部特性實(shí)現(xiàn)了光流的各向異性傳播，有效地抑制光流過傳播帶來的負(fù)影響.下面構(gòu)建適用于肺部4D-CT數(shù)據(jù)不同相位圖像間光流估計(jì)的計(jì)算模型.

首先，基于光流場(chǎng)基本理論，將肺部4D-CT圖像看作肺部圖像I隨時(shí)間發(fā)生形變而生成的圖像序列，也稱為圖像一致性假設(shè)，即呼吸運(yùn)動(dòng)過程中肺部圖像各個(gè)像素位置亮度沿著運(yùn)動(dòng)軌跡保持不變，即：

I(x(t),y(t),t)=C.

(1)

其次，引入對(duì)光流的l1范數(shù)約束和各項(xiàng)異性擴(kuò)散約束，以及對(duì)數(shù)據(jù)項(xiàng)的雙邊濾波操作，得到肺部4D-CT圖像的CLG變分光流估計(jì)模型：

(2)

(3)

其中，X=(x,y)為像素點(diǎn)坐標(biāo)，U0是光流場(chǎng)的初始估計(jì).數(shù)據(jù)殘差項(xiàng)的局部約束使得鄰域像素的光流會(huì)影響中心像素點(diǎn)光流的計(jì)算，為了確保在得到稠密光流場(chǎng)的同時(shí)弱化鄰域像素的影響，最有效的方法是引入對(duì)數(shù)據(jù)殘差項(xiàng)的濾波操作.濾波算法有雙邊濾波[18-19]以及基于奇異值分解的濾波算法[20-21]等，其中雙邊濾波算法實(shí)現(xiàn)簡(jiǎn)單且能滿足需求，故本文采用雙邊濾波對(duì)數(shù)據(jù)殘差項(xiàng)進(jìn)行約束.此外，在CLG變分光流估計(jì)模型的求解過程中，光流場(chǎng)的各向同性傳播會(huì)造成光流的過傳播，從而降低光流場(chǎng)的精度.本文中CLG變分光流模型通過加入光流的各向異性濾波操作，有效地抑制了光流的過傳播.

1.2求解光流模型

為了實(shí)現(xiàn)3D-CT序列圖像間光流場(chǎng)的快速求解，我們采用快速交替方向乘子法(alternating direc-tion method of multipliers, ADMM)算法求解本文中的變分光流模型.ADMM是一種求解優(yōu)化問題的計(jì)算框架，通過將全局問題分解成多個(gè)相互獨(dú)立、較小且易于求解的局部子問題，之后協(xié)調(diào)子問題的解得到全局問題的解.ADMM算法結(jié)合了增廣拉格朗日方法以及求解凸優(yōu)化問題的對(duì)偶分解策略的優(yōu)勢(shì)[13]，實(shí)現(xiàn)了對(duì)凸優(yōu)化問題的快速求解.定義Ψi(X)為張量X在第i維度上的融合算子，光流模型可轉(zhuǎn)化成等價(jià)問題：

s.t.U=Vi, 1≤i≤2，

(4)

其中，Vi(1≤i≤2)是松弛變量.為了求解式(4)中的條件極值問題，本文將該約束問題轉(zhuǎn)換成無約束的增廣拉格朗日函數(shù)：

(5)

為了求解式(5)，本文采用了ADMM算法，通過依次迭代更新U,Vi,Wi(i=1,2)來最小化式(5)中的目標(biāo)函數(shù)，得到光流U的局部最優(yōu)值.U,Vi,Wi(i=1,2)的更新規(guī)則：

具體迭代步驟為

Step1. 固定Vi和Wi，更新Uk+1：

(6)

最優(yōu)解為

(7)

Step2. 固定U和Wi，求解Vi,i=1,2

(8)

式(8)是可分解的，我們可以獨(dú)立地計(jì)算Vi,i=1,2，且可等價(jià)寫成:

(9)

(10)

其中，ν是張量Vi所包含的一維張量，t是Γi中與ν相對(duì)應(yīng)的一維張量，k是ν的維數(shù)，di是D中與ν相對(duì)應(yīng)一維張量的第i個(gè)元素.式(10)是一維全變分問題有很多算法可以快速求解.并行計(jì)算的引入使得Vi的計(jì)算可以在毫秒時(shí)間內(nèi)完成.

Step3. 計(jì)算Wi,i=1,2

(11)

該迭代求解過程在原始?xì)埐詈蛯?duì)偶?xì)埐頪13]滿足特定條件時(shí)結(jié)束.

文獻(xiàn)[14]將多維變分問題轉(zhuǎn)換成標(biāo)準(zhǔn)的l1范數(shù)問題，并證明了轉(zhuǎn)換后的問題滿足收斂的約束條件.而本文中光流求解算法僅是該算法中的特例(二維變分問題)，因此，可證明本文光流求解過程是收斂的.

2 改進(jìn)的非局部迭代反投影超分辨率重建

2.1迭代反投影算法

(12)

(13)

其中，Tk(·)是指圖像的仿射變換；h是點(diǎn)擴(kuò)散函數(shù)；↓s為下采樣算子.反投影重建過程可表示為

(14)

其中，T-1k是Tk的逆；p是反投影算子，p的取值將影響算法迭代時(shí)間；↑s表示上采樣算子.

2.2改進(jìn)的非局部迭代反投影算法

為了更好地利用不同幀圖像間的互補(bǔ)信息，我們結(jié)合迭代反投影算法和非局部迭代反投影算法提出了改進(jìn)的非局部迭代反投影，實(shí)現(xiàn)4D-CT圖像超分辨率重建.本文算法在2個(gè)方面具有創(chuàng)新性：1)對(duì)IBP算法中的投影重建誤差進(jìn)行非局部均值濾波，來保留圖像中的高頻細(xì)節(jié)；2)對(duì)反投影重建過程中的圖像融合策略進(jìn)行了改進(jìn).

IBP算法在反投影重建過程中采用了灰度平均的像素級(jí)圖像數(shù)據(jù)融合策略，該策略簡(jiǎn)單且便于實(shí)現(xiàn)，但在輸入低分辨率圖像間光流場(chǎng)估計(jì)誤差較大的情況下會(huì)嚴(yán)重影響超分辨率重建結(jié)果的質(zhì)量.基于此，本文利用數(shù)據(jù)加權(quán)平均法進(jìn)行融合，根據(jù)低分辨率圖像光流估計(jì)誤差確定該圖像的重建誤差在反投影過程中的權(quán)重.這樣便可以很好地抑制光流估計(jì)誤差對(duì)重建結(jié)果的影響，從而得到高質(zhì)量的高分辨率圖像.修改后的反投影過程可表示為

(15)

(16)

其中,t是控制核衰變速度的參數(shù).

我們將4D-CT圖像所包含的所有相位圖像作為已知的低分辨率圖像集合，從中選出某一相位圖像作為浮動(dòng)圖像.首先，用CLG變分光流估計(jì)算法得到浮動(dòng)圖像與其余相位圖像間的光流場(chǎng)；其次，結(jié)合該光流場(chǎng)以及圖像放大倍數(shù)，得到圖像間的仿射變換；最后，采用改進(jìn)的非局部迭代反投影算法重建出高分辨率圖像.算法1給出了本文方法的實(shí)現(xiàn)過程.

算法1. 基于變分光流估計(jì)的4D-CT圖像超分辨率重建算法.

輸出：高分辨率圖像Ih.

①U←利用式(2)得到Ik與參考幀Il間的光流場(chǎng)；

②Tk←由光流場(chǎng)U和倍數(shù)l得到仿射變換；

③ do

④e(n)←0； /*e(n)為第n次迭代中的誤差*/

⑤ fork=1 toK

⑨ end for

3 光流模型中參數(shù)問題

當(dāng)待配準(zhǔn)圖像與基準(zhǔn)間存在較大位移時(shí)，式(1)圖像灰度一致性假設(shè)會(huì)造成估計(jì)結(jié)果誤差較大.為此，我們?cè)谇蠼膺^程中，利用由粗到細(xì)的圖像金字塔方法提高光流估計(jì)精度.式(2)中各項(xiàng)異性擴(kuò)散因子D的參數(shù)設(shè)置為典型值α=5,β=1/2，權(quán)重系數(shù)λ=1 000.增廣拉格朗日函數(shù)中的懲罰項(xiàng)參數(shù)設(shè)置為ρ=10.由光流求解過程可知，雙邊濾波器的鄰域大小的選取對(duì)光流估計(jì)結(jié)果影響最大.因此，我們用6種不同大小濾波器并分別對(duì)2幅圖像間存在大位移和小位移的情況進(jìn)行了對(duì)比實(shí)驗(yàn)，如圖1所示.

Fig. 1 The influence of filter size on the accuracy of optical flow圖1 濾波器大小對(duì)光流精度的影響

圖1中，橫軸表示濾波器鄰域大小，縱軸表示圖像插值誤差[24].插值誤差是一種光流估計(jì)效果評(píng)價(jià)標(biāo)準(zhǔn)，其值為真實(shí)圖像與基于光流補(bǔ)償?shù)玫焦烙?jì)圖像間的均方差.由圖1可知，當(dāng)濾波器大小為3×3，5×5時(shí)，對(duì)于大位移圖像間的光流估計(jì)效果較差；當(dāng)濾波器大小為9×9,11×11,13×13時(shí)，雖然對(duì)圖像間存在大位移的情況下光流估計(jì)結(jié)果較精確，但計(jì)算時(shí)間較長(zhǎng)且對(duì)圖像間存在小位移情況估計(jì)偏差變大了，這是由于濾波器的鄰域變大，使得影響當(dāng)前中心像素的鄰域像素變多，從而影響估計(jì)結(jié)果；當(dāng)濾波器大小為7×7時(shí)，能夠很好地處理大位移光流估計(jì)問題，且處理時(shí)間也能達(dá)到要求.故實(shí)驗(yàn)中將濾波器的大小設(shè)置為7×7.

圖2顯示了使用基于CLG變分光流估計(jì)算法對(duì)不同相位中的3幅冠狀面圖像進(jìn)行光流估計(jì)的2個(gè)實(shí)例.這2幅冠狀面圖像是從同一幅4D-CT圖像中選取不同相位上的冠狀面圖像圖2(a)、圖2(b)和圖2(f)，其中圖2(a)作為參考幀圖像，圖2(b)和圖2(f)分別作為圖像，使用CLG變分光流估計(jì)算法分別得到圖2(a)和圖2(b)之間的光流場(chǎng)圖2(c)以及圖2(a)和圖2(f)之間的光流場(chǎng)圖2(g).利用得到的光流場(chǎng)圖2(c)對(duì)圖2(a)進(jìn)行運(yùn)動(dòng)補(bǔ)償后得到圖2(d)；利用得到的光流場(chǎng)圖2(g)對(duì)圖2(a)進(jìn)行運(yùn)動(dòng)補(bǔ)償后得到圖2(h).圖2(e)是圖2(b)和圖2(d)之間的差圖像，圖2(i)是圖2(f)和圖2(h)之間的差圖像.觀察可知，圖2(b)和圖2(f)中肺部相對(duì)于圖2(a)整體向下運(yùn)動(dòng)且圖2(a)和圖2(b)之間的位移較小、圖2(a)和圖2(f)之間的位移較大，這與得到光流場(chǎng)一致.結(jié)果可見，圖2(b)和圖2(d)以及圖2(f)和圖2(h)的差異很小，說明CLG變分光流估計(jì)能準(zhǔn)確估計(jì)圖像.間的位移場(chǎng)，有助于更好地進(jìn)行超分辨率重建.

Fig. 2 The solved optical flow field and the difference between the original image and the compensated according to the optical flow field圖2 光流場(chǎng)、補(bǔ)償圖像以及誤差圖像

4 實(shí)驗(yàn)結(jié)果對(duì)比與分析

實(shí)驗(yàn)所用數(shù)據(jù)采用的是德克薩斯大學(xué)安德森癌癥中心DIR-lab實(shí)驗(yàn)室*https://www.dir-lab.com/Referenceta.html提供的公開的標(biāo)準(zhǔn)肺部4D-CT數(shù)據(jù)集[25]，文獻(xiàn)[6-7]均采用了該數(shù)據(jù)集.我們從中選取了5組4D-CT數(shù)據(jù)(數(shù)據(jù)1～數(shù)據(jù)5).每組數(shù)據(jù)含有呼吸運(yùn)動(dòng)過程中的10個(gè)不同相位，包括呼氣末端相位和吸氣末端相位.數(shù)據(jù)的層內(nèi)像素尺寸范圍從(0.97×0.97)mm2到(1.16×1.16)mm2，層厚為2.5mm.下面將分別展示CLG變分光流估計(jì)和最終的超分辨率重建結(jié)果.

4.1實(shí)驗(yàn)對(duì)比

實(shí)驗(yàn)將本文算法與雙三次插值(Bicubic)算法、基于全搜索運(yùn)動(dòng)估計(jì)(Full Search)的超分辨算法以及非局部迭代反投影(NLIBP)重建算法進(jìn)行了對(duì)比，以評(píng)估本文算法在邊界保持方面的效果，并論證光流估計(jì)對(duì)超分辨重建效果的影響.

Bicubic算法利用待插值像素點(diǎn)鄰域內(nèi)的16個(gè)像素值得到插值點(diǎn)的像素值，實(shí)現(xiàn)了圖像的放大.但是，該方法未考慮圖像的各向異性特性，導(dǎo)致放大后的圖像輪廓和紋理比較模糊.Full Search算法通過基于圖像塊的全搜索算法對(duì)低分辨率序列進(jìn)行運(yùn)動(dòng)估計(jì)，再利用迭代反投影算法重建出高分辨率圖像，但是，由于全搜索運(yùn)動(dòng)估計(jì)精度較低，不能精確引導(dǎo)投影重建過程，造成重建結(jié)果質(zhì)量較差.NLIBP算法利用低分辨率圖像自身的非局部冗余信息引導(dǎo)重建誤差的反向投影，很好地抑制了高對(duì)比度邊緣的“振鈴”現(xiàn)象，但是，該方法僅利用了單幅低分辨率圖像的自身信息，導(dǎo)致重建圖像缺失部分細(xì)節(jié)信息.相比而言，本文算法在提高光流估計(jì)精度的基礎(chǔ)上，改進(jìn)了非局部反投影算法，更好地利用了不同相位圖像間的互補(bǔ)信息，從而得到具有更多細(xì)節(jié)信息且結(jié)構(gòu)清晰的高分辨率圖像.我們利用圖3和圖4中的對(duì)比結(jié)果對(duì)上述方法進(jìn)行驗(yàn)證.

Fig. 3 Visual comparison of reconstruction results of coronal images圖3 不同算法得到的冠狀圖超分辨率圖像視覺效果對(duì)比

Fig. 4 Visual comparison of reconstruction results of sagittal images圖4 不同算法得到的矢狀面超分辨率圖像視覺效果對(duì)比.

圖3展示了采用不同算法對(duì)數(shù)據(jù)2進(jìn)行冠狀面圖像重建的結(jié)果圖像.由于數(shù)據(jù)2中層內(nèi)像素尺寸為(1.16×1.16)mm2且層厚2.5mm，為了使軸向采樣率和層內(nèi)采樣率相同，本文將數(shù)據(jù)2冠狀圖上采樣倍數(shù)設(shè)置為2.15.圖3中列圖像分別是利用Bicubic算法、Full Search算法、NLIBP算法以及本文算法得到的冠狀面圖像超分辨率結(jié)果.前3行圖像是分別對(duì)應(yīng)數(shù)據(jù)2中不同相位(相位5、相位8和相位1)的超分辨率重建結(jié)果.為了更好地對(duì)比不同算法重建效果，第4行給出了第3行圖像中局部細(xì)節(jié)的放大圖.對(duì)比前3行圖像可以發(fā)現(xiàn)，當(dāng)以不同相位圖像作為基準(zhǔn)圖像時(shí)得到的超分辨率圖像間存在差異，但在圖像質(zhì)量方面相差無幾.對(duì)比放大后的肺實(shí)質(zhì)周圍的血管與組織邊緣處的細(xì)節(jié)可以看出，本文算法重建出的圖像具有更清晰的結(jié)構(gòu)，同時(shí)邊緣和細(xì)節(jié)信息也得到加強(qiáng).

圖4則給出采用不同算法對(duì)數(shù)據(jù)2進(jìn)行矢狀面圖像重建得到的結(jié)果圖像.圖4的圖片排列方式與圖3一致.對(duì)比放大后的血管與組織邊緣處的細(xì)節(jié)可以看出，Bicubic算法重建效果最差；較Full Search方法，NLIBP算法與本文算法重建出的圖像具有更清晰的結(jié)構(gòu)；較NLIBP算法，本文方法重建出的圖像中細(xì)節(jié)信息更加明顯.

為了更好地評(píng)價(jià)重建結(jié)果圖像的清晰度，我們采用平均梯度作為量化標(biāo)準(zhǔn).圖像邊緣附近的灰度值變化率較大，這種變化率大小可以用來表示圖像的清晰度.圖像平均梯度就是用圖像多維方向上灰度變化率之和來表征圖像的相對(duì)清晰度，值越大圖像越清晰：

(17)

表1列出了分別利用Bicubic算法、Full Search算法、NLIBP算法以及本文算法對(duì)5幅4D-CT圖像選擇不同相位圖像為浮動(dòng)圖像，對(duì)不同平面進(jìn)行超分辨率重建，然后計(jì)算得到所有重建圖像平均梯度值的平均值.對(duì)比表1中數(shù)據(jù)可知，較前3種方法，本文方法得到的平均梯度值更高，圖像清晰度明顯增強(qiáng).

Table 1 Average Gradient Using Four Different Algorithms

Fig. 5 Comparison of average gradient of super reconstruction results圖5 不同算法得到的高分辨率結(jié)果對(duì)比

此外，利用山東省千佛山醫(yī)院提供的10組4D-CT數(shù)據(jù)對(duì)相關(guān)算法進(jìn)行了對(duì)比分析.圖5給出了不同算法得到的高分辨率圖像平均梯度，可以看出，對(duì)于10組實(shí)驗(yàn)數(shù)據(jù)，本文算法得到的超分辨率重建結(jié)果均好于前3種算法.

4.2結(jié)果分析

從圖3和圖4中的細(xì)節(jié)放大圖中可以看出，全搜索算法由于得到的光流場(chǎng)不夠精確，導(dǎo)致重建圖像中出現(xiàn)一些噪聲.本文算法利用基于CLG變分光流估計(jì)算法計(jì)算得到比全搜索估計(jì)更加精確和稠密的光流場(chǎng)，所以超分辨結(jié)果要優(yōu)于全搜索算法.由于非局部迭代反投影只利用低分辨率圖像自身的冗余信息，導(dǎo)致超分辨率結(jié)果圖像在具有較少信息冗余的器官組織邊緣出現(xiàn)了模糊.而本文算法則是很好地利用不同幀圖像間的互補(bǔ)信息，在保留細(xì)節(jié)信息的同時(shí)提高了重建圖像的清晰度.對(duì)比表1中5份4D-CT數(shù)據(jù)重建結(jié)果圖像的平均梯度值可知，本文算法得到的平均梯度值都是最大的，這表明本文算法更好地保留了圖像的細(xì)節(jié)信息.

5.3算法時(shí)間效率分析

本文利用ADMM算法求解光流模型的過程中，主要工作是迭代更新U,Vi,Wi(i=1,2).設(shè)實(shí)驗(yàn)圖像包含的像素?cái)?shù)為n，迭代次數(shù)為k1，則在光流求解過程中，更新U的時(shí)間復(fù)雜度為O(n)，更新Vi(i=1,2)的時(shí)間復(fù)雜度為O(m)(其中m為實(shí)驗(yàn)圖像像素的列數(shù))，更新Wi(i=1,2)的時(shí)間復(fù)雜度為O(n).由于m

實(shí)驗(yàn)用C++語言實(shí)現(xiàn)了本文中光流求解算法，并在配備英特爾i5處理器(主頻為3.2 GHz)的主機(jī)上進(jìn)行了測(cè)試.對(duì)于分辨率為256×99的浮動(dòng)圖像和基準(zhǔn)圖像，平均光流計(jì)算時(shí)間為0.323 s，極大地縮短了光流計(jì)算時(shí)間.因此，對(duì)于10幅圖像組成的實(shí)驗(yàn)圖像序列，使用本文超分辨率算法可以在30 s內(nèi)得到超分辨結(jié)果.

5 總 結(jié)

本文提出了一種基于CLG變分光流估計(jì)的超分辨率重建算法，用于提高4D-CT冠狀面圖像分辨率.文中，我們將不同相位的3D-CT圖像看作不同“幀”圖像，基于CLG變分光流模型，并利用快速ADMM算法迅速收斂的性質(zhì)，實(shí)現(xiàn)了浮動(dòng)“幀”圖像與基準(zhǔn)“幀”圖像間光流場(chǎng)的快速計(jì)算；然后，基于得到的光流場(chǎng)，改進(jìn)了非局部迭代反投影算法，并利用不同相位間的互補(bǔ)信息增強(qiáng)重建了圖像的結(jié)構(gòu)細(xì)節(jié)信息.實(shí)驗(yàn)結(jié)果表明：本文算法不僅有效增強(qiáng)了圖像的紋理結(jié)構(gòu)，而且能夠更好地保留圖像中的輪廓信息.

在未來的研究工作中，我們將研究實(shí)現(xiàn)光流模型及其求解過程中參數(shù)的自適應(yīng)設(shè)定，以進(jìn)一步提高光流估計(jì)的精度和超分辨率重建質(zhì)量.此外，我們還將應(yīng)用GPU并行計(jì)算技術(shù)來縮短本文算法的運(yùn)行時(shí)間，用于滿足實(shí)時(shí)性需求.

[1] Liu Hui, Zhang Caiming, Deng Kai, et al. Research on fast FCM pulmonary nodule segmentation algorithm using improved self-adaption[J]. Journal of Computer-Aided Design and Computer Graphics, 2014, 26(10): 1727-1736 (in Chinese)(劉慧, 張彩明, 鄧凱, 等. 改進(jìn)局部自適應(yīng)的快速FCM肺結(jié)節(jié)分割方法[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2014, 26(10): 1727-1736)

[2] Khan F, Bell G, Antony J, et al. The use of 4D-CT to reduce lung dose: A dosimetric analysis[J]. Medical Dosimetry, 2010, 34(4): 273-278

[3] Li Tianfang, Schreibmann E, Thorndyke B, et al. Radiation dose reduction in four-dimensional computed tomography[J]. Medical Physics, 2005, 32(12): 3650-3660

[4] Carey W, Chuang D, Hemami S. Regularity-preserving image interpolation[J]. IEEE Trans on Image Processing, 1997, 8(9): 1293-1297

[5] Su Heng, Zhou Jie, Zhang Zhihao. Survey of super-resolution image reconstruction methods[J]. Acta Auto-matica Sinica, 2013, 39(8): 1202-1213 (in Chinese)(蘇衡, 周杰, 張志浩. 超分辨率圖像重建方法綜述[J]. 自動(dòng)化學(xué)報(bào), 2013, 39(8): 1202-1213)

[6] Zhang Yu, Wu Xiuxiu, Yang Wei, et al. Super-resolution reconstruction for 4D computed tomography of the lung via the projections onto convex sets approach[J]. Medical Physics, 2014, 41(11): No.111917

[7] Zhang Yu, Wu Guorong, Yap P, et al. Hierarchical patch-based sparse representation—A new approach for resolution enhancement of 4D-CT lung data[J]. IEEE Trans on Medical Imaging, 2012, 31(11): 1993-2005

[8] Li Min, Chen Jian, Le Xiang, et al. Super-resolution based on sparse dictionary coding[J]. Journal of Software, 2012(5): 1315-1324 (in Chinese)(李民, 程建, 樂翔, 等. 稀疏字典編碼的超分辨率重建[J]. 軟件學(xué)報(bào), 2012(5): 1315-1324)

[9] Yang Shuaifeng, Zhao Ruizhen. Image super-resolution reconstruction based on low-rank matrix and dictionary learning[J]. Journal of Computer Research and Development, 2016, 53(4): 884-891 (in Chinese)(楊帥鋒, 趙瑞珍. 基于低秩矩陣和字典學(xué)習(xí)的圖像超分辨率重建 DLOF[J]. 計(jì)算機(jī)研究與發(fā)展, 2016, 53(4): 884-891)

[10] Irani M, Peleg S. Motion analysis for image enhancement: resolution, occlusion, and transparency[J]. Journal of Visual Communication and Image Representation, 1993, 4(4): 324-347

[11] Stark H, Oskoui P. High-resolution image recovery from image-plane arrays using convex projections[J]. Journal of the Optical Society of America Optics and Image Science, 1989, 6(11): 1715-1726

[12] Levitan E, Herman G. A maximum a posteriori probability expectation maximization algorithm for image reconstruction in emission tomography[J]. IEEE Trans on Medical Imaging, 1987, 6(3): 185-192

[13] Boyd S, Parikh N, Chu E, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers[J]. Foundations and Trends in Machine Learning, 2011, 3(1): 89-97

[14] Yang Sen, Wang Jie, Fan Wei, et al. An efficient ADMM algorithm for multidimensional anisotropic total variation regularization problem[C] //Proc of the 19th ACM SIGKDD Int Conf on Knowledge Discovery and Data Mining. New York: ACM, 2013: 641-649

[15] Horn B, Schunck B. Determining optical flow[J]. Artificial Intelligence, 1981, 17(1/2/3): 185-203

[16] Marius D, Sergiu N. Total variation regularization of local-global optical flow[C] //Proc of the 14th Int IEEE Conf on Intelligent Transportation Systems. Piscataway, NJ: IEEE, 2011: 318-323

[17] Perona P, Malik J. Scale space and edge detection using anisotropic diffusion[J]. IEEE Trans on Pattern Analysis and Machine Intelligence, 1990, 7(12): 629-639

[18] Tomasi C, Manduchi R. Bilateral filtering for gray and color images[C] //Proc of the 6th Int Conf on Computer Vision(ICCV). Los Alamitos, CA: IEEE Computer Society, 1998: 839-846

[19] Durand F, Dorsey J. Fast bilateral filtering for the display of high-dynamic-range images[J]. ACM Trans on Graphics, 2002, 21(3): 257-266

[20] Guo Qiang, Zhang Caiming, Zhang Yunfeng. An Efficient SVD-Based Method for Image Denoising[J]. IEEE Trans on Circuits and Systems for Video Technology, 2016, 26(5): 868-880

[21] Guo Qiang, Zhang Caiming, Zhang Yunfeng, et al. Low-rank image denoising based on minimum variance estimator[J]. Journal of Computer-Aided Design and Computer Graphics, 2015, 27(12): 2237-2246 (in Chinese)(郭強(qiáng), 張彩明, 張?jiān)品? 等. 基于最小方差估計(jì)的圖像低秩去噪[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2015, 27(12): 2237-2246)

[22] Liu Jun, Yuan Lei, Ye Jieping. An efficient algorithm for a class of fused lasso problem[C] //Proc of the 19th ACM SIGKDD Int Conf on Knowledge Discovery and Data Mining. New York: ACM, 2010: 323-332

[23] Dong Weisheng, Zhang Lei, Shi Guangming, et al. Nonlocal back-projection for adaptive image enlargement[C] //Proc of the 16th IEEE Int Conf on Image. Piscataway, NJ: IEEE, 2009: 348-352

[24] Baker S, Scharstein D, Lewis J, et al. A dataset and evaluation methodology for optical flow[J]. International Journal of Computer Vision, 2011, 92(1): 1-31

[25] Castillo R, Castillo E, Guerra R, et al. A framework for evaluation of deformable image registration spatial accuracy using large landmark point sets[J]. Physical in Medicine and Biology, 2009, 54(7): 1849-1870

VariationalOpticalFlowEstimationBasedSuper-ResolutionReconstructionforLung4D-CTImage

Geng Fenghuan1,2, Liu Hui1,2, Guo Qiang1,2, and Yin Yilong1,3

1(CollegeofComputerScienceandTechnology,ShandongUniversityofFinanceandEconomics,Jinan250014)2(ShandongProvincialKeyLaboratoryofDigitalMediaTechnology,Jinan250014)3(CollegeofComputerScienceandTechnology,ShandongUniversity,Jinan250013)

The acquisition of lung 4D computed tomography (4D-CT) data is limited by the scanning time and radiation dose, which leads to the sampling rate in the axial direction is much less than that in the in-plane direction. In order to get better quality of 4D-CT images, based on the inherent self-similarity of medical images, a new method of image sequence super-resolution reconstruction is proposed in this paper. This method uses the local and global variational optical flow estimation to improve the quality of enlarged 4D-CT image. Firstly, we present a combined local and global variational optical flow model, in order to estimate the motion fields (i.e., the optical flow fields) between different phases in the corresponding positions. Then, the optical flow field is obtained by solving the model with the fast alternating direction method of multiplier. Finally, according to the calculated motion fields, we employ the improved non-local iterative back projection (NLIBP) algorithm to reconstruct high resolution lung images. The experimental results have shown that, in both quantification standard and visual perception, this method outperforms non-local iterative back projection algorithm and full search block matching based iterative back projection technique. Furthermore, our method can generate clear edges while enhancing the texture of images.

4D-CT image; super-resolution reconstruction; optical flow estimation; alternating direction method of multipliers; iterative back projection

Geng Fenghuan, born in 1992. Master candidate. His main research interests include image processing and machine learning.

Liu Hui, born in 1978. PhD, professor and master supervisor. Member of CCF. Her main research interests include image processing, computer vision, and machine learning.

Guo Qiang, born in 1979. PhD, associate professor and MSc supervisor. Member of CCF. His main research interests include image processing and computer vision.

Yin Yilong, born in 1972. PhD, professor and PhD supervisor. Senior member of CCF. His main research interests include machine learning, data mining, and computational medicine.

2017-05-23；

：2017-06-20

國(guó)家自然科學(xué)基金項(xiàng)目(61572286，61472220)；山東省重點(diǎn)研發(fā)項(xiàng)目(2014GGX101037)；濟(jì)南市高校自主創(chuàng)新項(xiàng)目(201401216)；山東省高校優(yōu)勢(shì)學(xué)科人才團(tuán)隊(duì)培育項(xiàng)目 This work was supported by the National Natural Science Foundation of China (61572286,61472220), the Project of Shandong Provincial Key Research and Development (2014GGX101037), the Project of Independent Innovation of Universities in Jinan (201401216), and the Fostering Project of Dominant Displine and Talent Team of Shandong Province Higher Education Institutions.

劉慧(liuh_lh@sdufe.edu.cn)

TP391