霍夫变换和分水岭图像分割

一、霍夫变换

1. 直线检测

对于一幅图像中的直线信息，人们的肉眼很容易就将其提取出来，但是对于计算机而言，这是难以实现的。因此我们需要对图像进行一些操作后提取出直线信息。前面提到的边缘检测或者阈值操作会因为噪声的存在而导致最后得到的结果存在一定的误差，比如说误判、错判和断线等。所以我们要使用基于边缘检测的结果使用霍夫变换来提取图像所有的直线信息。

原图和结果图如下：

霍夫线检测结果.PNG

边缘检测结果：

霍夫变换直线表示：

具体表示方法： $\begin{array}{l} y=\left(-\frac{\cos \theta}{\sin \theta}\right) x+\left(\frac{r}{\sin \theta}\right) , namely, y = kx + b\\ r=x \cos \theta+y \sin \theta \end{array}$
注意：
1. $\rho$ 和 $\theta$ 并不是极坐标表示方法，首先对于极坐标表示， $\rho$ 值会不断变换，而霍夫变换中的直线表示 $\rho$ 则不会变化。
2. 且 $\rho$ 是原点到直线的最短距离(既然最短则应为不变)
3. 一个 $\rho$ 和 $\theta$ 就能确定一条直线
直线检测的基本流程
1. 读入图像
2. 转化为灰度图，便于后续的阈值处理或边缘检测
3. 霍夫变换检测直线
  1. 设定 $\rho$ 的取值范围，对于一副图像而言，图中直线最长只能达到对角线长度。分为有限个取值区间。
  2. 将 $\theta$ 的取值范围分为若干个有限区间。
  3. 假设上面的 $\rho$ 和 $\theta$ 的区间个数分别为 8， 4。如下所示，最终会形成一个 8 * 4 的阵列。
  4. 遍历原图，遇到非边缘点则跳过，否则进入计算：将该点坐标 $(x_i,y_i)$ 代入 $r = xcos\theta + y sin\theta$ ，并将 $\theta$ 区间的所有值依次代入计算，会得到90个 r值，将其分别填入到 8 个区间之内。直到原图都被遍历完毕(想象为投票过程)。
  5. 得票数越多，则对应的 $\rho$ 和 $\theta$ 就越有可能确定一条直线。如下图所示：

对于cv.HoughLines得到的直线集，我们要从其每个元素的第一个元素的也即是 line in lines， line[0]中提取 $\rho$ 和 $\theta$ 信息。进而得到直线上距离远点最近的点的坐标 $x_0 = \rho cos\theta , y_0 = \rho sin\theta$ ，由于cv.line需要使用两个点才能绘制出一条直线，那么我们要根据 $(x_0,y_0)$ 和已知的 $\theta$ 来找出两个点确定此条直线。利用几何知识，我们可以选取一个足够大的沿直线的偏移量 bias，使 $x_1 = x_0 + bias * sin \theta, x_2 = x_0 - bias * sin \theta$ 以及 $y_1 = y_0 + bias * cos \theta, y_2 = y_0 - bias * cos \theta$ 。由此可以确定最终结果。，要注意霍夫变换的坐标轴表示方法：

实现代码：

#coding:utf-8

import numpy as np
import matplotlib.pyplot as plt
import os
import cv2

