Thursday, April 4, 2013

Wiener Filter + Inverse Filter, Contd.

In the last post i derived the formula for wiener filter. I am not getting enough time to write something here. This night i decided to write something.  If you carefully examine the wiener filter formula it can be seen that when the K is zero ( that is no noise),it act just an inverse filter.
     This means we can just divide the Fourier transform of the input signal(degraded) with the Fourier transform of the degrade function. This will produce the original data. You may be wondering how this works. It is based on convolution theorem.

In my pc mathematica stopped working. So i am switching to matlab from now on

Following mat lab snippet proves that convolution in spatial domain and frequency domain is same

%% in this example convolution signal is not centered ( so the source)
%% by this way we can avoid the one shifting(ifftshift) at the end.

data1 = [1,2,3,2,1,0,0];
kernel = [1,1,1];
kernel_padded = [1,1,1,0,0,0,0];
'convolution spatial'
conv(data1,kernel,'valid')
d = fft(data1) .* fft(kernel_padded);
'fourier'
ifft(d)


When you run it this will produce the result  1  3    6   7  6  3  1. (conv function will out two more zeros, but it is not a part of input data)

Using this idea inverse filter works. Since inverse filter uses a division problem can happen when the denominator is zero or say too high. because of this it is not stable as itself also with noise.

Here is the results and code

a = imread('f:\\ip\\imtests\\exotic.jpg');
d = rgb2gray(a);
degrade = [.1 .1 .1;.1 .5 .1;.1 .1 .1];
dim = conv2(d,degrade,'valid');
figure('Name','degraded image with some filter')
imshow(dim,[0,255])

img_size = size(d) ;
img_size = size(dim); 
conv_mat = zeros(img_size(1),img_size(2));

% create conv matrix for inverse.
pw = floor(img_size(1))-3
ph = floor(img_size(2))-3
conv_mat = padarray(degrade,[pw,ph],'post');

fo = fft2(dim);
fc = fft2(conv_mat);
% the final division.
myimg = fo ./ (fc )

corrected = ifft2(myimg);
figure('Name','corrected with inverse filter')
imshow(real(corrected),[0 255])

Oringal Image



Corrected with Inverse

Don't expect this kind of correction all cases , because of that division things can go totally wrong . So one idea is to add some constant with denominator before division and play with it until you find some satisfactory results Or you need to check it against zero or non-numeric value.
Ok its time for me to stop. I just made this because of me :) . I was worried about not writing anything for sometime. That's why this quick post.  More interesting this are coming(to my mind)