基于SPM的功能磁共振成像圖像頭動校正配準方法的分析與改進

陳 文 CHEN Wen 郁 蕓 YU Yun 劉宏毅 LIU Hongyi 胡新華 HU Xinhua 楊 坤 YANG Kun 趙 宇 ZHAO Yu

基于SPM的功能磁共振成像圖像頭動校正配準方法的分析與改進

陳 文 CHEN Wen 郁 蕓 YU Yun 劉宏毅 LIU Hongyi 胡新華 HU Xinhua 楊 坤 YANG Kun 趙 宇 ZHAO Yu

作者單位 南京醫(yī)科大學 江蘇 南京 210029

功能磁共振成像作為一種新興的神經(jīng)影像學成像技術，其重要性受到神經(jīng)、認知和心理學等領域的廣泛關注。本文針對功能磁共振成像數(shù)據(jù)預處理過程中的頭動校正環(huán)節(jié)，較為全面地闡述了國際知名腦功能成像軟件包SPM的配準原理，并在此基礎上提出了一種改進的配準方法，通過由粗到細的處理有效地提高了配準的準確度。此外，將差值平方和與歸一化互信息的比值作為配準測度的思想彌補了單用一種測度的局限性，增加了配準過程的魯棒性和可靠性。研究結果顯示該方法具有一定的實際意義和應用價值。

磁共振成像；圖像配準；圖像處理，計算機輔助；SPM；頭動校正

功能磁共振成像（functional magnetic resonance imaging，fMRI）是一種新的無侵入、無創(chuàng)傷的腦功能成像技術，為研究腦的認知活動、生理功能及病理狀態(tài)提供了有利的工具，已成為腦科學和生命科學研究的重要工具。

在fMRI技術中，影像圖像的后期處理具有重要意義，圖像的配準（registration）是其中的一個關鍵環(huán)節(jié)。在掃描采樣過程中，受試者不由自主的輕微頭動仍然不可避免，包括物理頭部運動（輕微的左右擺動或點頭等）和生理頭部運動（心跳、呼吸等），導致偽影出現(xiàn)［1］。為了提高腦功能圖像統(tǒng)計分析的準確性，需對圖像進行配準處理，以確保不同時間序列圖像空間位置的完全對應。

目前用于fMRI圖像處理的軟件包有AFNI、SPM、FSL、MRIcro等，其中SPM是一個專為腦功能成像數(shù)據(jù)分析而設計的通用軟件包，是國際上應用最普遍、最權威的腦功能成像處理分析軟件。本文基于SPM分析了fMRI圖像頭動校正的配準方法，并對其進行優(yōu)化和改進。

1 SPM圖像頭動校正配準的算法原理與結果分析

1.1 概述 為了適應圖像數(shù)據(jù)分析的需要，fMRI數(shù)據(jù)必須經(jīng)過多種圖像預處理過程，包括頭動校正、圖像融合、圖像分割、標準化、空間平滑等［1-3］。SPM中的頭動校正功能可以確定時間序列圖像與其參考圖像配準的參數(shù)，并在此基礎上將一個實驗序列中的每一幀圖像均與這個序列的參考圖像進行對齊，以矯正頭動［4］。本文討論的配準對應于SPM的Realign（estimate）過程。

1.2 配準方法 作為一種理想的單體（subject）單模態(tài)（modality）配準，頭部運動的校正常采用剛體運動模型。SPM配準算法的基本思想是通過迭代計算平移、旋轉參數(shù)，使參考圖像（通常為時間序列的第1幅）與后續(xù)序列圖像的不匹配度最小化，實現(xiàn)所有序列的配準［1］。

1.2.1 圖像配準原理 以2幅圖像為例，頭動校正是在相似性測度下尋找一個空間變換關系，使得經(jīng)過該變換后2幅圖像的相似性達到最大。例如，對參考圖像g和資源圖像f進行配準，則需要尋找一個變換Mr，使得：

其中，S是相似性測度，Mr為所求的最佳空間變換［5］。

因為需保證各向同性，需將體素坐標轉換為笛卡爾坐標。笛卡爾坐標系下的參考圖像ge和資源圖像fe為：

其中，Mg、Mf分別代表將參考圖像g和資源圖像f從體素坐標變換到笛卡爾坐標的轉換矩陣。

由公式（1）及公式（2）可得，變換過程為：