def lines_detector_hough(edge,ThetaDim = None,DistStep = None,threshold = None,halfThetaWindowSize = 2,halfDistWindowSize = None):
    '''
    :param edge: 经过边缘检测得到的二值图
    :param ThetaDim: hough空间中theta轴的刻度数量(将[0,pi)均分为多少份),反应theta轴的粒度,越大粒度越细
    :param DistStep: hough空间中dist轴的划分粒度,即dist轴的最小单位长度
    :param threshold: 投票表决认定存在直线的起始阈值
    :return: 返回检测出的所有直线的参数(theta,dist)
    @author: bilibili-会飞的吴克
    '''
    imgsize = edge.shape
    if ThetaDim == None:
        ThetaDim = 90
    if DistStep == None:
        DistStep = 1
    MaxDist = np.sqrt(imgsize[0]**2 + imgsize[1]**2)
    DistDim = int(np.ceil(MaxDist/DistStep))

    if halfDistWindowSize == None:
        halfDistWindowSize = int(DistDim/50)
    accumulator = np.zeros((ThetaDim,DistDim)) # theta的范围是[0,pi). 在这里将[0,pi)进行了线性映射.类似的,也对Dist轴进行了线性映射

    sinTheta = [np.sin(t*np.pi/ThetaDim) for t in range(ThetaDim)]
    cosTheta = [np.cos(t*np.pi/ThetaDim) for t in range(ThetaDim)]

    for i in range(imgsize[0]):
        for j in range(imgsize[1]):
            if not edge[i,j] == 0:
                for k in range(ThetaDim):
                    accumulator[k][int(round((i*cosTheta[k]+j*sinTheta[k])*DistDim/MaxDist))] += 1

    M = accumulator.max()

    if threshold == None:
        threshold = int(M*2.3875/10)
    result = np.array(np.where(accumulator > threshold)) # 阈值化
    temp = [[],[]]
    for i in range(result.shape[1]):
        eight_neiborhood = accumulator[max(0, result[0,i] - halfThetaWindowSize + 1):min(result[0,i] + halfThetaWindowSize, accumulator.shape[0]), max(0, result[1,i] - halfDistWindowSize + 1):min(result[1,i] + halfDistWindowSize, accumulator.shape[1])]
        if (accumulator[result[0,i],result[1,i]] >= eight_neiborhood).all():
            temp[0].append(result[0,i])
            temp[1].append(result[1,i])

    result = np.array(temp)    # 非极大值抑制

    result = result.astype(np.float64)
    result[0] = result[0]*np.pi/ThetaDim
    result[1] = result[1]*MaxDist/DistDim

    return result

def drawLines(lines,edge,color = (255,0,0),err = 3):
    if len(edge.shape) == 2:
        result = np.dstack((edge,edge,edge))
    else:
        result = edge
    Cos = np.cos(lines[0])
    Sin = np.sin(lines[0])

    for i in range(edge.shape[0]):
        for j in range(edge.shape[1]):
            e = np.abs(lines[1] - i*Cos - j*Sin)
            if (e < err).any():
                result[i,j] = color

    return result


if __name__=='__main__':
    pic_path = './HoughImg/'
    pics = os.listdir(pic_path)

    for i in pics:
        if i[-5:] == '.jpeg' or i[-4:] == '.jpg':
            img = plt.imread(pic_path+i)
            blurred = cv2.GaussianBlur(img, (3, 3), 0)
            plt.imshow(blurred,cmap='gray')
            plt.axis('off')
            plt.show()

            if not len(blurred.shape) == 2:
                gray = cv2.cvtColor(blurred, cv2.COLOR_RGB2GRAY)
            else:
                gray = blurred
            edge = cv2.Canny(gray, 50, 150)   #  二值图 (0 或 255) 得到 canny边缘检测的结果


            lines = lines_detector_hough(edge)
            final_img = drawLines(lines,blurred)
            

            plt.imshow(final_img,cmap='gray')
            plt.axis('off')
            plt.show()

2. 圆形检测

霍夫变换圆表示方法

与直线相比(只需要两个参数 $\rho$ 和 $\theta$ )，霍夫变换表示圆需要用到三个参数：圆心： $(x_{center}, y_{center})$ ，半径： $radius$ 。那么我们的投票阵列就是 3 维矩阵，每一个 $(x_i,y_i,z_i)$ 处记载了圆心为 $(x_i,y_i)$ ，半径为 $z_i$ 的圆的投票数。同属一个圆的边缘点都会为此圆进行投票。，如下所示：

注意：若是按照之前检测直线的方法来对圆形进行检测，那么我们最后的时间复杂度将会为 $O(n^4)$ 左右，一旦图像像素点比较多，那么所需时间是非常大的，所以我们就采用霍夫梯度法来对圆进行检测。
霍夫梯度法的基本思想

一个正圆周上面的所有点的梯度方向会交于一点，这个点就是圆心。确定圆心之后我们可以使用圆心获得半径的大小。而对于圆周边缘点的梯度方向，我们可以使用Sobel算子进行确定。
圆形检测的基本流程
1. 确定圆心：遍历图像，遇到边缘点，则使用Sobel算子计算出对 x 和 y 的偏导数，然后得出该点的梯度方向，然后可以利用这个方向和该定点(x,y)确定的直线将直线所在的所有像素点都投一票。最终票数最多的点一定是圆心。
2. 确定半径：利用筛选出来的圆心，所有可能的半径再次形成投票阵列。再次遍历所有边缘点，计算其到圆心的距离，然后根据距离为对应的半径长度进行投票。如下所示：

