楊雄,詹曙*,謝棟棟
結(jié)合峰值檢測(cè)的前列腺磁共振成像偏移場(chǎng)平滑擬合
楊雄1,詹曙1*,謝棟棟2
目的 研究前列腺磁共振圖像中灰度不均勻現(xiàn)象(偏場(chǎng))的校正方法。材料與方法 從幾組前列腺磁共振掃描數(shù)據(jù)中截取的橫斷面圖像。將真實(shí)圖像的分段常量特性與偏移場(chǎng)的平滑變化特性表達(dá)到圖像模型中,構(gòu)造一個(gè)能量函數(shù),通過(guò)能量函數(shù)的最小化實(shí)現(xiàn)偏移場(chǎng)評(píng)估和組織分割。利用峰值檢測(cè)技術(shù)自動(dòng)獲得能量函數(shù)的初始化參數(shù),并用結(jié)合三角函數(shù)與多項(xiàng)式函數(shù)的一組基函數(shù)實(shí)現(xiàn)對(duì)偏移場(chǎng)的平滑擬合。結(jié)果 定性的實(shí)驗(yàn)表明筆者的方法能對(duì)前列腺磁共振圖像中的偏場(chǎng)現(xiàn)象進(jìn)行有效的校正。另外通過(guò)與其他方法在變化系數(shù)、均方根、Jaccard相似度等指標(biāo)下的定量對(duì)比發(fā)現(xiàn),筆者的方法有更好的校正結(jié)果。結(jié)論 結(jié)合峰值檢測(cè)的偏移場(chǎng)校正方法能對(duì)前列腺磁共振圖像中的灰度不均勻現(xiàn)象有效改善。
磁共振成像;前列腺;峰值檢測(cè);偏場(chǎng)校正;平滑擬合
前列腺癌已成為嚴(yán)重危害老年男性健康的疾病,而且據(jù)相關(guān)報(bào)道,中國(guó)前列腺癌發(fā)病率正呈上升趨勢(shì)[1-2]。前列腺癌的早期診斷和檢測(cè)對(duì)疾病的治療預(yù)防非常重要??紤]到傳統(tǒng)前列腺癌早期檢測(cè)方法的局限性,一些借助于MRI技術(shù)的人工檢測(cè)方法,以及結(jié)合前列腺癌的病理學(xué)知識(shí)與MR影像特征的計(jì)算機(jī)輔助診斷方案也在研究之中[3-7]。
這些借助于MR影像的診斷方法對(duì)于影像的質(zhì)量有較高要求,然而醫(yī)學(xué)圖像中普遍存在著一種灰度不均勻現(xiàn)象,也即偏場(chǎng)。偏移場(chǎng)的存在會(huì)影響病灶區(qū)的視覺(jué)效果以及計(jì)算機(jī)對(duì)病灶組織的判斷。根據(jù)文獻(xiàn)[8-9],偏移場(chǎng)的產(chǎn)生是由于成像設(shè)備和患者特異性的綜合作用。目前對(duì)偏移場(chǎng)的校正方法大致可以分為前瞻性方法和回顧性方法[9]。前瞻性方法是從MR圖像獲得過(guò)程的角度來(lái)考慮校正方法,這類方法雖然也能取得一些效果,但是往往不能解決由于患者特異性所帶來(lái)的影響?;仡櫺苑椒ㄖ灰蕾囉讷@得圖像的信息,因而能更全面地解決偏場(chǎng),當(dāng)然也得到了更多的關(guān)注和研究。Johnston等[10]最先提出一種基于同態(tài)濾波的方法,這個(gè)方法假設(shè)偏移場(chǎng)具有低頻特性,但是圖像本身也具有很多低頻成分,因此高通濾波器也不能取得理想的結(jié)果。Miles等[11]提出一個(gè)結(jié)合直方圖信息和B樣條曲面擬合的方法。Li等[12]在乘法固有成分優(yōu)化的思想下,提出一個(gè)能量最小化模型,同時(shí)實(shí)現(xiàn)偏場(chǎng)矯正和多組織分割。
目前對(duì)于偏場(chǎng)矯正的方法很多,本文在已有算法的基礎(chǔ)上作了一些改進(jìn),提出一種針對(duì)前列腺組織MR圖像的偏場(chǎng)校正方法。首先根據(jù)前列腺磁共振圖像的灰度特征進(jìn)行自適應(yīng)的預(yù)處理,然后利用峰值檢測(cè)技術(shù)從圖像的直方圖信息中自動(dòng)獲得前列腺組織的類別參數(shù),接著根據(jù)圖像模型構(gòu)造能量最小化模型,最后提出一種由最簡(jiǎn)多項(xiàng)式函數(shù)與三角函數(shù)組合而成的曲面擬合方法,實(shí)現(xiàn)對(duì)偏場(chǎng)的平滑擬合。
2.1 MR圖像模型
目前一種普遍認(rèn)可的磁共振圖像模型如式(1)所示:
其中x代表圖像像素點(diǎn),v(x)表示獲得的圖像,u(x)代表沒(méi)有偏場(chǎng)與噪聲作用的真實(shí)圖像,b(x)是偏場(chǎng),n(x)則指圖像噪聲。通常認(rèn)為圖像中的背景與前景中的各個(gè)目標(biāo)在灰度上具有一致性,而圖像中的灰度不均勻變化是一個(gè)平滑過(guò)程。所以在這個(gè)圖像模型中,筆者認(rèn)為代表真實(shí)圖像的u(x)具有分段常量特性,即圖像中不同的組織可以近似為不同的灰度常量區(qū)域,而表征圖像不均勻現(xiàn)象的b(x)具有平滑特性。
2.2 能量函數(shù)
針對(duì)u(x)的分段常量特性,假設(shè)圖像中存在n個(gè)不同的組織,每個(gè)組織的灰度值為常量ci(i=1,...,n)。隸屬度函數(shù)μi(x)表示像素點(diǎn)x屬于第i類組織的概率。為了更真實(shí)地反映實(shí)像,這里參照模糊C均值(fuzzy c-means algorithm, FCM)算法采用一種模糊表示方法()。那么實(shí)像u(x)表示如下:
其中q是加權(quán)指數(shù)。對(duì)于偏場(chǎng),筆者用一系列連續(xù)函數(shù)的加權(quán)和表示:
m的取值決定著曲面的擬合精度。通常對(duì)于偏場(chǎng)的擬合會(huì)用多項(xiàng)式函數(shù)作為連續(xù)函數(shù),但是筆者發(fā)現(xiàn)三角函數(shù)有著等同高次函數(shù)的平滑特性,而且計(jì)算更簡(jiǎn)。所以用低次最簡(jiǎn)多項(xiàng)式函數(shù)與三角函數(shù)共同作為擬合曲面的基函數(shù)。
令a=(a1,...,am)T,G(x)=(g1(x),...,gm(x))T,得到偏場(chǎng)的向量表達(dá)形式:
將(3)、(5)代入(2)中得到能量的表達(dá):
根據(jù)文獻(xiàn)[11]中的能量?jī)?yōu)化過(guò)程,得到能量的變化形式:
2.3 能量?jī)?yōu)化
2.2 部分提出的能量函數(shù)實(shí)際是噪聲的表達(dá)。當(dāng)噪聲取得最小值時(shí),可以得到最接近真實(shí)前列腺磁共振圖像的校正圖像。由于能量函數(shù)是凸函數(shù),所以通過(guò)一種類似梯度下降法的迭代方法便可使能量達(dá)到全局最小。
對(duì)于變量c、u、a分別固定其中兩個(gè),然后讓能量函數(shù)對(duì)另外一個(gè)求導(dǎo),重復(fù)多次,直到收斂。計(jì)算得到變量的迭代式如下:
前列腺磁共振圖像偏移場(chǎng)校正的流程圖見(jiàn)圖1,詳細(xì)步驟將在后續(xù)部分介紹。
3.1 預(yù)處理
采用Pual等[13]開(kāi)發(fā)的ITK-SNAP軟件(www. itksnap.org)對(duì)前列腺磁共振圖像進(jìn)行分割獲得感興趣組織。然后通過(guò)一個(gè)直方圖線性拉伸,使感興趣組織獲得更好的視覺(jué)效果,預(yù)處理結(jié)果見(jiàn)圖2。
圖1 前列腺磁共振圖像偏場(chǎng)校正過(guò)程Fig. 1 Bias feld correction of prostate MR image.
圖2 預(yù)處理過(guò)程。A:原始圖像;B:分割結(jié)果;C:直方圖拉伸結(jié)果Fig. 2 Preprocessing. A: Original image; B: Segmentation results; C: Results of histogram stretch.
3.2 自動(dòng)獲取組織的類別參數(shù)
本研究筆者直接對(duì)感興趣組織的直方圖進(jìn)行處理。先對(duì)直方圖進(jìn)行均值濾波,消除波形中的不規(guī)則點(diǎn)(灰度值出現(xiàn)跳變的點(diǎn)),平滑波形。然后對(duì)波形進(jìn)行峰值檢測(cè)[14],根據(jù)檢測(cè)策略判斷圖像直方圖峰值的個(gè)數(shù),由于圖像中不同的組織會(huì)呈現(xiàn)不同的灰度直方圖分布,從而自動(dòng)獲得組織的類別參數(shù),波形平滑過(guò)程如圖3所示。
3.3 能量?jī)?yōu)化
能量?jī)?yōu)化過(guò)程是整個(gè)校正過(guò)程的核心。首先給定能量函數(shù)的初始化參數(shù)c、μ;然后分別對(duì)c、μ、a進(jìn)行迭代運(yùn)算,使能量函數(shù)沿著梯度下降的方向變化;接著檢查能量函數(shù)的值,通過(guò)比較判斷是否達(dá)到收斂;根據(jù)判斷結(jié)果再次執(zhí)行迭代運(yùn)算(不收斂)或者評(píng)估偏場(chǎng)(收斂)、校正輸出。
3.4 加權(quán)指數(shù)的選擇
據(jù)文獻(xiàn)[14]中描述,加權(quán)指數(shù)是模糊聚類算法中一個(gè)重要的參數(shù),它的選取決定著偏場(chǎng)矯正過(guò)程中能量?jī)?yōu)化的速度和程度。筆者通過(guò)圖4、5所示的兩組對(duì)比實(shí)驗(yàn)來(lái)確定指數(shù)q的值。
圖3 波形平滑過(guò)程。A為原圖;B為濾波后的圖形Fig. 3 Procedure of waveform smoothing. A: Original outline of ROI; B: Outline after ftting.
圖4 不同圖像在不同q下的優(yōu)化時(shí)間 圖5 能量函數(shù)在不同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.
通過(guò)以上兩組對(duì)比實(shí)驗(yàn)可以得出以下結(jié)論:(1)當(dāng)q值取整數(shù)時(shí)能量函數(shù)的迭代速度更快;(2)隨著q值的增大,能量函數(shù)能夠達(dá)到更好的收斂程度。另外結(jié)合文獻(xiàn)[15]中對(duì)于最佳q值的理論判斷,筆者使用q=2作為前列腺磁共振圖像偏場(chǎng)校正算法中的加權(quán)指數(shù)。
圖6 定性實(shí)驗(yàn)結(jié)果。第一行是截取的橫斷面圖像,第二行是預(yù)處理后的結(jié)果,第三行是偏場(chǎng)校正后的結(jié)果Fig. 6 Qualitative experiment result. The three row are original prostate axial image, results of preprocessing, results of bias feld correction.
圖7 定量實(shí)驗(yàn)結(jié)果。左、中、右分別為在rms、JS、BD下的對(duì)比實(shí)驗(yàn)結(jié)果Fig. 7 Quantitative experiment result. The left, middle, and right columns show the result on index rms, JS, BD respectively
表1 本文方法與Li等[12]的方法處理結(jié)果的CV值Tab.1 Contrast of CV between ours and Li.et al’s method
筆者將分別展示本研究的方法在定性和定量實(shí)驗(yàn)下的結(jié)果,所有的算法程序都運(yùn)行在2.53GHzCPU,4G內(nèi)存的 Matlab2012a平臺(tái)上。
4.1 定性評(píng)估
筆者選取6幅來(lái)自不同個(gè)體的前列腺磁共振掃描數(shù)據(jù)中的橫斷面圖像進(jìn)行實(shí)驗(yàn),定性的實(shí)驗(yàn)結(jié)果見(jiàn)圖6。
4.2 定量評(píng)估
為了定量評(píng)估算法的有效性,筆者將本研究的方法與Li等[12]的方法在一些公認(rèn)指標(biāo)下進(jìn)行對(duì)比。變異系數(shù)(coefficient of variation,CV)(式11)被定義為標(biāo)準(zhǔn)差與期望的比值,它反映了組織(T)的相對(duì)離散程度。通常認(rèn)為真實(shí)圖像的同種組織具有一致性,所以一個(gè)好的處理結(jié)果會(huì)得到一個(gè)小的CV值。本研究的方法與Li等[12]的方法對(duì)6幅前列腺橫斷面圖像處理的結(jié)果見(jiàn)表1。
均方根(root mean square,rms)(式12)反映了處理結(jié)果與真實(shí)圖像的差異。較小的rms值反映較好的處理結(jié)果。由于rms指標(biāo)需要無(wú)偏場(chǎng)圖像作為參考,所以筆者選取來(lái)自模擬腦數(shù)據(jù)庫(kù)(http:// brainweb.bic.mni.mcgill.ca/brainweb/)的10幅無(wú)噪聲、無(wú)不均勻現(xiàn)象的腦部磁共振圖像作為參考圖像,將對(duì)應(yīng)的10幅添加了3%噪聲、20%偏移場(chǎng)的圖像作為待處理圖像進(jìn)行實(shí)驗(yàn)。
JS (jaccard similarity)(式13)反映了處理結(jié)果與標(biāo)準(zhǔn)圖像的相似程度。S1、S2分別代表算法處理結(jié)果與參考圖像,較大的JS值反映了更好的處理結(jié)果。
巴氏距離(bhattacharyya distance)(式14)被認(rèn)為是測(cè)量直方圖相似性效果最好的度量工具。筆者通過(guò)測(cè)量校正圖像p(x)與參考圖像q(x)的直方圖之間的巴氏距離來(lái)衡量算法的處理效果。據(jù)算法定義,較高的巴氏距離表示了更好的處理效果。
本研究的算法與Li等[12]的算法在各項(xiàng)指標(biāo)下的實(shí)驗(yàn)情況見(jiàn)圖7。實(shí)驗(yàn)結(jié)果表明,筆者的方法相比Li等[12]的方法能取得更好的校正效果。
本文提出了一個(gè)針對(duì)前列腺磁共振圖像的偏移場(chǎng)校正方法。一系列實(shí)驗(yàn)和峰值檢測(cè)算法被用來(lái)自動(dòng)獲取能量?jī)?yōu)化過(guò)程的初始參數(shù),使得算法具有更好的準(zhǔn)確性和穩(wěn)定性。一組由三角函數(shù)和多項(xiàng)式函數(shù)組成的基函數(shù)保證了擬合偏場(chǎng)的平滑性。最后通過(guò)定性實(shí)驗(yàn)和與相關(guān)算法在定量指標(biāo)下的對(duì)比,驗(yàn)證了本文提出的算法的有效性。
[References]
[1] Sun YH. The research status of prostate cancer in china. Chin J Urol, 2004, 25(2): 77- 80.
孫穎浩. 我國(guó)前列腺癌的研究現(xiàn)狀. 中國(guó)泌尿外科雜志, 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.
胡夢(mèng)博, 姜昊文, 丁強(qiáng). 前列腺癌早期診斷新進(jìn)展. 上海醫(yī)藥, 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.
王建業(yè), 陳鑫, 劉明, 等. MRI在前列腺癌診斷和治療中的應(yīng)用. 磁共振成像, 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.
李拔森, 王良, 鄧明, 等. 多參數(shù)MRI前列腺影像報(bào)告和數(shù)據(jù)系統(tǒng)評(píng)分與經(jīng)直腸超聲引導(dǎo)下穿刺病理的相關(guān)性分析. 磁共振成像, 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.
智德波, 邱本勝. 正常前列腺組織磁共振彌散加權(quán)成像在不同b值下的最優(yōu)擬合模型研究. 磁共振成像, 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.
張冬, 吳婷, 史長(zhǎng)征, 等. 正常前列腺組織IVIM-DWI的可重現(xiàn)性研究及年齡相關(guān)性分析. 磁共振成像, 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.
王厚大, 梁棟, 楊恒新. 自動(dòng)閾值選取的兩種算法. 南京郵電學(xué)院學(xué)報(bào), 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.
高新波, 裴繼紅, 謝維信. 模糊均值聚類算法中加權(quán)指數(shù)的研究.電子學(xué)報(bào), 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
國(guó)家自然科學(xué)基金項(xiàng)目(編號(hào):61371156)
1.合肥工業(yè)大學(xué)計(jì)算機(jī)與信息學(xué)院,合肥 230009
2.安徽醫(yī)科大學(xué)第二附屬醫(yī)院,合肥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楊雄, 詹曙, 謝棟棟. 結(jié)合峰值檢測(cè)的前列腺M(fèi)RI偏移場(chǎng)平滑擬合. 磁共振成像, 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).