版权声明：本文为博主原创文章，未经博主允许不得转载。					https://blog.csdn.net/u011583927/article/details/53413825				</div>
							            <link rel="stylesheet" href="https://csdnimg.cn/release/phoenix/template/css/ck_htmledit_views-d7e2a68c7c.css">
					<div class="htmledit_views" id="content_views">

--------------------------------------------------------------------------------------------------------------------

附录A  教程【3】给出的matlab源码，附详细注释

function [  ] = MarkerControlled_Watershed_tutorial(  )
%标记符控制的分水岭算法教程
%本程序按照官方文档的步骤，展示了如何使用分水岭算法，并测试对于含手势图片的分割结果是否理想
%程序修改自matlab官方文档 Marker-Controlled Watershed Segmentation
%网址： http://cn.mathworks.com/help/images/examples/marker-controlled-watershed-segmentation.html?s_tid=gn_loc_drop

%%
%官方文档简介

%这个例子展示了如何使用分水岭分割算法来分割一副图像中的物体
%分水岭变换经常应用到这类问题中。
%分水岭变换通过将灰度图像看成拓扑表面（亮的像素看做较高的点，暗像素看做低点），以此在一幅图中寻找“汇水盆地”和“分水岭脊线”

%如果你可以识别或者说标记mark前景物体与背景的位置（提供一些先验信息），使用分水岭变换的分割将会表现的很好。
%当然，人为提供这些标记“mark”，算法就算有效，很大程度上也失去了意义，于是就有了下面可以自行标记前景背景的分水岭算法。

%Marker-Controlled Watershed Segmentation，即标记符控制的分水岭分割算法，主要遵循以下几个基本步骤：
%1.计算一个分割函数（最基本的分水岭分割）。你尝试去分割的图像需要满足：黑暗区域对应着待分割的目标物体
%2.计算前景标记。得到的前景标记是每个物体内部的联通像素斑块。
%3.计算背景标记。背景标记由那些不属于任何物体的的像素点组成。
%4.修改分割函数，使得它只能在前景标记和背景标记的位置处达到最小值。
%5.计算修改后分割函数的分水岭变换

close all;

%%
%步骤一：读取彩色图像并转化为灰度图
HandImage = imread(‘testImgae/pear.jpg’);
I = rgb2gray(HandImage);
imshow(I)

%%
%步骤二：使用梯度幅值作为分割函数********
%使用Sobel边缘检测算子，imfilter函数，和一些简单的四则运算来计算梯度幅值。
%梯度值总是在物体的边缘处高，而总是在物体内部低。
hy = fspecial(‘sobel’);%获取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)’)

%你能直接在梯度幅值矩阵上应用分水岭变换来分割图像吗？
L = watershed(gradmag);
Lrgb = label2rgb(L);
figure, imshow(Lrgb), title(‘Watershed transform of gradient magnitude (Lrgb)’)
%从结果可以看到，显然不能
%如果没有额外的预处理，比如添加标记符（如下所述），直接使用分水岭变换通常导致严重的过分割

%取消下面3行代码的注释可以观察label2rgb函数的作用
% testA=[1 1 2 2 2 2 ;1 2 3 3 3 4;0 1 1 3 4 4];
% testArgb = label2rgb(testA);
% figure, imshow(testArgb), title(‘测试label2rgb函数的作用，使用简单的标记函数testA’)

%%
%步骤三：标记前景物体******************
%我们可以应用多种方法来获得前景标记，这些标记必须是前景对象内部的连通像素斑块。
%这个例子中，将使用形态学技术“基于重建的开操作”和“基于重建的闭操作”来清理图像。
%这些操作将会在每个对象内部创建平坦的极大值小斑块，这些斑块可以使用imregionalmax函数来定位。

%我们的目标是

%开运算和闭运算：先腐蚀后膨胀称为开；先膨胀后腐蚀称为闭（所谓腐蚀还是膨胀是针对白色而言来理解的）
%开和闭这两种运算可以除去比结构元素小的特定图像细节，同时保证不产生全局几何失真。
%开运算可以把比结构元素小的突刺滤掉，切断细长搭接而起到分离作用；闭运算可以把比结构元素小的缺口或孔填充上，搭接短的间隔而起到连接作用。

%但是

%开操作的第一步腐蚀能去除白色小物体，随后的膨胀趋向于恢复保留下来的物体的形状，这种恢复是不精确的，其精确度取决于形状和结构体之间的相似性。
%基于重建进行的开操作能够准确的恢复腐蚀之后的物体形状。
%所以基于重建的开操作与单纯的开操作功能类似，但更好地保留了原物体的形状。

