分水岭算法的理解和应用

分水岭算法

• 主要思想
  • 图像的灰度空间很像地球表面的整个地理结构，每个像素的灰度值代表高度。分水岭就是灰度值较大的像素连成的线。二值化阈值可以理解为水平面，比灰度二值化阈值小的像素区域会被淹没。
  • 随着水位线的升高，被淹没的区域越来越多，通过“修建大坝”防止两个区域的水汇集。最终整张图像分成两个部分：一是水坝线（分水岭线像素集），二是水淹区域（像素集）。

1 预备知识

1.1 测地距离

• 首先，当根据经纬度计算地球上两点间的距离时，使用的是球面距离，而非欧氏距离
• 相似的，分水岭算法中用到的测地距离，是三维的，即[x,y,h]
• 测地距离计算方法为：两点间沿着三维表面的曲面走的最短路径

1.2 固定图像二值化：cv2.threshold函数解析

• retval, dst = cv2.threshold(src, thresh, maxval, type[, dst])
  • retval, dst分别表示二值化灰度阈值和二值化后的图像数据
  • scr:图像数据
  • thresh:自定义阈值
  • maxval:最大值,像素超过了阈值，所赋予的值
  • type:可选择的算法类型
• 例子
  • 超过阈值的设为maxval,其他为0：ret,thresh1 = cv2.threshold(img,127,255,cv2.THRESH_BINARY)
  • 超过阈值的设为0，其他未maxval：ret,thresh2 = cv2.threshold(img,127,255,cv2.THRESH_BINARY_INV)
  • 大于阈值部分设为阈值,其他不变:ret,thresh3 = cv2.threshold(img,127,255,cv2.THRESH_TRUNC)
  • 低于阈值的设置为0,其他保持不变：ret,thresh4 = cv2.threshold(img,127,255,cv2.THRESH_TOZERO)
  • 低于阈值的设置为maxval,其他保持不变ret,thresh5 = cv2.threshold(img,127,255,cv2.THRESH_TOZERO_INV)

1.3 自适应二值化：cv2.adaptiveThreshold

• dst = cv2.adaptiveThreshold(src, maxval, thresh_type, type, Block Size, C)
  • thresh_type:
    • cv2.ADAPTIVE_THRESH_MEAN_C：计算均值时每个像素的权值是相等的
    • cv2.ADAPTIVE_THRESH_GAUSSIAN_C：计算均值时每个像素的权值根据其到中心点的距离通过高斯方程得到
  • type:计算均值后，即生成阈值，而后采取固定图像二值化方法
  • block size:图片分块的大小
  • C:阈值计算中的所加减的常数项

1.4 图像形态学

• 原理介绍：需要在二值化后的图像中进行。相当于对图像进行池化操作，膨胀相当于最大池，图像的高亮区域或白色部分进行扩张,其运行结果图比原图的高亮区域更大；腐蚀相当于最小池，将图像中的高亮区域或白色部分进行缩减细化,其运行结果图比原图的高亮区域更小。
• 卷积核锚点（中心点）值被最亮/暗像素值替换

1. 图像膨胀

• dst = cv2.dilate(src, kernel, iterations)
  • kernel:卷积核，n*n维矩阵
  • interations:迭代次数
  • 作用：试图找到连通分支（具有相似颜色或强度的像素点的大块互相分离的区域）。因为多数情况下，一个大区域可能被噪声、阴影等分割成多个部分

2. 图像腐蚀

• dst = cv2.erode(src, kernel, iterations)
  • 作用：消除图像中“斑点”噪声，且能确保图像内较大区域依然存在

3. 创建核的方法

• 使用numpy构建结构化元素，都是正方形的核。
  • np.ones((5,5), np.unit8)
• 如果要构建圆形或者椭圆形的核，可以使用cv2.getStructuringElement()函数
  • kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (3,3))
  • MORPH_RECT:矩形
  • MORPH_CROSS:交叉十字
  • MORPH_ELLIPSE:椭圆
  • (n,n)：代表卷积核尺寸

