基于重采樣改進的多圖譜分割算法

夏 瑞，馬 瑜，王文娜，羅宇卓，尚夢玉

(寧夏大學 物理與電子電氣工程學院，寧夏 銀川 750021)

0 引 言

阿爾茲海默癥[1](Alzheimer’s disease，AD)的醫(yī)學診斷表現(xiàn)為顳中回萎縮，海馬萎縮是阿爾茲海默癥早期診斷的標志。研究人腦核磁共振(MR)圖像中海馬體的分割，有助于醫(yī)生了解患者海馬體病變狀況，同時益于醫(yī)生結(jié)合分割結(jié)果給出的定位信息精確地執(zhí)行手術計劃。

由于人腦結(jié)構的復雜性以及腦組織結(jié)構圖像低對比度等原因，傳統(tǒng)分割方法不易得到精確分割結(jié)果。基于多圖譜配準的分割方法利用多個圖譜的先驗知識，降低單一圖譜配準產(chǎn)生的不確定性以及誤差，得到的分割精確度更高。三維醫(yī)學圖像配準分割技術，得到的分割結(jié)果能為診斷提供豐富的信息，缺點是需要較長的計算時間。因此，本文提出基于重采樣改進的多圖譜配準方法，兩次提取感興趣區(qū)域(region of interest，ROI)，并改進標簽融合算法，有效提高配準時間效率，同時保持較高配準精度。

1 基于多圖譜配準的海馬體分割方法

1.1 圖譜圖像重采樣原理

傳統(tǒng)的多圖譜配準算法大多采用“粗精”混合的配準方法，先采用全局剛性配準進行圖譜與目標圖像的對齊，縮小兩者之間的差異，再通過精細的非剛性配準得到精確的配準結(jié)果。直接用精細配準無法使得兩幅圖像很好地對齊，而非剛性配準過程要求圖譜圖像與目標圖像具有相同的空間采樣點和采樣點空間距離，可對齊兩幅圖像。

由于3D醫(yī)學圖像具有巨大的數(shù)據(jù)量，圖像之間的配準過程消耗計算時間過長，因此“粗精”混合的配準方法不能同時達到高精度、高時間效率的要求。基于以上原因，本文提出利用重采樣的方法代替“粗”配準階段，重采樣過程包括輸入圖像、變換和校對機，圖像的空間坐標通過變換進行映射以便生成一個新的圖像。

一幅圖像是使用離散網(wǎng)格對連續(xù)場的采樣，重采樣就是基于灰度的圖像配準的本質(zhì)。重采樣具有和“粗”配準一樣可以達到圖譜圖像與目標圖像對齊的目的，使圖譜圖像和目標圖像有相同的空間采樣點和采樣點空間距離，同時減少了“粗”配準環(huán)節(jié)所需時間。在重采樣過程中，當圖像從一個空間到另一個空間的映射需要在非網(wǎng)格位置進行操作時，校對機用于計算非網(wǎng)格位置的亮度值。變換用來映射點坐標，包括旋轉(zhuǎn)、縮放、平移等，本文采用的變換是仿射變換。重采樣從參考圖像獲得輸出信息，根據(jù)提供的參考圖像，使用參考圖像的間距、原點和方向?qū)Υ蓸訄D像進行映射。進行重采樣后，調(diào)整了浮動圖像和參考圖像采樣點不一致的問題，并且圖像的尺寸大小得到矯正。

重采樣操作使參考圖像與浮動圖像具有同樣的各向同性采樣率，圖譜圖像與目標圖像大小、中心達到一致，為后續(xù)進行非剛性配準做好準備。圖1是一個簡單重采樣的原理圖，重采樣涉及從原始圖像中提取像素的位置、插值灰度級，并將其重定位到校正圖像中的近似矩陣坐標位置。內(nèi)插的方法包括最近鄰，雙線性插值和立方卷積。如圖1所示，第一行是原圖像素空間采樣點、采樣空間大小示意圖，第二行是對圖像像素進行重采樣，改變原點坐標且采樣空間大小為原來的兩倍的圖示。

