為了更好地利用相位信息補償血流造成的影響，本文對磁敏感加權成像（SWI）中存在的相位纏繞問題展開了研究。為提高解纏繞的準確性，本文提出了幅度圖像引導的磁敏感加權圖像相位解纏算法。基本思路如下：① 通過改進旋轉不變非局部主成分分析濾波（PRI-NL-PCA）降低噪聲影響；② 結合 C-V 模型水平集提取相位圖像中對應的實性組織區域，從而規避背景噪聲對相位解纏方法的影響；③ 采用相位補償的方法約束 K 空間重建出的相位圖像。最后，利用四種統計量作為量化指標，評價解纏繞方法的可靠性：相位誤差的突變點個數、均值（M）、方差（Var），以及陽性百分比（Pos）和陰性百分比（Neg）。通過對比仿真數據和 226 組真實頭部磁敏感數據，結果表明，本文算法相對于經典的枝切法和最小二乘法，解纏繞結果具有較高的準確性。
引用本文： 李心靈, 黃韞梔, 韓路易, 劉奇, 何凌, 張勁, 張軍鵬, 張江. 基于幅度圖像引導的磁敏感加權圖像相位解纏算法. 生物醫學工程學雜志, 2017, 34(5): 721-729. doi: 10.7507/1001-5515.201611011 復制
版權信息： ?四川大學華西醫院華西期刊社《生物醫學工程學雜志》版權所有，未經授權不得轉載、改編
引言
磁共振成像（magnetic resonance imaging，MRI）從 K 空間采集到的信號為復數形式，臨床應用中往往利用傅里葉變換重建出對應的幅度圖像，而忽略了相位信息。實際上相位信息包含了組織間磁場的變化，能很好地反映組織中如鐵成分等其它與臨床相關的生理參數，在磁敏感加權成像（susceptibility weighted imaging，SWI）[1-2]中尤為重要。但從復數信號中提取相位時，真實相位值會被纏繞在
之間，形成相位偽影，在相位圖像上表現為黑白相間的干涉條紋。Haacke 等[1]利用高通濾波器消除相位偽影，但由于高通濾波器只是濾除偽影相位處的低頻分量，并未恢復原始相位，使得相位圖像仍有殘留偽影。因此，為恢復出真實的平滑相位信息，消除相位偽影，需要對纏繞的相位進行解纏繞。
相位解纏（phase unwrapping）在 MRI 領域中已有廣泛的應用，如磁敏感映射（susceptibility mapping）[3-4]、Dixon 法水脂分離（water/fat separation）[5]和磁共振靜脈造影術（MR venography）[6]等。目前常用的相位解纏算法可大致分為三類：① 基于路徑的解纏；② 基于區域的解纏；③ 基于函數的解纏。
Goldstein 等[7]提出的枝切法是基于路徑解纏的代表算法，由于無迭代，運算速度較快，在殘差點較少的地方可以得到較好的解纏結果，是廣泛運用的路徑依賴法。但是，Goldstein 并未考慮殘差點的極性，可能會導致枝切線的重復連接，且當殘差點密度較大時，枝切線不可避免地出現纏繞，這些都會成為算法丟棄的孤島區；而且 Goldstein 并未區分邊界和殘差點，可能導致邊界的殘差點連接到邊界，導致錯誤的積分路徑。最近鄰算法[8]對 Goldstein 進行了兩方面的改進：A）區分了正負殘差點；B）積分路徑是以極性相反且距離最近的兩個殘差點的連線作為枝切線，能夠更好地處理殘差點集中的情況。但是該算法只考慮了局部最優的積分路徑，而這并不一定是全局最優解。另一種常用的路徑依賴法是相位質量圖引導法[9]，通過建立相位差相關模型，能夠較好地規避噪聲的影響，但是不同信噪比和分辨率的 MRI 相位圖像有較強的特異性。與路徑依賴法不同，區域生長解纏法能夠從不同的方向逼近待解纏區域，能夠在一定程度上減少誤差的傳播。最經典的是 Stephan 提出結合質量圖引導和區域生長結合的 ΦUN 算法[10]，但是解纏結果依賴于質量圖確定的種子點。另一種基于區域的區域融合法[11]則無需依賴于任何種子點，但是融合邊界是利用“貪心”優化算法，易陷入局部極值點。基于函數模型的解纏算法如加權、非加權最小二乘法[12]，通過求解整幅圖像解纏繞前后相位差梯度的最小值確定真實相位，但是全局解受噪聲影響較大。
為更好地規避背景噪聲的影響，本文提出一種基于幅度圖像引導的 SWI 圖像相位解纏算法（magnitude image-guided phase unwrapping algorithm，MGPU），該算法充分利用 K 空間重建出的幅值圖像信息，從其幅值和相位角度出發，利用 C-V 模型水平集[13-14]分割 K 空間重建出的相位圖像，分別從水平和垂直方向逐點檢測相位跳變，并以垂直方向的跳變作為迭代條件，在跳變處進行相位補償。其中，C-V 模型水平集在信號較強的部分能較好地從背景中區分出檢測物體。為克服檢測物體內部噪聲和邊界的影響，更好地恢復出真實相位，假設真實相位是連續變化的且相鄰相位差小于
。
1 MGPU 算法介紹
目前針對相位解纏所提出的算法均是直接對相位圖像進行處理，很大程度上受到背景噪聲的影響，進而影響最終解纏結果。填充在 K 空間的 SWI 原始數據經傅里葉變換重建得到 SWI 復數圖像，從其數據點陣中提取幅值和相角集合得到幅度圖像和相位圖像，因此，像素空間位置一一對應，MGPU 則利用幅度圖像提取目標圖像，利用 C-V 模型水平集分割 K 空間重建出的相位圖像，對其進行相位補償，恢復圖像真實相位。算法的流程如圖 1 所示。
圖1
MGPU 算法流程圖
PRI-NL-PCA：旋轉不變非局部主成分分析濾波
Figure1. Flowchart of MGPU algorithmPRI-NL-PCA: rotational invariant non-local principle component analysis denoising
1.1 有效相位
為了規避背景噪聲，本文首先利用幅度圖像提取有效相位。因為 SWI 圖像在采集過程中不可避免地受到外界的干擾，使圖像中存在噪聲，影響后續算法的運行，所以為得到有效的幅度圖像計算模板，首先需要對采集到的幅度圖像進行濾波。隨后分割濾波后的幅度圖像，精確提取目標區域。
1.1.1 改進的 PRI-NL-PCA 濾波 MR 幅度圖像的噪聲已被證明服從 Rician 分布[15-17]。最近 Manjón 等[17]提出一種 PRI-NL-PCA 濾波方法，其在非局部均值濾波[18]（non-local means，NLM）的基礎上，基于主成分分析[19]（principle component analysis，PCA），結合 MR 圖像的稀疏性和自相似性[20]實現濾波降噪。NLM 方法已被證明能夠有效消除 MR 幅度圖像的 Rician 噪聲，但其原始算法運算效率太低；PCA 能夠在圖像特征較多的情況下實現降維，大大減少了運算量。PRI-NL-PCA 將 NLM 與 PCA 結合起來，具有運算效率高、魯棒性強等優點。因此，本文選用 PRI-NL-PCA[17]對 SWI 幅度圖像進行濾波。為了減少原算法計算量和提高計算精度，本文在原算法的基礎上做了以下兩點改進：
① 利用矩陣相似度計算相似矩陣。矩陣的相似度通過計算兩個矩陣之間的相關系數來進行比較，相關系數的絕對值越高則兩個矩陣的相似度越高，計算公式見式（1）：
 |
