数字图像处理知识点梳理——第十章 图像分割

引言

1、什么是分割？

分割指把输入的图像的部分属性提取出来，它可将图像细分为构成它的子区域或物体，细分的程度要看具体要解决的问题。其中，异常图像的分割是图像处理中最困难的任务之一，分割的精细程度决定着计算分析过程最终的成败。

2、分割的原理是什么？

本章中的大多数分割算法均基于灰度值的两个基本性质之一：不连续性和相似性。

第一类中，方法是以灰度突变为基础分割一幅图像。如点检测、线检测、边缘检测。

第二类中，方法是根据一组预定义的准则将一幅图像分割为相似的区域。 如门限(阈值)分割、区域生长、区域分裂和合并。

3、分割的要求是什么？

（1）分割必须是完全的；

（2）区域中的点以某些预定义的方式来连接；

（3）各个区域必须是不相交的；

（4）分割后的区域中的像素必须满足属性；

（5）两个邻接区域R_i和R_j在属性Q的意义上必须是不同的。

分割图像的物体可以通过确定图像中的物体边界来完成，由此引出接下来的主题，不连续检测。

不连续检测

首先介绍点检测。

1、点检测

点的检测以二阶导数为基础。使用拉普拉斯
$${\nabla }^{2}f(x,y)=\frac{{\partial }^{2}f}{\partial x^2}+\frac{{\partial }^{2}f}{\partial y^2}$$

$${\nabla }^{2}f(x,y)=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y)$$

