MATLAB Answers

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

조회 수: 3(최근 30일)
Khai Hoang
Khai Hoang 2021년 5월 6일
답변: Bjorn Gustavsson 2021년 5월 6일
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.

채택된 답변

Bjorn Gustavsson
Bjorn Gustavsson 2021년 5월 6일
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

추가 답변(0개)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by