趙九龍，馬 瑜，李 爽，孟亞州，白 冰

（1.寧夏大學(xué) 研究生院，寧夏 銀川 750021；2.上海交通大學(xué) 電子與電氣工程學(xué)院 圖像處理與模式識(shí)別研究所，上海 200240）

1 引 言

計(jì)算機(jī)斷層掃描 （Computerized Tomography，CT）和核磁共振（Nuclear Magnetic Resonance，NMR）是目前獲取醫(yī)學(xué)圖像二維切片的主要技術(shù)手段，但二維圖像有很多局限性，不利于醫(yī)生及時(shí)有效地分析病情，因此近年來，三維醫(yī)學(xué)圖像的重構(gòu)可視化技術(shù)方興未艾［1－2］。利用醫(yī)學(xué)二維切片重構(gòu)出的三維醫(yī)學(xué)圖像可以清晰、準(zhǔn)確地顯示病患處的具體輪廓和細(xì)節(jié)信息，為醫(yī)生提供更逼真的顯示手段以及定量和定性分析工具，極大地提高了醫(yī)學(xué)分析的準(zhǔn)確性和有效性。但由于二維醫(yī)學(xué)圖像在采集和傳輸?shù)倪^程中難免會(huì)被噪聲污染，較常見的兩種噪聲是高斯噪聲和椒鹽噪聲。產(chǎn)生的主要原因是：醫(yī)學(xué)成像設(shè)備有很多是利用光學(xué)成像，如CT機(jī)為光學(xué)成像設(shè)備，其圖像傳感器的光電轉(zhuǎn)換過程中由于原件靈敏度的不均性，會(huì)產(chǎn)生高斯噪聲，數(shù)字化過程的量化誤差及人為因素干擾同樣會(huì)產(chǎn)生高斯噪聲；醫(yī)學(xué)成像設(shè)備的傳輸信道、解碼處理等過程會(huì)產(chǎn)生椒鹽噪聲。

圖像去噪是圖像預(yù)處理的重要部分，對(duì)含噪聲圖像有效地去噪［3］有利于后續(xù)對(duì)圖像特征的提取和分析等。若從醫(yī)學(xué)成像設(shè)備獲取的切片含噪聲且不能較好地去噪，很大程度上會(huì)影響三維醫(yī)學(xué)圖像的重構(gòu)精度和可視化質(zhì)量，因此能否有效地去除圖像在設(shè)備采集及傳輸中產(chǎn)生的噪聲直接影響三維重構(gòu)和可視化的效果，進(jìn)而會(huì)影響對(duì)病情分析的準(zhǔn)確性和高效性。

目前多數(shù)去噪算法只針對(duì)某一種特定噪聲且都應(yīng)用在二維圖像，本文將較新穎的分?jǐn)?shù)階積分去噪算法和經(jīng)典的非線性中值濾波算法結(jié)合，去除三維醫(yī)學(xué)圖像的混合噪聲。

2 典型去噪方法比較

（1）均值濾波。均值濾波算法是線性濾波算法的代表，其主要原理是用以待處理像素點(diǎn)為中心的若干鄰域像素的灰度算數(shù)平均值代替待處理像素點(diǎn)灰度值。其對(duì)高斯噪聲有不錯(cuò)的效果［4］，但是不適合去除突變的噪聲，對(duì)含細(xì)節(jié)較多的圖像失真較大。

（2）中值濾波。中值濾波有較長(zhǎng)的歷史，是一種非線性信號(hào)處理方法。其主要思想是基于排序理論，對(duì)待處理像素點(diǎn)的某一鄰域像素按灰度值排序，取排序后中間值替換當(dāng)前待處理像素。中值濾波可以一定程度地解決線性濾波后圖像細(xì)節(jié)失真的問題，常用于濾除幅值變化大的椒鹽噪聲，但需要選擇合適大小的濾波窗口，濾波窗過大會(huì)過多地破壞圖像信息，過小不能去除噪聲，文獻(xiàn)［5］提出了自適應(yīng)的中值濾波算法。但對(duì)于包含細(xì)節(jié)信息較多的噪聲圖像以及每一點(diǎn)均被污染的含高斯噪聲圖像，中值濾波效果不佳。

