? ? ? ? 从第九章图像形态学开始，研究方向便开始从输入和输出都是图像的图像处理方法，转变成了输入图像但输出为从这些图像中提取出来的属性的图像处理方法，这章中的分割是这一方向的另一个主要步骤。所谓分割，就是将图像分为构成它的子区域或物体，根据要解决的问题决定细分的程度，超过识别这些元素所需细节的分割是没有意义的。

? ? ? ? 在这个章节中的多数分割算法均基于灰度值的两个基本性质之一：不连续性和相似性。第一类中，方法是以灰度突变为基础分割一幅图像，比如图像的边缘等。第二类中，主要方法是根据一组预定义的准则将一幅图像分割为相似的区域。阈值处理、区域生长、区域分裂和区域聚合都是这类方法的例子。

一、基础知识?

? ? ? ? 假定令R表示一幅图像所占据的整个空间区域，由此就可以将图像分割为若干个子区域 $R_{1},R_{2},R_{3},R_{4},R_{5},R_{6}.....$ 且满足：

其中， $Q(R_{k})$ 是定义在集合 $R_{k}$ 的点上的一个逻辑属性?，如果 $R_{I}$ 和 $R_{j}$ 的并形成一个连通集，则称这两个区域是邻接的。

? ? ? ? 根据条件(a)我们知道，分割必须是完全的；即每个像素必须在一个区域内；条件(b)要求一个区域中的点以某些预定义的方式来连接(即这些点必须是4连接的或8连接的)；条件(c)指出，各个区域必须是不相交的。条件(d)涉及分割后的区域中的像素必须满足的属性；最后，条件(e)指出，两个邻接区域 $R_{I}$ 和 $R_{j}$ ,在属性Q的意义上必须是不同的。

? ? ? ? 所以分割中的基本问题就是把一幅图像分成满足前述条件的多个区域。下图说明了基于区域的分割。

从图中我们可以看出，图d与图a类似，但内部区域的灰度形成了一幅纹理模式，图e显示了计算该图像的边缘的结果?，灰度中大量的寄生变化使得识别原图像中的唯一边界很困难，因此基于边缘分割不是一种合适的方法，不过它们的外部区域都是恒定的，因此解决这个简单分割问题的全部就是区分纹理区域和恒定区域的一个属性。图c是将图像分割为两个区域后的结果，图e是基于区域特性的分割结果。

二、点、线和边缘检测

? ? ? ? 主要集中于以灰度局部剧烈变化检测为基础的分割方法上。涉及的三种图像特征是孤立点、线和边缘。

2.1 背景知识

? ? ? ? 正如利用局部平均平滑一幅图那样，假设平均处理类似于积分，对于灰度的突变，局部变化可以用微分来检测并不奇怪，由于变化非常短促，因此一阶微分和二阶微分适合于这种目的。

? ? ? ? 数字函数的导数可以用差分定义，对于一阶导数的任何近似的要求如下：

????????????????(1)在恒定灰度区域必须为零

????????????????(2)在灰度台阶或斜坡开始处必须不为零

????????????????(3)在沿灰度斜坡点处也必须不为零

????????类似地，对于所用的二阶导数的近似,要求如下：

? ? ? ? ? ? ? ? (1) 在恒定灰度区域必须为零

????????????????(2) 在灰度台阶或斜坡的开始处和结束处必须不为零

????????????????(3) 沿灰度斜坡必须为零

一维函数f(x)在x处的导数的近似：将函数f(x+ $\Delta x$ )展开为关于x的泰勒级数，令 $\Delta x$ =1，则：

二阶导数的表达式为：

二阶导数式中的第二行是关于点x+1的展开，变量减1得到点x的二阶导数：

????????一阶导数通常在图像中产生较粗的边缘。二阶导数对精细细节比如说细线、孤立点和噪声有较强的响应，在灰度斜坡和灰度台阶过滤处会产生双边缘响应，二阶导数的符号还可以用来确定边缘的过渡是从亮到暗还是从暗到亮。

? ? ? ? 空间滤波器用于计算图像中每个像素处的一阶导数和二阶导数，响应如下：

$z_{k}$ 是像素的灰度，该像素的空间位置对应于模板中第k个系数的位置。

3*3空间滤波器掩模?

2.2 孤立点检测

? ? ? ? 已知点的检测应以二阶导数为基础，则使用拉普拉斯：

?偏微分方程为：

将偏微分代入拉普拉斯公式可得：

????????使用上图a所示的拉普拉斯模板对图片进行检测，如果在某点处该模板的响应的绝对值超过了一个指定的阈值，那么就称在模板中心(x,y)处该点已经被检测到了。这时就可以用如下表达式产生一幅二值图像。

