李心靈 唐 銘 劉 奇

(四川大學(xué)電氣信息學(xué)院，成都 610065)

引言

磁共振相位信息在磁共振成像(magnetic resonance imaging，MRI)的臨床應(yīng)用中往往被忽略，實(shí)際上相位信息包含了大量組織間的局部敏感性變化，能夠用于測(cè)量例如鐵成分[1]等其他與臨床相關(guān)的物質(zhì)成分，在磁敏感加權(quán)成像(susceptibility weighted imaging, SWI)[1-3]中尤為重要。SWI是以T2*加權(quán)梯度回波序列作為序列基礎(chǔ)，根據(jù)不同組織間的磁敏感性差異對(duì)比增強(qiáng)目標(biāo)區(qū)域，因此相位信息在SWI中起著決定性作用。在相位圖像中，往往存在相位纏繞的問題，使得相位圖像存在偽影，不能真實(shí)反映組織的特性，所以需要對(duì)相位圖像進(jìn)行解纏繞處理。實(shí)際上，在合成孔徑雷達(dá)[4-5]、條紋投影[6-7]和光學(xué)干涉測(cè)量[8]等領(lǐng)域的學(xué)者對(duì)相位解纏繞技術(shù)的研究已有幾十年，隨著相位信息在MRI領(lǐng)域中的廣泛應(yīng)用，相位解纏繞算法對(duì)于減少磁共振圖像的相位偽影具有重要的研究?jī)r(jià)值。

近年來，研究學(xué)者針對(duì)MRI相位解纏繞提出了大量的算法，一般將MRI相位解纏繞算法分為兩大類：空域和時(shí)域。空域解纏繞算法假設(shè)相鄰像素之間相位變化緩慢，通常不超過半個(gè)周期，即π。然而，由于噪聲、欠采樣或者物體本身不連續(xù)等情況的存在，假設(shè)條件并不一定處處滿足，使得相位解纏繞面臨許多困難，這些不滿足假設(shè)條件的點(diǎn)稱為“殘差”。與此同時(shí)，時(shí)域解纏繞算法受到的限制較多，首先必須使用多次回波序列或一次回波序列采集多組信號(hào)，且所得到的解纏繞圖像并不是相位圖像，而是相位差圖像，這些原因都限制了時(shí)域解纏繞方法的廣泛應(yīng)用。本研究將空域解纏繞算法分為兩類：全局優(yōu)化算法和積分算法。

1)全局優(yōu)化算法。該類算法通常將解纏繞問題轉(zhuǎn)換為最小化代價(jià)函數(shù)進(jìn)行解纏繞，如最小范數(shù)法(minimum lp-norm)[9-10]、解纏繞離散估計(jì)的相位區(qū)域擴(kuò)張標(biāo)記法(phase region expanding labeler for unwrapping discrete estimations，PRELUDE)[11]、差值映射累加法(accumulation of residual maps，ARM)[12]等；或者通過在馬爾科夫隨機(jī)場(chǎng)中優(yōu)化能量函數(shù)實(shí)現(xiàn)解纏繞，如最大流/最小切相位解纏繞法(phase unwrapping max-flow/min-cut，PUMA)[13]、基于模型解纏繞[14]等；或者不通過優(yōu)化代價(jià)函數(shù)解纏繞，而是依賴于相鄰相位點(diǎn)的特征解纏繞，如排序多聚類法(sorted list，multi-clustering，SL-MC法)[15]等。

2)積分算法。該類算法以Itoh于1982年提出的積分理論[16]為基礎(chǔ)，采用積分的方式進(jìn)行解纏繞，如枝切法(branch-cut method)[17]、最優(yōu)路徑法(best-path method)[18]、基于傅里葉濾波的質(zhì)量圖引導(dǎo)法(windowed Fourier-filtered and quality-guided method，WFF-QG)[19]、PHUN法[20]、遞歸正交參考相位解纏繞法(phase unwrapping using recursive orthogonal referring，PUROR)[21]等。

時(shí)域解纏繞算法通常通過采集兩個(gè)或多個(gè)回波信號(hào)，恢復(fù)由時(shí)變效應(yīng)(time-dependent effects)引起的纏繞，如UMPIRE法[22](unwrapping multi-echo phase images with irregular echo spacings)、MAGPI法[23]等。