实现代码：

#coding:utf-8

import numpy as np
import matplotlib.pyplot as plt
import os
import cv2


def convlove(filter,mat,padding,strides):
    '''
    :param filter:卷积核，必须为二维(2 x 1也算二维) 否则返回None
    :param mat:图片
    :param padding:对齐
    :param strides:移动步长
    :return:返回卷积后的图片。(灰度图，彩图都适用)
    @author:bilibili-会飞的吴克
    '''
    result = None
    filter_size = filter.shape
    mat_size = mat.shape
    if len(filter_size) == 2:
        if len(mat_size) == 3:
            channel = []
            for i in range(mat_size[-1]):
                pad_mat = np.pad(mat[:,:,i], ((padding[0], padding[1]), (padding[2], padding[3])), 'constant')
                temp = []
                for j in range(0,mat_size[0],strides[1]):
                    temp.append([])
                    for k in range(0,mat_size[1],strides[0]):
                        val = (filter*pad_mat[j*strides[1]:j*strides[1]+filter_size[0],
                                      k*strides[0]:k*strides[0]+filter_size[1]]).sum()
                        temp[-1].append(val)
                channel.append(np.array(temp))

            channel = tuple(channel)
            result = np.dstack(channel)
        elif len(mat_size) == 2:
            channel = []
            pad_mat = np.pad(mat, ((padding[0], padding[1]), (padding[2], padding[3])), 'constant')
            for j in range(0, mat_size[0], strides[1]):
                channel.append([])
                for k in range(0, mat_size[1], strides[0]):
                    val = (filter * pad_mat[j * strides[1]:j * strides[1] + filter_size[0],
                                    k * strides[0]:k * strides[0] + filter_size[1]]).sum()
                    channel[-1].append(val)


            result = np.array(channel)


    return result


