劉紅亮,李 璟

(中國計(jì)量大學(xué) 機(jī)電工程學(xué)院，浙江 杭州 310018)

在醫(yī)學(xué)圖像處理過程中，去噪是很重要的一個預(yù)處理環(huán)節(jié).MRI(magnetic resonance imaging)可以提供詳細(xì)的身體內(nèi)部結(jié)構(gòu)和器官的斷層圖像，特別是可以很好地區(qū)分身體中不同軟組織結(jié)構(gòu)，所以MRI圖像越來越多的成為醫(yī)生診斷疾病的依據(jù)，以及科學(xué)工作者的研究材料.然而，在獲取圖像的過程中，總會受到外界噪聲的影響，而對于MRI來說其大部分噪聲來源于被測試的人體本身.除此之外，MRI成像系統(tǒng)自身的線圈和電子器件等也會成為噪聲源[1].早期的去噪算法一般是將MRI中的噪聲估計(jì)為高斯噪聲來計(jì)算的，實(shí)際上，MRI中的噪聲在后來的研究中被認(rèn)為是萊斯分布的噪聲.MRI中的這些噪聲極大地干擾了臨床醫(yī)生對于疾病的診斷，同時也限制了對其進(jìn)行的分析研究.因此，為了提高圖像質(zhì)量以便于正確的疾病診斷和分析，性能優(yōu)異的去噪算法及其實(shí)踐中的參數(shù)調(diào)校變得越來越重要.

在去噪過程中，相對于其他去噪算法來說，各向異性擴(kuò)散由于具有較好的保持邊緣和局部細(xì)節(jié)的能力，在目前的相關(guān)研究中已經(jīng)成為醫(yī)學(xué)圖像處理中廣泛應(yīng)用的去噪技術(shù).各向異性擴(kuò)散濾波是由Perona和Malik[2]于1990年提出的，其基本原理是：用偏微分方程作為擴(kuò)散方程，通過對圖像中不同方向上的梯度值進(jìn)行計(jì)算，進(jìn)而確定不同方向上的擴(kuò)散系數(shù)，最終使得在平滑噪聲的過程中同時也提升了保留局部細(xì)節(jié)的能力.

盡管各向異性擴(kuò)散濾波被廣泛應(yīng)用于醫(yī)學(xué)圖像去噪的預(yù)處理中，目前為止，對于在去噪過程中各向異性擴(kuò)散濾波的各個參數(shù)選擇優(yōu)化問題，還沒有一個系統(tǒng)完整的研究或驗(yàn)證.通過查閱相關(guān)文獻(xiàn)可以看出，在近年來的研究工作中，大部分研究者只著眼于算法本身的改進(jìn)，Hongjin Ma等人[3]提出了通過構(gòu)建合適的擴(kuò)散張量構(gòu)建新的各向異性擴(kuò)散模型使其在去噪過程中可以更好地保留原圖像信息；Justin Joseph[4]建立了通過圖像梯度均值來計(jì)算最佳梯度系數(shù)閾值的分析模型以改善各向異性擴(kuò)散濾波在2D圖像中的性能.Jianjun Yuan[5]提出了一個基于非局部均值理論的新模型，Mariem Ben Abdallah等人[6]引入了自動估計(jì)RGB噪聲模型的函數(shù)，Sanghun Kim 等人[7]提出了基于區(qū)域的自適應(yīng)平滑方法來提升ADF的性能，Jiangtao Xu 等人通過局部差值辨別出是否為噪聲點(diǎn)并用半自適應(yīng)閾值的方法來改善原來算法.而參數(shù)優(yōu)化方面，Ricardo J Ferrari 通過自動確定各向異性擴(kuò)散濾波的迭代次數(shù)來改善去噪性能，另外兩個參數(shù)并沒有優(yōu)化.然而，每個參數(shù)的改變都會對去噪效果產(chǎn)生影響，只有對所有參數(shù)進(jìn)行優(yōu)化，才能使得各向異性擴(kuò)散濾波在去噪過程中達(dá)到最優(yōu)性能.

從誕生以來，遺傳算法由于其強(qiáng)大的優(yōu)化性能，被越來越廣泛的運(yùn)用到各種算法優(yōu)化過程中，遺傳算法包括了計(jì)算簡單、更少的迭代次數(shù)和更少的數(shù)學(xué)復(fù)雜度等優(yōu)點(diǎn)[9].遺傳算法通?？梢钥焖僬业絾栴}的最優(yōu)解決方案，相對于傳統(tǒng)的優(yōu)化算法，其需要更少的先驗(yàn)知識.本文通過用改進(jìn)的遺傳算法(genetic algorithm, GA)來優(yōu)化ADF的各項(xiàng)參數(shù)，并且在不同噪聲條件下對各向異性擴(kuò)散濾波的參數(shù)進(jìn)行優(yōu)化選擇，而且還與未優(yōu)化參數(shù)的ADF的去噪性能進(jìn)行對比.

