在matlab中使用高斯滤波器进行图像去模糊，没有加性噪声

Posted 2023-04-14

【中文标题】在matlab中使用高斯滤波器进行图像去模糊，没有加性噪声

【英文标题】：Image deblurring using gaussian filter in matlab without additive noise

【发布时间】：2017-07-30 10:13:01

【问题描述】：

我必须使用逆滤镜来去除这张图片的模糊效果

。

不幸的是，我必须弄清楚成像的传递函数H 用于获得这些更清晰图像的系统，它应该是高斯的。因此，我应该通过在逆滤波器中尝试不同的高斯宽度并判断哪些结果图像看起来“最好”来确定高斯的近似宽度。

最好的结果是最清晰的——即边缘看起来很锐利，但不会有明显的振铃。

我尝试了 3 种方法：

我通过创建N 维度的网格，然后将高斯函数应用于此网格，创建了具有N 维度（奇数，为简单起见）的传递函数。之后，我们向这个传递函数添加零，以获得与原始图像相同的大小。但是，在将滤镜应用于原始图像后，我只看到了噪点（伪影太多）。 我通过创建与原始图像大小相同的网格来创建大小与原始图像一样高的传递函数。如果sigma 太小，则PSF FFT 幅度很宽。否则会变薄。如果sigma 很小，那么图像会更加模糊，但是如果我们将sigma 值设置得非常高，那么我们会得到相同的图像（一点也不更好）。 我使用了fspecial 函数，使用sigma 和h 的大小。但我仍然没有得到比原始模糊图像更清晰的东西。

有什么想法吗？

以下是方法一中用于创建传递函数的代码：

%Create Gaussian Filter
function h = transfer_function(N, sigma, I) %N is the dimension of the kernel
%create a 2D-grid that is the same size as the Gaussian filter matrix
grid = -floor(N/2) : floor(N/2);
[x, y] = meshgrid(grid, grid);
arg = -(x.*x + y.*y)/(2*sigma*sigma);
h = exp(arg); %gaussian 2D-function
kernel = h/sum(h(:)); %Normalize so that total weight equals 1

[rows,cols] = size(I);
add_zeros_w = (rows - N)/2;
add_zeros_h = (cols - N)/2;

h = padarray(kernel,[add_zeros_w  add_zeros_h],0,'both'); % h = kernel_final_matrix

end 

这是每种方法的代码：

I = imread('lena_blur.jpg');
I1 = rgb2gray(I);
figure(1),
I1 = double(I1);
%---------------Approach 1
% N = 5; %Dimension Assume is an odd number
% sigma = 20; %The bigger number, the thinner the PSF in FREQ
% H = transfer_function(N, sigma, I1);
%I1=I1(2:end,2:end); %To simplify operations
imagesc(I1); colormap('gray'); title('Original Blurred Image')

I_fft = fftshift(fft2(I1)); %Shift the image in Fourier domain to let its DC part in the center of the image



% %FILTER-----------Approach 2---------------
% N = 5; %Dimension Assume is an odd number
% sigma = 20; %The bigger number, the thinner the PSF in FREQ
% 
% 
% [x,y] = meshgrid(-size(I,2)/2:size(I,2)/2-1, -size(I,1)/2:size(I,1)/2-1);
% H = exp(-(x.^2+y.^2)*sigma/2);
% %// Normalize so that total area (sum of all weights) is 1
% H = H /sum(H(:));
% 
% %Avoid zero freqs
% for i = 1:size(I,2) %Cols
%     for j = 1:size(I,1) %Rows
%         if (H(i,j) == 0)
%             H(i,j) = 1e-8;
%         end
%     end
% end
% 
% [rows columns z] = size(I);
% G_filter_fft = fft2(H,rows,columns);
%FILTER---------------------------------


%Filter--------- Aproach 3------------
N = 21; %Dimension Assume is an odd number
sigma = 1.25; %The bigger number, the thinner the PSF in FREQ

H = fspecial('gaussian',N,sigma)
[rows columns z] = size(I);
G_filter_fft = fft2(H,rows,columns);

%Filter--------- Aproach 3------------



%DISPLAY FFT PSF MAGNITUDE
figure(2),
imshow(fftshift(abs(G_filter_fft)),[]); title('FFT PSF magnitude 2D');


% Yest = Y_blurred/Gaussian_Filter
I_restoration_fft = I_fft./G_filter_fft;
I_restoration = (ifft2(I_restoration_fft));
I_restoration = abs(I_restoration);




I_fft = abs(I_fft);

