44 views (last 30 days)

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.

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.

Opportunities for recent engineering grads.

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

Start Hunting!
## 0 Comments

Sign in to comment.