聶生東田曉迪錢黎俊陳增愛許建榮

1(上海理工大學(xué)醫(yī)學(xué)影像工程研究所，上海 200093)2(上海交通大學(xué)附屬仁濟醫(yī)院放射科，上海 200127)

基于相位展開的磁敏感加權(quán)圖像中相位偽影處理方法的研究

聶生東1*田曉迪1錢黎俊2陳增愛2許建榮2

1(上海理工大學(xué)醫(yī)學(xué)影像工程研究所，上海 200093)2(上海交通大學(xué)附屬仁濟醫(yī)院放射科，上海 200127)

針對Hanning濾波器難以有效解決磁敏感加權(quán)成像過程中產(chǎn)生的相位偽影問題，提出一種基于相位展開技術(shù)的相位偽影處理方法。采用兩個步驟有效消除磁敏感加權(quán)成像過程中產(chǎn)生的相位偽影。首先，采用傳統(tǒng)的Hanning濾波器對相位偽影進行初步處理;然后，采用相位展開技術(shù)，對由Hanning濾波器處理后的殘余偽影進行進一步處理，從而達到完全消除相位偽影的目的。對掃描層厚為0.6 mm的100層腦部SWI原始圖像進行相位偽影的處理，結(jié)果表明，該方法能有效消除磁敏感加權(quán)圖像中的相位偽影。采用該方法處理SWI圖像中的相位偽影，可為后續(xù)采用最小密度投影方法顯示清晰的血管影像奠定基礎(chǔ)。

磁敏感加權(quán)成像;相位偽影;Hanning濾波器;相位展開

Abstract:The phase artifact，produced during susceptibility weighted imaging(SWI)，can not be eliminated efficiently by a Hanning filter.To solve this kind of problem，a phase artifact processing method was presented in this paper based on phase unwrapping technology.The method includes two steps.Firstly，a traditional Hanning filter was adopted to process the phase artifact.Then，the remnant phase artifact was processed using phase unwrapping technique.The proposed method was evaluated using 100 thin-section brain SWI slices with the interval of 0.6 mm.The results showed that the phase artifact could be completely eliminated from SWI images.The application of our method could make a basis for displaying clear vein blood vessel image by using minimum intensity projection technique.

Key words:susceptibility weighted imaging;phase artifact;Hanning filter;phase unwrapping

引言

磁敏感加權(quán)成像(susceptibilityweighted imaging，SWI)是近年來發(fā)展起來的一種全新的磁共振成像方法［1-2］。不同于以往的質(zhì)子密度、T1或T2加權(quán)成像，這種方法可利用不同組織間磁敏感性的差異產(chǎn)生圖像對比，進而可對各組織顯影的新技術(shù)。與傳統(tǒng)磁共振成像方式相比，SWI在靜脈顯影方面具有獨特的優(yōu)勢，可應(yīng)用于腦腫瘤、腦出血或其他有靜脈參與的病灶研究中，從而有效地改善對這些疾病的診斷。靜脈血的主要成分為順磁性的去氧血紅蛋白，動脈血則是抗磁性的氧合血紅蛋白［3-5］。這種磁敏感性差異將最終導(dǎo)致兩種血管信號的相位信息強度不同，在與幅度信息相加權(quán)后，使靜脈能獨立于動脈清晰成像，這是SWI成像的理論基礎(chǔ)。在臨床醫(yī)學(xué)診斷中，由外傷所導(dǎo)致的微出血以及很多腫瘤病變都與靜脈有關(guān)，目前的成像手段對靜脈血管的成像效果都不甚理想，而利用SWI可獲得清晰的靜脈血管影像，因此SWI具有廣泛的臨床應(yīng)用價值。

圖2 相位圖像中的相位偽影及處理效果。(a)未經(jīng)處理的相位偽影;(b)經(jīng)Hanning濾波器處理后的殘留偽影;(c)偽影經(jīng)相位展開處理后的效果Fig.2 Phase artifact and the processed effect.(a)phase artifact without process;(b)the remnant phase artifact after processed by hanning filter;(c)processed result after phase unwrapping

