[Image Restoration, Wiener Filter] Ifft2 gives NaN complex matrix

Hello everyone

I am working on image restoration using Wiener Filtering, using the formula

from the book Digital Image Processing By Gonzales, Rafael C.Woods, Richard E ISBN-10: 9780133356724

However, I am getting NaN complex matrix when using ifft2 to get the restored image. Where did I do wrong?

img = imread('daisy.jpg');
img_values = double(rgb2gray(img));
[n1, n2] = size(img_values);

% Calculate the DFT of the image
f = fftshift(fft2(img_values));

% Get the degradation function
h1 = turbulence(0.01, n1,n2);
h1(h1 == 0) = 0.001; %Work around to avoid 0 values

%Apply degradation
f1 = f .* h1;

% Find the inverse FT
restore_img = ifft2(f1);
spectrum = abs(restore_img)/255;

% Generate Gaussian noise
mu= 0.1;
sigma = 0.15;
noise = normrnd(mu,sigma,size(img_values));
noise(noise == 0) = 0.001; %Work around to avoid 0 values

% Add noise to the spatial domain of the image
blurred_noisy_image = spectrum + noise;

% Display the images
subplot(1,3,1); imshow(img_values/255); title('original image');
subplot(1,3,2); imshow(spectrum); title('blurred image');
subplot(1,3,3);imshow(blurred_noisy_image); title('blurred and noisy image');


% Blurred and noisy image in Frequency domain
blur_noise_freq = fftshift(fft2(blurred_noisy_image));
% Wiener Filter

% H(u,v) = h1

%Calculate |H(u,v)|^2
h2 = h1.^2;

%Calculate Sn(u,V)
noise_freq = fft2(noise);
sn = abs(noise);

%Calculate Sf(u,v)
sf = abs(f);

wiener = (1./h1) .* (h2./(h2 + sn./sf));
result = wiener .* blur_noise_freq;
restored = ifft2(result); % This gives NaN Complex Matrix
figure;
imshow(abs(restored)/255);

The atmostpheric turbulence:

function h = turbulence(k, n1, n2)
[k1, k2] = meshgrid(-n2/2+1:n2/2, -n1/2+1:n1/2);
d = (k1.^2 + k2.^2);
h = exp(-k  .* d );
end

Any help would be greatly appreciated, please provide details explains as well.

Thanks in advance.

Answer 1

[Bjorn Gustavsson]

Just to get up to speed:

that then is F-transformed to:

which we can "solve" as:

or after multiplying the RHS with the complex conjugate of :

In order to avoid division-by-zero we add something like your to the RHS denominator:

which simplifies to:

I see no benefit from your last equation, what is gained by dividing by the complex-valued H?

I'd change this bit of your code:

%Calculate |H(u,v)|^2
h2 = abs(h1).^2; % I think this was the main error
%Calculate Sn(u,V)
noise_freq = fft2(noise);
sn = abs(noise);
%Calculate Sf(u,v)
sf = abs(f);
wiener = (conj(h1)./(h2 + sn./sf)); % simpler version than your

HTH