1醫學圖像配准方法
醫學影像學的發展為臨床診斷和治療提供了有效的輔助手段。目前,醫學成像模式可分解剖成像和功能成像,前者主要描述人體形態信息,後者主要描述人體代謝信息。基於各種原因,臨床上通常要將同一病人的多種模式成像結果結合起來進行分析,以提高醫學診斷和治療的水平。這就需要對不同模態的圖像進行融合,而融合首先要解決這幾幅圖像的嚴格對齊問題,即配准問題。醫學圖像配准,即通過尋找一種空間變換,使兩幅圖像的對應點達到空間位置和解剖位置的完全一致,配準的結果應使兩幅圖像上所有解剖點、或至少是所有具有診斷意義上的點都達到匹配。醫學圖像配准方法基本上可分為前瞻性和回溯性兩種[1]。前瞻性配準是基於外部特徵的圖像配准(有框架),它是利用外部基準點特徵,精度較高,但由於它是有創的,且操作較複雜,應用較少。回溯性配准即是基於內部圖像特徵的圖像配准(無框架),它是利用圖像內部本身的特徵,具有無創性和可回溯性,已成為配准算法研究的重點。回溯性配准技術可分為形狀特徵點的方法和像素(體素)相似性的方法。基於形狀的方法根據圖像間的共有特徵進行幾何配准,該方法原理簡單,應用面廣,但是需要較多的人工介入,且配准精度受特徵點(面)提取精度的限制;基於像素相似性的方法是基於圖像中所有體素的配准方法,不需要對圖像做預處理。近年來,引入了資訊理論中的一些概念,如聯合熵、相對熵和互信息[2],精度較高,可以達到亞像素級,已成為醫學圖像配准很重要的方法。互信息是資訊理論中的一個基本概念,是兩個隨機變量統計相關性的測度,當兩幅圖像基於共同解剖結構的圖像達到最佳配準時,它們對應像素的灰度互信息應達到最大。由於該測度不需要對不同成像模式下圖像灰度間的關係作任何假設,也不需要對圖像進行分割或任何預處理,所以被廣泛地用於多種醫學圖像配准,特別是當其中一個圖像的數據部分缺損時也能得到很好的配准效果。
2互信息基本概念
將兩幅待配準的圖像RF圖像和FA圖像的灰度值看作兩個隨機變量A和B,其灰度值為0～255,其機率密度函數分別為和,兩者的聯合機率密度函數為,則隨機變量A和B的互信息可以表示為:I(A,B)=H(A)+H(B)-H(A,B)(1)其中H(A)、H(B)和H(A,B)為隨機變量A與B的個體熵和聯合熵,其定義為:H(A)=-Σapr(a)logpr(a)H(B)=-Σbpf(b)logpf(b)a,b∈[0,255]H(A,B)=-ΣaΣbprf(a,b)logprf(a,b)(2)根據樣本集2個隨機變量之間互信息計算的Dobrushin公式[3],可以推導出互信息的計算公式為:I(A,B)=Σa,bprf(a,b)logprf(a,b)pr(a)pf(b)(a,b∈[0,255])(3)其中:pr(a)和pf(b)也稱為邊緣機率密度,可由聯合機率密度求得,即pr(a)=Σbprf(a,b)pf(b)=Σaprf(a,b)(4)對於聯合機率密度的估計,採用直方圖的方法進行。設表示隨機變量A和B的二維聯合直方圖,則其聯合機率密度prf(a,b)的估計為:prf(a,b)=h(a,b)Σa,bh(a,b)(5)在醫學圖像配准中,雖然兩幅圖像來源於不同的成像設備或來自於不同的時間,但它們都是基於共同個體的同一位置的解剖信息,所以當兩幅圖像在空間位置配準時,其重疊部分所對應的像素對的灰度的互信息達到最大值,以此時的變換參數作為空間變換的參數,通過空間變換達到圖像配準的目的。
3歸一化互信息
儘管互信息測度成功地應用於醫學圖像配准中,由於兩幅圖像重疊部分的大小對互信息的度量有很大影響,重疊部分減小,參與統計互信息的像素個數減小導致互信息值減小,互信息與兩個圖像重疊部分多少成正比,其次Studholme等[4]指出誤配數量增加可能導致互信息值增大。因此互信息值達到最大並不能保證得到正確的配准結果。為了解決這個問題,使目標函數能更加準確反映互信息量和配准參數之間的關係,Studholme等提出了一個歸一化互信息測度,歸一化互信息使配准函數更平滑,它能減少對圖像重疊部分敏感性,配准精度更高。NMI(A,B)=H(A)+H(B)H(A,B)(6)Maes等[1]利用熵相關係數(ECC)作為另一種歸一化相關係數。NMI和ECC有如下關係表示式:ECC=2-2/NMI(7)近年來歸一化互信息被許多學者廣泛應用,實踐證明,在剛性配准中,歸一化互信息比傳統互信息有更強的穩健型[5～7]。
4利用互信息配准常用方法
互信息的計算是最為關鍵的問題,計算結果對配准精度有很大的影響[8]。主要體現在插值運算、優化算法及提高配准速度等3個方面。Zhu和Cochoff[9]研究了採用不同的優化算法、插值方法、直方圖及多尺度逼近方法對配准結果的影響。目前,雖然提出了許多種方法,但還沒有一種非常有效的方法。
4.1空間變換
當兩幅圖像配準時,需要根據變換參數求出空間變換,浮動圖像的空間變換按照變換性質不同可分為剛體變換(rigid)、仿射變換(affine)、投影變換(projective)和非線性(curved)變換。剛體變換是指物體內部任意兩點間的距離保持不變;仿射變換是將直線映射為直線,並保持平行性;投影變換將直線映射為直線,但不再保持平行性質;非線性變換也稱做彎曲變換,它把直線變換為曲線。基於互信息的配准技術一般採用剛性變換或仿射變換。
4.2插值運算
在一般情況下,圖像經過幾何變換後,像素的坐標不會和原來的採樣網格完全重合,像素的灰度值也需要重新計算,這就需要對變換後的圖像進行重採樣和插值處理。常用的插值算法有最近鄰域法(Nearestneighboringinterpolation)、線性插值法(Trilinearinterpolation)和三線性部分體積分布(Trilinearpartialvolumedistribution)插值算法,簡稱PV插值方法。最近鄰域法具有計算量小、速度快的優點,但是存在質量不高的缺點。線性插值效果較好,運算量也不很大,故經常採用。三線性PV插值算法不是通過鄰居點確定所求像素的灰度,它是按照周圍8個像素和所求像素點的空間距離來分配權重,避免了一次插值運算,使互信息的計算更為精確,而且對於小的空間變換,增量,互信息的變化會更平滑,同時優化過程中的局部極值問題也會有所緩解[2]。文獻[10～12]討論了包括最近鄰居法、線性插值法、PV插值法在內的幾種插值算法,比較了各自的優缺點。文獻[10]詳細分析了線性插值法和PV插值法產生局部極值的原因,其實驗結果表明,非網格點上的PV插值操作造成的直方圖分散甚至比在位移更大的整數平移點上由於失配造成的直方圖分散還要嚴重,在整數平移點上產生局部極值。針對這個問題,文獻[11]提出了先驗聯合機率法和隨機重採樣法兩種改進的插值算法,先驗聯合機率法既保證了最大互信息量方法的有效性,又引入了與變換無關的先驗分布,增加了聯合分布的穩定性,使得互信息量隨配准參數變化更加平滑;隨機重採樣法在每個網格點施加一個輕微的擾動避免大量網格重合,這兩種方法在一定程度下抑制了局部極值的產生。文獻[12]提出帶擾動採樣的最近鄰居法(Nearestneighborwithjitteredsampling),網格點擾動後該坐標位置的灰度值由最近鄰居法確定,達到抑制局部極值的目的。
4.3優化策略
圖像配准在本質上是一個多參數優化問題[13],即尋找互信息達到最大時的幾個空間變換參數值。因此,配准問題實質是配准函數優化問題,但是,配准函數經常不是光滑的,存在許多局部極大值,給求解帶來很大難度。產生局部極大值主要有兩個原因[14]:一是兩幅圖像本身存在較好的局部匹配;二是在運算過程中產生的,如插值運算、圖像重疊部分的改變都有可能產生局部極值。避免局部極值的常用方法有:採用PV插值方法、圖像先濾波減少噪聲以及增大灰度直方圖窗口尺寸等方法。由於局部極大值的存在,優化算法的選取對配准結果有較大的影響,尤其對初始變換的魯棒性有很大影響。另一個關鍵問題是優化過程中參數變量的合適取值範圍[15]。優化過程中得到的值可能不是搜索空間中全局極大值,而是部分搜索空間中局部極大值。初始值偏離搜索區間太大很難使圖像配准正確。實踐證明,遺傳算法等隨機優化算法不適合求解此類問題,因為此類算法很容易跳出搜索區間。目前常用的方法是Powell優化算法[13],該算法輪流對變換參數進行優化,由於無需計算梯度,因而可以加快搜索最大互信息的速度,在每一維內使用Brent算法疊代搜索和估計配准參數,從而使互信息不斷增加。實踐表明Powell優化算法很容易受到局部極值的干擾[14]。另一個常用的方法是單純形算法,該方法也不要求計算梯度,與Powell算法只考慮單一變量相反,它同時考慮所有變量[9,13,16,17],但是該方法收斂速度不確定。Plattard等[18]使用Powell和單純形混合算法。Kagadis等[19]採用Powell和遺傳算法混合方法。Jenkinson和Smith用多解析度技術擴展了Powell算法。儘管爬山法是一種最簡單的優化方法,但是它在多解析度策略中有較好的效果,隨著圖像解析度的增加,爬山法搜索步長逐漸減小[20]。其他常規的算法如梯度上升法、Newton法、模擬退火法等也得到應用。為了提高配准函數全局最優值的搜索能力,Chen等[21]先將整幅圖像分成4個子塊,分別計算整幅圖像和4個子塊的互信息,並假定當整幅幅像互信息達到最大時,子塊圖像也應達到最大。Thevenaz等[22]同時利用3個不同窗口尺寸的直方圖計算的互信息來求解變換參數。採用多解析度方法是一種克服局部極值的一種有效方法[23]。
4.4優化速度加速策略
多解析度策略除了能夠避免出現局部極大值外,還能提高算法執行速度。配准按照從粗到精的方式執行,在低解析度下利用較少的時間進行初配,在高解析度下精配需要花銷較多的時間。Maes等[13]對多種優化算法在不同的解析度情況下結果進行了詳細比較。一般情況下,多解析度策略分為增加模型複雜度方法和增加數據複雜度方法[24],多解析度策略經常採用高斯金字塔方式。多解析度策略提高了運行速度、增加了魯棒性,不影響配准正確性。對於高解析度圖像採用多解析度策略是降低運算時間一種最有效的方式。最近提出的主要加速策略有:基於輪廓特徵點的互信息配准方法[25],先用小波變換或其它邊緣檢測運算元求出兩幅幅圖像的輪廓信息,利用聚類分析法求出輪廓特徵點,再求出特徵點對的互信息;Josien等[26]提出聯合互信息和梯度信息作為圖像配準的新標準,該算法不僅利用了圖像的灰度信息,而且還利用了圖像的空間信息(梯度信息),實驗表明,該方法比採用傳統互信息或歸一化互信息結果更準確,還能避免產生不正確的全局極大值以及由於插值產生的局部極大值。
5結束語
近幾年來,將資訊理論的互信息概念引入到醫學圖像配准中取得了很大進展。由於該方法不需要對圖像進行預處理、其精度高、適應性強、自動化程度高,能使多模態圖像的配准精度達到亞像素精度,因而在醫學圖像配准中得到廣泛應用。然而,互信息只是一個相似性測度,如何實現不同模態圖像間互信息的最大化到目前為止還沒有一套完全成熟的方法,選取什麼樣的優化策略、如何避免局部極值、如何提高優化速度、如何將互信息概念與其他概念結合起來形成一個更穩健的測度、怎樣將互信息應用到非剛性變換中等等問題均需要我們繼續深入研究。
醫學影像學的發展為臨床診斷和治療提供了有效的輔助手段。目前,醫學成像模式可分解剖成像和功能成像,前者主要描述人體形態信息,後者主要描述人體代謝信息。基於各種原因,臨床上通常要將同一病人的多種模式成像結果結合起來進行分析,以提高醫學診斷和治療的水平。這就需要對不同模態的圖像進行融合,而融合首先要解決這幾幅圖像的嚴格對齊問題,即配准問題。醫學圖像配准,即通過尋找一種空間變換,使兩幅圖像的對應點達到空間位置和解剖位置的完全一致,配準的結果應使兩幅圖像上所有解剖點、或至少是所有具有診斷意義上的點都達到匹配。醫學圖像配准方法基本上可分為前瞻性和回溯性兩種[1]。前瞻性配準是基於外部特徵的圖像配准(有框架),它是利用外部基準點特徵,精度較高,但由於它是有創的,且操作較複雜,應用較少。回溯性配准即是基於內部圖像特徵的圖像配准(無框架),它是利用圖像內部本身的特徵,具有無創性和可回溯性,已成為配准算法研究的重點。回溯性配准技術可分為形狀特徵點的方法和像素(體素)相似性的方法。基於形狀的方法根據圖像間的共有特徵進行幾何配准,該方法原理簡單,應用面廣,但是需要較多的人工介入,且配准精度受特徵點(面)提取精度的限制;基於像素相似性的方法是基於圖像中所有體素的配准方法,不需要對圖像做預處理。近年來,引入了資訊理論中的一些概念,如聯合熵、相對熵和互信息[2],精度較高,可以達到亞像素級,已成為醫學圖像配准很重要的方法。互信息是資訊理論中的一個基本概念,是兩個隨機變量統計相關性的測度,當兩幅圖像基於共同解剖結構的圖像達到最佳配準時,它們對應像素的灰度互信息應達到最大。由於該測度不需要對不同成像模式下圖像灰度間的關係作任何假設,也不需要對圖像進行分割或任何預處理,所以被廣泛地用於多種醫學圖像配准,特別是當其中一個圖像的數據部分缺損時也能得到很好的配准效果。
2互信息基本概念
將兩幅待配準的圖像RF圖像和FA圖像的灰度值看作兩個隨機變量A和B,其灰度值為0～255,其機率密度函數分別為和,兩者的聯合機率密度函數為,則隨機變量A和B的互信息可以表示為:I(A,B)=H(A)+H(B)-H(A,B)(1)其中H(A)、H(B)和H(A,B)為隨機變量A與B的個體熵和聯合熵,其定義為:H(A)=-Σapr(a)logpr(a)H(B)=-Σbpf(b)logpf(b)a,b∈[0,255]H(A,B)=-ΣaΣbprf(a,b)logprf(a,b)(2)根據樣本集2個隨機變量之間互信息計算的Dobrushin公式[3],可以推導出互信息的計算公式為:I(A,B)=Σa,bprf(a,b)logprf(a,b)pr(a)pf(b)(a,b∈[0,255])(3)其中:pr(a)和pf(b)也稱為邊緣機率密度,可由聯合機率密度求得,即pr(a)=Σbprf(a,b)pf(b)=Σaprf(a,b)(4)對於聯合機率密度的估計,採用直方圖的方法進行。設表示隨機變量A和B的二維聯合直方圖,則其聯合機率密度prf(a,b)的估計為:prf(a,b)=h(a,b)Σa,bh(a,b)(5)在醫學圖像配准中,雖然兩幅圖像來源於不同的成像設備或來自於不同的時間,但它們都是基於共同個體的同一位置的解剖信息,所以當兩幅圖像在空間位置配準時,其重疊部分所對應的像素對的灰度的互信息達到最大值,以此時的變換參數作為空間變換的參數,通過空間變換達到圖像配準的目的。
3歸一化互信息
儘管互信息測度成功地應用於醫學圖像配准中,由於兩幅圖像重疊部分的大小對互信息的度量有很大影響,重疊部分減小,參與統計互信息的像素個數減小導致互信息值減小,互信息與兩個圖像重疊部分多少成正比,其次Studholme等[4]指出誤配數量增加可能導致互信息值增大。因此互信息值達到最大並不能保證得到正確的配准結果。為了解決這個問題,使目標函數能更加準確反映互信息量和配准參數之間的關係,Studholme等提出了一個歸一化互信息測度,歸一化互信息使配准函數更平滑,它能減少對圖像重疊部分敏感性,配准精度更高。NMI(A,B)=H(A)+H(B)H(A,B)(6)Maes等[1]利用熵相關係數(ECC)作為另一種歸一化相關係數。NMI和ECC有如下關係表示式:ECC=2-2/NMI(7)近年來歸一化互信息被許多學者廣泛應用,實踐證明,在剛性配准中,歸一化互信息比傳統互信息有更強的穩健型[5～7]。
4利用互信息配准常用方法
互信息的計算是最為關鍵的問題,計算結果對配准精度有很大的影響[8]。主要體現在插值運算、優化算法及提高配准速度等3個方面。Zhu和Cochoff[9]研究了採用不同的優化算法、插值方法、直方圖及多尺度逼近方法對配准結果的影響。目前,雖然提出了許多種方法,但還沒有一種非常有效的方法。
4.1空間變換
當兩幅圖像配準時,需要根據變換參數求出空間變換,浮動圖像的空間變換按照變換性質不同可分為剛體變換(rigid)、仿射變換(affine)、投影變換(projective)和非線性(curved)變換。剛體變換是指物體內部任意兩點間的距離保持不變;仿射變換是將直線映射為直線,並保持平行性;投影變換將直線映射為直線,但不再保持平行性質;非線性變換也稱做彎曲變換,它把直線變換為曲線。基於互信息的配准技術一般採用剛性變換或仿射變換。
4.2插值運算
在一般情況下,圖像經過幾何變換後,像素的坐標不會和原來的採樣網格完全重合,像素的灰度值也需要重新計算,這就需要對變換後的圖像進行重採樣和插值處理。常用的插值算法有最近鄰域法(Nearestneighboringinterpolation)、線性插值法(Trilinearinterpolation)和三線性部分體積分布(Trilinearpartialvolumedistribution)插值算法,簡稱PV插值方法。最近鄰域法具有計算量小、速度快的優點,但是存在質量不高的缺點。線性插值效果較好,運算量也不很大,故經常採用。三線性PV插值算法不是通過鄰居點確定所求像素的灰度,它是按照周圍8個像素和所求像素點的空間距離來分配權重,避免了一次插值運算,使互信息的計算更為精確,而且對於小的空間變換,增量,互信息的變化會更平滑,同時優化過程中的局部極值問題也會有所緩解[2]。文獻[10～12]討論了包括最近鄰居法、線性插值法、PV插值法在內的幾種插值算法,比較了各自的優缺點。文獻[10]詳細分析了線性插值法和PV插值法產生局部極值的原因,其實驗結果表明,非網格點上的PV插值操作造成的直方圖分散甚至比在位移更大的整數平移點上由於失配造成的直方圖分散還要嚴重,在整數平移點上產生局部極值。針對這個問題,文獻[11]提出了先驗聯合機率法和隨機重採樣法兩種改進的插值算法,先驗聯合機率法既保證了最大互信息量方法的有效性,又引入了與變換無關的先驗分布,增加了聯合分布的穩定性,使得互信息量隨配准參數變化更加平滑;隨機重採樣法在每個網格點施加一個輕微的擾動避免大量網格重合,這兩種方法在一定程度下抑制了局部極值的產生。文獻[12]提出帶擾動採樣的最近鄰居法(Nearestneighborwithjitteredsampling),網格點擾動後該坐標位置的灰度值由最近鄰居法確定,達到抑制局部極值的目的。
4.3優化策略
圖像配准在本質上是一個多參數優化問題[13],即尋找互信息達到最大時的幾個空間變換參數值。因此,配准問題實質是配准函數優化問題,但是,配准函數經常不是光滑的,存在許多局部極大值,給求解帶來很大難度。產生局部極大值主要有兩個原因[14]:一是兩幅圖像本身存在較好的局部匹配;二是在運算過程中產生的,如插值運算、圖像重疊部分的改變都有可能產生局部極值。避免局部極值的常用方法有:採用PV插值方法、圖像先濾波減少噪聲以及增大灰度直方圖窗口尺寸等方法。由於局部極大值的存在,優化算法的選取對配准結果有較大的影響,尤其對初始變換的魯棒性有很大影響。另一個關鍵問題是優化過程中參數變量的合適取值範圍[15]。優化過程中得到的值可能不是搜索空間中全局極大值,而是部分搜索空間中局部極大值。初始值偏離搜索區間太大很難使圖像配准正確。實踐證明,遺傳算法等隨機優化算法不適合求解此類問題,因為此類算法很容易跳出搜索區間。目前常用的方法是Powell優化算法[13],該算法輪流對變換參數進行優化,由於無需計算梯度,因而可以加快搜索最大互信息的速度,在每一維內使用Brent算法疊代搜索和估計配准參數,從而使互信息不斷增加。實踐表明Powell優化算法很容易受到局部極值的干擾[14]。另一個常用的方法是單純形算法,該方法也不要求計算梯度,與Powell算法只考慮單一變量相反,它同時考慮所有變量[9,13,16,17],但是該方法收斂速度不確定。Plattard等[18]使用Powell和單純形混合算法。Kagadis等[19]採用Powell和遺傳算法混合方法。Jenkinson和Smith用多解析度技術擴展了Powell算法。儘管爬山法是一種最簡單的優化方法,但是它在多解析度策略中有較好的效果,隨著圖像解析度的增加,爬山法搜索步長逐漸減小[20]。其他常規的算法如梯度上升法、Newton法、模擬退火法等也得到應用。為了提高配准函數全局最優值的搜索能力,Chen等[21]先將整幅圖像分成4個子塊,分別計算整幅圖像和4個子塊的互信息,並假定當整幅幅像互信息達到最大時,子塊圖像也應達到最大。Thevenaz等[22]同時利用3個不同窗口尺寸的直方圖計算的互信息來求解變換參數。採用多解析度方法是一種克服局部極值的一種有效方法[23]。
4.4優化速度加速策略
多解析度策略除了能夠避免出現局部極大值外,還能提高算法執行速度。配准按照從粗到精的方式執行,在低解析度下利用較少的時間進行初配,在高解析度下精配需要花銷較多的時間。Maes等[13]對多種優化算法在不同的解析度情況下結果進行了詳細比較。一般情況下,多解析度策略分為增加模型複雜度方法和增加數據複雜度方法[24],多解析度策略經常採用高斯金字塔方式。多解析度策略提高了運行速度、增加了魯棒性,不影響配准正確性。對於高解析度圖像採用多解析度策略是降低運算時間一種最有效的方式。最近提出的主要加速策略有:基於輪廓特徵點的互信息配准方法[25],先用小波變換或其它邊緣檢測運算元求出兩幅幅圖像的輪廓信息,利用聚類分析法求出輪廓特徵點,再求出特徵點對的互信息;Josien等[26]提出聯合互信息和梯度信息作為圖像配準的新標準,該算法不僅利用了圖像的灰度信息,而且還利用了圖像的空間信息(梯度信息),實驗表明,該方法比採用傳統互信息或歸一化互信息結果更準確,還能避免產生不正確的全局極大值以及由於插值產生的局部極大值。
5結束語
近幾年來,將資訊理論的互信息概念引入到醫學圖像配准中取得了很大進展。由於該方法不需要對圖像進行預處理、其精度高、適應性強、自動化程度高,能使多模態圖像的配准精度達到亞像素精度,因而在醫學圖像配准中得到廣泛應用。然而,互信息只是一個相似性測度,如何實現不同模態圖像間互信息的最大化到目前為止還沒有一套完全成熟的方法,選取什麼樣的優化策略、如何避免局部極值、如何提高優化速度、如何將互信息概念與其他概念結合起來形成一個更穩健的測度、怎樣將互信息應用到非剛性變換中等等問題均需要我們繼續深入研究。
收藏