def AHTforCircles(edge,center_threhold_factor = None,score_threhold = None,min_center_dist = None,minRad = None,maxRad = None,center_axis_scale = None,radius_scale = None,halfWindow = None,max_circle_num = None):
    if center_threhold_factor == None:
        center_threhold_factor = 10.0
    if score_threhold == None:
        score_threhold = 15.0
    if min_center_dist == None:
        min_center_dist = 80.0
    if minRad == None:
        minRad = 0.0
    if maxRad == None:
        maxRad = 1e7*1.0
    if center_axis_scale == None:
        center_axis_scale = 1.0
    if radius_scale == None:
        radius_scale = 1.0
    if halfWindow == None:
        halfWindow = 2
    if max_circle_num == None:
        max_circle_num = 6
    min_center_dist_square = min_center_dist**2


    sobel_kernel_y = np.array([[-1.0, -2.0, -1.0], [0.0, 0.0, 0.0], [1.0, 2.0, 1.0]])
    sobel_kernel_x = np.array([[-1.0, 0.0, 1.0], [-2.0, 0.0, 2.0], [-1.0, 0.0, 1.0]])
    edge_x = convlove(sobel_kernel_x,edge,[1,1,1,1],[1,1])
    edge_y = convlove(sobel_kernel_y,edge,[1,1,1,1],[1,1])


    center_accumulator = np.zeros((int(np.ceil(center_axis_scale*edge.shape[0])),int(np.ceil(center_axis_scale*edge.shape[1]))))
    k = np.array([[r for c in range(center_accumulator.shape[1])] for r in range(center_accumulator.shape[0])])
    l = np.array([[c for c in range(center_accumulator.shape[1])] for r in range(center_accumulator.shape[0])])
    minRad_square = minRad**2
    maxRad_square = maxRad**2
    points = [[],[]]

    edge_x_pad = np.pad(edge_x,((1,1),(1,1)),'constant')
    edge_y_pad = np.pad(edge_y,((1,1),(1,1)),'constant')
    Gaussian_filter_3 = 1.0 / 16 * np.array([(1.0, 2.0, 1.0), (2.0, 4.0, 2.0), (1.0, 2.0, 1.0)])

    for i in range(edge.shape[0]):
        for j in range(edge.shape[1]):
            if not edge[i,j] == 0:
                dx_neibor = edge_x_pad[i:i+3,j:j+3]
                dy_neibor = edge_y_pad[i:i+3,j:j+3]
                dx = (dx_neibor*Gaussian_filter_3).sum()
                dy = (dy_neibor*Gaussian_filter_3).sum()
                if not (dx == 0 and dy == 0):
                    t1 = (k/center_axis_scale-i)
                    t2 = (l/center_axis_scale-j)
                    t3 = t1**2 + t2**2
                    temp = (t3 > minRad_square)&(t3 < maxRad_square)&(np.abs(dx*t1-dy*t2) < 1e-4)
                    center_accumulator[temp] += 1
                    points[0].append(i)
                    points[1].append(j)

    M = center_accumulator.mean()
    for i in range(center_accumulator.shape[0]):
        for j in range(center_accumulator.shape[1]):
            neibor = \
                center_accumulator[max(0, i - halfWindow + 1):min(i + halfWindow, center_accumulator.shape[0]),
                max(0, j - halfWindow + 1):min(j + halfWindow, center_accumulator.shape[1])]
            if not (center_accumulator[i,j] >= neibor).all():
                center_accumulator[i,j] = 0
                                                                        # 非极大值抑制

    plt.imshow(center_accumulator,cmap='gray')
    plt.axis('off')
    plt.show()

    center_threshold = M * center_threhold_factor
    possible_centers = np.array(np.where(center_accumulator > center_threshold))  # 阈值化


    sort_centers = []
    for i in range(possible_centers.shape[1]):
        sort_centers.append([])
        sort_centers[-1].append(possible_centers[0,i])
        sort_centers[-1].append(possible_centers[1, i])
        sort_centers[-1].append(center_accumulator[sort_centers[-1][0],sort_centers[-1][1]])

    sort_centers.sort(key=lambda x:x[2],reverse=True)

    centers = [[],[],[]]
    points = np.array(points)
    for i in range(len(sort_centers)):
        radius_accumulator = np.zeros(
            (int(np.ceil(radius_scale * min(maxRad, np.sqrt(edge.shape[0] ** 2 + edge.shape[1] ** 2)) + 1))),dtype=np.float32)
        if not len(centers[0]) < max_circle_num:
            break
        iscenter = True
        for j in range(len(centers[0])):
            d1 = sort_centers[i][0]/center_axis_scale - centers[0][j]
            d2 = sort_centers[i][1]/center_axis_scale - centers[1][j]
            if d1**2 + d2**2 < min_center_dist_square:
                iscenter = False
                break

        if not iscenter:
            continue

        temp = np.sqrt((points[0,:] - sort_centers[i][0] / center_axis_scale) ** 2 + (points[1,:] - sort_centers[i][1] / center_axis_scale) ** 2)
        temp2 = (temp > minRad) & (temp < maxRad)
        temp = (np.round(radius_scale * temp)).astype(np.int32)
        for j in range(temp.shape[0]):
            if temp2[j]:
                radius_accumulator[temp[j]] += 1
        for j in range(radius_accumulator.shape[0]):
            if j == 0 or j == 1:
                continue
            if not radius_accumulator[j] == 0:
                radius_accumulator[j] = radius_accumulator[j]*radius_scale/np.log(j) #radius_accumulator[j]*radius_scale/j
        score_i = radius_accumulator.argmax(axis=-1)
        if radius_accumulator[score_i] < score_threhold:
            iscenter = False

        if iscenter:
            centers[0].append(sort_centers[i][0]/center_axis_scale)
            centers[1].append(sort_centers[i][1]/center_axis_scale)
            centers[2].append(score_i/radius_scale)


    centers = np.array(centers)
    centers = centers.astype(np.float64)

    return centers








def drawCircles(circles,edge,color = (0,0,255),err = 600):
    if len(edge.shape) == 2:
        result = np.dstack((edge,edge,edge))
    else:
        result = edge
    for i in range(edge.shape[0]):
        for j in range(edge.shape[1]):
            dist_square = (circles[0]-i)**2 + (circles[1]-j)**2
            e = np.abs(circles[2]**2 - dist_square)
            if (e < err).any():
                result[i,j] = color
            if (dist_square < 25.0).any():
                result[i,j] = (255,0,0)
    return result


