数字旗手

电气化、自动化、数字化、智能化、智慧化

0%

ImagePy解析：15 -- 分水岭Watershed算法

Posted on Edited on In Disqus:

%%%%%%%%%%%%%%
2020-8-2更新：增加寻找邻居标识的过程解析
%%%%%%%%%%%%%%

参考文献：
OpenCV分水岭Watershed算法的前因后果
The Watershed Transformation

图像准备

此处分水岭算法所使用的样例图像为：
watershed-display
上述图像仅是为了便于展示，其实际尺寸为15乘以15像素，即像素矩阵为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
[[  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
[  0,   0, 255, 255, 255, 255, 255, 255, 255,   0,   0,   0,   0,    0,   0],
[  0,   0, 255, 255, 255, 255, 255, 255, 255,   0,   0,   0,   0,    0,   0],
[  0,   0, 255, 255, 255, 255, 255, 255, 255,   0,   0,   0,   0,    0,   0],
[  0,   0, 255, 255, 255, 255, 255, 255, 255,   0,   0,   0,   0,    0,   0],
[  0,   0, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,    0,   0],
[  0,   0, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,    0,   0],
[  0,   0, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,    0,   0],
[  0,   0,   0,   0,   0,   0, 255, 255, 255, 255, 255, 255, 255,     0,   0],
[  0,   0,   0,   0,   0,   0, 255, 255, 255, 255, 255, 255, 255,     0,   0],
[  0,   0,   0,   0,   0,   0, 255, 255, 255, 255, 255, 255, 255,     0,   0],
[  0,   0,   0,   0,   0,   0, 255, 255, 255, 255, 255, 255, 255,     0,   0],
[  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
[  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
[  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0]]

算法准备

分水岭算法需要配合其他的一些算法一块使用，才能获得较好的分割效果，即该算法需要一些前置算法的帮助，如获得种子标记等。
经过仔细对比，发现ImagePy的Process下的Binary下的Binary Watershed的流程较为简单，适合入手：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
class Watershed(Filter):
	"""Mark class plugin with events callback functions"""
	title = 'Binary Watershed'
	note = ['8-bit', 'auto_msk', 'auto_snap', 'preview']


	para = {'tor':2, 'con':False}
	view = [(int, 'tor', (0,255), 0, 'tolerance', 'value'),
			(bool, 'con', 'full connectivity')]

	def run(self, ips, snap, img, para = None):
		img[:] = snap>0
		dist = distance_transform_edt(snap, output=np.uint16)
		pts = find_maximum(dist, para['tor'], True)
		buf = np.zeros(ips.size, dtype=np.uint32)
		buf[pts[:,0], pts[:,1]] = img[pts[:,0], pts[:,1]] = 2
		markers, n = ndimg.label(buf, np.ones((3,3)))
		line = watershed(dist, markers, line=True, conn=para['con']+1, up=False)
		msk = apply_hysteresis_threshold(img, 0, 1)
		img[:] = snap * ~((line==0) & msk)

该算法的源程序在这里。
可以看出，该算法有两个参数tolerance和是否Full connectivity，以及需要前面的二值图。
另外，在Process的Hydrology下也有一个Find Watershed，它接收的参数多了sigma和thr，因为它接收一个灰度图，对其先做一个高斯滤波，然后根据阈值分割来确定标记种子。
这两处的分水岭原理都是一样的，都是调用了主文件夹下的ipyalg下的hydrology下的watershed算法，只是种子点的选取方式不同。

下面就是对二值分水岭Binary Watershed的run()函数的逐步解析。

二值化

1
img[:] = snap>0

前面的解析文章已经说过，snap是在图像被处理之前的一个快照。这一步是根据其是否大于0来对图像进行二值化操作，即0值保留，大于0的值都设为1。那么img经过二值化操作后，就变为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8)

注意，snap的值不变，最后一步还要用到它。

距离变换

1
dist = distance_transform_edt(snap, output=np.uint16)

这一步是得到距离变换图，具体距离变换的原理留待后续解析，这里直接给出距离变换的结果，即dist的值为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 2, 2, 2, 2, 2, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 2, 3, 3, 3, 2, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 2, 3, 3, 3, 2, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 2, 3, 3, 3, 2, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 2, 3, 3, 3, 2, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint16)

可以看出，两个白色方块的边缘上的距离值都是1，方块中心的距离值都是4，即计算结果为某一点与背景像素点（值为0）的棋盘距离（也许这里函数后缀为edt不合适？意味着是欧几里得距离？）。

寻找局部极值

1
pts = find_maximum(dist, para['tor'], True)

这一步是根据上面的距离变换图寻找局部极值（这里是极大值），具体的寻找原理见之前的解析，在这里。
这里的tolerance参数非常重要，它决定着每个局部极大值的“领地”，如果设置不当，就有可能找不到。针对这张图，因为其尺寸很小，这里只能设置为1。
经过计算后，pts的值为：

1
array([[4, 5], [8, 9]], dtype=int16)

即这两个白色方块的中心点的坐标位置。

标记局部极值的位置

1
2
buf = np.zeros(ips.size, dtype=np.uint32)
buf[pts[:,0], pts[:,1]] = img[pts[:,0], pts[:,1]] = 2

首先创建了一个与图像同样大小的缓冲区，初始值全部为0，然后根据前面得到的局部极大值的坐标位置，将缓冲区和二值化后的img的该位置处的像素值设为2。
那么，此时buf的值为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint32)

