6.1 圖像噪聲模型
1 %--------噪聲的MATLAB實現---------- 2 3 %圖像添加噪聲 4 J=imnoise(I,type,parameters);% type:噪聲類型 5 6 %高斯噪聲 7 %通過均值和方差產生高斯噪聲 8 J=imnoise(I,'gaussian',m,v);%m為均值，默認值0；v為方差默認值0.01 9 %通過位置信息產生高斯噪聲 10 J=imnoise(I,'localvar',V); 11 %根據亮度值產生高斯噪聲 12 J=imnoise(I,'localvar',h,v);%h表示圖像亮度值；v表示與h向量對應的高斯噪聲方差 13 14 %椒鹽噪聲 15 J=imnoise(I,'salt & pepper',d);%噪聲的密度為d，噪聲占整個像素總數的百分比，默認0.05 16 17 %泊松噪聲 18 J=imnoise(I,'possion');%從數據中產生泊松噪聲，不是人工噪聲添加到數據中 19 20 %乘性噪聲 21 J=imnoise(I, 'speckle', v);%J=I*n*I;n是均值為0，方差為V的均勻分布的隨機噪聲，V默認值0.04 22 23 %均勻分布的噪聲 24 m=256; n=256; 25 a=50; 26 b=180; 27 I=a+(b-a)*rand(m,n); 28 figure; 29 subplot(121); imshow(uint8(I)); 30 subplot(122); imhist(uint8(I)); 31 32 %指數分布的噪聲 33 m=256; n=256; 34 a=0.04; 35 k=-1/a; 36 I=k*log(1-rand(m, n)); 37 figure; 38 subplot(121); imshow(uint8(I)); 39 subplot(122); imhist(uint8(I));
6.2 空域內的濾波復原
1 %--------均值濾波-------- 2 3 %算術均值和幾何均值 4 I=imread('cameraman.tif'); 5 I=im2double(I); 6 I=imnoise(I, 'gaussian', 0.05); %添加高斯噪聲 7 PSF=fspecial('average', 3); %產生PSF 8 J=imfilter(I, PSF); %算術均值濾波 9 K=exp(imfilter(log(I), PSF)); %集合均值濾波 10 figure; 11 subplot(131); imshow(I); 12 subplot(132); imshow(J); 13 subplot(133); imshow(K); 14 15 %逆諧波均值濾波 16 I=imread('cameraman.tif'); 17 I=im2double(I); 18 I=imnoise(I, 'salt & pepper', 0.01); 19 PSF=fspecial('average', 3); 20 Q1=1.6; 21 Q2=-1.6; 22 j1=imfilter(I.^(Q1+1), PSF); %Q為正，去除椒噪聲 23 j2=imfilter(I.^Q1, PSF); 24 J=j1./j2; 25 k1=imfilter(I.^(Q2+1), PSF); %Q為負，去除鹽噪聲 26 k2=imfilter(I.^Q2, PSF); 27 K=k1./k2; 28 figure; 29 subplot(131); imshow(I); 30 subplot(132); imshow(J); 31 subplot(133); imshow(K); 32 33 %--------順序統計濾波-------- 34 35 %中值濾波器能很好的保留圖片邊緣，非常適合去除椒鹽噪聲，效果優於均值濾波 36 %最大值濾波器也能去除椒鹽噪聲，但會從黑色物體的邊緣去除一些黑色像素 37 %最小值濾波器會從白色物體的邊緣去除一些白色像素 38 39 %二維中值濾波 40 J=medfilt2(I); %窗口大小默認3*3 41 J=medfilt2(I,[m,n]); %窗口大小m*n 42 43 %二維排序濾波 44 J=ordfilt2(I, order, domain); %對矩陣domain中的非零值進行排序，order為選擇的像素位置 45 46 %最大值/最小值濾波 47 I=imread('cameraman.tif'); 48 I=im2double(I); 49 I=imnoise(I, 'salt & pepper', 0.01); 50 J=ordfilt2(I, 1, ones(4,4)); 51 K=ordfilt2(I, 9, ones(3)); 52 53 %--------自適應濾波-------- 54 %wiener2()根據圖像的噪聲進行自適應維納濾波，可以對噪聲進行估計。根據圖像的局部方差來調整濾波器的輸出 55 J=wiener2(I,[m,n],noise); %窗口大小m*n，默認3*3；noise為噪聲的能量 56 [J,noise]=wiener2(I,[m,n]); 57 58 RGB=imread('saturn.png'); 59 I=rgb2gray(RGB); 60 J=imnoise(I, 'gaussian', 0, 0.03); 61 [K, noise]=wiener2(J, [5, 5]); %自適應濾波 62 figure; 63 subplot(121); imshow(J); 64 subplot(122); imshow(K);
6.3 圖像復原方法
逆濾波復原
1 %逆濾波復原 2 clear all; close all; 3 I=imread('cameraman.tif'); 4 I=im2double(I); 5 [m, n]=size(I); 6 M=2*m; n=2*n; 7 u=-m/2:m/2-1; 8 v=-n/2:n/2-1; 9 [U, V]=meshgrid(u, v); 10 D=sqrt(U.^2+V.^2); 11 D0=130; %截至頻率 12 H=exp(-(D.^2)./(2*(D0^2))); %高斯低通濾波器 13 N=0.01*ones(size(I,1), size(I,2)); 14 N=imnoise(N, 'gaussian', 0, 0.001); %添加噪聲 15 J=fftfilter(I, H)+N; %頻域濾波並加入噪聲 16 figure; 17 subplot(121); imshow(I); %顯示原始圖像 18 subplot(122); imshow(J, [ ]); %顯示退化后的圖像 19 HC=zeros(m, n); 20 M1=H>0.1; %頻率范圍 21 HC(M1)=1./H(M1); 22 K=fftfilter(J, HC); %逆濾波 23 HC=zeros(m, n); 24 M2=H>0.01; 25 HC(M2)=1./H(M2); 26 L=fftfilter(J, HC); %進行逆濾波 27 figure; 28 subplot(121); imshow(K, [ ]); %顯示結果 29 subplot(122); imshow(L, [ ]); %頻率范圍較大，會放大噪聲的影響；頻率范圍較小，則達不到去除模糊的效果 30 31 function Z=fftfilter(X, H) 32 F=fft2(X, size(H,1), size(H, 2)); 33 Z=H.*F; 34 Z=ifftshift(Z); 35 Z=abs(ifft2(Z)); 36 Z=Z(1:size(X, 1), 1:size(X, 2)); 37 end
維納濾波復原
1 J=deconvwnr(I, PSF, NSR); %PSF為點擴散函數，NSR為信噪比 2 J=deconvwnr(I, PSF, NCORR,ICORR); %NCORR為噪聲的自相關函數，ICORR為原始圖像的自相關函數 3 4 %維納濾波對運動模糊的圖像進行復原 5 clear all; close all; 6 I=imread('onion.png'); 7 I=rgb2gray(I); 8 I=im2double(I); 9 LEN=25; %參數設置 10 THETA=20; 11 PSF=fspecial('motion', LEN, THETA); %產生PSF 12 J=imfilter(I, PSF, 'conv', 'circular'); %運動模糊 13 NSR=0; 14 K=deconvwnr(J, PSF, NSR); %維納濾波復原 15 figure; 16 subplot(131); imshow(I); %原圖像 17 subplot(132); imshow(J); %退化圖像 18 subplot(133); imshow(K); %復原圖像 19 20 %維納濾波對含有噪聲的運動模糊圖像進行復原 21 clear all; close all; 22 I=imread('cameraman.tif'); 23 I=im2double(I); 24 LEN=21; 25 THETA=11; 26 PSF=fspecial('motion', LEN, THETA); 27 J=imfilter(I, PSF, 'conv', 'circular'); 28 noise_mean=0; 29 noise_var=0.0001; 30 K=imnoise(J, 'gaussian', noise_mean, noise_var); %添加高斯噪聲 31 figure; 32 subplot(121); imshow(I); 33 subplot(122); imshow(K); 34 NSR1=0; 35 L1=deconvwnr(K, PSF, NSR1); %維納濾波復原 36 NSR2=noise_var/var(I(:)); 37 L2=deconvwnr(K, PSF, NSR2); %圖像復原 38 figure; 39 subplot(121); imshow(L1); 40 subplot(122); imshow(L2); 41 42 %通過圖像自相關信息進行復原 43 clear all; close all; 44 I=imread('rice.png'); 45 I=im2double(I); 46 LEN=20; 47 THETA=10; 48 PSF=fspecial('motion', LEN, THETA); 49 J=imfilter(I, PSF, 'conv', 'circular'); 50 figure; 51 subplot(121); imshow(I); 52 subplot(122); imshow(J); 53 noise=0.03*randn(size(I)); 54 K=imadd(J, noise); %添加噪聲 55 NP=abs(fft2(noise)).^2; 56 NPower=sum(NP(:))/prod(size(noise)); 57 NCORR=fftshift(real(ifft2(NP))); %噪聲的自相關函數 58 IP=abs(fft2(I)).^2; 59 IPower=sum(IP(:))/prod(size(I)); 60 ICORR=fftshift(real(ifft2(IP))); %圖像的自相關函數 61 L=deconvwnr(K, PSF, NCORR, ICORR); %圖像復原 62 figure; 63 subplot(121); imshow(K); 64 subplot(122); imshow(L);
約束最小二乘法復原
1 J=deconvreg(I,PSF); %PSF為點擴散函數 2 J=deconvreg(I,PSF,NOISEPOWER); %NOISEPOWER為噪聲強度，默認為0 3 J=deconvreg(I,PSF,NOISEPOWER,LRANGE); %LRANGE拉格朗日搜索范圍 4 J=deconvreg(I,PSF,NOISEPOWER,LRANGE，REGOP); %REGOP為約束算子 5 6 %約束最小二乘法進行圖像復原 7 clear all; close all; 8 I=imread('rice.png'); 9 I=im2double(I); 10 PSF=fspecial('gaussian', 8, 4); %產生PSF 11 J=imfilter(I, PSF, 'conv'); %圖像退化 12 figure; 13 subplot(121); imshow(I); 14 subplot(122); imshow(J); 15 v=0.02; 16 K=imnoise(J, 'gaussian', 0, v); %添加噪聲 17 NP=v*prod(size(I)); 18 L=deconvreg(K, PSF, NP); %圖像復原 19 figure; 20 subplot(121); imshow(K); 21 subplot(122); imshow(L);
Lucy-Richardson復原
1 J=deconvlucy(I, PSF); 2 J=deconvlucy(I, PSF, NUMIT); %NUMIT為算法的重復次數 3 J=deconvlucy(I, PSF, NUMIT, DAMPAR); %DAMPAR為偏差閾值 4 5 clear all; close all; 6 I=imread('rice.png'); 7 I=im2double(I); 8 LEN=30; 9 THETA=20; 10 PSF=fspecial('motion', LEN, THETA); 11 J=imfilter(I, PSF, 'circular', 'conv'); 12 figure; 13 subplot(121); imshow(I); 14 subplot(122); imshow(J); 15 K=deconvlucy(J, PSF, 5); %復原，5次迭代 16 L=deconvlucy(J, PSF, 15); %復原，15次迭代 17 figure; 18 subplot(121); imshow(K); 19 subplot(122); imshow(L);
不需要預先知道PSF，而且可以對PSF進行估計
優點：對退化圖像無先驗知識的情況下，仍然能進行復原
1 [J, PSF]=deconvblind(I, INITPSF);%INITPSF為PSF的估計值，返回PSF為算法實際采用的PSF值 2 [J, PSF]=deconvblind(I, INITPSF,NUMIT); %NUMIT為算法的重復次數，默認10 3 [J, PSF]=deconvblind(I, INITPSF,NUMIT,DAMPAR); %DAMPAR為偏移閾值 4 [J, PSF]=deconvblind(I, INITPSF,NUMIT,DAMPAR,WEIGHT); %WEIGHT為像素的加權值，默認為原始圖像的數值 5 [J, PSF]=deconvblind(I, INITPSF,NUMIT,DAMPAR,WEIGHT,READOUT); %READOUT為噪聲矩陣 6 7 %對運動模糊圖像采用盲解卷積算法進行復原 8 clear all; close all; 9 I=imread('cameraman.tif'); 10 I=im2double(I); 11 LEN=20; %設置參數 12 THETA=20; 13 PSF=fspecial('motion', LEN, THETA); %產生PSF 14 J=imfilter(I, PSF, 'circular', 'conv'); %運動模糊 15 INITPSF=ones(size(PSF)); 16 [K, PSF2]=deconvblind(J, INITPSF, 30); %圖像復原 17 figure; 18 subplot(121); imshow(PSF, []); %顯示原PSF 19 subplot(122); imshow(PSF2, []); %顯示估計PSF 20 axis auto; 21 figure; 22 subplot(121); imshow(J); %顯示退化圖像 23 subplot(122); imshow(K); %顯示復原圖像 24 25 %對退化圖像采用盲解卷積算法進行復原 26 clear all; close all; 27 I=checkerboard(8); %產生圖像 28 PSF=fspecial('gaussian', 7, 10); 29 v=0.001; 30 J=imnoise(imfilter(I, PSF), 'gaussian', 0, v); 31 INITPSF=ones(size(PSF)); 32 WT=zeros(size(I)); 33 WT(5:end-4, 5:end-4)=1; 34 [K, PSF2]=deconvblind(J, INITPSF, 20, 10*sqrt(v), WT); 35 figure; 36 subplot(131); imshow(I); 37 subplot(132); imshow(J); 38 subplot(133); imshow(K);