1.2.2 剛體變換（rigid transformation） 對于同一受試者，由于未產(chǎn)生結構間相對位置的改變，頭部運動即可分解為旋轉和平移兩部分的剛體變換。定義剛體變換向量=（q1，q2，q3，q4，q5，q6），其中q1，q2，q3依次代表頭部沿x，y，z 坐標軸的平移量大??；q4，q5，q6依次代表頭部繞x，y，z坐標軸的旋轉角度大小，則：

1.2.3 Gauss-Newton迭代算法 建立剛體變換模型之后，頭動校正的配準即可化歸為剛體變換向量= （q1，q2，q3，q4，q5，q6）的確定，在SPM中求解這6個獨立未知參數(shù)使用Gauss-Newton迭代算法來實現(xiàn)。

考慮將參考圖像g的坐標xi經(jīng)過剛體變換逆變換為資源圖像的坐標要使得兩者的相似性達到最大，SPM采用差值平方和作為測度，應使其最小化，定義目標函數(shù)為：

其中，xi表示第i個體素的坐標，為變換矩陣，f（）表示資源圖像的灰度值，g（）表示參考圖像的灰度值，u為2幅圖像的灰度平衡系數(shù)。

當差值平方和函數(shù)趨于平緩或達到指定的迭代次數(shù)Nth時迭代結束，此時的即為最佳的剛體變換參數(shù)。

1.3 實驗結果及分析 實驗所用為100張已做完時間校準的img圖像（a*.img），配準共用時27.2 s，所得結果為Graphics窗口中的2幅圖像，清晰地顯示了受試者在檢查時的平動和轉動情況（圖1）。

圖1 配準結果顯示的平動（A）和轉動（B）情況

經(jīng)分析發(fā)現(xiàn)，SPM的配準過程仍存在一些不足，有待進一步優(yōu)化：①SPM迭代可能出現(xiàn)過早收斂至局部最優(yōu)解的情況，具有一定的局限性［8］。伍亞軍等［9］提出了“安全區(qū)域”的概念，即配準目標函數(shù)最值的附近存在一個特殊區(qū)域，在此區(qū)域內(nèi)函數(shù)較為光滑，局部極值要緩和的多甚至消失，如果在迭代之前進行一次粗配準，使迭代的初始值進入安全區(qū)域，即可避開局部極值的影響。②SPM采用的配準測度是差值平方和最小化，然而這僅是一種相似性測度。目前基于像素灰度值的配準測度有很多，最常用的是圖像的差值平方和與互信息［10］。適當調(diào)整配準測度，將互信息最大化納入以彌補單用一種測度的局限性，將使配準的準度得到一定的提高。

2 SPM圖像頭動校正配準的算法改進與結果分析

2.1 配準方法的改進 從配準的精度和速度出發(fā)，本文在SPM的基礎上做了以下改進：為了避免收斂至局部極值，采用力矩主軸法在迭代之前進行一次粗配準；為了減少粗配準增加的時間，先提取圖像輪廓，省去不必要的計算負擔，相對加快了配準速度；為了減弱測度不全面性對配準精確度的影響，在原先配準測度的基礎上加入互信息測度，增加配準過程的魯棒性和可靠性。優(yōu)化后的算法流程見圖2。

圖2 優(yōu)化后的算法流程

2.1.1 輪廓提取 輪廓提取主要涉及形態(tài)學在邊緣提取、區(qū)域填充和移除對象方面的應用。首先對Canny算子提取的邊緣進行連接運算，然后進行區(qū)域填充、細化、腐蝕，用細化后圖像減去腐蝕過的圖像，移除干擾輪廓，即可得到去除內(nèi)部邊緣線的外輪廓圖［11］。

2.1.2 主軸矩法 力矩主軸法借用經(jīng)典力學中剛體質(zhì)量分布的概念，計算2幅圖像的質(zhì)心位置和主軸方向，進一步得到質(zhì)心和主軸對齊的平移、旋轉參數(shù)，從而限制位置偏差在一個較小范圍內(nèi)，達到粗略配準的目的。

首先，在笛卡爾坐標體系下，通過經(jīng)典力學方法計算質(zhì)心：

其中，x，y，z代表像素點的坐標，p（x，y，z）代表像素點的灰度值（視為像素的質(zhì)量）。

其次，根據(jù)“主軸是慣量矩陣的特征向量”，計算慣量矩陣：

