5.1 简介
该部分的学习内容是对经典的阈值分割算法进行回顾,图像阈值化分割是一种传统的最常用的图像分割方法,因其实现简单、计算量小、性能较稳定而成为图像分割中最基本和应用最广泛的分割技术。它特别适用于目标和背景占据不同灰度级范围的图像。它不仅可以极大的压缩数据量,而且也大大简化了分析和处理步骤,因此在很多情况下,是进行图像分析、特征提取与模式识别之前的必要的图像预处理过程。图像阈值化的目的是要按照灰度级,对像素集合进行一个划分,得到的每个子集形成一个与现实景物相对应的区域,各个区域内部具有一致的属性,而相邻区域不具有这种一致属性。这样的划分可以通过从灰度级出发选取一个或多个阈值来实现。
5.2 算法理论介绍
5.2.1 最大类间方差法(大津阈值法)
大津法(OTSU)是一种确定图像二值化分割阈值的算法,由日本学者大津于1979年提出。从大津法的原理上来讲,该方法又称作最大类间方差法,因为按照大津法求得的阈值进行图像二值化分割后,前景与背景图像的类间方差最大。
它被认为是图像分割中阈值选取的最佳算法,计算简单,不受图像亮度和对比度的影响,因此在数字图像处理上得到了广泛的应用。它是按图像的灰度特性,将图像分成背景和前景两部分。因方差是灰度分布均匀性的一种度量,背景和前景之间的类间方差越大,说明构成图像的两部分的差别越大,当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。
应用: 是求图像全局阈值的最佳方法,应用不言而喻,适用于大部分需要求图像全局阈值的场合。
优点: 计算简单快速,不受图像亮度和对比度的影响。
缺点: 对图像噪声敏感;只能针对单一目标分割;当目标和背景大小比例悬殊、类间方差函数可能呈现双峰或者多峰,这个时候效果不好。
原理非常简单,涉及的知识点就是均值、方差等概念和一些公式推导。为了便于理解,我们从目的入手,反推一下这著名的OTSU算法。
- 求类间方差:
OTSU算法的假设是存在阈值TH将图像所有像素分为两类C1(小于TH)和C2(大于TH),则这两类像素各自的均值就为m1、m2,图像全局均值为mG。同时像素被分为C1和C2类的概率分别为p1、p2。因此就有:
p
1
∗
m
1
+
p
2
∗
m
2
=
m
G
p_1*m_1+p_2*m_2=mG
p1∗m1+p2∗m2=mG
p
1
+
p
2
=
1
p_1+p_2=1
p1+p2=1
类间方差表达式:
σ
2
=
p
1
(
m
1
−
m
G
)
2
+
p
2
(
m
2
−
m
G
)
2
\sigma^2=p_1(m_1-mG)^2+p_2(m_2-mG)^2
σ2=p1(m1−mG)2+p2(m2−mG)2
将上式子化简,将
p
1
∗
m
1
+
p
2
∗
m
2
=
m
G
p_1*m_1+p_2*m_2=mG
p1∗m1+p2∗m2=mG代入
σ
2
=
p
1
(
m
1
−
m
G
)
2
+
p
2
(
m
2
−
m
G
)
2
\sigma^2=p_1(m_1-mG)^2+p_2(m_2-mG)^2
σ2=p1(m1−mG)2+p2(m2−mG)2,得
σ
2
=
p
1
p
2
(
m
1
−
m
2
)
2
\sigma^2=p_1 p_2(m_1-m_2)^2
σ2=p1p2(m1−m2)2
将上式最大化的灰度级k就是OTSU阈值。,其中:
p
1
=
∑
i
=
0
k
P
i
p_1=\sum_{i=0}^k P_i
p1=i=0∑kPi
m
1
=
1
/
p
1
∗
∑
i
=
0
k
i
p
i
m_1=1/p_1*\sum_{i=0}^k ip_i
m1=1/p1∗i=0∑kipi
m
2
=
1
/
p
2
∗
∑
i
=
k
+
1
L
−
1
i
p
i
m_2=1/p_2*\sum_{i=k+1}^{L-1} ip_i
m2=1/p2∗i=k+1∑L−1ipi
σ
2
=
p
1
p
2
(
m
1
−
m
2
)
2
\sigma^2=p_1 p_2(m_1-m_2)^2
σ2=p1p2(m1−m2)2还可以做变形:
首先灰度级K的累加均值m和图像全局均值mG分别:
m
=
∑
i
=
0
k
i
p
i
m=\sum_{i=0}^k ip_i
m=i=0∑kipi
m
G
=
∑
i
=
0
L
−
1
i
p
i
mG=\sum_{i=0}^{L-1}ip_i
mG=i=0∑L−1ipi
再看
m
1
=
1
/
p
1
∗
∑
i
=
0
k
i
p
i
m_1=1/p_1*\sum_{i=0}^k ip_i
m1=1/p1∗∑i=0kipi,m1、m2的变为:
m
1
=
1
/
p
1
∗
m
m_1=1/p_1*m
m1=1/p1∗m
m
2
=
1
/
p
2
∗
(
m
G
−
m
)
m_2=1/p2*(mG-m)
m2=1/p2∗(mG−m)
最终类间方差公式:
σ
=
(
m
G
∗
p
1
−
m
)
2
p
1
(
1
−
p
1
)
\sigma=\frac{(mG*p_1-m)^2}{p_1(1-p_1)}
σ=p1(1−p1)(mG∗p1−m)2
利用最后算到的p1,m,mG和上式:最大化的灰度级k就是OTSU阈值
- 分割
这个分割就是二值化。- THRESH_BINARY d s t ( x , y ) = { m a x v a l i f s r c ( x , y ) > t h r e s h 0 o h t e r w i s e dst(x,y)=\begin{cases}maxval &if src(x,y)>thresh\\0 &ohterwise\end{cases} dst(x,y)={maxval0ifsrc(x,y)>threshohterwise
- THRESH_BINARY_INV d s t ( x , y ) = { 0 i f s r c ( x , y ) > t h r e s h m a x v a l o t h e r w i s e dst(x,y)=\begin{cases}0 &if src(x,y)>thresh\\ maxval &otherwise\end{cases} dst(x,y)={0maxvalifsrc(x,y)>threshotherwise
- THRESH_TRUNC d s t ( x , y ) = { t h r e s h o l d i f s r c ( x , y ) > t h r e s h 0 o t h e r w i s e dst(x,y)=\begin{cases}threshold &if src(x,y)>thresh\\0 &otherwise\end{cases} dst(x,y)={threshold0ifsrc(x,y)>threshotherwise
- THRESH_TOZERO d s t ( x , y ) = { s r c ( x , y ) i f s r c ( x , y ) > t h r e s h 0 o t h e r w i s e dst(x,y)=\begin{cases}src(x,y) &if src(x,y)>thresh\\0 &otherwise\end{cases} dst(x,y)={src(x,y)0ifsrc(x,y)>threshotherwise
- THRESH_TOZERO_INV d s t ( x , y ) = { 0 s r c ( x , y ) > t h r e s h s r c ( x , y ) o t h e r w i s e dst(x,y)=\begin{cases}0 &src(x,y)>thresh\\src(x,y) &otherwise\end{cases} dst(x,y)={0src(x,y)src(x,y)>threshotherwise
5.2.2 自适应阈值
根据图像不同区域亮度分布,计算其局部阈值,所以对于图像不同区域,能够自适应计算不同的阈值,因此被称为自适应阈值法。(其实就是局部阈值法)
- 如何确定局部阈值呢?可以计算某个邻域(局部)的均值、中值、高斯加权平均(高斯滤波)来确定阈值。值得说明的是:如果用局部的均值作为局部的阈值,就是常说的移动平均法。
5.3 代码实践
5.3.1 简单阈值
对于每个像素,应用相同的阈值。如果像素值小于阈值,则将其设置为0,否则将其设置为最大值。
- 函数cv.threshold用于应用阈值。第一个参数是源图像,它应该是灰度图像。第二个参数是阈值,用于对像素值进行分类。第三个参数是分配给超过阈值的像素值的最大值。OpenCV提供了不同类型的阈值,这由函数的第四个参数给出。通过使用cv.THRESH_BINARY类型。所有简单的阈值类型为:
- cv.THRESH_BINARY
- cv.THRESH_BINARY_INV
- cv.THRESH_TRUNC
- cv.THRESH_TOZERO
- cv.THRESH_TOZERO_INV- 该方法返回两个输出。第一个是使用的阈值,第二个输出是阈值后的图像。
import cv2
import numpy as np
from matplotlib import pyplot as plt
img = cv.imread('C:/Users/83769/photo.jpg',0)
ret,thresh1 = cv2.threshold(img,127,255,cv.THRESH_BINARY)
ret,thresh2 = cv2.threshold(img,127,255,cv.THRESH_BINARY_INV)
ret,thresh3 = cv2.threshold(img,127,255,cv.THRESH_TRUNC)
ret,thresh4 = cv2.threshold(img,127,255,cv.THRESH_TOZERO)
ret,thresh5 = cv2.threshold(img,127,255,cv.THRESH_TOZERO_INV)
titles = ['Original Image','BINARY','BINARY_INV','TRUNC','TOZERO','TOZERO_INV']
images = [img, thresh1, thresh2, thresh3, thresh4, thresh5]
for i in range(6):
plt.subplot(2,3,i+1),plt.imshow(images[i],'gray')
plt.title(titles[i])
plt.xticks([]),plt.yticks([])
plt.show()
5.3.2 自适应阈值
根据图像上的每一个小区域计算与其对应的阀值。因此在同一幅图像上的不同区域采用的是不同的阀值,从而使我们能在亮度不同的情况下得到更好的结果。
这种方法需要我们指定三个参数,返回值只有一个。
Adaptive Method 指定计算阀值的方法
- cv2.ADAPTIVE_THRESH_MEAN_C:阀值取自相邻区域的平均值
- cv2.ADAPTIVE_THRESH_GAUSSIAN_C:阀值取自相邻区域的加权和,权重为一个高斯窗口
- Block Size 邻域大小(用来计算阀值的区域大小)
- C这就是一个常数,阀值就等于的平均值或者加权平均值减去这个常数
#中值滤波
img = cv2.medianBlur(img,5)
ret , th1 = cv2.threshold(img,127,255,cv2.THRESH_BINARY)
# 11为block size,2为C值
th2 = cv2.adaptiveThreshold(img,255,cv2.ADAPTIVE_THRESH_MEAN_C , cv2.THRESH_BINARY,11,2 )
th3 = cv2.adaptiveThreshold(img,255,cv2.ADAPTIVE_THRESH_GAUSSIAN_C , cv2.THRESH_BINARY,11,2)
titles = ['original image' , 'global thresholding (v=127)','Adaptive mean thresholding',
'adaptive gaussian thresholding']
images = [img,th1,th2,th3]
for i in range(4):
plt.subplot(2,2,i+1),plt.imshow(images[i],'gray')
plt.title(titles[i])
plt.xticks([]),plt.yticks([])
plt.show()
5.3.3 Otsu的二值化
在全局阈值化中,我们使用任意选择的值作为阈值。相反,Otsu的方法避免了必须选择一个值并自动确定它的情况。
考虑仅具有两个不同图像值的图像(双峰图像),其中直方图将仅包含两个峰。一个好的阈值应该在这两个值的中间。类似地,Otsu的方法从图像直方图中确定最佳全局阈值。
为此,使用了cv.threshold作为附加标志传递。阈值可以任意选择。然后,算法找到最佳阈值,该阈值作为第一输出返回。那么经过Otsu’s得到的那个阈值就是函数cv2.threshold的第一个参数了。
- 因为Otsu’s方法会产生一个阈值,那么函数cv2.threshold的的第二个参数(设置阈值)就是0了,并且在cv2.threshold的方法参数中还得加上语句cv2.THRESH_OTSU。
ret1,th1=cv2.threshold(img,127,255,cv2.THRESH_BINARY)
ret2,th2=cv2.threshold(img,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
#(5,5)为高斯核的大小,0为标准差
blur= cv2.GaussianBlur(img,(5,5),0)
#阀值一定要设为0
ret3,th3=cv2.threshold(blur,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
images=[img,0,th1,
img,0,th2,
img,0,th3]
titles =['original noisy image','histogram','global thresholding(v=127)',
'original noisy image','histogram',"otsu's thresholding",
'gaussian giltered image','histogram',"otus's thresholding"]
#这里使用了pyplot中画直方图的方法,plt.hist要注意的是他的参数是一维数组
#所以这里使用了(numpy)ravel方法,将多维数组转换成一维,也可以使用flatten方法
for i in range(3):
plt.subplot(3,3,i*3+1),plt.imshow(images[i*3],'gray')
plt.title(titles[i*3]),plt.xticks([]),plt.yticks([])
plt.subplot(3,3,i*3+2),plt.hist(images[i*3].ravel(),256)
plt.title(titles[i*3+1]),plt.xticks([]),plt.yticks([])
plt.subplot(3,3,i*3+3),plt.imshow(images[i*3+2],'gray')
plt.title(titles[i*3+2]),plt.xticks([]),plt.yticks([])
plt.show()