， ，

(浙江工業(yè)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，浙江 杭州 310014)

圖像配準(zhǔn)是圖像處理研究領(lǐng)域中的一個(gè)典型問題和技術(shù)難點(diǎn)，其目的在于對(duì)同一對(duì)象在不同條件下(不同的采集設(shè)備、不同的拍攝視角和不同的時(shí)間)獲取的圖像進(jìn)行對(duì)比、融合和拼接，圖像配準(zhǔn)是圖像融合的基礎(chǔ)，圖像配準(zhǔn)和圖像融合又是圖像拼接的關(guān)鍵技術(shù).例如，車牌識(shí)別，視頻全景拼接[1]，等等.具體地說，對(duì)于一組圖像數(shù)據(jù)集中的兩幅圖像，通過尋找一種空間變換把一幅圖像映射到另一幅圖像，使得兩幅圖像中對(duì)應(yīng)于空間同一位置的點(diǎn)一一對(duì)應(yīng)起來，從而達(dá)到信息融合的目的.其在醫(yī)學(xué)圖像處理、計(jì)算機(jī)視覺、計(jì)算機(jī)圖形學(xué)、運(yùn)動(dòng)分析和目標(biāo)檢測(cè)中都有著廣泛的應(yīng)用，Szeliski對(duì)其作了一定的回顧[2].本研究關(guān)注的是其在醫(yī)學(xué)圖像處理中的應(yīng)用.為了提高醫(yī)學(xué)診斷和治療水平，醫(yī)生需要對(duì)不同的醫(yī)學(xué)圖像進(jìn)行定性和定量分析，從而獲得更全面的信息，其中就會(huì)用到圖像配準(zhǔn)技術(shù)，圖像質(zhì)量也將影響到醫(yī)療診斷的準(zhǔn)確性[3].在醫(yī)學(xué)領(lǐng)域，諸多問題都會(huì)涉及到頸動(dòng)脈血管，其中一個(gè)典型的問題就是頸動(dòng)脈血管位移場(chǎng)的計(jì)算.在一個(gè)心動(dòng)周期內(nèi)，頸動(dòng)脈血管會(huì)隨著心臟的收縮和舒張發(fā)生細(xì)微的形變，由于形變量很小，其位移場(chǎng)很難被準(zhǔn)確計(jì)算出來.然而，通過配準(zhǔn)可以實(shí)現(xiàn)位移場(chǎng)的計(jì)算，但是已有的配準(zhǔn)方法鮮有針對(duì)頸動(dòng)脈血管的，其精度難以達(dá)到理想的效果，這就需要尋求新的方法來滿足要求.通過將頸動(dòng)脈血管圖像的形狀特征和紋理特征相結(jié)合，并進(jìn)行迭代運(yùn)算，可以將配準(zhǔn)精度提升到一個(gè)更高的層次，從而完成位移場(chǎng)的計(jì)算工作.這在醫(yī)學(xué)診療上具有很重要的現(xiàn)實(shí)意義，例如，動(dòng)脈粥樣硬化是很多心腦血管疾病中的主要誘因，其形成過程較為復(fù)雜，隨著年齡的增加或其他誘因，血管壁會(huì)逐漸增厚，這是動(dòng)脈粥樣硬化的早期表現(xiàn)，進(jìn)一步發(fā)展就有可能形成斑塊，斑塊嚴(yán)重到一定程度，就會(huì)導(dǎo)致血管腔狹窄，斑塊破裂也會(huì)導(dǎo)致血栓的形成.由此可見，本研究提出的方法具有較高的醫(yī)學(xué)價(jià)值.

1 相關(guān)工作

在圖像配準(zhǔn)過程中，主要涉及以下幾個(gè)方面：首先是特征空間，即圖像的特征提?。黄浯问撬阉骺臻g，即圖像變換的范圍和方式；第三是相似性度量；最后是搜索策略，即目標(biāo)函數(shù)的優(yōu)化.

在最初階段，研究人員們采用灰度不變性這一假設(shè)來建立圖像間的對(duì)應(yīng)關(guān)系，但是他們很快意識(shí)到，由于光照、視角和噪聲等因素的存在，圖像的匹配結(jié)果并不可靠.后來，互信息[4]、梯度[5]被用來進(jìn)行圖像匹配，但其在面對(duì)圖像劇烈變化時(shí)的表現(xiàn)并不是很好.之后，SIFT[6],shape context[7],HOG[8]等被應(yīng)用到圖像匹配當(dāng)中，并廣為流傳.