1 各向異性擴(kuò)散模型(PM模型)

在PM模型中，一般經(jīng)常采用梯度微分算子來識別圖像邊緣，從而可以將邊緣檢測和去除噪聲很好地統(tǒng)一起來，即統(tǒng)一到了基于變分法的偏微分方程中.Perona和Malik[2]提出的各向異性擴(kuò)散方程如下.

(1)

式(1)中，I表示需要處理的圖像，表示檢測邊緣的梯度算子，div用來計(jì)算散度，||·|| 表示圖像灰度值的幅度，表示擴(kuò)散方程，t作為方程中引入的時間算子.PM模型的基本原理是，通過計(jì)算的大小進(jìn)而有選擇地?cái)U(kuò)散平滑所要處理的圖像.Perona等人基于梯度與擴(kuò)散閾值的關(guān)系提出了具有兩種不同形式的擴(kuò)散方程.

(2)

(3)

離散化后的表達(dá)式為：

(4)

式(4)中，0≤λ≤1/4，N、S、E、W 分別代表北、南、東、西四個方向，并且有：

可以看出，在算法執(zhí)行過程中，需要設(shè)置三個參數(shù)，迭代次數(shù)t，參數(shù)λ，以及擴(kuò)散閾值k，每個參數(shù)的選擇都直接影響到濾波器的性能.比如，擴(kuò)散閾值參數(shù)k大的選擇就比較難以控制，設(shè)置不合適往往會出現(xiàn)明顯的“階梯”效應(yīng).單純的憑借經(jīng)驗(yàn)，或者用試湊法，并不能很好的發(fā)揮出濾波器的性能.

2 遺傳算法

2.1 遺傳算法簡介

1975年John Holland[10]提出了基于生物界自然選擇和自然遺傳的遺傳算法(genetic algorithm, GA)，這是一種在優(yōu)化過程中的啟發(fā)式隨機(jī)搜索算法.一方面，兩個個體之間通過基因交換來產(chǎn)生新的個體；另一方面，變異是基因隨機(jī)變化的結(jié)果[11].因此，變異的目的是為了在基因變化的過程中盡量提高數(shù)據(jù)的多樣性，從而可以避免在算法優(yōu)化過程中導(dǎo)致遺傳算法陷入局部最優(yōu)解.在整個過程中，產(chǎn)生的新個體有優(yōu)于舊個體的，同時也有比舊個體更差的，這樣在選擇的過程中就要剔除較差的個體，同時把更優(yōu)的個體保留下來[12]，如此就完成了選擇優(yōu)化的過程.

從本質(zhì)上來說，遺傳算法對于其他算法的優(yōu)化過程就是一個迭代過程，一般執(zhí)行過程如下：

1)確定所要優(yōu)化的參數(shù)，并將其轉(zhuǎn)換到二進(jìn)制串對應(yīng)的染色體中；

2)根據(jù)優(yōu)化目標(biāo)定義合適的適應(yīng)度函數(shù)；

3)選取合適的選擇、交叉和變異方法，并且確定交叉和變異概率；

4)設(shè)置初始種群大小，隨機(jī)初始種群，并且計(jì)算所有個體的適應(yīng)度值；

5)按照選擇、交叉、變異的順序?qū)ΨN群進(jìn)行操作得到下一代群體；

6)對新的種群進(jìn)行評估，如果達(dá)到預(yù)定迭代次數(shù)或者滿足某一指標(biāo)，則停止計(jì)算，否則返回第5)步.

在本文中，遺傳算法的實(shí)現(xiàn)及ADF的參數(shù)優(yōu)化過程，均在MATLAB-2016a的環(huán)境下進(jìn)行.

2.2 本文所使用的改進(jìn)的遺傳算法

傳統(tǒng)的遺傳算法中，種群中個體的選擇概率是按比例的適應(yīng)度來分配的，即個體的適應(yīng)度值和種群中所有個體適應(yīng)度值的和的比值，比值(概率)越大，被選中的概率就越高.對于本文的各向異性參數(shù)優(yōu)化來說，在試驗(yàn)過程中出現(xiàn)了明顯的震蕩現(xiàn)象，即優(yōu)劣個體并沒有完全按照預(yù)期的那樣被選擇出來，優(yōu)勢個體在下一代中可能被劣勢個體代替，最優(yōu)個體也不能保證在下一代中完全保存下來.

