dipexp6.m

这是一个维纳程序

% DIP experiment 6
% 逆滤波,维纳滤波
clc
clear all;

A = imread('Cameraman.bmp');
figure(1),subplot(221);imshow(A,[]);title('原图像');
figure(2),subplot(221);imshow(A,[]);title('原图像');
figure(3),subplot(221);imshow(A,[]);title('原图像');
A = double(A);
[M,N] = size(A);

%产生模糊冲击函数的模板,其中n为模板的大小
n = 7;
h = BlurTemplate(n);
sumh = sum(sum(h));
h = h/sumh;

%为了消除混叠效应引起的误差，需要对模板H及原图像的尺寸进行延拓
%对模板尺寸进行延拓
expandh = zeros(M+n-1,N+n-1);
expandh(1:n, 1:n) = h; 
%对图像尺寸进行延拓
expandA = zeros(M+n-1,N+n-1);
expandA(1:M, 1:N) = A; 


%在频域对图像进行模糊
%对延拓后的expandH进行离散傅立叶变换，为了提高运算效率，故调用matlab中自带的fft2函数
FH = fft2(expandh);
%对延拓后的图像进行离散傅立叶变换
FA = fft2(expandA);
%退化操作
FbA = FA .* FH;
%进行傅立叶反变换得到退化后的模糊图像
bA = abs(ifft2(FbA));
% figure(4),imshow(bA(1+ceil(n/2):M+floor(n/2),1+ceil(n/2):N+floor(n/2)),[]);
% title('与冲击函数卷积模糊后的图像')


%对模糊后的图像叠加高斯噪声
nbA1 = bA + sqrt(8)*randn([M+n-1,N+n-1]);%叠加均值为0，方差为8的高斯随机噪声
nbA2 = bA + sqrt(16)*randn([M+n-1,N+n-1]);%叠加均值为0，方差为16的高斯随机噪声
nbA3 = bA + sqrt(32)*randn([M+n-1,N+n-1]);%叠加均值为0，方差为32的高斯随机噪声

%画出模糊并添加高斯噪声后的的图像
figure(1),subplot(222);imshow(nbA1(1+floor(n/2):M+floor(n/2),1+floor(n/2):N+floor(n/2)),[]);title('模糊并添加均值为8高斯噪声后的图像');
figure(2),subplot(222);imshow(nbA2(1+floor(n/2):M+floor(n/2),1+floor(n/2):N+floor(n/2)),[]);title('模糊并添加均值为16高斯噪声后的图像');
figure(3),subplot(222);imshow(nbA3(1+floor(n/2):M+floor(n/2),1+floor(n/2):N+floor(n/2)),[]);title('模糊并添加均值为32高斯噪声后的图像');

%对模糊并添加噪声后的图像求离散傅立叶变换
FnbA1 = fft2(nbA1);
FnbA2 = fft2(nbA2);
FnbA3 = fft2(nbA3);

%********************************逆滤波********************************

%为了防止H(u,v)在UV平面上取0或很小且消除振铃效应，取恢复转移函数FrH(u,v)为
% 如果  H(u,v)<d,   FrH(u,v)=k
%       H(u,v)>=d,  FrH(u,v)=1/H(u,v)
% 其中d为小于1的常数，且选的较小为好
FrH = zeros(M+n-1,N+n-1);
for i = 1:M+n-1
    for j = 1:N+n-1
        if abs(FH(i,j)) < 0.1
            FrH(i,j) = 1;
        else
            FrH(i,j) = 1/FH(i,j);
        end
    end
end
%逆滤波过程
FrnbA1 = FnbA1 .* FrH;
FrnbA2 = FnbA2 .* FrH;
FrnbA3 = FnbA3 .* FrH;
%得到逆滤波复原后的图像
rnbA1 = abs(ifft2(FrnbA1));
rnbA2 = abs(ifft2(FrnbA2));
rnbA3 = abs(ifft2(FrnbA3));
%画出逆滤波复原后的图像
figure(1),subplot(223);imshow(rnbA1(1:M,1:N),[]);title('逆滤波复原的图像');
figure(2),subplot(223);imshow(rnbA2(1:M,1:N),[]);title('逆滤波复原的图像');
figure(3),subplot(223);imshow(rnbA3(1:M,1:N),[]);title('逆滤波复原的图像');


%***********************************维纳滤波**************************************

%为求得噪声的功率谱对噪声进行离散傅立叶变换
Fn1 = fft2(nbA1 - bA);
Fn2 = fft2(nbA2 - bA);
Fn3 = fft2(nbA3 - bA);
%求得维纳滤波器转移函数
W1 = conj(FH) ./ abs( FH.*conj(FH) + (Fn1.*conj(Fn1)) ./ (FA.*conj(FA)) );
W2 = conj(FH) ./ abs( FH.*conj(FH) + (Fn2.*conj(Fn2)) ./ (FA.*conj(FA)) );
W3 = conj(FH) ./ abs( FH.*conj(FH) + (Fn3.*conj(Fn3)) ./ (FA.*conj(FA)) );
%维纳滤波过程
FwrnbA1 =  FnbA1 .* W1;
FwrnbA2 =  FnbA2 .* W2;
FwrnbA3 =  FnbA3 .* W3;
%求得维纳滤波复原后的图像
wrnbA1 = abs(ifft2(FwrnbA1));
wrnbA2 = abs(ifft2(FwrnbA2));
wrnbA3 = abs(ifft2(FwrnbA3));
%画出维纳滤波复原后的图像
figure(1),subplot(224);imshow(wrnbA1(1:M,1:N),[]);title('维纳滤波复原的图像');
figure(2),subplot(224);imshow(wrnbA2(1:M,1:N),[]);title('维纳滤波复原的图像');
figure(3),subplot(224);imshow(wrnbA3(1:M,1:N),[]);title('维纳滤波复原的图像');