www.gusucode.com > images 案例代码 matlab源码程序 > images/WatershedSegmentationExample.m

    %% Marker-Controlled Watershed Segmentation
% This example shows how to use watershed segmentation to separate touching
% objects in an image. The watershed transform is often applied to this
% problem.  The watershed transform finds "catchment basins" and "watershed
% ridge lines" in an image by treating it as a surface where light pixels are
% high and dark pixels are low.
%
% Segmentation using the watershed transform works better if you can
% identify, or "mark," foreground objects and background locations.
% Marker-controlled watershed segmentation follows this basic procedure:
%
% 1. Compute a segmentation function.  This is an image whose dark
% regions are the objects you are trying to segment.
%
% 2. Compute foreground markers.  These are connected blobs of pixels
% within each of the objects.
%
% 3. Compute background markers.  These are pixels that are not part of
% any object.
%
% 4. Modify the segmentation function so that it only has minima at the
% foreground and background marker locations.
% 
% 5. Compute the watershed transform of the modified segmentation function.
%
% This example highlights many different Image Processing Toolbox(TM)
% functions, including |fspecial|, |imfilter|, |watershed|, |label2rgb|,
% |imopen|, |imclose|, |imreconstruct|, |imcomplement|, |imregionalmax|,
% |bwareaopen|, |graythresh|, and |imimposemin|.

% Copyright 1993-2015 The MathWorks, Inc.

%% Step 1: Read in the Color Image and Convert it to Grayscale

rgb = imread('pears.png');
I = rgb2gray(rgb);
imshow(I)

text(732,501,'Image courtesy of Corel(R)',...
     'FontSize',7,'HorizontalAlignment','right')

%% Step 2: Use the Gradient Magnitude as the Segmentation Function
% Use the Sobel edge masks, |imfilter|, and some simple arithmetic to
% compute the gradient magnitude.  The gradient is high at the borders of
% the objects and low (mostly) inside the objects.

hy = fspecial('sobel');
hx = hy';
Iy = imfilter(double(I), hy, 'replicate');
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);
figure
imshow(gradmag,[]), title('Gradient magnitude (gradmag)')

%% 
% Can you segment the image by using the watershed transform directly on
% the gradient magnitude?

L = watershed(gradmag);
Lrgb = label2rgb(L);
figure, imshow(Lrgb), title('Watershed transform of gradient magnitude (Lrgb)')

%%
% No.  Without additional preprocessing such as the marker computations below,
% using the watershed transform directly often results in
% "oversegmentation."

%% Step 3: Mark the Foreground Objects
% A variety of procedures could be applied here to find the foreground
% markers, which must be connected blobs of pixels inside each of the
% foreground objects.  In this example you'll use morphological
% techniques called "opening-by-reconstruction" and
% "closing-by-reconstruction" to "clean" up the image.  These
% operations will create flat maxima inside each object that can be
% located using |imregionalmax|.

%%
% Opening is an erosion followed by a dilation, while
% opening-by-reconstruction is an erosion followed by a morphological
% reconstruction.  Let's compare the two.  First, compute the opening using
% |imopen|.

se = strel('disk', 20);
Io = imopen(I, se);
figure
imshow(Io), title('Opening (Io)')

%%
% Next compute the opening-by-reconstruction using |imerode| and
% |imreconstruct|.

Ie = imerode(I, se);
Iobr = imreconstruct(Ie, I);
figure
imshow(Iobr), title('Opening-by-reconstruction (Iobr)')

%%
% Following the opening with a closing can remove the dark spots
% and stem marks.  Compare a regular morphological closing with a
% closing-by-reconstruction.  First try |imclose|:

Ioc = imclose(Io, se);
figure
imshow(Ioc), title('Opening-closing (Ioc)')

%%
% Now use |imdilate| followed by |imreconstruct|.  Notice you must 
% complement the image inputs and output of |imreconstruct|.

Iobrd = imdilate(Iobr, se);
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));
Iobrcbr = imcomplement(Iobrcbr);
figure
imshow(Iobrcbr), title('Opening-closing by reconstruction (Iobrcbr)')