由于時(shí)域解纏繞算法受到上述諸多限制，本研究?jī)H對(duì)空域解纏繞算法進(jìn)行分析與比較，并利用仿真數(shù)據(jù)和真實(shí)SWI數(shù)據(jù)對(duì)算法進(jìn)行測(cè)試，根據(jù)評(píng)價(jià)指標(biāo)完成對(duì)算法的評(píng)估。

1 方法

1.1 基本概念

1.1.1相位解纏繞

SWI相位圖像通過反正切計(jì)算，所得到的原始測(cè)量相位值的范圍被限制在(-π,π]區(qū)間內(nèi)，位于該區(qū)間外的相位通過加減整數(shù)倍2π被重新纏繞至該區(qū)間，對(duì)纏繞相位進(jìn)行一系列處理使其恢復(fù)至原始相位的過程稱為相位解纏繞。纏繞相位Φ(r)與真實(shí)相位φ(r)之間的關(guān)系如下：

Φ(r)=W[φ(r)]=φ(r)+2πn

(1)

式中：r是相位圖像中任意一點(diǎn)；W[·]是纏繞算子，其計(jì)算結(jié)果在(-π,π]區(qū)間內(nèi)；n是任意整數(shù)，同樣使得計(jì)算結(jié)果在(-π,π]區(qū)間內(nèi)。

1.1.2殘差

假設(shè)SWI相位圖像不受噪聲等因素影響，沿著任意閉合路徑對(duì)纏繞相位差進(jìn)行積分，得到的結(jié)果均為0，即積分與路徑無關(guān)[16]，其過程可以表示為

(2)

式中，ΔΦ(r)是相鄰兩點(diǎn)的相位差，m是閉合路徑的總像素點(diǎn)。

但是，真實(shí)的SWI相位圖像會(huì)受到噪聲或欠采樣等的影響。按順時(shí)針或逆時(shí)針的方向計(jì)算相鄰4個(gè)點(diǎn)的纏繞相位差之和，如果其值不為0，則該閉合環(huán)稱為殘差(residue)[17]。任何包含一個(gè)或多個(gè)殘差的閉合路徑在解纏繞圖像中會(huì)出現(xiàn)間斷區(qū)域，而對(duì)于存在殘差的纏繞相位圖現(xiàn)在還沒有研究出較好的相位解纏繞辦法。

1.2 空域相位解纏繞算法

本研究將空域解纏繞算法分為全局優(yōu)化算法和積分算法兩大類，并選取一共8種不同解纏繞的算法進(jìn)行分析與比較。

1.2.1全局優(yōu)化算法

全局優(yōu)化算法從整體相位出發(fā)，不處理殘差，通常通過將解纏繞問題轉(zhuǎn)換為優(yōu)化代價(jià)函數(shù)的方式進(jìn)行解纏繞。全局優(yōu)化算法中最具代表性的算法是由Ghiglia等于1994年提出的最小二乘法[9]，該算法分為加權(quán)最小二乘法(weighted least square，WLS)與無加權(quán)最小二乘法，相比之下，加權(quán)最小二乘法通過引入權(quán)重，減少了異常值對(duì)整體解纏繞結(jié)果的影響。Jenkinson于2003年提出的PRELUDE法[11]是目前魯棒性較強(qiáng)的一種算法，在之后提出的大量算法中，它均作為對(duì)比算法與結(jié)果進(jìn)行比較?；趫D割的相位解纏繞算法是近10年提出的新方法，Jose等于2007年提出PUMA算法[13]，在馬爾科夫隨機(jī)場(chǎng)中解決纏繞問題，魯棒性較好。Rivera等于2015年提出ARM算法[12]，它是魯棒性較強(qiáng)的新算法，能得到較好的解纏繞結(jié)果。Maier等于2015年提出SL-MC算法[15]，不是通過優(yōu)化代價(jià)函數(shù)進(jìn)行解纏繞，而是采用分組加權(quán)的方式解纏繞，該算法也是一種較新的全局算法。

1.2.1.1WLS

WLS將相位解纏繞問題轉(zhuǎn)化成范數(shù)形式函數(shù)的最小化問題，通過最小化解纏繞相位梯度與對(duì)應(yīng)的纏繞相位梯度之間的平方差總和，得到最終的解纏繞結(jié)果，其數(shù)學(xué)表達(dá)式為

(3)

式中，wij為權(quán)重，一般來自質(zhì)量圖或掩碼。

