# How can I lowpass filter by only reducing peak data?

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

• Can you explain a little more about restrictions 1 and 5? They appear to be (at first sight) contradictory.
– Peter K.
Aug 17 '11 at 10:00
• I am probably misunderstanding what you mean by "this algorithm", but 5 minutes for 5 MP seem like a lot for applying a lowpass filter. Aug 17 '11 at 11:59

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.

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.

• Hi, thanks for the pointer to grayscale morphological filtering. The description on wikipedia seems interesting and I will investigate that. However, in your link to the OpenCV documentation, I only see normal morphological filters, not grayscale ones. I will definately check this option and let you know the results. Thanks. Aug 17 '11 at 7:50
• Does rwong's suggestion of median filtering help at all? Explaining a little more about what you're trying to achieve by presenting a simple example of the data and a "fake" example of what you want to get out might help.
– Peter K.
Aug 17 '11 at 10:03
• I updated my question with sample data + results from various suggestions. I hope things are more clear now. Aug 18 '11 at 13:36

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')


• Thanks for your suggestion, I will investigate it. From your graph, it would seem like I need a much higher smoothingFactor than your example. The steep edge around x=700 is not removed and will be clearly visible. Also the initial bump in x = [0 ,400] is not removed at all. Don't you think this will have the same problem as any other (low-pass filter + move down) approaches? You can see the global offset between the two graphs, which will probably even increase when I use a higher smoothingFactor. Nov 2 '11 at 16:58
• @Pieter-JanBusschaert: Oh, I thought that the first peak was somehow useful to you. Anyway, all the low-pass-filter+move-down will have difficulties with the steep rise at ~650: They'll make this part flatter, and thus the curve has to be moved down a lot. Median Filter followed by a smoothing spline helps a bit. Nov 5 '11 at 2:53

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)