（3）高斯去噪。文獻(xiàn)［6］中提出了基于高斯去噪的三維邊緣曲面追蹤算法，將高斯去噪算法應(yīng)用到三維圖像重構(gòu)中，一定程度地解決了三維邊緣曲面追蹤算法中圖像受高斯噪聲污染的問題。

（4）分?jǐn)?shù)階積分去噪。分?jǐn)?shù)階微積分由周激流、蒲亦非等人引入數(shù)字圖像處理領(lǐng)域［7－9］，并驗(yàn)證了分?jǐn)?shù)階積分較很多傳統(tǒng)去噪算法有更好的結(jié)果［7－9］。文獻(xiàn)［7］和［9］已證明，由于分?jǐn)?shù)階積分可以削弱高頻信號(hào)并極大的非線性保留低頻信號(hào)，所以分?jǐn)?shù)階積分可以在去除高頻的噪聲同時(shí)保留低頻的圖像本身信息，因而其去噪效果較好。文獻(xiàn)［10］將分?jǐn)?shù)階積分理論推廣應(yīng)用至三維，并驗(yàn)證了基于分?jǐn)?shù)階積分去噪的邊緣曲面追蹤算法較文獻(xiàn)［6］有更好的效果。

以上去噪算法均針對(duì)某一種特定噪聲，并未考慮混合噪聲干擾，本文結(jié)合分?jǐn)?shù)階積分去噪和中值濾波的優(yōu)點(diǎn)，將兩種去噪算法結(jié)合應(yīng)用于三維醫(yī)學(xué)圖像混合噪聲去除。

3 混合去噪的邊緣曲面追蹤算法

3.1 三維分?jǐn)?shù)階積分定義

基于Riemann－Liouville三維分?jǐn)?shù)階積分的離散定義為［9－10］：

基于二維分?jǐn)?shù)階積分離散模板［9］和三維分?jǐn)?shù)階積分定義，結(jié)合三維圖像的空間信息，構(gòu)造三維分?jǐn)?shù)階積分離散模板如圖1所示。

圖1 三維分?jǐn)?shù)階離散模板Fig.1 3Dfractional discrete template

圖中有3層切片，分別為z層及其上方的z＋1層和下方的z－1層。如圖所示我們選取與中心體素（centre voxel）相鄰的18個(gè)方向體素為三維離散模板的有用體素。18個(gè)方向上濾波結(jié)果如下：

其中：VNε（i，j，k）為18個(gè)方向中某一方向的最終值，ε＝1，2，…18。fm（i，j，k）為從模板中心體素開始某一方向的第m 個(gè)體素值，m＝0，1，2，n。fe（m）為模板中濾波系數(shù)，由三維分?jǐn)?shù)階積分定義的離散表達(dá)式（1）～（3）可知取值為：

3.2 分?jǐn)?shù)階積分階數(shù)的自適應(yīng)選擇

考慮噪聲的高頻特性和局部灰度突變的特性，定義梯度判別因子選擇自適應(yīng)分?jǐn)?shù)階的階數(shù)。三維函數(shù)的梯度算子如下：

三維函數(shù)的梯度計(jì)算方式同樣為計(jì)算函數(shù)對(duì)不同分量的偏導(dǎo)數(shù)，求各偏導(dǎo)數(shù)的平方和后取算數(shù)平方根作為梯度模值。

在三維圖像中，各體素點(diǎn)最小距離為單位一，三維圖像的體素值即為離散的三維函數(shù)值，對(duì)于離散的三維函數(shù)采用雙向差分值的二分之一來計(jì)算梯度。如果三維圖像（信號(hào)）為I（i，j，k），則梯度的差分計(jì)算規(guī)則如下：