圖1 重采樣原理

1.2 微分同胚Demons配準算法

傳統(tǒng)Demons算法，空間變換的復合運算使用向量加法來近似，然而，圖像配準中用到的大多數(shù)空間變換形成的是李群，而不是向量空間。同時，傳統(tǒng)Demons算法不能到可逆變形場，易引起局部形變折疊。而微分同胚是可逆的光滑映射，微分同胚Demons算法[2-5]，能夠保證圖像的拓撲結(jié)構在配準前后保持不變，光滑且連續(xù)，并防止引入形變折疊，對大小形變都適用，在缺少可用的空間變換信息時，是很好的配準框架。

微分同胚Demons算法，利用公式c←s°exp(u)在李群中對幾何變換s進行更新。微分同胚Demons算法的目標能量函數(shù)為

(1)

其中，變形場u是一個稠密速度場

(2)

微分同胚Demons算法的實現(xiàn)過程：

步驟2 對u進行流體正則化，即u←Kfluid*u；

步驟3 利用牛頓方法計算李群得到的exp(u)，計算c←s°exp(u)；

步驟4 對s進行擴散正則化，s←Kdiff*s。

1.3 多分辨率配準

多分辨率配準[6]是在不同分辨率條件下用不同的比例進行配準，并且上一級的配準結(jié)果用作下一級配準的輸入，直到達到最好的比例范圍。多分辨率配準能夠提高配準的成功率，并且在粗糙比例時消除局部噪聲。多分辨率配準框架如圖2所示，參考圖像和浮動圖像被視為兩個圖像金字塔，金字塔根據(jù)用戶定義的收放因數(shù)來進行平滑和二次采樣圖像。

圖2 多分辨率配準處理概念

本文將多分辨率配準思想加入到微分同胚配準過程中，既提高了配準的精確度，又降低了配準的時間。

1.4 基于非局部圖像塊相似性的加權選擇算法

多圖譜配準得到的形變結(jié)果傳遞到圖譜相應的標記圖像，即可得到相應的分割結(jié)果，再將多個圖譜的分割結(jié)果進行融合得到一個綜合的分割結(jié)果。當前研究較為成熟的標簽融合算法有：加權選擇算法、STAPLE算法、COLLATE算法等。本文以K近鄰搜索算法對加權選擇算法進行改進，并在融合階段再次提取ROI，使融合算法運行更高效。

1.4.1 加權選擇標簽融合算法

簡單加權選擇算法[7,8](majority voting)考慮大多數(shù)分割結(jié)果在該采樣點是否存在，如果大多數(shù)分割結(jié)果在該點的像素值為1，那么判定該點為分割結(jié)果。該算法不需要先驗知識，對每個像素進行選擇，假設圖譜與目標圖像進行配準，并將形變向量傳遞給標記圖像后，圖譜表示為Ai,Li，對于空間中的每一個位置x，目標圖像的分割結(jié)果為

(3)

(4)

其中，i為圖譜的個數(shù)，Ltarget(x)表示在位置x處的目標圖像的標簽像素。

改進后的加權選擇算法，充分考慮到傳統(tǒng)融合算法中未用到圖譜灰度圖像中信息這個缺陷，得到如下的融合公式

(5)

(6)

其中，權重系數(shù)Wi是目標圖像與形變后圖譜灰度圖像的相似性，一般為歸一化互信息(NMI)、歸一化互相關系數(shù)(NCC)、平均誤差平方和(MSD)。融合算法中引入了圖譜灰度圖像的先驗信息，在一定程度上能夠減小誤差，得到更準確的融合結(jié)果。

為了更精確地應用目標圖像與圖譜灰度圖像的先驗信息，一些學者研究了基于局部塊相似性的加權融合算法，其思想是，在計算相似性時，將圖像中每個體素用以該體素為中心的一個圖像塊來代替，圖譜圖像塊與目標圖像塊之間的相似性越高，則其擁有相同標簽的概率越大

(7)

(8)