img的值为：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 2, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8)

标识种子（根据极值点位置进行标号）

1
markers, n = ndimg.label(buf, np.ones((3,3)))

这里特别注意scipy的ndimage包的label函数，它是为了标记矩阵中的特征，比如标识有多少个连通区域。
这里使用了一个3乘以3的、且值都是1的核，表示只要八邻域内有值即表示两者是同一个特征。
经过特征标注后，n的值为2，表示有两个种子点。
markers的值为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

分水岭

1
line = watershed(dist, markers, line=True, conn=para['con']+1, up=False)

这一步就是实际调用了通用的分水岭算法，算法在这里。
这里计算分水岭的算法基于浸没模拟（Immersion Simulation)思想：在种子点位置，刺穿一个小孔，或称打一个洞，然后把整个模型慢慢浸入水中，随着浸入的加深，种子点的影响区域慢慢向外扩展，在两个集水盆汇合处构筑大坝，即形成分水岭。由于修建的大坝将阻止聚合，即，浸没结束时，所建立的堤坝就对应于区域的轮廓，而集水盆则对应分割区域。

“高”和“低”是一个相对的概念，上面思想是在两个山谷处注水，从而找到“分水岭”；同样地，也可以将原图像进行“反向”，使得山顶变为山谷，然后在原“山顶”（即现”山谷“）处进行注水（注意，注水都是在山谷处进行），那么就得到了现图像的”分水岭“，效果就是找到了原图像的”河谷线“。

首先需要明确它的inputs是什么：距离变换图dist、种子标识markers、是否标识分水岭line、邻域范围conn、是否向上注水up。
具体解释一下：
距离变换图：之所以选择距离变换图作为分水岭算法的输入，是因为它是一个有高低起伏的地形图，必须要依据它进行“注水建大坝”；
种子标识：这里的种子是局部极大值点，种子也可以是一个区域；
是否标识分水岭：如果要标识出分水岭，就会在分水岭处设置一个大数，反之，两个连通域就硬碰硬；
邻域范围：这个参数决定了中心像素的邻域范围是多大，如果con为False，那么conn就是1，就是四邻域，如果con为True，那么conn就是2，就是八邻域；
是否向上注水：该参数与选择的种子点密切相关，即如果标记的是山峰位置，那么up就得是False，即向下下雨；如果标记的是盆地位置，那么up就是True，即向上注水。

变换数据范围

