I_o = imread('blur.jpg'); %图像文件读入
I_c1 = I_o(:,:,1);
I_c2 = I_o(:,:,2);
I_c3 = I_o(:,:,3); %将图像转为单层
H = fspecial('gaussian'); %高斯低通滤波作为退化函数
I_c1_0 = im2double(I_c1);
I_c2_0 = im2double(I_c2);
I_c3_0 = im2double(I_c3); %将unit8转为double
IH_c1 = imfilter(I_c1_0, H, 'conv');
IH_c2 = imfilter(I_c2_0, H, 'conv');
IH_c3 = imfilter(I_c3_0, H, 'conv'); %通过高斯滤波退化
IH(:,:,1) = IH_c1;
IH(:,:,2) = IH_c2;
IH(:,:,3) = IH_c3; %导出退化后图片
G_c1 = imnoise(IH_c1,'gaussian');
G_c2 = imnoise(IH_c2,'gaussian');
G_c3 = imnoise(IH_c3,'gaussian'); %加噪声,参数为默认参数0,0.01
G(:,:,1) = G_c1;
G(:,:,2) = G_c2;
G(:,:,3) = G_c3; %导出加噪后图片
[im_row, im_line, im_floor] = size(I_o); %读取图片尺寸
%基础参数设定OK
%--------------------------------------------------------------------------
%直接反变换法
I_c1_fft2 = fft2(I_c1_0); %基础图像傅里叶变换
G_c1_fft2 = fft2(G_c1); %加噪图像傅立叶变换
H_fft2 = fft2(H, im_row, im_line); %高斯滤波器傅里叶变换
F_c1_fft2 = G_c1_fft2 ./ H_fft2; %直接除掉,得到还原图像的频谱
F_c1 = ifft2(F_c1_fft2); %反傅里叶变换得到还原图像
I_c2_fft2 = fft2(I_c2_0);
G_c2_fft2 = fft2(G_c2);
F_c2_fft2 = G_c2_fft2 ./ H_fft2;
F_c2 = ifft2(F_c2_fft2);
I_c3_fft2 = fft2(I_c3_0);
G_c3_fft2 = fft2(G_c3);
F_c3_fft2 = G_c3_fft2 ./ H_fft2;
F_c3 = ifft2(F_c3_fft2);
F(:,:,1) = F_c1;
F(:,:,2) = F_c2;
F(:,:,3) = F_c3;
%--------------------------------------------------------------------------
%维纳滤波法
Noise_c1 = G_c1 - IH_c1; %噪声函数
Noise_c1_fft2 = fft2(Noise_c1); %噪声傅里叶变换
Noise_c1_Square = abs(Noise_c1_fft2).^2; %噪声功率
I_c1_Square = abs(I_c1_fft2).^2; %原图功率
H_Square = abs(H_fft2).^2; %滤波器的平方
F_c1_fftwiener = (H_Square .* G_c1_fft2) ./ (H_fft2 .* (H_Square + Noise_c1_Square ./ I_c1_Square));
F_c1_wiener = ifft2(F_c1_fftwiener); %反变换得到还原图像
Noise_c2 = G_c2 - IH_c2;
Noise_c2_fft2 = fft2(Noise_c2);
Noise_c2_Square = abs(Noise_c2_fft2).^2;
I_c2_Square = abs(I_c2_fft2).^2;
F_c2_fftwiener = (H_Square .* G_c2_fft2) ./ (H_fft2 .* (H_Square + Noise_c2_Square ./ I_c2_Square));
F_c2_wiener = ifft2(F_c2_fftwiener);
Noise_c3 = G_c3 - IH_c3;
Noise_c3_fft2 = fft2(Noise_c3);
Noise_c3_Square = abs(Noise_c3_fft2).^2;
I_c3_Square = abs(I_c3_fft2).^2;
F_c3_fftwiener = (H_Square .* G_c3_fft2) ./ (H_fft2 .* (H_Square + Noise_c3_Square ./ I_c3_Square));
F_c3_wiener = ifft2(F_c3_fftwiener);
F_wiener(:,:,1) = F_c1_wiener;
F_wiener(:,:,2) = F_c2_wiener;
F_wiener(:,:,3) = F_c3_wiener;
%--------------------------------------------------------------------------
%L_R迭代算法
f_c1_lr2 = G_c1;
f_c2_lr2 = G_c2;
f_c3_lr2 = G_c3;
for t = 1 : 2 %2次L_R迭代
f_c1_lr2 = f_c1_lr2 .* (imfilter((G_c1 ./ imfilter(f_c1_lr2, H, 'conv')), H, 'corr'));
f_c2_lr2 = f_c2_lr2 .* (imfilter((G_c2 ./ imfilter(f_c2_lr2, H, 'conv')), H, 'corr'));
f_c3_lr2 = f_c3_lr2 .* (imfilter((G_c3 ./ imfilter(f_c3_lr2, H, 'conv')), H, 'corr'));
end;
f_lr2(:,:,1) = f_c1_lr2;
f_lr2(:,:,2) = f_c2_lr2;
f_lr2(:,:,3) = f_c3_lr2;
f_c1_lr5 = G_c1;
f_c2_lr5 = G_c2;
f_c3_lr5 = G_c3;
for s = 1 : 5 %5次L_R迭代
f_c1_lr5 = f_c1_lr5 .* (imfilter((G_c1 ./ imfilter(f_c1_lr5, H, 'conv')), H, 'corr'));
f_c2_lr5 = f_c2_lr5 .* (imfilter((G_c2 ./ imfilter(f_c2_lr5, H, 'conv')), H, 'corr'));
f_c3_lr5 = f_c3_lr5 .* (imfilter((G_c3 ./ imfilter(f_c3_lr5, H, 'conv')), H, 'corr'));
end;
f_lr5(:,:,1) = f_c1_lr5;
f_lr5(:,:,2) = f_c2_lr5;
f_lr5(:,:,3) = f_c3_lr5;
%--------------------------------------------------------------------------
%图像显示并比较
subplot(2,4,1),imshow(im2uint8(I_o)); %原图
title('原图');
subplot(2,4,2),imshow(im2uint8(IH)); %滤波器后图像
title('滤波后图像');
subplot(2,4,3),imshow(im2uint8(G)); %加噪声后图像
title('加噪声后图像');
subplot(2,4,4),imshow(im2uint8(F)); %直接逆滤波
title('直接逆滤波');
subplot(2,4,5),imshow(im2uint8(F_wiener)); %维纳滤波
title('维纳滤波');
subplot(2,4,6),imshow(im2uint8(f_lr2)); %L_R迭代2次
title('LR迭代2次');
subplot(2,4,7),imshow(im2uint8(f_lr5)); %L_R迭代5次
title('LR迭代5次');
%--------------------------------------------------------------------------
%灰度平均梯度法
GMG_F = 0; GMG_wiener = 0; GMG_lr = 0;
for k = 1 : 1 : im_row - 1 ;
for j = 1 : 1 : im_line - 1 ;
GMG_F = GMG_F + sqrt(((F(k+1,j)-F(k,j))^2+(F(k,j+1)-F(k,j))^2)/2); %直接逆滤波的GMG
GMG_wiener = GMG_wiener + sqrt(((F_wiener(k+1,j)-F_wiener(k,j))^2+(F_wiener(k,j+1)-F_wiener(k,j))^2)/2); %维纳滤波法的GMG
GMG_lr = GMG_lr + sqrt(((f_lr2(k+1,j)-f_lr2(k,j))^2+(f_lr2(k,j+1)-f_lr2(k,j))^2)/2); %L_R迭代2次算法的GMG
end;
end;
GMG_F = GMG_F / (im_row - 1) / (im_line - 1);
GMG_wiener = GMG_wiener / (im_row - 1) / (im_line - 1);
GMG_lr = GMG_lr / (im_row - 1) / (im_line - 1);