g是输出图像，T是一个非负的阈值。

? ? ? ? 给出的例图中，图b是一个带有通孔的涡轮叶片的x射线图像，通孔含有一个黑色像素。图c是是拉普拉斯模板与图像卷积的结果。图d则是用g(x,y)后得到的结果图。

2.3 线检测

? ? ? ? 使用二阶导数将会产生比一阶导数更细的线，因此对于线检测，也可以使用上图给出的拉普拉斯模板，不过二阶导数会产生双线效应，这需要进行适当的处理。

? ? ? ? 上图为一个利用拉普拉斯进行线检测的例子，图a是一幅二值图像，图b为它的拉普拉斯图像，在拉普拉斯图像中包含有负值，图b中中等灰度值表示为0，较暗的灰色表示负值，较亮的色调表示正值。

? ? ? ? 负值可以通过取拉普拉斯图像的绝对值进行简单处理，但是就如图c所示，这个方法会使线的宽度加倍。仅使用拉普拉斯的正值时，会产生如图d的图像，对比使用拉普拉斯的绝对值，正值产生的线更细，也更有用。

? ? ? ? 对于不同方向上的线检测所用的拉普拉斯模板也不同，特定方向的拉普拉斯模板如下：

2.4 边缘检测?

? ? ? ? 边缘检测是基于灰度突变来分割图像的最常用的方法。台阶边缘是指在1个像素的距离上发生两个灰度级间理想的过渡。下图a显示了一个垂直台阶边缘的一部分和通过该边缘的一个水平剖面。

? ? ? ? 我们所知道的数字图像都存在被模糊且带有噪声的边缘，其模糊的程度主要取决于聚焦机理的限制，噪声水平主要取决于成像系统的电子元件。这里边缘被建模成一个更接近斜坡的模型，如上图图b。

? ? ? ? 第三种模型为屋顶边缘，如上图c。屋顶边缘是通过一个区域的线的模型，屋顶边缘的基底(宽度)由该线的宽度和尖锐度决定。在极限情形下，当其基底为1个像素宽时，屋顶边缘只不过是- 条穿过图像中一个区域的一条1像素宽的线。

? ? ? ? 图10.10的图a是从图10.8的图b中的一段截取出来的一幅图像。图b显示了一条水平灰度剖面线。该图形还显示了灰度剖面的一阶导数和二阶导数。当沿着灰度剖面从左到右移动时，我们注意到，在斜坡的开始处和在斜坡上的各个点处，一阶导数为正。而在恒定灰度区域的一阶导数为零.在斜坡的开始处，二阶导数为正；在斜坡的结束处，二阶导数为负；在斜坡上的各点处，二阶导数为零;在恒定灰度区域的各点处，二阶导数为零。对于从亮到暗过渡的边缘，刚刚讨论的导数的符号则正好相反。零灰度轴和二阶导数极值间的连线的交点称为该二阶导数的零交叉点。

? ? ? ? 下面给出执行边缘检测的三个基本步骤：

????????????????1.为降噪对图像进行平滑处理。

????????????????2.边缘点的检测。

????????????????3.边缘定位。

2.5 基本边缘检测、更先进的边缘检测

? ? ? ? 检测灰度变化可以使用一阶或二阶导数来完成，为了在一幅图像f的(x,y)位置处寻找边缘的强度和方向，因此选择的工具就是梯度，用向量定义如下：

梯度指出了f在(x,y)处的最大变换率的方向。

其长度可以用它的模来表示，是梯度向量方向变化率的值，即：

这里的 $g_{x},g_{y}$ ?和 $M(x,y)$ 都是于原图像大小相同的图像，是x和y在f中所有像素位置上变化时产生的。通常将其称为梯度图像或梯度。

其方向由对x轴度量的角度给出：

梯度算子

? ? ? ? 先前提及的图像梯度中存在着偏导数的计算，由于我们所处理的均为数字量，因此就要求关于一点的邻域上的偏导数的数字近似，由此就可以得到：

以及

下图中的两个一维模板可以用于实现上面给出的两个公式。只需用提供的模板对f(x,y)进行滤波即可。?

? ? ? ? 对于对角线方向的边缘进行研究时，则需要一个二维模板进行处理，这里引出罗伯特交叉梯度算子，这是最早尝试使用具有对角优势的二维模板之一。

罗伯特算子以求对角像素之差为基础：

这两个公司可以由下面给出的模板的b和c对图像进行滤波来实现。?

????????一个2*2的模板在概念上很简单，但它们对于用关于中心点对称的模板来计算边缘方向不是很有用。大小为3*3的模板近似偏导数的最简单的数字近似为：