%基于重建的开操作：I先腐蚀成Ie,后进行形态学重建。

%开操作是腐蚀后膨胀，基于重建的开操作是腐蚀后进行形态学重建。
%下面通过示例比较这两种方式，直观展示“基于重建的开操作”如何优于“开操作”
%首先，用imopen做开操作。
se = strel(‘disk’,20);%disk指定构建一个圆形的结构体，第二个参数指定结构体的半径
Io = imopen(I, se);%开操作
% figure
%imshow(Io), title(‘Opening (Io)’)

%接下来，进行基于重建的开操作。使用imerode和imreconstruct函数实现
Ie = imerode(I, se);%先腐蚀‘erosion’
Iobr = imreconstruct(Ie, I);%再重建
% figure
%imshow(Iobr), title(‘Opening-by-reconstruction (Iobr)’)

%对开操作后的结果进行闭操作，可以移除较暗的斑点和枝干标记。
%同样的，我们来对比常规的形态学闭操作和基于重建的闭操作的区别。

%首先，使用imclose对开操作的结果进行常规闭操作：
Ioc = imclose(Io, se);
% figure
%imshow(Ioc), title(‘Opening-closing 开操作+闭操作 (Ioc)’)

%然后我们来对开操作的结果进行基于重建的闭操作
%即首先使用imdilate函数进行膨胀，然后使用imreconstruct进行重建。
%这里要解释一下，重建操作是通过膨胀操作和交运算定义的，即 Hk+1=（Hk 膨胀 结构体B）∩G。
%重建定义的具体解释可以参考《数字图像处理的MATLAB实现（第2版）》P358页
%重建操作是对标记图像，在模板图像的限制下，进行膨胀。本质上是一种膨胀算法
%而这里需要的是模拟闭操作对膨胀后的图像进行腐蚀，怎么办呢？
%可以通过对模板图像和标记图像同时求补，然后进行重建操作，对补的膨胀相当于对原图的腐蚀。
%当然，重建的结果也要求补，才能得到实际的结果。
%IM2 = imcomplement(IM)计算图像IM的补集。IM可以是二值图像，或者RGB图像。IM2与IM有着相同的数据类型和大小。
Iobrd = imdilate(Iobr, se);%在基于重建的开操作的结果基础上，进行腐蚀
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));%重建，标记图像为腐蚀后图像取补，模板为腐蚀前原图取补。
Iobrcbr = imcomplement(Iobrcbr);%重建结果再取补，得到实际基于重建的闭操作的结果。
figure
imshow(Iobrcbr), title(‘Opening-closing by reconstruction 基于重建的开+闭操作 (Iobrcbr)’)
%如你所见，通过比较Iobrcbr和loc可以看到，在移除小污点同时不影响对象全局形状的应用下，基于开闭操作的重建要比标准的开闭操作更加有效。

%计算Iobrcbr的区域极大值来得到好的前景标记。
%得到的前景标记图fgm是二值图，白色对应前景区域
fgm = imregionalmax(Iobrcbr);
figure
imshow(fgm), title(‘Regional maxima of opening-closing by reconstruction (fgm)’)

%为了方便解释结果，将前景标记图叠加到原始图像上
I2 = I;
I2(fgm) = 255;%将fgm中的前景区域（像素值为1）标记到原图上（置白色）
figure
imshow(I2), title(‘Regional maxima superimposed on original image (I2)’)

%注意到一些大部分重合或被阴影遮挡的物体没有被标记出来。这意味着这些物体最终可能不会被正确的分割出来。
%并且，有些物体中前景标记正确的到达了物体的边缘。这意味着你应该清除掉标记斑块的边缘，向内收缩一点。
%你可以通过先闭操作，再腐蚀做到这点。
se2 = strel(ones(5,5));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);

%这个操作会导致遗留下一些离群的孤立点，这些是需要被移除的。
%你可以通过bwareaopen做到这点，函数将移除那些包含像素点个数少于指定值的区域。
fgm4 = bwareaopen(fgm3, 20);
I3 = I;
I3(fgm4) = 255;
figure
imshow(I3)
title(‘Modified regional maxima superimposed on original image (fgm4)’)

%%
%第四步：计算背景标记*************
%本例中设计出的标记背景算法的前提假设是：图像中相对亮的是物体，相对暗的区域是背景。
%如果不满足这条假设，标记结果可能不甚理想

