非剛性醫(yī)學(xué)圖像配準(zhǔn)算法的設(shè)計(jì)與實(shí)現(xiàn)

1、非剛性醫(yī)學(xué)圖像配準(zhǔn)算法的設(shè)計(jì)與實(shí)現(xiàn)【關(guān)鍵詞】醫(yī)學(xué)圖像；非剛性；圖像配準(zhǔn)；匹配矩陣；薄板樣條摘要：非剛性圖像匹配問題已成為醫(yī)學(xué)圖像分析中一個(gè)非常具有挑戰(zhàn)性的問題?；诒“鍢訔l插值方法,引入實(shí)匹配矩陣,并給出相應(yīng)配準(zhǔn)變換算法，該算法將薄板樣條參數(shù)表示成仿射分量和非仿射分量，并分別進(jìn)展求解。與其它非剛性匹配算法相比，該算法不僅保證了對(duì)應(yīng)特征點(diǎn)的雙向?qū)?yīng)，也實(shí)現(xiàn)了自動(dòng)特征點(diǎn)選擇，實(shí)驗(yàn)結(jié)果令人滿意。關(guān)鍵詞：醫(yī)學(xué)圖像；非剛性；圖像配準(zhǔn)；匹配矩陣；薄板樣條1引言在醫(yī)學(xué)診斷和治療過程中，常需要比照分析多幅圖像，以獲得更為準(zhǔn)確和全面的信息。圖像分析大都要求多幅圖像的幾何位置一致，因此，配準(zhǔn)是醫(yī)學(xué)圖像分析的一個(gè)重

2、大課題。醫(yī)學(xué)圖像配準(zhǔn)是指對(duì)于一幅醫(yī)學(xué)圖像尋求一種(或一系列)空間變換，使它與另一幅醫(yī)學(xué)圖像上的對(duì)應(yīng)點(diǎn)到達(dá)空間上的一致。這種一致是指人體上的同一解剖點(diǎn)在兩張匹配圖像上有一樣的空間位置。配準(zhǔn)的結(jié)果應(yīng)使兩幅圖像上所有的解剖點(diǎn)，或至少是所有具有診斷意義的點(diǎn)及手術(shù)感興趣的點(diǎn)都到達(dá)匹配。圖像配準(zhǔn)不僅可以校正病人屢次成像間的位置變化，也可以校正由于成像形式本身導(dǎo)致的畸變。對(duì)同一個(gè)病人的不同時(shí)間的圖像進(jìn)展配準(zhǔn)，可以理解發(fā)育過程及腫瘤病變的病情；對(duì)不同人的圖像進(jìn)展配準(zhǔn)，去除種族、年齡等臨床及遺傳差異，從而形成疾病或人群特異性圖譜，可用于正常與否的分析；對(duì)不同成像形式進(jìn)展配準(zhǔn)，可以獲得互補(bǔ)信息。醫(yī)學(xué)圖像配準(zhǔn)可分為

3、剛性配準(zhǔn)和非剛性配準(zhǔn)兩類。剛性配準(zhǔn)在許多情況下不能滿足臨床的需要，因?yàn)楹芏嘈巫兊男再|(zhì)是非剛體、非線性的。比方為了準(zhǔn)確定位R圖像左心室，常常伴有組織磁化系數(shù)差異、非水分子的化學(xué)位移以及血流流動(dòng)等因素導(dǎo)致的幾何畸變以及由于磁場(chǎng)不均勻、磁場(chǎng)梯度非線性及渦流等導(dǎo)致的探測(cè)畸變，因此在放療方案制定中，將R圖像配準(zhǔn)時(shí)，不能單純地使用剛性配準(zhǔn)，必須使用非剛性配準(zhǔn)。非剛性配準(zhǔn)算法可分為灰度驅(qū)動(dòng)、模型驅(qū)動(dòng)及混合算法三種13。灰度驅(qū)動(dòng)方法基于數(shù)學(xué)或統(tǒng)計(jì)尺度將一個(gè)灰度形式與另一個(gè)對(duì)準(zhǔn)。典型情況下，需要定義源系統(tǒng)與目的系統(tǒng)之間的灰度相似性的數(shù)學(xué)量度。灰度相似性測(cè)度包括象素灰度的均方差、相關(guān)或互信息。模型驅(qū)動(dòng)方法首先建立