%%
% As you can see by comparing |Iobrcbr| with |Ioc|, reconstruction-based
% opening and closing are more 
% effective than standard opening and closing at removing small blemishes
% without affecting the overall shapes of the objects.  Calculate the
% regional maxima of |Iobrcbr| to obtain good foreground markers.

fgm = imregionalmax(Iobrcbr);
figure
imshow(fgm), title('Regional maxima of opening-closing by reconstruction (fgm)')

%%
% To help interpret the result, superimpose the foreground marker image
% on the original image. 

I2 = I;
I2(fgm) = 255;
figure
imshow(I2), title('Regional maxima superimposed on original image (I2)')

%%
% Notice that some of the mostly-occluded and shadowed objects are not
% marked, which means that these objects will not be segmented properly
% in the end result.  Also, the foreground markers in some objects go right
% up to the objects' edge.  That means you should clean the edges of the
% marker blobs and then shrink them a bit.  You can do this by a closing
% followed by an erosion.

se2 = strel(ones(5,5));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);

%%
% This procedure tends to leave some stray isolated pixels that must be
% removed.  You can do this using |bwareaopen|, which removes all blobs
% that have fewer than a certain number of pixels. 

fgm4 = bwareaopen(fgm3, 20);
I3 = I;
I3(fgm4) = 255;
figure
imshow(I3)
title('Modified regional maxima superimposed on original image (fgm4)')

%% Step 4: Compute Background Markers
% Now you need to mark the background.  In the cleaned-up image, |Iobrcbr|,
% the dark pixels belong to the background, so you could start with a
% thresholding operation.

bw = imbinarize(Iobrcbr);
figure
imshow(bw), title('Thresholded opening-closing by reconstruction (bw)')

%%
% The background pixels are in black, but ideally we don't want the
% background markers to be too close to the edges of the objects we are
% trying to segment.  We'll "thin" the background by computing the
% "skeleton by influence zones", or SKIZ, of the foreground of |bw|.
% This can be done by computing the watershed transform of the distance
% transform of |bw|, and then looking for the watershed ridge lines (|DL ==
% 0|) of the result.

D = bwdist(bw);
DL = watershed(D);
bgm = DL == 0;
figure
imshow(bgm), title('Watershed ridge lines (bgm)')

%% Step 5: Compute the Watershed Transform of the Segmentation Function.
% The function |imimposemin| can be used to modify an image so that it has 
% regional minima only in certain desired locations.  Here you can use
% |imimposemin| to modify the gradient magnitude image so that its only
% regional minima occur at foreground and background marker pixels.

gradmag2 = imimposemin(gradmag, bgm | fgm4);

%%
% Finally we are ready to compute the watershed-based segmentation.

L = watershed(gradmag2);

%% Step 6: Visualize the Result
% One visualization technique is to superimpose the foreground
% markers, background markers, and segmented object boundaries on the
% original image.  You can use dilation as needed to make certain aspects,
% such as the object boundaries, more visible.  Object boundaries are
% located where |L == 0|.

I4 = I;
I4(imdilate(L == 0, ones(3, 3)) | bgm | fgm4) = 255;
figure
imshow(I4)
title('Markers and object boundaries superimposed on original image (I4)')

%%
% This visualization illustrates how the locations of the foreground and
% background markers affect the result.  In a couple of locations,
% partially occluded darker objects were merged with their brighter
% neighbor objects because the occluded objects did not have foreground
% markers.
%
% Another useful visualization technique is to display the label matrix
% as a color image.  Label matrices, such as those produced by
% |watershed| and |bwlabel|, can be converted to truecolor images for
% visualization purposes by using |label2rgb|.

Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');
figure
imshow(Lrgb)
title('Colored watershed label matrix (Lrgb)')

%%
% You can use transparency to superimpose this pseudo-color label matrix on
% top of the original intensity image.

figure
imshow(I)
hold on
himage = imshow(Lrgb);
himage.AlphaData = 0.3;
title('Lrgb superimposed transparently on original image')