import numpy as np
import cv2
from matplotlib import pyplot as plt

def imshow(imgname,img):
    h ,w = img.shape[:2]
    cv2.namedWindow(imgname, cv2.WINDOW_NORMAL)
    cv2.resizeWindow(imgname, int(w * 0.3), int(h * 0.3))
    cv2.imshow(imgname, img)

def waterSeg(img):
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

    ret, thresh = cv2.threshold(gray,0,255,cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)

    thresh = cv2.bitwise_not(thresh)
    imshow('thresh',thresh)
    # noise removal
    # kernel = np.ones((3,3),np.uint8)
    # opening = cv2.morphologyEx(thresh,cv2.MORPH_OPEN,kernel, iterations = 2) # 形态开运算
    # sure_bg = cv2.dilate(opening,kernel,iterations=3)
    # cv2.imshow('sure_bg',sure_bg)

    # Finding sure foreground area
    dist_transform = cv2.distanceTransform(thresh,cv2.DIST_L2,5)
    # imshow('dist_transform',dist_transform)

    # print(0.7*dist_transform.max())

    ret, sure_fg = cv2.threshold(dist_transform,0.7*dist_transform.max(),255,cv2.THRESH_BINARY)

    imshow('sure_fg',sure_fg)
    # Finding unknown region
    sure_fg = np.uint8(sure_fg)

    unknown = cv2.subtract(thresh,sure_fg)

    imshow('unknown',unknown)
    # Marker labelling
    ret, markers = cv2.connectedComponents(sure_fg)
    # Add one to all labels so that sure background is not 0, but 1
    markers = markers + 1

    # Now, mark the region of unknown with zero
    markers[unknown==255] = 0

    markersw = cv2.watershed(img,markers)

    img[markersw == -1] = [255,0,0]

    imshow('img',img)

    cv2.waitKey(0)
    cv2.destroyAllWindows()

if __name__ == '__main__':
    # img = cv2.imread('20200608\\1592124020(1).jpg')
    img = cv2.imread('handwrite\\1592629465(1).jpg')
    waterSeg(img)

    该篇文章是在老师指导下完成的作业，感谢老师的指导。
    分水岭分割借鉴了形态学理论，是一种新的基于区域的图形分割算法。在该方法中，将一幅图像看成一个地形图，灰度值对应地形的高度值，高度值对应着山峰，低灰度值对应着山谷。水总是朝地势低的地方流动，直到某个局部低洼处，这个低洼处就是盆地，最终所有的水都会处于不停地盆地，盆地之间的山脊称为分水岭。
1、使用距离变换的分水岭分割：
距离变换：每一个像素到最邻近非零值像素的距离。
程序：
clear all;close all;clc
fobrcbr=imread('F:\matlab\MATLAB上机操作\图形\f1.tif');
level=graythresh(uint8(fobrcbr));
g=im2bw(fobrcbr,level);  %用阈值法把平滑后的图像变为二值图像
subplot(221)
imshow(g,[])
title('二值图像')
gc=~g;
subplot(222)
imshow(gc,[])
title('二值图像的补图像')
D=bwdist(gc);   %求二值图像的距离变换
subplot(223)
imshow(D,[])
title('距离变换')
L=watershed(-D);  %对负距离变换图求分水岭变换
w=L==0;  %求出山脊线
subplot(224)
imshow(w,[])
title('负距离变换的分水岭山脊线')
figure(2)
g2=g&~w;  %把黑色山脊线叠加在原图上
imshow(g2,[])

2、使用梯度的分水岭分割：
梯度大小图沿着物体的边缘有很高的像素值，在其他地方有低的像素值。分水岭变换将导致沿目标边沿的分水岭山脊线。
3、标记控制的分水岭分割：
由于噪声和梯度的其他局域不规则性，直接把分水岭变换用于梯度图通常导致过度分割。
使用标记概念控制过度分割的方法：
                      内部标记：在每个感兴趣的目标内。
                      外部标记：在背景内。
 
L=watershed(D);  %对图求分水岭变换
bwdist函数用于计算元素之间的距离。
[D,L]=bwdist(a)
D表示零元素所在的位置靠近非零元素位置的最短距离；
L则表示在该元素所靠近的最近的非零元的位置。
imcomplement为一个函数的名称，其可以对图像数据进行取反运算（实现底片效果）。
函数imcomplement可近似的把RGB模型转换为CMY模型：cmy_image=imcomplement(rgb_image)
也可以使用该函数将CMY图像转换为RGB图像：rgb_image=imcomplement(cmy_image)

分水岭算法

Watershed Algorithm(分水岭算法)，顾名思义，就是根据分水岭的构成来考虑图像的分割。现实中我们可以或者说可以想象有山有湖的景象，那么那一定是水绕 山，山围水的情形。当然在需要的时候，要人工构筑分水岭，以防集水盆之间的互相穿透。而区分高山（plateaus）与水的界线，以及湖与湖之间的间隔或 都是连通的关系，就是我们可爱的分水岭（watershed）。

如果图像中的目标物体是连接在一起的，则分割起来会更困难，分水岭分割算法经常用于处理这类问题， 