定義G（i，j）為當(dāng)前體素點(diǎn)的梯度值，NεG（i，j）為以當(dāng)前體素為中心的相鄰18體素的梯度值，ε＝1，…18。首先定義梯度判別因子J（i，j）為：

式中：min為取最小值操作，‖為取模值，由上式可知，梯度判別因子J（i，j）為當(dāng)前體素的梯度模值與其三維離散模板內(nèi)18個(gè)相鄰點(diǎn)梯度模值做差的最小值。則定義自適應(yīng)分?jǐn)?shù)階階數(shù)表達(dá)式為：

其中：α為調(diào)節(jié)系數(shù)，MaxG（i，j）為18體素點(diǎn)的梯度最大值。T為依據(jù)不同圖像選擇的閾值。當(dāng)所處理體素的梯度值與周邊體素的梯度值發(fā)生較大變化時(shí)，即梯度判別因子大于給定閾值時(shí)，認(rèn)定此體素為高頻噪聲并進(jìn)行自適應(yīng)去噪處理。

式（8）自適應(yīng)分?jǐn)?shù)階的階數(shù)以當(dāng)前待處理像素的梯度值為分子，以18鄰域內(nèi)梯度最大值為分母，當(dāng)以梯度判別條件判定當(dāng)前體素為噪聲時(shí)，分子即噪聲的梯度模值與噪聲幅值變化成正比，對(duì)于高斯噪聲或椒鹽噪聲，其噪聲點(diǎn)的梯度即v的分子都會(huì)較大，則自適應(yīng)的階數(shù)v也會(huì)隨之增大，若噪聲幅值較小，則計(jì)算后的梯度模值也相對(duì)較小，與之匹配的是較小的分?jǐn)?shù)階階數(shù)；因此，符合文獻(xiàn)［7］中闡述的分?jǐn)?shù)階的幅頻相應(yīng)特性，即分?jǐn)?shù)階階數(shù)越高，削弱高頻噪聲的能力越強(qiáng)，保留低頻圖像信息能力越強(qiáng)。當(dāng)梯度判別后不是噪聲的點(diǎn)，采取保留原值的操作，可以最大程度保留圖像原有信息。

由式（4）得對(duì)當(dāng)前體素點(diǎn)濾波后的新的I（i，j，k）為：

式中的I（i，j，k）即為經(jīng)過自適應(yīng)分?jǐn)?shù)階積分去噪后的新的體素。在求得18個(gè)方向的卷積值后取算數(shù)平均即為新的體素值。按此規(guī)則用三維分?jǐn)?shù)階離散模板遍歷整個(gè)三維帶噪圖像（數(shù)據(jù)），即完成一次自適應(yīng)分?jǐn)?shù)階積分去噪。

3.3 中值濾波

自適應(yīng)分?jǐn)?shù)階積分可以更好地去除混合噪聲中的高斯噪聲，對(duì)幅值突變明顯的椒鹽噪聲也有一定削弱作用，但是對(duì)灰度值變化較明顯的噪聲，中值濾波更為有效，在經(jīng)自適應(yīng)分?jǐn)?shù)階積分濾波后，再用中值濾波進(jìn)行進(jìn)一步去噪，去除分?jǐn)?shù)階積分濾波后幅值仍較大或較小的噪聲。

圖2 中值濾波模板Fig.2 Median filter template

圖2所示為中值濾波的3×3×3模板，以當(dāng)前待處理體素為中心，其周圍緊鄰的26個(gè)體素為影響體素，對(duì)此27個(gè)體素進(jìn)行排序，取中值作為濾波結(jié)果。文獻(xiàn)［5］中改進(jìn)的自適應(yīng)中值濾波模板最多取到超過11×11，考慮到三維數(shù)據(jù)量較大，過大的濾波模板會(huì)將更多的無關(guān)體素考慮在內(nèi)，且計(jì)算量急劇增長(zhǎng)，因此兼顧計(jì)算量，圖中3×3×3模板已包含27個(gè)體素，可以達(dá)到中值濾波的效果。

