Fast Median Filter in C / C++ for `UINT16` 2D Array

Question

Does anyone know a fast median filter algorithm for 16-bit (unsigned short) arrays in c++?

http://nomis80.org/ctmf.html

This one seems quite promising, but it only seems to work with byte arrays. Does anyone know either how to modify it to work with shorts or an alternative algorithm?

Crossway

Did you try std::nth_element? It's O(n) compared to O(n log n) for a quicksort. – Municipality

You don't want to modify this algorithm to make it work with short since the running time per pixel is proportional to 2^n, where n is the number of bits in the datatype that is used. 256 for 8-bit arrays is already painful enough, you don't want to go to 65536 for 16-bit arrays. See my answer for a faster algorithm, even though it is O(log r) per pixel instead of O(1). – Nemesis

If you don't want to do median filtering, which is what you do in for example image processing where you find one median for each pixel, but just want to find one median, @smocking's comment is relevant. – Nemesis

Answer 1

The technique in the paper relies on creating a histogram with 256 bins for an 8 bit pixel channel. Converting to 16 bits per channel would require a histogram with 65536 bins, and a histogram is required for each column of the image. Inflating the memory requirements by 256 makes this a less efficient algorithm overall, but still probably doable with today's hardware.

Using their proposed optimization of breaking the histogram into coarse and fine sections should further reduce the runtime hit to only 16x.

For small radius values I think you'll find traditional methods of median filtering will be more performant.

Huntress

Answer 2

Fast Median Search - An ANSI C implementation (PDF) is something for C, it's a paper with the title "Fast median search: an ANSI C implementation". The author claims it's O(log(n)), he also provides some code, maybe it'll help you. It's not better than your suggested code, but maybe a look worth.

Theressa

The paper linked in the question is O(1) which is better than O(log n). – Huntress

Yep, but this one maybe gives some more impulses, but you are definitely right. I've made a small edit, clarifying my intention. – Theressa

@MarkRansom: O(1) doesn't automatically mean better than O(log(n)) as there is always a constant the O-notation hides. In this case, the O(1) algorithm it is much slower than the O(log(n)) algorithm (for other than 2-bit or maybe 4-bit channel values). The O(1) paper works with histograms, which makes the running time per pixel proportional to 2^b, where b is the number of bits per channel, which for 8 bits is 256 and for 16 bits is 65536 (even though these numbers are constants, hence O(1)). This quickly makes the O(1) algorithm extremely slow the more bits you add to the channel values. – Nemesis

Answer 3

std::vector<unsigned short> v{4, 2, 5, 1, 3};
std::vector<unsigned short> h(v.size()/2+1);
std::partial_sort_copy(v.begin(), v.end(), h.begin(), h.end());
int median = h.back();

Runs in O(N·log(N/2+1)) and does not modify your input.

Erubescence

Answer 4

This article describes a method for median filtering of images that runs in O(log r) time per pixel, where r is the filter radius, and works for any data type (be it 8 bit integers or doubles):

Fast Median and Bilateral Filtering

Nemesis

Is there a C / C++ implementation of Wiess' method anywhere? – Archduchy

@Archduchy I don't know. If you find out about any, please tell me! – Nemesis

Answer 5

See equations 4 and 5 in the following paper. The complexity is O(N*W) where W is the width of the filter and N is the number of samples.

See Noise Reduction by Vector Median Filtering.

Fatherly

Answer 6

I know this question is somewhat old but I also got interested in median filtering. If one is working with signals or images, then there will be a large overlap of data for the processing window. This can be taken advantage of.

I've posted some benchmark code here: 1D moving median filtering in C++

It's template based so it should work with most POD data types.

According to my results std::nth_element has poor performance for a moving median, as it must sort the window of values each time.

However, using a pool of values that is kept sorted, one can perform the median with 3 operation.

1. Remove oldest value out of the pool (calls std::lower_bound)
2. Insert new value into pool (calls std::lower_bound)
3. Store new value in history buffer

The median is now the middle value in the pool.

I hope someone finds this interesting and contributes their ideas!

Impropriety