非剛性醫學圖像配準算法的設計與實現

導語：非剛性醫學圖像配準算法的設計與實現一文來源于網友上傳，不代表本站觀點，若需要原創文章可咨詢客服老師，歡迎參考。

【關鍵詞】醫學圖像；非剛性；圖像配準；匹配矩陣；薄板樣條

摘要：非剛性圖像匹配問題已成為醫學圖像分析中一個非常具有挑戰性的問題。基于薄板樣條插值方法,引入實匹配矩陣,并給出相應配準變換算法，該算法將薄板樣條參數表示成仿射分量和非仿射分量，并分別進行求解。與其它非剛性匹配算法相比，該算法不僅保證了對應特征點的雙向對應，也實現了自動特征點選擇，實驗結果令人滿意。

關鍵詞：醫學圖像；非剛性；圖像配準；匹配矩陣；薄板樣條

1引言

在醫學診斷和治療過程中，常需要對比分析多幅圖像，以獲得更為精確和全面的信息。圖像分析大都要求多幅圖像的幾何位置一致，因此，配準是醫學圖像分析的一個重大課題。醫學圖像配準是指對于一幅醫學圖像尋求一種(或一系列)空間變換，使它與另一幅醫學圖像上的對應點達到空間上的一致。這種一致是指人體上的同一解剖點在兩張匹配圖像上有相同的空間位置。配準的結果應使兩幅圖像上所有的解剖點，或至少是所有具有診斷意義的點及手術感興趣的點都達到匹配。圖像配準不僅可以校正病人多次成像間的位置變化，也可以校正由于成像模式本身導致的畸變。對同一個病人的不同時間的圖像進行配準，可以了解發育過程及腫瘤病變的病情；對不同人的圖像進行配準，去除種族、年齡等臨床及遺傳差異，從而形成疾病或人群特異性圖譜，可用于正常與否的分析；對不同成像模式進行配準，可以獲得互補信息。

醫學圖像配準可分為剛性配準和非剛性配準兩類。剛性配準在許多情況下不能滿足臨床的需要，因為很多形變的性質是非剛體、非線性的。比如為了精確定位MR圖像左心室，常常伴有組織磁化系數差異、非水分子的化學位移以及血流流動等因素導致的幾何畸變以及由于磁場不均勻、磁場梯度非線性及渦流等導致的探測畸變，因此在放療計劃制定中，將MR圖像配準時，不能單純地使用剛性配準，必須使用非剛性配準。

非剛性配準算法可分為灰度驅動、模型驅動及混合算法三種［1~3］。灰度驅動方法基于數學或統計尺度將一個灰度模式與另一個對準。典型情況下，需要定義源系統與目標系統之間的灰度相似性的數學量度。灰度相似性測度包括象素灰度的均方差、相關或互信息。模型驅動方法首先建立明確的幾何模型，以此表示解剖標志。這些解剖標志包括有重要功能的表面、曲線和點。將源系統的解剖標志參數化，與目標系統的對應部分對準，以這種對應關系引導系統其余部分的變換。模型驅動算法包括點約束法、線約束法和面約束法。混合算法是結合使用以上兩種算法的方法。薄板樣條插值方法是非剛體變換中的一種特殊的變換，它允許局部調整，并符合某種連續性或平滑性要求。第2節討論剛性能量函數；第3節給出非剛性能量函數；第4節設計并實現一個非剛性配準算法；最后給出實驗結果。

2剛性能量函數

本研究之所以采用薄板樣條，是因為它的獨特性質，就是能夠將空間變換分解為一個全局仿射變換和一個局部非仿射變換。Booksteein［4］首先將薄板樣條函數應用于標志點的匹配，結果證明它是一個非常有用的形狀分析工具。假設在二維空間，已知兩個具有N對對應點的點集，Q={Qi，i=1，2，…，n}和P={Pi，i=1，2，…，n}，將點集Q，P表示為：

Q=1x1y1

1x2y2

………

1xnynP=1x1y1

1x2y2

………

1xnyn

下面我們建立從點集P到點集Q的薄板樣條映射f(Pi)，由于薄板樣條是不對稱的，因此從P到Q的映射不能簡單地反轉為從Q到P的映射。通過最小化下面的能量函數，可以得到一個剛性能量函數：

Etps(f)=∑ni=1‖Q-f(P)‖2+λJ(f)(1)

其中，J(f)=R22fx22+22fxy2+2f2y2dxdy

(1)式第一項代表經過變換的源標志點與目標標志點之間的距離和；第二項代表了獲得的變換的不平滑度，也叫懲罰函數。使該式最小化的變換既滿足變換后源標志點與目標標志點間接近(近似)的要求，同時也加入了足夠的平滑。系數λ(λ>0)表征了近似和平滑之間的相對關系:當λ較小時，獲得的變換表現了很好的近似效果;當λ較大時，就獲得了比較平滑的變換，對較大的局部畸變進行了調整。薄板函數計算如下：

