1 綜述

Separating touching objects in an image is one of the more difficult image processing operations. 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.

1.計算分割函式。影象中較暗的區域是要分割的物件。

2. Compute foreground markers. These are connected blobs of pixels within each of the objects.

2.計算前景標誌。這些是每個物件內部連線的斑點畫素。

3. Compute background markers. These are pixels that are not part of any object.

3.計算背景標誌。這些是不屬於任何物件的畫素。

4. Modify the segmentation function so that it only has minima at the foreground and background marker locations.

4.修改分割函式，使其僅在前景和後景標記位置有極小值。

5. Compute the watershed transform of the modified segmentation function.

5.對修改後的分割函式做分水嶺變換計算。

Use by Matlab Image Processing Toolbox

使用MATLAB影象處理工具箱

注：期間用到了很多影象處理工具箱的函式，例如fspecial、imfilter、watershed、label2rgb、imopen、imclose、imreconstruct、imcomplement、imregionalmax、bwareaopen、graythresh和imimposemin函式等。

2 步驟

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

第一步：讀入彩色影象，將其轉化成灰度影象

clc; clear all; close all;

rgb = imread('pears.png');

if ndims(rgb) == 3

    I = rgb2gray(rgb);

else

    I = rgb;

end

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(1, 2, 1); imshow(rgb); title('原圖');

subplot(1, 2, 2); imshow(I); title('灰度圖');

Step 2: Use the Gradient Magnitude as the Segmentation Function

第2步：將梯度幅值作為分割函式

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.

使用Sobel邊緣運算元對影象進行水平和垂直方向的濾波，然後求取模值，sobel運算元濾波後的影象在邊界處會顯示比較大的值，在沒有邊界處的值會很小。

hy = fspecial('sobel');

hx = hy';

Iy = imfilter(double(I), hy, 'replicate');

Ix = imfilter(double(I), hx, 'replicate');

gradmag = sqrt(Ix.^2 + Iy.^2);

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(1, 2, 1); imshow(I,[]), title('灰度影象')

subplot(1, 2, 2); imshow(gradmag,[]), title('梯度幅值影象')

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

可否直接對梯度幅值影象使用分水嶺演算法？

L = watershed(gradmag);

Lrgb = label2rgb(L);

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(1, 2, 1); imshow(gradmag,[]), title('梯度幅值影象')

subplot(1, 2, 2); imshow(Lrgb); title('梯度幅值做分水嶺變換')

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

第3步：標記前景物件

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.

有多種方法可以應用在這裡來獲得前景標記，這些標記必須是前景物件內部的連線斑點畫素。這個例子中，將使用形態學技術“基於開的重建”和“基於閉的重建”來清理影象。這些操作將會在每個物件內部建立單位極大值，使得可以使用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.

開操作是腐蝕後膨脹，基於開的重建（基於重建的開操作）是腐蝕後進行形態學重建。下面比較這兩種方式。首先，用imopen做開操作。

se = strel('disk', 20);

Io = imopen(I, se);

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(1, 2, 1); imshow(I, []); title('灰度影象');

subplot(1, 2, 2); imshow(Io), title('影象開操作')

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

接下來，通過腐蝕後重建來做基於開的重建計算。

Ie = imerode(I, se);

Iobr = imreconstruct(Ie, I);

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(1, 2, 1); imshow(I, []); title('灰度影象');

subplot(1, 2, 2); imshow(Iobr, []), title('基於開的重建影象')

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:

開操作後，接著進行閉操作，可以移除較暗的斑點和枝幹標記。對比常規的形態學閉操作和基於閉的重建操作。首先，使用imclose：

Ioc = imclose(Io, se);

Ic = imclose(I, se);

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(2, 2, 1); imshow(I, []); title('灰度影象');

subplot(2, 2, 2); imshow(Io, []); title('開操作影象');

subplot(2, 2, 3); imshow(Ic, []); title('閉操作影象');

subplot(2, 2, 4); imshow(Ioc, []), title('開閉操作');

Now use imdilate followed by imreconstruct. Notice you must complement the image inputs and output of imreconstruct. IM2 = imcomplement(IM) computes the complement(補集) of the image IM. IM can be a binary, intensity, or RGB image. IM2 has the same class and size as IM.

現在使用imdilate，然後使用imreconstruct。注意必須對輸入影象求補，對imreconstruct輸出影象求補。IM2 = imcomplement(IM)計算影象IM的補集。IM可以是二值影象，或者RGB影象。IM2與IM有著相同的資料型別和大小。

Iobrd = imdilate(Iobr, se);

Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));