3.4 邊緣曲面追蹤算法

三維邊緣曲面追蹤算法數(shù)學(xué)模型為：

式中：f（x，y，z）為三維體數(shù)據(jù)，.2f（x，y，z）為三維函數(shù)f（x，y，z）的Laplacian算子，‖.f（x，y，z）‖ 為三維函數(shù)f（x，y，z）的梯度模值。在三維邊緣曲面重構(gòu)算法中Laplacian算子可檢測(cè)并抽取出三維邊緣曲面的零交叉等值面，梯度算子可用來追蹤到高亮度的邊緣曲面，其中梯度算子和Laplacian算子均為離散值。三維邊緣曲面追蹤算法的詳細(xì)描述見文獻(xiàn)［11］和［12］。

3.5 基于混合去噪的邊緣曲面追蹤算法

基于自適應(yīng)三維分?jǐn)?shù)階積分和中值去噪的邊緣曲面追蹤算法步驟如下：

（1）利用自適應(yīng)分?jǐn)?shù)階積分和中值濾波對(duì)噪聲污染的三維切片進(jìn)行去噪；

（2）利用2.4節(jié)中介紹的算法追蹤三維邊緣曲面；

（3）計(jì)算三維邊緣曲面片并將其三角片化；

（4）通過計(jì)算機(jī)圖形學(xué)技術(shù)可視化三維邊緣曲面；

已被噪聲污染的三維圖像（數(shù)據(jù)）f（x，y，z）經(jīng)自適應(yīng)三維分?jǐn)?shù)階積分和中值去噪后變?yōu)镮（x，y，z），從去噪后的I（x，y，z）追蹤邊緣曲面，基于混合去噪的邊緣曲面追蹤算法如下：

4 實(shí)驗(yàn)結(jié)果與分析

實(shí)驗(yàn)以Visual C＋＋2010結(jié)合OpenGL為仿真環(huán)境。

為驗(yàn)證本文算法有效性，采用一組嬰兒頭部的CT圖像進(jìn)行驗(yàn)證。圖3為嬰兒頭部圖像的三維重建原圖以及不同去噪算法的處理結(jié)果。

如圖3所示，圖3（c）的均值去噪，圖3（d）的中值去噪和圖3（e）的經(jīng)分?jǐn)?shù)階階數(shù)為0.5的分?jǐn)?shù)階積分去噪對(duì)高斯噪聲和椒鹽噪聲的混合噪聲去噪效果均不是很好；圖3（c）均值去噪后仍有較明顯的黑色暗點(diǎn)塊和明顯的凸凹不平，圖3（d）中值去噪后整體仍有高斯噪聲引起的凸起等噪聲，圖3（e）為經(jīng)固定分?jǐn)?shù)階次為0.5的分?jǐn)?shù)階積分去噪結(jié)果，雖然取得了一定的去噪結(jié)果，但相比于本文算法仍有不足之處。圖（3）f為經(jīng)自適應(yīng)分?jǐn)?shù)階積分和中值結(jié)合去噪后的結(jié)果，相對(duì)于前幾種單一算法，明顯地減少了噪聲干擾。

圖3 嬰兒頭部圖像Fig.3 Denoising results of infant image

二維數(shù)字圖像的去噪結(jié)果常采用峰值信噪比（Peak Signal－to－Noise Ratio，PSNR）和均方誤差（Mean Square Error，MSE）來客觀評(píng)價(jià)，三維數(shù)據(jù)中我們定義并采用體數(shù)據(jù)峰值信噪比（Volumetric Peak Signal－to－Noise Ratio，VPSNR）和體數(shù)據(jù)均方誤差（Volumetric Mean Square Error，VMSE）［7，10，13］來客觀評(píng)價(jià)去噪效果，定義如下：

