1 算法简介

分水岭算法的原理很容易查到，但是很多文章都是直接用的opencv或matlab函数，看不到具体实现方法，这篇文章希望能对大家有点帮助。

分水岭算法就是往山谷中注水，把不同湖水接触的位置称作分水岭，这么做的结果有两个：
1.得到整个区域的分界线
2.把整个区域编号

这是一种很容易过分割的算法，需要一些预处理手段和一些优化，本文只解析基本原理
（下图的注水过程是一种基于人为标记的方法）

2 代码步骤说明

分水岭算法的步骤如下：

1.将图像的所有像素按像素值从小到大排序，这里可以利用直方图将像素信息塞入数组

2.开始按灰度级从小到大顺序遍历所有像素，先将该灰度级的全部像素标记为待计算点，若点的邻域内有已存在的水池，则放入一个队列（先将水池边缘像素入队）

3.开始遍历队列，直到队列为空

3.1 若当前像素周围没有已经编号过的像素，则是新的水池
给与一个新的编号，并将邻域内的同灰度级像素加入队列
一直遍历队列到队列为空，即将这个新水池底部填满3.2 若像素周围有已经计算过的像素则看情况：
3.2.1 如果周围有2个不同的水池编号，则是分水岭
3.2.2 如果周围只有1个水池编号，则该点也是这个编号
3.2.3 周围有同层（待计算点），入队

4.开始计算下一个灰度级

所以这个算法就是从每个水池的边界开始一圈一圈，一个灰度级一个灰度级地向外蔓延的过程

3 python实现

