結合峰值檢測的前列腺磁共振成像偏移場平滑擬合

楊雄，詹曙*，謝棟棟

結合峰值檢測的前列腺磁共振成像偏移場平滑擬合

楊雄1，詹曙1*，謝棟棟2

目的 研究前列腺磁共振圖像中灰度不均勻現象(偏場)的校正方法。材料與方法 從幾組前列腺磁共振掃描數據中截取的橫斷面圖像。將真實圖像的分段常量特性與偏移場的平滑變化特性表達到圖像模型中，構造一個能量函數，通過能量函數的最小化實現偏移場評估和組織分割。利用峰值檢測技術自動獲得能量函數的初始化參數，并用結合三角函數與多項式函數的一組基函數實現對偏移場的平滑擬合。結果 定性的實驗表明筆者的方法能對前列腺磁共振圖像中的偏場現象進行有效的校正。另外通過與其他方法在變化系數、均方根、Jaccard相似度等指標下的定量對比發現，筆者的方法有更好的校正結果。結論 結合峰值檢測的偏移場校正方法能對前列腺磁共振圖像中的灰度不均勻現象有效改善。

磁共振成像；前列腺；峰值檢測；偏場校正；平滑擬合

1 引言

前列腺癌已成為嚴重危害老年男性健康的疾病，而且據相關報道，中國前列腺癌發病率正呈上升趨勢[1-2]。前列腺癌的早期診斷和檢測對疾病的治療預防非常重要。考慮到傳統前列腺癌早期檢測方法的局限性，一些借助于MRI技術的人工檢測方法，以及結合前列腺癌的病理學知識與MR影像特征的計算機輔助診斷方案也在研究之中[3-7]。

這些借助于MR影像的診斷方法對于影像的質量有較高要求，然而醫學圖像中普遍存在著一種灰度不均勻現象，也即偏場。偏移場的存在會影響病灶區的視覺效果以及計算機對病灶組織的判斷。根據文獻[8-9]，偏移場的產生是由于成像設備和患者特異性的綜合作用。目前對偏移場的校正方法大致可以分為前瞻性方法和回顧性方法[9]。前瞻性方法是從MR圖像獲得過程的角度來考慮校正方法，這類方法雖然也能取得一些效果，但是往往不能解決由于患者特異性所帶來的影響。回顧性方法只依賴于獲得圖像的信息，因而能更全面地解決偏場，當然也得到了更多的關注和研究。Johnston等[10]最先提出一種基于同態濾波的方法，這個方法假設偏移場具有低頻特性，但是圖像本身也具有很多低頻成分，因此高通濾波器也不能取得理想的結果。Miles等[11]提出一個結合直方圖信息和B樣條曲面擬合的方法。Li等[12]在乘法固有成分優化的思想下，提出一個能量最小化模型，同時實現偏場矯正和多組織分割。

目前對于偏場矯正的方法很多，本文在已有算法的基礎上作了一些改進，提出一種針對前列腺組織MR圖像的偏場校正方法。首先根據前列腺磁共振圖像的灰度特征進行自適應的預處理，然后利用峰值檢測技術從圖像的直方圖信息中自動獲得前列腺組織的類別參數，接著根據圖像模型構造能量最小化模型，最后提出一種由最簡多項式函數與三角函數組合而成的曲面擬合方法，實現對偏場的平滑擬合。

2 相關闡述

2.1 MR圖像模型

目前一種普遍認可的磁共振圖像模型如式(1)所示：

其中x代表圖像像素點，v(x)表示獲得的圖像，u(x)代表沒有偏場與噪聲作用的真實圖像，b(x)是偏場，n(x)則指圖像噪聲。通常認為圖像中的背景與前景中的各個目標在灰度上具有一致性，而圖像中的灰度不均勻變化是一個平滑過程。所以在這個圖像模型中，筆者認為代表真實圖像的u(x)具有分段常量特性，即圖像中不同的組織可以近似為不同的灰度常量區域，而表征圖像不均勻現象的b(x)具有平滑特性。

2.2 能量函數

針對u(x)的分段常量特性，假設圖像中存在n個不同的組織，每個組織的灰度值為常量ci(i=1,...,n)。隸屬度函數μi(x)表示像素點x屬于第i類組織的概率。為了更真實地反映實像，這里參照模糊C均值(fuzzy c-means algorithm, FCM)算法采用一種模糊表示方法()。那么實像u(x)表示如下：