由于SWI是一種新的磁共振成像技術(shù)，現(xiàn)有的磁共振成像儀得到的僅僅是SWI原始圖像，其所配備的圖像處理工作站沒有對SWI原始圖像進行處理的相關(guān)軟件，目前國內(nèi)外也沒有成熟的軟件實現(xiàn)對SWI原始圖像的后處理［5］。所謂的 SWI原始圖像是指采用梯度回波序列掃描獲得的幅值圖像(magnitude image)和相位圖像(phase image)，如圖1中的(a)和(b)所示。常規(guī)磁共振成像中所使用的都是幅值圖像，它描述弛豫過程中質(zhì)子發(fā)出的信號強度，相位圖像則描述質(zhì)子在該過程中行經(jīng)的角度。幅值圖像中包含了絕大部分的組織對比信息，而相位圖像則從磁敏感性的角度反映組織對比，特別是磁化率差異較大的組織。這兩種圖像是在掃描過程中同時獲得的，所以，它們總是成對出現(xiàn)，每一對圖像所對應(yīng)的解剖位置都完全一致。但是，從圖1中的(a)和(b)很難獲得血管的診斷信息。要利用原始SWI圖像獲得連續(xù)清晰的血管影像，需要對原始圖像進行一系列復(fù)雜的后處理過程。

在對SWI原始圖像的處理中，需要解決的關(guān)鍵問題就是如何有效處理SWI相位圖像中的相位偽影(或 field inhomogeneity artifact，磁場不均勻性偽影)。相位偽影是由于磁場的不均勻性造成的，通常出現(xiàn)在具有不同磁化率物質(zhì)的交界面周圍，如副鼻竇、蝶鞍和顱底等位置［5］。在成像過程中由于磁場的不均勻性，導(dǎo)致了在SWI的相位圖像中出現(xiàn)相位模糊(或相位纏繞)問題，在相位圖上，表現(xiàn)為環(huán)繞在交界面周圍的擾動波紋，如圖2(a)所示(圖中圈出的部分)。近年來，對 SWI原始圖像進行處理的文獻報道非常少，且基本上都在國外文獻上發(fā)表。目前，已報道的對SWI原始圖像處理的文獻中通常采用 4 個步驟［1，5］:

1)幅值和相位圖像的重組;

2)消除相位圖像中磁場不均一性偽影;

3)制作相位蒙片并與幅值圖像加權(quán)獲得磁敏感加權(quán)圖像;

4)通過最小密度投影顯示連續(xù)層面的靜脈血管結(jié)構(gòu)。

其中，第2步是一個關(guān)鍵步驟，其對偽影的處理結(jié)果直接影響到后續(xù)的采用最小密度投影方法顯示血管影像的效果。在目前報道的文獻中，都是采用Hanning濾波器來處理SWI相位圖像中的偽影問題。利用Hanning濾波器處理SWI相位圖像中的相位偽影，存在的主要問題是:經(jīng)Hanning濾波器處理后，相位偽影不能完全被消除，仍然殘留部分偽影，如圖2(b)所示。如果不能完全消除偽影，殘留偽影將會保留在后續(xù)生成的蒙片上(如圖3(a)所示)，將會對最終血管顯示結(jié)果造成不利影響。為此，本研究提出基于相位展開技術(shù)的SWI圖像中相位偽影的處理方法，首次將廣泛應(yīng)用在合成孔徑雷達干涉測量、合成孔徑聲納、光學(xué)干涉、微波干涉等領(lǐng)域的相位展開技術(shù)應(yīng)用于處理 SWI圖像中的偽影［6，7］，成功解決了 SWI圖像中的相位偽影問題，如圖2(c)所示。由相位展開消除偽影后生成的蒙片，殘留的偽影明顯被消除，如圖3(b)所示。