通常会取得比较好的效果。分水岭分割算法把图像看成一幅“地形图”， 

其中亮度比较强的区域像素值较大，而比较暗的区域像素值较小，通过寻找“汇水盆地”和“分水岭界限”，对图像进行分割。



案例

案例参考matlab官网案例，添加了详细注释，做出一定的调整，更容易让读者理解和接受。 
实现功能： 

将明显的梨从一堆梨子中分离出来

最终结果： 


文末将回答一下几个问题： 

（1） 如果采用最大类间方差阈值分割方法进行分割，效果如何？为什么？ 

（2） 直接用分水岭分割把“pears.png”分割好么？为什么？ 

（3） 如何获得前景标记？ 

（4） Imregionalmax是什么作用，请举例说明。 

（5） bwareaopen是什么作用，请举例说明。它是不是用数学形态学算法实现？ 

（6） 如何获得背景标记？ 

（7） 最终如何用前景标记和背景标记实现标记分水岭分割？

第一步：读入彩色图像，将其转化成灰度图像



clc; clear ; close ;
rgb = imread('pears.png');
if ndims(rgb) == 3
I = rgb2gray(rgb);
else
I = rgb;
end
figure('units', 'normalized','Name','图像读取：原图及灰度图比较');
subplot(1, 2, 1); imshow(rgb); title('原图');
subplot(1, 2, 2); imshow(I); title('灰度图');





第2步：方法1：将梯度幅值作为分割函数

使用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','Name','直接计算梯度幅值作为分割函数的结果');
subplot(1, 3, 1); imshow(I,[]), title('灰度图像')
subplot(1, 3, 2); imshow(gradmag,[]), title('梯度幅值图像')
L = watershed(gradmag);
Lrgb = label2rgb(L);
subplot(1, 3, 3); imshow(Lrgb); title('梯度幅值做分水岭变换')



 直接使用梯度模值图像进行分水岭算法得到的结果往往会存在过度分割的现象。因此通常需要分别对前景对象和背景对象进行标记，以获得更好的分割效果。



第3步：使用形态学技术“基于开的重建”和“基于闭的重建”来清理图像。



se = strel('disk', 20);
Io = imopen(I, se);
% 通过腐蚀后重建来做基于开的重建计算。
Ie = imerode(I, se);
Iobr = imreconstruct(Ie, I);
% 开操作后，接着进行闭操作，可以移除较暗的斑点和枝干标记。对比常规的形态学闭操作和基于闭的重建操作。
Ioc = imclose(Io, se);
Ic = imclose(I, se);
% 现在使用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','Name','比较各种形态学运算后的重建结果');
subplot(2, 3, 1); imshow(I, []); title('灰度图像');
subplot(2, 3, 2); imshow(Io, []); title('开操作图像');
subplot(2, 3, 3); imshow(Ic, []); title('闭操作图像');
subplot(2, 3, 4); imshow(Ioc, []), title('开闭操作');
subplot(2, 3, 5); imshow(Iobr, []); title('基于开的重建图像-腐蚀原图作为Marker');
subplot(2, 3, 6); imshow(Iobrcbr, []), title('基于闭的重建图像-膨胀开重建（取反）作为Marker');
% 通过比较Iobrcbr和loc可以看到，在移除小污点同时不影响对象全局形状的应用下，
% 基于重建的开闭操作要比标准的开闭重建更加有效。计算Iobrcbr的局部极大来得到更好的前景标记。




处理结果发现：基于重建的开闭操作要比标准的开闭重建更加有效。所以计算Iobrcbr的局部极大来得到更好的前景标记。



第4步：标记前景对象

有多种方法可以应用在这里来获得前景标记，这些标记必须是前景对象内部的连接斑点像素。 这些操作将会在每个对象内部创建单位极大值，使得可以使用imregionalmax来定位。



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

figure('units', 'normalized','Name','前景标记及调整');
subplot(3, 3, 1); imshow(I, []); title('灰度图像');
subplot(3, 3, 2); imshow(Iobrcbr, []); title('基于开闭的重建操作');
subplot(3, 3, 3); imshow(fgm, []); title('局部极大图像');
subplot(3, 3, 4); imshow(rgb, []); title('原图像');
subplot(3, 3, 5); imshow(I2); title('局部极大叠加到原图像');
% 注意到大多闭塞处和阴影对象没有被标记，这就意味着这些对象在结果中将不会得到合理的分割。
% 而且，一些对象的前景标记会一直到对象的边缘。
% 这就意味着应该清理标记斑点的边缘，然后收缩它们。可以通过闭操作和腐蚀操作来完成。
se2 = strel(ones(5,5));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);
subplot(3, 3, 6); imshow(fgm2, []); title('闭操作后-局部极大图像');
subplot(3, 3, 7); imshow(fgm3, []); title('腐蚀操作后-局部极大图像');
% 这个过程将会留下一些偏离的孤立像素，应该移除它们。可以使用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);
subplot(3, 3, 8); imshow(fgm4, []); title('进一步去除小斑点操作');
subplot(3, 3, 9); imshow(I3, []); title('修改后局部极大叠加到原图像');