%现在你需要标记背景。在去除噪点后的图像Iobrcbr中，暗像素属于背景，所以你可以先进行一下阈值操作。
%bw = imbinarize(Iobrcbr);
bw=im2bw(Iobrcbr,graythresh(Iobrcbr));%我使用上一行代码报错了，故换了一种二值化方法
figure
imshow(bw), title(‘Thresholded opening-closing by reconstruction (bw)’)

%背景像素现在是黑的了，但是理想情况下我们不希望背景标记太接近我们想要分割的物体的边界。
%我们要使背景变瘦，通过计算“骨架影响范围”来“细化”背景，或者SKIZ，bw的前景。这句没翻通，原文如下
%We’ll “thin” the background by computing the “skeleton by influence zones”, or SKIZ, of the foreground of bw.
%这个可以通过对bw的距离变换进行分水岭变换来实现，然后寻找结果的分水岭脊线（DL==0）。
D = bwdist(bw);
%D = bwdist(BW)计算二值图像BW的欧几里得矩阵。对BW的每一个像素，距离变换指定像素和最近的BW非零像素的距离。
%bwdist默认使用欧几里得距离公式。BW可以由任意维数，D与BW有同样的大小。

%由于bw中目标物体是白色的1.所以D中对应的目标物体处均是0，随着进入背景越深，对应像素值越大。
%这时正好符合我们使用分水岭算法的假设（想分出的目标物体数值较低）
%于是得到的背景标记是 物体与背景间的一个圈，能够包住目标物体
DL = watershed(D);
bgm = DL == 0;%分水岭变换结果L中，同一区域用同一数字表示，区域间分界线同一由0标识
figure
imshow(bgm), title(‘Watershed ridge lines (bgm)’)

%%
%*第五步：计算分割函数（修改后）的分水岭变换
% 函数imimposemin可以被用来修改一副图片，使得其只在指定的位置处取得局部最小值
% 这里你可以使用imimposemin来修改梯度幅值图像，使得局部最小值只出现在前景标记和背景标记处。
% 从结果来看imimposemin会将指定区域置为-Inf，从而成为极小值
% 且图像变得相当“平整”，一块块的相同数值的区域
gradmag2 = imimposemin(gradmag, bgm | fgm4);

%这段测试imimposemin都做了什么
% gradmag2;
% temp=gradmag2-gradmag;
%
% maxValue=max(max(temp));
% minValue=min(min(temp));
% temp=temp/maxValue;
% figure
% imshow(temp)

% maxValue=max(max(gradmag));
% temp_gradmag=gradmag/maxValue;
% figure
% imshow(temp_gradmag)
%
% maxValue=max(max(gradmag2));
% temp_gradmag2=gradmag2/maxValue;
% figure
% imshow(temp_gradmag2)

%Finally we are ready to compute the watershed-based segmentation.
L = watershed(gradmag2);

%%
%第六步：可视化结果******************
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)’)

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

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

--------------------------------------------------------------------------------------------------------------------

附录B  提炼后的源码，去掉多余演示步骤，只可视化关键步骤

function [  ] = MarkerControlled_Watershed(  )
%标记符控制的分水岭算法
%与教程函数不同，本函数没有多余的步骤和注释，直接实现
%展示的结果图更少更精炼

close all;

%%
%步骤一：读取彩色图像并转化为灰度图
HandImage = imread(‘testImgae/hand1.jpg’);
I = rgb2gray(HandImage);

%如果原图中，目标物体是较暗的，即与假设相反，这里可以取反
% Fan=ones(size(I,1),size(I,2))*255;
% Fan=uint8(Fan);
% Fan=Fan-I;
% I=Fan;

imshow(I), title(‘原图I’)
%%
%步骤二：基于重建的开闭操作******************

se = strel(‘disk’,20);%disk指定构建一个圆形的结构体，第二个参数指定结构体的半径

%接下来，进行基于重建的开操作。使用imerode和imreconstruct函数实现
Ie = imerode(I, se);%先腐蚀‘erosion’
Iobr = imreconstruct(Ie, I);%再重建

Iobrd = imdilate(Iobr, se);%在基于重建的开操作的结果基础上，进行腐蚀
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));%重建，标记图像为腐蚀后图像取补，模板为腐蚀前原图取补。
Iobrcbr = imcomplement(Iobrcbr);%重建结果再取补，得到实际基于重建的闭操作的结果。
figure
imshow(Iobrcbr), title(‘基于重建的开+闭操作 (Iobrcbr)’)