設z(x,y)=-U(r)=-r2logr2，其中，r=x2+y2，U(r)是構建薄板樣條的基函數，設rij=|Pi-Pj|為點Pi與點Pj的歐幾里德距離。對分散點數據集Pi進行薄板樣條彈性插值后可以得到曲面。插值過程形象地模擬為一個薄金屬板在點約束下的扭曲變形，要使金屬板在點(xi，yi)處高度為zi，并且該板具有最小彎曲能量，即薄板函數f(x,y)使罰函數J(f)最小。定義n×n矩陣：

K=0U(r12)…U(r1n)

U(r21)0…U(r2n)

…………

U(rn1)U(rn2)…0

V=(z(x1,y1),z(x2,y2),…,z(xn,yn))T

通過解線性方程組（2）可以得到W=(w1,w2,…,wn)T和T=(a1,ax,ay)T

KW+PT=V

PTW=0(2)

W是n×3的非仿射變換形變參數矩陣，T是3×3的仿射形變參數矩陣，K是薄板樣條的核，為n×n矩陣。

然后構造函數：

f(x,y)=a1,axx+ayy+∑ni=1wiU(|(xi,yi)-(x,y)|)（3）

此時該函數對于所有i，有f(xi,yi)=zi，并使罰函數J(f)最小。

事實上，直接解方程組（2）是困難的，也不現實，我們將通過迭代求解點集之間的匹配矩陣來求方程（2）的參數W和T。

3非剛性能量函數

由剛性能量函數推導表明，只要已知兩個點集之間的對應點，就可以得到它們之間的薄板樣條映射參數。但是當對應點未知時，該如何處理呢?傳統的方法往往都是手動選點,這種方法費時費力,同時在結構不清的情況下,很難選擇到足夠多的精確對應點。而且其準確性也只是相對的,誤差是不可避免的。文獻［10］定義兩個點集之間的匹配矩陣M={Mij}：

Mij=1，若點Qi對應于點Pi

0，其他

由于兩個點集之間是雙向一一對應的，即一個點集中的每個點在另一個點集中至多有一個對應點，反之亦然。匹配矩陣{Mij}具有下面約束：

j，∑N1i=1Mij=1，i，∑N2j=1Mij=1（4）

N1和N2分別是兩個點集的點數，將匹配矩陣考慮到式(1)中，得到基于薄板樣條映射的非剛性匹配的能量函數為：

Etps(M，T，W)=∑N1i=1∑N2j=1Mij‖Q-PT-KW‖2+λJ(f)(5)

式(5)的第一項,可以使點集P中的點盡可能近地映射到Q中的點；第二項是平滑性約束,用于映射的調整，調整參數λ決定了映射的形變程度,當λ→0時，將得到對應點的精確匹配。

4非剛性配準算法的設計與實現

在保證式(4)的約束下,放松對匹配矩陣的約束,將二值的匹配矩陣轉化為連續實數矩陣,即Mij∈{0,1}→Mij∈［0,1］,允許部分匹配的存在,稱這樣的匹配矩陣為模糊匹配矩陣。由第下面的算法可以看到，隨著時間的推移，Mij的值逐漸變大，越來越接近二值矩陣，當時間足夠長時，就會得到最終的二值匹配矩陣。

根據匹配矩陣元素的性質，令：

Mij=1‖Q-PT-KW‖2+1(6)

當‖Q-PT-KW‖→0時，Mij→1。所以，非剛性能量函數（4）式可改寫為：

Etps(M，T，W)=∑N1i=1∑N2j=1(Mij2(‖Q-PT-KW‖2+1)+2Mij)+λJ(f)(7)

通過EMij=0能夠得到使能量函數(7)式極小的匹配矩陣元素Mij。

在第2節我們知道，直接解方程組（2）中參數W和T是困難的，對于仿射變換T的計算是獨立于（2）式。

因為二維歐氏空間上的仿射變換可寫為：S(Pj)=TPj+A，其中，A=(Δx,Δy)T為平移量，T=kcosθsinθ

-sinθcosθ,平移、旋轉、縮放及反射和剪切等是二維仿射變換的特例。此模型中的參數k、θ、Δx和Δy，即為兩圖像的配準參數。確定這幾個參數的步驟為：首先對需配準的兩幅圖像估計初始值k0、θ0、Δx0和Δy0，建立兩個點集的坐標對應關系。計算兩幅圖像對應點互信息，可得到Δx0和Δy0。對k0、θ0的選取可先確定一個大致范圍，然后設定一定的間隔μ作步長因子，設k=μ1k0，θ=μ2θ0，以互信息最大為原則進行迭代搜索，自適應取得最佳值。實際中圖像經過預處理后，圖像之間的旋轉角θ比較小，取值范圍為(-π/4，π/4)，就能保證找到正確的θ值。按照最大相關原則迭代搜索，以獲得最佳值θ，對k0的處理方法與此類似。