这个两个模板又称为prewitt算子。对应于下图的d、e。

当在中心系数上使用一个权值为2时，即：

模块就可以被称为Sobel算子。?对应于下图的f、g。

? ? ? ? prewitt模板实现要比Sobel模板更为简单一些，Sobel算子能够较好的抑制噪声。

高斯拉普拉斯算子(LoG)及Marr-Hildreth检测器

? ? ? ? 研究建议，用于边缘检验的算子应有两个显著的特点。第一个和最重要的特点是它应该是一个能计算图像中每一点处的一阶导数或二阶导数的数字近似的微分算子。第二，它应能被“调整”以便在任何期望的尺寸上起作用，因此，大的算子也可用于检测模糊边缘,小的算子可用于检测锐度集中的精细细节。

? ? ? ? 在众多算子中，滤波器 $\bigtriangledown ^{2}G$ 的表现较好， $\bigtriangledown ^{2}$ 是拉普拉斯算子，G为标准差为 $\sigma$ 的二维高斯函数：

$\bigtriangledown ^{2}G(x,y)$ 的表达式为：

? ? ? ? 上图为LoG的负函数的三维图、图像、剖面图和一个大小为5*5的模板。模板内的系数之和必须为，以保证模板的响应在恒定灰度区域为零。?

? ? ? ? 选取算子 $\bigtriangledown ^{2}G$ 的基础有两个基本概念，首先，算子的高斯部分会模糊图像，从而在尺寸上将结构的灰度降低到远小于 $\sigma$ 的程度。同时在滤波器的二阶导数方面，尽管一阶导数可用于检测灰度突变，但它们也是有方向的算子。其次，拉普拉斯有各向同性的优点，这不仅符合人的视觉系统特征而且对任何模板方向的灰度变化有相等的响应。

? ? ? ? Marr-Hildreth算法可以由LoG滤波器与f(x,y)卷积得到：

canny边缘检测器

? ? ? ? canny检测器是迄今为止边缘检测器中最为优秀的一种。基于三个基本目标：

1、低错误率.所有边缘都应被找到，并且应该没有伪响应.也就是检测到的边缘必须尽可能是真实的边缘。

2、边缘点应被很好地定位。已定位边缘必须尽可能接近真实边缘。也就是由检测器标记为边缘的点和真实边缘的中心之间的距离应该最小。

3、单一的边缘点响应。对于真实的边缘点，检测器仅应返回一个点.也就是真实边缘周围的局部最大数应该是最小的。这意味着在仅存一个单一边缘点的位置，检测器不应指出多个边缘像素。

?canny工作的本质是从数学上表达前面的三个准则，并尝试找出最佳解。对于一个由加性高斯白噪声污染的一维台阶边缘而言，最佳的台阶边缘检测器的一个较好近似是：

推广至二维时，令f(x.y)为输入图像，G(x,y)为高斯函数可以得到：

G与f的卷积为：

其梯度的幅值图像和角度图像计算为：

随后需对其赋值图像应用非最大抑制，并通过双阈值处理和连接分析检测并连接边缘。

利用python实现robert、prewitt、sobel、LoG、canny算子?

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

# 边缘检测的梯度算子 (Roberts 算子, Prewitt 算子, Sobel 算子,LoG算子，canny 算子)
img = cv2.imread("D:\\picture\\test2.jpg", 0)  

# Roberts 边缘算子
kernel_Roberts_x = np.array([[1, 0], [0, -1]])
kernel_Roberts_y = np.array([[0, -1], [1, 0]])
# Prewitt 边缘算子
kernel_Prewitt_x = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]])
kernel_Prewitt_y = np.array([[1, 1, 1], [0, 0, 0], [-1, -1, -1]])
# Sobel 边缘算子
kernel_Sobel_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
kernel_Sobel_y = np.array([[1, 2, 1], [0, 0, 0], [-1, -2, -1]])
# LoG算子
def LoG(image):
    kernel = np.array([[0, 0, 1, 0, 0], [0, 1, 2, 1, 0], [1, 2, -16, 2, 1], [0, 1, 2, 1, 0], [0, 0, 1, 0, 0]], dtype=int)
    img = cv2.filter2D(image, cv2.CV_16S, kernel)
    return cv2.convertScaleAbs(img)
