一種基于齒科序列圖像的超分辨率重建算法

程 斌，汪 偉

(上海理工大學(xué) 光電信息與計(jì)算機(jī)工程學(xué)院，上海 200093)

1 引 言

醫(yī)學(xué)影像在臨床上應(yīng)用廣泛，能夠?yàn)榧膊〉脑\斷提供科學(xué)、直觀的依據(jù)；同時(shí)也能夠結(jié)合臨床的癥狀、化驗(yàn)等結(jié)果，為最終病情的確定和后續(xù)的治療計(jì)劃提供重要參考.但是，由于物理分辨率的限制，成像的質(zhì)量遠(yuǎn)遠(yuǎn)不能滿足醫(yī)學(xué)上的需求.為了獲取更高分辨率影像，常采用減小傳感器像素尺寸或增大成像陣列芯片大小等硬件的方法.但當(dāng)像素尺寸減小時(shí)，加性噪聲幾乎維持不變，有效信號(hào)的能量將隨傳感器像素尺寸成比例減??；而后者高精度光學(xué)傳感器的增大又將會(huì)顯著增加成本，給醫(yī)療設(shè)備的普及帶來重要阻礙.因此，為了規(guī)避硬件發(fā)展的瓶頸，許多研究者深入地分析了圖像的成像過程，構(gòu)建出較為理想的圖像退化模型，并嘗試?yán)猛粓?chǎng)景下獲得的低分辨率序列圖像，通過軟件計(jì)算來重構(gòu)出高分辨率的影像.這一思想得到了快速發(fā)展，且很快成為了目前研究的熱點(diǎn)問題之一.

當(dāng)前，基于空域的超分辨率序列圖像重建算法不斷被提出，包括非均勻插值法、迭代反投影法(iterative back projection,IBP)、凸集投影法(projection onto convex sets,POCS)、最大后驗(yàn)概率法(maximum a posterior，MAP)，POCS-MAP混合法等[1,2].在醫(yī)學(xué)中，Wissam基于兩正交的低分辨率核磁共振(magnetic resonance imaging，MRI)數(shù)據(jù)集和電子計(jì)算機(jī)斷層掃描(computed tomography，CT)數(shù)據(jù)集，在上采樣的插值過程中，利用歐式距離計(jì)算出每個(gè)像素置信度，并結(jié)合區(qū)域和局部信息，構(gòu)造損失因子，修補(bǔ)了IBP算法中的高分辨率圖像初始估計(jì)[3]；Esben Plenge等人考慮不同的數(shù)據(jù)和噪聲模型，根據(jù)MAP法可靈活調(diào)整數(shù)據(jù)保真項(xiàng)(如l1,l2范數(shù))和正則項(xiàng)(如Tikhonov，總變分正則項(xiàng)等)的特點(diǎn)，將IBP、MAP+ Tikhonov正則項(xiàng)、MAP+Lambda權(quán)重正則項(xiàng)3種序列圖像超分辨率重建(multi-image super resolution reconstruction，MISR)方法計(jì)算結(jié)果與通過高分辨率硬件成像設(shè)備直接獲取的MRI圖像進(jìn)行比對(duì)分析，發(fā)現(xiàn)MISR方法能獲得更好的評(píng)價(jià)結(jié)果[4]；Jonghye Woo等結(jié)合Markov隨機(jī)場(chǎng)，提出了一種基于舌頭MRI序列圖像的重建方法[5]；除此之外，Hayit Greenspan對(duì)MISR方法在醫(yī)學(xué)圖像MRI，F(xiàn)MRI(functional MRI)，PET(positron emission tomography)上的應(yīng)用進(jìn)行了回顧與分析，并以大腦影像為例，比較了不同MISR方法的效果與局限性[6].

除了上述提到的MRI圖像外，M.Dirk Robinson等在前人應(yīng)用的基礎(chǔ)上，將MISR擴(kuò)展到乳房X光片(X-ray digital mammography)和光學(xué)相干斷層掃描(optical coherence tomography,OCT)應(yīng)用領(lǐng)域上[7]；Hong等利用光流法配準(zhǔn)后，基于IBP方法完成了超聲圖像的超分辨率重建[8].在2016年，Springer期刊上發(fā)表了一篇《智能優(yōu)化在生物醫(yī)學(xué)中的應(yīng)用》?？痆9]，詳細(xì)介紹了超分辨重建算法的發(fā)展，重點(diǎn)談到了最優(yōu)化方法在醫(yī)學(xué)圖像重建上的應(yīng)用.并結(jié)合癌癥成像數(shù)據(jù)庫，對(duì)采用bicubic插值、IBP、POCS、通過稀疏表示的梯度下降法這四種方法的重建結(jié)果進(jìn)行了分析.

