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

阈值化分割（二）OTSU法

本人邮箱：[e-mail]，欢迎交流讨论，
欢迎转载，转载请注明网址http://blog.csdn.net/u010128736/

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

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

二、单阈值OTSU法

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

N=N0+N1+⋯+NL−1

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

Pi=NiN

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

μT=∑i=0L−1iPi

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

ω0(T)=∑i=0TPiω1(T)=1−ω0(T)

c0和c1的均值为：

μ0(T)=∑Ti=0iPiω0(T)μ1(T)=μT−∑Ti=0iPiω1(T)

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

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

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

σ2B(T∗)=max0≤T≤L−1{σ2B(T)}

下面给出python源代码。

#coding:utf-8import cv2import numpy as npfrom matplotlib import pyplot as pltimage = 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_OTSUplt.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)}

其中，

σ2B(T1,T2,⋯,Tm)=∑i=0mωi(T1,T2,⋯,Tm)[μi(T1,T2,⋯,Tm)−μT]2μi(T1,T2,⋯,Tm)=∑i=TiTi+1iPiωi(T1,T2,⋯,Tm)ωi(T1,T2,⋯,Tm)=∑i=TiTi+1PiμT=∑i=0L−1iPi

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

四、遗传算法解OTSU

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

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

α={T1,T2,⋯,Tm}

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

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

  python代码如下所示：

#coding:utf-8import cv2import numpy as npfrom matplotlib import pyplot as pltimport 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 kdef 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# mainimagesrc = 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.0times = 0k = 0P, 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 + 1for j in range(len(Var)):    if Var[j] == max(Var):        k = getTh(Next[j], m)print kplt.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()

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