基于以上原因，本文采用了一種可以保證最優(yōu)個體保存下來的一種新的精英選擇法.具體操作是：首先，在當(dāng)代中找出適應(yīng)度值最大的個體，并且用它來取代最劣勢個體；然后，將原來最優(yōu)個體重新用一個隨機(jī)產(chǎn)生的個體取代，并且與原來次優(yōu)個體比較，如果優(yōu)于次優(yōu)個體則取代次優(yōu)個體，并且再將次優(yōu)個體用一個隨機(jī)產(chǎn)生的個體取代，若劣于次優(yōu)個體則保持不變.

Srinvivas等人提出了一種自適應(yīng)的遺傳算法，如果種群中的個體適應(yīng)度值趨于一致或者收斂域局部最優(yōu)，就會使交叉概率Pc和變異概率Pm增大；然而當(dāng)適應(yīng)度比較分散時，則會減小Pc和Pm的值.具體計(jì)算表達(dá)式如下.

其中，Pc1=0.9 ，Pc2=0.6 ，Pm1=0.1 ，Pm2=0.001.

最后，本文提出了將精英選擇策略與自適應(yīng)遺傳算法結(jié)合起來，并且作為一種改進(jìn)的遺傳算法，對各向異性擴(kuò)散濾波參數(shù)進(jìn)行優(yōu)化.這樣不僅保證了種群中每代最優(yōu)個體的保存，并且可以避免陷入局部最優(yōu)，在遺傳過程中保證了種群個體的多樣性，同時，也加快了收斂速度.

3 參數(shù)選擇及評價指標(biāo)

3.1 遺傳算法的參數(shù)選擇

通常來說，種群數(shù)目設(shè)置較大時可以同時處理更多的解，也更容易找到全局最優(yōu)解；但是，這也會增加每次迭代的時間，于是實(shí)際應(yīng)用中一般選取范圍在20～100，本文分別選取50和100進(jìn)行實(shí)驗(yàn).交叉概率和變異概率本文采用自適應(yīng)的方法，根據(jù)個體適應(yīng)度值的大小來自動調(diào)整.遺傳算法中染色體長度一般根據(jù)問題求解精度的要求來決定，精度越高，染色體的長度就要設(shè)置越長，本文中有三個基因段，染色體長度分別選取24，30，36來進(jìn)行實(shí)驗(yàn).最大遺傳代數(shù)一般作為模擬終止的條件，通常視具體情況而定，選取范圍在100～1 000，本文將分別選取50，100，200，300，500進(jìn)行實(shí)驗(yàn).具體如表1.

表1 遺傳算法的參數(shù)設(shè)置

3.2 各向異性擴(kuò)散濾波的相關(guān)參數(shù)選取范圍

Perona和Malik[2]明確給出了λ的取值范圍是[0,0.25]，而對于迭代次數(shù)t和擴(kuò)散門限k并沒有統(tǒng)一的標(biāo)準(zhǔn)或范圍，在后人的研究中，一般這兩個參數(shù)都是根據(jù)經(jīng)驗(yàn)或者試湊法來確定具體取值大小.在本試驗(yàn)中，將迭代次數(shù)t取值范圍限制在1～300，將擴(kuò)散門限k取值范圍限制在1～100.

3.3 相關(guān)評價指標(biāo)

為了評估通過遺傳算法選擇出的各向異性擴(kuò)散濾波器參數(shù)，本文選取了如下評價指標(biāo)：峰值信噪比(peak signal-to-noise ratio，PSNR)，結(jié)構(gòu)相似性指數(shù)(structural similarity index metric，SSIM)，均方差(mean squared error，MSE).

MSE(用來估計(jì)兩幅圖像相異程度)的計(jì)算方法為[13]

(9)

式(9)中，I0是無噪聲的原始圖像，In是添加噪聲以后的圖像，k、l分別代表行列數(shù).

(10)

式(10)中，N表示圖像中的灰度級數(shù).

SSIM是另一種評價兩幅圖像相似度的計(jì)算方法[14]

(11)

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

在遺傳算法的選擇過程中，為選出最佳參數(shù)組合，本實(shí)驗(yàn)先以最佳PSNR為標(biāo)準(zhǔn)選擇出相關(guān)參數(shù)組合，進(jìn)而再將MSE和SSIM的值的大小作為衡量所選濾波器參數(shù)性能的條件.