與基于空域的方法相對(duì)應(yīng)，研究者們也提出了基于頻域的超分辨率序列圖像重建算法.在頻域中，成像模型內(nèi)的卷積、幾何運(yùn)動(dòng)等操作都將被轉(zhuǎn)換為線性的計(jì)算，超分辨率重建問題等效為了線性逆向模型的求解，各種頻域變換如離散余弦變換、小波變換、Contourlet變換等在醫(yī)學(xué)上相繼得到了應(yīng)用.Kumar V.R.等人通過對(duì)小波變換后子帶系數(shù)的奇異值估計(jì)和插值法重建了大腦CT、胰腺CT和胸部X光片[10]；Yan Shen等基于Contourlet變換重建了乳房X射線圖，幫助乳腺癌診斷中的鈣化檢測(cè)[11].在2005年，Arthur L.等提出了非下采樣Contourlet變換(the non-subsampled contourlet transform,NSCT)[12]，該方法相對(duì)于小波變換、Contourlet變換具有平移不變性的特點(diǎn).國內(nèi)的楊貞萍，周靖鴻等人分別則采用該變換，選用光流法配準(zhǔn)和矩陣乘法配準(zhǔn)，通過系數(shù)融合策略重建出文本圖像[13]和遙感圖像[14]，這也給該變換在醫(yī)學(xué)圖像上的應(yīng)用提供了可能.

近年來，基于學(xué)習(xí)的超分辨率重建算法[15,16]得到了快速發(fā)展，其中單幀圖像的超分辨率重建(single-image super-resolution reconstruction,SISR)方法較多，比如Yang的稀疏表示模型、Dong的卷積神經(jīng)網(wǎng)絡(luò)SRCNN、Ledig等提出的生成對(duì)抗網(wǎng)絡(luò)SRGAN等[17]，并且有關(guān)研究人員都用在了醫(yī)學(xué)圖像上[18-20]，而基于學(xué)習(xí)的MISR方法相對(duì)較少，其大多在SISR基礎(chǔ)上，在前級(jí)輸入中加入運(yùn)動(dòng)補(bǔ)償處理，或設(shè)計(jì)網(wǎng)絡(luò)中就考慮多幀信息的融合問題，如VSRnet (video super-resolution with convolutional neural networks)[21],VESPCN (video efficient sub-pixel convolution network[22],Dynamic filter generate on network[17]等.

牙科疾病的多樣性和個(gè)體的差異性給學(xué)習(xí)類的重建方法帶來了重要麻煩，存在針對(duì)不同的病患群體和不同的病灶類型，需獲取大量的訓(xùn)練樣本和實(shí)時(shí)效率的問題；而在傳統(tǒng)的MISR方法中，文獻(xiàn)[13]、文獻(xiàn)[14]等提出的配準(zhǔn)方法又只能處理序列圖像間的平移情況；另外口腔圖像細(xì)節(jié)的損失將會(huì)嚴(yán)重影響牙醫(yī)的最終診斷，因此本文提出了一種針對(duì)醫(yī)學(xué)齒科序列圖像的超分辨率重建算法.該算法創(chuàng)新性地將點(diǎn)云中的配準(zhǔn)方法，即迭代最近鄰點(diǎn)法[23,24](iterative closest point，ICP)引入到齒科序列圖像的超分辨率重建過程中，并分別在配準(zhǔn)點(diǎn)集篩選和配準(zhǔn)策略上進(jìn)行了改進(jìn)，增強(qiáng)了算法的魯棒性與精確度，從而可以克服實(shí)際檢測(cè)過程中口腔內(nèi)部組織三維運(yùn)動(dòng)所帶來的配準(zhǔn)難題；其次，算法中采用基于多尺度的NSCT變換來融合齒科序列圖像的有效信息，NSCT變換具有平移不變性的特點(diǎn)，其非下采樣的過程不但可以彌補(bǔ)圖像分解過程中信息的丟失，使邊緣等細(xì)節(jié)特征得以保留，而且通過對(duì)子帶系數(shù)的處理，可以將噪聲信息從圖像信息中有效濾除，從而改善超分辨率齒科圖像的重建結(jié)果，更好地輔助牙醫(yī)診斷病情和制定治療計(jì)劃.

2 基于齒科序列圖像的超分辨率重建

2.1 圖像退化模型

圖像的退化過程是超分辨率重建算法得以實(shí)現(xiàn)的理論基礎(chǔ)和依據(jù).現(xiàn)常用的圖像退化(觀測(cè))模型如圖1所示.

圖1 圖像退化模型[1]Fig.1 Image observation model

上述模型數(shù)學(xué)表達(dá)為：

yk=DBkMkx+nk(1≤k≤p)

(1)

其中x表示高分辨率圖像HR；yk表示第k張低分辨率圖像LR；Mk表示幾何形變矩陣；Bk對(duì)應(yīng)著模糊矩陣，如傳感器、運(yùn)動(dòng)、光學(xué)等模糊，常被模型化成各種點(diǎn)擴(kuò)散函數(shù)(a point spread function, PSF)；D代表下采樣矩陣；nk、p分別表示加性噪聲和低分辨率序列圖像幀數(shù).

超分辨率算法的重建過程，即是圖像退化模型的逆過程.因此，要實(shí)現(xiàn)高分辨率圖像的重建，必須完成如下三個(gè)步驟：矯正低分辨率序列圖像之間的幾何運(yùn)動(dòng)、去除模糊和噪聲、圖像上采樣處理，這也構(gòu)成了本文算法的理論基礎(chǔ).

2.2 齒科序列圖像超分辨率重建算法流程

如圖2所示，齒科序列圖像超分辨率重建步驟可簡(jiǎn)述如下：

圖2 本文重建算法主流程Fig.2 Main flow chart of the proposed reconstruction algorithm

1)讀入從口腔診療設(shè)備獲得的低分辨率齒科序列圖像，并進(jìn)行預(yù)處理；

