基于多核的MCMC圖像去噪算法并行實現(xiàn)

成敏，張建州，于巖

1.四川大學(xué)計算機學(xué)院，成都 610065

2.中國人民解放軍78098部隊

基于多核的MCMC圖像去噪算法并行實現(xiàn)

成敏1，2，張建州1，于巖1

1.四川大學(xué)計算機學(xué)院，成都 610065

2.中國人民解放軍78098部隊

圖像在獲取、傳輸和存貯的過程中，總是不可避免地受到各種噪聲源的干擾。為抑制圖像噪聲，改善圖像質(zhì)量，從圖像中獲取更準確的信息，需要對圖像進行去噪處理。圖像去噪是圖像分割、復(fù)原等后續(xù)處理的重要前提條件[1]，是圖像處理領(lǐng)域中十分關(guān)鍵的技術(shù)，也一直是圖像處理領(lǐng)域的難題。

針對高斯噪聲，常見的去噪方法有：高斯濾波、小波變換域濾波、非局部均值濾波[2-3]、偏微分方程濾波等。隨著馬爾可夫鏈蒙特卡洛方法在統(tǒng)計學(xué)應(yīng)用領(lǐng)域的發(fā)展，馬爾可夫鏈蒙特卡洛方法也運用到了圖像去噪中。

馬爾可夫鏈蒙特卡洛去噪算法[4]是通過隨機抽取變化的圖像，根據(jù)樣本的條件概率分布，制定一個貝葉斯最小二乘法[5-6]，在S.Geman和D.Geman[7]對圖像的貝葉斯重建研究中，最大限度地得到原始圖像條件概率，然后運用后驗估計的迭代模擬方法解決圖像去噪，這是一種非參數(shù)的隨機方法。馬爾可夫鏈蒙特卡洛的去噪方法，以靈活的方式適應(yīng)了原始圖像的噪聲統(tǒng)計。

高斯噪聲圖像運用MCMC方法去噪后，圖像有明顯斑點，高噪聲情況下去噪效果不理想。Chatterjee和Milanfar[8-9]在高斯噪聲環(huán)境下的實驗說明，對圖像進行預(yù)去噪處理可以降低噪聲產(chǎn)生的影響，有利于進行分類和聚類。因此，本文在用MCMC方法進行去噪前先對噪聲圖像進行簡單的均值濾波預(yù)處理，既去除了圖像中部分噪聲，也有利于隨機聚類，減少斑點。

在MCMC方法中，需要知道噪聲方差，但在實際應(yīng)用中，噪聲方差顯然是不可知的，所以需要對噪聲方差進行估計。本文運用成熟的Donoho和Johnstone[10]頻域采樣法進行噪聲方差估計。

由于在MCMC方法中，每個像素點的計算都是獨立的，本文利用多核硬件計算資源，將算法進行并行處理[11]，提高了計算速度。在圖像規(guī)模較大的情況下，加速比更加明顯。

1 馬爾可夫鏈蒙特卡洛方法

1.1 問題模型

設(shè)S是圖像，s是圖像上的一個像素，s∈S。被噪聲污染的原始圖像是F={F(s)|s∈S}，無噪聲的圖像是G={G(s)|s∈S}，s上的隨機噪聲為N={N(s)|s∈S}。F，G和N之間的關(guān)系通?？梢员硎緸閇12]：

求的是G(s)，從公式（1）可看出，這其實是一個逆問題，最終只能求得G(s)的估計值，因此這被視為解決后驗估計的問題。運用貝葉斯最小二乘法來解決這個估計問題。

從公式（3）可以知道，G(s)的最佳估計值是E(G(s)|F(s))，G(s)是F(s)條件下的均值。但需要用一個非常復(fù)雜的非線性函數(shù)F(s)來計算。若要解決這個問題，計算式可作一個簡單的貝葉斯線性最小二乘式，建模后p(G(s)|F(s))使用參數(shù)統(tǒng)計模型。但使用這種方法往往導(dǎo)致噪聲強度大的圖像去噪更加復(fù)雜。

為使G(s)的估計值更接近于真實值，分析后，決定使用重要性加權(quán)馬爾可夫鏈蒙特卡洛后驗估計方法。這種方法在一個全圖的范圍內(nèi)考慮該像素的灰度值，將非參數(shù)近似計算p(G(s)|F(s))記為(G(s)|F(s))。在這種方法中本地空間強度的相互作用重在允許改善無噪聲圖像G(s)，噪聲N(s)的地方空間強度的相互作用在這里就會表現(xiàn)得非常弱。

1.2 馬爾可夫鏈蒙特卡洛去噪算法