圖3 后續(xù)生成的相位蒙片。(a)由圖2(b)生成的相位蒙片;(b)由圖2(c)生成的相位蒙片F(xiàn)ig.3 The different phase mask produced by Fig.2(b)and Fig.2(c).(a)the phase mask produced by Fig.2(b);(b)the phase mask produced by Fig.2(c)

1 相位展開技術(shù)及其應(yīng)用

1.1 關(guān)于相位纏繞與相位展開問題

關(guān)于相位纏繞(wrapping)問題，可以用圖4進行簡單說明。一個相位本應(yīng)在0～4π的信號因為某種原因，落在了-π到+π之間，這就是所謂的相位纏繞問題。將此纏繞信號重新恢復(fù)到原相位范圍，稱為相位展開(phase unwrapping)，有些文獻稱為相位解纏。

圖4 真實相位與纏繞相位Fig.4 The sketch map of the true phase and the intertwist phase

在合成孔徑雷達干涉測量(interferometric synthetic aperture radar，InSAR)中存在相位纏繞問題［6-13］。合成孔徑雷達干涉測量技術(shù)，是以合成孔徑雷達復(fù)數(shù)據(jù)提取的相位信息為信息源，獲取地表的三維信息和變化信息［6］。InSAR通過兩副天線，同時觀測(單軌模式)或兩次近平行觀測(重復(fù)軌道模式)獲取地面同一景觀的兩幅包含幅值與相位信息的二維復(fù)數(shù)據(jù)圖像［6-7，13］。因此，與 SWI 成像一樣，InSAR一次測量也得到同一景觀的一對圖像，即幅值圖像和相位圖像。由于相位展開技術(shù)能有效解決合成孔徑雷達干涉測量中的相位纏繞問題，因此應(yīng)用相位展開技術(shù)解決 SWI圖像中的相位纏繞(相位偽影)問題也是一個自然的選擇。

1.2 納入標準 （1）年齡≥60歲。（2）本次住院行冠脈造影，且造影結(jié)果顯示冠脈血管1支及（或）以上血管病變≥50%。（3）患者冠心病診斷符合我國2010年ST段抬高型心肌梗死診療指南［9］中或我國2012版非ST段抬高型急性冠脈綜合征［10］的標準。（4）患者糖尿病診斷均符合我國2型糖尿病防治指南（2013年版）［11］中的相關(guān)標準。（5）臨床資料及隨訪資料完整。（5）對本研究知情同意，且愿意按時接受隨訪。

與非相干信號處理只關(guān)心信號的幅度與強度不同，相干信號處理還同時涉及信號的時間和空間的相位信息［8］。SWI成像和干涉孔徑雷達成像是相干信號處理的典型代表，然而因為系統(tǒng)原因，所得到的與感興趣的物理量相聯(lián)系的相位均被纏繞在(-π，π］之間，存在著相位模糊(相位纏繞)的現(xiàn)象?；謴?fù)與磁場一致的真實相位，就是本研究所要討論的相位展開問題。

1.2 相位展開技術(shù)簡述

從數(shù)學(xué)上講［6，13］，給定一個二維相位矩陣 Ψi，j，為了對該矩陣進行相位展開，需要對每個像素點(i，j)加上或減去2π的整數(shù)倍，從而得到一個連續(xù)的函數(shù):

在解決相位纏繞的問題中，目前有很多方法［6，13］，但其中主要采用的方法有路徑跟蹤法和最小范數(shù)法。路徑跟蹤法(path-following method)是建立在殘差定理的基礎(chǔ)上，通過建立枝切線，阻止積分路徑穿過相位不連續(xù)區(qū)域，從而實現(xiàn)與路徑無關(guān)的相位積分，得到?jīng)]有誤差傳遞的相位展開結(jié)果。其代表算法包括Goldstein的枝切法(branch cut)［9］、Xu 和 Cumming的區(qū)域生長法(region growing)等。最小范數(shù)法則基于最小均方意義，其思想是最小化測量相位梯度與展開相位梯度的均方差，主要工作集中在 Romero、Franceschetti及Lanari的格林函數(shù)法，Ghigilia、Romero的先驗條件共軛梯度法(PCG)［10］，Pritt于 1996 年提出的加權(quán)多重網(wǎng)格算法［11］。此外，通過求相位數(shù)據(jù)的二次偏導(dǎo)數(shù)并設(shè)置閾值，從而生成二值質(zhì)量圖，這就是1991年Bone提出的利用質(zhì)量評價來指導(dǎo)相位展開的質(zhì)量圖引導(dǎo)法［14］。