由于無法獲取無噪聲的真實(shí)MRI圖像，為便于分析算法性能，本文選取BrainWeb database網(wǎng)站上的計(jì)算機(jī)模擬MRI圖像(T1-weighted)進(jìn)行實(shí)驗(yàn)，并且添加噪聲強(qiáng)度為3%、5%、7%和9%的萊斯噪聲，通過濾波后的噪聲圖像與無噪聲圖像對比，進(jìn)而評估濾波器性能.在本文中，所有的實(shí)驗(yàn)均在MATLAB R2016a中進(jìn)行.如圖1所示為參數(shù)優(yōu)化的系統(tǒng)框圖.

圖1 ADF的參數(shù)優(yōu)化系統(tǒng)框圖Figure 1 Optimization system block diagram for ADF parameters

如圖2所示，對于不同程度噪聲條件下MRI在遺傳算法優(yōu)化過程中圖像濾波后的最佳PSNR變化情況，從圖中可以看出，從300代以后變化已經(jīng)很小甚至幾乎沒有變化，所以經(jīng)過多次實(shí)驗(yàn)，在遺傳算法優(yōu)化各向異性擴(kuò)散濾波過程中，遺傳操作進(jìn)行到500 代即可.

圖2 不同噪聲條件下基于最佳PSNR的遺傳算法性能表現(xiàn)Figure 2 Performance of genetic algorithm based on optimal PSNR under different noise

(a)3%噪聲圖 (b)優(yōu)化后參數(shù)去噪效果 (c)其他參數(shù)去噪效果圖3 3%噪聲情況下去噪效果對比Figure 3 Comparison of denoising result for3% noised image

(a)5%噪聲圖 (b)優(yōu)化后參數(shù)去噪效果 (c)其他參數(shù)去噪效果圖4 5%噪聲情況下去噪效果對比Figure 4 Comparison of denoising result for 5% noised image

(a)7%噪聲圖 (b)優(yōu)化后參數(shù)去噪效果(c)其他參數(shù)去噪效果圖5 7%噪聲情況下去噪效果對比Figure 5 Comparison of denoising result for 7% noised image

(a)9%噪聲圖 (b)優(yōu)化后參數(shù)去噪效果(c)其他參數(shù)去噪效果圖6 9%噪聲情況下去噪效果對比Figure 6 Comparison of denoising result for 9% noised image

參數(shù)3%噪聲5%噪聲7%噪聲9%噪聲k53.343831.078191.234417.7969λ0.00100.00200.00300.0010t49434099

為了便于觀察結(jié)果，上面列出3%、5%、%7和9%噪聲的去噪效果對比。

從以上對比圖來看，相對于其他參數(shù)，經(jīng)過參數(shù)優(yōu)化的各向異性濾波的去噪效果更好，在去除噪聲的同時，更好的保持了邊緣，對MRI圖像中的的解剖結(jié)構(gòu)破壞較小，更有利于圖像處理的其他后續(xù)操作.

表3為從數(shù)據(jù)量化方面比較優(yōu)化過參數(shù)的濾波器去噪效果對比.

表3 不同噪聲條件下優(yōu)化后參數(shù)與其他參數(shù)去噪性能對比

從本次實(shí)驗(yàn)結(jié)果可以看出，不管是從視覺上還是從數(shù)據(jù)指標(biāo)上，相對于其他參數(shù)來說，經(jīng)過遺傳算法優(yōu)化過的各向異性擴(kuò)散濾波有更好的去噪效果.去噪是平滑的過程，在去噪過程中勢必會影響到圖像中的結(jié)構(gòu)，特別是一些細(xì)小結(jié)構(gòu).從實(shí)驗(yàn)結(jié)果中的圖像(圖3～6)來看，相對于優(yōu)化過參數(shù)的各向異性濾波來說，其他參數(shù)的濾波結(jié)果要么是看起來更加模糊，要么是噪聲去除得不夠理想.各向異性擴(kuò)散濾波本身的特性決定了它在保證去噪性能的前提下可以較好的保持原有的圖像結(jié)構(gòu)(即邊緣效果)，通過這種參數(shù)優(yōu)化方法可以讓它的去噪效果達(dá)到最優(yōu)，這樣既能最大限度的去除噪聲，同時又能保持邊緣不被破壞.

5 小 結(jié)