初始化：Ω是一個估計出來的目標(biāo)分配序列樣本集，Ω={s0，s1，…，sη}。初始序列樣本Ω，令s0=s。

第1步第k次迭代的候選樣本從公式（4）的Q(s′k|sk-1)分布中隨機抽取。

其中σs表示分布函數(shù)Q(s′|sk-1)的方差。

第2步由公式（5）計算接受概率α(s′k|sk-1)。

第3步由均勻分布U(0，1)生成一個隨機值u。

第4步如果u≤α(s′k|sk-1)，將候選樣本s′k納入樣本序列Ω中，否則，舍棄s′k。

第5步重復(fù)1～4步，迭代m次。

第6步將樣本序列Ω中的樣本權(quán)值，建立一個權(quán)值集[13]。

第7步將公式（7）中的權(quán)值代入公式（8）中，估計后驗概率密度函數(shù)(G(s)|F(s))。

1.3 噪聲圖像均值濾波預(yù)處理

均值濾波是一種平滑線性空間濾波，將包含在濾波掩模確定鄰域內(nèi)像素作簡單平均，用平均灰度值去代替圖像每個像素點的值。圖像灰度尖銳變化是圖像邊緣的特性，本文采用均值濾波作預(yù)處理后，圖像邊緣特性削弱，所以處理后存在著不希望的邊緣模糊的負面效應(yīng)。

1.4 噪聲方差估計

公式（6）需要知道噪聲方差，但在實際應(yīng)用中，處理含有噪聲的圖像時，噪聲方差是未知的，但噪聲方差又是算法中的重要參數(shù)，那就只有從噪聲圖像中估計得出。

噪聲方差的估計方法有很多，如K.Rank基于空域的噪聲標(biāo)準差估計，M.Jansen基于小波域的噪聲方差估計方法，魯棒中值絕對估計方法等。本文選用Donoho和Johnstone的頻域采樣法，這種方法是最簡便的。

Donoho噪聲估計方法是建立在理想的高斯噪聲模型基礎(chǔ)上，將圖像小波分解的第一級斜向高頻子圖像作為純子噪聲圖來估計原噪聲標(biāo)準方差。公式是MAD/0.674 5，其中MAD是HH帶小波系數(shù)幅度的絕對值的中值，理想高斯噪聲的絕對值中值和標(biāo)準差的比約是0.674 5。這種采樣方法適合高頻段含很少圖像信息的場合，得到的子圖像點數(shù)較多，保證了子噪聲對原噪聲特征的繼承。

2 馬爾可夫鏈蒙特卡洛算法的并行實現(xiàn)

2.1 并行計算

圖像處理通常存在著時間復(fù)雜度和結(jié)果圖像質(zhì)量相互矛盾的問題，為了提高圖像的處理效果，通常需要花費更多的時間。當(dāng)串行執(zhí)行的性能接近最大時，只有通過并行計算等方式提高算法的運算速度。在單機環(huán)境下，用多核計算機實現(xiàn)并行計算。在分布式環(huán)境下，通常由多處理器和多核計算機來構(gòu)成分布式的環(huán)境。

馬爾可夫鏈蒙特卡洛方法需要計算多次循環(huán)迭代，并且循環(huán)不依賴于其他迭代的結(jié)果，特別適用于并行計算，所以本文將循環(huán)迭代進行分組并行計算。因為循環(huán)內(nèi)有通信消耗，所以只有進行大量的復(fù)雜計算時才會大大降低時間復(fù)雜度。

2.2 并行算法實現(xiàn)

分析本文算法耗時原因發(fā)現(xiàn)，計算每個像素點的去噪過程是相同的，將這個圖像的大矩陣計算過程進行任務(wù)分配，可以在Matlab軟件環(huán)境下實現(xiàn)并行性[14]。

馬爾可夫鏈蒙特卡洛算法的并行過程如圖1所示。

由圖1可以看出，循環(huán)體的每次執(zhí)行都是一次迭代，Matlab可以以非特定順序相互獨立地計算每個迭代。因為每個像素點的計算都是獨立的，所以沒必要保證這些計算的同步性。如果CPU的數(shù)量與迭代次數(shù)相等，那么每個CPU執(zhí)行一次迭代。如果CPU數(shù)量少于迭代次數(shù)，那么有些CPU將執(zhí)行多次迭代；在這種情況下，為了減少通信的時間，一個CPU可能一次接到多次迭代。

圖1 并行計算流程

3 實驗

3.1 實驗環(huán)境及相關(guān)參數(shù)

