MATLAB图像复原【包括(自适应)均值滤波，(自适应)中值滤波，逆滤波，维纳滤波】

均值滤波

clc;clear;close all;
img=imread("ckt-board-orig.tif");
figure;
subplot(231);imshow(img);xlabel('a.原图');

%加噪声，密度0.1高斯白噪声
imgNoise=imnoise(img,'gaussian',0,0.01);
subplot(232);imshow(imgNoise);xlabel('b.加噪声');
[m,n]=size(imgNoise);

%选择mask大小，为可调奇数3、5、7、9...
mask=3;
imgNoise=double(imgNoise);

% varNoise=var(imgNoise(:));%无法获得准确值，一般采用估计的方法得到。
varNoise=0.01;

%进行拓展:因为要求mask扫到图像中的每一个像素点，所以要在图像周边加上一圈厚度为（mask-1）/2的边界
imgEx=[255*ones((mask-1)/2,n+(mask-1));
    255*ones(m,(mask-1)/2),imgNoise,255*ones(m,(mask-1)/2);
    255*ones((mask-1)/2,n+(mask-1))];
temp=zeros(mask,mask);
arithMean=zeros(m,n);
geomeMean=zeros(m,n);
adaptMean=zeros(m,n);
for x=1+(mask-1)/2:m+(mask-1)/2
    for y=1+(mask-1)/2:n+(mask-1)/2
        temp=imgEx(x-(mask-1)/2:x+(mask-1)/2,y-(mask-1)/2:y+(mask-1)/2);
        varLocal=var(temp(:));
        arithMean(x-(mask-1)/2,y-(mask-1)/2)=mean(mean(temp));
        geomeMean(x-(mask-1)/2,y-(mask-1)/2)=(prod(prod(temp),2))^(1/(mask*mask));
        adaptMean(x-(mask-1)/2,y-(mask-1)/2)=imgEx(x,y)-varNoise/varLocal*(imgEx(x,y)-arithMean(x-(mask-1)/2,y-(mask-1)/2));
    end
end
subplot(234);imshow(uint8(arithMean));xlabel('c.算术平均');
subplot(235);imshow(uint8(geomeMean));xlabel('d几何平均');
subplot(236);imshow(uint8(adaptMean));xlabel('e.自适应平均');

中值滤波

clc;clear;close all;
img=imread('ckt-board-orig.tif');
figure;
subplot(221);imshow(img);xlabel('a.原图');

%加噪声 密度为0.4椒盐噪声
imNoise=imnoise(img,'salt & pepper',0.4);
subplot(222);imshow(imNoise);xlabel('b.加噪');

%两次普通中值滤波
imMed=medfilt2(imNoise,[3 3]);
imMed=medfilt2(imMed,[3 3]);
subplot(223);imshow(imMed);xlabel('c.两次中值滤波');

%最大尺寸9*9的自适应中值滤波器
[w,h]=size(imNoise);
nmin=3;nmax=9; % 起始窗nmin*nmin,最大窗nmax*nmax
imAdaptMed=imNoise;% 定义复原后图像

%将iNoise扩展
imgEx=[zeros((nmax-1)/2,h+(nmax-1));
    zeros(w,(nmax-1)/2),imNoise,zeros(w,(nmax-1)/2);
    zeros((nmax-1)/2,h+(nmax-1))];
for x=1:w
    for y=1:h
        for n=nmin:2:nmax
            %iNoise中某点（x,y）的邻域
            Sxy=imgEx(x+(nmax-1)/2-(n-1)/2:x+(nmax-1)/2+(n-1)/2,y+(nmax-1)/2-(n-1)/2:y+(nmax-1)/2+(n-1)/2);
            Smax=max(max(Sxy));
            Smin=min(min(Sxy));
            Smed=median(median(Sxy));
            if Smed>Smin && Smed<Smax
                if imNoise(x,y)<=Smin || imNoise(x,y)>=Smax
                    imAdaptMed(x,y)=Smed;
                end
                break;
            end
        end
        imAdaptMed(x,y)=Smed; % 达到最大窗口都没满足上述条件姑且取最大窗中值
    end
end
subplot(224);imshow(imAdaptMed);xlabel('d.自适应中值滤波')
# 中值滤波逆滤波和维纳滤波

直接逆滤波和维纳滤波

clc;clear;close all;
% % %手打函数
% img=imread('original_DIP.tif');
% img=imread('aerial_view_no_turb.tif');
% img=imread('lena256.bmp');
% figure;
% subplot(231);imshow(img);xlabel('a.原图');
% f=double(img);
% F=fftshift(fft2(img));%原图傅里叶
% subplot(234);imshow(abs(F),[]);xlabel('b.傅里叶变换所的频谱');
% [m,n]=size(F);
% H=zeros(m,n);%传递函数
% B=zeros(m,n);%模糊傅里叶
% F1=zeros(m,n); % 逆滤波
% F2=zeros(m,n); % 维纳滤波
% a=0.1;b=0.1;T=1;K=0.00259;
% for u=1:m
%     for v=1:n
%         H(u,v)=(T/(pi*(a*u+b*v)))*sin(pi*(a*u+b*v))*exp(-1i*pi*(a*u+b*v));%自己随意设置的传递函数
%         B(u,v)=H(u,v)*F(u,v);
%     end
% end
% blur=ifft2(ifftshift(B));%模糊
% blur=256*blur/max(max(blur));
% blur=uint8(real(blur));
% subplot(232);imshow(blur);xlabel('c.傅里叶反变换所得模糊图像');
% 
% %模糊加高斯白噪声
% g=imnoise(blur,'gaussian',0,0.01);
% subplot(235);imshow(abs(g));xlabel('d.模糊加噪图像')
% g=double(g);
% G=fftshift(fft2(g));%模糊加高斯白噪声傅里叶
% 
% for u=1:m
%     for v=1:n
%         F1(u,v)=1/H(u,v)*G(u,v);
%         F2(u,v)=1/H(u,v)*(abs(H(u,v)))^2/((abs(H(u,v)))^2+K)*G(u,v);
%     end
% end
% f1=ifft2(ifftshift(F1));
% f1=256.*f1./max(max(f1));
% f1=uint8(real(f1));
% subplot(233);imshow(f1);xlabel('e.逆滤波图像');
% 
% f2=ifft2(ifftshift(F2));
% f2=256.*f2./max(max(f2));
% f2=uint8(real(f2));
% subplot(236);imshow(f2);xlabel('e.维纳滤波图像');