2)結(jié)合梯度算子計(jì)算邊緣圖像，形態(tài)學(xué)優(yōu)化后，提取出初步特征點(diǎn)集，然后經(jīng)過采樣和邊界點(diǎn)剔除，形成優(yōu)化后的配準(zhǔn)點(diǎn)集，如圖2中虛線框部分所示；

3)設(shè)定魯棒損失函數(shù)，進(jìn)行ICP算法的迭代配準(zhǔn)過程，求得最優(yōu)幾何變換矩陣，并對(duì)序列齒科圖像進(jìn)行矯正；

4)利用軟判決自適應(yīng)插值法，插值放大已配準(zhǔn)的低分辨率序列圖像；

5)采用NSCT變換，提取各圖像的低頻子帶和高頻子帶，對(duì)低頻子帶系數(shù)采用加權(quán)均值融合的策略，對(duì)高頻子帶系數(shù)用硬閾值去噪后，采用取最大絕對(duì)值的融合策略；

6)將NSCT各個(gè)子帶系數(shù)完成融合后，再結(jié)合NSCT反變換，得到最后的超分辨率重建齒科圖像.

2.3 圖像的預(yù)處理

為改善口腔視頻圖像存在的噪聲干擾及光照不均勻現(xiàn)象，首先需對(duì)源圖像進(jìn)行預(yù)處理，即圖像直方圖均衡化：

(2)

式(2)中rk代表歸一化后的離散灰度級(jí)，且0≤rk≤1,k=0,1,2,…,n-1；ni為圖像中第i個(gè)灰度級(jí)的像素?cái)?shù)；n為圖像中的像素總數(shù).

2.4 圖像的配準(zhǔn)

在序列圖像的超分辨率重建方法中，最基本的思想就是整合捕捉到的多幀低分辨率圖像信息來提高空間分辨率，而信息整合的好壞，很大程度上取決于圖像配準(zhǔn)的精度.

本文采用迭代最近鄰點(diǎn)法配準(zhǔn)低分辨率序列圖像，矯正序列圖像間存在的幾何運(yùn)動(dòng).經(jīng)典ICP配準(zhǔn)方法最早由Besl和McKay提出[23]，并在點(diǎn)云的配準(zhǔn)上得到了廣泛應(yīng)用.但該方法對(duì)異常點(diǎn)對(duì)不敏感、干擾點(diǎn)的存在會(huì)引起配準(zhǔn)效果的急劇惡化、點(diǎn)集之間需滿足包含的關(guān)系，而超分辨率重建結(jié)果嚴(yán)重依賴于配準(zhǔn)效果.所以，本文創(chuàng)新性地在點(diǎn)集的篩選和配準(zhǔn)策略上改進(jìn)了ICP算法，以提高低分辨率齒科序列圖像間運(yùn)動(dòng)估計(jì)的精度.

2.4.1 配準(zhǔn)點(diǎn)集提取與優(yōu)化

通過對(duì)多種已有特征點(diǎn)提取方法的比較與分析，本文利用預(yù)處理圖像中的梯度信息和邊緣信息，提取出牙體硬組織的邊緣特征點(diǎn)作為初始配準(zhǔn)點(diǎn)集.

1)牙體硬組織邊緣圖像LRE的計(jì)算

在常用邊緣檢測(cè)算法中，相比于Canny算子,Sobel算子雖然檢測(cè)精度較低，但具有較低的運(yùn)算復(fù)雜度、對(duì)噪聲不敏感.故本文采用Sobel算子提取低分辨率圖像LR的邊緣圖像:

LRE=Sobel(LR)

(3)

2)邊緣圖像的優(yōu)化

對(duì)所獲得的邊緣圖像采用數(shù)學(xué)形態(tài)學(xué)中的膨脹操作，既提高了配準(zhǔn)點(diǎn)集中有效點(diǎn)的數(shù)量，又減少了干擾點(diǎn)在點(diǎn)集中所占的比例，可有效提高配準(zhǔn)的精度.

(4)

(5)

3)初始配準(zhǔn)點(diǎn)集采樣

為了降低計(jì)算量，利用采樣函數(shù)對(duì)初始特征點(diǎn)集E進(jìn)行采樣：

Model=sample(ER)

(6)

Data=sample(ETi)(i=2,3,…,M)

(7)

式(6)和式(7)中，ER表示參考幀的初始配準(zhǔn)點(diǎn)集，ETi為序列中第i待配準(zhǔn)幀的初始配準(zhǔn)點(diǎn)集,M為低分辨率醫(yī)學(xué)序列圖像幀數(shù)，sample()為預(yù)設(shè)的采樣函數(shù)，一般可選為均勻采樣或者隨機(jī)采樣函數(shù).

4)配準(zhǔn)點(diǎn)集的優(yōu)化

