Local variance estimation for image filtering

Hi,

while implementing median filter, medfilt2(), why does medfilt2(I.^2) and medfilt2(I).^2 produce same result?

I am interested in calculating local variance , but using this methodology it is exactly 0.

How do I implement the formula correctly?

Answer 1

[Bjorn Gustavsson]

You get that result because medfilt2 will select the median value of the pixel-intensities of your region, and completely discard all the other values. Therefore when you square that median value you get medfilt2(I).^2. Since x^2 is a monotnous function the element that is median of x, will also be the median element of x.^2, and that's why you always get zero's. Simple illustrating example:

x = [1 3 4 412 -3]
x2 = x.^2
Q0 = median(x)
Q1 = median(x).^2
Q2 = median(x2)

Have a look at the code of wiener2 to see how it is done with local averaging filtering. If you want a different variance-type operation you could try mad instead of std (you have to figure out how you want to handle the conversion from mad similar to the conversion between std and var). Or you could try something like:

localMedian = medfilt2(x,nhood,'symmetric');
localMEDVar = medfilt2((x-localMedian).^2,nhood,'symmetric');
% or:
localMADVar = (filter2(ones(nhood), abs(x-localMedian) ) / prod(nhood)).^2;

You have so many options to calculate some measur of the local intensity variations. Which one suits your needs best is difficult to judge. At least you now know why you get the all-zeros result from your first-stab.

HTH

[DGM]

IPT has stdfilt(), if that would be of use...

[Albert Tagtum]

Thanks for your answer, it was helpful.

The thing I am trying to do is to modify the wiener2 filter so that instead of mean value I put median.

In original code (wiener2), lines

box_filter = fspecial ("average", N);
mean_im = convn (I, box_filter, "same");
variance_im = convn (I.^2, box_filter, "same") - mean_im.^ 2;

should be modified, but since median filter is not convolutional I have trouble with implementation. stdfilt() doesn't produce same output when I substitute it in third line. If anyone has some input, it would be of great help.

[Bjorn Gustavsson]

You got a couple of suggestions of how to modify the local variance estimates to somethings based on medians or me(di)an absolute deviations above:

localMedian = medfilt2(x,nhood,'symmetric');
localMEDVar = medfilt2((x-localMedian).^2,nhood,'symmetric');
% or:
localMADVar = (filter2(ones(nhood), abs(x-localMedian) ) / prod(nhood)).^2;

If you can explain in what ways they fail to give you the characteristics of the local variations you want to capture, please let us know. Because then we can figure out ways to come up with further modifications.

[Albert Tagtum]

I have an image (or better say set of images) that is uniform in intesity but corrupted with noise. I am testing different filters to see which one is best in increasing signal-to-noise ratio (defined as mean2(Im)/std2(Im)).

I am specially interested in median modified wiener filter, but I don't know what the end result should be so I tried to modify the wiener2 code to see which modification produced same SNR value, but none of the above solutions produced same SNR as wiener2() (I get worse SNR).

When I change the line

variance_im = convn (I.^2, box_filter, "same") - mean_im.^ 2;

in wiener2 code with stdfilt(I) or one of your solutions and compare SNR from (modified) wiener2 and my median modified code (using exact modification) I get better SNR with median modified version. So I try to implement the variance in median modified code so that it has same formalism as wiener2 in order to compare them properly.

[Bjorn Gustavsson]

To compare different filters on a uniform intensity seems rather dubious to me. If we have a uniform intensity the best estimate of that intensity is the average of all pixels. I've always understood filtering to be a compromise between noise-reduction and reduction of spatial resolution. In order to compare how different methods handle that compromise it is necessary to use images with relevant spatial intensity-variations.

Then when you make these comparisons that local variance estimate of wiener2 is not that sacrosanct. With the median-filtering you obviously cannot make use of the "litteral" translation, since that leads to all zeros. A similar enough transition might be to use the basic definition of variance:

Sure it is not exactly the same as in wiener2, but ought to be close enough, and from that you can use the translation to median-filtering.

Then you have to make "a similar enough" estimate of the local intensity-variation to make the same type of supression of the filtered image. The good measure of noise-removal is to compare the intensity distribution of the difference between the filtered image and the noisy image. If you manage to make that fit with what you know about the noise-distribution you're in the good region for the noise-supression. Amongst all filters that give you that you want to find the filter that reduces the spatial resolution the least, but to judge that you need an image with relevant intensity-variations.

[Albert Tagtum]

Thank you Bjorn for very detailed answers and proposals of solutions!

Much appreciated.

[Bjorn Gustavsson]

My pleasure, happy that it helped.

Answer 2

[Image Analyst]

Try stdfilt if you have the Image Processing Toolbox. It gives the standard deviation in a local window.