# Canny算子
kSize = (5, 5)
imgGauss1 = cv2.GaussianBlur(img, kSize, sigmaX=1.0)
t1, t2 = 50, 150
Roberts_x = cv2.filter2D(img, -1, kernel_Roberts_x)
Roberts_y = cv2.filter2D(img, -1, kernel_Roberts_y)
Roberts = np.uint8(cv2.normalize(abs(Roberts_x) + abs(Roberts_y), None, 0, 255, cv2.NORM_MINMAX))
Prewitt_x = cv2.filter2D(img, -1, kernel_Prewitt_x)
Prewitt_y = cv2.filter2D(img, -1, kernel_Prewitt_y)
Prewitt = np.uint8(cv2.normalize(abs(Prewitt_x) + abs(Prewitt_y), None, 0, 255, cv2.NORM_MINMAX))
Sobel_x = cv2.filter2D(img, -1, kernel_Sobel_x)
Sobel_y = cv2.filter2D(img, -1, kernel_Sobel_y)
Sobel = np.uint8(cv2.normalize(abs(Sobel_x) + abs(Sobel_y), None, 0, 255, cv2.NORM_MINMAX))
LoG_img = LoG(img)
Canny = cv2.Canny(imgGauss1, t1, t2)
plt.figure(figsize=(12, 8))
plt.subplot(231), plt.title('Origin'), plt.imshow(img, cmap='gray'), plt.axis('off')
plt.subplot(232), plt.title('Roberts'), plt.imshow(Roberts, cmap='gray'), plt.axis('off')
plt.subplot(233), plt.title('Prewitt'), plt.imshow(Prewitt, cmap='gray'), plt.axis('off')
plt.subplot(234), plt.title('Sobel'), plt.imshow(Sobel, cmap='gray'), plt.axis('off')
plt.subplot(235), plt.title('Log'), plt.imshow(LoG_img, cmap='gray'), plt.axis('off')
plt.subplot(236), plt.title('Canny'), plt.imshow(Canny, cmap='gray'), plt.axis('off')
plt.tight_layout()
plt.show()

三、阈值处理

? ? ? ? 在先前的讨论中，采用的是寻找边缘线段，并通过将这些线段连接为边界的方法来识别区域。这里将从灰度值以及灰度值的特性把图像直接划分为区域进行讨论。

3.1 基础知识?

? ? ? ? 将上方的灰度直方图对应位图像f(x,y)，图像由暗色背景上的较亮物体组成，物体的像素和背景像素所具有的灰度值组合成了两种支配模式。分割这两种模式的阈值设为T，则分割后的图像就可以由下面的公式给出：

????????当T适用于整个图像时，上式给出的出来就称为全局阈值处理。当T值在一幅图像上改变时，称为可变阈值处理。局部阈值处理或区域阈值处理有时也可用于表示可变阈值处理。可变阈值处理的T值取决于(x,y)的邻域的特性。

? ? ? ? 上图中的图b就展示了一个可变阈值处理，包含有三个支配的直方图，利用表达式：

a，b，c是任意的三个不同的灰度值。?

? ? ? ?当图像中存在噪声时，其直方图也会发生不同的变化，下图给出了不同噪声值情况下图像直方图的例子。

从上图可以看出当噪声的灰度级越来越大时，对于图像直方图的影响也会发生巨大变化。?

? ? ? ? 除了噪声外，光照和反射同样会对直方图产生影响，同样给出例子。

不难看处，光照和反射对于直方图的影响也很大。

3.2 基本的全局阈值处理?

? ? ? ? 通常情况下，图像之间有较大的变化，纵使使用了全局阈值的方法计算其阈值也是需要的，这里就需要一个能够自动估计阈值的算法：

????????1、为全局阈值T选择一个初始估计值。

? ? ? ? 2、在公式中用 T分割该图像。这将产生两组像素: $G_{1}$ 由灰度值大于T的所有像素组成， $G_{2}$ 由所有小于等于T的像素组成。

? ? ? ? 3、对 $G_{1}$ 和 $G_{2}$ 的像素分别计算平均灰度值(均值) m和m2。

? ? ? ? 4、计算一个新的阈值：? ? ? ?

???????? $T = \frac{1}{2}(m_{1}+m_{2})$

? ? ? ? 5、重复步骤2到步骤4,直到连续迭代中的T值间的差小于一个预定义的参数△T为止。?

3.3 用Otsu方法的最佳全局阈值处理

? ? ? ? 阈值处理可视为一种统计决策理论问题，目的是在把像素分配给两个或多个组的过程中引入的平均误差最小。Otsu方法便是一种方案，该方法在类间方差最大的情况下是最佳的，也就是众所周知的统计鉴别分析中所用到的度量。当然除了最佳性外Otsu方法还有一个重要的特性，即可以完全以在一幅图像的直方图上执行计算为基础，直方图是可以容易得到的一维阵列。