進(jìn)一步地，圖像間對(duì)應(yīng)關(guān)系可以分為稀疏表示和稠密表示.僅僅利用圖像的特征點(diǎn)，可以獲得一幅圖像的稀疏表示，通常情況下，基于稀疏表示的算法是比較有效的.本研究也可以利用圖像中每一個(gè)像素點(diǎn)的信息獲得一個(gè)圖像的稠密表示，利用圖像間的稠密對(duì)應(yīng)關(guān)系，可以很容易地將一幅圖像變形到另一幅圖像.

在圖像配準(zhǔn)領(lǐng)域，研究人員們提出了多種多樣的方法.Belongie等[7]提出了shape context描述符，其為利用目標(biāo)的形狀輪廓特征進(jìn)行配準(zhǔn)提供了可能；Liu等[9]提出了SIFT Flow，可以使用目標(biāo)圖像逐像素的SIFT特征進(jìn)行配準(zhǔn)；Ying等[10]提出了一種利用圖像的分布指導(dǎo)配準(zhǔn)的方法；Ye等[11]提出了一種基于流形的可變形的配準(zhǔn)方法；Seiler等[12]利用多尺度的局部仿射變換進(jìn)行配準(zhǔn)，其在仿射變換中引入先驗(yàn)知識(shí)，避免了局部極小值，并用生成的統(tǒng)計(jì)模型來重新表示仿射變換，通過最優(yōu)化后驗(yàn)概率的最大值來找到最佳變換；Wachinger等[13]提出了一種解決多幅圖像同時(shí)配準(zhǔn)問題的方法；Cordero-Grande等[14]提出了一個(gè)針對(duì)時(shí)間序列圖像聯(lián)合對(duì)齊的方法理論，該理論基于一個(gè)群組配準(zhǔn)過程，該過程使用預(yù)期的稀疏的時(shí)間強(qiáng)度曲線來匹配被對(duì)齊的點(diǎn)，該方法已經(jīng)被應(yīng)用在了強(qiáng)對(duì)比度的初次灌注的心臟磁共振圖像的運(yùn)動(dòng)校正中.本研究旨在尋求一種適合于頸動(dòng)脈血管的高精度的配準(zhǔn)方法.

2 基于形狀特征和紋理特征的配準(zhǔn)方法

對(duì)于頸動(dòng)脈血管，其內(nèi)壁和外壁具有可見的輪廓，因此可以提取其形狀特征，利用形狀上下文信息進(jìn)行匹配和配準(zhǔn).同時(shí)，由于組織結(jié)構(gòu)的不同，其成像后具有區(qū)別于周圍組織的紋理特征，我們可以基于逐像素的SIFT特征對(duì)其進(jìn)行配準(zhǔn).

首先把一個(gè)目標(biāo)的形狀用其輪廓上的一組離散的點(diǎn)集來表示.對(duì)于一個(gè)含有n個(gè)點(diǎn)的集合P，當(dāng)n越大時(shí)，可以獲得越好的形狀估計(jì)，即

P={p1,p2,…,pn}pi∈IR2

(1)

對(duì)于輪廓上的任意一點(diǎn)pi，其形狀上下文信息被定義為一個(gè)點(diǎn)其關(guān)于其余n-1個(gè)點(diǎn)相關(guān)位置的直方圖hi.通過在極坐標(biāo)中表示其相互間的位置關(guān)系可以使其對(duì)距離關(guān)系更敏感，即

hi(k)=#{q≠pi:(q-pi)∈bin(k)}

(2)

對(duì)于形狀1的任意一點(diǎn)pi和形狀2中的任意一點(diǎn)qj，定義這兩個(gè)點(diǎn)的匹配代價(jià)Cij，即

(3)

式中：hi(k)，hj(k)分別為點(diǎn)pi和qj的K-bin歸一化直方圖.進(jìn)而，對(duì)于所有的點(diǎn)對(duì)pi和qj，我們要最小化其總的匹配代價(jià)H(π)，即

(4)

式中π為一個(gè)排列.

Shape context作為一個(gè)描述符本身具有平移不變性，以n2對(duì)點(diǎn)的平均徑向距離來歸一化這n2對(duì)點(diǎn)可以使其具有尺度不變性，把每個(gè)點(diǎn)處的切向量看作是與x軸正向的切角可以使描述符具有旋轉(zhuǎn)不變性.

SIFT作為一個(gè)局部描述符，可以描述局部梯度信息[6].圖像中每一個(gè)像素的4×4×8=128維向量作為一個(gè)像素的SIFT描述，進(jìn)而可以得到一幅逐像素SIFT描述符表示的圖像.