在低分辨率序列圖像間的幾何運(yùn)動(dòng)，造成圖像的四周邊界處差異最大，而從這些邊界處提取出的配準(zhǔn)點(diǎn)集間很少存在包含關(guān)系，即會(huì)對(duì)配準(zhǔn)迭代的過程產(chǎn)生不利影響.因此，需要去除圖像邊界處的邊緣點(diǎn)，從而優(yōu)化配準(zhǔn)點(diǎn)集：

Model′=Model(|x-xmid|≤T1,|y-ymid|≤T2)

(8)

Data′=Data(|x-xmid|≤T3,|y-ymid|≤T4)

(9)

其中，T1，T2，T3，T4分別為設(shè)定的參考幀點(diǎn)集與待配準(zhǔn)幀點(diǎn)集坐標(biāo)變化閾值，xmid、ymid分別為圖像中心像素的橫坐標(biāo)與縱坐標(biāo).

2.4.2 ICP算法迭代

基于低分辨率參考幀和待配準(zhǔn)幀完成特征點(diǎn)集提取和優(yōu)化后，得到優(yōu)化點(diǎn)集Model′和Data′.在經(jīng)典的ICP配準(zhǔn)策略中，所有最近鄰點(diǎn)對(duì)的權(quán)重都是相同的；但是，對(duì)于匹配點(diǎn)間歐式距離過大的情況，極大概率是屬于異常干擾點(diǎn).因此，為增加配準(zhǔn)過程的魯棒性，引入魯棒損失函數(shù)來調(diào)整不同匹配點(diǎn)對(duì)所對(duì)應(yīng)的權(quán)值，從而減少異常匹配點(diǎn)對(duì)的干擾.其改進(jìn)后的ICP算法迭代步驟可簡(jiǎn)要描述為：

1)設(shè)置魯棒損失函數(shù).常見的魯棒損失函數(shù)[25]及對(duì)應(yīng)的權(quán)值函數(shù)如下：

-Huber函數(shù)

(10)

(11)

-Tukey′s Biweigh函數(shù)

(12)

(13)

-Cauchy函數(shù)

(14)

(15)

其中r指的是殘差，kHu、kTu和kCa為預(yù)設(shè)值.

2)輸入低分辨率參考幀點(diǎn)集Model′和待配準(zhǔn)幀點(diǎn)集Data′，根據(jù)魯棒損失函數(shù)調(diào)用相應(yīng)的權(quán)值函數(shù)，并設(shè)置初始旋轉(zhuǎn)矩陣和平移向量為單位矩陣E和零向量；

4)以配對(duì)點(diǎn)對(duì)之間的歐式距離為依據(jù)，結(jié)合權(quán)值函數(shù)，確定所有配準(zhǔn)點(diǎn)對(duì)的權(quán)值ωj；

5)根據(jù)上述權(quán)值大于0的配對(duì)點(diǎn)，求得使目標(biāo)函數(shù)Dmin最小化的剛體變換，得出平移參數(shù)和旋轉(zhuǎn)量；

(16)

(17)

其中R(k)，T(k)分別為迭代計(jì)算的旋轉(zhuǎn)矩陣和平移向量.

6)如果變換點(diǎn)集與源點(diǎn)集之間的均方根誤差RMSE小于某一給定閾值,或達(dá)到設(shè)定的最大配準(zhǔn)迭代次數(shù)時(shí)，停止迭代計(jì)算；否則，將變換后的點(diǎn)集作為新的待配準(zhǔn)點(diǎn)集，按照步驟3)至步驟5)繼續(xù)與源點(diǎn)集迭代計(jì)算，求取最優(yōu)的變換參數(shù).

其完整的配準(zhǔn)迭代過程偽代碼[25]如下：

算法.The Modified ICP algorithm

輸出:R(k),T(k)

1.Repeat

2.k=k+1

5. if anyωj>0 then

8.R(k)=R*R(k-1)

9.T(k)=R*T(k-1)+T*

10. else

12.R(k)=R(k-1)

13.T(k)=T(k-1)

14. Terminate iteration

15. end if

16. until convergence

17.ReturnR(k),T(k)

2.5 軟判決自適應(yīng)插值放大

在改進(jìn)ICP配準(zhǔn)法完成齒科序列圖像的幀間配準(zhǔn)后，本文選用在保持圖像空間一致性上有優(yōu)異表現(xiàn)的軟判決自適應(yīng)插值法(soft-decision adaptive interpolation，SAI)[26]提高圖像空間分辨率，有效降低插值后的模糊、振鈴和鋸齒現(xiàn)象.

如圖3所示，黑點(diǎn)代表已知的低分辨率像素點(diǎn)，白點(diǎn)和灰點(diǎn)表示待插值的高分辨率像素點(diǎn).SAI插值算法將先建立圖像的分段自回歸過程(PAR)模型，其中模型參數(shù)分別由圖3(b)中八鄰域像素值、圖3(c)中四鄰域像素值得出.然后，基于指定滑窗中的最小均方誤差估計(jì)出圖3(a)中灰點(diǎn)的像素值；最后，將模型旋轉(zhuǎn)45°，估計(jì)出白點(diǎn)的像素值，從而得到完整放大后的高分辨率圖像.

圖3 SAI插值法Fig.3 Soft-decision adaptive interpolation