4、明確的幾何模型，以此表示解剖標(biāo)志。這些解剖標(biāo)志包括有重要功能的外表、曲線和點(diǎn)。將源系統(tǒng)的解剖標(biāo)志參數(shù)化，與目的系統(tǒng)的對(duì)應(yīng)部分對(duì)準(zhǔn)，以這種對(duì)應(yīng)關(guān)系引導(dǎo)系統(tǒng)其余部分的變換。模型驅(qū)動(dòng)算法包括點(diǎn)約束法、線約束法和面約束法?；旌纤惴ㄊ墙Y(jié)合使用以上兩種算法的方法。薄板樣條插值方法是非剛體變換中的一種特殊的變換，它允許部分調(diào)整，并符合某種連續(xù)性或平滑性要求。第2節(jié)討論剛性能量函數(shù)；第3節(jié)給出非剛性能量函數(shù)；第4節(jié)設(shè)計(jì)并實(shí)現(xiàn)一個(gè)非剛性配準(zhǔn)算法；最后給出實(shí)驗(yàn)結(jié)果。2剛性能量函數(shù)本研究之所以采用薄板樣條，是因?yàn)樗莫?dú)特性質(zhì)，就是可以將空間變換分解為一個(gè)全局仿射變換和一個(gè)部分非仿射變換。BksteEin4首先將薄板

5、樣條函數(shù)應(yīng)用于標(biāo)志點(diǎn)的匹配，結(jié)果證明它是一個(gè)非常有用的形狀分析工具。假設(shè)在二維空間，兩個(gè)具有N對(duì)對(duì)應(yīng)點(diǎn)的點(diǎn)集，Q=Qi，i=1，2，n和P=Pi，i=1，2，n，將點(diǎn)集Q，P表示為：Q=1x1y11x2y21xnynP=1x1y11x2y21xnyn下面我們建立從點(diǎn)集P到點(diǎn)集Q的薄板樣條映射f(Pi)，由于薄板樣條是不對(duì)稱的，因此從P到Q的映射不能簡(jiǎn)單地反轉(zhuǎn)為從Q到P的映射。通過最小化下面的能量函數(shù)，可以得到一個(gè)剛性能量函數(shù)：Etps(f)=ni=1Q-f(P)2+J(f)(1)其中，J(f)=R22fx22+22fxy2+2f2y2dxdy(1)式第一項(xiàng)代表經(jīng)過變換的源標(biāo)志點(diǎn)與目的標(biāo)志點(diǎn)之間