%%
%步骤三：标记前景物体******************
%计算Iobrcbr的区域极大值来得到好的前景标记。
%得到的前景标记图fgm是二值图，白色对应前景区域
fgm = imregionalmax(Iobrcbr);

%为了方便解释结果，将前景标记图叠加到原始图像上
I2 = I;
I2(fgm) = 255;%将fgm中的前景区域（像素值为1）标记到原图上（置白色）

%注意到一些大部分重合或被阴影遮挡的物体没有被标记出来。这意味着这些物体最终可能不会被正确的分割出来。
%并且，有些物体中前景标记正确的到达了物体的边缘。这意味着你应该清除掉标记斑块的边缘，向内收缩一点。
%你可以通过先闭操作，再腐蚀做到这点。
se2 = strel(ones(5,5));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);

%这个操作会导致遗留下一些离群的孤立点，这些是需要被移除的。
%你可以通过bwareaopen做到这点，函数将移除那些包含像素点个数少于指定值的区域。
fgm4 = bwareaopen(fgm3, 20);
I3 = I;
I3(fgm4) = 255;

%%
%第四步：计算背景标记*************
%本例中设计出的标记背景算法的前提假设是：图像中相对亮的是物体，相对暗的区域是背景。
%如果不满足这条假设，标记结果可能不甚理想

%现在你需要标记背景。在去除噪点后的图像Iobrcbr中，暗像素属于背景，所以你可以先进行一下阈值操作。
%bw = imbinarize(Iobrcbr);
bw=im2bw(Iobrcbr,graythresh(Iobrcbr));%我使用上一行代码报错了，故换了一种二值化方法

%背景像素现在是黑的了，但是理想情况下我们不希望背景标记太接近我们想要分割的物体的边界。

D = bwdist(bw);
%D = bwdist(BW)计算二值图像BW的欧几里得矩阵。对BW的每一个像素，距离变换指定像素和最近的BW非零像素的距离。
%bwdist默认使用欧几里得距离公式。BW可以由任意维数，D与BW有同样的大小。

%由于bw中目标物体是白色的1.所以D中对应的目标物体处均是0，随着进入背景越深，对应像素值越大。
%这时正好符合我们使用分水岭算法的假设（想分出的目标物体数值较低）
%于是得到的背景标记是 物体与背景间的一个圈，能够包住目标物体
DL = watershed(D);
bgm = DL == 0;%分水岭变换结果L中，同一区域用同一数字表示，区域间分界线同一由0标识

%%
%*第五步：计算分割函数（修改后）的分水岭变换
% 函数imimposemin可以被用来修改一副图片，使得其只在指定的位置处取得局部最小值
% 这里你可以使用imimposemin来修改梯度幅值图像，使得局部最小值只出现在前景标记和背景标记处。
% 从结果来看imimposemin会将指定区域置为-Inf，从而成为极小值
% 且图像变得相当“平整”，一块块的相同数值的区域

%使用Sobel边缘检测算子，imfilter函数，和一些简单的四则运算来计算梯度幅值。
%梯度值总是在物体的边缘处高，而总是在物体内部低。
hy = fspecial(‘sobel’);%获取sobel算子模板，计算纵向梯度
hx = hy’;%横向模板
Iy = imfilter(double(I), hy, ‘replicate’);%计算纵向梯度
Ix = imfilter(double(I), hx, ‘replicate’);
gradmag = sqrt(Ix.^2 + Iy.^2);%计算梯度幅值

gradmag2 = imimposemin(gradmag, bgm | fgm4);

%Finally we are ready to compute the watershed-based segmentation.
L = watershed(gradmag2);

%%
%第六步：可视化结果******************
I4 = I;
I4(imdilate(L == 0, ones(3, 3)) | bgm | fgm4) = 255;
figure
imshow(I4)
title(‘在原图上绘制前景、背景标记，以及分割边界’)

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

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

end

• 算法介绍

原始的分水岭算法对于存在噪声以及梯度不规则的图像极易造成过度分割（存在大量分割区域），解决该问题的一种方案是通过融入预处理步骤来限制允许存在的区域数目，因此基于标记的分水岭算法应运而生。标记是属于一幅图像的一个连通分量，与感兴趣物体相联系的标记称为内部标记，与背景相关联的标记称为外部标记。

• Opencv函数介绍

void watershed( 
                 InputArray image,         //Input 8-bit 3-channel image
                 InputOutputArray markers  //Input/output 32-bit single-channel 
              );