因此，對(duì)已配準(zhǔn)的低分辨率齒科序列圖像采用SAI插值法放大factor倍，其計(jì)算公式可表示為：

(18)

2.6 NSCT域圖像融合

插值放大后的低分辨率序列圖像像素增加，但圖像間的信息仍然獨(dú)立，需結(jié)合NSCT變換[12]實(shí)現(xiàn)齒科序列圖像信息的有效融合.相對(duì)于Contourlet變換，NSCT去掉了金字塔分解和方向?yàn)V波器分解中的下采樣及綜合濾波前的插值過程，而改為對(duì)濾波器進(jìn)行上采樣、再對(duì)信號(hào)進(jìn)行分析濾波和綜合濾波，有效減少了信息的丟失.其融合過程可表述為：

(19)

(20)

2)HRi序列圖像經(jīng)過NSCT變換分解后，本文采用的融合策略為：低頻采用加權(quán)均值的融合方法，主要由于超分辨率重建的重點(diǎn)是重建參考幀，所以會(huì)改變參考幀在均值融合過程中所占的比例；高頻采用比較經(jīng)典的取絕對(duì)值最大的融合方法，同時(shí)為有效減少圖像噪聲，會(huì)在高頻融合前，進(jìn)行一次硬閾值去噪，即去噪閾值設(shè)為對(duì)應(yīng)子帶系數(shù)的均方差倍比.

NSCT變換后低頻系數(shù)的融合：

(21)

NSCT變換后高頻系數(shù)的融合：

(22)

(23)

(24)

(25)

3)融合后便可得到優(yōu)化的低頻分量和高頻分量，低分辨率齒科序列圖像間的信息進(jìn)行了整合與互補(bǔ)，采用NSCT反變換，得到最后的超分辨率齒科重建圖像HR.

(26)

其中c-1為NSCT反變換算子.

3 實(shí)驗(yàn)結(jié)果

3.1 建立圖像數(shù)據(jù)庫

本文所采用圖像數(shù)據(jù)庫均來源于通過口腔診療設(shè)備采集的齒科視頻文件，共截取20,000余幀，分辨率皆為640*480.按照實(shí)驗(yàn)的要求，從中選取200幀較為清晰的圖片作為超分辨率重建結(jié)果的參考圖像.根據(jù)這200張圖像，會(huì)根據(jù)圖1的圖像退化模型模擬生成800張低分辨率圖像，即1張高清圖像對(duì)應(yīng)于4張低分辨率圖像.因退化圖像相對(duì)于參考圖像會(huì)存在相對(duì)的幾何運(yùn)動(dòng)，考慮到序列圖像間差異度較小，其平移的像素范圍預(yù)設(shè)為[-10,+10]，角度旋轉(zhuǎn)范圍為內(nèi)[-10°,+10°]；除此之外，生成的退化圖像中還添加了高斯噪聲，其方差范圍為(0.0001～0.0005).最終生成的800張退化圖像構(gòu)成了算法測(cè)試數(shù)據(jù)庫，其中平移情況占30%，旋轉(zhuǎn)情況占30%，旋轉(zhuǎn)與平移兩種情況皆有的例子占40%.

為充分驗(yàn)證本文算法的有效性，我們?cè)O(shè)計(jì)了相關(guān)實(shí)驗(yàn)，所采用的仿真平臺(tái)為：Windows10操作系統(tǒng)、64位Intel Core i5處理器、主機(jī)頻率2.5GHZ、8G運(yùn)行內(nèi)存、Matlab2016b.

3.2 算法評(píng)價(jià)體系及實(shí)驗(yàn)設(shè)定

在超分辨率的重建過程中，配準(zhǔn)精度根據(jù)均方根誤差(root mean square error,RMSE)[25]這一通用方法來評(píng)價(jià),如式(27)，n為配準(zhǔn)點(diǎn)數(shù)，di為配準(zhǔn)后點(diǎn)集與目標(biāo)點(diǎn)集之間歐式距離的測(cè)量值與平均值的偏差:

(27)

對(duì)于超分辨率重建結(jié)果的評(píng)價(jià)方法有主觀的評(píng)價(jià)標(biāo)準(zhǔn)和客觀評(píng)價(jià)的標(biāo)準(zhǔn)，本文中采用峰值信噪比(peak signal-to-noise,PSNR)和結(jié)構(gòu)相似性度量(structural similarity,SSIM)作為最后重建的評(píng)價(jià)指標(biāo)[27].

(28)

(29)

式(28)中，m,n分別表示圖像高和寬的像素?cái)?shù)；P為圖像像素灰度的最大可能取值，對(duì)于一般的灰度圖，P=255.PSNR取值范圍為[0,+)，其值越高，則認(rèn)定圖像質(zhì)量越好.式(29)中，μ1，μ2分別為滑動(dòng)窗口內(nèi)像素的均值；σ1，σ2分別為滑窗內(nèi)像素的方差；σ12為兩幅圖像對(duì)應(yīng)滑窗內(nèi)像素的協(xié)方差；k1，k2為常數(shù)，一般設(shè)定k1=0.01，k2=0.03.SSIM取值區(qū)間為[-1,+1)，其值越大表示圖像質(zhì)量越高.

