劉紅亮,李 璟

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

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

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

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

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

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

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

(1)

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

(2)

(3)

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

(4)

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

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

2 遺傳算法

2.1 遺傳算法簡(jiǎn)介

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 參數(shù)選擇及評(píng)價(jià)指標(biāo)

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

通常來(lái)說(shuō)，種群數(shù)目設(shè)置較大時(shí)可以同時(shí)處理更多的解，也更容易找到全局最優(yōu)解；但是，這也會(huì)增加每次迭代的時(shí)間，于是實(shí)際應(yīng)用中一般選取范圍在20～100，本文分別選取50和100進(jìn)行實(shí)驗(yàn).交叉概率和變異概率本文采用自適應(yīng)的方法，根據(jù)個(gè)體適應(yīng)度值的大小來(lái)自動(dòng)調(diào)整.遺傳算法中染色體長(zhǎng)度一般根據(jù)問(wèn)題求解精度的要求來(lái)決定，精度越高，染色體的長(zhǎng)度就要設(shè)置越長(zhǎng)，本文中有三個(gè)基因段，染色體長(zhǎng)度分別選取24，30，36來(lái)進(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]，而對(duì)于迭代次數(shù)t和擴(kuò)散門(mén)限k并沒(méi)有統(tǒng)一的標(biāo)準(zhǔn)或范圍，在后人的研究中，一般這兩個(gè)參數(shù)都是根據(jù)經(jīng)驗(yàn)或者試湊法來(lái)確定具體取值大小.在本試驗(yàn)中，將迭代次數(shù)t取值范圍限制在1～300，將擴(kuò)散門(mén)限k取值范圍限制在1～100.

3.3 相關(guān)評(píng)價(jià)指標(biāo)

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

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

(9)

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

(10)

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

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

(11)

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

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

由于無(wú)法獲取無(wú)噪聲的真實(shí)MRI圖像，為便于分析算法性能，本文選取BrainWeb database網(wǎng)站上的計(jì)算機(jī)模擬MRI圖像(T1-weighted)進(jìn)行實(shí)驗(yàn)，并且添加噪聲強(qiáng)度為3%、5%、7%和9%的萊斯噪聲，通過(guò)濾波后的噪聲圖像與無(wú)噪聲圖像對(duì)比，進(jìn)而評(píng)估濾波器性能.在本文中，所有的實(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所示，對(duì)于不同程度噪聲條件下MRI在遺傳算法優(yōu)化過(guò)程中圖像濾波后的最佳PSNR變化情況，從圖中可以看出，從300代以后變化已經(jīng)很小甚至幾乎沒(méi)有變化，所以經(jīng)過(guò)多次實(shí)驗(yàn)，在遺傳算法優(yōu)化各向異性擴(kuò)散濾波過(guò)程中，遺傳操作進(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%噪聲情況下去噪效果對(duì)比Figure 3 Comparison of denoising result for3% noised image

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

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

(a)9%噪聲圖 (b)優(yōu)化后參數(shù)去噪效果(c)其他參數(shù)去噪效果圖6 9%噪聲情況下去噪效果對(duì)比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%噪聲的去噪效果對(duì)比。

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

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

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

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

5 小 結(jié)

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

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