第一个参数image是CV_8UC3格式的原始图像；

第二个参数markers是CV_32S格式的标记图像，即该图像必须包含注水点（即种子点）。注水点的获取有以下三种方式：

        第一种：使用findContours函数自动标记。

        第二种：connectedComponents函数标记连通域。

        第三种：用鼠标划线标记（opencv官方例程即使用该方法），该方法需要手动操作。

算法执行后，区域间的边界像素会设为-1，所有0值像素会被设置为给定标记中的一个。

• Opencv实例

• 效果图

                                                                                                   原图

                                                                                                  结果图

• 代码

//
//opencv4.1.0
//

#include <opencv2/opencv.hpp>
#include <iostream>

using namespace std;
using namespace cv;

void findContours_watershed(Mat& dist_binary, Mat& src);        //第一种方式使用findContours函数标记
void connectedComponents_watershed(Mat& dist_binary, Mat& src); //第二种方式：使用connectedComponents函数标记
void DisplayRegion(int num, Mat& src, Mat& markers, const cv::String& winname);//使用不同的颜色显示分割后的各个区域,便于观察结果

int main() {
	Mat src, gray, binary, open_img, dilate_img, dist_img, dist_binary;
	src = imread("12.png");
	cvtColor(src, gray, COLOR_BGR2GRAY);
	threshold(gray, binary, 0, 255, THRESH_BINARY_INV | THRESH_OTSU);

	//形态学操作
	Mat k = getStructuringElement(MORPH_RECT, Size(3, 3));
	morphologyEx(binary, open_img, MORPH_OPEN, k, Point(-1,-1), 2);
	dilate(open_img, dilate_img, k, Point(-1, -1), 3);

	//距离变换
	distanceTransform(open_img, dist_img, DIST_L2, 5);

	//归一化方法一
	//normalize(dist_img, dist_img, 0, 1,NORM_MINMAX);
	
	//归一化方法二
	double minVal,maxVal;
	minMaxLoc(dist_img, &minVal, &maxVal,0,0);
	Mat dst;
	dist_img.convertTo(dst, CV_8UC1, 1 / (maxVal-minVal)*255);

	//阈值分割提取每个区域的标记图像
	threshold(dst, dist_binary, 180, 255, 0);	

	//方式一
	findContours_watershed(dist_binary, src);
	//方式二
	connectedComponents_watershed(dist_binary, src);

	imshow("原图", src);
	waitKey(0);
	return 0;
}

void connectedComponents_watershed(Mat& dist_binary, Mat& src) {

	//定义一个标记图像
	Mat markers(src.size(), CV_32S);
	
	//计算连通域
	int num = connectedComponents(dist_binary, markers);

	//基于标记的分水岭分割
	watershed(src, markers);

	//使用不同的颜色显示各个区域
	DisplayRegion(num, src, markers, "result2");
}

void findContours_watershed(Mat& dist_binary, Mat& src) {

	//定义一个标记图像，并初始化为0
	Mat markers(src.size(), CV_32S);
	markers = Scalar::all(0);

	//寻找轮廓
	vector<vector<Point>> contours;
	vector<Vec4i> hierarchy;
	findContours(dist_binary, contours, hierarchy, RETR_CCOMP, CHAIN_APPROX_SIMPLE);

	//绘制轮廓
	for (int index = 0; index < contours.size(); index++)
		drawContours(markers, contours, index, Scalar::all(index), -1, 8, hierarchy, INT_MAX);

	//基于标记的分水岭分割
	watershed(src, markers);

	//使用不同的颜色显示各个区域
	int num = contours.size();
	DisplayRegion(num, src, markers, "result1");
}

void DisplayRegion(int num,Mat& src, Mat& markers,const cv::String& winname) {
	RNG rng(123456);
	vector<Vec3b> color(num);
	//color[0] = Vec3b(0, 0, 0);                                                           //设置背景颜色
	for (size_t i = 1; i < num; i++) {
		color[i] = Vec3b(rng.uniform(0, 256), rng.uniform(0, 256), rng.uniform(0, 256)); //设置目标颜色
	}

	Mat markersColor = Mat::zeros(src.size(), src.type());
	int w = src.cols;
	int h = src.rows;
	for (size_t row = 0; row < h; row++) {
		for (size_t col = 0; col < w; col++) {
			int label = markers.at<int>(row, col);
			if (label == -1) continue;
			markersColor.at<Vec3b>(row, col) = color[label];
		}
	}
	imshow(winname, markersColor);
}