import cv2
import numpy as np
from queue import Queue#队列maxhigh = 255 #水位
mask = -2       #用于标记每次涨水时，将会被水淹没的像素
watershed = 0   #用于标记分水岭边缘
#fixmark = np.array([-1,-1])#用于隔开队列的每个部分def checkedge(inputidx,w,h):if inputidx[0] >= h:return -1if inputidx[1] >= w:return -1if inputidx[0] < 0:return -1if inputidx[1] < 0:return -1return 0#def mark_end_que(que):# que.put(fixmark) def water(inputimg, size):imsize = size[0]*size[1]w = size[1]h = size[0]pixarray = np.zeros([imsize,2])#用于储存所有像素，以及坐标labelout = np.zeros(size)-1#distance = np.zeros(size)putquemask = np.zeros(size)labelsize = np.zeros(imsize)cnt = np.zeros(257)#累积直方图currenLabel = 0#湖区标号que = Queue(maxsize=0)#创建队列neighborbias4 = np.array([[0,-1],[-1,0],[0,1],[1,0]])#邻域偏移#neighborbias8 = np.array([[0,-1],[-1,-1],[-1,0],[-1,1],[0,1],[1,1],[1,0],[1,-1]])#邻域偏移#计算直方图，知道每个灰度级的位置hist = cv2.calcHist([inputimg],[0],None,[256],[0,256])#累计直方图的点就是每个灰度级在pixarray中的起始位置，但是要去除空的点hist[0] = hist[0] - 1#由于累积直方图是作为下标使用，因此最终求和结果不可等于imsizehistsum = 0for i in range(0,256):histsum = histsum + hist[i]cnt[i+1] = histsumcnt = cnt.astype(np.int32)cntidx = cnt[1:257].copy()#遍历图片记录全部像素，根据累积直方图排列像素顺序，到时候按顺序从低到高一个个涨水for y in range(0,h):for x in range(0,w):pix = inputimg[y,x]pixarray[(cntidx[pix]),:] = y,xcntidx[pix] = cntidx[pix] - 1#准备完成，开始涨水for nowgraylevel in range(0,maxhigh+1):print(nowgraylevel)#标记当前层灰度for idx in range(cnt[nowgraylevel],cnt[nowgraylevel+1]+1):nowpixaxis = pixarray[idx].astype(np.int32)labelout[nowpixaxis[0],nowpixaxis[1]] = mask#把在水池边缘的当前灰度级入队for nei in range(0,4):neighboraxis  = nowpixaxis + neighborbias4[nei]if checkedge(neighboraxis,w,h) < 0:continueif labelout[neighboraxis[0],neighboraxis[1]] >= 0:#周围有已经计算过的putquemask[neighboraxis[0],neighboraxis[1]] > 0#标记已经入队过了que.put(nowpixaxis)#入队break#开始遍历队列while(True):if que.empty():breaknowpixaxis = que.get()#蔓延过程#1 周围有1个label，加入其中#2 周围两个不同label，分水岭#3 周围无label，新水池#4 周围有同层，入队蔓延for nei in range(0,4):neighboraxis  = nowpixaxis + neighborbias4[nei]if checkedge(neighboraxis,w,h) < 0:continuelabelnow = labelout[nowpixaxis[0],nowpixaxis[1]]labelneighbor = int(labelout[neighboraxis[0],neighboraxis[1]])if labelnow == -2 and labelneighbor > 0:labelout[nowpixaxis[0],nowpixaxis[1]] = labelneighborlabelsize[labelneighbor] = labelsize[labelneighbor] + 1elif labelnow > 0 and labelneighbor > 0 and labelnow != labelneighbor:labelout[nowpixaxis[0],nowpixaxis[1]] = 0if labelneighbor == -2 and putquemask[neighboraxis[0],neighboraxis[1]] == 0:que.put(neighboraxis)putquemask[neighboraxis[0],neighboraxis[1]] = 1#标记这个像素已经入队了不要重复使用#找到了新的水坑（邻域没有水）for idx in range(cnt[nowgraylevel],cnt[nowgraylevel+1]+1):nowpixaxis = pixarray[idx,:].astype(np.int32)if(labelout[nowpixaxis[0],nowpixaxis[1]] == -2):#经过之前的赋值，仍然没有被标记的是新水池#print("new pool %d",)currenLabel = currenLabel + 1#配置序号labelout[nowpixaxis[0],nowpixaxis[1]] = currenLabelque.put(nowpixaxis)#将新坑底所有同灰度级像素都填上水while not que.empty():nowpixaxis = que.get()for nei in range(0,4):neighboraxis = nowpixaxis + neighborbias4[nei]#防越界if checkedge(neighboraxis,w,h) < 0:continueif  putquemask[neighboraxis[0],neighboraxis[1]] == 0 and labelout[neighboraxis[0],neighboraxis[1]] == -2:labelout[neighboraxis[0],neighboraxis[1]] = currenLabelque.put(neighboraxis)putquemask[neighboraxis[0],neighboraxis[1]] = 1labelsize[currenLabel] = labelsize[currenLabel] + 1return labelout    image = cv2.imread('F:/tf_learn--------/impo/coin3.png')  #F:/tf_learn--------/impixiv/1606.jpg #F:/tf_learn--------/impo
gray = cv2.cvtColor(image,cv2.COLOR_BGR2GRAY, cv2.CV_16S)size = gray.shape  # h w c#二值化
ret,gray = cv2.threshold(gray,40,255,cv2.THRESH_BINARY)
cv2.imshow("ret",gray)
cv2.waitKey(0)
#膨胀
kernel = np.ones((5,5),np.uint8)
gray = cv2.dilate(gray,kernel)#腐蚀
gray = cv2.erode(gray,kernel)#边缘检测
#edgeimage = np.uint8(np.abs(cv2.Sobel(gray,cv2.CV_16S, 1, 1, ksize=3)))#这个sobel不太方便
edgeimage = cv2.Canny(gray, 10, 100) cv2.imshow("gray",edgeimage)
cv2.waitKey(0)label = water(edgeimage,size)
label = cv2.resize(label,(size[1],size[0]))#标记
b=image[:,:,0]
g=image[:,:,1]
r=image[:,:,2]
r[label==0] = 255
g[label==2] = 255
b[label==3] = 255
b[label==4] = 125
g[label==4] = 125largeimage = cv2.resize(image,(size[1],size[0]),cv2.INTER_LINEAR)
cv2.imshow("lab",largeimage)
cv2.waitKey(0)

4 效果