其中 X0 表示以原噪聲圖像 I 中一點 x0 為中心建立 K×K 大小的降噪窗口；Xj 表示 L×L 大小的訓練窗口（L>K）中大小為 K×K 的子訓練窗口；
、
分別為 X0、Xj 的均值。
② 參考文獻[21-23]中的方法對權值中的參數 h 實現自動調參。
改進后的算法能夠更好地保留圖像的邊界，且有效提高計算效率。
1.1.2 基于 C-V 模型水平集分割目標圖像 為了精確地提取濾波后的幅度圖像的有效相位，本文采用 C-V 模型水平集分割目標圖像。C-V 模型水平集將幅度圖像按邊界分成目標區域與背景區域，隨后通過最小化能量函數的方式來演化曲線[13-14]，達到分割目標與背景的目的。C-V 模型水平集具有可有效避免邊界泄露現象、抗噪能力強、對初始輪廓曲線位置不敏感等優點，在信號較強的部分能較好地從背景中區分出檢測物體。
1.2 相位補償
對提取出有效實性組織的相位圖像，需要對其纏繞的相位進行恢復，本文采用相位補償的方法實現解纏繞。
由于傅里葉變換的特性，傅里葉反變換只能恢復出
內的相位，因此從 K 空間重建出的相位圖像已產生周期性纏繞。假設實際相鄰空間位置的相位連續變化，即相鄰像素相位差小于
；反之，則表示存在較大跳變。如圖 2 所示，在點 A 附近的小空間范圍
內，真實相位理論上是連續的，但纏繞后的相位存在相位突變。
圖2
真實相位與纏繞相位關系示意圖
Figure2.
Illustration of the relationship between real and wrapped phase
為恢復出空間范圍內的實際相位，根據相位周期性變化原理，需利用式（2）確定相位圖像中每個像素點對應的纏繞周期 n。
 |
其中，
為真實相位，
為纏繞相位，
。
為確定真實相位纏繞的周期數，在相位連續的約束條件下，本文利用 Goldstein 等[7]提出的殘差點的概念，由局部到全局依次修正相位的纏繞，具體的設計流程如圖 3 所示。
圖3
相位修正流程圖
Figure3.
Flowchart of phase correction
1.2.1 局部周期修正 實性組織內部的物理和化學特性相近，在真實相位圖中理論上呈漸變趨勢，在相位纏繞結果中滿足周期連續。為保證小范圍內相位連續變化的條件，在有效相位區域內首先進行局部相位修正。根據記錄局部相位的跳變情況，進行相位周期的補償。本文考慮最近相鄰的相位作為局部修正的范圍，以保證相位連續。
① 水平方向最近鄰相位補償
計算按式（3）對水平方向上相鄰相位進行周期補償：
 |
其中
表示相位圖上點
修正后的相位；
表示對應纏繞的相位；T 表示補償單個周期；
。
② 垂直方向最近鄰相位補償
局部補償相位圖像中的行方向相位后，在沒有噪聲干擾的情況下，可認為行內的相位已正確解纏。接下來就是實現行間的周期補償，但受到有效相位的限制，相鄰行所包含的像素數存在一定的差異，不能直接利用水平方向的補償方案。為減少這一因素的影響，當相鄰行間有效相位差總和超過特定閾值范圍時，根據差異結果進行形如式（4）的列方向的周期補償。
 |