其中q是加權指數。對于偏場，筆者用一系列連續函數的加權和表示：

m的取值決定著曲面的擬合精度。通常對于偏場的擬合會用多項式函數作為連續函數，但是筆者發現三角函數有著等同高次函數的平滑特性，而且計算更簡。所以用低次最簡多項式函數與三角函數共同作為擬合曲面的基函數。

令a=(a1,...,am)T，G(x)=(g1(x),...,gm(x))T，得到偏場的向量表達形式：

將(3)、(5)代入(2)中得到能量的表達：

根據文獻[11]中的能量優化過程，得到能量的變化形式：

2.3 能量優化

2.2 部分提出的能量函數實際是噪聲的表達。當噪聲取得最小值時，可以得到最接近真實前列腺磁共振圖像的校正圖像。由于能量函數是凸函數，所以通過一種類似梯度下降法的迭代方法便可使能量達到全局最小。

對于變量c、u、a分別固定其中兩個，然后讓能量函數對另外一個求導，重復多次，直到收斂。計算得到變量的迭代式如下：

3 實施

前列腺磁共振圖像偏移場校正的流程圖見圖1，詳細步驟將在后續部分介紹。

3.1 預處理

采用Pual等[13]開發的ITK-SNAP軟件(www. itksnap.org)對前列腺磁共振圖像進行分割獲得感興趣組織。然后通過一個直方圖線性拉伸，使感興趣組織獲得更好的視覺效果，預處理結果見圖2。

圖1 前列腺磁共振圖像偏場校正過程Fig. 1 Bias feld correction of prostate MR image.

圖2 預處理過程。A：原始圖像；B：分割結果；C：直方圖拉伸結果Fig. 2 Preprocessing. A: Original image; B: Segmentation results; C: Results of histogram stretch.

3.2 自動獲取組織的類別參數

本研究筆者直接對感興趣組織的直方圖進行處理。先對直方圖進行均值濾波，消除波形中的不規則點(灰度值出現跳變的點)，平滑波形。然后對波形進行峰值檢測[14]，根據檢測策略判斷圖像直方圖峰值的個數，由于圖像中不同的組織會呈現不同的灰度直方圖分布，從而自動獲得組織的類別參數，波形平滑過程如圖3所示。

3.3 能量優化

能量優化過程是整個校正過程的核心。首先給定能量函數的初始化參數c、μ；然后分別對c、μ、a進行迭代運算，使能量函數沿著梯度下降的方向變化；接著檢查能量函數的值，通過比較判斷是否達到收斂；根據判斷結果再次執行迭代運算(不收斂)或者評估偏場(收斂)、校正輸出。

3.4 加權指數的選擇

據文獻[14]中描述，加權指數是模糊聚類算法中一個重要的參數，它的選取決定著偏場矯正過程中能量優化的速度和程度。筆者通過圖4、5所示的兩組對比實驗來確定指數q的值。

圖3 波形平滑過程。A為原圖；B為濾波后的圖形Fig. 3 Procedure of waveform smoothing. A: Original outline of ROI; B: Outline after ftting.

圖4 不同圖像在不同q下的優化時間 圖5 能量函數在不同q值下的收斂程度Fig. 4 Show the time of four group different MR image in different exponent q. Fig. 5 Show the different degree of convergence in different exponent q.

通過以上兩組對比實驗可以得出以下結論：(1)當q值取整數時能量函數的迭代速度更快；(2)隨著q值的增大，能量函數能夠達到更好的收斂程度。另外結合文獻[15]中對于最佳q值的理論判斷，筆者使用q=2作為前列腺磁共振圖像偏場校正算法中的加權指數。

圖6 定性實驗結果。第一行是截取的橫斷面圖像，第二行是預處理后的結果，第三行是偏場校正后的結果Fig. 6 Qualitative experiment result. The three row are original prostate axial image, results of preprocessing, results of bias feld correction.

圖7 定量實驗結果。左、中、右分別為在rms、JS、BD下的對比實驗結果Fig. 7 Quantitative experiment result. The left, middle, and right columns show the result on index rms, JS, BD respectively

表1 本文方法與Li等[12]的方法處理結果的CV值Tab.1 Contrast of CV between ours and Li.et al’s method

4 實驗

