6. 图像复原技术
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);