1.2.1.2PRELUDE

PRELUDE是一種區(qū)域融合的多維相位解纏繞技術(shù)，基于區(qū)域融合的方法來最優(yōu)化代價(jià)函數(shù)，使得邊界處的相位差保持在半個(gè)周期內(nèi)。PRELUDE按照3個(gè)步驟完成相位解纏繞：

1)設(shè)置相位區(qū)間，將多維空間數(shù)據(jù)分成多個(gè)子區(qū)域，并認(rèn)為每個(gè)區(qū)域內(nèi)不存在相位纏繞，例如以π/3為一個(gè)間隔，將(-π，π]區(qū)間劃分為6個(gè)小區(qū)間；

2)將每個(gè)區(qū)域視為一個(gè)整體并對(duì)其進(jìn)行編號(hào)，根據(jù)式(4)計(jì)算兩兩交界區(qū)域的代價(jià)函數(shù)；

3)根據(jù)代價(jià)函數(shù)的結(jié)果決定融合順序，通過迭代的方法融合子區(qū)域直至所有區(qū)域融合，有

(4)

式中：MAB=MA-MB，MA與MB分別為交界區(qū)域A與B的整數(shù)偏移量；j和k分別為A和B區(qū)域中的兩個(gè)相鄰體素，N(j)為j的一個(gè)簡(jiǎn)單連通域。

1.2.1.3PUMA

PUMA實(shí)際上是一種圖割的最大流/最小割算法，通過建立馬爾科夫能量模型，推導(dǎo)出能量函數(shù)如下：

(5)

然后利用圖割理論求解最優(yōu)化的問題。

1.2.1.4ARM

ARM是一種不依賴路徑的全局優(yōu)化算法，在算法中嵌套多重網(wǎng)格法進(jìn)行迭代，最小化代價(jià)函數(shù)實(shí)現(xiàn)解纏繞。ARM不是直接對(duì)纏繞相位圖進(jìn)行優(yōu)化，而是優(yōu)化冗余相位圖實(shí)現(xiàn)解纏繞。假設(shè)現(xiàn)在存在初始解纏繞相位圖f0，則有

φ(r)=f0(r)+δ(r)

(6)

式中，δ為冗余相位。

λ(δr-δs)2]+μ(1-ωrs)2}

(7)

1.2.1.5SL-MC

SL-MC首先將MR復(fù)數(shù)信號(hào)的像素按照幅值大小遞減的順序排列建立一個(gè)list，再按照下列方式進(jìn)行處理：

1)如果像素ri沒有被修正過的相鄰點(diǎn)，則它被單獨(dú)分為一組，ri+1以ri為標(biāo)準(zhǔn)進(jìn)行修正并分為一組，且新組的權(quán)重為ri的幅值。

2)如果ri的相鄰點(diǎn)已被修正，則ri被加至相鄰點(diǎn)的分組，該組的權(quán)重為該組所有像素對(duì)應(yīng)的幅值。

3)如果ri有多個(gè)不同分組的相鄰點(diǎn)，這些分組之間可能存在相位偏移，那么組別之間將進(jìn)行修正，修正的方式為權(quán)重較大的組融合權(quán)重較小的組。當(dāng)list中的所有像素均被修正時(shí)，SL-MC則完成解纏繞。

1.2.2積分算法

積分算法以Itoh的積分理論為基礎(chǔ)，從殘差的分布出發(fā)，設(shè)法繞過殘差解纏繞。處理殘差的方式?jīng)Q定了算法的分類，目前使用較多的方法有枝切線法和質(zhì)量引導(dǎo)法。枝切線法根據(jù)殘差的分布將殘差連起來形成枝切線，隨后避開枝切線對(duì)圖像積分，Goldstein等于1988年提出的Branch-cut法[17]是枝切線法的代表算法，之后此類的大量算法均是在枝切線的設(shè)置上進(jìn)行改進(jìn)。質(zhì)量引導(dǎo)法不直接處理殘差，通常采用不同的方式獲得質(zhì)量圖，根據(jù)質(zhì)量圖中各個(gè)像素點(diǎn)的質(zhì)量引導(dǎo)出合理的積分路徑來解纏繞，Qian等于2008提出的WFF-QG法未采用相位圖像生成質(zhì)量圖，而是將濾波后的幅度圖像作為質(zhì)量圖引導(dǎo)解纏繞[19]。此外，有的路徑積分算法并不通過尋找路徑或者質(zhì)量引導(dǎo)的方法來解纏繞，如Liu等于2012年提出的PUROR法從全局出發(fā)對(duì)積分圖像進(jìn)行約束解纏繞[21]。