Otsu算法的步骤：

1、计算输入图像的归一化直方图，并使用 $p_{i}$ 表示直方图上的各个分量，所有分量之和为1。

?2、对于0<k<L-1，计算雷击和 $P_{1}(k)$

3、同样的范围，计算累积均值m(k)

4、计算全局灰度均值 $m_{G}$

5、计算类间方差 $\sigma _{B}^{2}(k)$

6、得到阈值 $k^{*}$ ，若最大值不唯一，用相应检测到的各个最大值k的平均求得 $k^{*}$ ，其最大化 $\sigma _{B}^{2}(k^{*})$ 为：

7、最后计算得到可分性度量 $\eta ^{*}$

下面给出使用该方法处理图像的例子：

实验?

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

img = cv2.imread("D:\\picture\\test4.jpg", 0)

deltaT = 1  
histCV = cv2.calcHist([img], [0], None, [256], [0, 256])  # 灰度直方图
grayScale = range(256)  # 灰度级 [0,255]
totalPixels = img.shape[0] * img.shape[1]  # 像素总数
totalGray = np.dot(histCV[:, 0], grayScale)  
T = round(totalGray / totalPixels)  # 平均灰度
while True:
    numC1, sumC1 = 0, 0
    for i in range(T):  # 计算 C1: (0,T) 平均灰度
        numC1 += histCV[i, 0]  # C1 像素数量
        sumC1 += histCV[i, 0] * i  # C1 灰度值总和
    numC2, sumC2 = (totalPixels - numC1), (totalGray - sumC1)  # C2 像素数量, 灰度值总和
    T1 = round(sumC1 / numC1)  # C1 平均灰度
    T2 = round(sumC2 / numC2)  # C2 平均灰度
    Tnew = round((T1 + T2) / 2)  # 计算新的阈值
    print("T={}, m1={}, m2={}, Tnew={}".format(T, T1, T2, Tnew))
    if abs(T - Tnew) < deltaT:  
        break
    else:
        T = Tnew

# 阈值处理
ret2, img_Otsu = cv2.threshold(img, T, 255, cv2.THRESH_OTSU)  # 阈值分割
print(ret2)

plt.figure(figsize=(9,3))
plt.subplot(131), plt.axis('off'), plt.title("Origin"), plt.imshow(img, 'gray')
plt.subplot(132, yticks=[]), plt.title("Gray Hist")  
histNP, bins = np.histogram(img.flatten(), bins=255, range=[0, 255], density=True)
plt.bar(bins[:-1], histNP[:])
plt.subplot(133), plt.title("OTSU binary(T={})".format(round(ret2))), plt.axis('off')
plt.imshow(img_Otsu, 'gray')
plt.tight_layout()
plt.show()

3.4 图像平滑、边缘改善全局阈值处理?

噪声

? ? ? ? 噪声对于图像处理而言会产生巨大的影响，它会将一个简单的阈值处理问题变成一个不可解决的问题，所以在选择阈值处理进行分割处理时，需要在阈值处理之间先进行平滑图像。

从上图可以看出，当不对图像中的噪声加以处理时，使用Otsu方法进行分割时会出现很多黑点，分割很不成功。而进行平滑后的图像其直方图改进明显，利用Otsu阈值处理法处理后的图像也更加完美。

边缘

? ? ? ? 除了噪声外，波峰图的形状对于阈值处理的是否成功也会造成影响，对于大背景区域上一个小物体组成的图像的直方图将会有一个较大的波峰，因此需要对其进行改进。

????????如果仅用位于或接近物体和背景间的边缘的像素，则得到的直方图将有几个高度近似相同的波峰。此外，像素位于物体上的概率将近似等于其位于背景中的概率，从而改进了直方图模式的对称性最后，如下段中指出的那样，用满足某些基于梯度和拉普拉斯算子的简单度量的像素有加深直方图波峰间的波谷的倾向。

边缘实验?

img = cv2.imread("D:\\picture\\tupian.jpg", flags=0)


histCV1 = cv2.calcHist([img], [0], None, [256], [0, 256])  # 灰度直方图
totalPixels = img.shape[0] * img.shape[1]  # 像素总数
totalGray = np.dot(histCV1[:, 0], range(256))
meanGray = round(totalGray / totalPixels)  # 平均灰度
ret, imgBin = cv2.threshold(img, meanGray, 255, cv2.THRESH_BINARY)  # thresh=meanGray