υ(x)是目標圖像在位置x處的標簽值，W(x,xs)是以x為中心的目標圖像塊與第s個圖譜中以x為中心的圖譜圖像塊的相似性，Lxs是第s個圖譜中在位置x處的標簽值，h是局部適應參數(shù)。

考慮到配準可能產(chǎn)生誤差，為了減弱這種誤差對融合結(jié)果的影響，Pierrick Coupe等[9]研究了基于非局部塊相似性的加權選擇算法，克服了局部塊加權的缺點，通過在對應塊的鄰域進行搜索，增加了圖像信息，來計算相似性

(9)

(10)

V是圖譜在位置x處的一個搜索鄰域塊，Lxs,j是第s個圖譜在搜素鄰域V中位置j處的標簽，W(x,xs,j)代表以x中心的目標圖像塊與以xs,j為中心的圖譜圖像塊的相似性。

對于目標圖像的同一個體素x，對應的圖譜圖像增加了塊的搜索范圍，則搜素鄰域內(nèi)圖像塊的個數(shù)變多，將相似度較低的塊的權重直接置為0，可以排除不相似塊對整體權重的影響。Pierrick Coupe等[9]在計算權重時，選擇如下相似性公式對圖譜塊進行預選

(11)

其中，μ和σ是在位置x處目標圖像塊的均值和方差，μs,j和σs,j是第s個圖譜在x的搜素鄰域中位置j處的圖譜圖像塊的均值和方差。s值越大，兩個塊越相似，當s小于預先設定的閾值時，舍棄與該s對應的圖譜塊，不進行權重計算。

1.4.2 基于K近鄰搜索的加權選擇融合算法

為了提高基于非局部塊加權融合算法的效率，本文進行了更深入的研究，通過引入K近鄰搜索算法，快速搜索圖像塊的k個近鄰，只用k個非局部塊來進行局部加權。

K近鄰搜索算法是一種有監(jiān)督的學習分類算法。對于一個樣本，尋找與它最相近的k個鄰居樣本。實現(xiàn)K近鄰搜索算法可以用kd-tree或VP-tree。

k-d樹(k-dimensional樹的縮寫)是一種對K維空間中的實例點進行存儲以方便對齊進行快速檢索的樹形數(shù)據(jù)結(jié)構，k-d樹可以用于多維空間關鍵數(shù)據(jù)的搜索，例如范圍搜索和最近鄰搜索。

最簡單的最近鄰搜索，也稱為窮舉搜索，依次計算樣本集中每個實例點到目標點的距離，然后取最小距離的那個點。當樣本集較大時，這種策略非常耗時。對于大規(guī)模的高維數(shù)據(jù)空間，kd-tree是一種有效的快速搜素方法。利用kd-tree搜索最相近的k個數(shù)據(jù)點，減少了搜索全部數(shù)據(jù)點的計算量。

因此，本文將權重公式改寫如下

(12)

xs,j∈K表示xs,j在x的搜索鄰域V中，并且是x的k個近鄰，這樣這個圖像塊才被進行權重計算。由于非局部塊加權融合算法增加了搜索范圍，導致圖像塊的數(shù)量增加，加入K近鄰搜素算法，只使用最近的k個非局部塊來進行局部加權，大大減少了計算量，提高了融合算法的效率。

下面是基于K近鄰搜素的非局部塊加權融合算法實現(xiàn)過程：

步驟1 標注所有配準形變之后的圖譜標記圖像的連通區(qū)域，獲得包含所有標注區(qū)域的一個區(qū)域塊

Vo=V1∪V2∪…Vn

步驟2 根據(jù)步驟1中計算得到的立方區(qū)域塊Vo的大小裁剪目標圖像T得到TROI，以及圖譜圖像Ai,Li其,中i=1,2,…,N為參與配準的圖譜數(shù)量，得到AROI,LROI；

步驟3 K近鄰搜素算法對圖譜塊進行預選，剔除與目標圖像塊相差較大的圖譜塊，減少更多的計算量；

步驟4 構建概率圖譜，獲得最終的分割結(jié)果。

1.5 分割結(jié)果的評價