1.2.2.1Branch-cut

Branch-cut是枝切線法的代表算法，通過計(jì)算殘差點(diǎn)，連接正負(fù)殘差點(diǎn)，設(shè)置枝切線，使電荷達(dá)到平衡，然后繞過枝切線對(duì)纏繞相位進(jìn)行積分解纏繞。Branch-cut最大的一個(gè)缺點(diǎn)是：在殘差點(diǎn)較多的地方，算法直接放棄這部分區(qū)域形成“孤島”。此外，Branch-cut連接正負(fù)殘差點(diǎn)時(shí)未考慮最優(yōu)連接方式，因此枝切線往往設(shè)置得較復(fù)雜。

1.2.2.2WFF-QG

WFF-QG是一種質(zhì)量圖引導(dǎo)的解纏繞算法，其首先利用窗口傅里葉變換(windowed Fourier transform，WFT)對(duì)幅度圖像進(jìn)行濾波，然后將經(jīng)過濾波后的幅度圖像作為實(shí)值質(zhì)量圖進(jìn)行解纏繞。

1.2.2.3PUROR

PUROR在Itoh理論的基礎(chǔ)上，分別對(duì)纏繞圖像從行、列兩個(gè)方向積分，然后分別設(shè)定兩個(gè)不同的閾值來約束兩幅初始積分圖像，利用遞歸正交參考的方法使兩幅初始積分圖像最終得到相同的解纏繞結(jié)果。

1.3 驗(yàn)證方法

1.3.1實(shí)驗(yàn)材料

在Matlab中分別生成兩組不同的仿真圖：一組為參照文獻(xiàn)[15]中的方法生成的復(fù)數(shù)環(huán)形仿真圖，SNR按順時(shí)針方向從6開始每30°遞減0.5，如圖1所示；另一組為存在相位不連續(xù)的仿真圖，并對(duì)其添加不同方差的隨機(jī)噪聲，使之產(chǎn)生5組不同數(shù)量的殘差，其殘差所占百分比分別為0%、0.5%、1%、1.5%、2%，通過式(1)計(jì)算得到纏繞后的相位圖如圖1所示。兩組仿真纏繞相位圖的相位范圍均為(-π,π]，即每一個(gè)黑白相間的相位圈內(nèi)相位跨越范圍均為|2π|。

除了仿真圖片，兩組由Alltech公司提供的真實(shí)SWI圖像被用于測(cè)試算法，一組為ACR體膜數(shù)據(jù)，一組為真實(shí)人腦數(shù)據(jù)。兩組數(shù)據(jù)的格式為.dat文件，采自MRI Comfort 1.5T機(jī)型，8通道頭線圈。掃描參數(shù)設(shè)置為：TE/TR=40 ms/50 ms，反轉(zhuǎn)角20°，激勵(lì)次數(shù)1次，層厚2 mm，層間距0 mm，層數(shù)66，采集矩陣200×1024。本研究采用VRC法[25]，將每一層8個(gè)通道的幅度圖像和相位圖像分別進(jìn)行重建，得到共264幅相位圖像和幅度圖像。

1.3.2評(píng)價(jià)方法

為有效評(píng)價(jià)算法的魯棒性，在Matlab 2016 a環(huán)境中分別對(duì)仿真圖片和真實(shí)SWI圖像進(jìn)行測(cè)試，其中PRELUDE算法在由英國(guó)牛津大學(xué)發(fā)布的FSL(FMRIB Software Library)軟件中進(jìn)行測(cè)試。計(jì)算解纏繞相位與真實(shí)相位之間的相位誤差和再纏繞相位與原始纏繞相位之間的相位誤差，本研究將誤差大于|0.001π|的點(diǎn)定義為突變點(diǎn)，即在規(guī)定誤差范圍內(nèi)未完全恢復(fù)原始相位的點(diǎn)，并以突變點(diǎn)的統(tǒng)計(jì)結(jié)果作為算法的量化指標(biāo)。

1)計(jì)算突變點(diǎn)的均值M。M表示突變點(diǎn)的值與真實(shí)值之間的大小差異，值越接近于零表明算法解纏繞的結(jié)果越精確。