# (1) 计算 Sobel 梯度算子
SobelX = cv2.Sobel(img, cv2.CV_32F, 1, 0)  # 计算 x 轴方向
SobelY = cv2.Sobel(img, cv2.CV_32F, 0, 1)
grad = np.sqrt(SobelX ** 2 + SobelY ** 2)
gradMax = np.int(np.max(grad))
# (2) 设置阈值 T=0.3*gradMax，对梯度图像进行阈值处理，作为遮罩模板
_, maskBW = cv2.threshold(np.uint8(grad), 0.3 * gradMax, 255, cv2.THRESH_BINARY)
# (3) 计算基于遮罩模板的直方图分布，以排除无效背景像素的影响
histCV2 = cv2.calcHist([img], [0], maskBW, [256], [0, 256])
histCV2[0] = 0
# (4) 排除无效背景像素影响后，进行阈值处理
Tmask = 120 
_, result = cv2.threshold(img, Tmask, 255, cv2.THRESH_BINARY)

print(gradMax, meanGray)
plt.figure(figsize=(9, 3))
plt.subplot(231), plt.axis('off'), plt.title("Origin"), plt.imshow(img, 'gray')
plt.subplot(232, yticks=[]), plt.axis([0, 255, 0, np.max(histCV1)])
plt.bar(range(256), histCV1[:, 0]), plt.title("Hist")
plt.subplot(233), plt.axis('off'), plt.title("Binary"), plt.imshow(imgBin, 'gray')
plt.subplot(234), plt.axis('off'), plt.title("Sobel")
plt.imshow(maskBW, cmap='gray')
plt.subplot(235, yticks=[]), plt.title("Hist of boundries")  
plt.bar(range(256), histCV2[:, 0])
plt.subplot(236), plt.axis('off'), plt.title("result")
plt.imshow(result, 'gray')
plt.tight_layout()
plt.show()

3.5 多阈值处理

? ? ? ? 在前面的讨论中虽然有提及多阈值情况的阈值处理，但处理的图像仍旧用的是单一的全局阈值方法，这里将对任意数量的阈值进行讨论。在K个类 $C_{1},C_{2},C_{3},C_{4}....C_{K}$ 的情况下，类间方差可写为：

其中

$m_{G}$ 是全局均值，在之前的讨论中有过推导，K类由K-1个阈值来分离。

? ? ? ? 在实际处理中，当相信两个阈值就可以有效解决问题的时候，使用多个全局阈值处理就被是为一种可行的方法。对于由三个灰度间隔组成的三个类，类间方差为：

且

$P_{1}m_{1}+P_{2}m_{2}+P_{3}m_{3}=m_{G},P_{1}+P_{2}+P_{3}=1$

其最佳阈值为：

由此就可以得到阈值处理过后的图像：

下图为多阈值处理的示例图?

实验?

双Otsu阈值算法的python实现

def double(img):
    histCV = cv2.calcHist([img], [0], None, [256], [0, 256])
    grayScale = np.arange(0, 256, 1)
    totalPixels = img.shape[0] * img.shape[1]
    totalGray = np.dot(histCV[:, 0], grayScale)
    mG = totalGray / totalPixels
    varG = sum(((i - mG) ** 2 * histCV[i, 0] / totalPixels) for i in range(256))

    T1, T2, varMax = 1, 2, 0.0
    for k1 in range(1, 254):
        n1 = sum(histCV[:k1, 0])  # C1 像素数量
        s1 = sum((i * histCV[i, 0]) for i in range(k1))
        P1 = n1 / totalPixels  # C1 像素数占比
        m1 = (s1 / n1) if n1 > 0 else 0  # C1 平均灰度

        for k2 in range(k1 + 1, 256):  # k2: [2,254], k2>k1
            n3 = sum(histCV[k2 + 1:, 0])  # C3 像素数量
            s3 = sum((i * histCV[i, 0]) for i in range(k2 + 1, 256))
            P3 = n3 / totalPixels  # C3 像素数占比
            m3 = (s3 / n3) if n3 > 0 else 0  # C3 平均灰度
            P2 = 1.0 - P1 - P3  # C2 像素数占比
            m2 = (mG - P1 * m1 - P3 * m3) / P2 if P2 > 1e-6 else 0  # C2 平均灰度
            var = P1 * (m1 - mG) ** 2 + P2 * (m2 - mG) ** 2 + P3 * (m3 - mG) ** 2
            if var > varMax:
                T1, T2, varMax = k1, k2, var

    epsT = varMax / varG  #
    print(totalPixels, mG, varG, varMax, epsT, T1, T2)
    return T1, T2, epsT

结果图：

3.6 可变阈值处理

? ? ? ? 可变阈值处理常用的方法有三种：图像分块、基于局部图像特性、使用平均移动。