1、和分水岭算法相关的资料

[1] 标记符控制的分水岭算法原理及matlab实现_CSDN
[2] 图像处理——分水岭算法
[3]

2、《数字图像处理：原理与实践》

基于标记控制的分水岭分割方法有以下基本步骤：
1.计算分割函数。图像中较暗的区域是要分割的对象。
2.计算前景标志。这些是每个对象内部连接的斑点像素。
3.计算背景标志。这些是不属于任何对象的像素。
4.修改分割函数，使其仅在前景和后景标记位置有极小值。
5.对修改后的分割函数做分水岭变换计算。

（1）读入彩色图，转为为灰度图

rgb = imread('tudou500.jpg');
I = rgb2gray(rgb);

（2）使用Sobel算子边缘滤波

hy = fspecial('sobel');  % 水平 y 轴
hx = hy';  % 竖直 x
Iy = imfilter(double(I), hy, 'replicate');
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);
figure;
subplot(221);imshow(rgb);title('original');
subplot(222);imshow(gradmag, []);title('梯度幅值图像');
subplot(223);imshow(Iy, []);title('水平边缘检测');
subplot(224);imshow(Ix, []);title('垂直边缘检测');

（2-1）直接对梯度幅值图像gradmag进行watershed()操作

L = watershed(gradmag);
Lrgb = label2rgb(L);
figure;
subplot(211);imshow(I);title('original');
subplot(212);imshow(Lrgb);title('对梯度幅值直接做分水岭分割');

结果，过度分割达不到我们想要的效果。

（3）对图像进行增强

se = strel('disk',13);
Ie = imerode((I),se);
% Ie = imerode(imcomplement(I),se);  % 腐蚀
% Ie = imcomplement(Ie);
Iobr = imreconstruct(Ie,I);  % 重建
figure
subplot(231);imshow(I);title('original I');
subplot(232);imshow(Ie);title('腐蚀 Ie');
subplot(233);imshow(Iobr);title('基于开的重建 Iobr');

Iobrd = imdilate(Iobr, se);
Iobrdr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));
Iobrdr = imcomplement( Iobrdr );
subplot(234);imshow(I);title('original I');
subplot(235);imshow(Iobrd);title('Iobr + dilate = Iobrd');
subplot(236);imshow(Iobrdr);title('基于重建的开闭操作 Iobrd+r=Iobrdr');

（4）获取前景标记

fgm = imregionalmax(Iobrdr);
It1 = rgb(:,:,1);
It2 = rgb(:,:,2);
It3 = rgb(:,:,3);
% 将fgm区域变成红色，其他区域不变
It1(fgm) = 255; It2(fgm) = 0;  It3(fgm) = 0;  
I2 = cat(3,It1, It2, It3);
figure
subplot(121);imshow(fgm,[]);title('局部极大值图像');
subplot(122);imshow(I2,[]);title('局部极大值叠加图像');

（5）极大值结果修正

se2 = strel(ones(12,12));
fgm2 = imcomplement(imclose(imcomplement(fgm), se2));
% fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);
% fgm3 = imerode(fgm3, se2);
fgm3 = imdilate(fgm3, se2);
fgm3 = imdilate(fgm3, se2);
fgm3 = imerode(fgm3, se2);
% fgm3 = imopen(fgm2, se2);
fgm4 = bwareaopen(fgm3, 300);
figure
subplot(221);imshow(fgm);title('fgm');
subplot(222);imshow(fgm2);title('fgm2');
subplot(223);imshow(fgm3);title('fgm3');
subplot(224);imshow(fgm4);title('fgm4');

（6）使用分水岭分割watershed()

BW = im2bw(Iobrdr, graythresh(Iobrdr));
D = bwdist(BW);
DL = watershed(D);
bgm = DL == 0;  % DL == 0 确定矩阵DL等于0的值所在位置，=0的位置置为1
gradmag2 = imimposemin(gradmag, bgm | fgm4);


%% 分水岭分割
L = watershed(gradmag2);
% first way show
Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');
figure
subplot(121);imshow(Lrgb);title('分水岭结果显示');

% second way show
subplot(122);imshow(rgb, []);title('分水岭结果显示（叠加）');
hold on;
himage = imshow(Lrgb);
set(himage, 'AlphaData', 0.3)

（上周就做完的实验，但是换方向，拖到了周六/(ㄒoㄒ)/~~）
end