Re: median filter that preserves peaks

Posted by Kenneth Sloan-2 on
URL: http://imagej.273.s1.nabble.com/median-filter-that-preserves-peaks-tp5004693p5004717.html

Just so.

More generally (I guess):

a) produce a boolean mask identifying "good pixels" and "pixels in need of correction"
b) replace all the "pixels in need of correction" by a value computed from the "good pixels"

Again (I really think this is a point worth emphasizing) - this is the kind of freedom you have (use it, or abuse it) when you give up the idea that you are *only* doing linear signal processing, or that the proper definition of "what you are doing" is best stated in terms of the Fourier Domain.

Just, please note the "abuse it" clause.  Not everything that is possible is a good idea!

for example…one reason that people like Median instead of Mean is that Median *tends* to completely discount outliers, while Mean can be significantly affected by certain types of noise.  The Median computed from the "good pixels" is not terribly different than the Median of all the pixels (but the Mean might be).

If you implement the idea of a mask to identify the pixels that need to be changed AND the ones to NOT use in the computation, you can get this effect with nearly any filter.  The only rub is that it's not linear anymore.  To compensate you, note that you might get away with doing the filtering IN-PLACE (since you only change pixels that do not participate in the computation!)

I think someone might have a lot of fun with this idea.  Please steal it, and let me know if it gains any traction.
 --
Kenneth Sloan
[hidden email]


On Sep 7, 2013, at 12:04 , Gabriel Landini <[hidden email]> wrote:

> On Saturday 07 Sep 2013 06:19:18 Kenneth Sloan wrote:
>> The nice thing about filters like the MedianFilter is that once you have
>> committed to that style, you can implement ANY function of the values
>> inside the window.  This is what separates Median from more
>> signal-processing-inspired methods like Mean.  Median is a frank "spatial
>> domain" operation, with no clear partner in the Fourier domain.  Given a
>> 5x5 window, say, you have 25 values - and are free to compute "anything you
>> want to" from those 25 values.
>>
>> To completely solve this problem, I propose a small suite of programs, all
>> modified in roughly the same way.  The idea is that small, local
>> fluctuations are "noise" but macroscopic major features are "signal".
>>
>> The basic idea is: if the central value in a window is the MIN or MAX within
>> the window, LEAVE IT ALONE.  Otherwise - do the usual thing.
>>
>> So, XXXFilter becomes:
>>
>> for (x = 0, width)
>> for (y = 0, height)
>>   {
>>     gatherWindowAndCompute(x,y,windowW, windowH)
>>     if (MinOrMax(input(x,y)) output(x,y) = input(x,y)
>>     else output(x,y) =  XXX(window);
>>   }
>>
>> I hope that communicates the idea.
>>
>> Now…to fully implement what the OP wants, I suggest implementing
>> MedianFilter and also MeanFilter.
>
> Very interesting. I guess that this would be similar to:
> 1. compute the local minima and maxima for a kernel size and their grey
> values,
> 2. filter the original image with a median filter
> 3. replace, in the filtered image, the local minima and maxima of the
> original.
>
> This macro implements that:
> //------------------8<---------------
> // MedianLocalMinMax.txt
> // Based on the original idea of Kenneth Sloan
> // posted to the IJ mailing list on Sept 7, 2013.
> // G. Landini at bham. ac. uk
>
> a=getNumber("Local radius",3);
> setBatchMode(true);
> b=getTitle();
> setOption("BlackBackground", true);
> run("Duplicate...", "title=Median");
> run("Duplicate...", "title=Maxima");
> run("Maximum...", "radius="+a);
> run("Image Calculator...", "image1=Maxima operation=Subtract image2=Median  
> create");
> setThreshold(0, 0);
> run("Convert to Mask");
> selectWindow("Maxima");
> run("Close");
> selectWindow("Result of Maxima");
> rename("LocalMaxima");
>
> selectWindow("Median");
> run("Duplicate...", "title=Minima");
> run("Minimum...", "radius="+a);
> run("Image Calculator...", "image1=Median operation=Subtract image2=Minima
> create");
> setThreshold(0, 0);
> run("Convert to Mask");
> selectWindow("Minima");
> run("Close");
> selectWindow("Result of Median");
> rename("LocalMinima");
>
> imageCalculator("Add" , "LocalMaxima","LocalMinima");
> selectWindow("LocalMinima");
> close();
>
> selectWindow("Median");
> run("Median...", "radius="+a);
> imageCalculator("Subtract", "Median","LocalMaxima");
> imageCalculator("AND", "LocalMaxima", b);
> imageCalculator("Add", "Median","LocalMaxima");
>
> selectWindow("LocalMaxima");
> close();
>
> selectWindow("Median");
> rename("MedianLocalMinMax");
>
> setBatchMode(false);
> //------------------8<---------------
>
>
> Furthermore, it is also  possible to replace with the *regional* minima and
> maxima instead, as those are regions guaranteed to be surrounded by brighter
> and darker pixels respectively and more likely to be peaks and wells in the
> image.
> Interestingly the regional minima and maxima do not depend on kernel sizes.
>
> The following macro does that (the Domes plugin is part of the morphology
> collection, available in my page)
>
> //------------------8<---------------
> // MedianRegionalMinMax.txt
> // Based on the original idea of Kenneth Sloan
> // posted to the IJ mailing list on Sept 7, 2013.
> //  G. Landini at bham. ac. uk
>
> a=getNumber("Local radius",3);
> setBatchMode(true);
> b=getTitle();
> setOption("BlackBackground", true);
> run("Duplicate...", "title=Median");
> run("Duplicate...", "title=RMaxima");
> run("Domes ", "height=1");
> setThreshold(1, 255);
> run("Convert to Mask");
>
> selectWindow("Median");
> run("Duplicate...", "title=RMinima");
> run("Domes ", "height=1 basins");
> setThreshold(1, 255);
> run("Convert to Mask");
>
> imageCalculator("Add" , "RMaxima","RMinima");
> selectWindow("RMinima");
> close();
>
> selectWindow("Median");
> run("Median...", "radius="+a);
> imageCalculator("Subtract", "Median","RMaxima");
> imageCalculator("AND", "RMaxima", b);
> imageCalculator("Add", "Median","RMaxima");
>
> selectWindow("RMaxima");
> close();
>
> selectWindow("Median");
> rename("MedianRegionalMinMax");
>
> setBatchMode(false);
> //------------------8<---------------
>
>
> Hope it is useful. Let me know of any bugs.
> Cheers
>
> Gabriel
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html