實(shí)驗(yàn)中，設(shè)定LR低分辨率序列圖像幀數(shù)M為4，重建倍數(shù)factor為4；在配準(zhǔn)點(diǎn)集優(yōu)化的過程中，根據(jù)實(shí)際過程中幾何運(yùn)動(dòng)的程度和圖像大小，閾值T1，T2，T3，T4分別設(shè)定為70、55、70和55；ICP迭代過程中，RMSE閾值設(shè)為0.000001，最大迭代次數(shù)設(shè)為100，魯棒損失函數(shù)選擇Tukey bi-weight函數(shù)；NSCT采用1級(jí)金字塔分解及4次方向分解，分別采用“maxflat”多尺度濾波器和“dmaxflat7”菱形最大平坦濾波器；融合過程中，低頻權(quán)值向量ω設(shè)為(3/6,1/6,1/6,1/6)，高頻硬閾值去噪中的閾值α參考文獻(xiàn)[12]設(shè)為3.

3.3 實(shí)驗(yàn)結(jié)果

1)齒科序列圖像配準(zhǔn)算法驗(yàn)證實(shí)驗(yàn)

SIFT (scale invariant feature transform)[28]，即尺度不變特征變換，是一種基于圖像局部特征的描述子，通過在空間尺度中尋找極值點(diǎn)，并提取出其位置、尺度、旋轉(zhuǎn)不變量.為驗(yàn)證ICP配準(zhǔn)方法的有效性與準(zhǔn)確性，基于測(cè)試樣本庫，將ICP配準(zhǔn)法與SIFT特征點(diǎn)配準(zhǔn)法進(jìn)行實(shí)驗(yàn)比對(duì).

在測(cè)試的200組源圖片中，我們發(fā)現(xiàn)其中存在11幅圖像由于SIFT特征點(diǎn)集的缺少，導(dǎo)致采用SIFT法，不能得出相應(yīng)的單應(yīng)變換矩陣；在其他成功配準(zhǔn)的齒科圖像中，因?yàn)辇X科圖像灰度變化緩慢，分布均勻，造成SIFT特征點(diǎn)提取的不穩(wěn)定，易收斂于局部最小值.表1中結(jié)果說明：SIFT配準(zhǔn)法的RMSE明顯大于ICP算法求出的結(jié)果，而ICP算法穩(wěn)定性較好，計(jì)算復(fù)雜度低.

表1 SIFT配準(zhǔn)與ICP配準(zhǔn)相關(guān)統(tǒng)計(jì)結(jié)果
Table 1 Statistical results related to SIFT and ICP registration

RMSE平均時(shí)間RMSEmaxRMSEminSIFT配準(zhǔn)3.2460.053s15.7130.224ICP配準(zhǔn)0.4340.045s1.4650.012

2)超分辨率重建算法驗(yàn)證實(shí)驗(yàn)

以圖4中低分辨率序列圖像為例，采用本文算法得到的實(shí)驗(yàn)結(jié)果如下所示：圖4(a)是從低分辨率序列圖像中選取的參考幀；圖4(b)是在低分辨率序列中，與參考幀存在幾何相對(duì)運(yùn)動(dòng)的第i幀；經(jīng)過ICP配準(zhǔn)算法的迭代計(jì)算后，可求出位移變換矩陣，圖4(c)給出了對(duì)序列圖像中第i幀配準(zhǔn)后的結(jié)果.

在ICP配準(zhǔn)的過程中，特征點(diǎn)集的提取與變換過程如圖4(d)至圖4(f)所示：圖4(d)是未優(yōu)化的配準(zhǔn)點(diǎn)集，存在很多干擾點(diǎn)；圖4(e)中點(diǎn)集經(jīng)過優(yōu)化后，ICP算法的包含關(guān)系得到滿足；圖4(f)顯示了點(diǎn)集的配準(zhǔn)結(jié)果.結(jié)合SAI插值法對(duì)圖4(f)進(jìn)行處理，得到高分辨率圖4(g).對(duì)序列中的其他幀采用同樣的方法進(jìn)行配準(zhǔn)和插值處理后，得到高分辨率的序列圖像，對(duì)該序列圖像中的每一幀依次進(jìn)行NSCT變換，并采用基于NSCT變換域的改進(jìn)融合策略，得到重建后的低頻分量和高頻分量如圖4(h)和圖4(i)所示.圖4(j)對(duì)應(yīng)于NSCT分解后的頻譜圖；圖4(k)為重建結(jié)果；圖4(l)為真值圖.在該例中，ICP配準(zhǔn)結(jié)果的均方根誤差RMSE為0.4634，重建結(jié)果的PSNR和SSIM分別為37.029和0.909.

圖4 典型案例重建算法示意圖Fig.4 Reconstruction steps for a typical case

圖5 3組實(shí)例對(duì)應(yīng)低分辨率序列與Ground Truth圖Fig.5 Corresponding LR sequence and Ground Truth diagram

3)超分辨率重建算法比對(duì)實(shí)驗(yàn)