????????当感兴趣物体和背景占据合理的可比大小的区域是，通常选用图像分块进行可变阈值处理；基于局部图像特性的可变阈值处理是以一个或多个在(x,y)邻域计算的特性为基础的。用一幅图像中每个点的一个邻域内的像素的标准差和均值来说明局部阈值处理的基本方法。这两个量对于确定局部阈值非常有用，因为它们是局部对比度和平均灰度的描述子；平均移动则一般使用在速度是一种基本要求的文档处理中。

图像分块示例：

基于局部图像特性的示例：

平均移动的示例：

四、基于区域的分割?

? ? ? ? 分割的目的就是将一幅图像划分为多个区域。在上面的讨论中主要是基于灰度级、像素特性等来研究分割。这里将讨论以直接寻找区域为基础的分割技术。

4.1 区域生长

? ? ? ? 根据预先定义的生长准则将像素或子区域组合为更大区域的过程就称为区域生长。利用区域生长来处理问题时通常可以基于问题的性质来选择一组或多组开始点，这一过程是在每个像素处计算相同特性的集合。最终，在生长过程期间，根据这些特性集合把像素分配给各个区域。若这些计算的结果显示了一族值，则特性靠近这些族的中心的像素可以作为种子。

? ? ? ? 若区域生长过程中没有使用连通属性，那么单独的描绘子会产生错误的结果。且当不再有像素满足加入某个区域的准则是，区域生长需要停止，这一中止的表示法也需要考虑。

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

# 用来计算两个像素的距离
def getGraydis(image, currentPoint, tmpPoint):
    return abs(int(image[currentPoint[0], currentPoint[1]]) - int(image[tmpPoint[0], tmpPoint[1]]))

# 区域生长算法
def growth(img, seeds, thresh=5):
    height, weight = img.shape
    seedMark = np.zeros(img.shape)
    seedList = []
    for seed in seeds:
        if (0 < seed[0] < height and 0 < seed[1] < weight): seedList.append(seed)
    label = 1  # 种子位置标记
    connects = [(-1, -1), (0, -1), (1, -1), (1, 0), (1, 1), (0, 1), (-1, 1), (-1, 0)]  # 8 邻接连通
    while (len(seedList) > 0):  # 如果列表里还存在点
        currentPoint = seedList.pop(0)  # 将最前面的那个抛出
        seedMark[currentPoint[0], currentPoint[1]] = label  # 将对应位置的点标记为 1
        for i in range(8):  # 对这个点周围的8个点一次进行相似性判断
            tmpX = currentPoint[0] + connects[i][0]
            tmpY = currentPoint[1] + connects[i][1]
            if tmpX < 0 or tmpY < 0 or tmpX >= height or tmpY >= weight:  # 是否超出限定阈值
                continue
            grayDiff = getGraydis(img, currentPoint, (tmpX, tmpY))  # 计算灰度差
            if grayDiff < thresh and seedMark[tmpX, tmpY] == 0:
                seedMark[tmpX, tmpY] = label
                seedList.append((tmpX, tmpY))
    return seedMark

# 区域生长 主程序
img = cv2.imread("D:\\picture\\test5.jpg", 0)
# 区域生长图像分割
imgBlur = cv2.blur(img, (3, 3))
_, imgTop = cv2.threshold(imgBlur, 250, 255, cv2.THRESH_BINARY)
nseeds, labels, stats, centroids = cv2.connectedComponentsWithStats(imgTop)  
seeds = centroids.astype(int)  
result = growth(img, seeds, 8)
plt.figure(figsize=(8, 4))
plt.subplot(121), plt.axis('off'), plt.title("Origin")
plt.imshow(img, 'gray')
plt.subplot(122), plt.axis('off'), plt.title("grow")
plt.imshow(255 - result, 'gray')
plt.tight_layout()
plt.show()

4.2 区域分裂与聚合

????????除了区域生长外，还有另一种可供选择的方法，先将一幅图像细分为一组任意的不相交区域，然后聚合或分裂这些区域。下面为示意图：

步骤为：

1、对满足 $Q(R_{i})=FALSE$ 的任何区域 $R_{i}$ 分裂为4个不相交的象限区域。

2、当不可能进一步分裂时，对满足条件 $Q(R_{j}\bigcup R_{k})=TRUE$ 的任意两个邻接区域 $R_{j}$ 和 $R_{k}$ 进行聚合。

3、当无法进一步操作时，停止操作。?

五、形态学分水岭的分割

