MATLAB Answers

How to build a wiener deconvolution manually? Or how to control the K value using wiener2 if possible?

27 views (last 30 days)
Joan B
Joan B on 29 Nov 2019
Commented: Joan B on 17 Dec 2019
Hi,
I'm trying to do Wiener deconvolution on a 256x256 image, where I have to be able to specify the K value.
g = noisy image matrix
G = fftshift(fft2(g))
x = clean/original image
for u = 1:256 % width of image
for v = 1:256 % height of image
M = 256;
N = 256;
D0 = 4; % pixel width is 4
D(u,v) = (((u - M./2).^2) + ((v - N./2).^2)).^(1./2);
H(u,v) = exp(-((D(u,v)).^2)./(2.*(D0.^2))); % Gaussian deblurring function
end
end
H2 = H*H; % is this the correct way to do |H|^2 with a matrix?
for k = -20:20 % trying different k values to see what works best (will try a larger range but keeping it small until I get the code working)
K = k.*(ones(256)); % matrix of K values, where K = |N|^2/|F|^2
F = (H2*G)./(H*(H2 + K)); % I don't think this is correct with the '.' but I get a singular matrix when I don't use that
f = ifft2(F);
difference = abs(f - x); % finding difference between element value in f and in x
S = (sum(difference,"all"))./(256.*256); % averaging the difference
good = min(S); % trying to find the k value which gives the minimum average difference
end
Any help is really appreciated.

  0 Comments

Sign in to comment.

Answers (1)

Raunak Gupta
Raunak Gupta on 4 Dec 2019
Hi,
You may use deconvwnr as it’s used for deconvolution (deblurring) where the psf can be the gaussian point spread function which can be computed using fspecial. The k value that is intended to be included can be done by setting the nsr value which stands for noise to signal power ratio.
wiener2 is locally adaptive filter that estimate noise mean and variance in a neighborhood of a pixel so including a global k value is not possible.
As per the manual implementation I see that after calculating the H matrix the frequency domain equivalent is missing so there it need some changes. Also I would suggest using H2 = H.*H where H is frequency domain filter because H*H return something different then the required |H(u,v)|^2. Also while looping over the k values it should not be negative as it represent noise to signal power ratio or even greater than 1 as then reconstructed image will not contain meaningful information.
Hope it helps.

  1 Comment

Joan B
Joan B on 17 Dec 2019
Hi Raunak,
Thank you for your help. Sorry for the late reply.
I should have explained better: I was required to try k values from -2000 to +2000 in increments of 1, and then use results from that and decide on a method to choose the best k value from those. That was why I was running it as a for loop and storing values for the difference between the original/clean and filtered images.
Anyway I am no longer working on that so no matter; thanks again for your reply!
Joan

Sign in to comment.

Sign in to answer this question.