獲得了最佳值k和θ后，再對Δx0和Δy0按最大互信息原則沿圖像兩個正交方向逐像素搜索，以取得最佳值Δx和Δy。

確定了配準參數K、θ、Δx和Δy，就可對圖像進行平移、旋轉和縮放。將得到的T代入（2）式，求得參數W。W的作用是將圖像經T變換后坐標值不落在像素點上的點調整到像素點上。

本研究的算法主要包括以下幾個步驟：

①給定特征點集Q和P；

②初始化：Mij=1(全1矩陣)，W=1(全1矩陣),T=0（零矩陣）,A=0（零向量），λ=λ0，N=100，構造初始薄板樣條；

③根據仿射坐標最大相關原則迭代計算T，A；

④根據（4）式計算W；

⑤根據（6）式計算Mij；

⑥根據（3）式，構造薄板樣條；

⑦如果M滿足（4）式約束或迭代次數大于>N，則轉⑧，否則轉③；

⑧end。

需要指出的是,由于互信息量是統計量,因此我們對標準圖確定的點數不能太少，保證互信息量的統計有效性。

5實驗結果

圖1顯示了對一對96×96,灰度級為256級的MRI圖像進行非剛性配準的實例。(a)是一幅標準的MRI心臟長軸的一個切片；(b)是一幅有變形的MRI圖像心臟長軸同一周期另一個切片,此圖像是用圖像處理軟件Photoshop手動變形處理獲得；(c)是標準圖像中特征點的選取,圖中的小紅圓圈為選擇的特征點；(d)是最大互信息搜索方法在變形圖像上搜索到的對應特征點；(e)是配準結果。由圖1(a)、(b)、(e)我們可以看出,經過配準變換后,圖1(b)變形為圖1(e),與圖1(a)已基本完全一致了。通過計算整幅圖像之間的互信息量我們得到:圖1(a)和圖1(b)的互信息量為1.037,圖1(a)和圖1(e)之間的互信息量為3.44,而圖1(a)與自身計算得到的互信息量為3.46。可見,經過非剛性性配準變換,圖像達到了較高的匹配精度。

abcde

圖1MRI非剛性配準的實例

6結束語

一直以來，學者們提出了各種非剛性圖像匹配方法。文獻［5］將圖像分解為許多子圖像，估計每個子圖像之間仿射變換，這樣用多個仿射變換近似反映整個圖像之間的非剛性映射。文獻［6］首先從圖像中提取出特征點，然后將特征點擬合成曲線或曲面進行匹配。這種方法在擬合的曲線或曲面比較光滑的情況下效果好，但是當圖像所包含的形狀很復雜時，曲線或曲面的擬合就變得非常困難，這種方法的好處是曲線的匹配相對于點集的匹配要容易些。

最近還有許多非剛性匹配的研究主要集中于非剛性形狀的統計特性的學習上［7］，其方法類似于迭代最近點(ICP)算法，是基于局部的啟發式搜索。以上這些方法大都存在魯棒性差，而且通常不能保證圖像之間的一一對應和特征點的自動確定。本研究提出了一種通過薄板樣條函數來表征特征點集之間的非剛性映射，把該映射分解為仿射變換和非仿射變換，并分別計算求解薄板樣條的參數并滿足雙向對應的約束。

參考文獻

1LikarB,pernuiF.Ahierarchicalapproachtoelasticregistrationbasedonmutualinformation.ImagVisionComput，2001，19:33~44.

2KyriacouSK，DavatzikosC.Nonlinearelasticregistrationofbrainimageswithtumorpathologyusingabiomechanicalmodel［MRI］.IEEEtransMedImag，1999，18:580~592.

3ButtA，AcharyaR，SibataC，putMedImagGraph，1998，22:13~23.

4BooksteinFL.Principalwarps:Thinplatesplinesandthedecompositionofdeformations.IEEETarans.1989，PAMI11(6):567~585.

5PappuS,GoldS.AFrameworkforNonRigidMatchingandCorrespondence.AdvanceinNeuralinformationProcessingSystems,MITpress,Cambridge:1996,8.795~801.

6SzeliskiR,LavalleeS.Matching3Danatomicsurfaceswithnonrigiddeformationsusingoctreesplines.Int1JComputerVision,1996,18:171~186.

7DutaN,JainAK,etal.Learning2Dshapemodels.IEEEConf.OnComputerVisionandPatternRecognition(CVPR).Colombia:FortColins,1999,2：8~14.

8SunDongmei,QiuZhengding.ANewNonRigidImageMatchingAlgorithmUsingThinPlatesSpline.ACTAElectronicSinica(inChina).2002,30(8):1104~1108.

9RenHaihua，etal.MedicalImageElasticRegistrationBasedonThinPlateSplines.BeijingBiomedicalEngineering(inChnia).2003,22(4):241~245.

10孫冬梅,裘正定.利用薄板樣條函數實現非剛性圖像匹配算法.電子學報，2002,30(8):1104~1107.