利用兩幅圖像中逐像素的描述符，將這些描述符進(jìn)行匹配并建立一個(gè)稠密對(duì)應(yīng)關(guān)系，進(jìn)而進(jìn)行坐標(biāo)變換，最后完成配準(zhǔn).

參數(shù)模型的選擇是我們下面要解決的問題.首先需要估計(jì)一個(gè)平面變換T：IR2→IR，將一個(gè)點(diǎn)集映射到另一個(gè)點(diǎn)集.

在此選擇薄板樣條函數(shù)模型(TPS)[15-16].Bookstein[17]發(fā)現(xiàn)其在應(yīng)對(duì)生物形態(tài)變化的建模上很有效.構(gòu)建坐標(biāo)變換模型T(x,y)，即

T(x,y)=(fx(x,y),fy(x,y))

(5)

其中f(x,y)表示薄板樣條插值函數(shù).

在得到了初步配準(zhǔn)的圖像之后，進(jìn)一步用逐像素SIFT信息進(jìn)行配準(zhǔn)以達(dá)到更精確的結(jié)果.其能量函數(shù)定義為

(α|u(p)-u(q)|,d)+min(α|v(p)-

v(q)|,d)

(6)

式中：p=(x,y)為圖像中p點(diǎn)的坐標(biāo)；w(p)=(u(p)，v(p))為p點(diǎn)的流向量，u(p)為水平流，v(p)為垂直流；s1和s2分別為兩幅待匹配的SIFT圖像；ε為空間鄰域；α,t,d分別為閾值.能量函數(shù)中包含了數(shù)據(jù)項(xiàng)，位移項(xiàng)和正則項(xiàng).其中數(shù)據(jù)項(xiàng)用來約束SIFT描述符沿著流向量w(p)的方向進(jìn)行匹配，位移項(xiàng)用來約束流向量盡可能的小，正則項(xiàng)用來約束鄰接像素的流向量盡可能相似.

在優(yōu)化目標(biāo)函數(shù)時(shí)，使用了雙層環(huán)路置信傳播.因?yàn)槟繕?biāo)函數(shù)中含有L1范數(shù)，可以使用距離變換函數(shù)[18]進(jìn)一步降低復(fù)雜度，連續(xù)的置信傳播(BP-S)[19]也可以使其更好地收斂.

最后將上述兩部分進(jìn)行迭代，通過手動(dòng)設(shè)置迭代次數(shù)，使得配準(zhǔn)結(jié)果達(dá)到最佳，從而進(jìn)一步優(yōu)化我們的配準(zhǔn)結(jié)果，即

(7)

式中：Rij為參考圖像對(duì)應(yīng)位置像素的灰度值；Fij為浮動(dòng)圖像對(duì)應(yīng)位置像素的灰度值；M·N為圖像的大小.

通過式(7)計(jì)算參考圖像和配準(zhǔn)后圖像的均方誤差(Mean square error，MSE)，并以此來評(píng)估配準(zhǔn)結(jié)果.本研究方法算法流程如圖1所示.

圖1 算法流程圖Fig.1 The flowchart of algorithm

在本方法中，我們先采用形狀上下文特征進(jìn)行了初步的配準(zhǔn)，之后采用逐像素SIFT特征進(jìn)行精確的配準(zhǔn).選用這樣的順序是考慮到原始圖像具有更為自然的形狀特征，先采用形狀上下文特征可以更充分地發(fā)揮目標(biāo)形狀特征在初步配準(zhǔn)的效力，而這一步處理對(duì)后續(xù)利用逐像素SIFT特征進(jìn)行精確配準(zhǔn)影響較小.對(duì)于先采用形狀上下文特征還是逐像素SIFT特征，通過實(shí)驗(yàn)對(duì)比發(fā)現(xiàn)先采用形狀上下文特征要優(yōu)于先采用逐像素SIFT特征.

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

在實(shí)驗(yàn)中，選取了心動(dòng)周期內(nèi)一個(gè)序列的人體頸部橫切面圖像(圖2)，該序列圖像包含16幅大小為512×512的圖像.通過截取等處理，獲得頸動(dòng)脈血管的圖像(圖3).

圖2 人體頸部橫切面圖像Fig.2 The image of human neck transection

圖3 頸動(dòng)脈橫切面圖像Fig.3 The image of carotid artery transection

