阈值化分割（二）OTSU法-附Python实现

最新推荐文章于 2024-03-25 22:57:18 发布

xiaopihaierletian

最新推荐文章于 2024-03-25 22:57:18 发布

阅读量6.2k

点赞数 8

分类专栏： cv

本文链接：https://blog.csdn.net/xiaopihaierletian/article/details/76340380

版权

cv 专栏收录该内容

119 篇文章 1 订阅

订阅专栏

一、OTSU法（大津阈值分割法）介绍

OTSU算法是由日本学者OTSU于1979年提出的一种对图像进行二值化的高效算法，是一种自适应的阈值确定的方法，又称大津阈值分割法，是最小二乘法意义下的最优分割。

二、单阈值OTSU法

设图像包含L个灰度级，灰度值为i的像素点个数为Ni，像素总点数为：

N = N 0 + N 1 + \dots + N L - 1

则灰度值为i的点的概率为：

P i = N i N

根据期望公式，图像灰度的均值为：

μ T = \sum i = 0 L - 1 i P i

按图像的灰度特性，使用阈值T将图像分成目标c0和背景c1两类，则ω0(T)和ω1(T)分别表示阈值为T时，c0和c1发生的概率，即：

ω 0 (T) = \sum i = 0 T P i ω 1 (T) = 1 - ω 0 (T)

c0和c1的均值为：

μ 0 (T) = \sum T i = 0 i P i ω 0 ( T ) μ 1 (T) = μ T - \sum T i = 0 i P i ω 1 ( T )

σ^2(T)表示直方图中阈值为T的类间方差，定义为：

σ 2 B (T) = ω 0 (T) [μ 0 (T) - μ T] 2 + ω 1 (T) [μ 1 (T) - μ T] 2

最优阈值定义为类间方差最大时对应的T值，即：

σ 2 B (T *) = max 0 \leq T \leq L - 1 {σ 2 B (T)}

下面给出Python源代码。

#coding:utf-8
import cv2
import numpy as np
from matplotlib import pyplot as plt

image = cv2.imread("E:/python/cv/OTSU/test.jpg")
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

plt.subplot(131), plt.imshow(image, "gray")
plt.title("source image"), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.hist(image.ravel(), 256)
plt.title("Histogram"), plt.xticks([]), plt.yticks([])
ret1, th1 = cv2.threshold(gray, 0, 255, cv2.THRESH_OTSU)  #方法选择为THRESH_OTSU
plt.subplot(133), plt.imshow(th1, "gray")
plt.title("OTSU,threshold is " + str(ret1)), plt.xticks([]), plt.yticks([])
plt.show()
 
 2
3

结果如下所示。可以看到，使用OTSU法计算出来的阈值为165。
单阈值OTSU

三、多阈值OTSU法

将单阈值的OTSU推广到多阈值的图像分割中，假设将图像直方图分为m+1类，对应的阈值为T1，T2，···，Tm。则最大类间方差为：

σ 2 B (T * 1, T * 2, \dots, T * m) = max 0 \leq T 1 \leq T 2 \leq \dots \leq L - 1 {σ 2 B (T 1, T 2, \dots, T m)}

其中，

σ 2 B (T 1, T 2, \dots, T m) = \sum i = 0 m ω i (T 1, T 2, \dots, T m) [μ i (T 1, T 2, \dots, T m) - μ T] 2 μ i (T 1, T 2, \dots, T m) = \sum i = T i T i + 1 i P i ω i ( T 1 , T 2 , \dots , T m ) ω i (T 1, T 2, \dots, T m) = \sum i = T i T i + 1 P i μ T = \sum i = 0 L - 1 i P i

为求得最优阈值，需要使用穷举搜索，随着m增大，计算量骤增。若使用牛顿迭代等优化搜索方法，容易陷入局部最优解。

四、遗传算法解OTSU

遗传算法是一种基于自然选择和群体遗传机理的搜索算法。它模拟了自然选择和自然遗传过程中发生的繁殖、交配和突变现象，将每一个可能的解看作是群体的一个个体，并将每个个体编码，根据设定的目标函数对每个个体进行评价，给出一个适应度值。开始时随机产生一些个体，利用遗传算子产生新一代的个体，新个体继承上一代的优良性状，逐步向更优解进化。由于遗传算法在每一代同时搜索参数空间的不同区域，从而能够使找到全局最优解的可能性大大增加。遗传算法属于启发式算法，无限趋紧最优解并收敛。

那么怎么将图像分割问题抽象成遗传问题，即怎么将问题编码成基因串，如何构造适应度函数来度量每条基因的适应度值。假设如上述三所示，将图像分为m+1类，则m个阈值按顺序排列起来构成一个基因串：

α = {T 1, T 2, \dots, T m}

由于灰度为0~255，所以可以使用8位二进制代码表示每个阈值，此时每个基因串由长度为8*m个比特位的传组成。

将类间方差作为其适应度函数，类间方差越大，适应度函数值就越高。

python代码如下所示：

#coding:utf-8
import cv2
import numpy as np
from matplotlib import pyplot as plt
import random


#将不足8*m位的染色体扩展为8*m位
def expand(k, m):
    for i in range(len(k)):
        k[i] = k[i][:2] + '0'*(8*m+2 - len(k[i])) + k[i][2:len(k[i])]
    return k


