# Find the stomata in a plant microscopy image

Here is a question for image processing experts.

I am working on a difficult computer vision problem. The task is to count the stomata (marked below) in DIC microscopy images. These images are resistant to most superficial image processing techniques like morphological operations and edge detection. It is also different from other cell counting tasks.

I am using OpenCV. My plan is to review potentially useful features for stomata discrimination.

• Texture classifiers
• DCT (Discrete cosine transform/frequency-domain analysis)
• LBP (Local binary patterns)
• HOG (Histogram of oriented gradients)
• Robust feature detectors (I am skeptical)
• Harris corners
• SIFT, SURF, STAR, etc.

And possibly design a novel feature descriptor. I am leaving out the selection of a classifier for now.

What have I missed? How would you solve this? Solutions for similar object detection problems would be very helpful.

Sample images here.

After bandpass filter:

Canny edge detection is not promising. Some image areas are out of focus:

• Maybe instead of trying to find the stomata, you could try to remove the mazy lines? – endolith Nov 22 '11 at 1:26
• How many images do you have to process? How fast does it need to be? How automated does it have to be? – endolith Nov 22 '11 at 20:36
• It does not have to be very fast. We are processing on the order of 1000 images. It should be automatic - dump images into a directory and go. – Matt M. Nov 23 '11 at 3:03

Sorry I don't know OpenCV, and this is more a pre-processing step than a complete answer:

First, you don't want an edge detector. An edge detector converts transitions (like this dark-to-light):

into ridges (bright lines on dark) like this:

It performs a differentiation, in other words.

But in your images, there is a light shining down from one direction, which shows us the relief of the 3D surface. We perceive this as lines and edges, because we're used to seeing things in 3D, but they aren't really, which is why edge detectors aren't working, and template matching won't work easily with rotated images (a perfect match at 0 degrees rotation would actually cancel out completely at 180 degrees, because light and dark would line up with each other).

If the height of one of these mazy lines looks like this from the side:

then the brightness function when illuminated from one side will look like this:

This is what you see in your images. The facing surface becomes brighter and the trailing surface becomes darker. So you don't want to differentiate. You need to integrate the image along the direction of illumination, and it will give you the original height map of the surface (approximately). Then it will be easier to match things, whether through Hough transform or template matching or whatever.

I'm not sure how to automate finding the direction of illumination. If it's the same for all your images, great. Otherwise you'd have to find the biggest contrast line and assume the light is perpendicular to it or something. For my example, I rotated the image manually to what I thought was the right direction, with light coming from the left:

You also need to remove all the low-frequency changes in the image, though, to highlight only the quickly-changing line-like features. To avoid ringing artifacts, I used 2D Gaussian blur and then subtracted that from the original:

The integration (cumulative sum) can runaway easily, which produces horizontal streaks. I removed these with another Gaussian high-pass, but only in the horizontal direction this time:

Now the stomata are white ellipses all the way around, instead of white in some places and black in others.

Original:

Integrated:

from pylab import *
import Image
from scipy.ndimage import gaussian_filter, gaussian_filter1d

filename = 'rotated_sample.jpg'
I = Image.open(filename).convert('L')
I = asarray(I)

# Remove DC offset
I = I - average(I)

close('all')
figure()
imshow(I)
gray()
show()
title('Original')

# Remove slowly-varying features
sigma_2d = 2
I = I - gaussian_filter(I, sigma_2d)

figure()
imshow(I)
title('2D filtered with %s' % sigma_2d)

# Integrate
summed = cumsum(I, 1)

# Remove slowly-changing streaks in horizontal direction
sigma_1d = 5
output = summed - gaussian_filter1d(summed, sigma_1d, axis=1)

figure()
imshow(output)
title('1D filtered with %s' % sigma_1d)


The Hough transform can be used to detect ridge ellipses like this, made of "edge pixels", though it's really expensive in computation and memory, and they are not perfect ellipses so it would have to be a bit of a "sloppy" detector. I've never done it, but there are a lot of Google results for "hough ellipse detection". I'd say if you detect one ellipse inside the other, within a certain size search space, it should be counted as a stoma.

Also see:

• P.S. Does what I did here have a name? Is it a common filter type? – endolith Apr 16 '12 at 21:25
• +1 - Great answer! About the automation of light source angle - you could use edge detector that computes both magnitude and gradient and then compute the weighted (by mag.) average of the gradient. The strongest responses should be in the direction of the illumination. – Andrey Rubshtein Oct 5 '12 at 18:33

First thing I would try is template matching, with templates rotated for all the angles with some step. Rotating template essential here. Also choice of template could be non-trivial - could be several with different lighting, and it could be blurred to allow for difference in shapes.

http://en.wikipedia.org/wiki/Template_matching#Template-based_matching_and_convolution

Next - HOG looks promising here. Another solution could be using strong corner detector like Moravec or Shi-Tomasi (with non-maximum suppression) and look for groups of 2-corners or 3-4 corners on the same line as candidates. After finding candidates you can apply active contour for verification (not sure if it would really help, but that is possibility)

http://en.wikipedia.org/wiki/Corner_detector

http://en.wikipedia.org/wiki/Active_contour

Yet another possibility is to use Hough transform for ellipses, possibly with not 2 but 3-4 free parameters.

Partial answer. Finding candidates with Mathematica:

p = ColorConvert[Import@"http://i.stack.imgur.com/38Ysw.jpg",