1
2
3
4
def watershed(img, mark, conn=1, line=False, up=True):
    maxv = img.max(); minv = img.min();
    if img.dtype != np.uint8:
        img = ((img-minv)*(255/(maxv-minv))).astype(np.uint8）

注意，上面的img形参就是传入的距离图dist，它的dtype是uint16，它最大值是4，最小值是0，上面的代码将其范围转换为0到255之间。具体做法就是先取得最大值和最小值，然后将其归一化，再乘以255，那么img的值变为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
array([[  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,    0,   0],
       [  0,   0,  63,  63,  63,  63,  63,  63,  63,   0,   0,   0,   0,    0,   0],
       [  0,   0,  63, 127, 127, 127, 127, 127,  63,   0,   0,   0,   0,    0,   0],
       [  0,   0,  63, 127, 191, 191, 191, 127,  63,   0,   0,   0,   0,    0,   0],
       [  0,   0,  63, 127, 191, 255, 191, 127,  63,   0,   0,   0,   0,    0,   0],
       [  0,   0,  63, 127, 191, 191, 191, 127,  63,  63,  63,  63,  63,   0,   0],
       [  0,   0,  63, 127, 127, 127, 127, 127, 127, 127, 127, 127,  63,    0,   0],
       [  0,   0,  63,  63,  63,  63,  63, 127, 191, 191, 191, 127,  63,    0,   0],
       [  0,   0,   0,   0,   0,   0,  63, 127, 191, 255, 191, 127,  63,    0,   0],
       [  0,   0,   0,   0,   0,   0,  63, 127, 191, 191, 191, 127,  63,    0,   0],
       [  0,   0,   0,   0,   0,   0,  63, 127, 127, 127, 127, 127,  63,    0,   0],
       [  0,   0,   0,   0,   0,   0,  63,  63,  63,  63,  63,  63,  63,    0,   0],
       [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,    0,   0],
       [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,    0,   0],
       [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,    0,   0]], dtype=uint8)

四周填充

1
2
3
4
5
6
img = buffer(img, np.uint8)

def buffer(img, dtype):
    buf = np.zeros(tuple(np.array(img.shape)+2), dtype=dtype)
    buf[tuple([slice(1,-1)]*buf.ndim)] = img
    return buf

这一步是将img的四周用0填充，即原来img是15乘以15大小，填充后的img是17乘以17大小，且具体数值为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
array([[  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,    0,   0,   0,   0],
       [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,    0,   0,   0,   0],
       [  0,   0,   0,  63,  63,  63,  63,  63,  63,  63,   0,   0,   0,    0,   0,   0,   0],
       [  0,   0,   0,  63, 127, 127, 127, 127, 127,  63,   0,   0,   0,    0,   0,   0,   0],
       [  0,   0,   0,  63, 127, 191, 191, 191, 127,  63,   0,   0,   0,    0,   0,   0,   0],
       [  0,   0,   0,  63, 127, 191, 255, 191, 127,  63,   0,   0,   0,    0,   0,   0,   0],
       [  0,   0,   0,  63, 127, 191, 191, 191, 127,  63,  63,  63,  63,  63,   0,   0,   0],
       [  0,   0,   0,  63, 127, 127, 127, 127, 127, 127, 127, 127, 127, 63,   0,   0,   0],
       [  0,   0,   0,  63,  63,  63,  63,  63, 127, 191, 191, 191, 127, 63,   0,   0,   0],
       [  0,   0,   0,   0,   0,   0,   0,  63, 127, 191, 255, 191, 127, 63,   0,   0,   0],
       [  0,   0,   0,   0,   0,   0,   0,  63, 127, 191, 191, 191, 127, 63,   0,   0,   0],
       [  0,   0,   0,   0,   0,   0,   0,  63, 127, 127, 127, 127, 127, 63,   0,   0,   0],
       [  0,   0,   0,   0,   0,   0,   0,  63,  63,  63,  63,  63,  63, 63,   0,   0,   0],
       [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
       [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
       [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
       [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0]], dtype=uint8)

同理，对于种子标识点也做同样的填充：
1
mark = buffer(mark, np.uint32)

mark即变为：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint32)

但注意两者的数据类型dtype，前者是uint8类型，后者是uint32类型。
然后对mark来说，将四周填充的0改为2的32次方减1（用十六进制表示就是0xffffffff），这样就明确标识出了四周边界（后面还会用4294967294，即0xfffffffe，来标识出内部边界，即分水岭）：
1
2
3
4
ndim = img.ndim
for n in range(ndim):
    idx = [slice(None) if i!=n else [0,-1] for i in range(ndim)]
    mark[tuple(idx)] = 0xffffffff

那么，mark变为：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
array([[4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          1,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            2,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295,          0,          0,          0,          0,           0,          0,          0,          0,          0,            0,          0,          0,          0,          0,              0,             4294967295],
       [4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295]], dtype=uint32)

寻找邻居标识并压平图像

1
nbs = neighbors(img.shape, conn)

该步是计算某个像素的邻居点的标识，比如它左侧的邻居像素标识就是-1，右侧的邻居标识就是1。neighbors函数的具体写法之所以看着繁琐，是为了适配任意维度，保证程序的鲁棒性。
另外，这里用到了整个分水岭算法的第四个参数conn，即Full Connectivity，实测该参数是控制计算四邻域还是八邻域。默认full connectivity是False。
寻找邻居标识的具体过程为：

1
block = generate_binary_structure(dim, conn)

即先生成一个二元的结构单元，这个结构单元的大小和内容由维度和连通性来决定。如果连通性为1，即寻找的是4邻域，另外，若dim为2，即图像的shape有二维，那么该结构单元为：
1
2
3
[[False  True False]
[ True  True  True]
[False  True False]]

若dim为3，且连通性仍为1，则：
1
2
3
4
5
6
7
8
9
[[[False False False]
  [False  True False]
  [False False False]]
[[False  True False]
  [ True  True  True]
  [False  True False]]
[[False False False]
  [False  True False]
  [False False False]]]

这里可以将其联想成一个3乘3的魔方，二维时候的4邻域就变成了三维时的6邻域，此时中心像素与邻居的距离都是1个棋盘距离；同理，二维时的8邻域就变成了三维时的18邻域，这时中心像素与邻居的距离都是2个棋盘距离以内；同时，三维时还有魔方的八个角点也作为邻居时的情形，此时中心像素与邻居的距离都是3个棋盘距离以内，即此时三维的邻域就是26邻域。

下面的解析都以二维且4邻域为例：

1
block[tuple([1]*dim)] = 0

这一步是将中心像素的标识置为0，因此，此时block变为：
1
2
3
[[False  True False]
[ True  Flase  True]
[False  True False]]

接下来取得这些邻居所在的索引指标：
1
idx = np.where(block>0)

idx的输出为：
1
idx =  (array([0, 1, 1, 2], dtype=int64), array([1, 0, 2, 1], dtype=int64))

因此idx的行标号和列标号是分离的，还需要将它们组合起来：
1
idx = np.array(idx, dtype=np.uint8).T

一个简单的转置操作即可实现，此时idx为：
1
2
3
4
idx =  [[0 1]
[1 0]
[1 2]
[2 1]]

然后：
1
idx = np.array(idx-[1]*dim)

这一步是取得邻居与中心像素的相对距离，此时idx变为：
1
2
3
4
idx =  [[-1  0]
[ 0 -1]
[ 0  1]
[ 1  0]]

注意，此时这个相对距离还是以“图像”为载体，即是一个矩阵上的行列的相对距离，实际在使用时，是将图像压平为一个很长的一维链，所以下一步是根据图像的宽度将这里的相对距离转换为一维链上的相对距离：
1
acc = np.cumprod((1,)+shape[::-1][:-1])

这一步是最重要的一步：
（1）对shape的操作是先将它翻转，比如原来是(10, 110)，先翻转为(110, 10)，然后排除掉最后一个元素，即只取110，这是因为在一维链上的相对位置与图像宽度有关，而与高度无关；如果shape是(10, 110, 3)，即是一个三通道图像，那么排除最后一个元素后，就是(3, 110)。
（2）接着是另一个精髓的操作，将(1,)这个元组与shape截取后的元组相加，注意元组相加其实是组合效果，即比如(1,)+(110, )等于(1, 110)，以及(1, )+(3, 110)等于(1, 3, 110)，这里之所以用(1,)相加，是为了对应与中心像素相同高度上的邻居。
（3）最后是另一个神来之笔，即用numpy的累乘，这一步是用来将通道考虑进去。

经过一系列操作，如果是二维的shape(10, 110)，那么acc为：

1
acc =  [1 110]

如果是三维的shape(10, 110, 3)，那么acc为：
1
acc =  [1  3 330]

至此，就知道了在一维链上邻居点相对中心像素的绝对距离，那么，针对于具体的邻居点，得到其具体的绝对距离：
1
np.dot(idx, acc[::-1])

比如二维时结果为：
1
array([-110,   -1,    1,  110])

三维时结果为：
1
array([-330,   -3,   -1,    1,    3,  330])

以上就是寻找邻居标识的全部过程。

在此例中，nbs就是：

1
array([-17,  -1,   1,  17])

如果full connectivity设为True，那么nbs的值就是：
1
array([-18, -17, -16,  -1,   1,  16,  17,  18])

接下来就是将img和mark进行压平，这样上面的邻居标识就可以配合压平后的图像进行使用：
1
2
img = img.ravel()
mark1d = mark.ravel()

这样img和mark1d的shape都是：
1
(289, )

统计灰度直方图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
pts = np.zeros(img.size//3, dtype=np.int64)
s, bins = collect(img, mark1d, nbs, pts)

def collect(img, mark, nbs, pts):
    bins = np.zeros(img.max()+1, dtype=np.uint32)
    cur = 0
    for p in range(len(mark)):
        bins[img[p]] += 1
        if mark[p]==0xffffffff: continue # edge
        if mark[p]==0: continue      # zero
        for dp in nbs:
            if mark[p+dp]!=mark[p]:
                pts[cur] = p
                cur += 1
                break
    return cur, bins

这一步是统计了图像的灰度直方图bins，即每一个灰度阶有多少个像素。s存储的是后面注水时一步步水漫时的中心像素的个数。
bins的数值是：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
array([200,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,  40,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,  31,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,  16,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   2], dtype=uint32)

即灰度为0的像素有200个，灰度为63的像素有40个，灰度为127的像素有31个，灰度为191的像素有16个，灰度为255的像素有2个。
cur就是种子数，为2。
同时，该函数也改变了pts的数值，它存储了后面注水时一步步水漫时的中心像素的位置，初始值是两个种子点的位置：
1
2
3
4
5
6
7
8
array([ 91, 163,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0], dtype=int64)

注水

1
2
3
4
for level in range(len(bins))[::1 if up else -1]:
    if bins[level]==0:continue
    s, c = clear(pts, s, 0)
    s = step(img, mark1d, line, pts, s, level, up, nbs)

首先，这里面用到了最开始调用分水岭算法时的最后一个参数up，up控制的是灰度的遍历顺序，即从0到255，还是从255到0，即：

1
range(len(bins))[::1 if up else -1]

若up为True，则1为步长，该range就是range(0,256)；若up为False，则-1为步长，该range就是range(255, -1, -1)。

然后对所有的灰度阶（即水平面）进行判断：
（1）如果该灰度阶没有对应的像素，那么就跳过下面的动作，继续对下一个灰度阶进行判断；
（2）对于有对应像素的灰度阶，首先清除pts中无效标识，即clear函数所做的动作，如果pts中的标识为-1，则清除，且将后面的标识往前推；
（3）接着就是具体的注水过程，执行step()函数，它接收的参数比较多：img是地形图，mark1d是标识点，line是是否标识分水岭，pts中存储了水漫时的中心位置点，s是中心位置点个数，level是当前水平面，up是是否注水还是下雨，nbs是邻居像素的标识。
（3.1）对水漫时的中心像素进行循环，依次从pts中取出该像素的位置标识，首先判断此处的标记值是否为0xfffffffe，或者是否向上注水up和此处img灰度值是否大于该灰度阶，或者是否向下注水not up和此处img灰度值是否小于该灰度阶。
这里需要明确一下Python的逻辑运算符的优先级，这个链接是非常清晰的一个解释。摘抄如下：

首先，‘and’、‘or’和‘not’的优先级是not大于and大于or；
其次，逻辑操作符and 和or 也称作短路操作符（short-circuitlogic）或者惰性求值（lazy evaluation）：它们的参数从左向右解析，一旦结果可以确定就停止。
例如，如果A 和C 为真而B 为假， A and B and C 不会解析C 。作用于一个普通的非逻辑值时，短路操作符的返回值通常是最后一个变量。
因此，and运算符必须所有的运算数都是true才会把所有的运算数都解析，并且返回最后一个变量。而或逻辑（or），即只要有一个是true，即停止解析运算数，返回最近为true的变量。

（3.2）如果上述条件满足，那么就跳过此处的像素，继续对pts中的位置进行循环，这一步实现的效果就是：（a）如果是向上注水，即up为True，那么代表在集水盆地注水，此时在这个盆地内与注水点相同水平面（即同一个level值）的像素点的标记值会被置为与注水点相同的标记值，而当某一像素值大于此时的level时，就不会被淹没，因为它的水平面更高，这是为了保证分水岭两侧的水平面是同步上涨的，否则就像是水会沿着山体爬上去，然后就会越过山体与另一侧汇合，这样也就无法找到正确的分水岭；（b）如果是向下注水，即up为False，那么代表在山顶开始下雨，水位也是在同一个level时停止，保证水在两个山体下滑的速度相等；（c）如果该像素已经是分水岭，即标记值为0xfffffffe，那么也跳过。

1
2
3
if msk[p] == 0xfffffffe or up and img[p]>level or not up and img[p]<level:
    cur += 1
    continue

（3.3）如果上述条件不满足，则对该中心像素的邻居像素进行循环：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
for dp in nbs:
    cp = p+dp
    if msk[cp]==0:
        msk[cp] = msk[p]
        if s == len(pts):
            s, cur = clear(pts, s, cur)
        pts[s] = cp
        s+=1

    elif msk[cp] == msk[p]:
        continue
    elif msk[cp] == 0xffffffff:
        continue
    elif msk[cp] == 0xfffffffe:
        continue

    elif line : msk[cp] = 0xfffffffe

首先判断该邻居的标记值是否为0，如果是，则将该中心像素的标记值赋给它，同时将该邻居位置记录进入pts中，后面将以该点为中心继续填充；如果不为0，再判断该邻居像素与中心像素的标记值是否相等，若相等，则跳过下面步骤；若两者不相等，则再判断该邻居像素的标记值是否为0xffffffff，即是否在边界上，若在边界上，即跳过下面步骤；若不在边界上，则继续判断该像素的标记值是否为0xfffffffe（该值是为了标识分水岭位置），若是分水岭，则跳过以下步骤；若不是分水岭，则继续判断line是否为True，该参数控制是否将分水岭标识出来，若为True，则此时将该邻居像素的标记值置为0xfffffffe，即明确显示出来分水岭。
因此，这里面的逻辑就是：
0代表还没占领的，0xffffffff代表边界，0xfffffffe代表已经确定是分水岭的，msk[p]代表自己已经占领的，如果这几个都不满足或命中，则言外之意是到了他人占领的地方了，此时如果设置line=True，那么就划定界线，将该像素标为0xfffffffe，即分水岭。

（3.4）对该中心像素的邻居遍历完后，还将该中心像素在pts中的标识转为-1，方便后面在clear函数中进行删除。

经过上面的注水过程，得到的标识图像mark为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
array([[4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295],
       [4294967295,          1,          1,          1,          1,          1,          1,          1,          1,          1,         1,          1,          1,          1, 4294967294, 2, 4294967295],
       [4294967295,          1,          1,          1,          1,          1,          1,          1,          1,          1,         1,          1,          1,          4294967294,  2,  2, 4294967295],
       [4294967295,          1,          1,          1,          1,          1,          1,          1,          1,          1,         1,          1, 4294967294,          2,          2,   2, 4294967295],
       [4294967295,          1,          1,          1,          1,          1,          1,          1,          1,          1,         1, 4294967294,          2,          2,          2,   2, 4294967295],
       [4294967295,          1,          1,          1,          1,          1,          1,          1,          1,          1,         4294967294,          2,          2,          2,      2,  2, 4294967295],
       [4294967295,          1,          1,          1,          1,          1,          1,          1,          1, 4294967294,   2,          2,          2,          2,          2,    2, 4294967295],
       [4294967295,          1,          1,          1,          1,          1,          1,          1, 4294967294,          2,   2,          2,          2,          2,          2,   2, 4294967295],
       [4294967295,          1,          1,          1,          1,          1,          1, 4294967294,          2,          2,   2,          2,          2,          2,          2,   2, 4294967295],
       [4294967295,          1,          1,          1,          1,          1, 4294967294,          2,          2,          2,   2,          2,          2,          2,          2,   2, 4294967295],
       [4294967295,          1,          1,          1,          1,  4294967294,          2,          2,          2,          2,        2,          2,          2,          2,          2,    2, 4294967295],
       [4294967295,          1,          1,          1, 4294967294,          2,          2,          2,          2,          2,     2,          2,          2,          2,          2,        2, 4294967295],
       [4294967295,          1,          1, 4294967294,          2,          2,          2,          2,          2,          2,     2,          2,          2,          2,          2,        2, 4294967295],
       [4294967295,          1, 4294967294,          2,          2,          2,          2,          2,          2,          2,     2,          2,          2,          2,          2,        2, 4294967295],
       [4294967295, 4294967294,          2,          2,          2,          2,          2,          2,          2,          2,     2,          2,          2,          2,          2,        2, 4294967295],
       [4294967295,          2,          2,          2,          2,      2,          2,          2,          2,          2,          2,          2,          2,          2,          2,     2, 4294967295],
       [4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295, 4294967295]], dtype=uint32)

蚀刻边界

1
if line:erose(mark1d)

erose函数将边界和分水岭边界给去掉，此时mark变为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 2, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 2, 2, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 2, 2, 2, 0],
       [0, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0],
       [0, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0],
       [0, 1, 1, 1, 1, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0],
       [0, 1, 1, 1, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0],
       [0, 1, 1, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0],
       [0, 1, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0],
       [0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0],
       [0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint32)

去除四周边界填充

最后，将之前加入的四周边界填充去掉：

1
mark[tuple([slice(1,-1)]*ndim)]

那么，mark即变为：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
array([[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 2],
       [1, 1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 2, 2],
       [1, 1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 2, 2, 2],
       [1, 1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 2, 2, 2, 2],
       [1, 1, 1, 1, 1, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2],
       [1, 1, 1, 1, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
       [1, 1, 1, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
       [1, 1, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
       [1, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
       [0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
       [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]], dtype=uint32)

至此，整个通用的分水岭算法就结束了。

后处理显示

这一步就是为了将分水岭可视化出来，根据不同的问题，有不一样的处理方式，比如这里的二值分水岭和另一个灰度图分水岭的处理方式就不同。但目的是一致的，就是可视化。
还是看这里的二值分水岭的处理方法：

创建原图掩膜

1
msk = apply_hysteresis_threshold(img, 0, 1)

这一步是基于已经标识了极大值位置的img来创建掩膜，注意不是上面的mark，而是img，即：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 2, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8)


apply_hysteresis_threshold()函数的用法是：
1
2
3
4
skimage.filters.apply_hysteresis_threshold(image, low, high)

Apply hysteresis thresholding to image.
This algorithm finds regions where image is greater than high OR image is greater than low and that region is connected to a region greater than high.

即找到图像中的这样的区域：（1）比high值要大的区域，或者（2）比low值要大，但该区域要连接比high值大的区域。
那么，掩膜msk的值就是：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
array([[False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False,  True,  True,  True,  True,  True,  True,  True, False, False, False, False, False, False],
       [False, False,  True,  True,  True,  True,  True,  True,  True, False, False, False, False, False, False],
       [False, False,  True,  True,  True,  True,  True,  True,  True, False, False, False, False, False, False],
       [False, False,  True,  True,  True,  True,  True,  True,  True, False, False, False, False, False, False],
       [False, False,  True,  True,  True,  True,  True,  True,  True, True,  True,  True,  True, False, False],
       [False, False,  True,  True,  True,  True,  True,  True,  True, True,  True,  True,  True, False, False],
       [False, False,  True,  True,  True,  True,  True,  True,  True, True,  True,  True,  True, False, False],
       [False, False, False, False, False, False,  True,  True,  True, True,  True,  True,  True, False, False],
       [False, False, False, False, False, False,  True,  True,  True, True,  True,  True,  True, False, False],
       [False, False, False, False, False, False,  True,  True,  True, True,  True,  True,  True, False, False],
       [False, False, False, False, False, False,  True,  True,  True, True,  True,  True,  True, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False]])

与分水岭做与操作

1
(line==0) & msk

将上述掩膜与分水岭边界进行与操作，就可以得到分水岭的一个子集，即两个矩形区块相碰的区域，注意看中间的三个True值：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
array([[False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False,  True, False, False, False, False, False, False],
       [False, False, False, False, False, False, False,  True, False, False, False, False, False, False, False],
       [False, False, False, False, False, False,  True, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False],
       [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False]])

然后再取反，再与原先的截图snap进行相乘：
1
img[:] = snap * ~((line==0) & msk

从而得到原图被分水岭分割后的效果，即：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
array([[  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
       [  0,   0, 255, 255, 255, 255, 255, 255, 255,   0,   0,   0,   0,   0,   0],
       [  0,   0, 255, 255, 255, 255, 255, 255, 255,   0,   0,   0,   0,   0,   0],
       [  0,   0, 255, 255, 255, 255, 255, 255, 255,   0,   0,   0,   0,   0,   0],
       [  0,   0, 255, 255, 255, 255, 255, 255, 255,   0,   0,   0,   0,   0,   0],
       [  0,   0, 255, 255, 255, 255, 255, 255,   0, 255, 255, 255, 255, 0,   0],
       [  0,   0, 255, 255, 255, 255, 255,   0, 255, 255, 255, 255, 255, 0,   0],
       [  0,   0, 255, 255, 255, 255,   0, 255, 255, 255, 255, 255, 255, 0,   0],
       [  0,   0,   0,   0,   0,   0, 255, 255, 255, 255, 255, 255, 255,  0,   0],
       [  0,   0,   0,   0,   0,   0, 255, 255, 255, 255, 255, 255, 255,  0,   0],
       [  0,   0,   0,   0,   0,   0, 255, 255, 255, 255, 255, 255, 255,  0,   0],
       [  0,   0,   0,   0,   0,   0, 255, 255, 255, 255, 255, 255, 255,  0,   0],
       [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
       [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
       [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0]], dtype=uint8)

就是这个样子：
watershed-after-display

山脊线

与分水岭不同，山脊线可以形成孤立的线，线段基本是沿着地理意义上的山脊行走，除了主脉，也有余脉。
具体寻找山脊线时，类似分水岭算法，仍然采用“涨水法”，但只用一个集水盆，并且制定新的终止规则：
（1）不能淹没终端像素；
（2）不能将原本相连的区域分成两个孤岛。
image