為了從眾多相位展開方法中選擇能夠有效處理SWI圖像中相位偽影的方法，分別采用Goldstein枝切法、質(zhì)量圖引導(dǎo)法和最小二乘法等不同算法，對SWI相位數(shù)據(jù)進行對比處理實驗，結(jié)果發(fā)現(xiàn)，以Goldstein的枝切法為代表的路徑積分法對SWI相位數(shù)據(jù)的處理效率最高［15］。因此，下面主要介紹Goldstein的枝切法在SWI相位偽影處理中的應(yīng)用。

1.3 基于路徑跟蹤的相位展開算法

Goldstein等提出的枝切法是基于路徑跟蹤的相位展開方法的典型代表［9］。其基本思想是:首先，通過識別正負殘差點，并連接鄰近的正負殘差點，實現(xiàn)殘差點的電荷平衡;然后，生成最優(yōu)枝切線來確定積分路徑，防止誤差沿積分路徑傳遞;最后，使相位展開的積分路徑不穿過枝切線，從而得到真實的相位信息。算法基本步驟包括識別殘差點、生成枝切線和繞過枝切線積分，以下對算法進行詳細描述。

首先，遍歷圖像尋找殘差點，當(dāng)掃描到殘差點時，將一個3像素×3像素的窗口放置在該殘差點上，在窗口內(nèi)尋找其他的殘差點，當(dāng)尋找到新的殘差點時將它們連接起來;接著，判斷極性，如果極性相反，則表明通過連接枝切線實現(xiàn)了殘差點的極性平衡，如果極性相同且此點未與其他殘差點相連，則將其加入到極性表中記錄其極性，并將其電荷值累計;然后，將窗口移到新的殘差點上繼續(xù)搜索，直到達到極性平衡，否則將窗口放大到5像素×5像素、7像素×7像素或更大范圍來搜索。當(dāng)搜索到邊界時，將殘差點與邊界相連。生成枝切線后，通過積分相鄰像素上的相位差分進行相位展開。其中的關(guān)鍵是:無論在窗口中找到的殘差點是否與其他殘差點相連，都將它與窗口中心的殘差點相連。另外，當(dāng)搜索窗口到達圖像邊界時，將殘差點與邊界相連，阻止積分路徑穿過。圖5是路徑跟蹤法框圖，圖6描述了枝切線設(shè)置的基本過程。

圖5 路徑跟蹤法框圖Fig.5 The flow chart of the path-following method

2 實驗結(jié)果及討論

對100層SWI相位圖片中的第1層運用Goldstein方法進行搜索，產(chǎn)生殘差點和枝切線結(jié)果。圖7是第一次對圖2(a)搜索時產(chǎn)生的殘差點結(jié)果圖像，其中被圈起及劃線的區(qū)域是被消去的對偶極子較多的區(qū)域;圖8是消除對偶極子后所剩的殘差點;圖9為消除對偶極子后產(chǎn)生的枝切線。

圖2(c)是對圖2(a)進行識別殘差點、消除對偶極子、設(shè)置枝切線后相位展開的結(jié)果圖像。圖10(a)是利用Hanning濾波器對相位偽影處理后，再利用最小密度投影獲得的圖像。圖10(b)是利用相位展開技術(shù)，對經(jīng)Hanning濾波器處理后的殘留偽影進行處理，再利用最小密度投影獲得的圖像。從圖10可以明顯看出，經(jīng)相位展開技術(shù)處理后，可以獲得更清晰的血管顯示影像。