if __name__=='__main__':
    pic_path = './HoughCircleImg/'
    pics = os.listdir(pic_path)


    params = {
        '1.jpeg':{
        'center_threhold_factor': 3.33,
        'score_threhold':15.0,
         'min_center_dist':80.0,
         'minRad':0.0,
         'maxRad':1e7*1.0,
         'center_axis_scale':1.0,
         'radius_scale':1.0,
         'halfWindow':2,
        'max_circle_num':1
         },
        '4.jpg':{
        'center_threhold_factor': 2.0,
        'score_threhold': 15.0,
         'min_center_dist': 80.0,
         'minRad': 0.0,
         'maxRad': 1e7 * 1.0,
         'center_axis_scale': 1.0,
         'radius_scale': 1.0,
         'halfWindow': 2,
         'max_circle_num': 6
         },
        '2.jpeg':{
        'center_threhold_factor': 3.33,
        'score_threhold': 50.0,
         'min_center_dist': 80.0,
         'minRad': 0.0,
         'maxRad': 1e7 * 1.0,
         'center_axis_scale': 0.9,
         'radius_scale': 1.0,
         'halfWindow': 2,
         'max_circle_num': 1
         },
        '3.jpeg':{
        'center_threhold_factor': 1.5,
        'score_threhold': 56.0,
         'min_center_dist': 80.0,
         'minRad': 0.0,
         'maxRad': 1e7 * 1.0,
         'center_axis_scale': 0.8,
         'radius_scale': 1.0,
         'halfWindow': 2,
         'max_circle_num': 1
         },
        '0.jpg':{
        'center_threhold_factor': 1.5,
        'score_threhold': 30.0,
         'min_center_dist': 80.0,
         'minRad': 0.0,
         'maxRad': 1e7 * 1.0,
         'center_axis_scale': 1.0,
         'radius_scale': 1.0,
         'halfWindow': 2,
         'max_circle_num': 1
         }
    }
    for i in pics:
        if i[-5:] == '.jpeg' or i[-4:] == '.jpg':
            img = plt.imread(pic_path+i)
            blurred = cv2.GaussianBlur(img, (3, 3), 0)
            plt.imshow(blurred)
            plt.axis('off')
            plt.show()

            if not len(blurred.shape) == 2:
                gray = cv2.cvtColor(blurred, cv2.COLOR_RGB2GRAY)
            else:
                gray = blurred
            edge = cv2.Canny(gray, 50, 150)   #  二值图 (0 或 255) 得到 canny边缘检测的结果

            circles = AHTforCircles(edge,
                                    center_threhold_factor=params[i]['center_threhold_factor'],score_threhold=params[i]['score_threhold'],min_center_dist=params[i]['min_center_dist'],minRad=params[i]['minRad'],
                                    maxRad=params[i]['maxRad'],center_axis_scale=params[i]['center_axis_scale'],radius_scale=params[i]['radius_scale'],
                                    halfWindow=params[i]['halfWindow'],max_circle_num=params[i]['max_circle_num'])
            final_img = drawCircles(circles,blurred)

            plt.imshow(final_img)
            plt.axis('off')
            plt.show()

结果：

可以看出误差还是很大。

二、分水岭算法

1. 概念

处理对象：分水岭算法针对的是灰度图。
下雨思想：将灰度图视为一张三维的地形图，俯视图自然就是灰度图本身，高度的话就是灰度值。灰度值较低的区域会形成山谷，灰度值较高的区域会形成山巅。给每一个山谷区域的“积水”涂上不同的颜色，然后我们不断的模拟“下雨”(水不会溢出到图像外)，那么这些山谷中的水会逐渐汇集，汇集的同时，这些颜色会变黑，也就是边缘会变黑形成屏障，同时双方的水就不会互相交融。最终，直到水位超过最高的山巅，得到的图像就会有各种颜色以及黑色边界。，如图：
洪泛思想：将灰度图的山谷下面打穿留一个洞，然后将这些地形图放入无穷深的大海中，那么海水就会从这些洞中涌出，接下来就和下雨一致了，如下图所示：

2. 算法流程

找出所有的极小值区域
被水淹没

三、分水岭算法具体说明(转自知乎@程序员阿德)

分水岭（Watershed）是基于地理形态的分析的图像分割算法，模仿地理结构（比如山川、沟壑，盆地）来实现对不同物体的分类。

