Wiener Filtering: Mean Square Error Problems

Question

I am attempting to write a code to deconvolve an Instrument Response function, IRF (also known as a Point Spread Function, PSF), from an image, convolved with this unknown function, using the original 'true' image. The original and the convoluted image are both greyscale images which have been centred and scaled appropriately. The true image does not have noise but the other does which has lead me to Wiener filtering to try and remove noise whilst acquiring the IRF.

The Parameters are as such:

t: true image

c: convoluted image

T=fft2(t);
C=fft2(C);

K is scalar value, which is the noise-to-signal ratio (or 1/signal-to-noise, whichever way suits)

The formula I have found for deconvolution using this filter is:

L=conj(T)./(T*conj(T)+K)
IRF=L.*C

What I am attempting to do is estimate the signal-to-noise ratio (which is unknown) in order to find an ideal K by computing the forumla above and then calculating the mean square error for various values of K between 0.001 and 0.1. I have the following code:

K=linspace(0.001,0.1,100);
for i=1:length(K);
    L=conj(T)./(T*conj(T)+K(i));
    IRF=L.*C;
    Cv=uint8(ifft2(C));
    %Mean Square Error is Here:
    e=uint8(t)-Cv;      
    Ev(i)=mean(e(:))^2;
end

[minErrorValue minErrorPos]=min(Ev);
idealK=K(minErrorPos);
L=conj(T)./(T*conj(T)+idealK);
IRF=L.*C;
IRFfinal=ifft2(IRF);

I found the code for the mean square error online but I believe its wrong cause all I'm getting is a vector whose entries are all the same value. Plus I'm getting absolute nonsense in my final image which leads me to believe that the noise is still present within the data. I thought how it should be calculated would lead to a mean square error for each individual K(i). Is there another way to calculate this?

lennon310 lennon310 · Accepted Answer · 2014-02-21T03:37:08

In wiener deconvolution L=conj(H)./(H*conj(H)+K);, H = fft2(IRF) is the instrument response. And the deconvolution is used to estimate the 'true' image given the blurred image and PSF. I don't think you can use the formula to estimate the IRF;
All your Ev elements in your code are constant because they are just the difference between your true image and convoluted image. Trying to find out K doesn't make much sense, because what you need to compare is the estimated IRF and true IRF, and the latter is just unknown.

There are two possible methods to determine the blurred kernel.

The 1st one is to optimize the objective function:

min(sum([c(x,y)-IRF(x,y)*t(x,y)]^2) + lambda*sum(IRF(x,y)^2)); for each pixel (x,y), where regularization term is to constrain the norm of the blurred kernel.

You have:

gradient_pixel(x,y) = -[c(x,y) - IRF(x,y)*t(x,y)]*t(x,y) + lambda*IRF(x,y);
IRF(x,y) = IRF(x,y) - eta * gradient_pixel(x,y);

Iteratively estimate the IRF.

The 2nd method is to use the cross correlation, suppose the size of IRF is sIRF:

t1 = im2col(t, sIRF, 'sliding')';
IRF1 = t1 \ c(:);
IRF = rot90(reshape(IRF1, sIRF), 2);

This is pretty easy.

Wiener Filtering: Mean Square Error Problems

2 Answers