Iobrcbr = imcomplement(Iobrcbr);

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(2, 2, 1); imshow(I, []); title('灰度影象');

subplot(2, 2, 2); imshow(Ioc, []); title('開閉操作');

subplot(2, 2, 3); imshow(Iobr, []); title('基於開的重建影象');

subplot(2, 2, 4); imshow(Iobrcbr, []), title('基於閉的重建影象');

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.

通過比較Iobrcbr和loc可以看到，在移除小汙點同時不影響物件全域性形狀的應用下，基於重建的開閉操作要比標準的開閉重建更加有效。計算Iobrcbr的區域性極大來得到更好的前景標記。

fgm = imregionalmax(Iobrcbr);

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(1, 3, 1); imshow(I, []); title('灰度影象');

subplot(1, 3, 2); imshow(Iobrcbr, []); title('基於重建的開閉操作');

subplot(1, 3, 3); imshow(fgm, []); title('區域性極大影象');

To help interpret the result, superimpose(疊加) the foreground marker image on the original image.

為了幫助理解這個結果，疊加前景標記到原圖上。

It1 = rgb(:, :, 1);

It2 = rgb(:, :, 2);

It3 = rgb(:, :, 3);

It1(fgm) = 255; It2(fgm) = 0; It3(fgm) = 0;

I2 = cat(3, It1, It2, It3);

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(2, 2, 1); imshow(rgb, []); title('原影象');

subplot(2, 2, 2); imshow(Iobrcbr, []); title('基於重建的開閉操作');

subplot(2, 2, 3); imshow(fgm, []); title('區域性極大影象');

subplot(2, 2, 4); imshow(I2); title('區域性極大疊加到原影象');

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);

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(2, 2, 1); imshow(Iobrcbr, []); title('基於重建的開閉操作');

subplot(2, 2, 2); imshow(fgm, []); title('區域性極大影象');

subplot(2, 2, 3); imshow(fgm2, []); title('閉操作');

subplot(2, 2, 4); imshow(fgm3, []); title('腐蝕操作');

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. BW2 = bwareaopen(BW,P) removes from a binary image all connected components (objects) that have fewer than P pixels, producing another binary image, BW2.

這個過程將會留下一些偏離的孤立畫素，應該移除它們。可以使用bwareaopen，用來移除少於特定畫素個數的斑點。BW2 = bwareaopen(BW,P)從二值影象中移除所以少於P畫素值的連通塊，得到另外的二值影象BW2。

fgm4 = bwareaopen(fgm3, 20);

It1 = rgb(:, :, 1);

It2 = rgb(:, :, 2);

It3 = rgb(:, :, 3);

It1(fgm4) = 255; It2(fgm4) = 0; It3(fgm4) = 0;

I3 = cat(3, It1, It2, It3);

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(2, 2, 1); imshow(I2, []); title('區域性極大疊加到原影象');

subplot(2, 2, 2); imshow(fgm3, []); title('閉腐蝕操作');

subplot(2, 2, 3); imshow(fgm4, []); title('去除小斑點操作');

subplot(2, 2, 4); imshow(I3, []); title('修改區域性極大疊加到原影象');

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.

第4步：計算背景標記

現在，需要標記背景。在清理後的影象Iobrcbr中，暗畫素屬於背景，所以可以從閾值操作開始。

bw = im2bw(Iobrcbr, graythresh(Iobrcbr));

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(1, 2, 1); imshow(Iobrcbr, []); title('基於重建的開閉操作');

subplot(1, 2, 2); imshow(bw, []); title('閾值分割');

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) computes the Euclidean distance transform of the binary image BW. For each pixel in BW, the distance transform assigns a number that is the distance between that pixel and the nearest nonzero pixel of BW. bwdist uses the Euclidean distance metric by default. BW can have any dimension. D is the same size as BW.

背景畫素在黑色區域，但是理想情形下，不必要求背景標記太接近於要分割的物件邊緣。通過計算“骨架影響範圍”來“細化”背景，或者SKIZ，bw的前景。這個可以通過計算bw的距離變換的分水嶺變換來實現，然後尋找結果的分水嶺脊線（DL==0）。D = bwdist(BW)計算二值影象BW的歐幾里得矩陣。對BW的每一個畫素，距離變換指定畫素和最近的BW非零畫素的距離。bwdist預設使用歐幾里得距離公式。BW可以由任意維數，D與BW有同樣的大小。

D = bwdist(bw);

DL = watershed(D);

bgm = DL == 0;

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(2, 2, 1); imshow(Iobrcbr, []); title('基於重建的開閉操作');

subplot(2, 2, 2); imshow(bw, []); title('閾值分割');

