Median filter (Python)

From LiteratePrograms

This is a Python-implementation of the median image processing filter for 8-bit greyscale images. It uses the Python Imaging Library (PIL) for loading/displaying images and Psyco for performance improvements (but the latter is optional), which are not part of the standard Python distribution:

<<median.py>>=
import Image, array, math, time
import psyco   # use optimizer - comment out these two
psyco.full()   # lines to deactivate the optimizer

Psyco is an incredible piece of Python code: it compiles python bytecode into native-code on the fly, which leads to huge performance increases. However, it can lead to instabilities, especially in some Python IDEs written Python or in debuggers, so you can just turn it off by uncommenting these two lines. Also, I think it's only available for x86 platforms.

Now, to the algorithm: the median of a list is the "middle" element in the sorted version of the list. More generally, a quantile picks out a value at a specified position in the sorted version of a list (a 0%-quantile is the minimum, 50%-quantile is the median, and 100%-quantile is the maximum). So a straight-forward implementation could look like this:

<<median.py>>=
# simple implementation:
#def quantile(lst, f=0.5):
#    """Returns the f%-quantile of a list. Returns the median for f=0.5"""
#    return sorted(lst)[min(len(lst)-1, int(len(lst)*f))]

Unfortunately, this function is O=n*log(n), which can get pretty expensive for bigger filter masks. Also, sorting is an expensive operation in python (due to dynamic typing).

Luckily, we don't need the list fully sorted, we only need the nth-biggest value of it. Even more luckily, we know the range of the values in advance (image pixels are in the range 0-256). So, we can instead do a binary search on the byte range, looking for the value that has the specified number of values smaller than it in the list:

<<median.py>>=
# fast implementation:
def quantile(lst, f=0.5):
    """Returns the f%-quantile of a list. Returns the median for f=0.5. 
    Works only if all list members are integers in the range 0..255!"""
    l,r = 0,256
    targetCount = int(len(lst)*f)
    while l+1 < r:
        m = (l+r)/2
        count = 0
        for x in lst:
            if x < m: count += 1
        if count <= targetCount:
            l = m
        else:
            r = m
    return l

Note: There is a binary-search function in the Python standard library (bisect), but it would have to do lots of C-to-Python callbacks, which can't be optimized by Psyco, so I decided to implement it "by hand". Also, I didn't know about bisect at the time I wrote that code...

Here's another version using a histogram that's about 30% faster:

<<median.py>>=
# version using histogram
def quantile2(lst, f=0.5):
    # build histogram    
    h = [0]*256
    for x in lst:
        h[x] = h[x]+1
    # find median value in histogram
    sum = 0
    targetSum = int(len(lst)*f)
    for i in range(0, 255):
        sum += h[i]
        if (sum > targetSum):
            return i;

Now, the actual image filtering function: The first thing we have to do is get the pixels out of the image objects for faster access (Image.GetPixel is convenient, but slow, and worse, results in a Python-to-C-call which, you guessed it, can't be optimized by Psyco):

<<median.py>>=
def filter(img):
    """filters a region of an image with a circular-shaped median/quantile filter.
    Takes a PIL-Image as parameter and returns the filtered image."""
    # extract access variables from image
    src = array.array('B', img.getdata())
    dst = array.array('B', img.getdata())
    w,h = img.size

Considering that this is actually all it takes to get the "raw data" of an image, Python's PIL can be considered pretty image-processing friendly: Getting to the naked data of a bitmap image in C++, C# or Java is a lot more complex, which was actually the main reason why I chose Python when I did this: Trading performance for developer time seemed like a good idea for an experimental filter. And most of the optimizations made and tested in this code (like the quantile-function above or the offest-calculations below) could probably be translated 1:1 to C++ if needed.

Now, the next step is to prepare the mask, that is, the relative pixel positions that will be median'ed for every processed pixel:

<<median.py>>=
    # prepare circular filter mask
    r = 10
    mask = []
    for dy in range(-r, r+1):
        dx = int(math.sqrt(r*r-dy*dy))
        for x in range(-dx,dx+1):
            mask += [(x,dy)]

The next line converts these X/Y-offsets into raw-address-offsets in the image. This saves multiplications, and, more importantly, a tuple-unpacking in the image processing loop:

<<median.py>>=
    # calculare address offsets for mask coordinates
    maskOfs = [x+w*y for x,y in mask]

The image processing loop will go through each pixel, loop through the maskOfs list for them, and put the pixel values at (pixel address+mask offset) into this list:

<<median.py>>=
    pixelsInMask = [0]*len(mask)

and then calculate the median of that list.

And here it is:

<<median.py>>=
    starttime = time.clock()
    pixelsProcessed = 0
    # image processing loop
    for y in range(0+r,h-r):
        elapsed = time.clock()-starttime
        if elapsed > 0.1 and y % 10 == 0:
            print "Applying filter: %3i / %3i: %0.2f %%; %0.1f p/s\r" % (y, h, (y-r)*100.0/(h-r*2), pixelsProcessed/elapsed),
        for x in range(0+r,w-r):

As you can see, I don't process all of the pixels, but leave the borders unprocessed. This is because the mask offset would point "out of" the image for the border pixels, leading to invalid index-exceptions. I was only interested in the middle of the image anyway, so this didn't bother me. If you want to, you can insert tests for each pixel if x+mask.x s in the range 0..w, plus the same for y, and process the whole image (just testing if address+offset is in the array's bounds won't do!).

Now, the function gathers the pixel values for every pixel in the mask:

<<median.py>>=
            # filter kernel
            i = 0
            # gather all the pixels in the filter mask into to pixelsInMask list
            base = x+y*w
            for ofs in maskOfs:
                pixelsInMask[i] = src[base+ofs]
                i += 1

Calculates the median, and stores the result in the destination array

<<median.py>>=
            dst[x+y*w] = quantile(pixelsInMask, 0.5)
            pixelsProcessed += 1

That's it, more or less, only thing left to do is creating an image from the dst-array:

<<median.py>>=
    # create another image from results
    imgDest = img.copy()
    imgDest.putdata(dst)
    print
    return imgDest

And this is how to call it:

<<median.py>>=
if __name__=="__main__":
    img = Image.open("your-bitmap-here.bmp")
    imgDest = filter(img)
    imgDest.show()

Remember, it only works for 8-bit-greyscale images. Adjusting it for 24-bit RGB shouldn't be too hard though, if you understand what the code does and how the PIL packs the pixels.

Download code