其中，Ixx=μ020+μ002，Iyy=μ200+μ002，Izz=μ200+μ020，Ixy=Iyx=μ110，Ixz=Izx=μ101，Iyz=Izy=μ011。μpqr=∑x，y，z（x－xc）p（y－yc）q（z－zc）r，表示p+q+r階中心矩。

由于慣量矩陣I的標準化特征向量E等于旋轉矩陣R。

其中，α，β，γ分別為圖像繞x，y，z軸旋轉的角度。

最后求得：

由公式（8）、公式（11）可分別算出2幅圖像的xc，yc，zc，α，β，γ，即質(zhì)心和主軸方向得到確定。經(jīng)過作差，可粗略地得到平移和旋轉的配準參數(shù)［11-13］。

2.1.3 SPM迭代 完成粗配準后，將得到的配準參數(shù)作為SPM迭代的初始值，以達到避免誤差函數(shù)過早收斂的目的，同時可一定程度地減少迭代次數(shù)和縮短配準時間。

2.1.4 尋優(yōu) 經(jīng)過SPM迭代，在差值平方和最小化的原則下獲取了較為準確的配準參數(shù)，但這只是相似性測度的一種表現(xiàn)形式。為了提高準確度，彌補單用一種配準測度的局限性，考慮將互信息（mutual information）作為測度的一部分。

互信息具有精度高、適應性強、自動化程度高的優(yōu)點，在配準多模圖像和具有不同成像參數(shù)的同模圖像領域均有較好的效果［10，14］。在信息統(tǒng)計學中，互信息通常用于描述兩個系統(tǒng)之間的統(tǒng)計相關性，把2幅圖像的灰度值看作隨機變量A、B，其概率密度分布函數(shù)分別為pA（a）和pB（b），兩者的聯(lián)合概率密度分布函數(shù)為pAB（a，b），則隨機變量A和B的互信息可表示為：

其中，H（A），H（B）和H（A，B）為隨機變量A和B的個體熵和聯(lián)合熵，其定義為：

既往研究表明，當兩組基于共同解剖結構的圖像的空間位置達到一致時，由于圖像之間存在的相關性最大，灰度聯(lián)合概率密度分布p（a，b）最集中，聯(lián)合熵H（A，B）最小，互信息量I（A，B）最大?；バ畔y度對于2幅圖像的重疊覆蓋程度較為敏感，Studholme等提出了歸一化互信息（normalized mutual information）的概念［5］，并將其作為一種更為準確的配準測度，其定義為：

綜合兩種測度，本文選擇差值平方和與歸一化互信息的比值作為相關性測度：

當S值取得最小值時，則近似地認為相關程度達到最大。

該步驟在SPM所得配準參數(shù)的基礎上運用交替方向法進行進一步尋優(yōu)。

2.2 實驗結果分析 本文針對前面所述的100幅圖像進行頭動校正處理。根據(jù)實驗結果，100幅圖像中76%的圖像數(shù)據(jù)提示該方法較改進前具有較好的改善，其中17%顯著改善，其余24%的數(shù)據(jù)表示算法無改進。表1中為幾組典型的數(shù)據(jù)。

由表1可見，經(jīng)過配準后第9、26、43、65、87幅圖像的SSD減小、NMI增大、SSD/NMI減小，其配準的準確度得到較為全面的提高；第55幅和第73幅圖像盡管SSD有所增加，但NMI增大、SSD/NMI減小，即互信息水平得到提高，且互信息對配準效果的影響較大；第35幅圖像NMI有一定的增大，但意義不如前兩者明顯。

表1 改進前和改進后的相關典型數(shù)據(jù)

以第43幅圖像的第15層為例，改進前后的圖像見圖3。

圖3 改進前后配準結果比較

本文算法中，輪廓提取有效地省去了不必要的計算負擔，相對加快了粗配準的速度；基于主軸矩法的粗配準可限制位置偏差在一個較小范圍內(nèi)，避免迭代時收斂至局部極值，同時減少了迭代次數(shù)，一定程度地提高了運算速度；互信息測度的加入及其交替方向法判優(yōu)，彌補了單用一種測度的局限性，因而增加了配準過程的魯棒性和可靠性，但由于計算量較大，又延長了配準時間。

因此，此算法對SPM配準算法的改進有一定的意義，對fMRI圖像預處理中的配準操作具有一定的實用價值。然而，在配準時間上，此算法仍存在缺陷，有待進一步優(yōu)化。

3 結束語

