otsu阈值分割

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

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

二、单阈值OTSU法

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

N=N0+N1+⋯+NL−1N=N0+N1+⋯+NL−1

下面给出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()

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

三、多阈值OTSU法

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

σ2B(T∗1,T∗2,⋯,T∗m)=max0≤T1≤T2≤⋯≤L−1{σ2B(T1,T2,⋯,Tm)}σB2(T1∗,T2∗,⋯,Tm∗)=max0≤T1≤T2≤⋯≤L−1{σB2(T1,T2,⋯,Tm)}

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

四、遗传算法解OTSU

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

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

α={T1,T2,⋯,Tm}α={T1,T2,⋯,Tm}

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

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

  python代码如下所示：

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

#将不足8m位的染色体扩展为8m位
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=hw
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

<span class="hljs-keyword">for</span> i <span class="hljs-keyword">in</span> range(<span class="hljs-number">256</span>):
    g_muT = g_muT + i * p[i]

<span class="hljs-keyword">for</span> i <span class="hljs-keyword">in</span> range(len(seed)):
    th = getTh(seed[i], m)
    <span class="hljs-keyword">for</span> j <span class="hljs-keyword">in</span> range(len(th)-<span class="hljs-number">1</span>):
        w = [<span class="hljs-number">0.0</span>] * (len(th)-<span class="hljs-number">1</span>)
        muT = [<span class="hljs-number">0.0</span>] * (len(th)-<span class="hljs-number">1</span>)
        mu = [<span class="hljs-number">0.0</span>] * (len(th)-<span class="hljs-number">1</span>)
        <span class="hljs-keyword">for</span> k <span class="hljs-keyword">in</span> range(th[j], th[j+<span class="hljs-number">1</span>]):
            w[j] = w[j] + p[k]
            muT[j] = muT[j] +  + p[k] * k
        <span class="hljs-keyword">if</span> w[j] &gt; <span class="hljs-number">0</span>:
            mu[j] = muT[j] / w[j]
            Var[i] = Var[i] + w[j] * pow(mu[j] - g_muT, <span class="hljs-number">2</span>)
<span class="hljs-keyword">return</span>
 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:20m] #产生随机x坐标
y = items_y[0:20m] #产生随机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[ij]][y[ij]] << 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()

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