def Hist(image):
    a=[0]*256
    h=image.shape[0]
    w=image.shape[1]
    MN=h*w
    average=0.0
    for i in range(w):
        for j in range(h):
            pixel=int(image[j][i])
            a[pixel]=a[pixel]+1
    for i in range(256):
        a[i]=a[i]/float(MN)
        average=average+a[i]*i
    return a, average


#解析多阈值基因串
def getTh(seed, m):
    th = [0, 256]
    seedInt = long(seed, 2)
    for i in range(0, m):
        tmp = seedInt & 255
        if tmp != 0:
            th.append(tmp)
        seedInt = seedInt >> 8
    th.sort()
    return th


#适应度函数 Ostu全局算法
def fitness(seed, p, average, m):
    Var = [0.0] * len(seed)
    g_muT = 0.0

    for i in range(256):
        g_muT = g_muT + i * p[i]

    for i in range(len(seed)):
        th = getTh(seed[i], m)
        for j in range(len(th)-1):
            w = [0.0] * (len(th)-1)
            muT = [0.0] * (len(th)-1)
            mu = [0.0] * (len(th)-1)
            for k in range(th[j], th[j+1]):
                w[j] = w[j] + p[k]
                muT[j] = muT[j] +  + p[k] * k
            if w[j] > 0:
                mu[j] = muT[j] / w[j]
                Var[i] = Var[i] + w[j] * pow(mu[j] - g_muT, 2)
    return Var


#选择算子 轮盘赌选择算法
def wheel_selection(seed, Var):
    var = [0.0]*len(Var)
    s = 0.0
    n = ['']*len(seed)
    sumV = sum(Var)
    for i in range(len(Var)):
        var[i] = Var[i]/sumV
    for i in range(1, len(Var)):
        var[i] = var[i] + var[i-1]
    for i in range(len(seed)):
        s = random.random()
        for j in range(len(var)):
            if s <= var[j]:
                n[i] = seed[j]
    return n


#单点交叉算子
def Cross(Next, m):
    for i in range(0, len(Next) - 1, 2):
        if random.random() < 0.7:
            if m > 2:
                tmp = Next[i][10:]
                Next[i] = Next[i][:10] + Next[i+1][10:]
                Next[i+1] = Next[i+1][:10] + tmp
            else:
                tmp = Next[i][6:]
                Next[i] = Next[i][:6] + Next[i+1][6:]
                Next[i+1] = Next[i+1][:6] + tmp
    return Next


#变异算子
def Variation(Next):
   for i in range(len(Next)):
        if random.random()<0.06:
            Next[i]=bin(long(Next[i],2)+2)
   return Next


#多阈值分割
def genetic_thres(image, k, m):
    th = image
    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            for t in range(1, len(k)-1):
                if k[t-1] <= image[i][j] < k[t]:
                    th[i][j] = int(k[t-1])
    return th


# main
imagesrc = cv2.imread("E:/python/cv/OTSU/test2.jpg")
gray = cv2.cvtColor(imagesrc, cv2.COLOR_BGR2GRAY)

m = 3  #阈值数
items_x = range(0, imagesrc.shape[0])
items_y = range(0, imagesrc.shape[1])
random.shuffle(items_x)
random.shuffle(items_y)
x = items_x[0:20*m]  #产生随机x坐标
y = items_y[0:20*m]  #产生随机y坐标
seed = []
Var = 0.0
times = 0
k = 0
P, average = Hist(gray)  #计算直方图，P为各灰度的概率的数组，average为均值
for i in range(0, 20):
    code = long(0)
    for j in range(0, m):
        code = code + gray[x[i*j]][y[i*j]] << j*8  #将阈值连起来组成一个8*m比特的基因串
    seed.append(bin(code))  #生成第一代

while times < 2000:
    Var = fitness(seed, P, average, m)
    Next = wheel_selection(seed, Var)
    Next = Cross(Next, m)
    Next = expand(Variation(Next), m)
    seed = Next
    times = times + 1

for j in range(len(Var)):
    if Var[j] == max(Var):
        k = getTh(Next[j], m)
print k

plt.subplot(131), plt.imshow(imagesrc, "gray")
plt.title("source image"), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.hist(imagesrc.ravel(), 256)
plt.title("Histogram"), plt.xticks([]), plt.yticks([])
th1 = genetic_thres(gray, k, m)
plt.subplot(133), plt.imshow(th1, "gray")
titleName = ''
for i in range(1, len(k)-1):
    titleName = titleName + str(k[i]) + ', '
titleName = titleName[:len(titleName)-2]
plt.title("threshold is " + titleName), plt.xticks([]), plt.yticks([])
plt.show()

 
 1
2
3
4
5
6
13
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
3
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
1
152
153
154
155
1157
158160
161

这里使用的是标准二进制，若使用格雷码，应该能收敛得更好。结果如下所示：
遗传算法OTSU

xiaopihaierletian

关注

8
点赞
踩
41

收藏

觉得还不错? 一键收藏
1
评论
阈值化分割（二）OTSU法-附Python实现

一、OTSU法（大津阈值分割法）介绍 OTSU算法是由日本学者OTSU于1979年提出的一种对图像进行二值化的高效算法，是一种自适应的阈值确定的方法，又称大津阈值分割法，是最小二乘法意义下的最优分割。二、单阈值OTSU法设图像包含L个灰度级，灰度值为i的像素点个数为Ni，像素总点数为： N=N0+N1+⋯+NL−1则灰度值为i的点的概率为： Pi=Ni
复制链接

扫一扫