實驗環(huán)境為：Intel?CoreTMi5 CPU M430@2.27 GHz內(nèi)存2 GB，Win7旗艦版，Matlab R2010（b）。參數(shù)的選取為，均值預(yù)處理窗口為3×3，空間方差σs=21，迭代次數(shù)m=200，搜索窗口為7×7。本文只對該算法的并行計算有效性和可行性進行測試，所以，只采用簡單的單機多核進行實驗。在以后的實驗中將會由多處理器和多核計算機來構(gòu)成分布式環(huán)境。

在Matlab中，首先要用matlabpool[15]啟動并行運行環(huán)境，設(shè)置計算機CPU數(shù)量，為了最大地達到計算機的運行速度，一般按照計算機CPU的核數(shù)進行設(shè)置，在實驗中設(shè)置為雙核。雙核運算功能啟動后，Windows任務(wù)管理器里有三個Matlab.exe進程，其中一個占有內(nèi)存較多，另兩個占有內(nèi)存一樣多。在設(shè)計MCMC圖像去噪并行處理時，需要使用parfor語句，要注意將不變參數(shù)設(shè)置在循環(huán)體外，生成的隨機數(shù)等要設(shè)置在循環(huán)體內(nèi)，否則運算出來的結(jié)果會發(fā)生錯誤。

3.2 本文算法與傳統(tǒng)去噪算法效果對比

首先，驗證本文去噪方法的有效性，分別選取大小為512×512的“Lena”圖像、512×512的“Peppers”圖像添加噪聲進行實驗，并對高斯平滑濾波、維納濾波[16]進行去噪效果的比較，實驗具有很強代表性，實驗結(jié)果如圖2，圖3所示。

圖2 σn=50圖像去噪

圖3 σn=50圖像去噪

從圖2、圖3可以看出，常用的一些圖像去噪算法在一定程度上削弱了圖像噪聲，提高了圖像的質(zhì)量，但用高斯平滑去噪方法去噪后，圖像變得模糊，邊緣特征不清晰；用維納去噪后圖像上仍存在許多明顯噪點。本文算法去噪后的圖像，噪點基本被去除，邊緣特征較清晰。

從圖中看出，用MCMC方法可以獲得比較好的去噪效果，下面定量分析圖像去噪效果。用幾種去噪算法對Lena圖像和Peppers圖像的不同噪聲方差圖像去噪后所得圖像峰值信噪比如表1所示。從表1看出，當(dāng)噪聲水平很小時，維納去噪方法的效果最好，本文算法略低于該方法，但隨著噪聲方差的增加，維納去噪方法的去噪效果越來越不理想，本文提出的算法去噪效果是最好的。

3.3 算法運行時間比較分析

從表2看出，矩陣規(guī)模較大（512×512）比計算規(guī)模較小（256×256）多核加速明顯，這是因為將任務(wù)分配給CPU時耗時造成，在計算規(guī)模很小的情況下并行得來的性能提升有限，不能彌補這種消耗。

表1 PSNR比較表dB

表2 單核和雙核運行效率對比s

4 結(jié)論與展望

本文用MCMC的并行方法對圖像去噪，既保證了去噪效果也提高了運行速度。在這項研究中，通過對噪聲圖像的預(yù)處理，估計噪聲方差，使用馬爾可夫鏈蒙特卡洛抽樣方法進行隨機優(yōu)化，達到圖像去噪的目的。實驗是采用添加噪聲然后進行去噪來測試該算法信噪比，用數(shù)據(jù)說明了其去噪性能的優(yōu)越性。該算法還可擴展為多通道圖像，視頻和3D數(shù)據(jù)量。

為了提高MCMC算法的運行速度，特別是該算法運用到視頻和3D上時為了解決其計算量的龐大，本文針對MCMC算法的可并行性進行了研究與實現(xiàn)，提高了大規(guī)模矩陣的運算速度和多核CPU資源的利用率。本文的并行設(shè)計方法為其他類似的數(shù)據(jù)處理提供了一種將串行任務(wù)進行分割的新途徑。Matlab豐富的庫函數(shù)和操作系統(tǒng)上的集群環(huán)境、多核的硬件條件可更簡捷有效地解決串行程序性能提升方面的瓶頸問題。目前，Matlab并行計算應(yīng)用發(fā)展迅速，將圖像放在分布式的環(huán)境下將大大提升運行速度。

[1]Gonzalez C.Digital image processing[M].2nd ed.Beijing：Publishing House of Electronics Industry，2011.

[2]Buades A，Coll B，Morel J M.A review of image denoising algorithms，with a new one[J].Multiscale Model Simul，2005，4（2）：490-530.

