Detecting Touching Objects Using Watershed Segmentation

Object detection in an image is an example of image segmentation. To segment touching objects, the watershed transform is often used. This transform finds intensity valleys in an image if you view an image as a surface with mountains (high intensity) and valleys (low intensity).

Key concepts:	Contrast enhancement, regional minima, watershed transform, label matrix, segmentation
Key functions	strel, imtophat, imbothat, imadd, imsubtract, imcomplement, imextendedmin, imimposemin, regionprops, label2rgb

Overview of Demo

The demo includes these steps:

• Step 1: Read Image
• Step 2: Maximize Contrast
• Step 3: Subtract Images
• Step 4: Convert Objects of Interest
• Step 5: Detect Intensity Valleys
• Step 6: Watershed Segmentation
• Step 7: Extract Features from Label Matrix

Step 1: Read Image

Read in the 'afmsurf.tif' image, which is an atomic force microscope image of a surface coating.

afm = imread('afmsurf.tif');
figure, imshow(afm), title('surface image');

Step 2: Maximize Contrast

We see many objects of different sizes that are touching each other. To minimize the number of valleys found by the watershed transform, we maximize the contrast of our objects of interest. A common technique for contrast enhancement is the combined use of the top-hat and bottom-hat transforms.

The top-hat transform is defined as the difference between the original image and its opening. The opening of an image is the collection of foreground parts of an image that fit a particular structuring element. The bottom-hat transform is defined as the difference between the closing of the original image and the original image. The closing of an image is the collection of background parts of an image that fit a particular structuring element.

The common structuring elements are squares, rectangles, disks, diamonds, balls, and lines. Since the objects of interest in our image look like disks, we create a disk structuring element having a radius of 15 pixels with the strel function. This disk size was an estimation of the average object radius in the image.

se = strel('disk', 15);

The imtophat and imbothat functions return the top-hat and bottom-hat transforms, respectively, of our original image.

Itop = imtophat(afm, se);
Ibot = imbothat(afm, se);
figure, imshow(Itop, []), title('top-hat image');


figure, imshow(Ibot, []), title('bottom-hat image');

Step 3: Subtract Images

We see that the top-hat image contains the "peaks" of objects that fit the structuring element. In contrast, the bottom-hat image shows the gaps between the objects of interest. To maximize the contrast between the objects and the gaps that separate them from each other, the "bottom-hat" image is subtracted from the "original + top-hat" image. The imadd function adds two uint8 images; similarly, the imsubtract function subtracts two uint8 images.

Ienhance = imsubtract(imadd(Itop, afm), Ibot);
figure, imshow(Ienhance), title('original + top-hat - bottom-hat');

Step 4: Convert Objects of Interest

Recall that watershed transform detects intensity "valleys" in an image. We use the imcomplement function on our enhanced image to convert our objects of interest to intensity valleys.

Iec = imcomplement(Ienhance);
figure, imshow(Iec), title('complement of enhanced image');

Step 5: Detect Intensity Valleys

We detect all intensity valleys below a particular threshold with the imextendedmin function. The output of the imextendedmin function is a binary image. The location rather than the size of the regions in the imextendedmin image is important. The imimposemin function will modify the image to contain only those valleys found by the imextendedmin function. The imimposemin function will also change a valley's pixel values to zero (deepest possible valley for uint8 images). All regions containing an imposed minima will be detected by the watershed transform.

Iemin = imextendedmin(Iec, 22);
Iimpose = imimposemin(Iec, Iemin);
figure, imshow(Iemin), title('extended minima image');


figure, imshow(Iimpose), title('imposed minima image');

Step 6: Watershed Segmentation

Watershed segmentation of the imposed minima image is accomplished with the watershed function.

wat = watershed(Iimpose);

The watershed function returns a label matrix containing non-negative numbers that correspond to watershed regions. Pixels that do not fall into any watershed region are given a pixel value of zero.

We view the label matrix as an image using label2rgb.

rgb = label2rgb(wat);
figure, imshow(rgb);
title('watershed segmented image');

Step 7: Extract Features from Label Matrix

Features can be extracted from the label matrix with the regionprops function. For example, we can calculate two measurements (area and orientation) and view them as a function of one another.

stats = regionprops(wat, 'Area', 'Orientation');
area = [stats(:).Area];
orient = [stats(:).Orientation];
figure, plot(area, orient, 'b*');
title('Relationship of Particle Orientation to Area');
xlabel('particle area (pixels)');
ylabel('particle orientation (degrees)');