4. 更加通用的形态学

• opening = cv2.morphologyEx(data, type, kernel, iterations)

• 开运算

  • 先腐蚀再膨胀：向上的孤立点被消除，去除小的明亮区域，并且剩余的明亮区域被隔绝，但大小不变
  • type:cv2.MORPH_OPEN

• 闭运算

  • 先膨胀后腐蚀：消除低亮度值的孤立点，亮的区域合在一起，但他们的基本大小不变
  • type:cv2.MORPH_CLOSE

• 形态学梯度

  • 膨胀后的图像-腐蚀后的图像
  • 可以突出边缘
  • type:cv2.MORPH_GRADIENT

• 礼帽和黑帽

  • 礼帽：原图-开运算后的图像
    • 作用：突出比原图中周围区域更明亮的部分，与核的大小相关
  • 黑帽：原图-闭运算后的图像
    • 作用：突出比原图中周围区域更暗的部分，与核的大小相关

2 算法流程

• 把梯度图像中的所有像素按照灰度值分类，并设定一个测地距离阈值。
• 找到灰度值最小的像素点（默认标记为灰度值最低点），让threshold从最小值开始增长，这些点为起始点。
• 水平面在增长的过程中，会碰到周围的邻域像素，测量这些像素到起始点（灰度值最低点）的测地距离，如果小于设定阈值，则将这些像素淹没，否则在这些像素上设置大坝，这样就对这些邻域像素进行了分类。
• 思考：这个阈值应该如何设定？若设定偏低，则会出现较多杂乱的边缘线，若设定偏高，则往往“防洪”效果不好。

3 opencv实例操作

3.1 导入相关库并读取数据

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

path = r''
origin = cv2.imread(path)

3.2 图像预处理

(一) 转化为灰度图片，并进行二值化处理
gray = cv2.cvtColor(origin, cv2.COLOR_BGR2GRAY)
ret, binary = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY|cv2.THRESH_OTSU)

(二) 去除杂点：三次膨胀+两次腐蚀+一次膨胀
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (2,2))
opening1 = cv2.dilate(binary, kernel, iterations = 3)
opening2 = cv2.erode(opening1, kernel, iterations = 2)
opening3 = cv2.dilate(opening2, kernel, iterations=1)

3.3 确定前景和背景和未知区域

• 下面是我自己的理解啦
• 前景：确定属于目标的部分
• 背景：一定包含全部目标的部分，背景包含前景
• 未知区域：背景-前景，因此包含边缘线

(一) 确定前景
sure_bg = cv2.dilate(opening3, kernel, iterations = 3) # 通过膨胀获取背景区域

(二)确定背景
dist = cv2.distanceTransform(opening3, cv2.DIST_L2, 5) # 通过距离变换获取前景区域
ret, sure_fg = cv2.threshold(dist, dist.max()*0.45, 255, cv2.THRESH_BINARY) # 0.45根据实际情况来调

(三)确定未知区域
surface_fg = np.uint8(sure_fg)  # 背景空间为整型空间，前景为浮点型空间，这里做一下统一
unknown = cv2.subtract(sure_bg, surface_fg)
plt.imshow(unknown, cmap = plt.cm.gray)

3.4 创建markers并实现watershed

• cv2.watershed函数需要将包含边缘的未知区域设置为0

ret, markers = cv2.connectedComponents(surface_fg) # 所有的前景区域非0
markers = markers + 1 # 此时所有的区域非零，其中包含边缘的未知区域和非背景区域为1
markers[unknown == 255] = 0  # 令包含边缘的位置区域为0

# 注:watershed会将找到的栅栏在markers中设置为-1

markers = cv2.watershed(origin, markers = markers) # 获取栅栏
origin[markers == -1] = [0,0,255]  # 根据栅栏，我们对原图像进行操作，对栅栏区域设置为红色

• 分水岭算法的简单应用就先到这里啦，想要得到更精细的“水坝线”，还需要自己多多调试
• 欢迎邮件交流：hu.yf@outlook.com