1、實(shí)驗(yàn)三 圖像恢復(fù)技術(shù)一、 實(shí)驗(yàn)?zāi)康恼莆胀嘶瘓D像中常見噪聲模型及參數(shù)估計(jì)方法；加深對(duì)幾種常用的圖像復(fù)原方法的理解；通過Matlab開發(fā)環(huán)境實(shí)現(xiàn)參數(shù)維納濾波器復(fù)原圖像。二、 實(shí)驗(yàn)內(nèi)容1 對(duì)各種常見噪聲模型進(jìn)行仿真實(shí)現(xiàn)；2 把仿真噪聲加入純凈圖像上形成退化圖像，讓后對(duì)噪聲進(jìn)行參數(shù)估計(jì)，并與原設(shè)定參數(shù)進(jìn)行比較以衡量估計(jì)的準(zhǔn)確性；3 仿真幾種常見噪聲，并用參數(shù)維納濾波器方法進(jìn)行圖像復(fù)原；4 觀察復(fù)原效果，并進(jìn)行仔細(xì)討論。三、 知識(shí)要點(diǎn)與范例1. 圖像恢復(fù)退化模型f(x,y)退化函數(shù)H+g(x,y)n(x,y)恢復(fù)濾波器f(x,y)的最優(yōu)估計(jì)退化過程恢復(fù)過程圖像恢復(fù)問題可以分為兩大類：1) 假定系統(tǒng)退化函

2、數(shù)H為恒等算子，而只考慮環(huán)境噪聲引起退化的情況；2) 同時(shí)考慮H和環(huán)境噪聲h。2. 噪聲模型對(duì)噪聲的行為和效果進(jìn)行仿真的能力是圖像恢復(fù)處理中的重點(diǎn)。除了周期噪聲，通常假定噪聲與象素坐標(biāo)無關(guān)。Matlab噪聲仿真函數(shù)為：g = imnoise(f, type, parameters);常見的噪聲（高斯、椒鹽噪聲）的仿真實(shí)例如下：g = imnoise(f, gaussian, m, var);g = imnoise(f,salt & pepper, d); (where d represents the noise density).3. 產(chǎn)生具有某種固定分布的空間隨機(jī)噪聲.MATLAB

3、中兩個(gè)常用的隨機(jī)分布產(chǎn)生函數(shù)是：rand（均勻分布噪聲）, randn（標(biāo)準(zhǔn)正態(tài)分布噪聲）空間隨機(jī)噪聲是有各自的PDF 或 CDF所決定的。除了有上面所講的immoise噪聲產(chǎn)生函數(shù)外， 具有某種規(guī)定 CDF的大部分隨機(jī)噪聲都可由均勻分布的隨機(jī)量（這里用w表示）通過某個(gè)隨機(jī)數(shù)產(chǎn)生方程來生成。也就是說，給定w (0,1), 具有某個(gè)確定分布Fz 的隨機(jī)變量z 可以通過解如下方程獲得：z = Fz-1(w).如：變量z 滿足如下Rayleigh分布：則：, 這個(gè)方程也通常叫做隨機(jī)數(shù)產(chǎn)生方程 (RNGE)。對(duì)其他類型的 RNGEs, 請(qǐng)見參考教材2的146頁。產(chǎn)生兩種噪聲的Matlab程序?qū)嵗?1)

4、 generate the gaussian noise of 500*500 size, mean=0, variance=1 R = imnoise2(gaussian, 500,500, 0, 1)Whose histogram is below, generated by:hist(R:), 50);(2) rayleigh noise (a and b is by default)r= imnoise2(rayleight, 100000, 1)4. 從灰度直方圖估計(jì)噪聲的方法所要進(jìn)行估計(jì)的參數(shù)主要有兩個(gè)：均值和方差；參考教材2中的Matlab實(shí)現(xiàn)函數(shù)為：v, unv = statm