為了充分顯示本文算法的有效性與優(yōu)越性，現(xiàn)結(jié)合圖5實(shí)例分析對(duì)比4種算法重建結(jié)果：1)直接采用cubic插值法；2)對(duì)序列低分辨率圖像采用本文中的改進(jìn)ICP算法進(jìn)行幀間配準(zhǔn)，配準(zhǔn)后的序列圖像采用SAI插值算法，然后在Contourlet變換域(contourlet transform,CT)中，采用與本文算法相同的融合策略完成序列圖像的超分辨率重建；3)對(duì)序列低分辨率圖像采用SIFT算法進(jìn)行幀間配準(zhǔn)，配準(zhǔn)后的序列圖像采用SAI插值算法，然后在NSCT變換域中，采用與本文算法相同的融合策略完成序列圖像的超分辨率重建；4)本文所提出的超分辨率重建算法.

實(shí)驗(yàn)結(jié)果如表2所示，表2中每一例的第一行對(duì)應(yīng)PSNR值、第二行對(duì)應(yīng)SSIM值.

表2 不同算法實(shí)驗(yàn)結(jié)果
Table 2 Evaluation results of different algorithms

CubicICP+SAI+CTSIFT+SAI+NSCTICP+SAI+NSCT例132.72335.98820.81536.0880.8110.8860.7590.890例234.09237.24319.15737.9830.8380.9110.8210.930例334.22838.42838.48138.8090.8200.8840.8920.901

由實(shí)驗(yàn)數(shù)據(jù)可知，3例LR圖像經(jīng)cubic插值，重建后的HR圖像PSNR和SSIM都較低；而采用方法2后，通過ICP配準(zhǔn)和在Contourlet變換域中完成序列圖像的融合，評(píng)價(jià)指標(biāo)明顯上升；方法4即本文算法，將Contourlet變換改成NSCT后，PSNR和SSIM得到了進(jìn)一步提高，圖像效果更佳；如果將方法4中的改進(jìn)ICP配準(zhǔn)算法變更為SIFT特征點(diǎn)配準(zhǔn)算法，即方法3，當(dāng)配準(zhǔn)精度與ICP配準(zhǔn)相當(dāng)時(shí)，如例3重建結(jié)果表現(xiàn)良好，例1和例2則由于配準(zhǔn)誤差較大，使得重建的指標(biāo)急劇下降，PSNR只有20.815和19.157.

通過實(shí)驗(yàn)比對(duì)分析，本文所提出的算法取得了最好的實(shí)驗(yàn)效果.其中，改進(jìn)ICP算法和基于NSCT變換域的融合策略較其它算法，在齒科序列圖像的超分辨率重建過程中，具有一定的優(yōu)越性.

4)測(cè)試圖庫統(tǒng)計(jì)分析結(jié)果

為定量分析本文算法的有效性，使用超分辨率重建算法比對(duì)實(shí)驗(yàn)中所提及的4種不同算法，對(duì)200張由源圖像生成的800張測(cè)試圖像進(jìn)行了超分辨率重建，并統(tǒng)計(jì)了重建后的評(píng)價(jià)指標(biāo).齒科序列圖像重建結(jié)果的峰值信噪比PSNR如圖6所示：其中橫坐標(biāo)是PSNR峰值信噪比；縱坐標(biāo)是對(duì)應(yīng)PSNR值下的樣本個(gè)數(shù)；白灰條形圖表示cubic雙三次插值重建樣本數(shù)(方法1)；淺灰條形圖表示使用ICP配準(zhǔn)、SAI插值、Contourlet變換域融合重建的樣本數(shù)(方法2)；深灰條形圖為SIFT特征點(diǎn)配準(zhǔn)、SAI插值然后采用NSCT分解融合的重建樣本數(shù)(方法3)；黑色條形圖表示的是本文提出算法的重建樣本數(shù)(方法4).

由圖6可以看出，cubic插值法重建結(jié)果(方法1)相對(duì)于方法2和方法4的重建結(jié)果，其PSNR值分布在較低數(shù)值段；而比對(duì)方法2和方法4，即基于Contourlet變換域融合與基于NSCT變換域融合的實(shí)驗(yàn)對(duì)比中，即使采用相同的融合策略，在較高PSNR值區(qū)間內(nèi)，方法4所對(duì)應(yīng)的重建樣本數(shù)明顯多于方法2所對(duì)應(yīng)的重建樣本數(shù)；另一方面，比對(duì)方法3和方法4，即都在NSCT變換域中進(jìn)行融合時(shí)，雖然有些通過SIFT特征點(diǎn)配準(zhǔn)算法所得的重建結(jié)果與通過ICP配準(zhǔn)算法所得重建結(jié)果在PSNR值上相當(dāng)，但是由于SIFT配準(zhǔn)在齒科圖像應(yīng)用上效果的不穩(wěn)定性，使得方法3所對(duì)應(yīng)的PSNR值在圖6上呈現(xiàn)出零散分布的狀態(tài).因此，綜合穩(wěn)定性與峰值信噪比來看，本文算法有效提高了齒科圖像的重建效果.

圖6 四種不同超分辨率重建算法重建結(jié)果PSNR分布Fig.6 PSNR statistical distribution of 4 different SR algorithms