6、的間隔 和；第二項(xiàng)代表了獲得的變換的不平滑度，也叫懲罰函數(shù)。使該式最小化的變換既滿足變換后源標(biāo)志點(diǎn)與目的標(biāo)志點(diǎn)間接近(近似)的要求，同時(shí)也參加了足夠的平滑。系數(shù)(0)表征了近似和平滑之間的相對(duì)關(guān)系:當(dāng)較小時(shí)，獲得的變換表現(xiàn)了很好的近似效果;當(dāng)較大時(shí)，就獲得了比擬平滑的變換，對(duì)較大的部分畸變進(jìn)展了調(diào)整。薄板函數(shù)計(jì)算如下：設(shè)z(x,y)=-U(r)=-r2lgr2，其中，r=x2+y2，U(r)是構(gòu)建薄板樣條的基函數(shù)，設(shè)rij=|Pi-Pj|為點(diǎn)Pi與點(diǎn)Pj的歐幾里德間隔 。對(duì)分散點(diǎn)數(shù)據(jù)集Pi進(jìn)展薄板樣條彈性插值后可以得到曲面。插值過程形象地模擬為一個(gè)薄金屬板在點(diǎn)約束下的扭曲變形，要使金屬板在點(diǎn)(

7、xi，yi)處高度為zi，并且該板具有最小彎曲能量，即薄板函數(shù)f(x,y)使罰函數(shù)J(f)最校定義nn矩陣：K=0U(r12)U(r1n)U(r21)0U(r2n)U(rn1)U(rn2)0V=(z(x1,y1),z(x2,y2),z(xn,yn)T通過解線性方程組2可以得到=(1,2,n)T和T=(a1,ax,ay)TK+PT=VPT=0(2)是n3的非仿射變換形變參數(shù)矩陣，T是33的仿射形變參數(shù)矩陣，K是薄板樣條的核，為nn矩陣。然后構(gòu)造函數(shù)：f(x,y)=a1,axx+ayy+ni=1iU(|(xi,yi)-(x,y)|)3此時(shí)該函數(shù)對(duì)于所有i，有f(xi,yi)=zi，并使罰函數(shù)J(f

8、)最校事實(shí)上，直接解方程組2是困難的，也不現(xiàn)實(shí)，我們將通過迭代求解點(diǎn)集之間的匹配矩陣來求方程2的參數(shù)和T。3非剛性能量函數(shù)由剛性能量函數(shù)推導(dǎo)說明，只要兩個(gè)點(diǎn)集之間的對(duì)應(yīng)點(diǎn)，就可以得到它們之間的薄板樣條映射參數(shù)。但是當(dāng)對(duì)應(yīng)點(diǎn)未知時(shí)，該如何處理呢?傳統(tǒng)的方法往往都是手動(dòng)選點(diǎn),這種方法費(fèi)時(shí)費(fèi)力,同時(shí)在構(gòu)造不清的情況下,很難選擇到足夠多的準(zhǔn)確對(duì)應(yīng)點(diǎn)。而且其準(zhǔn)確性也只是相對(duì)的,誤差是不可防止的。文獻(xiàn)10定義兩個(gè)點(diǎn)集之間的匹配矩陣=ij：ij=1，假設(shè)點(diǎn)Qi對(duì)應(yīng)于點(diǎn)Pi0，其他由于兩個(gè)點(diǎn)集之間是雙向一一對(duì)應(yīng)的，即一個(gè)點(diǎn)集中的每個(gè)點(diǎn)在另一個(gè)點(diǎn)集中至多有一個(gè)對(duì)應(yīng)點(diǎn)，反之亦然。匹配矩陣ij具有下面約束：j，N1

9、i=1ij=1，i，N2j=1ij=14N1和N2分別是兩個(gè)點(diǎn)集的點(diǎn)數(shù)，將匹配矩陣考慮到式(1)中，得到基于薄板樣條映射的非剛性匹配的能量函數(shù)為：Etps(，T，)=N1i=1N2j=1ijQ-PT-K2+J(f)(5)式(5)的第一項(xiàng),可以使點(diǎn)集P中的點(diǎn)盡可能近地映射到Q中的點(diǎn)；第二項(xiàng)是平滑性約束,用于映射的調(diào)整，調(diào)整參數(shù)決定了映射的形變程度,當(dāng)0時(shí)，將得到對(duì)應(yīng)點(diǎn)的準(zhǔn)確匹配。4非剛性配準(zhǔn)算法的設(shè)計(jì)與實(shí)現(xiàn)在保證式(4)的約束下,放松對(duì)匹配矩陣的約束,將二值的匹配矩陣轉(zhuǎn)化為連續(xù)實(shí)數(shù)矩陣,即ij0,1ij0,1,允許部分匹配的存在,稱這樣的匹配矩陣為模糊匹配矩陣。由第下面的算法可以看到，隨著時(shí)間的

10、推移，ij的值逐漸變大，越來越接近二值矩陣，當(dāng)時(shí)間足夠長時(shí)，就會(huì)得到最終的二值匹配矩陣。根據(jù)匹配矩陣元素的性質(zhì)，令：ij=1Q-PT-K2+1(6)當(dāng)Q-PT-K0時(shí)，ij1。所以，非剛性能量函數(shù)4式可改寫為：Etps(，T，)=N1i=1N2j=1(ij2(Q-PT-K2+1)+2ij)+J(f)(7)通過Eij=0可以得到使能量函數(shù)(7)式極小的匹配矩陣元素ij。在第2節(jié)我們知道，直接解方程組2中參數(shù)和T是困難的，對(duì)于仿射變換T的計(jì)算是獨(dú)立于2式。因?yàn)槎S歐氏空間上的仿射變換可寫為：S(Pj)=TPj+A，其中，A=(x,y)T為平移量，T=kssin-sins,平移、旋轉(zhuǎn)、縮放及反射和剪