subplot(2, 2, 3); imshow(label2rgb(DL), []); title('分水嶺變換示意圖');

subplot(2, 2, 4); imshow(bgm, []); title('分水嶺變換脊線圖');

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.

第5步：計算分割函式的分水嶺變換

函式imimposemin可以用來修改影象，使其只是在特定的要求位置有區域性極小。這裡可以使用imimposemin來修改梯度幅值影象，使其只在前景和後景標記畫素有區域性極小。

gradmag2 = imimposemin(gradmag, bgm | fgm4);

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(2, 2, 2); imshow(fgm4, []); title('前景標記');

subplot(2, 2, 3); imshow(gradmag, []); title('梯度幅值影象');

subplot(2, 2, 4); imshow(gradmag2, []); title('修改梯度幅值影象');

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

最後，可以做基於分水嶺的影象分割計算。

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.

第6步：檢視結果

一個視覺化技術是疊加前景標記、背景標記、分割物件邊界到初始影象。可以使用膨脹來實現某些要求，比如物件邊界，更加清晰可見。物件邊界定位於L==0的位置。

It1 = rgb(:, :, 1);

It2 = rgb(:, :, 2);

It3 = rgb(:, :, 3);

fgm5 = imdilate(L == 0, ones(3, 3)) | bgm | fgm4;

It1(fgm5) = 255; It2(fgm5) = 0; It3(fgm5) = 0;

I4 = cat(3, It1, It2, It3);

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(1, 2, 1); imshow(rgb, []); title('原影象');

subplot(1, 2, 2); imshow(I4, []); title('標記和物件邊緣疊加到原影象');

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.

另外一個有用的視覺化技術是將標記矩陣作為彩色影象進行顯示。標記矩陣，比如通過watershed和bwlabel得到的，可以使用label2rgb轉換到真彩圖像來顯示。

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(1, 2, 1); imshow(rgb, []); title('原影象');

subplot(1, 2, 2); imshow(Lrgb); title('彩色分水嶺標記矩陣');

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

可以使用透明度來疊加這個偽彩色標記矩陣在原亮度影象上進行顯示。

figure('units', 'normalized', 'position', [0 0 1 1]);

subplot(1, 2, 1); imshow(rgb, []); title('原影象');

subplot(1, 2, 2); imshow(rgb, []); hold on;

himage = imshow(Lrgb);

set(himage, 'AlphaData', 0.3);

title('標記矩陣疊加到原影象');

總結

程式碼：

clc; clear all; close all;

rgb = imread('pears.png');

if ndims(rgb) == 3

    I = rgb2gray(rgb);

else

    I = rgb;

end

hy = fspecial('sobel');

hx = hy';

Iy = imfilter(double(I), hy, 'replicate');

Ix = imfilter(double(I), hx, 'replicate');

gradmag = sqrt(Ix.^2 + Iy.^2);

L = watershed(gradmag);

Lrgb = label2rgb(L);

se = strel('disk', 20);

Io = imopen(I, se);

Ie = imerode(I, se);

Iobr = imreconstruct(Ie, I);

Ioc = imclose(Io, se);

Iobrd = imdilate(Iobr, se);

Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));

Iobrcbr = imcomplement(Iobrcbr);

fgm = imregionalmax(Iobrcbr);

It1 = rgb(:, :, 1);

It2 = rgb(:, :, 2);

It3 = rgb(:, :, 3);

It1(fgm) = 255; It2(fgm) = 0; It3(fgm) = 0;

I2 = cat(3, It1, It2, It3);

se2 = strel(ones(5,5));

fgm2 = imclose(fgm, se2);

fgm3 = imerode(fgm2, se2);

fgm4 = bwareaopen(fgm3, 20);

It1 = rgb(:, :, 1);

It2 = rgb(:, :, 2);

It3 = rgb(:, :, 3);

It1(fgm4) = 255; It2(fgm4) = 0; It3(fgm4) = 0;

I3 = cat(3, It1, It2, It3);

bw = im2bw(Iobrcbr, graythresh(Iobrcbr));

D = bwdist(bw);

DL = watershed(D);

bgm = DL == 0;

gradmag2 = imimposemin(gradmag, bgm | fgm4);

L = watershed(gradmag2);

It1 = rgb(:, :, 1);

It2 = rgb(:, :, 2);

It3 = rgb(:, :, 3);

fgm5 = imdilate(L == 0, ones(3, 3)) | bgm | fgm4;

It1(fgm5) = 255; It2(fgm5) = 0; It3(fgm5) = 0;

I4 = cat(3, It1, It2, It3);

Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');