I wanted to implement an extended depth of field algorithm last week for visualisation.

The idea of extended depth of focus is that you have a stack of images, taken at different focal points, and you build a single artificial image so that you get everything in focus. Some simple methods work OK, some very advanced methods work better. We were happy with OK for our setting. Here’s how to do it with mahotas:

Start with standard imports:

import numpy as np import mahotas as mh

We are going to assume that you have an `image` object, which has three dimensions: the stack, height, and width:

stack,h,w = image.shape

We use `mh.sobel` as the measure of infocusness for each pixel [1]:

focus = np.array([mh.sobel(t, just_filter=True) for t in image])

Now, we select the best slice at each pixel location:

best = np.argmax(focus, 0)

So far, very easy. The next part is the hard part. We want to do the following:

r = np.zeros((h,w))-1 for y in xrange(h): for x in xrange(w): r[y,x] = image[best[y,x], y, x]

But this is very slow (never run nested loops in Python if you can avoid it). We get the same result with a slightly less legible, but faster manipulation [2]:

image = image.reshape((stack,-1)) # image is now (stack, nr_pixels) image = image.transpose() # image is now (nr_pixels, stack) r = image[np.arange(len(image)), best.ravel()] # Select the right pixel at each location r = r.reshape((h,w)) # reshape to get final result

Et voilà!

§

Here is an example, from a stack of plankton images. This is is the maximum intensity projection:

This is the *most in-focus slice* (using the sobel operator):

And this is the extended depth of field result:

It is clearly sharper (look at some of the organelles on the top-left: they are a blur in the other two images),at the expense of some possible noise. I actually played around with blurring the image a little bit and it did improve things ever so slightly.

§

Cite the mahotas paper if using this in a scientific publication.

[1] | Other methods simply use a different measure here. |

[2] | I am not 100% convinced that this is the best. After all we create an array of size len(image) just to index. I would be happy to find an alternative. |

###### Related articles

- Mahotas software paper published (metarabbit.wordpress.com)