其中
，表示相鄰行間有效相位總差異的歸一化結果；
表示將總差異向上取整的結果作為補償的周期；th 表示發生行間補償的閾值條件，
，th 越小有效相位解纏繞范圍越大，但造成的誤差也隨之增大，因此其值為數據訓練的結果，本文中訓練結果為 0.75。
1.2.2 全局周期修正 經過局部修正后，由于噪聲的影響，相位圖像仍存在部分相位的纏繞，主要存在列方向的相位跳變。為克服噪聲對解纏的影響，本文從全局出發，提出了減少局部噪聲引起相位修正誤差的方案。具體思路如下：
① 根據局部解纏后的結果，確定在有效相位范圍內的最大連通域，如圖 4 所示，并以該連通域的中心線（Cen）作為基準線（Base Line）；
圖4
全局相位修正示意圖
Figure4.
Illustration of global phase correction
② 利用相位梯度檢測相位跳變邊界，并檢測跳變邊界線內部是否存在相位跳變，若有，則按局部相位補償方案進行補償；
③ 檢測距離基準線最近的行方向的相位跳變線；
④ 計算跳變邊界與基準線所在區域的相位差，并加以補償；
⑤ 重復 ②—④，直到整幅圖片不存在跳變線。
2 實驗材料與方法
2.1 仿真圖片
在 Matlab 中生成 512×512 大小的仿真圖片，使用半徑為 200 像素的圓形模板，對其添加方差為 50 的隨機噪聲，并加至仿真圖片中作為背景。將目標區域原始相位范圍設置為
，對其添加方差分別為 0.4、1 和 2 的乘性噪聲，通過式（2）計算得到纏繞后的相位圖，分別如圖 5 所示。纏繞相位圖的相位范圍為
，即每一個黑白相間的相位圈內相位跨越范圍均為
。
圖5
512×512 不同噪聲方差仿真圖
Figure5.
Simulation images with size of 512×512 and different noise
2.2 真實 SWI 圖像
真實 SWI 人體數據由 Alltech 公司提供，數據格式為.dat 文件，采自 MRI Comfort 1.5T 機型，8 通道頭線圈，掃描參數設置為：TE/TR=40 ms/50 ms，反轉角 20°，激勵次數 1 次，層厚 2 mm，層間距 0 mm，層數 32，采集矩陣 200×1 024。
真實 SWI 圖像共 226 組原始 K 空間數據，對每一組數據均分別采用枝切法、最小二乘法和 MGPU 法進行實驗，實驗中將 SWI 原始 K 空間數據重建為 256×256 大小的幅度圖像和相位圖像。文中選取分別來自 8 層 1 通道、12 層 7 通道、20 層 2 通道和 27 層 8 通道的四組不同層數、不同通道的數據對實驗結果進行簡要闡述。
2.3 評價方法
為驗證本文提出的 MGPU 法的有效性及魯棒性，在 MATLAB 環境中分別對仿真圖片和真實腦部 SWI 圖像進行測試，并與枝切法和基于離散余弦變換的最小二乘法進行對比，計算解纏相位與真實相位之間的相位誤差和再纏繞相位與纏繞相位之間的相位誤差，本文認為相位誤差在
，
之間的點近似為零，超過該范圍的點定義為突變點，以突變點的統計結果作為算法的量化指標：
① 計算突變點的均值（M）：值越接近于零表明算法越優越；
② 計算突變點的方差（Var）：值越接近于零表明相位誤差變化越小，算法越穩定；
③ 計算突變點的陰、陽性百分比 Pos（%）、Neg（%）：定義解纏相位小于真實相位的點為陰性點，解纏相位大于真實相位的點為陽性點，分別計算兩種突變點占全部像素點個數的百分比，百分比之和（P+N）越小表明解纏繞結果越可靠、算法越有效。
3 實驗結果與分析
3.1 仿真實驗結果與分析
為了評價算法的有效性，分別用枝切法、最小二乘法和 MGPU 法對圖 5 中已知相位分布的不同乘性噪聲方差纏繞相位仿真圖進行實驗，并統計量化指標值。
3.1.1 無噪聲影響 為了討論噪聲的存在對解纏算法的影響，本文首先對無噪聲影響的仿真圖進行實驗，并統計其量化指標值，實驗結果和指標值如圖 6 和表 1 所示。
圖6
仿真圖解纏繞結果
Figure6.
The unwrapping results of simulation images
表1
無噪聲影響的量化指標結果
Table1.
Indicator values of noise-free simulation image
圖 6 可直觀說明，枝切法解纏相位圖能夠恢復大部分真實相位，但邊界由于受到背景噪聲的影響出現了一些“孤島”區域，并且再纏繞相位圖與原始纏繞相位圖有一定的偏差；最小二乘法由于受到背景噪聲的影響，其邊界和背景不易區分，幾乎融合在一起，其再纏繞相位圖不能得到理想的結果；MGPU 法能夠恢復真實相位并再纏繞得到纏繞相位。
為了統計有效的評價指標，實驗中將前文所描述的相位圖像目標區域計算模板 Mask 用于枝切法的計算中，最小二乘法由于其全局優化的本質，故考慮全局的影響。表 1 數據表明，在無噪聲影響下，枝切法的 P+N 為 47.81%，突變點個數為 125 324，值得注意的是，計算模板 Mask 的像素點個數即為 125 324，即枝切法未能對無噪聲圖像實現正確解纏繞。最小二乘法的 P+N 為 100%，突變點個數為 262 144，意味著所有點均未正確解纏繞。而 MGPU 法的 M、Var、P+N 均為 0，即說明在誤差允許的范圍內 MGPU 法能夠對無噪聲圖像實現完全解纏繞。
3.1.2 有噪聲影響 對噪聲方差分別為 0.4、1 和 2 的乘性噪聲仿真相位圖進行實驗，并統計其量化指標值，實驗結果如表 2 和表 3 所示，其中真實相位是未添加噪聲時的仿真圖片原始相位，解纏相位是對添加噪聲仿真圖的解纏結果，再纏繞相位是對解纏相位使用式（2）所得結果，纏繞相位是添加噪聲仿真圖的相位。表 2 記錄了突變點個數，由于最小二乘法的突變點占 100%，故只統計了枝切法與 MGPU 法在不同噪聲方差下的突變點個數。表 3 統計了不同噪聲影響下的量化指標結果。
表2
突變點的個數
Table2.
Number of error points
表3
不同噪聲影響下的量化指標結果
Table3.
Indicator values of simulation image with different noise
數據結果表明，無論是解纏相位與真實相位，還是再纏繞相位與纏繞相位，枝切法的突變點個數均為 125 324，P+N 為 47.81%，即為計算模板 Mask 的大小，表明枝切法不能正確恢復真實相位和纏繞相位。值得注意的是，前文提到枝切法由于受到背景噪聲的影響，導致邊界出現“孤島”區域，表 2 和表 3 的數據結果未能證明目標區域的解纏相位是否受到背景噪聲的影響，但不能排除其影響。MGPU 法的解纏繞相位與真實相位的突變點 P+N 也為 47.81%，說明算法不能恢復受到了噪聲污染的相位，但是已經被污染的噪聲相位通過 MGPU 則可以完全恢復出來，不受背景噪聲的影響。隨著噪聲方差的增大，對比三種算法的 M、Var、P+N 指標可得，MGPU 法解纏繞準確性最高，枝切法次之，最小二乘法最差。
3.2 真實 SWI 圖像實驗結果與分析
3.2.1 去噪結果 如圖 7 所示分別為原始幅度圖像、經改進的 PRI-NL-PCA 法處理后得到的濾波圖像，以及計算模板 Mask。其中目標區域能夠被有效并完整地提取，規避了背景噪聲。
圖7
原始幅度圖像處理結果
Figure7.
Results of original magnitude images
3.2.2 三種方法實驗結果對比 如圖 8 第二列所示為枝切法處理結果，圖中目標區域邊界附近殘差點較多，噪聲較大，導致枝切法不能有效地解纏繞，這部分區域亦稱之為“孤島”區域。此外，枝切法的運算過程中需要選擇一個已知點的真實相位值，并且積分路徑與枝切線的設置密切相關，這樣使得得到的結果必然有誤差。第三列為基于最小二乘法的處理結果，由于最小二乘法是一種全局最優算法，其局部誤差易傳遞至整幅圖像，其背景區域與目標區域幾乎融合在了一起，受背景噪聲影響較大。最后一列為 MGPU 算法處理結果，在幅度圖像的引導下消除了背景噪聲對解纏結果的影響，從整體出發依次由局部到全局進行周期修正，避免了“孤島”現象的產生，并能得到較理想的解纏繞結果。
圖8
相位圖像解纏繞結果
Figure8.
Unwrapping results of real phase images
為了評價算法的有效性，現將圖 8 所示的解纏結果通過式（2）計算得到再纏繞結果，兩者相減統計再纏繞相位與纏繞相位之間的相位誤差評價指標，如表 4 所示。
表4
真實 SWI 數據的量化指標統計結果
Table4.
The indicators values of true SWI data
表 4 統計了上述四幅圖像的量化指標結果以及 226 組數據的平均評價指標值，結果表明，針對不同的圖像枝切法在不同程度上殘留突變點，在誤差允許的范圍內其 P+N 為 0.10%，即完全解纏繞程度可達 99.90%，而 MGPU 法可達 100%，優于枝切法。
經統計 226 組數據經枝切法處理后的突變點相關指標，其所有數據突變點總數為 3 401，平均每組數據殘留約 15 個突變點，均值 M 為 15.048 7，方差 Var 為 31.939 4。突變點大于 M 的數據有 60 組，占 26.55%；大于 Var 的數據有 30 組，占 13.37%；突變點最多的一組為 2 層 5 通道，有 443 個突變點，P+N 達 3.37%，遠超突變點的 M 和 Var；突變點為零的數據有 10 組，占 4.42%，即有 10 組數據在誤差允許的范圍內經枝切法完全正確解纏繞。
上述數據結果表明，在相位誤差允許的范圍內，MGPU 法無殘留突變點，即完全正確解纏繞的能力可達 100%。而枝切法仍不同程度地殘留突變點，完全正確解纏繞的可能性僅為 4.42%，其原因很大程度上是受背景噪聲的影響及其算法本身的影響，算法把整幅相位圖像作為輸入，未規避背景噪聲。
4 討論與總結
針對 SWI 相位圖像存在的相位偽影以及相位纏繞問題，本文提出一種基于幅度圖像引導的 SWI 圖像相位解纏算法，并對 PRI-NL-PCA 濾波算法做了一些改進，簡化運算復雜度并保留圖像邊緣信息；對幅度圖像濾波后利用 C-V 模型水平集分割幅度圖像目標圖像，提取相位圖像實性組織，規避背景噪聲的干擾；在相位連續的約束條件下，由局部到全局依次修正相位的纏繞，保證整體解纏結果有效且可靠。實驗結果表明 MGPU 能夠在規避背景噪聲的同時有效消除明暗相間的干涉條紋狀相位偽影，恢復出原始相位值。
為了準確量化實驗結果，本文提出基于相位誤差的突變點的四組量化指標，分別是 M、Var、Pos 和 Neg，并用于仿真實驗和真實 SWI 數據中，其值越接近于零表明突變點越少，解纏繞效果越好。仿真實驗表明：① 在無噪聲影響下，MGPU 法能夠完全實現正確解纏繞；枝切法對目標區域的點的解纏繞存在偏差；而最小二乘法由于其函數優化的本質，所有點均未能正確解纏繞。② 隨著噪聲方差的增大，MGPU 法對目標區域的解纏繞仍能有效完成，枝切法和最小二乘法的解纏繞結果偏差增大。③ 在不同程度噪聲的影響下，MGPU 法仍然能夠完全恢復噪聲相位，表明噪聲僅僅是污染了真實相位，但并未對 MGPU 算法本身造成影響。真實 SWI 數據表明，在相位誤差允許的范圍內，MGPU 法能夠 100% 實現正確解纏繞；枝切法仍然不同程度地殘留突變點，對一組數據正確解纏繞的可能性為 4.42%，每組數據突變點出現的可能性為 0.10%；最小二乘法則未能正確解纏繞。
MGPU 法具有以下幾個顯著優點：① 有效規避背景噪聲，減少噪聲對相位解纏繞的影響；② 進入目標區域的補償周期初值可控，不存在局部最優；③ 由局部到全局依次修正周期避免產生“孤島”現象。總地來說，MGPU 具有較強魯棒性，能有效且可靠地對 SWI 相位圖像進行解纏繞。
致謝：感謝成都奧泰醫療系統有限責任公司（AllTech，http://www.alltechmed.com/）為本文研究提供 1.5T SWI 原始數據。
引言
磁共振成像（magnetic resonance imaging，MRI）從 K 空間采集到的信號為復數形式，臨床應用中往往利用傅里葉變換重建出對應的幅度圖像，而忽略了相位信息。實際上相位信息包含了組織間磁場的變化，能很好地反映組織中如鐵成分等其它與臨床相關的生理參數，在磁敏感加權成像（susceptibility weighted imaging，SWI）[1-2]中尤為重要。但從復數信號中提取相位時，真實相位值會被纏繞在
之間，形成相位偽影，在相位圖像上表現為黑白相間的干涉條紋。Haacke 等[1]利用高通濾波器消除相位偽影，但由于高通濾波器只是濾除偽影相位處的低頻分量，并未恢復原始相位，使得相位圖像仍有殘留偽影。因此，為恢復出真實的平滑相位信息，消除相位偽影，需要對纏繞的相位進行解纏繞。
相位解纏（phase unwrapping）在 MRI 領域中已有廣泛的應用，如磁敏感映射（susceptibility mapping）[3-4]、Dixon 法水脂分離（water/fat separation）[5]和磁共振靜脈造影術（MR venography）[6]等。目前常用的相位解纏算法可大致分為三類：① 基于路徑的解纏；② 基于區域的解纏；③ 基于函數的解纏。
Goldstein 等[7]提出的枝切法是基于路徑解纏的代表算法，由于無迭代，運算速度較快，在殘差點較少的地方可以得到較好的解纏結果，是廣泛運用的路徑依賴法。但是，Goldstein 并未考慮殘差點的極性，可能會導致枝切線的重復連接，且當殘差點密度較大時，枝切線不可避免地出現纏繞，這些都會成為算法丟棄的孤島區；而且 Goldstein 并未區分邊界和殘差點，可能導致邊界的殘差點連接到邊界，導致錯誤的積分路徑。最近鄰算法[8]對 Goldstein 進行了兩方面的改進：A）區分了正負殘差點；B）積分路徑是以極性相反且距離最近的兩個殘差點的連線作為枝切線，能夠更好地處理殘差點集中的情況。但是該算法只考慮了局部最優的積分路徑，而這并不一定是全局最優解。另一種常用的路徑依賴法是相位質量圖引導法[9]，通過建立相位差相關模型，能夠較好地規避噪聲的影響，但是不同信噪比和分辨率的 MRI 相位圖像有較強的特異性。與路徑依賴法不同，區域生長解纏法能夠從不同的方向逼近待解纏區域，能夠在一定程度上減少誤差的傳播。最經典的是 Stephan 提出結合質量圖引導和區域生長結合的 ΦUN 算法[10]，但是解纏結果依賴于質量圖確定的種子點。另一種基于區域的區域融合法[11]則無需依賴于任何種子點，但是融合邊界是利用“貪心”優化算法，易陷入局部極值點。基于函數模型的解纏算法如加權、非加權最小二乘法[12]，通過求解整幅圖像解纏繞前后相位差梯度的最小值確定真實相位，但是全局解受噪聲影響較大。
為更好地規避背景噪聲的影響，本文提出一種基于幅度圖像引導的 SWI 圖像相位解纏算法（magnitude image-guided phase unwrapping algorithm，MGPU），該算法充分利用 K 空間重建出的幅值圖像信息，從其幅值和相位角度出發，利用 C-V 模型水平集[13-14]分割 K 空間重建出的相位圖像，分別從水平和垂直方向逐點檢測相位跳變，并以垂直方向的跳變作為迭代條件，在跳變處進行相位補償。其中，C-V 模型水平集在信號較強的部分能較好地從背景中區分出檢測物體。為克服檢測物體內部噪聲和邊界的影響，更好地恢復出真實相位，假設真實相位是連續變化的且相鄰相位差小于
。
1 MGPU 算法介紹
目前針對相位解纏所提出的算法均是直接對相位圖像進行處理，很大程度上受到背景噪聲的影響，進而影響最終解纏結果。填充在 K 空間的 SWI 原始數據經傅里葉變換重建得到 SWI 復數圖像，從其數據點陣中提取幅值和相角集合得到幅度圖像和相位圖像，因此，像素空間位置一一對應，MGPU 則利用幅度圖像提取目標圖像，利用 C-V 模型水平集分割 K 空間重建出的相位圖像，對其進行相位補償，恢復圖像真實相位。算法的流程如圖 1 所示。
圖1
MGPU 算法流程圖
PRI-NL-PCA：旋轉不變非局部主成分分析濾波
Figure1. Flowchart of MGPU algorithmPRI-NL-PCA: rotational invariant non-local principle component analysis denoising
1.1 有效相位
為了規避背景噪聲，本文首先利用幅度圖像提取有效相位。因為 SWI 圖像在采集過程中不可避免地受到外界的干擾，使圖像中存在噪聲，影響后續算法的運行，所以為得到有效的幅度圖像計算模板，首先需要對采集到的幅度圖像進行濾波。隨后分割濾波后的幅度圖像，精確提取目標區域。
1.1.1 改進的 PRI-NL-PCA 濾波 MR 幅度圖像的噪聲已被證明服從 Rician 分布[15-17]。最近 Manjón 等[17]提出一種 PRI-NL-PCA 濾波方法，其在非局部均值濾波[18]（non-local means，NLM）的基礎上，基于主成分分析[19]（principle component analysis，PCA），結合 MR 圖像的稀疏性和自相似性[20]實現濾波降噪。NLM 方法已被證明能夠有效消除 MR 幅度圖像的 Rician 噪聲，但其原始算法運算效率太低；PCA 能夠在圖像特征較多的情況下實現降維，大大減少了運算量。PRI-NL-PCA 將 NLM 與 PCA 結合起來，具有運算效率高、魯棒性強等優點。因此，本文選用 PRI-NL-PCA[17]對 SWI 幅度圖像進行濾波。為了減少原算法計算量和提高計算精度，本文在原算法的基礎上做了以下兩點改進：
① 利用矩陣相似度計算相似矩陣。矩陣的相似度通過計算兩個矩陣之間的相關系數來進行比較，相關系數的絕對值越高則兩個矩陣的相似度越高，計算公式見式（1）：
|
其中 X0 表示以原噪聲圖像 I 中一點 x0 為中心建立 K×K 大小的降噪窗口；Xj 表示 L×L 大小的訓練窗口（L>K）中大小為 K×K 的子訓練窗口；
、
分別為 X0、Xj 的均值。
② 參考文獻[21-23]中的方法對權值中的參數 h 實現自動調參。
改進后的算法能夠更好地保留圖像的邊界，且有效提高計算效率。
1.1.2 基于 C-V 模型水平集分割目標圖像 為了精確地提取濾波后的幅度圖像的有效相位，本文采用 C-V 模型水平集分割目標圖像。C-V 模型水平集將幅度圖像按邊界分成目標區域與背景區域，隨后通過最小化能量函數的方式來演化曲線[13-14]，達到分割目標與背景的目的。C-V 模型水平集具有可有效避免邊界泄露現象、抗噪能力強、對初始輪廓曲線位置不敏感等優點，在信號較強的部分能較好地從背景中區分出檢測物體。
1.2 相位補償
對提取出有效實性組織的相位圖像，需要對其纏繞的相位進行恢復，本文采用相位補償的方法實現解纏繞。
由于傅里葉變換的特性，傅里葉反變換只能恢復出
內的相位，因此從 K 空間重建出的相位圖像已產生周期性纏繞。假設實際相鄰空間位置的相位連續變化，即相鄰像素相位差小于
；反之，則表示存在較大跳變。如圖 2 所示，在點 A 附近的小空間范圍
內，真實相位理論上是連續的，但纏繞后的相位存在相位突變。
圖2
真實相位與纏繞相位關系示意圖
Figure2.
Illustration of the relationship between real and wrapped phase
為恢復出空間范圍內的實際相位，根據相位周期性變化原理，需利用式（2）確定相位圖像中每個像素點對應的纏繞周期 n。
|
其中，
為真實相位，
為纏繞相位，
。
為確定真實相位纏繞的周期數，在相位連續的約束條件下，本文利用 Goldstein 等[7]提出的殘差點的概念，由局部到全局依次修正相位的纏繞，具體的設計流程如圖 3 所示。
圖3
相位修正流程圖
Figure3.
Flowchart of phase correction
1.2.1 局部周期修正 實性組織內部的物理和化學特性相近，在真實相位圖中理論上呈漸變趨勢，在相位纏繞結果中滿足周期連續。為保證小范圍內相位連續變化的條件，在有效相位區域內首先進行局部相位修正。根據記錄局部相位的跳變情況，進行相位周期的補償。本文考慮最近相鄰的相位作為局部修正的范圍，以保證相位連續。
① 水平方向最近鄰相位補償
計算按式（3）對水平方向上相鄰相位進行周期補償：
|
其中
表示相位圖上點
修正后的相位；
表示對應纏繞的相位；T 表示補償單個周期；
。
② 垂直方向最近鄰相位補償
局部補償相位圖像中的行方向相位后，在沒有噪聲干擾的情況下，可認為行內的相位已正確解纏。接下來就是實現行間的周期補償，但受到有效相位的限制，相鄰行所包含的像素數存在一定的差異，不能直接利用水平方向的補償方案。為減少這一因素的影響，當相鄰行間有效相位差總和超過特定閾值范圍時，根據差異結果進行形如式（4）的列方向的周期補償。
|
其中
，表示相鄰行間有效相位總差異的歸一化結果；
表示將總差異向上取整的結果作為補償的周期；th 表示發生行間補償的閾值條件，
，th 越小有效相位解纏繞范圍越大，但造成的誤差也隨之增大，因此其值為數據訓練的結果，本文中訓練結果為 0.75。
1.2.2 全局周期修正 經過局部修正后，由于噪聲的影響，相位圖像仍存在部分相位的纏繞，主要存在列方向的相位跳變。為克服噪聲對解纏的影響，本文從全局出發，提出了減少局部噪聲引起相位修正誤差的方案。具體思路如下：
① 根據局部解纏后的結果，確定在有效相位范圍內的最大連通域，如圖 4 所示，并以該連通域的中心線（Cen）作為基準線（Base Line）；
圖4
全局相位修正示意圖
Figure4.
Illustration of global phase correction
② 利用相位梯度檢測相位跳變邊界，并檢測跳變邊界線內部是否存在相位跳變，若有，則按局部相位補償方案進行補償；
③ 檢測距離基準線最近的行方向的相位跳變線；
④ 計算跳變邊界與基準線所在區域的相位差，并加以補償；
⑤ 重復 ②—④，直到整幅圖片不存在跳變線。
2 實驗材料與方法
2.1 仿真圖片
在 Matlab 中生成 512×512 大小的仿真圖片，使用半徑為 200 像素的圓形模板，對其添加方差為 50 的隨機噪聲，并加至仿真圖片中作為背景。將目標區域原始相位范圍設置為
，對其添加方差分別為 0.4、1 和 2 的乘性噪聲，通過式（2）計算得到纏繞后的相位圖，分別如圖 5 所示。纏繞相位圖的相位范圍為
，即每一個黑白相間的相位圈內相位跨越范圍均為
。
圖5
512×512 不同噪聲方差仿真圖
Figure5.
Simulation images with size of 512×512 and different noise
2.2 真實 SWI 圖像
真實 SWI 人體數據由 Alltech 公司提供，數據格式為.dat 文件，采自 MRI Comfort 1.5T 機型，8 通道頭線圈，掃描參數設置為：TE/TR=40 ms/50 ms，反轉角 20°，激勵次數 1 次，層厚 2 mm，層間距 0 mm，層數 32，采集矩陣 200×1 024。
真實 SWI 圖像共 226 組原始 K 空間數據，對每一組數據均分別采用枝切法、最小二乘法和 MGPU 法進行實驗，實驗中將 SWI 原始 K 空間數據重建為 256×256 大小的幅度圖像和相位圖像。文中選取分別來自 8 層 1 通道、12 層 7 通道、20 層 2 通道和 27 層 8 通道的四組不同層數、不同通道的數據對實驗結果進行簡要闡述。
2.3 評價方法
為驗證本文提出的 MGPU 法的有效性及魯棒性，在 MATLAB 環境中分別對仿真圖片和真實腦部 SWI 圖像進行測試，并與枝切法和基于離散余弦變換的最小二乘法進行對比，計算解纏相位與真實相位之間的相位誤差和再纏繞相位與纏繞相位之間的相位誤差，本文認為相位誤差在
，
之間的點近似為零，超過該范圍的點定義為突變點，以突變點的統計結果作為算法的量化指標：
① 計算突變點的均值（M）：值越接近于零表明算法越優越；
② 計算突變點的方差（Var）：值越接近于零表明相位誤差變化越小，算法越穩定；
③ 計算突變點的陰、陽性百分比 Pos（%）、Neg（%）：定義解纏相位小于真實相位的點為陰性點，解纏相位大于真實相位的點為陽性點，分別計算兩種突變點占全部像素點個數的百分比，百分比之和（P+N）越小表明解纏繞結果越可靠、算法越有效。
3 實驗結果與分析
3.1 仿真實驗結果與分析
為了評價算法的有效性，分別用枝切法、最小二乘法和 MGPU 法對圖 5 中已知相位分布的不同乘性噪聲方差纏繞相位仿真圖進行實驗，并統計量化指標值。
3.1.1 無噪聲影響 為了討論噪聲的存在對解纏算法的影響，本文首先對無噪聲影響的仿真圖進行實驗，并統計其量化指標值，實驗結果和指標值如圖 6 和表 1 所示。
圖6
仿真圖解纏繞結果
Figure6.
The unwrapping results of simulation images
表1
無噪聲影響的量化指標結果
Table1.
Indicator values of noise-free simulation image
圖 6 可直觀說明，枝切法解纏相位圖能夠恢復大部分真實相位，但邊界由于受到背景噪聲的影響出現了一些“孤島”區域，并且再纏繞相位圖與原始纏繞相位圖有一定的偏差；最小二乘法由于受到背景噪聲的影響，其邊界和背景不易區分，幾乎融合在一起，其再纏繞相位圖不能得到理想的結果；MGPU 法能夠恢復真實相位并再纏繞得到纏繞相位。
為了統計有效的評價指標，實驗中將前文所描述的相位圖像目標區域計算模板 Mask 用于枝切法的計算中，最小二乘法由于其全局優化的本質，故考慮全局的影響。表 1 數據表明，在無噪聲影響下，枝切法的 P+N 為 47.81%，突變點個數為 125 324，值得注意的是，計算模板 Mask 的像素點個數即為 125 324，即枝切法未能對無噪聲圖像實現正確解纏繞。最小二乘法的 P+N 為 100%，突變點個數為 262 144，意味著所有點均未正確解纏繞。而 MGPU 法的 M、Var、P+N 均為 0，即說明在誤差允許的范圍內 MGPU 法能夠對無噪聲圖像實現完全解纏繞。
3.1.2 有噪聲影響 對噪聲方差分別為 0.4、1 和 2 的乘性噪聲仿真相位圖進行實驗，并統計其量化指標值，實驗結果如表 2 和表 3 所示，其中真實相位是未添加噪聲時的仿真圖片原始相位，解纏相位是對添加噪聲仿真圖的解纏結果，再纏繞相位是對解纏相位使用式（2）所得結果，纏繞相位是添加噪聲仿真圖的相位。表 2 記錄了突變點個數，由于最小二乘法的突變點占 100%，故只統計了枝切法與 MGPU 法在不同噪聲方差下的突變點個數。表 3 統計了不同噪聲影響下的量化指標結果。
表2
突變點的個數
Table2.
Number of error points
表3
不同噪聲影響下的量化指標結果
Table3.
Indicator values of simulation image with different noise
數據結果表明，無論是解纏相位與真實相位，還是再纏繞相位與纏繞相位，枝切法的突變點個數均為 125 324，P+N 為 47.81%，即為計算模板 Mask 的大小，表明枝切法不能正確恢復真實相位和纏繞相位。值得注意的是，前文提到枝切法由于受到背景噪聲的影響，導致邊界出現“孤島”區域，表 2 和表 3 的數據結果未能證明目標區域的解纏相位是否受到背景噪聲的影響，但不能排除其影響。MGPU 法的解纏繞相位與真實相位的突變點 P+N 也為 47.81%，說明算法不能恢復受到了噪聲污染的相位，但是已經被污染的噪聲相位通過 MGPU 則可以完全恢復出來，不受背景噪聲的影響。隨著噪聲方差的增大，對比三種算法的 M、Var、P+N 指標可得，MGPU 法解纏繞準確性最高，枝切法次之，最小二乘法最差。
3.2 真實 SWI 圖像實驗結果與分析
3.2.1 去噪結果 如圖 7 所示分別為原始幅度圖像、經改進的 PRI-NL-PCA 法處理后得到的濾波圖像，以及計算模板 Mask。其中目標區域能夠被有效并完整地提取，規避了背景噪聲。
圖7
原始幅度圖像處理結果
Figure7.
Results of original magnitude images
3.2.2 三種方法實驗結果對比 如圖 8 第二列所示為枝切法處理結果，圖中目標區域邊界附近殘差點較多，噪聲較大，導致枝切法不能有效地解纏繞，這部分區域亦稱之為“孤島”區域。此外，枝切法的運算過程中需要選擇一個已知點的真實相位值，并且積分路徑與枝切線的設置密切相關，這樣使得得到的結果必然有誤差。第三列為基于最小二乘法的處理結果，由于最小二乘法是一種全局最優算法，其局部誤差易傳遞至整幅圖像，其背景區域與目標區域幾乎融合在了一起，受背景噪聲影響較大。最后一列為 MGPU 算法處理結果，在幅度圖像的引導下消除了背景噪聲對解纏結果的影響，從整體出發依次由局部到全局進行周期修正，避免了“孤島”現象的產生，并能得到較理想的解纏繞結果。
圖8
相位圖像解纏繞結果
Figure8.
Unwrapping results of real phase images
為了評價算法的有效性，現將圖 8 所示的解纏結果通過式（2）計算得到再纏繞結果，兩者相減統計再纏繞相位與纏繞相位之間的相位誤差評價指標，如表 4 所示。
表4
真實 SWI 數據的量化指標統計結果
Table4.
The indicators values of true SWI data
表 4 統計了上述四幅圖像的量化指標結果以及 226 組數據的平均評價指標值，結果表明，針對不同的圖像枝切法在不同程度上殘留突變點，在誤差允許的范圍內其 P+N 為 0.10%，即完全解纏繞程度可達 99.90%，而 MGPU 法可達 100%，優于枝切法。
經統計 226 組數據經枝切法處理后的突變點相關指標，其所有數據突變點總數為 3 401，平均每組數據殘留約 15 個突變點，均值 M 為 15.048 7，方差 Var 為 31.939 4。突變點大于 M 的數據有 60 組，占 26.55%；大于 Var 的數據有 30 組，占 13.37%；突變點最多的一組為 2 層 5 通道，有 443 個突變點，P+N 達 3.37%，遠超突變點的 M 和 Var；突變點為零的數據有 10 組，占 4.42%，即有 10 組數據在誤差允許的范圍內經枝切法完全正確解纏繞。
上述數據結果表明，在相位誤差允許的范圍內，MGPU 法無殘留突變點，即完全正確解纏繞的能力可達 100%。而枝切法仍不同程度地殘留突變點，完全正確解纏繞的可能性僅為 4.42%，其原因很大程度上是受背景噪聲的影響及其算法本身的影響，算法把整幅相位圖像作為輸入，未規避背景噪聲。
4 討論與總結
針對 SWI 相位圖像存在的相位偽影以及相位纏繞問題，本文提出一種基于幅度圖像引導的 SWI 圖像相位解纏算法，并對 PRI-NL-PCA 濾波算法做了一些改進，簡化運算復雜度并保留圖像邊緣信息；對幅度圖像濾波后利用 C-V 模型水平集分割幅度圖像目標圖像，提取相位圖像實性組織，規避背景噪聲的干擾；在相位連續的約束條件下，由局部到全局依次修正相位的纏繞，保證整體解纏結果有效且可靠。實驗結果表明 MGPU 能夠在規避背景噪聲的同時有效消除明暗相間的干涉條紋狀相位偽影，恢復出原始相位值。
為了準確量化實驗結果，本文提出基于相位誤差的突變點的四組量化指標，分別是 M、Var、Pos 和 Neg，并用于仿真實驗和真實 SWI 數據中，其值越接近于零表明突變點越少，解纏繞效果越好。仿真實驗表明：① 在無噪聲影響下，MGPU 法能夠完全實現正確解纏繞；枝切法對目標區域的點的解纏繞存在偏差；而最小二乘法由于其函數優化的本質，所有點均未能正確解纏繞。② 隨著噪聲方差的增大，MGPU 法對目標區域的解纏繞仍能有效完成，枝切法和最小二乘法的解纏繞結果偏差增大。③ 在不同程度噪聲的影響下，MGPU 法仍然能夠完全恢復噪聲相位，表明噪聲僅僅是污染了真實相位，但并未對 MGPU 算法本身造成影響。真實 SWI 數據表明，在相位誤差允許的范圍內，MGPU 法能夠 100% 實現正確解纏繞；枝切法仍然不同程度地殘留突變點，對一組數據正確解纏繞的可能性為 4.42%，每組數據突變點出現的可能性為 0.10%；最小二乘法則未能正確解纏繞。
MGPU 法具有以下幾個顯著優點：① 有效規避背景噪聲，減少噪聲對相位解纏繞的影響；② 進入目標區域的補償周期初值可控，不存在局部最優；③ 由局部到全局依次修正周期避免產生“孤島”現象。總地來說，MGPU 具有較強魯棒性，能有效且可靠地對 SWI 相位圖像進行解纏繞。
致謝：感謝成都奧泰醫療系統有限責任公司（AllTech，http://www.alltechmed.com/）為本文研究提供 1.5T SWI 原始數據。