在對(duì)頸動(dòng)脈血管切面進(jìn)行配準(zhǔn)的過程當(dāng)中，選取序列中第一幅圖像作為參考圖像，其余圖像作為浮動(dòng)圖像進(jìn)行配準(zhǔn).首先要得到的就是血管內(nèi)外壁的輪廓，圖像分割技術(shù)目前已經(jīng)實(shí)現(xiàn)血管內(nèi)腔的分割和建模[20-21].對(duì)于血管內(nèi)壁輪廓來說，由于其與血管腔具有明顯的灰度差異，正如圖2中看到的那樣，是可以進(jìn)行邊緣檢測(cè)的.然而，血管外壁盡管具有可見的輪廓，一般的邊緣檢測(cè)方法卻無能為力，因此選擇手動(dòng)標(biāo)出血管內(nèi)外壁輪廓，進(jìn)而在其上進(jìn)行點(diǎn)采樣，從而獲得一個(gè)邊緣點(diǎn)集.為了更好地表示輪廓以及點(diǎn)與點(diǎn)之間的位置關(guān)系，應(yīng)盡可能均勻的在輪廓線上進(jìn)行點(diǎn)采樣.最終在頸動(dòng)脈血管內(nèi)外壁邊緣上選取了400個(gè)點(diǎn)(圖4)，圖4(a)為血管內(nèi)外壁采樣點(diǎn)在配準(zhǔn)前的位置和匹配關(guān)系，圖4(b)為血管內(nèi)外壁采樣點(diǎn)在配準(zhǔn)后的位置關(guān)系.在比較配準(zhǔn)結(jié)果時(shí)，為了更具說明性，選擇一個(gè)圓形模板突出血管部分，并計(jì)算其均方誤差來評(píng)估配準(zhǔn)結(jié)果.本研究提出的方法對(duì)頸動(dòng)脈血管的配準(zhǔn)結(jié)果如圖5(a)所示，圖5(b)為使用sift flow進(jìn)行配準(zhǔn)后的結(jié)果，圖5(c)為使用shape context進(jìn)行配準(zhǔn)后的結(jié)果.

圖4 血管壁采樣點(diǎn)配準(zhǔn)前后圖像Fig.4 The image of carotid artery wall samples before and after registration

圖5 不同方法配準(zhǔn)后圖像Fig.5 The registered image of different method

表1比較了不同方法配準(zhǔn)后的均方誤差(MSE)，MSE(mRn)中m表示參考圖像序號(hào)，n表示浮動(dòng)圖像序號(hào)，隨著迭代次數(shù)的增加，整體看來，配準(zhǔn)精度在增加，配準(zhǔn)時(shí)間也隨之增加，而在3～4次迭代后，配準(zhǔn)結(jié)果趨于穩(wěn)定并達(dá)到最優(yōu).從圖6中可以直觀地看出不同配準(zhǔn)方法之間的比較，其中“·”表示配準(zhǔn)之前兩幅圖像的均方誤差；“*”表示利用SIFT Flow進(jìn)行配準(zhǔn)之后的均方誤差；“×”表示利用shape context進(jìn)行配準(zhǔn)之后的均方誤差；“■”表示應(yīng)用本研究所提出的配準(zhǔn)方法配準(zhǔn)之后的均方誤差.在單獨(dú)使用SIFT Flow或shape context進(jìn)行配準(zhǔn)時(shí)，其表現(xiàn)并不穩(wěn)定，會(huì)出現(xiàn)無法配準(zhǔn)的情況，而本研究所提出的方法在頸動(dòng)脈血管配準(zhǔn)的應(yīng)用中，相較于單一的配準(zhǔn)方法，具有更高的配準(zhǔn)精度.

表1 采用不同方法配準(zhǔn)后圖像的均方誤差Table 1 The comparison of mean squared error after registration using different method

圖6 不同方法配準(zhǔn)后的MSEFig.6 The MSE of different method after registration

4 結(jié) 論

基于目標(biāo)形狀特征和紋理特征的迭代配準(zhǔn)方法，為解決人體頸動(dòng)脈血管的配準(zhǔn)問題提供了一個(gè)有效途徑，對(duì)于頸動(dòng)脈血管位移場(chǎng)的計(jì)算也提供了重要的參考依據(jù)，這些對(duì)于醫(yī)學(xué)診斷和醫(yī)學(xué)治療都具有非常重要的意義.通過利用圖像的不同特征進(jìn)行迭代配準(zhǔn)，筆者方法達(dá)到了更精確的配準(zhǔn)結(jié)果，具有較好的魯棒性，在運(yùn)行時(shí)間上也具有了可接受的水平.可以預(yù)見，對(duì)于處理形狀特征和紋理特征明顯的圖像的配準(zhǔn)問題，這也將是一個(gè)不錯(cuò)的選擇.

致謝：劍橋大學(xué)放射系騰忠照老師研究團(tuán)隊(duì)為本研究提供了實(shí)驗(yàn)數(shù)據(jù)，在此表示衷心的感謝.

