Размытие изображения с использованием фильтра Гаусса в Matlab без аддитивного шума

Я должен использовать обратный фильтр, чтобы удалить размытие с этого изображения.

.

К сожалению, мне нужно вычислить передаточную функцию H системы визуализации, используемой для получения этих более четких изображений. Она должна быть гауссовой. Итак, я должен определить приблизительную ширину Гаусса, пробуя разные ширины Гаусса в обратном фильтре и оценивая, какие результирующие изображения выглядят «лучшими».

Наилучший результат будет оптимально четким, т. е. края будут выглядеть четкими, но не будут иметь видимого звона.

Я пробовал, используя 3 подхода:

1. Я создал передаточную функцию с N измерениями (нечетное число для простоты), создав сетку N измерений, а затем применив к этой сетке функцию Гаусса. После этого мы добавляем к этой передаточной функции нули, чтобы получить тот же размер, что и исходное изображение. Однако после применения фильтра к исходному изображению я вижу только шум (слишком много артефактов).
2. Я создал передаточную функцию размером с исходное изображение, создав сетку того же размера, что и исходное изображение. Если sigma слишком мало, то величина PSF FFT широкая. В противном случае он становится тоньше. Если sigma мало, то изображение будет еще более размытым, но если мы установим очень высокое значение sigma, то получим такое же изображение (совсем не лучше).
3. Я использовал функцию fspecial, играя с размерами sigma и h. Но все же у меня не получается ничего более резкого, чем исходное размытое изображение.

Любые идеи?

Вот код, используемый для создания передаточной функции в подходе 1:

%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')