2)計(jì)算突變點(diǎn)的百分比P(%)。P表示突變點(diǎn)數(shù)量百分比，即未正確解纏繞的點(diǎn)占整幅相位圖像的百分比，P越小表明解纏繞結(jié)果越可靠，算法越有效。

3)運(yùn)行時(shí)間time(s)。

2 結(jié)果

仿真圖實(shí)驗(yàn)結(jié)果和評(píng)價(jià)指標(biāo)統(tǒng)計(jì)結(jié)果分別如圖1和表1所示。對(duì)復(fù)數(shù)圖像來說，WLS出現(xiàn)明顯的解纏繞誤差，表1數(shù)據(jù)顯示有72.58%的像素點(diǎn)未有效解纏繞。PRELUDE根據(jù)幅度圖像提取目標(biāo)相位，在SNR較低的區(qū)域無法提取有效目標(biāo)相位，使得部分相位缺失，但已解纏繞的區(qū)域不存在突變點(diǎn)，即該區(qū)域被有效解纏繞。PUMA、SL-MC和PUROR 3種方法均能有效解纏繞，在低SNR區(qū)域的解纏繞結(jié)果各不相同，其中PUROR殘留明顯的誤差“碎片”，但表1中其P±|M|值為零，說明此“碎片”區(qū)域與真實(shí)解纏繞相位之間相差整數(shù)倍周期。圖1中ARM和WFF-QG的解纏繞效果較好，但其突變點(diǎn)百分比均在99%以上。Branch-cut缺失低SNR區(qū)域，且有53.61%的突變點(diǎn)。從運(yùn)行時(shí)間上來說，PRELUDE運(yùn)行速度最快，PUMA最慢。

圖1 仿真圖的解纏繞結(jié)果。(a)復(fù)數(shù)環(huán)形仿真圖，幅值圖像和相位圖像由復(fù)數(shù)圖像計(jì)算得到；(B)存在不同殘差數(shù)

表1 仿真圖片評(píng)價(jià)指標(biāo)

對(duì)存在不同殘差數(shù)的仿真圖片來說，當(dāng)不存在殘差時(shí)，PRELUDE、PUMA和SL-MC 3種方法的結(jié)果相似；WLS、Branch-cut和PUROR的突變點(diǎn)為0，但其結(jié)果與原始圖像存在差異；ARM和WFF-QG均存在突變點(diǎn)，前者為99.91%，后者為14.27%，且WFF-QG的解纏繞圖像與PRELUDE的相似。隨著殘差數(shù)增大，PUMA的實(shí)驗(yàn)結(jié)果最穩(wěn)定，其余算法的結(jié)果均與原始仿真圖之間存在一定的差異，其中PRELUDE的時(shí)間代價(jià)最大，且殘差數(shù)越多運(yùn)行時(shí)間越長(zhǎng)。PUROR的運(yùn)行時(shí)間最短。

圖2為部分真實(shí)SWI圖像的實(shí)驗(yàn)結(jié)果，第1行與第2行分別為第1層和第61層的phantom圖像，第3行與第4行分別為第13層和第54層的真實(shí)人腦圖像。表2統(tǒng)計(jì)了圖2和全部真實(shí)SWI圖像實(shí)驗(yàn)的評(píng)價(jià)指標(biāo)，其中第1行和第4行的P±|M|為單層圖像平均值，T為66層圖像總運(yùn)行時(shí)間。結(jié)果顯示，PRELUDE運(yùn)行時(shí)間最短，但其受本身幅度圖像模板的影響，導(dǎo)致目標(biāo)相位部分丟失，在剩余的目標(biāo)相位中，其解纏繞準(zhǔn)確性較高。PUMA、SL-MC和PUROR的解纏繞圖像誤差較少，但PUMA和SL-MC的時(shí)間代價(jià)太大，PUROR的運(yùn)行時(shí)間比PRELUDE慢20 s左右。WLS運(yùn)行時(shí)間較快，但解纏繞圖像中的誤差“碎片”較多，結(jié)果不理想。ARM和WFF-QG的突變點(diǎn)幾乎在99%以上，運(yùn)行時(shí)間在200 s以上，且WFF-QG的目標(biāo)與背景難以區(qū)分，邊界有重影出現(xiàn)。Branch-cut也存在突變點(diǎn)，且圖像中殘留有“孤島”區(qū)域。