式中：m，n和l為三維圖像（數(shù)據(jù)）規(guī)模大小，I（i，j，k）為源圖像，I＊（i，j，k）為去噪后的圖像。顯然體數(shù)據(jù)峰值信噪比越大且體數(shù)據(jù)均方誤差越小說明去噪效果越好。

表1為不同去噪方法的體數(shù)據(jù)均方誤差和體數(shù)據(jù)峰值信噪比。由體數(shù)據(jù)均方誤差和體數(shù)據(jù)峰值信噪比的定義可知，三維圖像去噪后的體數(shù)據(jù)均方誤差越小，體數(shù)據(jù)峰值信噪比越大說明去噪效果越好。從表中可以看出，對(duì)于混合噪聲，傳統(tǒng)的去噪方法體數(shù)據(jù)均方誤差均較大，本文自適應(yīng)分?jǐn)?shù)階積分與中值結(jié)合的方法可以獲得更小的均方誤差和更大的峰值信噪比。

表1 不同去噪方法的體數(shù)據(jù)均方誤差和體數(shù)據(jù)峰值信噪比Tab.1 VMSE and VPSNR of different denoising methods

5 結(jié) 論

目前去噪算法很多，各有利弊。對(duì)于椒鹽噪聲，已經(jīng)出現(xiàn)了很多改進(jìn)的中值濾波算法，分?jǐn)?shù)階積分去噪歷史較短，但對(duì)于高斯噪聲去噪效果不錯(cuò)。對(duì)于高斯噪聲與椒鹽噪聲的混合噪聲，相關(guān)的去噪方法較少，本文提出了針對(duì)不同梯度信息的自適應(yīng)分?jǐn)?shù)階積分，并與中值濾波相結(jié)合，應(yīng)用于三維醫(yī)學(xué)圖像的混合噪聲去除。實(shí)驗(yàn)結(jié)果中本文算法對(duì)圖像去噪后的 VMSE降低至243.95，VPSNR提高到23.70，實(shí)驗(yàn)結(jié)果表明，本文方法可以取得相對(duì)不錯(cuò)的去噪效果。

本文將較新穎的分?jǐn)?shù)階積分和經(jīng)典的中值濾波結(jié)合用于混合噪聲去除。分?jǐn)?shù)階積分核心思想是加權(quán)平均的平滑濾波過程，中值對(duì)于灰度值極高或極低的椒鹽噪聲有良好的效果，本文處理后三維圖像雖有相對(duì)好的效果，但細(xì)節(jié)保留不夠好，整體略顯平滑。如何更好地利用分?jǐn)?shù)階積分和中值的特性，完全實(shí)現(xiàn)混合噪聲的自適應(yīng)去除，將是下一步研究的方向。

［1］ 武君勝，楊紅遠(yuǎn)，諶洪初.醫(yī)學(xué)TPS中放射線劑量分布的三維可視化方法［J］.計(jì)算機(jī)應(yīng)用，2010，30（3）：582－584.Wu J S，Yang H Y，Chen H C.3Dvisualization method of radiotherapy dose distribution in medical TPS ［J］.Journal of Computer Applications，2010，30（3）：582－584.（in Chinese）

［2］ 李祥，傅俊瓊.基于圖像融合的微表面快速三維重構(gòu)算法研究［J］.計(jì)算機(jī)應(yīng)用研究，2009，26（10）：3992－3994.Li X，F(xiàn)u J Q.3Dreconstruction algorithm for micro－surface based on image fusion ［J］.Application Research of Computer，2009，26（10）：3992－3994.（in Chinese）

［3］ 董雪，林志賢，郭太良.基于LoG算子改進(jìn)的自適應(yīng)閾值小波去噪算法［J］.液晶與顯示，2014，29（2）：275－280.Dong X，Lin Z X，Guo T L.Improved self－adaptive threshold wavelet denoising analysisbased on LoG operator［J］.Chinese Journal of Liquid Crystals and Displays，2014，29（2）：275－280.（in Chinese）