筆者將分別展示本研究的方法在定性和定量實驗下的結果，所有的算法程序都運行在2.53GHzCPU，4G內存的 Matlab2012a平臺上。

4.1 定性評估

筆者選取6幅來自不同個體的前列腺磁共振掃描數據中的橫斷面圖像進行實驗，定性的實驗結果見圖6。

4.2 定量評估

為了定量評估算法的有效性，筆者將本研究的方法與Li等[12]的方法在一些公認指標下進行對比。變異系數(coefficient of variation，CV)(式11)被定義為標準差與期望的比值，它反映了組織(T)的相對離散程度。通常認為真實圖像的同種組織具有一致性，所以一個好的處理結果會得到一個小的CV值。本研究的方法與Li等[12]的方法對6幅前列腺橫斷面圖像處理的結果見表1。

均方根(root mean square，rms)(式12)反映了處理結果與真實圖像的差異。較小的rms值反映較好的處理結果。由于rms指標需要無偏場圖像作為參考，所以筆者選取來自模擬腦數據庫(http:// brainweb.bic.mni.mcgill.ca/brainweb/)的10幅無噪聲、無不均勻現象的腦部磁共振圖像作為參考圖像，將對應的10幅添加了3%噪聲、20%偏移場的圖像作為待處理圖像進行實驗。

JS (jaccard similarity)(式13)反映了處理結果與標準圖像的相似程度。S1、S2分別代表算法處理結果與參考圖像，較大的JS值反映了更好的處理結果。

巴氏距離(bhattacharyya distance)(式14)被認為是測量直方圖相似性效果最好的度量工具。筆者通過測量校正圖像p(x)與參考圖像q(x)的直方圖之間的巴氏距離來衡量算法的處理效果。據算法定義，較高的巴氏距離表示了更好的處理效果。

本研究的算法與Li等[12]的算法在各項指標下的實驗情況見圖7。實驗結果表明，筆者的方法相比Li等[12]的方法能取得更好的校正效果。

5 總結

本文提出了一個針對前列腺磁共振圖像的偏移場校正方法。一系列實驗和峰值檢測算法被用來自動獲取能量優化過程的初始參數，使得算法具有更好的準確性和穩定性。一組由三角函數和多項式函數組成的基函數保證了擬合偏場的平滑性。最后通過定性實驗和與相關算法在定量指標下的對比，驗證了本文提出的算法的有效性。

[References]

[1] Sun YH. The research status of prostate cancer in china. Chin J Urol, 2004, 25(2): 77- 80.

孫穎浩. 我國前列腺癌的研究現狀. 中國泌尿外科雜志, 2004, 25(2): 77-80.

[2] Hu MB, Jiang HW, Ding Q. Recent progress in early detection of prostate cancer. Shanghai Medical & Pharmaceutical Journal, 2013, 34(20): 12-15.

胡夢博, 姜昊文, 丁強. 前列腺癌早期診斷新進展. 上海醫藥, 2013, 34(20): 12-15.

[3] Wang JY, Chen X, Liu M, et al. MRI in clinical diagnosis and treatment of prostate cancer. Chin J Magn Reson Imaging, 2010, 1(4): 253-256.

王建業, 陳鑫, 劉明, 等. MRI在前列腺癌診斷和治療中的應用. 磁共振成像, 2010, 1(4): 253-256.

[4] Li BS, Wang L, Deng M, et al. The correlation between multiparametric MRI of prostate imaging reporting and data system score and transrectal ultrasound guided needle biopsy. Chin J Magn Reson Imaging, 2016, 7(5): 321-326.

李拔森, 王良, 鄧明, 等. 多參數MRI前列腺影像報告和數據系統評分與經直腸超聲引導下穿刺病理的相關性分析. 磁共振成像, 2016, 7(5): 321-326.

[5] Geert L, Oscar D, Jelle B, et al. Computer-Aided Detection of Prostate Cancer in MRI. IEEE Trans Med Image, 2014, 33(5): 1083-1092.

[6] Zhi DB, Qiu BS. A study on the best fitting model for diffusionweighted magnetic resonance imaging of normal prostate tissue at different b-values. Chin J Magn Reson Imaging, 2015, 6(8): 631-635.

智德波, 邱本勝. 正常前列腺組織磁共振彌散加權成像在不同b值下的最優擬合模型研究. 磁共振成像, 2015, 6(8): 631-635.

[7] Zhang D, Wu T, Shi CZ, et al. The reproducibility of IVIM-DWI in normal prostate and age-related analysis. Chin J Magn Reson Imaging, 2015, 6(11): 848-854.

張冬, 吳婷, 史長征, 等. 正常前列腺組織IVIM-DWI的可重現性研究及年齡相關性分析. 磁共振成像, 2015, 6(11): 848-854.

[8] Andrew S, Paul ST, Gareth JB, et al. Sources of intensity nonuniformity in spin echo images at 1.5 T. Magn Reson Med, 1994, 32(1): 121-128.

[9] Uros V, Franjo P, Bostjan L. A review of methods for correction of intensity inhomogeneity in MRI. IEEE Trans Med Imag, 2007, 17(3): 405-421.

[10] Johnston B, Atkins MS, Mackiewich B, et al. Segmentation of multiple sclerosis lesions in intensity corrected multispec-tral MRI. IEEE Trans Med Imag, 1996, 15(2): 154-169.

[11] Milles J, Zhu YM, Gimenez G, et al. MRI intensity nonuniformity correction using simultaneously spatial and gray-level histogram information. Medical imaging, 2004, 31(2): 81-90.

[12] Li CM, John CG, Christos D. Multiplicative intrinsic component optimization (MICO) for MRI bias field estimation and tissue segmentation. Magnetic Resonance Imaging, 2014, 32(7): 913-923.

[13] Paul AY, Joseph P, Heather CH, et al. User-guided 3D active contour segmentation of anatomical structures: Significantly improved effciency and reliability. Neuroimage, 2006, 31(3): 1116-1128.

[14] Wang HD, Liang D, Yang HX. Two algorithms for auto thresholding selection. J Nanjing Univ Posts Telecommun, 2002, 22(4): 85-88.

王厚大, 梁棟, 楊恒新. 自動閾值選取的兩種算法. 南京郵電學院學報, 2002, 22(4): 85-88.

[15] Gao XB, Pei JH, Xie WX. A study of weighting exponent m in a fuzzy c-means algorithm. Acta Electron Sin, 2000, 28(4): 1-4.

高新波, 裴繼紅, 謝維信. 模糊均值聚類算法中加權指數的研究.電子學報, 2000, 28(4): 1-4.

Smooth ftting of bias feld in prostate MRI with peak detection

YANG Xiong1, ZHAN Shu1*, XIE Dong-dong21School of Computer and Information, Hefei University of Technology, Hefei 230009, China
2Second Affliated Hospital, Anhui Medical University, Hefei 230601, China

Objective: To study correction of the inhomogeneity of grayscale (Bias Field) in prostate MR image. Materials and Methods: Several transverse images derived from magnetic resonance scanning data of prostate. The piecewise constant property of the real image and the smooth change characteristic of the bias feld are expressed in the image model. An energy function is constructed and the bias field estimation and tissue segmentation are realized by minimizing the energy function. The initial parameters of the energy function are obtained automatically by using the peak detection technique, and the smoothing ftting of the offset feld is realized by using a set of basis functions combined with trigonometric functions and polynomial functions. Results: Some qualitative evaluations showed the signifcant improvement of prostate MR image with severe intensity inhomogeneity by using our method. The comparison with other methods in some quantitative evaluation indexes (Coeffcient of variation, Root mean square and Jaccard similarity) is shown to demonstrate the better result of our method. Conclusion: Peak detection based bias correction method can perfect the intensity inhomogeneity in prostate MR image．

Magnetic resonance imaging; Prostate; Peak detection; Bias field correction; Smooth ftting

國家自然科學基金項目(編號：61371156)

1.合肥工業大學計算機與信息學院，合肥 230009

2.安徽醫科大學第二附屬醫院，合肥230601

詹曙，E-mail：shu_zhan@hfut.edu.cn

2016-08-01

接受日期：2016-09-23

R445.2；R737.25

A

10.12015/issn.1674-8034.2016.10.011楊雄, 詹曙, 謝棟棟. 結合峰值檢測的前列腺MRI偏移場平滑擬合. 磁共振成像, 2016, 7(10): 775-779.

*Correspondence to: Zhan S, E-mail: shu_zhan@hfut.edu.cn

Received 1 Aug 2016, Accepted 23 Sep 2016

ACKNOWLEDGEMENTS This work was supported by National Natural Science Foundation of China (No. 61371156).