基于構(gòu)建的齒科低分辨率圖像庫，表3分別列出了4種不同重建算法重建結(jié)果的平均指標(biāo)，包括峰值信噪比PSNR、均方根誤差MSE和結(jié)構(gòu)相似性度量SSIM.可以看出，cubic插值算重建效果最差，PSNR為33.275，SSIM為0.824；而當(dāng)采用方法2，即采用ICP配準(zhǔn)、SAI插值和Contourlet變換融合時(shí)，相應(yīng)指標(biāo)明顯上升，低分辨率序列圖像信息得到了有效利用；但是相對(duì)于NSCT，CT變換過程中由于存在下采樣的過程，不可避免會(huì)損失圖像信息，細(xì)節(jié)等，NSCT就表現(xiàn)出更優(yōu)的效果.

表3 不同超分辨率重建方法的平均精度對(duì)比
Table 3 Average indicators comparison of different SR reconstruction methods

PSNRMSESSIMCubic33.27531.0130.824ICP+SAI+CT37.28512.8120.907SIFT+SAI+NSCT35.72591.6100.882ICP+SAI+NSCT37.53912.1060.917

另外，重建結(jié)果的好壞還嚴(yán)重依賴于配準(zhǔn)的精度，SIFT在齒科圖像重建過程中，由于特征點(diǎn)的缺少和誤匹配的存在，導(dǎo)致了方法3整體上劣于方法4，即本文提出算法.本文所提出針對(duì)齒科低分辨率序列圖像重建的方法，無論是在精度上還是穩(wěn)定性上都有較大的優(yōu)勢(shì).

3.4 實(shí)驗(yàn)結(jié)果分析與討論

通過對(duì)實(shí)驗(yàn)結(jié)果的分析，本文所提出的超分辨率齒科序列圖像重建算法，能夠有效提高重建高分辨率圖像的表現(xiàn)及客觀指標(biāo)，其主要原因在于：

1)經(jīng)典ICP配準(zhǔn)法要求目標(biāo)點(diǎn)集與源點(diǎn)集存在包含關(guān)系，本文則提出了改進(jìn)ICP的配準(zhǔn)策略，引入了魯棒損失函數(shù)，又優(yōu)化了點(diǎn)集的篩選，使得低分辨率序列圖像的配準(zhǔn)精度得到了有效保證，同時(shí)跟SIFT特征點(diǎn)配準(zhǔn)法相比，明顯魯棒性較強(qiáng)；

2)軟判決自適應(yīng)插值法的應(yīng)用，有效保護(hù)了圖像的空間一致性，邊緣以及牙齒圖像細(xì)節(jié)等信息；

3)NSCT變換相對(duì)于傅里葉變換、小波變換等，利用類似輪廓段的結(jié)構(gòu)表征圖像，具有良好的稀疏特性，同時(shí)由于沒有下采樣過程，圖像保留細(xì)節(jié)的能力加強(qiáng)；

4)對(duì)NSCT變換后的系數(shù)進(jìn)行了有效融合，低分辨率序列圖像重建主要重建參考幀，低頻融合時(shí)加強(qiáng)了參考幀所占比重；高頻融合時(shí)，硬閾值的去噪可以降低噪聲，然后又采用取絕對(duì)值最大的方法融合，實(shí)現(xiàn)信息的有效互補(bǔ)，整體提高超分辨率的重建效果.

但是從測(cè)試樣本庫的重建結(jié)果來看，基于本文所提出的超分辨率齒科序列圖像重建算法仍存在不足之處，主要如下：

1)基于點(diǎn)集配對(duì)的ICP配準(zhǔn)法，仍然無法做到完全的正確配對(duì)，同時(shí)干擾點(diǎn)的存在會(huì)影響配準(zhǔn)結(jié)果；

2)多分辨率多尺度的變換不可避免的會(huì)損失一些圖像信息，另外由于既要去除噪聲，又要充分保留圖像信息，而且序列間的互補(bǔ)信息需要有效融合，之間的權(quán)衡關(guān)系無法做到完全精準(zhǔn)，因此本文所提出的融合策略需要進(jìn)一步加強(qiáng).

4 結(jié)束語

本文提出了一種基于齒科序列圖像的超分辨率重建算法，該算法首先對(duì)低分辨率的齒科序列圖像進(jìn)行預(yù)處理，減少光照的影響；然后選取牙體硬組織的邊緣點(diǎn)作為特征點(diǎn)集，通過改進(jìn)的ICP配準(zhǔn)法，在點(diǎn)集篩選和迭代策略上進(jìn)行創(chuàng)新，引入魯棒損失函數(shù)，以此矯正低分辨率齒科序列幀間的幾何運(yùn)動(dòng)，有效減少了偏差較大點(diǎn)的干擾；其次，已配準(zhǔn)的序列圖像選取有邊緣保護(hù)作用的SAI插值法放大后，基于超分辨率重建參考幀的特點(diǎn)，通過調(diào)整NSCT變換域的去噪與融合策略，實(shí)現(xiàn)信息的有效互補(bǔ)，最后利用NSCT的反變換得到超分辨率齒科圖像的重建結(jié)果.實(shí)驗(yàn)結(jié)果表明，不論是主觀視覺效果上、還是峰值信噪比、結(jié)構(gòu)相似性等客觀指標(biāo)上，本文算法都能夠有效提高齒科序列圖像的超分辨率重建效果.