% Display of Frequency domain (To compare with the slides) 
figure(3),
subplot(1,3,1); 
imagesc(I_fft);colormap('gray');title('|DFT Blurred Image|')
subplot(1,3,2)
imshow(log(fftshift(abs(G_filter_fft))+1),[]) ;title('| Log DFT Point Spread Function + 1|');
subplot(1,3,3)
imagesc(abs(I_restoration_fft));colormap('gray'); title('|DFT Deblurred|')
% imshow(log(I_restoration+1),[])

%Display PSF FFT in 3D
figure(4)
hf_abs = abs(G_filter_fft);
%270x270
surf([-134:135]/135,[-134:135]/135,fftshift(hf_abs));
% surf([-134:134]/134,[-134:134]/134,fftshift(hf_abs));
shading interp, camlight, colormap jet
xlabel('PSF FFT magnitude')

%Display Result (it should be the de-blurred image)
figure(5),
%imshow(fftshift(I_restoration));
imagesc(I_restoration);colormap('gray'); title('Deblurred Image')

%Pseudo Inverse restoration
% cam_pinv = real(ifft2((abs(G_filter_fft) > 0.1).*I_fft./G_filter_fft));
% imshow(fftshift(cam_pinv));
% xlabel('pseudo-inverse restoration')

【问题讨论】：

Matlab网页中有解决方案：blogs.mathworks.com/steve/2007/08/13/…，你看看这些吗？？

看看deconwnr也可能很有趣

@monakons 是的，实际上我的解决方案部分基于该博客。问题是我根本没有噪音。图像只是模糊了（不是附加噪声）。

@m7913d 我没有任何噪音...

您可以将其用于无噪点图像，只需将 NSR 设为零即可。请注意，您总是有量化噪声。

【参考方案1】：

一个可能的解决方案是deconvwr。我将首先从未失真的 lena 图像开始展示它的性能。所以，我确切地知道高斯模糊函数。请注意，将estimated_nsr 设置为零会由于量化噪声而完全破坏性能。

 I_ori = imread('lenaTest3.jpg'); % Download an original undistorted lena file
 N = 19;
 sigma = 5;
 H = fspecial('gaussian',N,sigma) 
 estimated_nsr = 0.05;

 I = imfilter(I_ori, H)

 wnr3 = deconvwnr(I, H, estimated_nsr); 
 figure
 subplot(1, 4, 1);
 imshow(I_ori) 
 subplot(1, 4, 2);
 imshow(I) 
 subplot(1, 4, 3);
 imshow(wnr3) 
 title('Restoration of Blurred, Noisy Image Using Estimated NSR'); 
 subplot(1, 4, 4);
 imshow(H, []);

我通过反复试验为您的问题找到的最佳参数。

 N = 19; 
 sigma = 2; 
 H = fspecial('gaussian',N,sigma) 
 estimated_nsr = 0.05;

编辑：精确计算使用的模糊滤镜

如果你下载的是未失真的 lena (I_original_fft)，你可以计算使用的模糊滤镜如下：

G_filter_fft = I_fft./I_original_fft

【讨论】：

感谢您的回答。但我需要使用我在第一张图片中发布的这张确切的图片。正如我所写：我应该通过在 INVERSE 滤波器中尝试不同的高斯宽度并判断哪些结果图像看起来“最好”来确定高斯的近似宽度。这意味着如果我写： I_restoration_fft = I_fft./G_filter_fft;那么逆FFT应该是我想要的图像。我通过找到正确的高斯参数来实现这一点。 deconwnr 没用（我没有被要求处理噪音或使用维纳滤波器）

如上例所示。假设图像上有一些噪点是个好主意。如果你运行上面的例子，如果你将estimated_nsr 设置为零，你会得到非常糟糕的结果，即使高斯模糊滤波器是完全已知的。如果我正确地查看了您的代码，那么您基本上是在以零噪声实现 deconwnr。我的第二部分代码是稍微锐化图像的参数。

无论如何，如果你使用：I = imfilter(I_ori, H)，你就是在空间域中进行卷积。而我应该做的是频域中的逐点除法，即I_restoration_fft = I_fft./G_filter_fft，这与卷积un空间域不同。

请注意，我使用 imfilter 来模糊未失真的 lena 图像，因此我知道确切的模糊功能并且可以测试 deconvwnr。我用deconvwnr来恢复镜像，不用自己写一堆代码。其次，请注意您说您的方法基于this post，这与该帖子的 cmets 中讨论的 Wiener 方法非常相似。

哦，好的。好吧，您可能会得到一个不错的结果，因为您已经知道高斯（当您模糊图像时）并且您使用相同的高斯对模糊图像进行去模糊。这里的重点是猜测高斯的参数，因为你得到了一个模糊的图像（它已经被具有特定参数的高斯模糊；你必须猜测相同的参数）。所以，做我的客人，用我的形象尝试 deconwnr，和我在第一篇文章中发布的一样。你会看到 deconwnr 没有给出任何好的结果。