5、oments(p, n)選取無特征背景區(qū)域 (ROI)的函數(shù)為：B = roipoly(f, c, r), or B, r, c = roipoly() for chose interactively.實(shí)例：仿真某種噪聲加入到某個(gè)無噪聲的圖像，然后把上面方法得到的參數(shù)估計(jì)與原參數(shù)進(jìn)行比較。f=imread('Fig0504(a)(noisy_image).tif');B, c, r = roipoly(f)subplot(221);imshow(f);subplot(222);imshow(B);p,npix = histroi(f,c,r);subplot(223);bar(

6、p,1);v,unv=statmoments(p,2);x = imnoise2(gaussian, npix, 1, floor(unv), floor(sqrt(v);v,unv=statmoments(p,2)x = imnoise2('gaussian', npix, 1, floor(unv(1), floor(sqrt(unv(2);subplot(224), hist(x, 130);axis(0 300 0 50)結(jié)果：5. 只有空間加性噪聲時(shí)的圖像恢復(fù)問題當(dāng)退化僅僅是由上式所示的加性噪聲造成時(shí)，可以采用教材第三章中所講的圖像強(qiáng)技術(shù)。Matlab函數(shù)spfilt

7、（見附錄）集成了差不多所有空間濾波的功能，包括 including: arithmetic mean, geometric mean, harmonic mean, contraharmonic mean, median, max, min, midpoint, alpha-trimmed mean；用上面的函數(shù)進(jìn)行圖像復(fù)原的實(shí)例：原圖像是一幅電路板的x-image，圖像的相關(guān)信息為： unit8, TIFF, 448*464 (heigh*width); 對(duì)圖像加入椒鹽噪聲后再用各種方法進(jìn)行恢復(fù)：% adding pepper-salt noiseM, N=size(f);R=imnoise

8、2('salt & pepper',M, N, 0.1,0);c = find(R=0);gp = f;gp(c)=0;% the following codes generate only salt noise only: (see upper-middle figure)R=imnoise2('salt & pepper',M, N, 0,0.1);c = find(R=1);gs = f;gs(c)=255;%the upper-right figure is generated by filtering pepper noise usin

9、g contraharmonic filter with Q>0:fp = spfilt(gp, 'chmean',3,3,1.5);%the bottom-left figure is generated by filtering salt noise using contraharmonic filter with Q<0:fs = spfilt(gs, 'chmean',3,3,-1.5);%the bottom-middle figure is generated by filtering pepper noise using max fil

10、ter:fpmax = spfilt(gp,'max',3,3);%the bottom-right figure is generated by filtering salt noise using min filter:fsmin = spfilt(gs,'min',3,3);=6. 空間退化函數(shù)建模有先驗(yàn)知識(shí)的退化函數(shù)建模函數(shù)：PSF = fspecial(motion, len , theta); %it approximates the effects of linear motion of a camera.為比較各種恢復(fù)方法而設(shè)計(jì)的一種合適的測(cè)試模

11、式：C = checkerboard(NP, M, N);一個(gè)產(chǎn)生“checkerboard” 圖像和其測(cè)試圖像的實(shí)例：+f = checkerboard(8);PSF = fspecial ('motion', 7, 45); % generate degradation functiongb = imfilter(f, PSF, 'circular'); % producing degraded imagenoise = imnoise(zeros(size(f), 'gaussian', 0, 0.001); % simulate Gaus

12、sian noiseg = gb + noise; % generate blurred image by degradation and noisefr = pixeldup(f, 8); % pixel replication to form 512*512 imagegbr = pixeldup(gb, 8);nr = pixeldup(noise,8);gr = pixeldup(g, 8);subplot(221);imshow(fr,);title(' original test pattern of size 512*512');subplot(222);imsh

13、ow(gbr, );title(' degreded image by PSF');subplot(223);imshow(nr, );title(' gaussian noise image');subplot(224);imshow(gr, );title('blurred image by degradation and noise');+7. Wiener filtering基本原理：Matlab實(shí)現(xiàn)：fr = deconvwnr(g, PSF, NACORR, FACORR);四、 參數(shù)維納濾波參考程序兩種維納濾波方法與直接逆濾波方法間

14、的比較： +f = checkerboard(8); % the original undegraded imagePSF = fspecial ('motion', 7, 45); % degradation functiongb = imfilter(f, PSF, 'circular'); % producing degraded imagenoise = imnoise(zeros(size(f), 'gaussian', 0, 0.001); % simulate Gaussian noiseg = gb + noise; % degr

15、aded image mixed with noise % direct verse filteringfr1 = deconvwnr(g, PSF); % parametric Wiener filtering (assume noise-to-signal ratio is average power ratio)Sf = abs(fft2(f) . 2; % clear image power spectrumfA = sum(Sf(:)/prod(size(f); % the image average power; Sn = abs(fft2(noise) . 2; % noise

16、power spectrumnA = sum(Sn(:)/prod(size(noise); % the average noise power;R = nA/fA;fr2 = deconvwnr(g, PSF, R) % parametric Wiener filtering by autocorrelation functionsNCORR = fftshift(real(ifft2(Sn);ICORR = fftshift(real(ifft2(Sf);fr3 = deconvwnr(g, PSF, NCORR, ICORR); subplot(221); imshow(g, );title('degraded image');subplot(222); im