% % %系统自带函数
% img = imread('original_DIP.tif');
% img = im2double(img);%为什么要double一下？
% figure;
% subplot(231);
% imshow(img);
% title('Original image');
% 
% % % 函数说明
% % ①FSPECIAL('motion',LEN,THETA)为运动模糊算子，表示摄像物体逆时针方向以
% % theta 角度运动了 len 个像素， len 的默认值为 9， theta 的默认值为 0；% ②imfilter(f, w, filtering_mode, boundary_options, size_options)， f 为输入
% % 图像， w 为滤波掩模，
% % filtering_mode 用于指定在滤波过程中是使用“相关”‘corr’还是“卷积”‘conv’，
% % boundary_options 用于处理边界充零问题，边界的大小由滤波器的大小确定
% % (‘circular’图像大小通过将图像看成是一个二维周期函数的一个周期来扩展)
% 
% %模糊图像
% PSF = fspecial('motion', 50, -45);
% img1 = imfilter(img, PSF, 'conv', 'circular');
% subplot(232);
% imshow(img1);
% title('Blurred image');
% % 添加高斯噪声 noise_var 可变，按照书上分别设为 650,65,0.0065
% noise_var = 0.0065;
% img2 = imnoise(img1, 'gaussian', 0, noise_var);
% subplot(233);
% imshow(img2);
% title(['add Gaussian noise with variance is ',num2str(noise_var)]);
% 
% % % 函数说明
% % deconvwnr(I，PSF，NSR)其中，I值退化的图像，是原图像卷积一个点扩散函数 PSF
% % 然后加上加性噪声而得到的;
% % NSR噪信比，可以为一标量，也可以为与同样大小的矩阵，默认值为 0。
% 
% % 逆滤波， NSR设为0(噪信比为0，参数维纳滤波退化为理想逆滤波)
% % Specifying 0 for the NSR is equivalent to creating an ideal inverse filter.
% img3 = deconvwnr(img2, PSF, 0.0);
% subplot(223);
% imshow(img3);
% title('Result of invense filtering');
% % 参数维纳滤波，计算噪信比带入
% estimated_NSR = noise_var / var(img(:));
% img4 = deconvwnr(img2, PSF, estimated_NSR);
% subplot(224);
% imshow(img4);
% title('Result of Wiener filtering');

%%两种方式对比
img=imread('lena256.bmp');
img = im2double(img);
figure;
subplot(331);
imshow(img);
title('原图');
%模糊图像
PSF = fspecial('motion', 30, -45);
img1 = imfilter(img, PSF, 'conv', 'circular');
subplot(332);
imshow(img1);
title('原图模糊掉');
%添加高斯噪声
noise_var = 0.000001;
img2 = imnoise(img1, 'gaussian', 0, noise_var);
subplot(333);imshow(img2);title(['原图模糊掉并加上方差',num2str(noise_var)]);
%逆滤波
img3 = deconvwnr(img2, PSF, 0.0);
subplot(323);imshow(img3);title('matlab自带逆滤波');
%维纳滤波
estimated_NSR = noise_var / var(img(:));
img4 = deconvwnr(img2, PSF, estimated_NSR);
subplot(324);imshow(img4);title('matlab自带维纳滤波');
imwrite(img4,strcat('C:\Users\WangHao\Desktop\pic_here\matlab_WienerFiltering.jpg'),'jpg');
%手打逆滤波和维纳滤波
g=double(img1);
G=fftshift(fft2(g)); %模糊加高斯白噪声傅里叶
[m,n]=size(G);
H=zeros(m,n); %传递函数
F1=zeros(m,n); %逆滤波
F2=zeros(m,n); %维纳滤波
a=0.1;b=0.1;T=1;K=0.06;
for u=1:m
    for v=1:n
        H(u,v)=(T/(pi*(a*u+b*v)))*sin(pi*(a*u+b*v))*exp(-1i*pi*(a*u+b*v));
        F1(u,v)=1/H(u,v)*G(u,v);
        F2(u,v)=1/H(u,v)*(abs(H(u,v)))^2/((abs(H(u,v)))^2+K)*G(u,v);
    end
end
f1=ifft2(ifftshift(F1));
f1=256*f1/max(max(f1));
f1=uint8(real(f1));
subplot(325);imshow(f1);title('手打逆滤波');

f2=ifft2(ifftshift(F2));
f2=256*f2/max(max(f2));
f2=uint8(real(f2));
subplot(326);imshow(f2);title(['手打维纳滤波K=',num2str(K)]);
imwrite(f2,strcat('C:\Users\WangHao\Desktop\pic_here\Mine_WienerFiltering.jpg'),'jpg');

参考链接与图片

链接: link1.
link2.

图片: