李心靈 唐 銘 劉 奇

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

引言

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

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

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

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

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

由于時域解纏繞算法受到上述諸多限制，本研究僅對空域解纏繞算法進行分析與比較，并利用仿真數(shù)據(jù)和真實SWI數(shù)據(jù)對算法進行測試，根據(jù)評價指標完成對算法的評估。

1 方法

1.1 基本概念

1.1.1相位解纏繞

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

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

(1)

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

1.1.2殘差

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

(2)

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

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

1.2 空域相位解纏繞算法

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

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

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

1.2.1.1WLS

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

(3)

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

1.2.1.2PRELUDE

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

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

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

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

(4)

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

1.2.1.3PUMA

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

(5)

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

1.2.1.4ARM

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

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

(6)

式中，δ為冗余相位。

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

(7)

1.2.1.5SL-MC

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

1)如果像素ri沒有被修正過的相鄰點，則它被單獨分為一組，ri+1以ri為標準進行修正并分為一組，且新組的權重為ri的幅值。

2)如果ri的相鄰點已被修正，則ri被加至相鄰點的分組，該組的權重為該組所有像素對應的幅值。

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

1.2.2積分算法

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

1.2.2.1Branch-cut

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

1.2.2.2WFF-QG

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

1.2.2.3PUROR

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

1.3 驗證方法

1.3.1實驗材料

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

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

1.3.2評價方法

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

1)計算突變點的均值M。M表示突變點的值與真實值之間的大小差異，值越接近于零表明算法解纏繞的結果越精確。

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

3)運行時間time(s)。

2 結果

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

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

表1 仿真圖片評價指標

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

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

圖2 真實SWI圖像實驗結果

實驗圖像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的結果受到其幅度圖像的模板影響較大，為了更公平地評價實驗結果，精確地提取目標相位的模板，本研究利用簡單的PM擴散模型[26]和Ostu閾值分割法[27]提取幅度圖像模板，改進的模板較PRELUDE傳統(tǒng)模板更為精確，能有效并完整地提取目標圖像邊界，并應用在PRELUDE、PUMA、SL-MC和PUROR這4種方法中，實驗結果和評價指標分別如圖3和表3所示。

圖3 真實SWI圖像實驗結果。(a)Phantom;(b)人腦圖像

表3 改進模板后4種算法的運行時間(s)

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

3 討論

為了找到能有效解決SWI相位解纏繞問題的方法，本研究在MRI空域解纏繞算法中選取8種不同類型的方法，分別對兩組仿真圖像和兩組真實SWI圖像進行測試，以突變點和運行時間作為算法的評價指標，對8種算法進行評價。

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

表4 8種算法對比分析

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

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

4 結論

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

(致謝：感謝成都奧泰醫(yī)療系統(tǒng)有限責任公司(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.