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

24 ビュー (過去 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
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 );
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

その他の回答 (0 件)

Community Treasure Hunt

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

Start Hunting!

Translated by