圖11～圖13是利用枝切法對40層圖像處理的結(jié)果。

對比圖11和圖12可以看出，利用本方法可以有效地消除SWI圖像中的相位偽影，經(jīng)最小密度投影處理后，可得到清晰的血管影像(見圖13)。

另外，本研究從算法的處理速度和處理效果兩個方面，對常用的兩種算法——枝切法和質(zhì)量圖引導(dǎo)法的處理結(jié)果進行了對比實驗。從視覺效果上看，兩種算法對SWI圖像中相位偽影的處理效果基本相同。在處理速度方面，利用枝切法處理每幅SWI圖像的平均時間是476 ms，而利用質(zhì)量圖引導(dǎo)法處理每幅SWI圖像的平均處理時間為736 ms。因此，對SWI圖像中相位偽影的處理速度，枝切法比質(zhì)量圖引導(dǎo)法要高得多。

圖6 枝切線設(shè)置的基本過程。(a)正負殘差點;(b)連接正負殘差點;(c)正負殘差點;(d)3像素×3像素搜索窗;(e)連接殘差點;(f)搜索窗擴大至3像素×3像素;(g)3個殘差點;(h)連接3個殘差點;(i)將殘差點接地Fig.6 The basic process of set up branch cut line.(a)the negative＆positive residues;(b)to connect the negative＆positive residues;(c)the negative＆positive residues;(d)the 3 pixels×3 pixels searching window;(e)to connect the negative＆positive residues;(f)extending the searching window to 5 pixels×5 pixels;(g)three residues;(h)to connect the three residues;(i)grounding the residues.

圖7 第1次生成的殘差點Fig.7 The residues produced at the first time

圖8 消去對偶極子后的殘差點Fig.8 The residues after the dipole was deleted

最小二乘相位展開算法是求解二維相位展開問題最穩(wěn)健的方法之一，也是一種在雷達干涉孔徑成像中常用的相位展開技術(shù)［16］，筆者也試圖在 SWI圖像的相位偽影處理中應(yīng)用這一技術(shù)。但通過實驗發(fā)現(xiàn)，使用最小二乘法進行相位展開時，算法平均運行時間在60 s以上，運行耗時太高。又因所處理的SWI圖像為180層連續(xù)圖像，如果只進行相位展開處理就需要60×180 s的時間，而 SWI后處理還包含其他步驟。因此，在處理效果基本相同的情況下，運行時耗的多少是決定SWI的后處理算法選取的一個重要因素。由于最小二乘法相對較低的處理效率，研究中沒有選取最小二乘法對SWI圖像進行相位展開處理。

圖9 消除對偶極子后的枝切線Fig.9 The branch cut line after the dipole was deleted

圖10 經(jīng)最小密度投影處理后的圖像顯示結(jié)果。(a)有殘留偽影的血管顯示圖像;(b)消除殘留偽影后的血管顯示圖像Fig.10 The blood vessel image display results after processed by minimum intensity projection.(a)blood vessel display image with remnant phase artifact;(b)blood vessel display image after remnant phase artifact was eliminated

圖11 處理前的相位圖像Fig.11 The phase image before processed

圖12 消除偽影后的相位圖像Fig.12 The phase image after artifact be eliminated

圖13 經(jīng)最小密度投影處理后的血管圖像Fig.13 The blood vessel image display results after processed by minimum intensity projection

3 結(jié)論