11、切等是二維仿射變換的特例。此模型中的參數(shù)k、x和y，即為兩圖像的配準(zhǔn)參數(shù)。確定這幾個(gè)參數(shù)的步驟為：首先對(duì)需配準(zhǔn)的兩幅圖像估計(jì)初始值k0、0、x0和y0，建立兩個(gè)點(diǎn)集的坐標(biāo)對(duì)應(yīng)關(guān)系。計(jì)算兩幅圖像對(duì)應(yīng)點(diǎn)互信息，可得到x0和y0。對(duì)k0、0的選取可先確定一個(gè)大致范圍，然后設(shè)定一定的間隔作步長因子，設(shè)k=1k0，=20，以互信息最大為原那么進(jìn)展迭代搜索，自適應(yīng)獲得最正確值。實(shí)際中圖像經(jīng)過預(yù)處理后，圖像之間的旋轉(zhuǎn)角比擬小，取值范圍為(-/4，/4)，就能保證找到正確的值。按照最大相關(guān)原那么迭代搜索，以獲得最正確值，對(duì)k0的處理方法與此類似。獲得了最正確值k和后，再對(duì)x0和y0按最大互信息原那么沿圖像兩個(gè)

12、正交方向逐像素搜索，以獲得最正確值x和y。確定了配準(zhǔn)參數(shù)K、x和y，就可對(duì)圖像進(jìn)展平移、旋轉(zhuǎn)和縮放。將得到的T代入2式，求得參數(shù)。的作用是將圖像經(jīng)T變換后坐標(biāo)值不落在像素點(diǎn)上的點(diǎn)調(diào)整到像素點(diǎn)上。本研究的算法主要包括以下幾個(gè)步驟：給定特征點(diǎn)集Q和P；初始化：ij=1(全1矩陣)，=1(全1矩陣),T=0零矩陣,A=0零向量，=0，N=100，構(gòu)造初始薄板樣條；根據(jù)仿射坐標(biāo)最大相關(guān)原那么迭代計(jì)算T，A；根據(jù)4式計(jì)算；根據(jù)6式計(jì)算ij；根據(jù)3式，構(gòu)造薄板樣條；假如滿足4式約束或迭代次數(shù)大于N，那么轉(zhuǎn)，否那么轉(zhuǎn)；end。需要指出的是,由于互信息量是統(tǒng)計(jì)量,因此我們對(duì)標(biāo)準(zhǔn)圖確定的點(diǎn)數(shù)不能太少，保證互信息