圖2 真實(shí)SWI圖像實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)圖像WLSPRELUDEPUMAARMSL-MCBranch-cutWFF-QGPURORP±|M|T/sP±|M|T/sP±|M|T/sP±|M|T/sP±|M|T/sP±|M|T/sP±|M|T/sP±|M|T/sPhantom—97.2—40.6—877.799.92±0.05245.9—590.235.61±0.07385.896.04±0.002480.2—75.2第1層—1.8—0.1—11.799.95±0.124.1—8.429.80±0.936.297.48±0.0026.9—1.4第61層—1.9—0.1—12.199.92±0.023.9—8.530.34±0.796.197.95±0.0066.3—1.4人腦—95.4—77.1—748.699.91±0.02248.6—572.563.80±0.06215.299.36±0.008487.5—93.7第13層—1.8—0.2—9.799.93±0.153.7—7.960.25±0.043.699.24±0.0226.4—1.7第54層—1.8—0.3—12.899.83±0.013.7—8.153.23±1.024.499.22±0.0136.3—1.6

由于PRELUDE的結(jié)果受到其幅度圖像的模板影響較大，為了更公平地評(píng)價(jià)實(shí)驗(yàn)結(jié)果，精確地提取目標(biāo)相位的模板，本研究利用簡(jiǎn)單的PM擴(kuò)散模型[26]和Ostu閾值分割法[27]提取幅度圖像模板，改進(jìn)的模板較PRELUDE傳統(tǒng)模板更為精確，能有效并完整地提取目標(biāo)圖像邊界，并應(yīng)用在PRELUDE、PUMA、SL-MC和PUROR這4種方法中，實(shí)驗(yàn)結(jié)果和評(píng)價(jià)指標(biāo)分別如圖3和表3所示。

圖3 真實(shí)SWI圖像實(shí)驗(yàn)結(jié)果。(a)Phantom;(b)人腦圖像

表3 改進(jìn)模板后4種算法的運(yùn)行時(shí)間(s)

使用本研究提取的模板后，PRELUDE對(duì)phantom和人腦實(shí)驗(yàn)的運(yùn)行時(shí)間分別增加了2/3和1/2，對(duì)比原算法提取的模板，改進(jìn)后的模板對(duì)缺失的相位區(qū)域進(jìn)行了補(bǔ)充，即目標(biāo)相位區(qū)域增大，說明目標(biāo)圖像的大小會(huì)影響其運(yùn)行時(shí)間；圖3表明，PRELUDE存在一定的誤差， 但大部分實(shí)驗(yàn)情況良好，魯棒性依然很高。PUMA的運(yùn)行時(shí)間均減少了約1/3，運(yùn)行結(jié)果未受到影響。SL-MC的運(yùn)行時(shí)間均減少了約100 s，運(yùn)行結(jié)果與圖2中的有所差異，相比之下圖2的結(jié)果誤差“碎片”更少。PUROR的運(yùn)行時(shí)間均減少了約1/2，同時(shí)其誤差“碎片”增多。

3 討論

為了找到能有效解決SWI相位解纏繞問題的方法，本研究在MRI空域解纏繞算法中選取8種不同類型的方法，分別對(duì)兩組仿真圖像和兩組真實(shí)SWI圖像進(jìn)行測(cè)試，以突變點(diǎn)和運(yùn)行時(shí)間作為算法的評(píng)價(jià)指標(biāo)，對(duì)8種算法進(jìn)行評(píng)價(jià)。

根據(jù)8種算法的實(shí)驗(yàn)結(jié)果，表4列出了對(duì)比分析結(jié)果。WLS由于其范數(shù)最小化的本質(zhì)[9]，使得存在較多誤差“碎片”；PUMA和SL-MC受背景噪聲影響較大；ARM、Branch-cut和WFF-QG的P±|M|值不為0，真實(shí)相位未被恢復(fù)；PUMA、ARM、SL-MC和WFF-QG這4種方法的運(yùn)行時(shí)間較長(zhǎng)；PRELUDE和PUROR的P±|M|值為0，算法運(yùn)行時(shí)間相對(duì)較短，解纏繞效果較好。

表4 8種算法對(duì)比分析

