非剛性醫學圖像配準是醫學圖像研究領域的熱門專題之一,具有重要的臨床應用價值。本文提出一種改進的Demons算法,將灰度守恒模型和局部結構張量守恒模型結合起來,構造一個新的能量函數處理多模態配準問題,然后采用L-BFGS算法優化能量函數,解決復雜三維數據的優化問題,并采用多尺度分層細化的思想解決大形變配準。實驗表明,本文算法對大形變和多模態三維醫學圖像配準都有較好的效果。
引用本文： 郝培博, 陳震, 江少鋒, 汪洋. 基于Demons算法的多模態醫學圖像非剛性配準研究. 生物醫學工程學雜志, 2014, 31(1): 161-165. doi: 10.7507/1001-5515.20140032 復制
引言
醫學圖像的配準已經成為研究領域的熱門之一。醫學圖像配準是指通過尋找某種空間變換，使兩幅圖像的對應點達到空間位置和解剖結構上的完全一致，要求配準的結果能使兩幅圖像上所有的解剖點，或至少是所有具有診斷意義以及手術區域的點都達到匹配[1]。醫學圖像配準具有很重要的臨床應用價值。
隨著醫學圖像學的不斷發展，剛性配準已經不能滿足臨床的需要，越來越多的學者開始研究非剛性醫學圖像配準。研究人員提出了大量的非剛性醫學圖像配準算法。這些算法可分為二大類：基于特征點的算法、基于物理模型的算法[2]。
基于特征點的算法首先要在每個需要配準的數據集上建立代表不同解剖元素的幾何模型，然后與目標圖像中的對應解剖元素進行配準，最后利用解剖元素之間配準得到的對應關系引導圖像其余部分之間的配準。Likar等[3]提出將圖像劃分為若干子塊，取每個子塊的中心作為標記對應點，但沒有旋轉變換，提取的標記點有偏差，配準精度有待提高。馮林等[4]提出了一種基于分層互信息和薄板樣條自動確定標記對應點的選取方法，但是由于該算法依賴灰度統計相關性，無法取得足夠多的標記點。基于物理模型的算法通常采用來自數學或統計學的相似性準則，以匹配圖像每個區域的灰度分布模式為目標，通常會在形變的模板圖像和目標圖像之間定義一個灰度相似性測度，最后通過物理模型變換得到配準圖像。D’Agostino等[5]用粘滯流體模型控制形變，該算法可以實現任何復雜形狀的形變，但是所需計算量大。Lester等[6]在流體模型中，把粘滯常數看成可以變化的，因此對于圖像不同的部分可以有不同的形變特性。
Demons算法是一種小形變無參數非剛性配準方法[7]。它具有較高的運算效率和良好的配準效果。Hellier等[8]對常用的六種配準算法的回顧性對比研究表明，Demons算法在多個評價指標中都優于其他五種算法。但是Demons算法在處理大形變配準和多模態配準時候，會產生較大的誤配準和物理上的不合理變形。因此，本研究旨在對Demons算法進行一定的改進，使其能在大形變和多模態下具有較高的準確性。
1 Demons算法原理及改進
Demons配準驅動力是以光流方程式驅動圖像小形變變化。對靜態圖像F一個給定點p，用f表示靜態圖像的p點的強度，m表示浮動圖像的p點的強度，那么浮動圖像到靜態圖像配準點p的估計偏移量u為
| $~u=\frac{\left( m-f \right)\nabla f}{{{\left| \nabla f \right|}^{2}}+{{\left( m-f \right)}^{2}}},$ |
其中u=(ux，uy)是浮動圖像上任意點p到靜態圖像上對應點的偏移量，?f是靜態圖像的梯度值，(m-f)2是為了梯度方程更加穩定地添加到方程中。Wang 等[9]增加了一個方程使用浮動圖像邊緣作為內在驅動力，改進了算法的收斂速度和穩定性，
| $u=\frac{\left( m-f \right)\nabla f}{{{\left| \nabla f \right|}^{2}}+{{\partial }^{2}}{{\left( m-f \right)}^{2}}}+\frac{\left( m-f \right)\nabla m}{{{\left| \nabla m \right|}^{2}}+{{\partial }^{2}}{{\left( m-f \right)}^{2}}}$ |
為了調整驅動力的強度,Cachier等[10]提出了歸一化因子?。
Vercauteren等[11]描繪了標準配準模型；一個基于能量的相似性測度，求得一個配準誤差函數，當誤差函數為0時求得最優解。相似性測度是衡量變換后圖像m(x-u)好壞的重要標準。基于灰度守恒的配準中最常用的相似性測度是像素距離平方：
| $E\left( m,f,u \right)=\frac{1}{2}\int{{{\left| m\left( x-u \right)-f\left( x \right) \right|}^{2}}dx,}$ |
所以可以計算出誤差函數：
| $\nabla E\left( u \right)=2{{\left( \nabla m \right)}^{T}}\left( f-m+u\nabla m \right)$ |
若直接對式(4) 求最小化會產生不惟一解，那優化問題就存在不適定性，因此需要對式(4) 增加正則項來解決不適定性，但是Thirion提出的Demons算法只是通過簡單的高斯平滑來正則化，使圖像結構丟失，圖像出現變形。
雖然圖像的梯度方向也包含了局部結構信息，但是在計算圖像的局部方向時需要在每個像素的鄰域做平均，如果在某個像素的局部鄰域中同時包含上升邊緣的梯度和下降邊緣的梯度(如一條細邊緣的兩側)，由于此時梯度出現正負抵消，因此無法反映該點的局部方向。局部結構張量很好地避免了梯度方向的這一不足，同時能夠區分圖像中的各種不同結構，近年來被廣泛應用于圖像處理領域。局部結構張量守恒假設是根據圖像局部結構在運動方向上基本不變的特性提出的，這樣就可以避免運動方向性的約束。 形式為
| $T\left( X+w,t+1 \right)-T\left( X,t \right)=0,$ |
式(5) 中T為圖像的局部結構張量，包含圖像的局部運動信息(任意方向運動信息)。
| $T=\left[ \begin{matrix} I_{x}^{2} & {{I}_{x}}{{I}_{y}} \\ {{I}_{x}}{{I}_{y}} & I_{y}^{2} \\ \end{matrix} \right]$ |
在設計能量函數時，可以根據圖像中運動場景的不同和各守恒假設適用的范圍，合理地選擇其中一種或者幾種守恒假設結合的方式，需注意的是數據項中選擇的守恒假設越多，越有可能導致模型計算復雜度的提高和效率的降低，因此通常選擇兩種假設模型的結合方式構成數據項模型。
綜合考慮灰度守恒假設與圖像局部結構守恒假設的適用范圍和優缺點，將兩者結合建立一個新的能量函數[12]，得到的能量函數(相似性測度)為
| $\eqalign{ & {E_{data}} = {E_g} + {\lambda _1}{E_t} = \cr & \int_\Omega {\left( {{{\left| {I\left( {X + W} \right) - I\left( X \right)} \right|}^2} + {\lambda _1}*{{\left| {T\left( {X + W} \right) - T\left( X \right)} \right|}^2}} \right)dX} \cr & {\int_\Omega {\left( {\left| {{m_x}\left( {x - u} \right) - f\left( x \right)} \right|} \right.} ^2} + {\lambda _1}* \cr & \left[ \matrix{ {\left| {{m_x}\left( {x - u} \right) - f\left( x \right)} \right|^2} + {\left| {{m_y}\left( {x - u} \right) - {f_y}\left( x \right)} \right|^2} + \hfill \cr {\left| {{m_x}\left( {x - u} \right)*{m_y}\left( {x - u} \right) - {f_x}\left( x \right){f_y}\left( x \right)} \right|^2} \hfill \cr} \right]\left. {} \right)dX, \cr} $ |
其中λ1表示數據項中兩守恒假設的權重系數。則新的誤差函數為
| $\begin{align} & \nabla E\left( u \right)=2{{\left( \nabla m \right)}^{T}}\left( f-m+u\nabla m \right)+2{{\left( \nabla {{m}_{x}} \right)}^{T}} \\ & \left( {{f}_{x}}-{{m}_{x}}+u\nabla {{m}_{x}} \right)+2{{\left( \nabla {{m}_{y}} \right)}^{T}}\left( {{f}_{y}}-{{m}_{y}}+u\nabla {{m}_{y}} \right) \\ & +2{{\left( \nabla {{m}_{x}} \right)}^{T}}{{\left( \nabla {{m}_{y}} \right)}^{T}}\left( {{f}_{x}}{{f}_{y}}-{{m}_{x}}{{m}_{y}}+2u\nabla {{m}_{x}}\nabla {{m}_{y}} \right) \\ \end{align}$ |
其中mx、my為浮動圖像在x和y方向的梯度，fx、fy為參考圖像在x和y方向的梯度。
為了將算法應用到三維圖像中，可以把二維局部結構張量擴展到三維，三維結構張量包含了三維空間中的典型結構。設I(X)為圖像在點X(x,y,z)的灰度，此時的局部結構張量為
| $T=\left[ \begin{matrix} I_{x}^{2} & {{I}_{x}}{{I}_{y}} & {{I}_{x}}{{I}_{z}} \\ {{I}_{x}}{{I}_{y}} & I_{y}^{2} & {{I}_{y}}{{I}_{z}} \\ {{I}_{x}}{{I}_{z}} & {{I}_{y}}{{I}_{z}} & I_{z}^{2} \\ \end{matrix} \right]$ |
三維圖像數據量大，如果不對三維數據進行優化處理，配準過程相當耗時，本文采用L-BFGS算法優化能量函數，此算法將BFGS方法與有限內存技術相結合，用來求解大規模優化問題[13]。其特點是不需要明確地形成或存儲Hessian矩陣，只保存x處前m次的更新歷史以及梯度，是一種有限內存擬牛頓法。
Demons算法對于大形變非剛性配準不能得到良好效果，為了改善這種情況，采用多尺度分層細化思想，在同一幅圖像中，較高分辨率下可以觀察到尺寸很小或者灰度級不高的物體，而在較低的分辨率下卻只能觀察到尺寸很大或者灰度級較高的物體，對圖像采用多分辨率處理具有明顯的優勢。因此，這樣的從粗到精的多分辨率思想在大位移非剛性運動時特別適用。
2 改進算法的實現
改進算法的具體步驟：
(1) 將浮動圖像m和靜態圖像f進行多尺度金字塔分層。
(2) 從金字塔最高層k層開始，將mk,fk代入式(2) 中，計算出形變場uk。
(3) 由浮動圖像mk和形變場uk得到新的浮動圖像mk(x-u)。
(4) 將mk(x-u)和fk代入能量函數式(7) 中，判斷是否滿足相似性測度。若不滿足需將mk(x-u)和fk重新代入誤差函數式(8) 中，計算出新的形變場uk，直到滿足相似性測度；若滿足則將形變場uk插值得到下一層的初始形變場uk-1。
(5) 將mk-1，fk-1和uk-1重新代入步驟3中迭代，直到計算到金字塔底層。
(6) 將底層得到的形變場u0代入浮動圖像m，得到最終配準后圖像。
3 實驗結果及討論
實驗目的是為了驗證本算法在不同模態情況下和大形變情況下的配準效果優于Demons算法。
3.1 配準精度的評價方法
目前醫學圖像配準還沒有金標準來衡量配準結果的好壞，在這里采用主觀評價，最小均方誤差(MSE)和相關系數(Rcc)來評價算法的優劣。
MSE的計算公式為
| $MSE = \sqrt {\sum\limits_{x = 0}^{N - 1} {\sum\limits_{y = 0}^{N - 1} {{{\left[ {M\left( {x,y} \right) - F\left( {x,y} \right)} \right]}^2}} } } ,$ |
MSE 越小表示配準效果越好。
Rcc的計算公式為
| ${{R}_{cc}}=\frac{\frac{1}{N-1}\sum\limits_{i=0}^{N-1}{\left( {{M}_{i}}-\bar{M} \right)\left( {{F}_{i}}-\bar{F} \right)}}{\sqrt{\frac{1}{N-1}\sum\limits_{i=0}^{N-1}{{{\left( {{M}_{i}}-\bar{M} \right)}^{2}}}}\sqrt{\frac{1}{N-1}\sum\limits_{i=0}^{N-1}{{{\left( {{F}_{i}}-\bar{F} \right)}^{2}}}}},$ |
其中N為像素點個數，
3.2 實驗結果
實驗一：驗證算法在二維圖像配準中也能取得良好結果，采用256×256像素的顱腦MR圖像，用質子密度加權成像和T1加權成像代表不同模態的圖像。實驗結果如圖 1所示，表 1為本文算法與Demons算法等的比較。
圖1
二維圖像配準結果
(a)浮動圖像;(b)參考圖像;(c)互信息算法結果;(d)B樣條算法結果;(e)Demons算法結果;(f)本文算法結果
Figure1. 2D image registration results(a) moving image; (b) static image; (c) mutual information algorithm result; (d) B-spline algorithm result; (e) Demons algorithm result; (f) the proposed algorithm result
表1
不同算法的配準精度比較
Table1.
Registration accuracy comparison among different algorithms
從圖 1中可知，(a)為質子密度加權成像圖像偏移了13x17y，(b)為T1加權成像圖像，(c)為互信息算法配準后圖像，(d)為B樣條算法配準后圖像，(e)為Demons算法配準后圖像，(f)為本文算法配準后
的圖像。可以直觀看到，B樣條算法和Demons算法在處理不同模態圖像配準時，出現了嚴重的變形和失真；而本文算法和互信息算法能很好地處理多模態情況下的配準。從表 1中可知，本文算法的MSE小于其他幾種算法，而且Rcc都高于其他幾種算法，說明本文算法優于其他幾種算法，能夠很好地處理多模態圖像配準。
實驗二：驗證算法在三維配準中的應用，采用模擬MRT腦數據網站(http://brainweb.bic.mni.mcgill.ca/brainweb/)中的模擬MR腦數據，選取mr_T1數據為181×217×181，切片厚度為1 mm，含有3%噪聲的圖像序列繞垂直軸旋轉10°作為浮動圖像；選取質子密度加權成像mr_PD數據為181×217×181，切片厚度為1 mm，含有3%噪聲的圖像序列作為參考圖像。配準后的結果如
圖2
三維多模態配準
(a)浮動圖像;(b)參考圖像;(c)配準后圖像;(d)配準前后第90張切片圖像差值
Figure2. 3D image registration results(a) moving image; (b) static image; (c) registration result; (d) image difference of the 90th slice before and after registration
從圖 2中可知，(a)為T1數據181×217×181繞垂直軸旋轉10°并偏移5x5y，(b)為質子密度加權成像數據181×217×181，(c)為配準后T1數據，從直觀上可以看到矢狀面和冠狀面與參考圖像基本一致，水平面上的偏移量和旋轉角度都得到了較好的配準，本文算法對多模態三維圖像配準也能夠得到良好的效果。目前還沒有衡量圖像三維配準的金標準，故采用配準前后與參考圖像的差值比較來表示配準結果的好壞。如圖 2(d)所示，(d)中左邊是配準前圖像差值，配準前有明顯的偏移和誤差；右邊是配準后圖像差值，偏移和誤差有明顯的改善。本文算法對不同模態的三維醫學圖像也能取得較好的配準結果。
實驗三：選取mr_T1數據為181×217×181，切片厚度為1 mm，含有3%噪聲的圖像序列繞垂直軸旋轉10°作為浮動圖像；選取mr_T1數據為181×217×60，切片厚度為3 mm，含有3%噪聲的圖像序列作為參考圖像。配準結果如圖 3所示。
圖3
不同切片厚度的三維配準
(a)浮動圖像;(b)參考圖像;(c)配準后圖像;(d)配準前后第90張切片圖像差值
Figure3. 3D registration with different slice thicknessesl(a) moving image; (b) static image; (c) registration result; (d) image difference of the 90th slice before and after registration
如圖 3所示，(a)為T1數據181×217×181繞垂直軸旋轉10°并偏移5x5y，(b)為T1數據181×217×60，切片厚度為3 mm，由于切片厚度不同，在配準前先進行部分體積差值，使圖像尺寸相同再進行配準。(c)為配準后的數據，從直觀上可以看到冠狀面和矢狀面的灰度值有所變化，這是由部分體積差值引起的，水平面上選擇和偏移都得到了良好的改善，配準前后的差異如(d)所示，左邊是配準前圖像差值，右邊是配準后圖像差值。從(d)中可以看到，配準前兩幅圖像差異非常明顯，有較多的偽影和偏移；配準后圖像得到了較大的改善，偽影和偏移明顯減少。本文算法在處理不同尺寸的體素配準時也能取得較好的結果。
4 結論
本文從多種能量守恒相結合的思想出發，把局部結構張量守恒假設和灰度守恒相結合，得到新的能量函數，使算法能適應不同模態圖像的配準。將結構張量應用到能量函數中，使圖像的結構得以保持，以適應大形變的圖像配準。在算法實現中，采用多尺度分層細化策略，使圖像運動信息更加準確。實驗表明，本文算法對大形變配準、多模態配準都能得到較好的結果，并能得到比Demons算法精確度更好的結果。
引言
醫學圖像的配準已經成為研究領域的熱門之一。醫學圖像配準是指通過尋找某種空間變換，使兩幅圖像的對應點達到空間位置和解剖結構上的完全一致，要求配準的結果能使兩幅圖像上所有的解剖點，或至少是所有具有診斷意義以及手術區域的點都達到匹配[1]。醫學圖像配準具有很重要的臨床應用價值。
隨著醫學圖像學的不斷發展，剛性配準已經不能滿足臨床的需要，越來越多的學者開始研究非剛性醫學圖像配準。研究人員提出了大量的非剛性醫學圖像配準算法。這些算法可分為二大類：基于特征點的算法、基于物理模型的算法[2]。
基于特征點的算法首先要在每個需要配準的數據集上建立代表不同解剖元素的幾何模型，然后與目標圖像中的對應解剖元素進行配準，最后利用解剖元素之間配準得到的對應關系引導圖像其余部分之間的配準。Likar等[3]提出將圖像劃分為若干子塊，取每個子塊的中心作為標記對應點，但沒有旋轉變換，提取的標記點有偏差，配準精度有待提高。馮林等[4]提出了一種基于分層互信息和薄板樣條自動確定標記對應點的選取方法，但是由于該算法依賴灰度統計相關性，無法取得足夠多的標記點。基于物理模型的算法通常采用來自數學或統計學的相似性準則，以匹配圖像每個區域的灰度分布模式為目標，通常會在形變的模板圖像和目標圖像之間定義一個灰度相似性測度，最后通過物理模型變換得到配準圖像。D’Agostino等[5]用粘滯流體模型控制形變，該算法可以實現任何復雜形狀的形變，但是所需計算量大。Lester等[6]在流體模型中，把粘滯常數看成可以變化的，因此對于圖像不同的部分可以有不同的形變特性。
Demons算法是一種小形變無參數非剛性配準方法[7]。它具有較高的運算效率和良好的配準效果。Hellier等[8]對常用的六種配準算法的回顧性對比研究表明，Demons算法在多個評價指標中都優于其他五種算法。但是Demons算法在處理大形變配準和多模態配準時候，會產生較大的誤配準和物理上的不合理變形。因此，本研究旨在對Demons算法進行一定的改進，使其能在大形變和多模態下具有較高的準確性。
1 Demons算法原理及改進
Demons配準驅動力是以光流方程式驅動圖像小形變變化。對靜態圖像F一個給定點p，用f表示靜態圖像的p點的強度，m表示浮動圖像的p點的強度，那么浮動圖像到靜態圖像配準點p的估計偏移量u為
| $~u=\frac{\left( m-f \right)\nabla f}{{{\left| \nabla f \right|}^{2}}+{{\left( m-f \right)}^{2}}},$ |
其中u=(ux，uy)是浮動圖像上任意點p到靜態圖像上對應點的偏移量，?f是靜態圖像的梯度值，(m-f)2是為了梯度方程更加穩定地添加到方程中。Wang 等[9]增加了一個方程使用浮動圖像邊緣作為內在驅動力，改進了算法的收斂速度和穩定性，
| $u=\frac{\left( m-f \right)\nabla f}{{{\left| \nabla f \right|}^{2}}+{{\partial }^{2}}{{\left( m-f \right)}^{2}}}+\frac{\left( m-f \right)\nabla m}{{{\left| \nabla m \right|}^{2}}+{{\partial }^{2}}{{\left( m-f \right)}^{2}}}$ |
為了調整驅動力的強度,Cachier等[10]提出了歸一化因子?。
Vercauteren等[11]描繪了標準配準模型；一個基于能量的相似性測度，求得一個配準誤差函數，當誤差函數為0時求得最優解。相似性測度是衡量變換后圖像m(x-u)好壞的重要標準。基于灰度守恒的配準中最常用的相似性測度是像素距離平方：
| $E\left( m,f,u \right)=\frac{1}{2}\int{{{\left| m\left( x-u \right)-f\left( x \right) \right|}^{2}}dx,}$ |
所以可以計算出誤差函數：
| $\nabla E\left( u \right)=2{{\left( \nabla m \right)}^{T}}\left( f-m+u\nabla m \right)$ |
若直接對式(4) 求最小化會產生不惟一解，那優化問題就存在不適定性，因此需要對式(4) 增加正則項來解決不適定性，但是Thirion提出的Demons算法只是通過簡單的高斯平滑來正則化，使圖像結構丟失，圖像出現變形。
雖然圖像的梯度方向也包含了局部結構信息，但是在計算圖像的局部方向時需要在每個像素的鄰域做平均，如果在某個像素的局部鄰域中同時包含上升邊緣的梯度和下降邊緣的梯度(如一條細邊緣的兩側)，由于此時梯度出現正負抵消，因此無法反映該點的局部方向。局部結構張量很好地避免了梯度方向的這一不足，同時能夠區分圖像中的各種不同結構，近年來被廣泛應用于圖像處理領域。局部結構張量守恒假設是根據圖像局部結構在運動方向上基本不變的特性提出的，這樣就可以避免運動方向性的約束。 形式為
| $T\left( X+w,t+1 \right)-T\left( X,t \right)=0,$ |
式(5) 中T為圖像的局部結構張量，包含圖像的局部運動信息(任意方向運動信息)。
| $T=\left[ \begin{matrix} I_{x}^{2} & {{I}_{x}}{{I}_{y}} \\ {{I}_{x}}{{I}_{y}} & I_{y}^{2} \\ \end{matrix} \right]$ |
在設計能量函數時，可以根據圖像中運動場景的不同和各守恒假設適用的范圍，合理地選擇其中一種或者幾種守恒假設結合的方式，需注意的是數據項中選擇的守恒假設越多，越有可能導致模型計算復雜度的提高和效率的降低，因此通常選擇兩種假設模型的結合方式構成數據項模型。
綜合考慮灰度守恒假設與圖像局部結構守恒假設的適用范圍和優缺點，將兩者結合建立一個新的能量函數[12]，得到的能量函數(相似性測度)為
| $\eqalign{ & {E_{data}} = {E_g} + {\lambda _1}{E_t} = \cr & \int_\Omega {\left( {{{\left| {I\left( {X + W} \right) - I\left( X \right)} \right|}^2} + {\lambda _1}*{{\left| {T\left( {X + W} \right) - T\left( X \right)} \right|}^2}} \right)dX} \cr & {\int_\Omega {\left( {\left| {{m_x}\left( {x - u} \right) - f\left( x \right)} \right|} \right.} ^2} + {\lambda _1}* \cr & \left[ \matrix{ {\left| {{m_x}\left( {x - u} \right) - f\left( x \right)} \right|^2} + {\left| {{m_y}\left( {x - u} \right) - {f_y}\left( x \right)} \right|^2} + \hfill \cr {\left| {{m_x}\left( {x - u} \right)*{m_y}\left( {x - u} \right) - {f_x}\left( x \right){f_y}\left( x \right)} \right|^2} \hfill \cr} \right]\left. {} \right)dX, \cr} $ |
其中λ1表示數據項中兩守恒假設的權重系數。則新的誤差函數為
| $\begin{align} & \nabla E\left( u \right)=2{{\left( \nabla m \right)}^{T}}\left( f-m+u\nabla m \right)+2{{\left( \nabla {{m}_{x}} \right)}^{T}} \\ & \left( {{f}_{x}}-{{m}_{x}}+u\nabla {{m}_{x}} \right)+2{{\left( \nabla {{m}_{y}} \right)}^{T}}\left( {{f}_{y}}-{{m}_{y}}+u\nabla {{m}_{y}} \right) \\ & +2{{\left( \nabla {{m}_{x}} \right)}^{T}}{{\left( \nabla {{m}_{y}} \right)}^{T}}\left( {{f}_{x}}{{f}_{y}}-{{m}_{x}}{{m}_{y}}+2u\nabla {{m}_{x}}\nabla {{m}_{y}} \right) \\ \end{align}$ |
其中mx、my為浮動圖像在x和y方向的梯度，fx、fy為參考圖像在x和y方向的梯度。
為了將算法應用到三維圖像中，可以把二維局部結構張量擴展到三維，三維結構張量包含了三維空間中的典型結構。設I(X)為圖像在點X(x,y,z)的灰度，此時的局部結構張量為
| $T=\left[ \begin{matrix} I_{x}^{2} & {{I}_{x}}{{I}_{y}} & {{I}_{x}}{{I}_{z}} \\ {{I}_{x}}{{I}_{y}} & I_{y}^{2} & {{I}_{y}}{{I}_{z}} \\ {{I}_{x}}{{I}_{z}} & {{I}_{y}}{{I}_{z}} & I_{z}^{2} \\ \end{matrix} \right]$ |
三維圖像數據量大，如果不對三維數據進行優化處理，配準過程相當耗時，本文采用L-BFGS算法優化能量函數，此算法將BFGS方法與有限內存技術相結合，用來求解大規模優化問題[13]。其特點是不需要明確地形成或存儲Hessian矩陣，只保存x處前m次的更新歷史以及梯度，是一種有限內存擬牛頓法。
Demons算法對于大形變非剛性配準不能得到良好效果，為了改善這種情況，采用多尺度分層細化思想，在同一幅圖像中，較高分辨率下可以觀察到尺寸很小或者灰度級不高的物體，而在較低的分辨率下卻只能觀察到尺寸很大或者灰度級較高的物體，對圖像采用多分辨率處理具有明顯的優勢。因此，這樣的從粗到精的多分辨率思想在大位移非剛性運動時特別適用。
2 改進算法的實現
改進算法的具體步驟：
(1) 將浮動圖像m和靜態圖像f進行多尺度金字塔分層。
(2) 從金字塔最高層k層開始，將mk,fk代入式(2) 中，計算出形變場uk。
(3) 由浮動圖像mk和形變場uk得到新的浮動圖像mk(x-u)。
(4) 將mk(x-u)和fk代入能量函數式(7) 中，判斷是否滿足相似性測度。若不滿足需將mk(x-u)和fk重新代入誤差函數式(8) 中，計算出新的形變場uk，直到滿足相似性測度；若滿足則將形變場uk插值得到下一層的初始形變場uk-1。
(5) 將mk-1，fk-1和uk-1重新代入步驟3中迭代，直到計算到金字塔底層。
(6) 將底層得到的形變場u0代入浮動圖像m，得到最終配準后圖像。
3 實驗結果及討論
實驗目的是為了驗證本算法在不同模態情況下和大形變情況下的配準效果優于Demons算法。
3.1 配準精度的評價方法
目前醫學圖像配準還沒有金標準來衡量配準結果的好壞，在這里采用主觀評價，最小均方誤差(MSE)和相關系數(Rcc)來評價算法的優劣。
MSE的計算公式為
| $MSE = \sqrt {\sum\limits_{x = 0}^{N - 1} {\sum\limits_{y = 0}^{N - 1} {{{\left[ {M\left( {x,y} \right) - F\left( {x,y} \right)} \right]}^2}} } } ,$ |
MSE 越小表示配準效果越好。
Rcc的計算公式為
| ${{R}_{cc}}=\frac{\frac{1}{N-1}\sum\limits_{i=0}^{N-1}{\left( {{M}_{i}}-\bar{M} \right)\left( {{F}_{i}}-\bar{F} \right)}}{\sqrt{\frac{1}{N-1}\sum\limits_{i=0}^{N-1}{{{\left( {{M}_{i}}-\bar{M} \right)}^{2}}}}\sqrt{\frac{1}{N-1}\sum\limits_{i=0}^{N-1}{{{\left( {{F}_{i}}-\bar{F} \right)}^{2}}}}},$ |
其中N為像素點個數，
3.2 實驗結果
實驗一：驗證算法在二維圖像配準中也能取得良好結果，采用256×256像素的顱腦MR圖像，用質子密度加權成像和T1加權成像代表不同模態的圖像。實驗結果如圖 1所示，表 1為本文算法與Demons算法等的比較。
圖1
二維圖像配準結果
(a)浮動圖像;(b)參考圖像;(c)互信息算法結果;(d)B樣條算法結果;(e)Demons算法結果;(f)本文算法結果
Figure1. 2D image registration results(a) moving image; (b) static image; (c) mutual information algorithm result; (d) B-spline algorithm result; (e) Demons algorithm result; (f) the proposed algorithm result
表1
不同算法的配準精度比較
Table1.
Registration accuracy comparison among different algorithms
從圖 1中可知，(a)為質子密度加權成像圖像偏移了13x17y，(b)為T1加權成像圖像，(c)為互信息算法配準后圖像，(d)為B樣條算法配準后圖像，(e)為Demons算法配準后圖像，(f)為本文算法配準后
的圖像。可以直觀看到，B樣條算法和Demons算法在處理不同模態圖像配準時，出現了嚴重的變形和失真；而本文算法和互信息算法能很好地處理多模態情況下的配準。從表 1中可知，本文算法的MSE小于其他幾種算法，而且Rcc都高于其他幾種算法，說明本文算法優于其他幾種算法，能夠很好地處理多模態圖像配準。
實驗二：驗證算法在三維配準中的應用，采用模擬MRT腦數據網站(http://brainweb.bic.mni.mcgill.ca/brainweb/)中的模擬MR腦數據，選取mr_T1數據為181×217×181，切片厚度為1 mm，含有3%噪聲的圖像序列繞垂直軸旋轉10°作為浮動圖像；選取質子密度加權成像mr_PD數據為181×217×181，切片厚度為1 mm，含有3%噪聲的圖像序列作為參考圖像。配準后的結果如
圖2
三維多模態配準
(a)浮動圖像;(b)參考圖像;(c)配準后圖像;(d)配準前后第90張切片圖像差值
Figure2. 3D image registration results(a) moving image; (b) static image; (c) registration result; (d) image difference of the 90th slice before and after registration
從圖 2中可知，(a)為T1數據181×217×181繞垂直軸旋轉10°并偏移5x5y，(b)為質子密度加權成像數據181×217×181，(c)為配準后T1數據，從直觀上可以看到矢狀面和冠狀面與參考圖像基本一致，水平面上的偏移量和旋轉角度都得到了較好的配準，本文算法對多模態三維圖像配準也能夠得到良好的效果。目前還沒有衡量圖像三維配準的金標準，故采用配準前后與參考圖像的差值比較來表示配準結果的好壞。如圖 2(d)所示，(d)中左邊是配準前圖像差值，配準前有明顯的偏移和誤差；右邊是配準后圖像差值，偏移和誤差有明顯的改善。本文算法對不同模態的三維醫學圖像也能取得較好的配準結果。
實驗三：選取mr_T1數據為181×217×181，切片厚度為1 mm，含有3%噪聲的圖像序列繞垂直軸旋轉10°作為浮動圖像；選取mr_T1數據為181×217×60，切片厚度為3 mm，含有3%噪聲的圖像序列作為參考圖像。配準結果如圖 3所示。
圖3
不同切片厚度的三維配準
(a)浮動圖像;(b)參考圖像;(c)配準后圖像;(d)配準前后第90張切片圖像差值
Figure3. 3D registration with different slice thicknessesl(a) moving image; (b) static image; (c) registration result; (d) image difference of the 90th slice before and after registration
如圖 3所示，(a)為T1數據181×217×181繞垂直軸旋轉10°并偏移5x5y，(b)為T1數據181×217×60，切片厚度為3 mm，由于切片厚度不同，在配準前先進行部分體積差值，使圖像尺寸相同再進行配準。(c)為配準后的數據，從直觀上可以看到冠狀面和矢狀面的灰度值有所變化，這是由部分體積差值引起的，水平面上選擇和偏移都得到了良好的改善，配準前后的差異如(d)所示，左邊是配準前圖像差值，右邊是配準后圖像差值。從(d)中可以看到，配準前兩幅圖像差異非常明顯，有較多的偽影和偏移；配準后圖像得到了較大的改善，偽影和偏移明顯減少。本文算法在處理不同尺寸的體素配準時也能取得較好的結果。
4 結論
本文從多種能量守恒相結合的思想出發，把局部結構張量守恒假設和灰度守恒相結合，得到新的能量函數，使算法能適應不同模態圖像的配準。將結構張量應用到能量函數中，使圖像的結構得以保持，以適應大形變的圖像配準。在算法實現中，采用多尺度分層細化策略，使圖像運動信息更加準確。實驗表明，本文算法對大形變配準、多模態配準都能得到較好的結果，并能得到比Demons算法精確度更好的結果。