［4］ 陳明舉.基于統(tǒng)計(jì)特性的非局部均值去噪算法［J］.液晶與顯示，2014，29（3）：450－454.Chen M J.Non－local means image denoising algorithm based on statistical property［J］.Chinese Journal of Liquid Crystals and Displays，2014，29（3）：450－454.（in Chinese）

［5］ 郭海霞，解凱.一種改進(jìn)的自適應(yīng)中值濾波算法［J］.中國(guó)圖象圖形學(xué)報(bào)，2007，12（7）：1185－1188.Guo H X，Xie K.An improved method of adaptive median filter［J］.Journal of Image and Graphics，2007，12（7）：1185－1188.（in Chinese）

［6］ Ma Y，Zhang Y N，Wang Y G，et al.An improved denoised 3Dedge extraction operator within biomedical images［C］.Berling：CMSP2012，Springer Communications in Computer and Information Science，2012，346：309－316.

［7］ Hu J R，Pu Y F，Zhou J L.A novel image denoising algorithm based on Riemann－Liouville definition［J］.Journal of Computers，2011，6（7）：1332－1338.

［8］ 趙健.分?jǐn)?shù)階微分在圖像紋理增強(qiáng)中的應(yīng)用［J］.液晶與顯示，2012，27（1）：121－124.Zhao J.Fractional differential and its application in image texture enhancement［J］.Chinese Journal of Liquid Crystals and Displays，2012，27（1）：121－124.（in Chinese）

［9］ 劉彥，蒲亦非，周激流.一種基于分?jǐn)?shù)階積分的數(shù)字圖像去噪算法［J］.四川大學(xué)學(xué)報(bào)：工程科學(xué)版，2011，43（3）：90－95.Liu Y，Pu Y F，Zhou J L.A digital image denoising method based on fractional calculus［J］.Journal of Sichuan University：Engineering Science Edition，2011，43（3）：90－95.（in Chinese）

［10］ 馬瑜，張艷寧，王利生.噪音工業(yè)計(jì)算機(jī)斷層圖像中追蹤三維邊緣曲面［J］.光子學(xué)報(bào)，2013，42（8）：916－923.Ma Y，Zhang Y N，Wang L S.Three dimensional edge surfaces tracked from noisy industrial CT slice images［J］.Acta Photonica Sinica，2013，42（8）：916－923.（in Chinese）

［11］ Wang L S，Bai J，Ying K.Adaptive approximation of the boundary surface of a neuron in confocal microscopy volumetric images［J］.Medical and Biological Engineering and Computing，2003，41（5）：601－607.

［12］ 馬瑜，王利生，唐淵圓.三維圖像中階梯型邊緣曲面的追蹤算法［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)報(bào)，2007，19（3）：329－333.Ma Y，Wang L S，Tang Y Y.A novel algorithm for tracking step－like edge surfaces within 3Dimages［J］.Journal of Computer－Aided Design ＆ ComputerGraphics，2007，19（3）：329－333.（in Chinese）

［13］ 黃果，蒲亦非，陳慶利，等.基于分?jǐn)?shù)階積分的圖像去噪［J］.系統(tǒng)工程與電子技術(shù)，2011，33（4）：925－932.Huang G，Pu Y F，Cheng Q L，et al.Research on image denoising based on fractional order integral［J］.Systems Engineering and Electronics，2011，33（4）：925－932.（in Chinese）

［14］ Tao R，Deng B，Wang Y.Research progress of the fractional Fourier transform in signal processing［J］.Science in China Series F，2006，49（1）：1－25.

［15］ Oldham K B，Spanier J.The Fractional Calculus：Theory and Applications of Differentiation and Integration to Arbitrary Order ［M］.New York：Academic Press，1974.

［16］ Podlubny I.Fractional Differential equations：An introduction to Fractional Derivatives，F(xiàn)ractional Differential Equations，To methods of Their Solution and Some of Their Applications ［M］.Academic press，1998.