? ? ? ? 分水岭的概念是以三维形象化一幅图像为基础的，两个空间坐标作为灰度的函数。在地形学中，考虑三种类型的点：(1)属于一个区域最小值的点(2)把一点看成是一个水滴，将点放在任意位置上，水滴一定会下落到一个单一的最小值点(c)处在该点的水会等可能性的流向不止一个这样的最小值点。

? ? ? ? 基于这样的概念的分割算法的主要目标就是要找出分水线，假设在每个区域的最小值上打一个洞，并且让水通过洞以均匀的速率上升，从低到高淹没整个地形。当不同的汇水盆地中上升的水聚集时，修建一个水坝来阻止这种聚合。水将达到在水线上只能见到各个水坝的顶部这样一个程度。这些大坝的边界对应于分水岭的分割线。因此，它们是由分水岭算法提取出来的(连接的)边界。可以借助下图进行理解。

分水岭分割算法

? ? ? ? 令 $M_{1},M_{2},M_{3}.....$ 是表示图像g(x,y)的区域最小值点的坐标的集合。这时一幅典型的梯度图像。令 $C(M_{i})$ 是表示与区域最小值 $M_{i}$ 相联系的汇水盆地中的点的坐标集合。符号min和max表示图像的最小值和最大值，令T[n]表示满足g(s,t)<n的坐标(s,t)的集合：

当水位上升时，地形将会被淹没，在这个过程中算法需要知道淹没深度之下的点的数量，由此就将处在n平面之下的坐标标记为黑色，其他为白色，当向下观察xy平面时就会得到一幅二值图像。令 $C_{n}(M_{i})$ 表示汇水盆地中与淹没阶段n的最小值 $M_{i}$ 相关联的点的坐标集。由此可以将其看成是如下公式给出的二值图像：

? ? ? ? 接着，令C[n]表示在阶段n中已被水淹没的汇水盆地的“并”：

C[max+1]则表示所有汇水盆地的并：

通过证明可以知道在执行过程中? $C_{n}(M_{i})$ 和T[n]中的元素是不会被替换的，n增大时集合中的元素不是增加就是保持相同。

示例图：

python实现形态学梯度的分割算法：

img = cv2.imread("D:\\picture\\tupian.jpg", flags=1)  # 读取彩色图像(BGR)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)  # 转为灰度图像

# 图像的形态学梯度
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (3, 3))  # 生成 3*3 结构元
grad = cv2.morphologyEx(gray, cv2.MORPH_GRADIENT, kernel)  #梯度

# 阈值分割，将灰度图像分为黑白二值图像
_, thresh = cv2.threshold(np.uint8(grad), 0.1 * grad.max(), 255, cv2.THRESH_BINARY)
# 形态学操作，生成 "确定背景" 区域
opening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel, iterations=2)  # 开运算，消除噪点
img_back = cv2.dilate(opening, kernel, iterations=3)  # 膨胀操作，生成 "确定背景" 区域
# 距离变换，生成 "确定前景" 区域
distance = cv2.distanceTransform(opening, cv2.DIST_L2, 5)  
_, img_front = cv2.threshold(distance, 0.1 * distance.max(), 255, 0)  # 阈值选择 0.1*max 效果较好
img_front = np.uint8(img_front)
# 连通域处理
ret, component = cv2.connectedComponents(img_front, connectivity=8)  # 对连通区域进行标号，序号为 0-N-1
markers = component + 1  # OpenCV 分水岭算法设置标注从 1 开始，而连通域编从 0 开始
kinds = markers.max()  # 标注连通域的数量
maxKind = np.argmax(np.bincount(markers.flatten()))  # 出现最多的序号，所占面积最大，选为底色
markersBGR = np.ones_like(img) * 255
for i in range(kinds):
    if (i != maxKind):
        colorKind = [np.random.randint(0, 255), np.random.randint(0, 255), np.random.randint(0, 255)]
        markersBGR[markers == i] = colorKind
# 去除连通域中的背景区域部分
unknown = cv2.subtract(img_back, img_front)
markers[unknown == 255] = 0
# 分水岭算法标注目标的轮廓
markers = cv2.watershed(img, markers)

# 把轮廓添加到原始图像上
img_yingshe = img.copy()
img_yingshe[markers == -1] = [0, 0, 255]
print(img.shape, markers.shape, markers.max(), markers.min(), ret)

plt.figure(figsize=(9, 3))
plt.subplot(131), plt.axis('off'), plt.title("Origin")
plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.subplot(132), plt.axis('off'), plt.title("tidu")
plt.imshow(grad, 'gray') 
plt.subplot(133), plt.axis('off'), plt.title("yingshe")
plt.imshow(cv2.cvtColor(img_yingshe, cv2.COLOR_BGR2RGB))
plt.tight_layout()
plt.show()