第十章 图像分割
一、基础知识
令R表示一幅图像占据的整个空间区域。我们可以将图像分割视为把R分为n个子区域
R
1
,
R
2
,
.
.
.
,
R
n
R_{1},R_{2},...,R_{n}
R1,R2,...,Rn的过程,满足:
(a)
⋃
i
=
1
n
R
i
=
R
\bigcup_{i=1}^{n}R_{i}=R
⋃i=1nRi=R。
(b)
R
i
R_{i}
Ri是一个连通集,
i
=
1
,
2
,
.
.
.
,
n
i=1,2,...,n
i=1,2,...,n。
(c)
R
i
∩
R
j
=
∅
,
R_{i}\cap R_{j}=\varnothing ,
Ri∩Rj=∅,对于所有
i
和
j
,
i
≠
j
i和j,i\neq j
i和j,i=j。
(d)
Q
(
R
i
)
=
T
R
U
E
,
i
=
1
,
2
,
.
.
.
,
n
Q(R_{i})=TRUE,i=1,2,...,n
Q(Ri)=TRUE,i=1,2,...,n。
(e)
Q
(
R
i
∪
R
j
)
=
F
A
L
S
E
Q(R_{i}\cup R_{j})=FALSE
Q(Ri∪Rj)=FALSE,对于任何
R
i
和
R
j
R_{i}和R_{j}
Ri和Rj的邻接区域。
条件(a)指出,分割必须是完全的;也就是说每个像素都必须在一个区域内。条件(b)要求一个区域中的点以某些预定义的方式来连接(即这些点必须是4连接的或8连接的)。条件(c)指出,各个区域必须是不相交的。条件(d)涉及分割后的区域中的像素必须满足的属性。最后,条件(e)指出,两个邻接区域
R
i
R_{i}
Ri和
R
j
R_{j}
Rj在属性Q的意义上必须是不同的。
分割的目的是将图像划分为不同区域。三大分割方法为:根据区域间灰度不连续搜寻区域之间的边界;以像素性质的分布进行阈值处理;直接搜寻区域进行分割。分别在点、线和边缘检测;阈值处理;基于区域的分割中介绍。
二、点、线和边缘检测
2.1点检测
点的检测以二阶导数为基础。使用拉普拉斯算子
▽
2
f
(
x
,
y
)
=
∂
2
f
∂
2
x
+
∂
2
f
∂
2
y
\triangledown ^{2}f(x,y)=\frac{\partial ^{2}f}{\partial ^{2}x}+\frac{\partial ^{2}f}{\partial ^{2}y}
▽2f(x,y)=∂2x∂2f+∂2y∂2f
▽
2
f
(
x
,
y
)
=
f
(
x
+
1
,
y
)
+
f
(
x
−
1
,
y
)
+
f
(
x
,
y
+
1
)
+
f
(
x
,
y
−
1
)
−
4
f
(
x
,
y
)
\triangledown ^{2}f(x,y)=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y)
▽2f(x,y)=f(x+1,y)+f(x−1,y)+f(x,y+1)+f(x,y−1)−4f(x,y)
输出
g
(
x
,
y
)
=
{
1
,
∣
R
(
x
,
y
)
∣
≥
T
0
,
其他
g(x,y)=\left\{\begin{matrix} 1, &\left | R(x,y) \right | \geq T\\ 0, & 其他 \end{matrix}\right.
g(x,y)={1,0,∣R(x,y)∣≥T其他
如下图所示的模板:
如果
∣
R
∣
≥
T
\left | R \right |\geq T
∣R∣≥T,则模板中心位置检测到一个点。其中,T是阈值,R是模板计算值。
如果一个孤立点与它周围的点不 同,则可以使用上述模板进行检测。如果模板响应为0,则表示在灰度级为常数的区域。
2.2线检测
四个线检测模板:
第一个模板对水平线有最大响应,第二个模板对
4
5
∘
45^{\circ}
45∘方向线有最大响应,第三个模板对垂直线有最大响应,第四个模板对
−
4
5
∘
-45^{\circ}
−45∘方向线有最大响应。
用
R
1
,
R
2
,
R
3
和
R
4
R_{1},R_{2},R_{3}和R_{4}
R1,R2,R3和R4分别代表水平、
4
5
∘
45^{\circ}
45∘、垂直和
−
4
5
∘
-45^{\circ}
−45∘方向线的模板响应,在图像中心的点,如过
∣
R
i
∣
>
∣
R
j
∣
,
j
≠
i
\left | R_{i} \right | >\left | R_{j}\right |,j≠i
∣Ri∣>∣Rj∣,j=i,则此点被认为与在模板i方向上的线更相关。在灰度恒定的区域,上述4个模板的响应为零。
2.3边缘检测
2.3.1梯度
图像f(x,y)在位置(x,y)的梯度用向量定义为
▽
f
=
[
g
x
g
y
]
=
[
∂
f
∂
x
∂
f
∂
y
]
\triangledown f=\begin{bmatrix} g_{x}\\g_{y} \end{bmatrix}=\begin{bmatrix} \frac{\partial f}{\partial x}\\ \frac{\partial f}{\partial y} \end{bmatrix}
▽f=[gxgy]=[∂x∂f∂y∂f]
向量
▽
f
\triangledown f
▽f的大小表示为
M
(
x
,
y
)
M(x,y)
M(x,y),即
M
(
x
,
y
)
=
m
a
g
(
▽
f
)
=
g
x
2
+
g
y
2
M(x,y)=mag(\triangledown f)=\sqrt{g_{x}^{2}+g_{y}^{2}}
M(x,y)=mag(▽f)=gx2+gy2
梯度向量的方向由下列对于x轴量度的角度给出:
α
(
x
,
y
)
=
a
r
c
t
a
n
[
g
y
g
x
]
\alpha (x,y)=arctan\begin{bmatrix} \frac{g_{y}}{g_{x}} \end{bmatrix}
α(x,y)=arctan[gxgy]
2.3.2梯度算子
边缘检测算子模板如下图所示:
Roberts交叉梯度算子:
G
x
=
Z
9
−
Z
5
G_{x}=Z_{9}-Z_{5}
Gx=Z9−Z5
G
y
=
Z
8
−
Z
6
G_{y}=Z_{8}-Z_{6}
Gy=Z8−Z6
Prewitt边缘检测算子:
G
x
=
(
Z
7
+
Z
8
+
Z
9
)
−
(
Z
1
+
Z
2
+
Z
3
)
G_{x}=(Z_{7}+Z_{8}+Z_{9})-(Z_{1}+Z_{2}+Z_{3})
Gx=(Z7+Z8+Z9)−(Z1+Z2+Z3)
G
y
=
(
Z
3
+
Z
6
+
Z
9
)
−
(
Z
1
+
Z
4
+
Z
7
)
G_{y}=(Z_{3}+Z_{6}+Z_{9})-(Z_{1}+Z_{4}+Z_{7})
Gy=(Z3+Z6+Z9)−(Z1+Z4+Z7)
Sobel算子:
G
x
=
(
Z
7
+
2
Z
8
+
Z
9
)
−
(
Z
1
+
2
Z
2
+
Z
3
)
G_{x}=(Z_{7}+2Z_{8}+Z_{9})-(Z_{1}+2Z_{2}+Z_{3})
Gx=(Z7+2Z8+Z9)−(Z1+2Z2+Z3)
G
y
=
(
Z
3
+
2
Z
6
+
Z
9
)
−
(
Z
1
+
2
Z
4
+
Z
7
)
G_{y}=(Z_{3}+2Z_{6}+Z_{9})-(Z_{1}+2Z_{4}+Z_{7})
Gy=(Z3+2Z6+Z9)−(Z1+2Z4+Z7)
2.3.3Python实现边缘检测的梯度算子
代码:
#边缘检测的梯度算子 (Roberts 算子, Prewitt 算子, Sobel 算子)
import numpy as np
import cv2
from matplotlib import pyplot as plt
img = cv2.imread("JMU.jpg", flags=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]])
# 卷积运算
imgBlur = cv2.blur(img, (3, 3)) # Blur 平滑后再做 Laplacian 变换
imgRoberts_x = cv2.filter2D(img, -1, kernel_Roberts_x)
imgRoberts_y = cv2.filter2D(img, -1, kernel_Roberts_y)
imgRoberts = np.uint8(cv2.normalize(abs(imgRoberts_x) + abs(imgRoberts_y), None, 0, 255, cv2.NORM_MINMAX))
imgPrewitt_x = cv2.filter2D(img, -1, kernel_Prewitt_x)
imgPrewitt_y = cv2.filter2D(img, -1, kernel_Prewitt_y)
imgPrewitt = np.uint8(cv2.normalize(abs(imgPrewitt_x) + abs(imgPrewitt_y), None, 0, 255, cv2.NORM_MINMAX))
imgSobel_x = cv2.filter2D(img, -1, kernel_Sobel_x)
imgSobel_y = cv2.filter2D(img, -1, kernel_Sobel_y)
imgSobel = np.uint8(cv2.normalize(abs(imgSobel_x) + abs(imgSobel_y), None, 0, 255, cv2.NORM_MINMAX))
plt.figure(figsize=(12, 8))
plt.subplot(341), plt.title('Origin'), plt.imshow(img, cmap='gray'), plt.axis('off')
plt.subplot(342), plt.title('Roberts'), plt.imshow(imgRoberts, cmap='gray'), plt.axis('off')
plt.subplot(346), plt.title('Roberts_X'), plt.imshow(imgRoberts_x, cmap='gray'), plt.axis('off')
plt.subplot(3, 4, 10), plt.title('Roberts_Y'), plt.imshow(imgRoberts_y, cmap='gray'), plt.axis('off')
plt.subplot(343), plt.title('Prewitt'), plt.imshow(imgPrewitt, cmap='gray'), plt.axis('off')
plt.subplot(347), plt.title('Prewitt_X'), plt.imshow(imgPrewitt_x, cmap='gray'), plt.axis('off')
plt.subplot(3, 4, 11), plt.title('Prewitt_Y'), plt.imshow(imgPrewitt_y, cmap='gray'), plt.axis('off')
plt.subplot(344), plt.title('Sobel'), plt.imshow(imgSobel, cmap='gray'), plt.axis('off')
plt.subplot(348), plt.title('Sobel_X'), plt.imshow(imgSobel_x, cmap='gray'), plt.axis('off')
plt.subplot(3, 4, 12), plt.title('Sobel_Y'), plt.imshow(imgSobel_y, cmap='gray'), plt.axis('off')
plt.tight_layout()
plt.show()
(代码参考链接:http://t.csdn.cn/hxkoF)
结果:
分析:
x方向上的梯度分量,水平细节更清楚;y方向上的梯度分量,垂直细节更清楚;梯度图像则水平和垂直细节都很清楚。
Roberts对边缘定位不是很准确,Prewitt和Sobel算子是计算数字梯度时常用的算子,Prewitt模板比Sobel模板简单,但Sobel模板能够有效抑制噪声。
2.4更先进的边缘检测技术
2.4.1Marr-Hildreth边缘检测器
高斯拉普拉斯(LoG):
G
(
x
,
y
)
=
e
−
x
2
+
y
2
2
σ
2
G(x,y)=e^{-\frac{x^{2}+y^{2}}{2\sigma ^{2}}}
G(x,y)=e−2σ2x2+y2
▽
2
G
(
x
,
y
)
=
[
x
2
+
y
2
−
2
σ
2
σ
4
]
e
−
x
2
+
y
2
2
σ
2
\triangledown ^{2}G(x,y)=[\frac{x^{2}+y^{2}-2\sigma ^{2}}{\sigma ^{4}}]e^{-\frac{x^{2}+y^{2}}{2\sigma ^{2}}}
▽2G(x,y)=[σ4x2+y2−2σ2]e−2σ2x2+y2
g
(
x
,
y
)
=
▽
2
[
G
(
x
,
y
)
★
f
(
x
,
y
)
]
g(x,y)=\triangledown ^{2}[G(x,y)\bigstar f(x,y)]
g(x,y)=▽2[G(x,y)★f(x,y)]
Marr-Hildreth边缘检测算法:
1、用对式
G
(
x
,
y
)
=
e
−
x
2
+
y
2
2
σ
2
G(x,y)=e^{-\frac{x^{2}+y^{2}}{2\sigma ^{2}}}
G(x,y)=e−2σ2x2+y2取样得到的一个n×n的高斯低通滤波器对输入图像滤波;
2、计算由第一步得到的图像的拉普拉斯;
3、找到步骤2所的图像的零交叉。
高斯差分(DoG):
D
o
G
(
x
,
y
)
=
1
2
π
σ
1
2
e
−
x
2
+
y
2
2
σ
1
2
−
1
2
π
σ
2
2
e
−
x
2
+
y
2
2
σ
2
2
,其中
σ
1
>
σ
2
DoG(x,y)=\frac{1}{2\pi \sigma _{1}^{2}}e^{-\frac{x^{2}+y^{2}}{2\sigma _{1}^{2}}}-\frac{1}{2\pi \sigma _{2}^{2}}e^{-\frac{x^{2}+y^{2}}{2\sigma _{2}^{2}}},其中\sigma _{1}>\sigma _{2}
DoG(x,y)=2πσ121e−2σ12x2+y2−2πσ221e−2σ22x2+y2,其中σ1>σ2
2.4.2坎尼边缘检测器
坎尼检测算法基本步骤:
1、用一个高斯滤波器平滑输入图像;
2、计算梯度幅值图像和角度图像;
3、对梯度幅值图像应用非最大抑制;
4、用双阈值处理和连接分析来检测并连接边缘。
2.4.3Python实现Marr-Hildreth、Canny检测算法
代码:
import cv2
import numpy as np
import matplotlib.pyplot as plt
# 读入图像并灰度化
image = cv2.imread("JMU.jpg", cv2.IMREAD_UNCHANGED)
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
# 拉普拉斯卷积核
kernel_Laplacian_1 = np.array([
[0, 1, 0],
[1, -4, 1],
[0, 1, 0]])
kernel_Laplacian_2 = np.array([
[1, 1, 1],
[1, -8, 1],
[1, 1, 1]])
# 下面两个卷积核不具有旋转不变性
kernel_Laplacian_3 = np.array([
[2, -1, 2],
[-1, -4, -1],
[2, 1, 2]])
kernel_Laplacian_4 = np.array([
[-1, 2, -1],
[2, -4, 2],
[-1, 2, -1]])
# 5*5 LOG卷积模板
kernel_LoG = 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]])
# Canny边缘检测 k为高斯核大小,t1,t2为阈值大小
def Canny(image, k, t1, t2):
img = cv2.GaussianBlur(image, (k, k), 0)
canny = cv2.Canny(img, t1, t2)
return canny
# 卷积
output_1 = cv2.filter2D(image, -1, kernel_Laplacian_1)
output_2 = cv2.filter2D(image, -1, kernel_LoG)
output_3 = Canny(image, 3, 50, 150)
# 显示处理后的图像
plt.subplot(2, 3, 1), plt.title('Original Image')
plt.imshow(image, cmap='gray'), plt.axis('off')
plt.subplot(2, 3, 4), plt.title('Laplacian')
plt.imshow(output_1, cmap='gray'), plt.axis('off')
plt.subplot(2, 3, 5), plt.title('LOG')
plt.imshow(output_2, cmap='gray'), plt.axis('off')
plt.subplot(2, 3, 6), plt.title('Canny')
plt.imshow(output_3, cmap='gray'), plt.axis('off')
plt.show()
(代码参考链接:http://t.csdn.cn/uDPhA)
结果:
分析:
对比2.3.3中实验结果的图像,可以看出Laplacian对图像中的阶跃性边缘点定位准确,对噪声敏感。LOG图像检测出了大多数的主要边缘,并且滤掉一些无关特征。Canny图像检测到了弱边缘,不易受噪声的干扰,对比其他图像有了明显的改进。
2.5边缘连接和边界检测
由于噪声、照明等产生边缘间断,使得一组像素难以完整形成边缘。因此,在边缘检测算法后,使用连接过程将间断的边缘像素组合成完整边缘。
2.5.1局部处理
分析图像中每个边缘点(x,y)的一个邻域内的像素,根据某种准则将相似点进行连接,由满足该准则的像素连接形成边缘。通过边缘像素梯度算子的响应强度和方向来确定边缘像素的相似性。
边缘像素梯度算子的响应强度:如果
∣
▽
f
(
x
,
y
)
−
▽
f
(
x
0
,
y
0
)
∣
≤
E
\left | \triangledown f(x,y)-\triangledown f(x_{0},y_{0}) \right |\leq E
∣▽f(x,y)−▽f(x0,y0)∣≤E
则
(
x
,
y
)
(x,y)
(x,y)领域内坐标为
(
x
0
,
y
0
)
(x_{0},y_{0})
(x0,y0)的边缘要素在幅度上相似于
(
x
,
y
)
(x,y)
(x,y)的像素。
边缘像素梯度算子的方向:如果
∣
α
(
x
,
y
)
−
α
(
x
0
−
y
0
)
∣
<
A
,
α
(
x
,
y
)
=
a
r
c
t
a
n
(
G
y
G
x
)
\left | \alpha (x,y)-\alpha (x_{0}-y_{0})\right |< A,\alpha (x,y)=arctan(\frac{G_{y}}{G_{x}})
∣α(x,y)−α(x0−y0)∣<A,α(x,y)=arctan(GxGy)
则
(
x
,
y
)
(x,y)
(x,y)邻域内坐标为
(
x
0
,
y
0
)
(x_{0},y_{0})
(x0,y0)的边缘像素,在角度 上相似于
(
x
,
y
)
(x,y)
(x,y)的像素。
2.5.2使用霍夫变换的全局处理
Hough变换的基本思想:
对于边界上的n个点的点集,找出共线的点集和直线方程。对于任意两点的直线方程:
y
=
a
x
+
b
y=ax+b
y=ax+b,构造一个参数a,b的平面,从而有结论:xy平面上的任意一条直线
y
=
a
x
+
b
y=ax+b
y=ax+b,对应在参数a、b平面上都有一个点过xy平一个点(x,y)的所有直线,构成参数ab平面上的一条直线。
如果点
(
x
1
,
y
1
)
与点
(
x
2
,
y
2
)
(x_{1},y_{1})与点(x_{2},y_{2})
(x1,y1)与点(x2,y2)共线,那么这两点在参数ab平面上的直线将有一个交点,具有相同的a和b。在参数ab平面上相交直线最多的点,对应的xy平面上的直线就是我们的解。如下图所示:
Hough变换算法实现:
由于垂直直线a为无穷大,我们改用极坐标形式:
x
c
o
s
θ
+
y
s
i
n
θ
=
ρ
。
xcos\theta + ysin\theta = \rho 。
xcosθ+ysinθ=ρ。参数平面为
θ
,
\theta ,
θ,
ρ
,
\rho ,
ρ,对应不是直线而是正弦曲线。使用交点累加器,或交点统计直方图,找出相交线段最多的参数空间的点,然后找出该点对应的xy平面的直线线段。
三、阈值处理
3.1基础知识
阈值处理后的图像
g
(
x
,
y
)
g(x,y)
g(x,y)定义为
g
(
x
,
y
)
=
{
1
f
(
x
,
y
)
>
T
0
f
(
x
,
y
)
≤
T
g(x,y)=\left\{\begin{matrix}1 &f(x,y)> T \\ 0&f(x,y)\leq T\end{matrix}\right.
g(x,y)={10f(x,y)>Tf(x,y)≤T
当T时一个适用于整个图像的常数时,上式给出的处理称为全局阈值处理;当T值在一幅图像上改变时,我们使用可变阈值处理这一术语;术语局部阈值处理或区域阈值处理有时用于表示可变阈值处理,这时图像中任何点(x,y)处的T值取决于(x,y)的领域的特性;若T取决于空间坐标(x,y)本身,则可变阈值处理通常称为动态阈值处理或自适应阈值处理。
3.2基本的全局阈值处理
3.2.1计算基本全局阈值算法:
1、选择一个T的初始估计值;
2、用T分割图像,生成两组像素:
G
1
G_{1}
G1由所有灰度值大于T的像素组成,而
G
2
G_{2}
G2由所有灰度值小于或等于T的像素组成;
3、对区域
G
1
和
G
2
G_{1}和G_{2}
G1和G2中的所有像素计算平均灰度值
μ
1
和
μ
2
\mu _{1}和\mu _{2}
μ1和μ2;
4、计算新的阈值
T
=
1
2
(
μ
1
+
μ
2
)
T=\frac{1}{2}(\mu _{1}+\mu _{2})
T=21(μ1+μ2);
5、重复步骤2到4,直到逐次迭代所得的T值之差小于事先定义的参数
T
0
T_{0}
T0。
3.2.2Python实现基本的全局阈值处理
代码:
import numpy as np
import cv2
import matplotlib.pyplot as plt
img = cv2.imread('JMU1.jpg', 0)
# 精度
eps = 1
iry = np.array(img)
r, c = img.shape
avg = 0
for i in range(r):
for j in range(c):
avg += iry[i][j]
T =int(avg/(r*c))
while 1:
G1, G2, cnt1, cnt2 = 0, 0, 0, 0
for i in range(r):
for j in range(c):
if iry[i][j] >= T: G1 += iry[i][j]; cnt1 += 1
else: G2 += iry[i][j]; cnt2 += 1
u1 = int(G1 / cnt1)
u2 = int(G2 / cnt2)
T2 = (u1 + u2) / 2
dis = abs(T2 - T)
if(dis <= eps): break
else :T = T2
new_img = np.zeros((r, c),np.uint8)
for i in range(r):
for j in range(c):
if iry[i][j] >= T: new_img[i][j] = 255
else: new_img[i][j] = 0
cv2.imshow('1', img)
cv2.imshow('2', new_img)
cv2.waitKey()
plt.hist(img.ravel(),256,[0,256])
plt.show()
(代码参考链接:http://t.csdn.cn/lt05n)
结果:
分析:
直方图显示有明显的波谷,从直方图清晰的模式分离中,物体和背景的分割效果很好。
3.3用Otsu方法的最佳全局阈值处理
3.3.1Otsu算法
1、计算输入图像的归一化直方图。使用
p
i
,
i
=
0
,
1
,
2
,
.
.
.
,
L
−
1
p_{i},i=0,1,2,...,L-1
pi,i=0,1,2,...,L−1表示该直方图的各个分量;
2、用式
P
1
(
k
)
=
∑
i
=
0
k
P
i
P_{1}(k)=\sum_{i=0}^{k}P_{i}
P1(k)=∑i=0kPi,对于
k
=
0
,
1
,
2
,
.
.
.
,
L
−
1
k=0,1,2,...,L-1
k=0,1,2,...,L−1,计算累积和
P
1
(
k
)
P_{1}(k)
P1(k);
3、用式
m
(
k
)
=
∑
i
=
0
k
i
P
i
m(k)=\sum_{i=0}^{k}iP_{i}
m(k)=∑i=0kiPi,对于
k
=
0
,
1
,
2
,
.
.
.
,
L
−
1
k=0,1,2,...,L-1
k=0,1,2,...,L−1,计算累积均值
m
(
k
)
m(k)
m(k);
4、用式
m
G
=
∑
i
=
0
L
−
1
i
P
i
m_{G}=\sum_{i=0}^{L-1}iP_{i}
mG=∑i=0L−1iPi计算全局灰度均值
m
G
m_{G}
mG;
5、用式
σ
B
2
(
k
)
=
[
m
G
P
1
(
k
)
−
m
(
k
)
]
2
P
1
(
k
)
[
1
−
P
1
(
k
)
]
\sigma _{B}^{2}(k)=\frac{[m_{G}P_{1}(k)-m(k)]^{2}}{P_{1}(k)[1-P_{1}(k)]}
σB2(k)=P1(k)[1−P1(k)][mGP1(k)−m(k)]2,对于
k
=
0
,
1
,
2
,
.
.
.
,
L
−
1
k=0,1,2,...,L-1
k=0,1,2,...,L−1,计算类间方差
σ
B
2
(
k
)
\sigma _{B}^{2}(k)
σB2(k);
6、得到Orsu阈值
k
∗
k^{\ast }
k∗,即使得
σ
B
2
(
k
)
\sigma _{B}^{2}(k)
σB2(k)最大的k值。如果最大值不唯一,用相应检测到的各个最大值k的平均得到
k
∗
k^{\ast }
k∗。
7、在
k
=
k
∗
k=k^{\ast }
k=k∗处计算式
η
(
k
)
=
σ
B
2
(
k
)
σ
G
2
\eta (k)=\frac{\sigma _{B}^{2}(k)}{\sigma _{G}^{2}}
η(k)=σG2σB2(k),得到可分性测度
η
∗
\eta ^{\ast }
η∗。
3.3.2Python实现Otsu
代码:
import cv2
import os
import numpy as np
import pandas as pd
from tqdm import tqdm
import matplotlib.pyplot as plt
### 首先将图片转化为灰度图像
image = cv2.imread("JMU1.jpg")
def rgb2gray(image):
h = image.shape[0]
w = image.shape[1]
grayimage = np.zeros((h,w),np.uint8)
for i in tqdm(range(h)):
for j in range(w):
grayimage [i,j] = 0.144*image[i,j,0]+0.587*image[i,j,1]+0.299*image[i,j,1]
return grayimage
### 大津法
def otsu(image):
### 高和宽
h = image.shape[0]
w = image.shape[1]
### 求总像素
m = h*w
otsuimg = np.zeros((h, w), np.uint8)
##初始阈值
initial_threshold = 0
### 最终阈值
final_threshold = 0
# 初始化各灰度级个数统计参数
histogram = np.zeros(256, np.int32)
# 初始化各灰度级占图像中的分布的统计参数
probability = np.zeros(256, np.float32)
### 各个灰度级的个数统计
for i in tqdm(range(h)):
for j in range(w):
s = image[i,j]
histogram[s] = histogram[s] +1
### 各灰度级占图像中的分布的统计参数
for i in tqdm(range(256)):
probability[i] = histogram[i]/m
for i in tqdm(range(255)):
w0 = w1 = 0 ## 前景和背景的灰度数
fgs = bgs = 0 # 定义前景像素点灰度级总和背景像素点灰度级总和
for j in range(256):
if j <= i: # 当前i为分割阈值
w0 += probability[j] # 前景像素点占整幅图像的比例累加
fgs += j * probability[j]
else:
w1 += probability[j] # 背景像素点占整幅图像的比例累加
bgs += j * probability[j]
u0 = fgs / w0 # 前景像素点的平均灰度
u1 = bgs / w1 # 背景像素点的平均灰度
G = w0*w1*(u0-u1)**2
if G >= initial_threshold:
initial_threshold = G
final_threshold = i
print(final_threshold)
for i in range(h):
for j in range(w):
if image[i, j] > final_threshold:
otsuimg[i, j] = 255
else:
otsuimg[i, j] = 0
return otsuimg
grayimage = rgb2gray(image)
otsuimage = otsu(grayimage)
cv2.imshow("grayimage",grayimage)
cv2.imshow("otsuimage",otsuimage)
cv2.waitKey()
plt.hist(image.ravel(),256,[0,256])
plt.show()
(代码参考链接:http://t.csdn.cn/fNjGK)
结果:
分析:
对比3.2中基本全局算法得到的分割结果,低值中的细节更突出,更有效地将目标与背景分割。
3.4自适应阈值处理
在前面的部分我们使用是全局阈值,整幅图像采用同一个数作为阈值。当时这种方法并不适应与所有情况,尤其是当同一幅图像上的不同部分的具有不同亮度时,这种情况下我们需要采用自适应阈值。此时的阈值是根据图像上的每一个小区域计算与其对应的阈值。因此在同一幅图像上的不同区域采用的是不同的阈值,从而使我们能在亮度不同的情况下得到更好的结果。
Python实现自适应阈值处理
代码:
import cv2
import matplotlib.pyplot as plt
img = cv2.imread('cat.jpg', cv2.IMREAD_GRAYSCALE) # 读取图像,转化为单通道灰度图像
cv2.imshow('img', img)
img2 = cv2.adaptiveThreshold(img, 255, cv2.ADAPTIVE_THRESH_MEAN_C,
cv2.THRESH_BINARY, 5, 10) # 阈值处理
cv2.imshow('img2', img2)
cv2.waitKey(0)
plt.hist(img.ravel(),256,[0,256])
plt.show()
使用自适应阈值处理结果:
使用Otsu算法的结果:
分析:
从结果中两种处理得到的图像及直方图对比可以看出,对于色彩均衡的图像来说,直接使用一个阈值就能完成对图像的阈值化处理。但是,有时图像的色彩是不均衡的,如果只用一个阈值,就无法得到清晰有效的阈值分割结果图像。例如图像左侧花朵的区域灰度高于背景灰度,Otsu放大的最佳全局阈值处理的结果更好,但是图像右侧猫的区域效果则不是很理想。因此当同一幅图像上的不同部分的具有不同亮度时,需要采用自适应阈值。
四、基于区域的分割
4.1基本概念
将区域R划分为若干个子区域
R
1
,
R
2
,
.
.
.
,
R
n
R_{1},R_{2},...,R_{n}
R1,R2,...,Rn,这些子区域要满足五个条件:
1、完备性:
⋃
i
=
1
n
=
R
\bigcup_{i=1}^{n}=R
⋃i=1n=R
2、连通性:每个
R
i
都是一个连通区域
R_{i}都是一个连通区域
Ri都是一个连通区域
3、独立性:对于任意
i
≠
j
,
R
i
∩
R
j
=
∅
i≠j,R_{i}\cap R_{j}=\varnothing
i=j,Ri∩Rj=∅
4、单一性:每个区域内的灰度级相等
5、互斥性:任意两个区域的灰度级不等
4.2区域增长
区域增长的算法实现:
1、根据图像的不同应用选择一个或一组种子,它或者是最亮或最暗的点,或者是位于点簇中心的点;
2、选择一个描述符(条件);
3、从该种子开始向外扩张,首先把种子像素加入结果集合,然后不断将与集合中各个像素连通、且满足描述符的像素加入集合;
4、上一过程进行到不再有满足条件的新结点加入集合为止。
Python实现区域增长的分割
代码:
import numpy as np
import cv2
class Point(object):
def __init__(self, x, y):
self.x = x
self.y = y
def getX(self):
return self.x
def getY(self):
return self.y
def getGrayDiff(img, currentPoint, tmpPoint):
return abs(int(img[currentPoint.x, currentPoint.y]) - int(img[tmpPoint.x, tmpPoint.y]))
def selectConnects(p):
if p != 0:
connects = [Point(-1, -1), Point(0, -1), Point(1, -1), Point(1, 0), Point(1, 1), \
Point(0, 1), Point(-1, 1), Point(-1, 0)]
else:
connects = [Point(0, -1), Point(1, 0), Point(0, 1), Point(-1, 0)]
return connects
def regionGrow(img, seeds, thresh, p=1):
height, weight = img.shape
seedMark = np.zeros(img.shape)
seedList = []
for seed in seeds:
seedList.append(seed)
label = 1
connects = selectConnects(p)
while (len(seedList) > 0):
currentPoint = seedList.pop(0)
seedMark[currentPoint.x, currentPoint.y] = label
for i in range(8):
tmpX = currentPoint.x + connects[i].x
tmpY = currentPoint.y + connects[i].y
if tmpX < 0 or tmpY < 0 or tmpX >= height or tmpY >= weight:
continue
grayDiff = getGrayDiff(img, currentPoint, Point(tmpX, tmpY))
if grayDiff < thresh and seedMark[tmpX, tmpY] == 0:
seedMark[tmpX, tmpY] = label
seedList.append(Point(tmpX, tmpY))
return seedMark
img = cv2.imread('cat.jpg', 0)
seeds = [Point(10, 10), Point(82, 150), Point(20, 300)]
binaryImg = regionGrow(img, seeds, 10)
cv2.imshow('1', img)
cv2.imshow('2', binaryImg)
cv2.waitKey(0)
(代码参考链接:http://t.csdn.cn/RNJRQ)
结果:
分析:
区域生长的基本思想是将具有相似性质的像素集合起来构成区域。具体先对每个需要分割的区域找一个种子像素作为生长起点,然后将种子像素和周围邻域中与种子像素有相同或相似性质的像素合并到种子像素所在的区域中。将这些新像素当作新的种子继续上面的过程,直到没有满足条件的像素可被包括进来。
该实验分割出了图像中花朵的区域。
4.3区域分裂与聚合
区域分裂与聚合的算法实现:
1、对图像中灰度级不同的区域,均分为四个子区域;
2、如果相邻的子区域所有像素的灰度级相同,则将其合并;
3、反复进行上两步操作,直至不再有新的分裂与合并为止。
五、总结
图像分割是一种基本的预处理步骤,本章学习了多种图像分割技术,在应用中选择何种技术需要考虑问题的特殊性。