[3]Buades A，Coll B，Morel J M.On image denoising methods[J].SIAM Review，2010，52（1）：113-147.

[4]Wong A，Mishra A.Stochastic image denoising based on Markov-chain Monte Carlo sampling[J].Signal Processing 2011，91：2112-2120.

[5]Koutrounbss K.Pattern recognition[M].Beijing：Publishing House of Electronics Industry，2009.

[6]Pryce J，Bruce A.Statistical mechanics of image restoration[J].Journal of Physics A：Mathematical and General，1995，28（3）：511-532.

[7]Geman S，Geman D.Stochastic relaxation，Gibbs distributions，and the Bayesian restoration of images[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，1984，6：721-741.

[8]Chatterjee P，Milanfar P.Is denoising dead?[J].IEEE Transactions on Image Pocessing，2010，19（4）：895-911.

[9]Chatterjee P，Milanfar P.Practical bounds on image denoising：from estimation to information[J].IEEE Transactions on Image Pocessing，2011，20（5）：1221-1233.

[10]Donoho D L，Johnstone I M.Ideal spatial adaptation by wavelet shrinkage[J].Biometrika，1994，81：425-455.

[11]Downton A，Crookes D.Parallel architectures for image processing[J].Electronics&Communication Engineering Journal，1998，10（3）.

[12]Chan TonyF，Shen Jianhon.Image processing and analysis[M].Beijing：Publishing House of Science，2009.

[13]Chen M.Importance-weighted marginal Bayesian posterior density estimation[J].Journal of the American Statistical Association，1994，89（427）：818-824.

[14]MathWorks.Matlab distributed computing server[EB/OL].（2009-10-01）.http：//www.mathworks.com/help/toolbox/ mdce/index.html.

[15]Gonzalez C.Digital image processing using MATLAB[M]. Beijing：Publishing House of Electronics Industry，2009.

[16]Nowak R D.Wavelet-based rician noise removal for magnetic resonance imaging[J].IEEE Transactions on Image Pocessing，1999，8（10）：1408-1419.

CHENG Min1，2,ZHANG Jianzhou1,YU Yan1

1.College of Computer,Sichuan University,Chengdu 610065,China
2.Unit 78098 of PLA,China

Image denoising is a prerequisite for the many processing tasks of image.Markov Chain Monte Carlo algorithm is an important method of image denoising.However,the method has some problems such that the denoised image has obvious spots,the denoising image corrupted by heavy noise is not satisfactory,the noise variance needs to be estimated in practical application,and the operation speed of this method is slow.This paper proposes a two-step denoising method. It preprocesses the noise image using the mean filter.It estimates the pretreated image noise variance.It uses the MCMC image denoising method.To take full advantage of multi-core processor resources,this paper studies the parallel programming of MCMC algorithm.The multi-core program increases the speed of MCMC algorithm.The experiments show that the denoising method given in this paper reduces spots and improves the signal-to-noise ratio.Parallel processing can make the algorithm more efficient.

image denoising;Markov Chain Monte Carlo（MCMC）;variance estimation;pretreatment;parallel processing

圖像去噪是許多圖像處理任務(wù)的前提。馬爾可夫鏈蒙特卡洛圖像去噪算法是很重要的一種圖像去噪方法，但去噪后圖像存在明顯斑點，在高噪聲情況下去噪效果不理想，實際應(yīng)用中需要進行噪聲方差估計，運算速度慢。提出兩步去噪方法，用均值濾波對噪聲圖像進行預(yù)處理，估計預(yù)處理后圖像噪聲方差，進行MCMC圖像去噪；為充分利用多核處理器的硬件資源，研究了將MCMC算法進行并行編程，提高了程序的運行速度。實驗表明兩步去噪方法減少了斑點、提高了信噪比；并行實現(xiàn)提高了運算效率。

圖像去噪；馬爾可夫鏈蒙特卡洛方法（MCMC）；方差估計；預(yù)處理；并行處理

A

TP391.9

10.3778/j.issn.1002-8331.1211-0155

CHENG Min,ZHANG Jianzhou,YU Yan.Parallel image denoising algorithm based on multi-core MCMC.Computer Engineering and Applications,2014,50（18）：152-155.

成敏（1976—），女，工程師，主要研究領(lǐng)域為計算機視覺。E-mail：scsdcm@sohu.com

2012-11-14

2013-01-10

1002-8331（2014）18-0152-04

CNKI網(wǎng)絡(luò)優(yōu)先出版：2013-01-29，http://www.cnki.net/kcms/detail/11.2127.TP.20130129.1539.011.html