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 -1
    if inputidx[1] >= w:
        return -1
    if inputidx[0] < 0:
        return -1
    if inputidx[1] < 0:
        return -1
    
    return 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#由于累积直方图是作为下标使用，因此最终求和结果不可等于imsize
    histsum = 0
    for i in range(0,256):
        histsum = histsum + hist[i]
        cnt[i+1] = histsum
    
    cnt = 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,x
 
            cntidx[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:
                    continue
                
                
                if labelout[neighboraxis[0],neighboraxis[1]] >= 0:#周围有已经计算过的

                    putquemask[neighboraxis[0],neighboraxis[1]] > 0#标记已经入队过了
                    que.put(nowpixaxis)#入队
                    break
            

        #开始遍历队列
        while(True):
            if que.empty():
                break
            
            nowpixaxis = 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:
                    continue
                labelnow = labelout[nowpixaxis[0],nowpixaxis[1]]
                labelneighbor = int(labelout[neighboraxis[0],neighboraxis[1]])
                
                if labelnow == -2 and labelneighbor > 0:
                    labelout[nowpixaxis[0],nowpixaxis[1]] = labelneighbor
                    labelsize[labelneighbor] = labelsize[labelneighbor] + 1
                    
                elif labelnow > 0 and labelneighbor > 0 and labelnow != labelneighbor:
                    labelout[nowpixaxis[0],nowpixaxis[1]] = 0
                    
                if 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]] = currenLabel
                
                que.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:
                            continue
                        if  putquemask[neighboraxis[0],neighboraxis[1]] == 0 and labelout[neighboraxis[0],neighboraxis[1]] == -2:
                          
                            labelout[neighboraxis[0],neighboraxis[1]] = currenLabel
                            que.put(neighboraxis)
                            putquemask[neighboraxis[0],neighboraxis[1]] = 1
                            labelsize[currenLabel] = labelsize[currenLabel] + 1

    return 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] = 125

largeimage = cv2.resize(image,(size[1],size[0]),cv2.INTER_LINEAR)
cv2.imshow("lab",largeimage)
cv2.waitKey(0)
    

4 效果