欢迎关注 『youcans 的 OpenCV 例程 200 篇』 系列,持续更新中
OpenCV 例程200篇 总目录-202205更新
【youcans 的 OpenCV 例程200篇】165.多阈值 OTSU 处理方法
3.5 多阈值处理方法
OTSU 方法使用最大化类间方差(intra-class variance)作为评价准则,基于对图像直方图的计算,可以给出类间最优分离的最优阈值。
OTSU 方法可以扩展到任意数量的阈值。假设有 K 个分类
c
1
,
c
2
,
.
.
.
,
c
k
c_1, c_2,...,c_k
c1,c2,...,ck 时,可以定义类间方差为:
σ
B
2
=
∑
k
=
1
K
P
k
(
m
k
−
m
G
)
2
\sigma^2_B = \sum^K_{k=1} P_k(m_k - m_G)^2
σB2=k=1∑KPk(mk−mG)2
常用地,考虑由 3个灰度区间组成的 3个类,可以由 2个阈值分割,定义类间方差为:
σ
B
2
=
P
1
(
m
1
−
m
G
)
2
+
P
2
(
m
2
−
m
G
)
2
+
P
3
(
m
3
−
m
G
)
2
P
1
+
P
2
+
P
3
=
1
P
1
∗
m
1
+
P
2
∗
m
2
+
P
3
∗
m
3
=
m
G
P
1
=
∑
i
=
0
k
1
p
i
,
P
2
=
∑
i
=
k
1
+
1
k
2
p
i
,
P
3
=
∑
i
=
k
2
+
1
L
−
1
p
i
m
1
=
∑
i
=
0
k
1
i
∗
p
i
/
P
1
,
m
2
=
∑
i
=
k
1
+
1
k
2
i
∗
p
i
/
P
2
,
m
3
=
∑
i
=
k
2
+
1
L
−
1
i
∗
p
i
/
P
3
,
\sigma^2_B = P_1(m_1 - m_G)^2 + P_2(m_2 - m_G)^2 + P_3(m_3 - m_G)^2 \\ P_1 + P_2 + P_3 = 1\\ P_1*m_1 + P_2*m_2 + P_3*m_3 = m_G \\ P_1 = \sum^{k_1}_{i=0} p_i ,\ P_2 = \sum^{k_2}_{i=k_1+1} p_i ,\ P_3 = \sum^{L-1}_{i=k_2+1} p_i \\ m_1 = \sum^{k_1}_{i=0} i*p_i /P_1 ,\ m_2 = \sum^{k_2}_{i=k_1+1} i*p_i /P_2 ,\ m_3 = \sum^{L-1}_{i=k_2+1} i*p_i /P_3 ,\\
σB2=P1(m1−mG)2+P2(m2−mG)2+P3(m3−mG)2P1+P2+P3=1P1∗m1+P2∗m2+P3∗m3=mGP1=i=0∑k1pi, P2=i=k1+1∑k2pi, P3=i=k2+1∑L−1pim1=i=0∑k1i∗pi/P1, m2=i=k1+1∑k2i∗pi/P2, m3=i=k2+1∑L−1i∗pi/P3,
使类间方差 ICV 最大化的灰度值
k
1
,
k
2
k_1, k_2
k1,k2 就是最优阈值。
阈值处理后的图像由下式给出:
g
(
x
,
y
)
=
{
a
,
f
(
x
,
y
)
≤
k
1
∗
b
,
k
1
∗
≤
f
(
x
,
y
)
≤
k
2
∗
c
,
f
(
x
,
y
)
≥
k
2
∗
g(x,y) = \begin{cases} a, & f(x,y) \le k_1^* \\ b, & k_1^* \le f(x,y) \le k_2^* \\ c, & f(x,y) \ge k_2^* \\ \end{cases}
g(x,y)=⎩⎪⎨⎪⎧a,b,c,f(x,y)≤k1∗k1∗≤f(x,y)≤k2∗f(x,y)≥k2∗
需要指出的是,双阈值处理涉及两重循环,算法耗时较长,但仍是基于直方图进行运算的,因此用时与图像无关。对于双重循环可以进行优化,就不在本节介绍了。
如果需要使用更多个阈值进行分割,通常不再用循环遍历的方法,而是采用聚类或启发式方法来获得分割阈值。
例程 11.22:阈值处理之多阈值 OTSU
# 11.22 阈值处理之多阈值 OTSU
def doubleThreshold(img):
histCV = cv2.calcHist([img], [0], None, [256], [0, 256]) # 灰度直方图
grayScale = np.arange(0, 256, 1) # 灰度级 [0,255]
totalPixels = img.shape[0] * img.shape[1] # 像素总数
totalGray = np.dot(histCV[:,0], grayScale) # 内积, 总和灰度值
mG = totalGray / totalPixels # 平均灰度,meanGray
varG = sum(((i-mG)**2 * histCV[i,0]/totalPixels) for i in range(256))
T1, T2, varMax = 1, 2, 0.0
# minGary, maxGray = np.min(img), np.max(img) # 最小灰度,最大灰度
for k1 in range(1, 254): # k1: [1,253], 1<=k1<k2<=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
# n2 = sum(histCV[k1+1:k2,0]) # C2 像素数量
# s2 = sum( (i * histCV[i,0]) for i in range(k1+1,k2) )
# P2 = n2 / totalPixels # C2 像素数占比
# m2 = (s2/n2) if n2>0 else 0 # C2 平均灰度
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
img = cv2.imread("../images/Fig1043a.tif", flags=0)
# img = cv2.imread("../images/Fig1045a.tif", flags=0)
histCV = cv2.calcHist([img], [0], None, [256], [0, 256]) # 灰度直方图
T1, T2, epsT = doubleThreshold(img)
print("T1={}, T2={}, esp={:.4f}".format(T1, T2, epsT))
binary = img.copy()
binary[binary<T1] = 0
binary[binary>T2] = 255
ret, imgOtsu = cv2.threshold(img, 127, 255, cv2.THRESH_OTSU) # OTSU 阈值分割
ret1, binary1 = cv2.threshold(img, T1, 255, cv2.THRESH_TOZERO) # 小于阈值置 0,大于阈值不变
ret2, binary2 = cv2.threshold(img, T2, 255, cv2.THRESH_TOZERO)
plt.figure(figsize=(9, 6))
plt.subplot(231), plt.axis('off'), plt.title("Origin"), plt.imshow(img, 'gray')
plt.subplot(232,yticks=[]), plt.axis([0,255,0,np.max(histCV)])
plt.bar(range(256), histCV[:,0]), plt.title("Gray Hist")
plt.subplot(233), plt.title("OTSU binary(T={})".format(round(ret))), plt.axis('off')
plt.imshow(imgOtsu, 'gray')
plt.subplot(234), plt.axis('off'), plt.title("Threshold(T={})".format(T1))
plt.imshow(binary1, 'gray')
plt.subplot(235), plt.axis('off'), plt.title("Threshold(T={})".format(T2))
plt.imshow(binary2, 'gray')
plt.subplot(236), plt.axis('off'), plt.title("DoubleT({},{})".format(T1,T2))
plt.imshow(binary, 'gray')
plt.show()
运行结果:
Fig1043a.tif:T1=35, T2=101, esp=0.8733
Fig1045a.tif:T1=81, T2=177, esp=0.9540
(本节完)
版权声明:
youcans@xupt 原创作品,转载必须标注原文链接:(https://blog.csdn.net/youcans/article/details/124400439)
Copyright 2022 youcans, XUPT
Crated:2022-4-25
欢迎关注 『youcans 的 OpenCV 例程 200 篇』 系列,持续更新中
欢迎关注 『youcans 的 OpenCV学习课』 系列,持续更新中【youcans 的 OpenCV 例程200篇】158. 阈值处理之固定阈值法
【youcans 的 OpenCV 例程200篇】159. 图像分割之全局阈值处理
【youcans 的 OpenCV 例程200篇】160. 图像处理之OTSU 方法
【youcans 的 OpenCV 例程200篇】161. OTSU 阈值处理算法的实现
【youcans 的 OpenCV 例程200篇】162. 全局阈值处理改进方法
【youcans 的 OpenCV 例程200篇】163. 基于边缘信息改进全局阈值处理
【youcans 的 OpenCV 例程200篇】164.使用 Laplace 边缘信息改进全局阈值处理
【youcans 的 OpenCV 例程200篇】165.多阈值 OTSU 处理方法
【youcans 的 OpenCV 例程200篇】166.自适应阈值处理
【youcans 的 OpenCV 例程200篇】167.基于移动平均的可变阈值处理
更多内容,请见:
【OpenCV 例程200篇 总目录-202206更新】