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

参考文献:
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)

这一步就是实际调用了通用的分水岭算法,算法在这里
首先需要明确它的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,那么在此例中,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,即明确显示出来分水岭。
(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
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