第5步：计算背景标记



% 现在，需要标记背景。在清理后的图像Iobrcbr中，暗像素属于背景，所以可以从阈值操作开始。
bw =  imbinarize(Iobrcbr);
% 背景像素在黑色区域，但是理想情形下，不必要求背景标记太接近于要分割的对象边缘。
% 通过计算“骨架影响范围”来“细化”背景，或者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','Name','背景标记');
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('分水岭变换脊线图');





第6步：计算分割函数的分水岭变换



% 函数imimposemin可以用来修改图像，使其只是在特定的要求位置有局部极小。
% 这里可以使用imimposemin来修改梯度幅值图像，使其只在前景和后景标记像素有局部极小。
gradmag2 = imimposemin(gradmag, bgm | fgm4);
figure('units', 'normalized','Name','分割函数的分水岭变换对比');
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('修改后梯度幅值图像');





第7步：查看结果



% 一个可视化技术是叠加前景标记、背景标记、分割对象边界到初始图像。可
% 以使用膨胀来实现某些要求，比如对象边界，更加清晰可见。对象边界定位于L==0的位置。
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);
figure('units', 'normalized','Name','分割结果');
subplot(2, 2, 1); imshow(rgb, []); title('原图像');
subplot(2, 2, 2); imshow(I4, []); title('标记和对象边缘叠加到原图像');
% 另外一个有用的可视化技术是将标记矩阵作为彩色图像进行显示。标记矩阵，
% 比如通过watershed和bwlabel得到的，可以使用label2rgb转换到真彩图像来显示。
Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');
subplot(2, 2, 3); imshow(Lrgb); title('彩色分水岭标记矩阵');
% 可以使用透明度来叠加这个伪彩色标记矩阵在原亮度图像上进行显示。
subplot(2, 2, 4); imshow(rgb, []); hold on;
himage = imshow(Lrgb);
set(himage, 'AlphaData', 0.3);
title('标记矩阵叠加到原图像');



问题解答： 

（1） 如果采用最大类间方差阈值分割方法进行分割，效果如何？为什么？ 

最大类间方差阈值分割方法实际上是当做双峰分布，分割结果不理想。不适合做背景比较复杂的图像的分割



clc,clear
rgb=imread('pears.png');
I = rgb2gray(rgb);
imshow(I)

T=graythresh(I);%通过graythresh选择阈值
BW=im2bw(I,T);%用Otus阈值对图像进行分割
figure,imshow(BW);

（2） 直接用分水岭分割把“pears.png”分割好么？为什么？ 

直接分割会出现分割过度 

（3） 如何获得前景标记？ 

1.首先使用形态学技术“基于开的重建”和“基于闭的重建”来清理图像。发现基于开+闭的重建效果最好 

2.对重建后的图像在每个对象内部创建单位极大值，使得可以使用imregionalmax来定位 

3.这个过程将会留下一些偏离的孤立像素，应该移除它们。可以使用bwareaopen，用来移除少于特定像素个数的斑点。 

（4） Imregionalmax是什么作用，请举例说明。 

使得可以使用imregionalmax来定位极大值和极小值



A = 10*ones(10,10);
A(2:4,2:4) = 22; 
A(6:8,6:8) = 33; 
A(2,7) = 44;
A(3,8) = 45;
A(4,9) = 44





regmax = imregionalmax(A)



（5） bwareaopen是什么作用，请举例说明。它是不是用数学形态学算法实现？  

BW2 = bwareaopen(BW,P)从二值图像中移除所以少于P像素值的连通块，得到另外的二值图像BW2。 

感觉是形态学的方法 
官方Example 

（6） 如何获得背景标记？

% 现在，需要标记背景。在清理后的图像Iobrcbr中，暗像素属于背景，所以可以从阈值操作开始。
bw = im2bw(Iobrcbr, graythresh(Iobrcbr));
% 背景像素在黑色区域，但是理想情形下，不必要求背景标记太接近于要分割的对象边缘。
% 通过计算“骨架影响范围”来“细化”背景，或者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','Name','背景标记');
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('分水岭变换脊线图');

（7） 最终如何用前景标记和背景标记实现标记分水岭分割？ 

 函数imimposemin可以用来修改图像，使其只是在特定的要求位置有局部极小。 

 这里可以使用imimposemin来修改梯度幅值图像，使其只在前景和后景标记像素有局部极小。



gradmag2 = imimposemin(gradmag, bgm | fgm4);
L = watershed(gradmag2);

完整的代码：https://code.csdn.net/snippets/2603174


  
参考 
http://blog.sina.cn/dpool/blog/s/blog_725866260100rz7x.html

讲解以及实现
对分水岭的讲解

下面有对分水岭算法进行的实现的博客，可能两篇有出入，懂得大概，能用即可。

对分水岭的实现1

对分水岭的实现2

中间有一些函数的解释链接

如bwdist函数
下面是对三篇的引用

https://blog.csdn.net/qq_36808245/article/details/101562927

https://www.codeleading.com/article/5970803031/

https://blog.csdn.net/abcjennifer/article/details/6669074