SWI是一種新型的磁共振成像方式，其主要特點就是利用不同物質(zhì)間的磁敏感性差異進行靜脈的顯影成像。目前，SWI在臨床上應(yīng)用的瓶頸問題是圖像后處理技術(shù)，特別是如何消除磁場不均一性偽影問題。近年來，有關(guān)處理SWI圖像的文獻都是采用Hanning濾波器，濾除部分磁場不均一性偽影。由于經(jīng)Hanning濾波器處理后仍殘留部分偽影，嚴重影響了最終的血管顯示結(jié)果。因此，本研究首次在SWI圖像后處理中引入了在雷達干涉孔徑成像中常用的相位展開技術(shù)，專門來解決SWI相位圖像中的磁場不均一性偽影問題。實驗結(jié)果表明，相位展開技術(shù)可有效處理磁場不均一性偽影，為最終顯示清晰的血管影像奠定了基礎(chǔ)。

［1］Thomas B，Somasundaram S，Thamburaj K.Clinical applications of susceptibility weighted MR imaging of the brain-a pictorial review［J］.Neuroradiology，2008，50:105-116.

［2］Shrinivas B.SWI，a new MRI sequence—How useful it is?［J］.Ind J Radil Imag，2006，16(1):13-14.

［3］Ogawa S，Lee TM，Nayak AS，et al.，Oxygenation sensitive contrast in magnetic resonance image of rodent brain at high magnetic fields［J］.Magn Reson Med，1990，14(1):68-78.

［4］Ogawa S，Lee TM，Kay AR，et al.Brain magnetic resonance imaging with contrast dependent on blood oxygenation［J］.Proc Natl Acad Sci，1990，87(24):9868-9872.

［5］錢黎俊.磁敏感加權(quán)成像后處理算法實現(xiàn)與臨床應(yīng)用［D］.上海:上海交通大學(xué)，2007.

［6］孟智強，夏雨人.相位質(zhì)量圖在二維相位解纏中的應(yīng)用探討［J］.計算機應(yīng)用與軟件，2006，23(2):63-65.

［7］岑小林.InSAR相位解纏算法研究［D］.長沙:湖南大學(xué)，2008.

［8］李笑郁，毛士藝.干涉SAR與MRI中相位展開算法研究［J］.中國體視學(xué)與圖像分析，2001，4(6):193-197.

［9］Goldstein RM，ZebkerHA，WernerCL.Satellite radar interferometr:two-dimensional phase unwrapping［J］.Radio Sci，1988，23(4):713-720.

［10］Ghiglia DC，Romero LA.Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods［J］.J Opt Soc Amer A，1994，11(1):107-117.

［11］Pritt MD.Phase unwrapping by means of multigrid techniques for Interferometric SAR［J］.IEEE Trans on GRS，1996，34(3):728-738.

［12］孫龍.偏微分方程在合成孔徑雷達干涉測量中的應(yīng)用研究［D］.合肥:安徽大學(xué)，2005.

［13］左臣婕.分布式小衛(wèi)星SAR地面目標檢測及速度估計［D］.北京:中國科學(xué)院電子學(xué)研究所，2007.

［14］Bone DJ.Fourier fringe analysis:The two dimensional phase unwrapping problem［J］.Applied Optics，1991，30(25):3627-3632.

［15］田曉迪.基于相位展開的磁敏感加權(quán)圖像后處理方法的研究［D］.上海:上海理工大學(xué)，2009.

［16］Xu Wei，Cumming.I A region growing algorithm for InSAR phase unwrapping［J］.IEEE Transactions on Geoscience and Remote Sensing，1999，37(1):124-134.

Study on Processing Method for Phase Artifact in SWI Based on Phase Unwrapping Technique

NIE Sheng-Dong1*TIAN Xiao-Di1QIAN Li-Jun2CHEN Zeng-Ai2XU Jian-Rong2
2(Institute of Medical Imaging Engineering，University of Shanghai for Science and Technology，Shanghai 200093，China)2(Renji Hospital，Shanghai Jiaotong University，Shanghai 200127，China)

TP391

A

0258-8021(2010)04-0517-07

10.3969/j.issn.0258-8021.2010.04.007

2009-12-31，

2010-04-12

國家自然科學(xué)基金資助項目(60972122);上海市教委科研創(chuàng)新項目(09YZ216)

*通訊作者。 E-mail:nsd4647@sohu.com