%读取原RGB图象并将其转化为灰度图象
I = imread('lena.bmp');
figure
imshow(I)
I = imfill(I,'holes');

figure
imshow(I)
title('Original gray image(I)');%原始灰度图像

%%计算灰度梯度（Gradient Magnitude）
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)');%梯度幅度

%%标记前景目标（内部约束）
%基于重建（reconstruct）的开—闭运算相对基本的开—闭运算来说，在不影响目标图象的前提下，去除图象中的碎片效果更好
%对灰度图象进行“opening-by-reconstruction”运算
I = imcomplement(I);
se = strel('disk', 20);
Ie = imerode(I, se);
Iobr = imreconstruct(Ie, I);
figure, imshow(Iobr), title('Opening-by-reconstruction (Iobr)');%开幕重建
%对进行开运算后的图象进行“closing-by-reconstruction”运算
%对重建前后的图象都要进行求补运算
%Following the opening with a closing can remove the dark spots and stem marks
Iobrd = imdilate(Iobr, se);
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));
Iobrcbr = imcomplement(Iobrcbr);
figure, imshow(Iobrcbr), title('Opening-closing by reconstruction (Iobrcbr)');%开合重建
%计算Iobrcbr的局部最大值，以得到更好的前景目标
%fgm = imextendedmax(Iobrcbr,30);%30是自定义的数值，明显不具有普遍性，原来是：
fgm = imregionalmax(Iobrcbr);
figure, imshow(fgm), title('Regional maxima of opening-closing by reconstruction (fgm)');%区域极大的开放关闭的重建
%将前景目标叠加到原图上(在实际应用中可以去除这步）
I2 = I;
I2(fgm) = 255;
figure, imshow(I2), title('Regional maxima superimposed on original image (I2)');%区域最大值叠加在原始图像
% clean the edges of the marker blobs and then shrink them a bit
se2 = strel(ones(5,5));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);
%去除所有小于特定像素数（此处设定为20）的斑点污渍
fgm4 = bwareaopen(fgm3, 20);
I3 = I;
I3(fgm4) = 255;
figure, imshow(I3)
title('Modified regional maxima superimposed on original image (fgm4)')
%极大改进的区域叠加在原始图像
%%计算背景标记(外部约束)
%将基于重建的开闭运算图象转化为二值图象？？？？？？？？？？？？？？？？key
bw = im2bw(Iobrcbr, graythresh(Iobrcbr));
figure, imshow(bw), title('Thresholded opening-closing by reconstruction (bw)') %阈值开合的重建
%计算距离函数，产生外部约束
D = bwdist(bw);
%将不属于目标的像素置为—inf
%D = -D;
%D(bw) = -Inf;
DL = watershed(D);
bgm = DL == 0;
figure, imshow(bgm), title('Watershed ridge lines (bgm)')
%分水岭基线
%%分水岭算法分割
%重构梯度图
gradmag2 = imimposemin(gradmag, bgm | fgm4);
%分割
L = watershed(gradmag2);
figure,imshow(L),title('Wate