输出为
$$g(x,y)=\{\begin{array}{l}1,|R(x,y)|⩾T\\ 0,\text{其他}\end{array}，T是一个非负的阈值$$
实验：

经过[-1 -1 -1 ; -1 8 -1 ; -1 -1 -1]，得到点检测的图像（右图）

import numpy as np
import cv2

img = cv2.imread('C:/Users/13121/Desktop/8.png',0)
cv2.imshow("original",img)
G = np.array([[-1, -1, -1],[-1, 8, -1],[-1, -1, -1]])
def spots(img,threshold):
    rows = np.size(img, 0)
    columns = np.size(img, 1)
    mag = np.zeros(img.shape)
    for i in range(0, rows - 2):
        for j in range(0, columns - 2):
            v = sum(sum(G * img[i:i+3, j:j+3])) 
            mag[i+1, j+1] = np.sqrt((v ** 2))       
    for p in range(0, rows):
        for q in range(0, columns):
            if mag[p, q] < threshold:
                mag[p, q] = 0
    return mag
img1 = spots(img,0)
cv2.imshow("spots",img1)     

cv2.waitKey(0)
cv2.destroyAllWindows() 

2、线检测

预期二阶导数将导致更强的响应。二阶导数的双线效应应必须做适当的处理。

拉普拉斯测子是各向同性的，因此其响应与方向无关。

如果对检测图像中由给定模板定义的方向上的所有线感兴趣，则只需简单地对该图像运行这个模板，并对结果的绝对值进行阈值处理，留下的点是有最强响应的点，对于1个像素宽度的线来说，相应的点最接近于模板定义的方向。

3、边缘检测

介绍边缘检测之前，要先了解到一组基础知识，一阶导数和二阶导数。

• 一阶导数为基础

$$\frac{\partial f}{\partial x}=f^{\prime }(x)=f(x+1)-f(x)$$

​ 恒定灰度区域为0，灰度台阶或斜坡处不为0

• 二阶导数为基础

$$\begin{array}{rl}\frac{{\partial }^{2}f}{\partial x^2}& =\frac{\partial f^{\prime }(x)}{\partial x}=f^{\prime }(x+1)-f^{\prime }(x)\\ & =f(x+2)-f(x+1)-f(x+1)+f(x)\\ & =f(x+2)-2f(x+1)+f(x)\end{array}$$
恒定灰度为0，灰度台阶或斜坡开始处、结束处不为0，沿灰度斜坡为0

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

接下来进入边缘检测。

边缘检测是根据图像的灰度剖面来分类的。包括台阶模型、斜坡模型、屋顶边缘模型。

边缘模型	内容	灰度剖面
台阶模型	一个像素的距离发生两个灰度级间理想的过度
斜坡模型	斜坡斜度与边缘模糊程度成反比。一个边缘点现在是斜坡中包含的任何点，而一条边缘线段将是一组一连接起来的这样的点
屋顶边缘模型	屋顶边缘的基地（宽度）由该线的宽度和尖锐度决定（当其基地为1个像素宽时，屋顶边缘时一条穿过图像中一个区域的1像素宽的线）

基本的边缘检测方法

首先来了解一些基础知识：
$$梯度图像：M(x,y)=\mathrm{mag}(\nabla f)=\sqrt{g_x^2+g_y^2}$$

$$梯度向量的方向：\alpha (x,y)=\arctan \left[\frac{g_y}{g_x}\right]$$

（后续处理中，因M(x,y)计算不便，改成取绝对值相加的方式计算，即$M(x,y)=|g_x|+|g_y|$）

其中，
$$\begin{array}{l}{g}_x=\frac{\partial f(x,y)}{\partial x}=f(x+1,y)-f(x,y)\\ g_y=\frac{\partial f(x,y)}{\partial y}=f(x,y+1)-f(x,y)\end{array}$$
（α(x,y)、g_x、g_y、M(x,y)图像尺寸相同）

接下来学习基本的边缘检测方法（基于一阶导数）：

边缘检测算子名称	相应算子（水平方向）	相应算子（对角方向）
Roberts
Prewitt
Sobel

我们以Sobel算子为例进行分析。

（原图）

（左图为x方向，右图为y方向，设阈值均为20）

（综合的Sobel算子，设阈值为40）

可见,图像检测中，对于精细的细节（如图中的墙上的砖）会引起噪声。对于这种现象，我们要先对图像进行平滑处理，再进行边缘检测。

例如：对原图进行一个5×5的中值滤波平滑，得到的结果为

平滑图像后进行边缘检测的结果为


可见，经过平滑后再边缘检测的图像的噪声减少了，边缘更加清晰。

import numpy as np
import cv2

img = cv2.imread('C:/Users/13121/Desktop/9.png',0)
cv2.imshow("original",img)

img = cv2.medianBlur(img, 5)
cv2.imshow("img_MedianBlur_5", img)

def sobel_v(img, threshold):
    G_x = np.array([[-1, 0, 1],[-2, 0, 2],[-1, 0, 1]])
    rows = np.size(img, 0)
    columns = np.size(img, 1)
    mag = np.zeros(img.shape)
    for i in range(0, rows - 2):
        for j in range(0, columns - 2):
            v = sum(sum(G_x * img[i:i+3, j:j+3]))  # vertical
            mag[i+1, j+1] = v           
    for p in range(0, rows):
        for q in range(0, columns):
            if mag[p, q] < threshold:
                mag[p, q] = 0
    return mag
def sobel_h(img, threshold):
    G_y = np.array([[-1, -2, -1],[0, 0, 0],[1, 2, 1]])
    rows = np.size(img, 0)
    columns = np.size(img, 1)
    mag = np.zeros(img.shape)
    for i in range(0, rows - 2):
        for j in range(0, columns - 2):
            h = sum(sum(G_y * img[i:i+3, j:j+3]))  # horizon
            mag[i+1, j+1] = h
            
    for p in range(0, rows):
        for q in range(0, columns):
            if mag[p, q] < threshold:
                mag[p, q] = 0
    return mag
def sobel(img, threshold):
    G_x = np.array([[-1, 0, 1],[-2, 0, 2],[-1, 0, 1]])
    G_y = np.array([[-1, -2, -1],[0, 0, 0],[1, 2, 1]])
    rows = np.size(img, 0)
    columns = np.size(img, 1)
    mag = np.zeros(img.shape)
    for i in range(0, rows - 2):
        for j in range(0, columns - 2):
            v = sum(sum(G_x * img[i:i+3, j:j+3]))  # vertical
            h = sum(sum(G_y * img[i:i+3, j:j+3]))  # horizon
            mag[i+1, j+1] = np.sqrt((v ** 2) + (h ** 2))       
    for p in range(0, rows):
        for q in range(0, columns):
            if mag[p, q] < threshold:
                mag[p, q] = 0
    return mag

edge_img1 = sobel_v(img,20)
cv2.imshow("edge1_img",edge_img1)
edge_img2 = sobel_h(img,20)
cv2.imshow("edge2_img",edge_img2)
edge_img3 = sobel(img,40)
cv2.imshow("edge_img",edge_img3)

cv2.waitKey(0)
cv2.destroyAllWindows() 

这一小节涉及的三种算子以用一个或多个模板对图像进行滤波为基础，而未对图像特性和噪声内容采取预防措施。接下来考虑诸如图像噪声和边缘本身特性的因素改进简单的边缘检测技术。

更先进的边缘检测技术

1、Marr-Hildreth边缘检测器

• 高斯拉普拉斯（LoG）

$$G(x,y)={\mathrm{e}}^{-\frac{x^2+y^2}{2{\sigma }^2}}$$

$${\nabla }^{2}G(x,y)=\left[\frac{x^2+y^2-2{\sigma }^2}{{\sigma }^4}\right]{\mathrm{e}}^{-\frac{x^2+y^2}{2{\sigma }^2}}$$

$$g(x,y)={\nabla }^2[G(x,y)\text{ ★ }f(x,y)]$$

该检测算法小结：1.用一个对式$G(x,y)={\mathrm{e}}^{-\frac{x^2+y^2}{2{\sigma }^2}}$取样得到的n×n的高斯低通滤波器对输入图像滤波。2.计算由第一步得到的图像的拉普拉斯3.找到步骤2所的图像的零交叉

• 高斯差分（DoG）

$$\mathrm{DoG}(x,y)=\frac{1}{2\pi {\sigma }_1^2}{\mathrm{e}}^{\frac{x^2+y^2}{2{\sigma }_1^2}}-\frac{1}{2\pi {\sigma }_2^2}{\mathrm{e}}^{-\frac{x^2+y^2}{2{\sigma }_2^2}}，其中{\sigma }_1>{\sigma }_2$$

2、坎尼边缘检测器

坎尼检测算法基本步骤：1.用一个高斯滤波器平滑输入图像2.计算梯度幅值图像和角度图像3.对梯度幅值图像应用非最大抑制4.用双阈值处理和连接分析来检测并连接边缘

但当目的是突出主要边缘并尽可能维护连接性时，事件中通常平滑处理和阈值处理两者都使用。

接下来学习阈值处理

阈值处理

先了解基本知识：
$$g(x,y)=\{\begin{array}{llc}1,& f(x,y)>T& 此时(x,y)称为对象点\\ 0,& f(x,y)⩽T& 此时(x,y)称为背景点\end{array}$$
当T是一个适用于整个图像的常数时，该公式给出的处理称为全局阈值处理。

当T值在一幅图像上改变时，我们称为可变阈值处理。

局部阈值处理或区域阈值处理有时用于表示可变阈值处理。如果T取决于空间坐标(x,y)本身，则可变阈值处理通常称为动态阈值处理或自适应阈值处理

灰度阈值的成功与否直接关系到可区分的直方图模式的谷的宽度和深度，影响波谷特性的关键因素是：（1）波峰间的间隔（波峰离得越远，分离这些模式的机会越好）（2）图像中的噪声内容（模式随噪声的增加而展宽）（3）物体和背景的相对尺寸（4）光源的均匀性（5）图像反射特性的均匀性。

接下来学习阈值处理的具体内容：

1、基本的全局阈值处理

当物体和背景像素的灰度分布十分明显时，可以用适用于整个图像的单个（全局）阈值。

计算流程为：

1. 为全局阈值T选择一个初始估计值
2. 用式$g(x,y)=\{\begin{array}{ll}1,& f(x,y)>T\\ 0,& f(x,y)⩽T\end{array}$来分割图像，
   G₁由灰度值大于T的所有像素组成，G₂由所有小于等于T的像素组成
3. 对G₁和G₂的像素分别计算平均灰度值（均值）m₁和m₂
4. 计算新阈值$T=\frac{1}{2}\left(m_1+m_2\right)$
5. 重复步骤2到步骤4，直到连续迭代种的T值间的差小于一个预定义的参数ΔT为止。

原图为

使用全局阈值后为：

import cv2 as cv
import numpy as np
# 转灰
from matplotlib.pyplot import hist
def rgb2gray(img):
    h=img.shape[0]
    w=img.shape[1]
    img1=np.zeros((h,w),np.uint8)
    for i in range(h):
        for j in range(w):
            img1[i,j]=0.144*img[i,j,0]+0.587*img[i,j,1]+0.299*img[i,j,2]
    return img1

# 计算新阈值
def threshold(img,T):
    h=img.shape[0]
    w=img.shape[1]
    G1=G2=0
    g1=g2=0
    for i in range (h):
        for j in range (w):
            if img[i,j]>T:
                G1+=img[i,j]
                g1+=1
            else:
                G2+=img[i,j]
                g2+=1
    m1=int(G1/g1)
    m2=int(G2/g2)   # m1，m2计算两组像素均值
    T0=int((m1+m2)/2)   # 据公式计算新的阈值
    return T0

def decide(img,T):
    h=img.shape[0]
    w=img.shape[1]
    img1=np.zeros((h,w),np.uint8)
    T0=T
    T1=threshold(img,T0)
    for k in range (100):   # 迭代次数为经验值，可据实际情况选定
        if abs(T1-T0)==0:   # 若新阈值减旧阈值差值为零，则为二值图最佳阈值
            for i in range (h):
                for j in range (w):
                    if img[i,j]>T1:
                        img1[i,j]=255
                    else:
                        img1[i,j]=0
            break
        else:
            T2=threshold(img,T1)
            T0=T1
            T1=T2   # 变量转换，保证if条件为新阈值减旧阈值
    return img1

image=cv.imread('C:/Users/13121/Desktop/6.jpg')
image=cv.resize(image,(224,224),interpolation=cv.INTER_CUBIC)
grayimage=rgb2gray(image)
thresholdimage=decide(grayimage,127)
cv.imshow("image",image)
cv.imshow("grayimage",grayimage)
cv.imshow("thresholdimage",thresholdimage)
cv.waitKey(0)
cv.destroyAllWindows()

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

阈值处理可视为一种统计决策理论问题，其目的是在把像素分配给两个或多个组（也称为分类）的过程中引入的平均误差最小。

Otsu方法的基本原理是以最佳阈值将图像的灰度直方图分割成两部分，使两部分之间的方差取最大值，即分离性最大。

原图为：

使用Otsu方法后为：

ret1, dst1 = cv.threshold(grayimage, 0 , 255, cv.THRESH_OTSU)
cv.imshow("otsu", dst1)

3、可变阈值

t, dst = cv.threshold(grayimage, 127, 255, cv.THRESH_TOZERO)
cv.imshow("dst", dst)

全局阈值的改善

1、用图像平滑改善全局阈值处理

噪声可将一个简单的阈值处理问题变为一个不可解决的问题，当噪声不能在源头减少，并且阈值处理又是所选择的分割方法时，那么通常可增强性能的一种技术是在阈值处理之前平滑图像。

我们平滑一幅图像侵蚀的越多，分割后的结果中的边界误差就越大。

区域太小，则该区域对直方图的贡献与由噪声引起的灰度扩散相比无足轻重。

2、利用边缘改进全局阈值处理

改进直方图形状的一种方法是仅考虑那些位于或靠近物体和背景间的边缘的像素，一种直接且明显的改进是直方图很少依赖物体和背景的相对大小。

边缘连接

在进行边缘检测时，会出现细节断裂的情况，因此就引入了本小节的内容，边缘连接。边缘连接又分为局部的连接处理、区域的连接处理、和使用霍夫变换的全局连接处理。

局部处理：
连接边缘点的最简单的方法之一是在每个点(x,y)处的一个小邻域内分析像素的特点。
用于确定边缘像素相似性的两个主要性质是：（1）梯度向量的强度（幅度）（2）梯度向量的方向

区域处理：
通常图像中感兴趣区域的位置已知或可以确定，这种情况下使用区域的基础上连接像素的技术，期望的结果是该区域边界的近似。
必须指定两个起始点，并且所有的点都必须排序（即按顺时针方向或逆时针方向）

使用霍夫变换的全局处理：
Hough变换可以用于将边缘像素连接起来得到边界曲线。
在已知曲线形状的条件下，Hough变换实际上是利用分散的边缘点进行曲线逼近，它也可看成是一种聚类分析技术。

利用Hough变换检测直线

1. 由图像中直线斜率可能的范围，确定好参数空间中a的范围(a_min, a_max)；
2. 由图像中直线截距可能的范围，确定好参数空间中b的范围(b_min, b_max)；
3. 将(a_min, a_max)和(b_min,bmax)离散，建立一个二维网格数组A，始化为零 。
4. 计算图像空间中每个点(x_i, y_i)的Hough变换：固定a_i∈(a_min,a_max)，计算b_i= −x_ia_i+y_i，对数组A中对应于(a_i , b_i)的位置加1。
5. 计算结束后, A(a,b)的值就是图像空间中落在以a为斜率,b为截距的直线上点的数目。
6. 设定阈值，点数高于阈值的为检测到的直线。

图像分割

1、形态学分水岭分割

分水岭阈值选择算法可以看成是一种自适应的多阈值分割算法。

分水岭算法不是简单地将图像在最佳灰度级进行阈值处理，而是从一个偏低但仍然能正确分割各个物体的阈值开始。然后随着阈值逐渐上升到最佳值，使各个物体不会被合并。这个方法可以解决那些由于物体靠得太近而不能用全局阈值解决的问题。

原图为：

用分水岭分割后的结果图为：

import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt

img = cv.imread('C:/Users/13121/Desktop/11.png')
cv.imshow("original",img)
gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
ret, thresh = cv.threshold(gray, 0, 255, cv.THRESH_BINARY_INV + cv.THRESH_OTSU)
# 去噪处理
kernel = np.ones((3, 3), np.uint8)
opening = cv.morphologyEx(thresh, cv.MORPH_OPEN, kernel, iterations=2)  # 开运算
sure_bg = cv.dilate(opening, kernel, iterations=2)  # 膨胀操作

# Finding sure foreground area
dist_transform = cv.distanceTransform(opening, cv.DIST_L2, 3)
ret, sure_fg = cv.threshold(dist_transform, 0.6 * dist_transform.max(), 255, 0)  # 距离背景点足够远的点认为是确定前景

# Finding unknown region
sure_fg = np.uint8(sure_fg)
unknown = cv.subtract(sure_bg, sure_fg)  # 确定未知区域：减法运算
# Marker labelling
ret, markers = cv.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 = cv.watershed(img, markers)
img[markers == -1] = [0, 0, 255]

cv.imshow('final',img)
cv.waitKey(0)
cv.destroyAllWindows()

2、基于区域的分割

阈值分割法由于没有或很少考虑空间关系，使得阈值分割方法受到限制，基于区域的分割方法可以弥补这点不足，它利用的是图像的空间性质，该方法认为分割出来的属于同一区域的像素应具有相似的性质。

传统的区域分割算法有区域生长法和区域分裂合并法

区域生长法使得图像由小变大，区域分裂合并法使得图像由大变小

1、区域生长
区域生长是一种根据事先定义的准则将像素或子区域聚合成更大区域的过程，主要考虑像素及其空间邻域像素之间的关系。

区域生长的步骤：
1、选择合适的种子点
2、确定相似性准则(生长准则)
3、确定生长停止条件

2、区域分裂合并

算法流程：

1. 设整幅图像为初始区域。
2. 对每一区域R，如果P( R )=FLASE，则把该区域分裂成四个子区域。
3. 重复上一步，直到没有区域可以分裂。
4. 对图像中任意两个相邻的R1和R2，如果P(R1UR2)=TRUE，则把这两个区域合 并成一个区域。
5. 重复上一步，直到没有相邻区域可以合并，算法结束。