分水岭算法中会用到一个重要的概念——测地线距离。

1. 测地线距离(Geodesic Distance)

测地线距离就是地球表面两点之间的最短路径（可执行路径）的距离，在图论中，Geodesic Distance 就是图中两节点的最短路径的距离，这与平时在几何空间通常用到的 Euclidean Distance（欧氏距离），即两点之间的最短距离有所区别。前者是沿着圆走，后者两点之间直线最短。

在下图中，两个黑点的 Euclidean Distance 是用虚线所表示的线段的长度 $d_{15}$ ，而 Geodesic Distance 作为实际路径的最短距离，其距离应为沿途实线段距离之和的最小值，即 $d_{12}+d_{23}+d_{34}+d_{45}$ 。

在三维曲面空间中两点间的测地距离就是两点间沿着三维曲面的表面走的最短路径。

2. 分水岭算法

图像的灰度空间很像地球表面的整个地理结构，每个像素的灰度值代表高度。其中的灰度值较大的像素连成的线可以看做山脊，也就是分水岭。其中的水就是用于二值化的gray threshold level，二值化阈值可以理解为水平面，比水平面低的区域会被淹没，刚开始用水填充每个孤立的山谷(局部最小值)。

当水平面上升到一定高度时，水就会溢出当前山谷，可以通过在分水岭上修大坝，从而避免两个山谷的水汇集，这样图像就被分成2个像素集，一个是被水淹没的山谷像素集，一个是分水岭线像素集。最终这些大坝形成的线就对整个图像进行了分区，实现对图像的分割。

在该算法中，空间上相邻并且灰度值相近的像素被划分为一个区域。

3. 分水岭算法的整个过程

把梯度图像中的所有像素按照灰度值进行分类，并设定一个测地距离阈值。
找到灰度值最小的像素点（默认标记为灰度值最低点），让threshold从最小值开始增长，这些点为起始点。
水平面在增长的过程中，会碰到周围的邻域像素，测量这些像素到起始点（灰度值最低点）的测地距离（沿着山坡走），如果小于设定阈值，则将这些像素淹没，否则在这些像素上设置大坝，这样就对这些邻域像素进行了分类。

随着水平面越来越高，会设置更多更高的大坝，直到灰度值的最大值，所有区域都在分水岭线上相遇，这些大坝就对整个图像像素的进行了分区。

整个过程可以查看下面这个动图：

用上面的算法对图像进行分水岭运算，由于噪声点或其它因素的干扰，可能会得到密密麻麻的小区域，即图像被分得太细（over-segmented，过度分割），这因为图像中有非常多的局部极小值点，每个点都会自成一个小区域。

其中的解决方法：

对图像进行高斯平滑操作，抹除很多小的最小值，这些小分区就会合并。
不从最小值开始增长，可以将相对较高的灰度值像素作为起始点（需要用户手动标记），从标记处开始进行淹没，则很多小区域都会被合并为一个区域，这被称为基于图像标记(mark)的分水岭算法。

下面三个图分别是原图，分水岭过分割的图以及基于标记的分水岭算法得到的图：

其中标记的每个点就相当于分水岭中的注水点，从这些点开始注水使得水平面上升，但是如上图所示，图像中需要分割的区域太多了，手动标记太麻烦，我们可使用距离转换的方法进行标记，OpenCV中就是使用的这种方法。

四、 OpenCV中分水岭算法

在OpenCV中，我们需要给不同区域贴上不同的标签。用大于1的整数表示我们确定为前景或对象的区域，用1表示我们确定为背景或非对象的区域，最后用0表示我们无法确定的区域。然后应用分水岭算法，我们的标记图像将被更新，更新后的标记图像的边界像素值为-1。

下面对相互接触的硬币应用距离变换和分水岭分割。

先使用 Otsu's 二值化对图像进行二值化。

import cv2
import numpy as np

img = cv2.imread('coins.png')
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
ret, thresh = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)

先使用开运算去除图像中的细小白色噪点(此处采用中值滤波也可以)，然后通过腐蚀运算移除边界像素，得到的图像中的白色区域肯定是真实前景，即靠近硬币中心的区域（下面左边的图）；膨胀运算使得一部分背景成为了物体的边界，得到的图像中的黑色区域肯定是真实背景，即远离硬币的区域（下面中间的图）。