13、量的統(tǒng)計(jì)有效性。5實(shí)驗(yàn)結(jié)果圖1顯示了對(duì)一對(duì)9696,灰度級(jí)為256級(jí)的RI圖像進(jìn)展非剛性配準(zhǔn)的實(shí)例。(a)是一幅標(biāo)準(zhǔn)的RI心臟長軸的一個(gè)切片；(b)是一幅有變形的RI圖像心臟長軸同一周期另一個(gè)切片,此圖像是用圖像處理軟件Phtshp手動(dòng)變形處理獲得；()是標(biāo)準(zhǔn)圖像中特征點(diǎn)的選取,圖中的小紅圓圈為選擇的特征點(diǎn)；(d)是最大互信息搜索方法在變形圖像上搜索到的對(duì)應(yīng)特征點(diǎn)；(e)是配準(zhǔn)結(jié)果。由圖1(a)、(b)、(e)我們可以看出,經(jīng)過配準(zhǔn)變換后,圖1(b)變形為圖1(e),與圖1(a)已根本完全一致了。通過計(jì)算整幅圖像之間的互信息量我們得到:圖1(a)和圖1(b)的互信息量為1.037,圖1(a)和

14、圖1(e)之間的互信息量為3.44,而圖1(a)與自身計(jì)算得到的互信息量為3.46?？梢?經(jīng)過非剛性性配準(zhǔn)變換,圖像到達(dá)了較高的匹配精度。abde圖1RI非剛性配準(zhǔn)的實(shí)例6完畢語一直以來，學(xué)者們提出了各種非剛性圖像匹配方法。文獻(xiàn)5將圖像分解為許多子圖像，估計(jì)每個(gè)子圖像之間仿射變換，這樣用多個(gè)仿射變換近似反映整個(gè)圖像之間的非剛性映射。文獻(xiàn)6首先從圖像中提取出特征點(diǎn)，然后將特征點(diǎn)擬合成曲線或曲面進(jìn)展匹配。這種方法在擬合的曲線或曲面比擬光滑的情況下效果好，但是當(dāng)圖像所包含的形狀很復(fù)雜時(shí)，曲線或曲面的擬合就變得非常困難，這種方法的好處是曲線的匹配相對(duì)于點(diǎn)集的匹配要容易些。最近還有許多非剛性匹配的研究主

15、要集中于非剛性形狀的統(tǒng)計(jì)特性的學(xué)習(xí)上7，其方法類似于迭代最近點(diǎn)(IP)算法，是基于部分的啟發(fā)式搜索。以上這些方法大都存在魯棒性差，而且通常不能保證圖像之間的一一對(duì)應(yīng)和特征點(diǎn)的自動(dòng)確定。本研究提出了一種通過薄板樣條函數(shù)來表征特征點(diǎn)集之間的非剛性映射，把該映射分解為仿射變換和非仿射變換，并分別計(jì)算求解薄板樣條的參數(shù)并滿足雙向?qū)?yīng)的約束。參考文獻(xiàn)1LikarB,pernuiF.Ahierarhialapprahtelastiregistratinbasednutualinfratin.IagVisinput，2001，19:3344.2KyriauSK，Davatziks.Nnlinearelast

16、iregistratinfbrainiagesithturpathlgyusingabiehanialdelRI.IEEEtransedIag，1999，18:580592.3ButtA，AharyaR，Sibata，etal.Surfaeathingfultidalityiagevluesbyafuzzyelastiregistratintehnique.putedIagGraph，1998，22:1323.4BkstEinFL.Prinipalarps:Thinplatesplinesandthedepsitinfdefratins.IEEETarans.1989，PAI11(6):567

17、585.5PappuS,GldS.AFraerkfrNnRigidathingandrrespndene.AdvaneinNeuralinfratinPressingSystes,ITpress,abridge:1996,8.795801.6SzeliskiR,LavalleeS.athing3Danatisurfaesithnnrigiddefratinsusingtreesplines.Int1JputerVisin,1996,18:171186.7DutaN,JainAK,etal.Learning2Dshapedels.IEEEnf.nputerVisinandPatternRegnitin(VPR).lbia:Frtlins,1999,2：814.8SunDngei,QiuZhen