圖像分割結(jié)果的好壞需要評價標準來進行衡量。主觀評價即觀察者主觀視覺觀察，易受環(huán)境等因素影響和觀察者的主觀認知影響。常用的客觀評價方法為相似性測度(Dice)，計算兩個二值分割圖像的重疊程度，式(13)來評價目標圖像的專家手動分割結(jié)果SF和浮動圖像分割結(jié)果SM重疊的部分

(13)

∩表示兩個集合的交集。相似性測度因子d取值范圍是0到1，0表示兩分割區(qū)域完全不重疊。1表示兩分割區(qū)域完全重疊，即配準分割結(jié)果與金標準一致，分割結(jié)果十分精確。

2 實驗結(jié)果分析

本文采用Visual Studio 2010軟件平臺結(jié)合ITK開源代碼庫進行MR人腦圖像海馬體的自動分割。數(shù)據(jù)來源于倫敦帝國理工學院醫(yī)學腦部研究數(shù)據(jù)庫[10]，共20個T1-MR灰度圖像和對應的標記圖像樣本，灰度圖像標號為a01到a20，對應的標記圖像標號為a01-seg到a20-seg。其中，標記圖像中標記了人腦中67個結(jié)構。

算法流程如圖3所示。

圖3 本文算法流程

圖譜的選取在一定程度上影響著分割結(jié)果的好壞，當與目標圖像差別較大時，配準結(jié)果越差，得到分割結(jié)果越偏離金標準。P.Aljabar等[11]在文章中描述了相似性測度因子和圖像的數(shù)目存在如下關系

(14)

a、b是常數(shù)，滿足0≤a≤1,b>0，n是圖譜的數(shù)目。圖譜數(shù)目的增多反而使分割結(jié)果與金標準的相似度下降。Sabuncu等[7]的研究表明，若隨機選擇圖譜，Majority Vo-ting的表現(xiàn)不是單調(diào)的，且圖譜數(shù)量超過10精確度開始下降。若選擇最優(yōu)圖譜，圖譜在15個左右時Majority Voting的結(jié)果變得平穩(wěn)?；谏厦鎺讞l結(jié)論，綜合本文的實驗，同時考慮標簽融合的效率問題，本文選取了與目標圖像互信息最高的10個樣本作為圖譜圖像進行實驗，對分割精度的影響較小。

2.1 圖像預處理

2.1.1 顱骨剔除

人腦MRI中背景和非腦組織占有較大比重，容易造成配準的誤差。剔除腦殼是重要的預處理工作，本文利用Stefan Bauer等[12]提出的算法剔除腦殼，克服了傳統(tǒng)BET算法運算時間長的缺點，并可以有效剔除影響配準精確度的腦殼與背景部分。

圖4 顱骨剔除結(jié)果

2.1.2 重采樣并提取待分割組織ROI

本文以灰度圖a04、標記圖a04-seg為參考圖像，對圖譜進行以參考圖像為基準的重采樣操作。參考圖像的尺寸大小為195×198×170，體素間距是0.937×0.937×0.937，圖像中心為(90.89, 92.29, 79.18)。例如，當圖譜a01的尺寸大小、體素間距、圖像中心分別為176×198×160、0.937×0.937×0.937、(81.99, 92.29, 74.49)，則重采樣后a01具有和參考圖像一樣的大小、間距、圖像中心。

對目標圖像以及圖譜圖像進行初次提取感興趣區(qū)域(ROI)。如圖5所示，為灰度圖譜及對應的標記圖像初次提取ROI的結(jié)果。

圖5 以左、右海馬體為中心提取ROI

2.2 配準及標記融合結(jié)果

本文以重采樣結(jié)果取代剛性配準環(huán)節(jié)，并以提取的ROI作為微分同胚Demons配準對象，大幅度縮短配準時間。微分同胚Demons配準產(chǎn)生的一個形變場的示意圖，如圖6所示。

圖6 待配準圖像映射到目標圖像的變形域

