How can I lowpass filter by only reducing peak data?

Question

I have a 2D image, which I want to lowpass filter, with these constraints / quality metrics:

1. I can not "add" light to the image, so each pixel in the result should be <= the corresponding pixel in the input.
2. The lowpass cutoff frequency should be a parameter, to experiment with
3. Applying this filter repeatedly should not change the result in a significant way.
4. The time it takes to run this algorithm (5 minutes for a 5MPix image seems reasonable)
5. Minimizing the amount of light that is filtered away.

Below are some approaches I have tried, together with their shortcomings:

1. Gaussian filter like normal, then pull the result down to comply with constraint 1. This complies very well with the first 3 points, but reduces a lot more light than necessary.

2. Fitting "upwards" parabolas through the "low" points and "downwards" parabolas between them to smooth out. This works great in 1D, but applying it first horizontally, then vertically produces bad results in 2D. It takes a lot longer, but not too long for my application. However, repeatably applying this filter will drastically change the result. If the (1D) input is a perfect "downward" parabola (which should not be filtered at all), it will be replaced by 2 "upward" parabolas sitting at the start/end.

3. Using some other form of 2D "basic" functions and linear solving to find the optimal parameters. This is an idea only currently, not implemented / tested yet.

My domain of experience in signal processing is almost exclusively image processing, so I hope to find alternatives to this problem with the input of experts active in other areas of signal processing.

update 2011/08/18

Based on the current reactions, I decided to make things a bit more clear by adding graphs of a typical input and the results of the 3 approaches I described originally + the suggestions I received so far. For easy comparison, I used just 1D filtering in these examples.

Input data:

Gaussian filter + bring it down to comply with requirement (1).
You can see that bringing it down results in unnecessary light reductions on the right hand side.

Parabolas
As far as I am concerned, this is pretty much excellent, sadly it does not translate perfectly into 2D by applying first horizontal, then vertical. In this case, you also see I can evaluate the fitted parabolas in floating point resolution, which is a small benefit, but not absolutely required.

Grayscale erosion
Based on the suggestion from rwong, I tried grayscale erosion. I used a structuring element with the same parabolic shape as my "fitted" parabolas. The result is almost exactly the same, so this looks promising. However, there are still a few problems: 1. My structuring element was not "big enough" (although it was already 801 pixels wide) 1. I only have "upwards" parabolas, no "downward parabolas to smooth the transition from one parabola to the next.

Median filtering
Only included for completeness, it is not really what I want.

raw data
I pasted the raw input data + the various python commands onto pastebin, so you can experiment with the same data too.
http://pastebin.com/ASnJ9M0p

Answer 1

There is indeed a 2D version for your attempt #2 - it is similar in theory, but it cannot be decomposed into two 1D operations. Please read about "2D grayscale morphological filtering". It is faster than curve fitting.

• http://opencv.willowgarage.com/documentation/image_filtering.html
• http://en.wikipedia.org/wiki/Dilation_(morphology)

Median filtering might also be useful if you are trying to remove speckles. A more advanced form of median filtering is "ordinal filtering".

In all cases, requirement #1 can be met trivially by taking the pixelwise minimum between the output and the input. It is an important quality criteria, but it will not limit the choice of algorithms.

---

Gaussian filtering (and a number of other useful filters) can be decomposed (first from 2D to 1D operations, then via Fourier transform), but there are many other useful image processing techniques that are not decomposable, which make them slow but does not diminish their usefulness.

Answer 2

I suggest using a smoothing spline.

Here's how you can do this using Matlab with the robust spline smoothing function SMOOTHN from the Matlab File Exchange (which contains the full source code, so that you can re-implement it somewhere else if needed). Note that it functions with n-dimensional data as well:

%# - get inputlist from pastebin

%# - smoothen data. Lower factor means less smooth
smoothingFactor = 1000;
smoothData = smoothn(inputlist,smoothingFactor);

%# - shift down
smoothData = smoothData - max(inputlist-smoothData);

%# - show results
plot(inputlist,'b'),hold on,plot(smoothData,'r')

Answer 3

One simple way to ensure that a filtered signal is always larger (or smaller) than the original is to use a moving maximum (or minimum respectively) window, followed by a non-negative filter with the same support and unity gain at DC. In other words, you can first perform gray-scale dilation with a flat structuring element, then filter the result using a regular blur kernel of the same size.

The dilation process will force all points within the support of the blur to be at least as large (or small) as the original point, so the blurred value will be at least (or at most) as large as the original point times the DC gain of the blur (which we set to unity by design). On the other hand, compared to a pure dilation the result will be smoother and typically closer to the original since the blur will smooth out the flat areas created by the dilation.

If implemented on the CPU, then in the 1D case the rectangular moving maximum/minimum window can be structured as a binary tree where the leaves act like a circular buffer, with $O(\log N)$ update cost per sample for a window size of N. In the case of a rectangular support the operation is separable. The blur can also be chosen to be separable.

Now, I have no idea what this actually looks like when applied to images, but the basic technique is rather well-known in the context of computing envelopes for brickwall lookahead limiting of audio and the extension to 2D is rather straight-forward. I know it's an old question, but figure it might still be fun.

Here is a (ridiculously inefficient; this is just to see how it looks) Shadertoy implementation: https://www.shadertoy.com/view/3ttcDM (including a variation that first does dilate+blur horizontally, then dilate + blur vertically)