imatest 中噪声SNR_BW 主要的噪声模型

转载

编程梦想编织者 2024-05-07 17:12:13

文章标签 imatest 中噪声SNR_BW 2d 图像复原维纳滤波 文章分类 架构后端开发

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);

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。