# Tutorial: Extended Depth of Field¶

This is an example of how to use mahotas to implement an algorithm that it does
not have built-in: *extended depth of field*.

The idea is that you have a stack of images, taken at different focal points, and you build a single image so that you get everything in focus.

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 microbes imaged. This is is the maximum intensity projection:

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

And this is the extended depth of field result:

It is clearly sharper, perhaps 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.

[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. |