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.
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, 1); imshow(bgm, []); title('分水嶺變換脊線圖');
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轉換到真彩圖像來顯示。
Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');
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('標記矩陣疊加到原圖像');
參考：http://blueben-chong.blog.sohu.com/141881444.html
總結
代碼：
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');