剩下的区域（硬币的边界附近）还不能确定是前景还是背景。可通过膨胀图减去腐蚀图得到，下图中的白色部分为不确定区域（下面右边的图）。

# noise removal
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (3,3))
opening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel, iterations=2)

sure_bg = cv2.dilate(opening, kernel, iterations=2)  # sure background area
sure_fg = cv2.erode(opening, kernel, iterations=2)  # sure foreground area
unknown = cv2.subtract(sure_bg, sure_fg)  # unknown area

剩下的区域不确定是硬币还是背景，这些区域通常在前景和背景接触的区域(或者两个不同硬币接触的区域)，我们称之为边界。通过分水岭算法应该能找到确定的边界。

注意：为什么硬币之间相互接触就不能直接使用分水岭算法呢？

首先硬币之间相互接触会造成不确定区域的非联通，如果我们直接使用分水岭算法则会出现漏选情况，也就是分水岭算法只能基于未确定的区域进行操作，若是未确定区域没有囊括可能的所有边缘，那么分水岭算法也不会得出对应的边缘。

由于硬币之间彼此接触，我们使用另一个确定前景的方法，就是带阈值的距离变换。

下面左边的图为得到的距离转换图像，其中每个像素的值为其到最近的背景像素（灰度值为0）的距离，可以看到硬币的中心像素值最大最亮（中心离背景像素最远）。对其进行二值处理就得到了分离的前景图（下面中间的图），白色区域肯定是硬币区域，而且还相互分离，下面右边的图为之前的膨胀图减去中间这个表示前景的图。

# Perform the distance transform algorithm
# cv2.DIST_L2表明使用欧式距离计算
dist_transform = cv2.distanceTransform(opening, cv2.DIST_L2, 5)
# Normalize the distance image for range = {0.0, 1.0}
cv2.normalize(dist_transform, dist_transform, 0, 1.0, cv2.NORM_MINMAX)

# Finding sure foreground area
ret, sure_fg = cv2.threshold(dist_transform, 0.5*dist_transform.max(), 255, 0)

# Finding unknown region
sure_fg = np.uint8(sure_fg)
unknown = cv2.subtract(sure_bg,sure_fg)

现在我们可以确定哪些是硬币区域，哪些是背景区域。然后创建标记(marker，它是一个与原始图像大小相同的矩阵，int32数据类型)，表示其中的每个区域。分水岭算法将标记的0的区域视为不确定区域，将标记为1的区域视为背景区域，将标记大于1的正整数表示我们想得到的前景。

我们可以使用 cv2.connectedComponents() 来实现这个功能，它是用0标记图像的背景，用大于0的整数标记其他对象。所以我们需要对其进行加一，用1来标记图像的背景。

cv2.connectedComponents() 将传入图像中的白色区域视为组件（前景）。

# 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

注意：得到的markers矩阵的元素类型为 int32，要使用 imshow() 进行显示，需要将其转换为 uint8 类型（ markers=np.uint8(markers) ）。

我们对得到的markers进行显示：

markers_copy = markers.copy()
markers_copy[markers==0] = 150  # 灰色表示背景
markers_copy[markers==1] = 0    # 黑色表示背景
markers_copy[markers>1] = 255   # 白色表示前景

markers_copy = np.uint8(markers_copy)

标记图像已经完成了，最后应用分水岭算法。然后标记图像将被修改，边界区域将被标记为-1。

# 使用分水岭算法执行基于标记的图像分割，将图像中的对象与背景分离
markers = cv2.watershed(img, markers)
img[markers==-1] = [0,0,255]  # 将边界标记为红色

经过分水岭算法得到的新的标记图像和分割后的图像如下图所示：

任何两个相邻连接的组件不一定被分水岭边界(-1的像素)分开；例如在传递给 watershed 函数的初始标记图像中的物体相互接触。

1. 总结

我们通过一个例子介绍了分水岭算法的整个过程，主要分为以下几步：

对图进行灰度化和二值化得到二值图像
通过膨胀得到确定的背景区域，通过距离转换得到确定的前景区域，剩余部分为不确定区域
对确定的前景图像进行连接组件处理，得到标记图像
根据标记图像对原图像应用分水岭算法，更新标记图像

OpenCV Extra 01 - 分水岭算法与霍夫变换