分水岭算法——MATLAB

% These operations will create flat maxima inside each object that can be located using imregionalmax.
% 有多种方法可以应用在这里来获得前景标记，这些标记必须是前景对象内部的连接斑点像素。这个例子中，将使用形态学技术“基于开的重建”和“基于闭的重建”来清理图像。
% 接下来，通过腐蚀后重建来做基于开的重建计算。
Ie = imerode(I, se);
Iobr = imreconstruct(Ie, I);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(1, 2, 1); imshow(I, []); title('灰度图像');
% 使用MATLAB图像处理工具箱
% 注：期间用到了很多图像处理工具箱的函数，例如fspecial、imfilter、watershed、label2rgb、imopen、
% imclose、imreconstruct、imcomplement、imregionalmax、bwareaopen、graythresh和imimposemin函数等。
% The gradient is high at the borders of the objects and low (mostly) inside the objects.
% 使用Sobel边缘算子对图像进行水平和垂直方向的滤波，然后求取模值，sobel算子滤波后的图像在边界处会显示比较大的值，在没有边界处的值会很小。
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(1, 2, 1); imshow(gradmag,[]), title('梯度幅值图像')
subplot(1, 2, 2); imshow(Lrgb); 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.
% 开操作后，接着进行闭操作，可以移除较暗的斑点和枝干标记。对比常规的形态学闭操作和基于闭的重建操作。首先，使用imclose：
Ioc = imclose(Io, se);
Ic = imclose(I, se);
figure('units', 'normalized', 'position', [0 0 1 1]);
% No. Without additional preprocessing such as the marker computations below, using the watershed transform directly often results in "oversegmentation."
subplot(1, 2, 1); imshow(I, []); title('灰度图像');
subplot(1, 2, 2); imshow(Io), title('图像开操作')
% Next compute the opening-by-reconstruction using imerode and imreconstruct.
% 现在使用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]);
% 4.修改分割函数，使其仅在前景和后景标记位置有极小值。
% 5. Compute the watershed transform of the modified segmentation function.
% 5.对修改后的分割函数做分水岭变换计算。
% Use by Matlab Image Processing Toolbox
% Can you segment the image by using the watershed transform directly on the gradient magnitude?
% 可否直接对梯度幅值图像使用分水岭算法？
L = watershed(gradmag);
Lrgb = label2rgb(L);
% 开操作是腐蚀后膨胀，基于开的重建（基于重建的开操作）是腐蚀后进行形态学重建。下面比较这两种方式。首先，用imopen做开操作。
se = strel('disk', 20);
Io = imopen(I, se);
figure('units', 'normalized', 'position', [0 0 1 1]);
% 这些操作将会在每个对象内部创建单位极大值，使得可以使用imregionalmax来定位。
% 开运算和闭运算：先腐蚀后膨胀称为开；先膨胀后腐蚀称为闭。开和闭这两种运算可以除去比结构元素小的特定图像细节，同时保证不产生全局几何失真。
% 开运算可以把比结构元素小的突刺滤掉，切断细长搭接而起到分离作用；闭运算可以把比结构元素小的缺口或孔填充上，搭接短的间隔而起到连接作用。
% 直接使用梯度模值图像进行分水岭算法得到的结果往往会存在过度分割的现象。因此通常需要分别对前景对象和背景对象进行标记，以获得更好的分割效果。
% Step 3: Mark the Foreground Objects
% 第3步：标记前景对象
% A variety of procedures could be applied here to find the foreground markers,
% 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
% 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.
% 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.
hy = fspecial('sobel');
hx = பைடு நூலகம்y';
Iy = imfilter(double(I), hy, 'replicate');
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);
% 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.
figure('units', 'normalized', 'position', [0 0 1 1]);