由實(shí)驗(yàn)結(jié)果可知，PRELUDE和PUROR兩種方法的魯棒性較強(qiáng)，但均有各自的優(yōu)勢(shì)及問題存在。PRELUDE的解纏繞過程依賴幅度圖像，若不利用幅度圖像提取模板、減少計(jì)算量，則運(yùn)行時(shí)間會(huì)大大增加；幅度圖像模板的精確度同樣也會(huì)影響PRELUDE的運(yùn)行結(jié)果與時(shí)間，若模板不精確，則會(huì)出現(xiàn)圖2中缺失目標(biāo)相位的問題。此外，PRELUDE尚存在一個(gè)相位區(qū)間的擇優(yōu)問題[11]，若相位區(qū)間設(shè)置較大，則有可能存在纏繞相位，產(chǎn)生計(jì)算誤差；若相位區(qū)間設(shè)置較小，則會(huì)增加運(yùn)行時(shí)間。因此，如何設(shè)置一個(gè)合適的相位區(qū)間，是PRELUDE需要考慮的問題。PUROR法是8種算法中運(yùn)行時(shí)間最短的方法，其計(jì)算時(shí)間與圖像數(shù)量成正相關(guān)，從效率上來說占據(jù)較大的優(yōu)勢(shì)；此外，PUROR不需要依賴集合高質(zhì)量像素點(diǎn)的模板，只需要從全局出發(fā)對(duì)圖像進(jìn)行解纏繞，考慮到了圖像中所有的像素點(diǎn)。但PUROR仍存在一定的問題，比如對(duì)于SNR極低的圖像，PUROR將失效[21]。

本研究將MRI空域解纏繞算法分為全局優(yōu)化算法和積分算法兩類，從突變點(diǎn)的指標(biāo)來看，全局優(yōu)化算法通過優(yōu)化代價(jià)函數(shù)的方式，最小化殘差點(diǎn)的冗余值，使得P±|M|值能保持為0；積分算法由于Itoh積分理論的本質(zhì)是設(shè)法繞過殘差點(diǎn)[16]，使得殘差點(diǎn)未被處理而P±|M|不為0。但是，全局優(yōu)化算法容易出現(xiàn)誤差“碎片”。從運(yùn)行時(shí)間來看，由于每種算法的具體計(jì)算過程不同，所以不能簡(jiǎn)單地用這一個(gè)指標(biāo)評(píng)價(jià)兩類算法。

4 結(jié)論

目前在實(shí)際臨床應(yīng)用中，SWI圖像相位處理仍普遍采用Haccke提出的零差濾波法[1]。該方法最大的特點(diǎn)是能在一定程度上濾除相位偽影且運(yùn)行時(shí)間短，但其對(duì)于相位纏繞的問題不能有效解決，仍存在誤差。因此，在臨床應(yīng)用中可參考相位解纏繞算法，對(duì)相位圖像的纏繞問題進(jìn)行解決。本研究對(duì)8種空域解纏繞算法進(jìn)行了對(duì)比，結(jié)果表明PUROR法魯棒性強(qiáng)，運(yùn)行時(shí)間短。實(shí)際臨床應(yīng)用中，在考慮準(zhǔn)確性和運(yùn)行時(shí)間的情況下，可嘗試采用PUROR法。但與零差濾波法相比，相位解纏繞算法的時(shí)間代價(jià)相對(duì)較高。因此，找出解纏繞效果好且運(yùn)行時(shí)間短的相位解纏繞方法，是今后相位解纏繞算法的研究重點(diǎn)。