本文較為全面地闡述了SPM頭動校正配準的算法原理，客觀地分析了SPM處理的優(yōu)、缺點，并在此基礎上提出了一種改進的配準方法，加入了輪廓提取、主軸矩法和基于SSD/NMI新測度的交替方向法判優(yōu)的步驟，由粗到精地進行處理，避免了局部極值問題，有效地提高了配準過程的準確度。后續(xù)研究將針對算法運算時間較長的問題進一步深入分析，探討并提出優(yōu)于SSD/NMI的配準測度。

［1］ 張志強，王世杰，盧光明.功能磁共振數(shù)據(jù)處理分析的原理及應用.中國醫(yī)學影像技術，2004，20（10）：1632-1635.

［2］ 趙健儀，尹建忠，紀東旭，等.嗅覺腦功能MRI對阿爾茨海默病的早期診斷.中國醫(yī)學影像學雜志，2013，21（5）：321-325.

［3］ 曾統(tǒng)軍，賴新生，黃泳，等.針刺外關穴和非穴fMRI腦功能成像的比較研究.中國醫(yī)學影像學雜志，2010，18（6）：532-537.

［4］ 劉永，張莉，馬竟鋒，等.基于fMRI的立體視覺數(shù)據(jù)分析輔助診斷系統(tǒng).中國醫(yī)療器械信息，2015，（7）：66-70.

［5］ 王麗，孫豐榮，王奕琨，等.基于互信息的顱腦MR影像序列的三維配準.計算機工程與應用，2011，47（31）：160-163.

［6］ 呂曉琪，李娜，張寶華，等.基于體素相似性的三維多模腦圖像配準研究.中國醫(yī)學影像學雜志，2013，21（2）：146-151.

［7］ 李京娜，王國宏，孫少燕，等.基于修正的結構相似度為測度的三維腦圖像配準.中國醫(yī)學影像學雜志，2013，21（8）：618-623.

［8］ 王超，杜劍平，閆鑌，等.SPM軟件包頭動校正原理的分析與評價.中國組織工程研究與臨床康復，2010，14（52）：9772-9776.

［9］ 伍亞軍，周正東，戴耀東.基于特征點互信息預配準的醫(yī)學圖像配準技術.南京航空航天大學學報，2007，39（5）：622-627.

［10］ 孫少燕，李斌，唐一源.fMRI時間序列圖像的配準測度及其相關優(yōu)化方法.應用基礎與工程科學學報，2004，12（3）：241-249.

［11］ 錢慧，胡志忠.基于粗配準和混合互信息的醫(yī)學圖像配準算法.生物醫(yī)學工程研究，2007，3（3）：266-270.

［12］ 吳鋒，錢宗才，杭洽時，等.基于輪廓的力矩主軸法在醫(yī)學圖像配準中的應用.第四軍醫(yī)大學學報，2001，22（6）：567-569.

［13］ 胡久鄉(xiāng)，丁國治.主軸法在三維圖像匹配中的應用.華中科技大學學報（自然科學版），2002，30（1）：38-40.

［14］ 杜俊俐，郭清宇.互信息醫(yī)學圖像配準的實時性研究.計算機工程與應用，2011，47（28）：162-165.

（本文編輯 張春輝）

Analysis and Improvement of Functional MRI Image Registration Method Based on SPM

As a new neuroimaging method，functional magnetic resonance （fMRI） with its importance of data processing has been widely recognized by neurology and cognitive psychology.Focusing on the realignment section in fMRI image preprocessing，this paper comprehensively describes the registration principle of SPM，an internationally-known software package which is specially designed for cerebral function imaging.An improved registration method is presented which effectively increases the accuracy.In addition，choosing the ratio of SSD and NMI as the registration measure can compensate for the limitation of using single measurement，which improves the robustness and reliability of the registration process.Experimental results prove the feasibility of this method.

Magnetic resonance imaging； Image registration； Image processing，computer-assisted； SPM； Realignment

10.3969/j.issn.1005-5185.2015.10.017

郁 蕓

Nanjing Medical University，Nanjing 210029，China

Address Correspondence to：YU Yun E-mail：yuyun@njmu.edu.cn

江蘇省大學生實踐創(chuàng)新訓練計劃項目（201410312057X）；南京市醫(yī)學科技發(fā)展項目（YKK12125）；南京市醫(yī)學科技發(fā)展項目（201303009）。

R445.2

2015-03-17

2015-04-20

中國醫(yī)學影像學雜志2015年 第23卷10期：780-784

Chinese Journal of Medical Imaging 2015 Volume 23（10）：780-784