Example: Automatic color image segmentation.

This example illustrates automatic segmentation of color images using several of the new ImageProcessing functions added in Version 2. The scheme employed in this example partitions the image by computing watershed zones in a surface formed from a computation of the second-order (joint) histogram of its chromatic components.

$HistoryLength = 1;

This loads the package and an example color image image.

<< ImageProcessing`

Show[Graphics[img = ImageRead["peppers.tif"]]];

Lab domain (peppers.tif)

It is well-known that the CIE \(TraditionalForm\`\(L\^*\) \(a\^*\) \(b\^*\)\) color is close to optimal in separating the color information in an image from its brigthness. The color  channels \(TraditionalForm\`\(a\^*\)\) and \(TraditionalForm\`\(b\^*\)\) roughly represent the color differences red minus green and green minus blue, respectively. The following transforms the example image to the  \(TraditionalForm\`\(L\^*\) \(a\^*\) \(b\^*\)\) color space and separates the channels.

\(\((imgLab = ToCIELABColor[img/255. ]; 
{L, a, b} = 
        Round[ToChannels[imgLab]];)\)\)

This calculates and displays the joint histogram of the two color channels.

\(\((h = JointImageHistogram[a, b, ImageRange[a], ImageRange[b]]; 

      ListPlot3D[h[\([\)\(All, All, 3\)\(]\)], Mesh -> False, 
        PlotRange -> All, 
        MeshRange -> {ImageRange[a], ImageRange[b]}];)\)\)

Here, the histogram data is rescaled and inverted to form a topographic surface suitable for partitioning by the watershed transformation.

data = 
        1 - ScaleLinear[
            h[[All, All, 3]], {0, 
              1}];

However, the surface exhibits many small local minima which leads to oversegmentation of the domain. This problem has been widely studied with numerous suggested solutions to choose from. As our interest lies in selecting only the most significant peaks in the histogram data we use recursive smoothing with the nonlinear median filter to obtain the so-called root median of the surface.

tmp = FixedPoint[BlockMedian[#, {3, 3}] &, data];

This shows the original and smoothed surfaces representing the joint color histogram of the example image.

DisplayTogetherArray[{ListPlot3D[data, Mesh -> False, 
            PlotRange -> All], 
          ListPlot3D[tmp, Mesh -> False, 
            PlotRange -> All]}];

The following shows the watershed segmentation results for successive iterates of the root median calculation.

tmp = data;
Do[(tmp = 
          BlockMedian[tmp, {3, 3}]; 
ListDensityPlot[
          WatershedTransform[tmp], Mesh -> False, AspectRatio -> 1, 
          ImageSize -> 120]), {12}]

The partitioning of the root median surface results in two distinct regions in image's chromatic space.

label = WatershedTransform[tmp];

Next, the locations of all the pixels belonging to each of the two regions must be calculated. First, the following returns all the color pairs belonging to each of the partitions.

minvals = {Min[a], Min[b]};
Do[
        region[i] = minvals + # & /@ Where[label, i], {i, 
          2}];

In the final step, the pixel locations for each of the regions are found and the corresponding image masks are formed.

{h, w} = ImageDimensions[img];
Do[
        tmp = Flatten[Intersection[Where[a, #[[1]]], 
                  Where[b, #[[2]]]] & /@ region[i], 
            1]; 
mask[i] = 
          1 - ToGrayLevel[FromSet[{h, w}, tmp, {0, 0}]], {i, 
          2}];

This shows the two masks.

DisplayTogetherArray[{Graphics[mask[1]], Graphics[mask[2]]}, 
        ImageSize -> 400];

Here are the final segmentation results shown in the original RGB color space.

DisplayTogetherArray[{Graphics[img ToRGBColor[mask[1]]], 
          Graphics[img ToRGBColor[mask[2]]]}, 
        ImageSize -> 400];