(致謝：感謝成都奧泰醫(yī)療系統(tǒng)有限責(zé)任公司(AllTech，http://www.alltechmed.com/)為本研究提供1.5T SWI原始數(shù)據(jù))

[1] Haacke EM. Susceptibility weighted imaging (SWI) [J]. Magnetic Resonance in Medicine, 2004, 52(3):612-618.

[2] Di IA, Le RP, Carsin-Nicol B, et al. Diagnostic value of fractal analysis for the differentiation of brain tumors using 3-tesla magnetic resonance susceptibility-weighted imaging [J]. Neurosurgery, 2016, 79(6): 839-846.

[3] Heyn C, Alcaide-Leon P, Bharatha A, et al. Susceptibility-weighted imaging in neurovascular disease [J]. Top Magn Reson Imaging, 2016, 25(2): 63-71.

[4] Pritt MD. Phase unwrapping by means of multigrid techniques for interferometric SAR [J]. IEEE Transactions on Geoscience and Remote Sensing, 1996, 34(3): 728-738.

[5] Chen CW, Zebker HA. Phase unwrapping for large SAR interferograms: Statistical segmentation and generalized network models [J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(8): 1709-1719.

[6] Judge TR, Bryanston-Cross PJ. A review of phase unwrapping techniques in fringe analysis [J]. Optics and Lasers in Engineering, 1994, 21(4): 199-239.

[7] Zhang H, Lalor MJ, Burton DR. Spatiotemporal phase unwrapping for the measurement of discontinuous objects in dynamic fringe-projection phase-shifting profilometry [J]. Applied Optics, 1999, 38(16): 3534-3541.

[8] Schofield MA, Zhu Y. Fast phase unwrapping algorithm for interferometric applications [J]. Optics Letters, 2003, 28(14): 1194-1196.

[9] Ghiglia DC, Romero LA. Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods [J]. Journal of the Optical Society of America A, 1994, 11(1):107-117.

[10] Ghiglia DC, Romero LA. Minimum LP-norm two-dimensional phase unwrapping [J]. J Opt Soc Am A, 1999, 13(10):1999-2013.

[11] Mark J. Fast, automated, N-dimensional phase-unwrapping algorithm [J]. Magnetic resonance in medicine, 2003, 49(1): 193-197.

[12] Rivera M, Hernandez-Lopez FJ, Gonzalez A. Phase unwrapping by accumulation of residual maps [J]. Optics & Lasers in Engineering, 2015, 64(64):51-58.

[13] Bioucasdias JM, Valad?o G. Phase unwrapping via graph cuts [J]. IEEE Transactions on Image Processing, 2007, 16(3):698-709.

[14] Ying L, Liang ZP, David C, et al. Unwrapping of MR phase images using a Markov random field model [J]. IEEE Trans Med Imaging, 2006, 25(1):128-136.

[15] Maier F, Fuentes D, Weinberg JS, et al. Robust phase unwrapping for MR temperature imaging using a magnitude-sorted list, multi-clustering algorithm [J]. Magnetic Resonance in Medicine, 2015, 73(4):1662-1668.

[16] K. Itoh. Analysis of the phase unwrapping algorithm [J]. Appl Opt, 1982, 21(14): 2470-2470.

[17] Goldstein RM, Zebker HA, Werner CL. Satellite radar interferometry: Two-dimensional phase unwrapping [J]. Radio science, 1988, 23(4):713-720.

[18] Abdul-Rahman HS, Gdeisat MA, Burton DR, et al. Fast and robust three-dimensional best path phase unwrapping algorithm [J]. Appl Opt, 2007, 46(26): 6623-6635.

[19] Kemao Q, Gao W, Wang H. Windowed Fourier-filtered and quality-guided phase-unwrapping algorithm [J]. Applied Optics, 2008, 47(29): 5420-5428.

[20] Witoszynskyj S, Rauscher A, Reichenbach JR, et al. Phase unwrapping of MR images using PHUN-A fast and robust region growing algorithm [J]. Medical Image Analysis, 2009, 13(2): 257-268.

[21] Junmin L,Maria D. Intervention-based multidimensional phase unwrapping using recursive orthogonal referring [J]. Magnetic Resonance in Medicine, 2012, 68(4): 1303-1316.

[22] Robinson S, Schodl H, Trattnig S. A method for unwrapping highly wrapped multi-echo phase images at very high field: UMPIRE [J]. Magn Reson Med, 2014, 72(1): 80-92.

[23] Dagher J, Reese T. Bilgin A. High-resolution, large dynamic range field map estimation [J]. Magn Reson Med, 2014, 71(1):105-117.

[24] Rivera M, Marroquin JL. Half-quadratic cost functions for phase unwrapping [J]. Optics Letters, 2004, 29(5):504-506.

[25] Parker DL, Payne A, Todd N, et al. Phase reconstruction from multiple coil data using a virtual reference coil [J]. Magnetic Resonance in Medicine, 2014, 72(2): 563-569.

[26] Perona P, Malik J. Scale-space and edge detection using anisotropic diffusing [J]. IEEE Transaction on Pattern Analysis and Machine Intelligence, 1990, 12(7): 629-639.

[27] Otsu N. A threshold selection method from gray-level histograms [J]. IEEE Transactions on Systems Man and Cybernetics, 1975, 9(1): 63-66.