[1] 范菁，吳佳敏，葉陽，等.帶顯著性區(qū)域約束的高效視頻全景拼接方法[J].浙江工業(yè)大學(xué)學(xué)報(bào)，2015，43(5)：479-486.

[2] SZELISKI R. Image alignment and stiching:a tutorial[J]. Foundations and trends in computer graphics and computer vision,2006,2(1):1-104.

[3] 張劍華，張自然，汪曉妍，等.基于結(jié)構(gòu)顯著性的醫(yī)學(xué)圖像質(zhì)量評(píng)價(jià)[J].浙江工業(yè)大學(xué)學(xué)報(bào)，2015，43(6)：636-641.

[4] VIOLA P, WELLS W M I. Alignment by maximization of mutual information[J]. IEEE computer society,1995,24(2):16-23.

[5] BROX T, BRUHN A, PAPENBERG N, et al. High accuracy optical flow estimation based on a theory for warping[J]. Lecture notes in computer science,2004,3024(10):25-36.

[6] LOWE D G. Object recognition from local scale-invariant features[C]//ICCV. IEEE International Conference on Computer Vision. Kerkyra, Greece: IEEE,1999:1150-1157.

[7] BELONGIE S, MALIK J, PUZICHA J. Shape matching and object recognition using shape contexts[J]. IEEE transactions

on pattern analysis and machine intelligence,2002,24(4):509-522.

[8] DALAL N, TRIGGS B. Histograms of oriented gradients for human detection[C]//CVPR. IEEE Conference on Computer Vision and Pattern Recognition. San Diego, CA, USA: IEEE,2005:886-893.

[9] LIU C, YUEN J, TORRALBA A. SIFT flow: dense correspondence across scenes and its applications[J]. IEEE transcations on pattern analysis and machine intelligence,2011,33(5):978-994.

[10] YING S, WU G, WANG Q, et al. Groupwise registration via graph shrinkage on the image manifold[C]//CVPR. IEEE Conference on Computer Vision and Pattern Recognition. Oregon, Portland: IEEE,2013:2323-2330.

[11] YE D H, HAMM J, KWON D, et al. Regional manifold learning for deformable registration of brain MR images[C]//MICCAI. International Conference on Medical Image Computing and Computer-Assisted Intervention. Nice, France: Springer Berlin Heidelberg,2012:131-138.

[12] SEILER C, PENNEC X, REYES M. Simultaneous multiscale polyaffine registration by incorporating deformation statistics[C]//MICCAI. International Conference on Medical Image Computing and Computer-Assisted Intervention. Nice, France: Springer Berlin Heidelberg,2012:130-137.

[13] WACHINGER C, NAVAB N. Simultaneous registration of multiple images: similarity metrics and efficient optimization[J]. IEEE transactions on pattern analysis and machine intelligence,2013,35(5):1221-1233.

[14] CORDERO-GRANDE L, MERINO-CAVIEDES S, AJA-FERNNDEZ S, et al. Groupwise elastic registration by a new sparsity-promoting metric: application to the alignment of cardiac magnetic resonance perfusion images[J]. IEEE transactions on pattern analysis and machine intelligence,2013,35(11):2638-2650.

[15] DUCHON J. Splines minimizing rotation-invariant semi-norms in Sobolev spaces[J]. Lecture notes in mathematics,2006(571):85-100.

[16] MEINGUET J. Multivariate interpolation at arbitrary points made simple[J]. Zeitschrift für angewandte mathematik und physik ZAMP,1979,30(2):292-304.

[17] BOOKSTEIN F L. Morphometric tools for landmark data: geometry and biology[M]. Cambridge: Cambridge University Press,1997.

[18] FELZENSZWALB P F, HUTTENLOCHER D P. Efficient belief propagation for early vision[J]. International journal of computer vision,2006,1(1):41-54.

[19] SZELISKI R, ZABIH R, SCHARSTEIN D, et al. A comparative study of energy minimization methods for markov random fields with smoothness-based priors[J]. IEEE transactions on pattern analysis and machine intelligence,2008,30(6):1068-1080.

[20] 鄭河榮，劉家好，何玲娜.腦血管計(jì)算機(jī)模型的建立與有限元分析[J].浙江工業(yè)大學(xué)學(xué)報(bào)，2014，42(3)：253-256.

[21] 尚巖峰，汪輝，汪寧，等.管狀特性和主動(dòng)輪廓的三位血管自動(dòng)提取[J].中國圖像圖形學(xué)報(bào)，2013，18(3)：290-298.