標簽融合的結(jié)果即分割結(jié)果。圖7為目標圖像右海馬體的金標準(專家手動分割結(jié)果)與本文方法分割結(jié)果比較圖：第一行至第三行分別為軸向、矢狀面、冠狀面海馬體切片圖。每行左圖為目標圖像金標準；右圖為分割結(jié)果與金標準的覆蓋比較圖，其中白色區(qū)域為兩者重合區(qū)域，灰色邊緣為兩者非重合區(qū)域。從圖中可以看到本文分割方法得到的海馬體覆蓋了大部分金標準結(jié)構，即分割結(jié)果較為精確。

圖7 海馬體切片比較

如圖8所示，第一行是目標圖像右海馬體的金標準；第二行是傳統(tǒng)方法分割結(jié)果(左圖)以及本文方法分割結(jié)果(右圖)的三維顯示效果?？梢钥吹皆诜指罴毠?jié)上，本文方法更接近于專家手工分割，傳統(tǒng)方法分割結(jié)果表面過于光滑，丟失了分割對象的部分表面紋理信息。

圖8 右海馬體分割結(jié)果三維顯示效果

如表1所列，是以a04為目標圖像的海馬體分割精度比較。改進的加權選擇算法比簡單加權選擇算法精度略有提高，而本文方法得到的左右海馬體的分割精度均高于簡單加權選擇算法、改進的加權選擇算法，Dice值有10%～20%的提高。

表1 分割結(jié)果與目標圖像的相似性測度

以右海馬體分割過程耗時來進行分割效率分析，表2、表3中Resamp表示重采樣、ROI表示提取感興趣區(qū)域、GReg表示剛性配準、DReg表示微分同胚Demons配準、Wrap表示形變映射、Fuse 表示標簽融合、Total是總時間。如表2所列，是使用本文方法分割海馬體時各操作過程所耗時間，表3所列，是使用傳統(tǒng)“粗精”配準方法分割海馬體所用時間。

表2 本文方法分割右海馬體時間/s

表3 傳統(tǒng)“粗精”配準法分割右海馬體時間/s

從表1、表2可以看到，分割以a04為目標圖像分割右海馬體時，本文方法重采樣耗時僅4 s～5 s,相比于傳統(tǒng)“粗精”混合算法中粗配準時間387 s～781 s，效率提高了95～155倍；并且本文方法提取ROI后，微分同胚Demons配準時間為2 s～3 s，而傳統(tǒng)方法中微分同胚Demons配準用時147 s～151 s，精配準的效率提高了49～72倍；在形變域映射過程，相較于傳統(tǒng)方法中7 s～8 s的用時，本文方法僅為0.5 s～1 s；最后，標簽融合階段，改進之后的加權選擇算法比傳統(tǒng)方法減少了40.1 s。整個多圖譜分割過程，本文方法所用時間共117.7 s，而傳統(tǒng)方法用時為6852.3 s，大大提高了分割的效率問題，實現(xiàn)了在兩分鐘之內(nèi)完成人腦MR圖像中海馬體的分割。

3 結(jié)束語

三維醫(yī)學圖像配準分割是醫(yī)學圖像研究的熱點，大多數(shù)研究者在分割精度上進行了深入的研究，并取得了較好的結(jié)果。處理三維醫(yī)學圖像耗時較長，人們在專注研究分割精度時往往忽略分割時間效率，甚至以犧牲時間效率來提高分割精度，很難適應于臨床應用對于快、精、準的要求。基于此，本文對基于多圖譜配準的海馬體分割算法進行研究，針對分割效率問題，改變傳統(tǒng)“粗精”混合配準方法，兩次提取人腦海馬體ROI，在重采樣基礎上進行微分同胚Demons配準，并改進標簽融合算法。實驗結(jié)果表明，本文方法在分割精度上比傳統(tǒng)方法提高了10%～20%左右，在分割時間效率上是傳統(tǒng)“粗精”配準分割算法速度的58～60倍，精度與效率均有提升。后續(xù)將考慮用實際數(shù)據(jù)來對本文算法進行驗證，繼續(xù)完善算法，以提高精度。