為了解決實(shí)際應(yīng)用中各向異性擴(kuò)散濾波參數(shù)的選擇問題，本文采用了一種改進(jìn)的遺傳算法對其參數(shù)進(jìn)行優(yōu)化，提出了一種新的精英選擇策略，并結(jié)合自適應(yīng)的變異和交叉概率，在優(yōu)化過程中保證了最佳個體在每代中的延續(xù)從而避免了陷入局部最優(yōu).同時，將PSNR作為目標(biāo)函數(shù)來選擇最佳個體，進(jìn)而得到最終的不同噪聲條件下的最佳濾波器參數(shù).

本次實(shí)驗(yàn)數(shù)據(jù)是基于電腦模擬的頭部MRI圖像，可以很好地去測試算法性能或者驗(yàn)證算法的有效性，在研究算法方面可以更加方便、快速、有效.但是對于實(shí)際的臨床上的MRI圖像來說，并不存在所謂無噪聲圖像，因?yàn)樵讷@取圖像的過程中勢必會引入噪聲，即臨床上獲取的真實(shí)人體MRI圖像本身就含有噪聲.由于MSE、PSNR都是針對無噪聲圖像和噪聲圖像來說的，所以本文提出的這種方法不能直接用PSNR作為評價指標(biāo).在實(shí)際情況中，對于這種只能獲取帶噪聲的MRI圖像來說，我們可以通過噪聲估計(jì)算法(比如常用的基于統(tǒng)計(jì)評價或者算術(shù)平均的算法)，將噪聲估計(jì)結(jié)果作為一個評價指標(biāo)，并且作為一個遺傳算法過程中的目標(biāo)函數(shù)來對各向異性擴(kuò)散濾波進(jìn)行優(yōu)化.即通過噪聲估計(jì)參數(shù)來代替PSNR，作為遺傳算法的篩選標(biāo)準(zhǔn)，進(jìn)而選擇出最合適的各向異性濾波參數(shù).

[1] NISHIMURA D G.PrinciplesofMagneticResonanceImaging[M]. Stanford:Stanford University, 2010:23-31.

[2] PERONA P, MALIK J. Scale-space and edge detection usinganisotropic diffusion[J].IEEETransactionsonPatternAnalysisandMachineIntelligence, 1990,12(7):629-639.

[3] MA H J, NIE Y F. An edge fusion scheme for image denoising based on anisotropic diffusion models[J].JournalofVisualCommunicationandImageRepresentation, 2016, 40:406-417.

[4] JOSEPH J, PERIYASAMY R. An analytical method for the adaptive computation of threshold of gradient modulus in 2D anisotropic diffusion filter[J].BiocyberneticsandBiomedicalEngineering, 2017, 37(1):1-10.

[5] YUAN J J. Improved anisotropic diffusion equation based on new non-local information scheme for image denoising[J].IETComputerVision, 2015, 9(6): 864-870.

[6] ABDALLAH M B, MALEK J, AZAR A T, et al. Adaptive noise-reducing anisotropic diffusion filter[J].NeuralComputing&Applications, 2016, 27(5): 1273-1300.

[7] KIM S, KANG S J, KIM Y H. Anisotropic diffusion noise filtering using region adaptive smoothing strength[J].JournalofVisualCommunicationandImageRepresentation, 2016,40: 384-391.

[8] XU J T, JIA Y Y, SHI Z F, et al. An improved anisotropic diffusion filter with semi-adaptive threshold for edge preservation[J].SignalProcessing, 2016 ,119(C):80-91.

[9] MISRA D, SARKER S, DHABAL S, et al. Effect of using genetic algorithm to denoise MRI Images corrupted with Rician noise[C] //IEEEInternationalConferenceonEmergingTrendsinComputing,CommunicationandNanotechnology. Tirunelveli, India: IEEE, 2013: 146-151.

[10] HOLLAND J H.AdaptationinNaturalandArtificialSystems:anIntroductoryAnalysiswithApplicationstoBiology,Control,andArtificialIntelligence[M]. Massachusetts, America: MIT Press, 1992:37-43.

[11] PREBYS E K. The genetic algorithm in computer science[J].MITUndergraduateJournalofMathematics, 2007, 3:165-170.

[12] RADWAN A A A, LATEF B A A, MGEID A, et al. Using genetic algorithm to improve information retrieval systems[J].ProceedingsofWorldAcademyofScienceEngineering&Technology, 2008, 17(2):6-12.

[14] WANG Z, BOVIK A C, SHEIKH H R, et al. Image quality assessment: from error visibility to structural similarity[J].IEEETransImageProcess, 2004, 13(4):600-612.

[15] COCOSCO C A, KOLLOKIAN V, KWAN K S, et al. BrainWeb: online interface to a 3D MRI simulated brain database[J].Neuroimage,1997, 5:425.