基于遗传算法的图像多阈值分割

基于遗传算法对图像进行多阈值分割

一、遗传算法基本原理

0、遗传学算法

遗传算法由Holland提出来的最初的目的是研究自然系统的自适应行为并设计具有自适应功能的软件系统。它的特点是对参数进行编码运算不需要有关体系的任何先验知识，沿多种路线进行平行搜索不会落人局部较优的陷阱。

进化算法最初是借鉴了进化生物学中的一些现象而发展起来的，这些现象包括遗传、突变.自然选择以及杂交等。遗传算法在适应度函数选择不当的情况下有可能收敛于局部最优,而不能达到全局最优。

遗传算法是从代表问题可能潜在的解集的一个种群(population)开始的,而一个种群则由经过基因(gene)编码的一定数目的个体(individual)组成。

每个个体实际上是染色体(chromosome)带有特征的实体。染色体作为遗传物质的主要载体，即多个基因的集合,其内部表现(即基因型)是某种基因组合,它决定了个体的形状的外部表现，如黑头发的特征是由染色体中控制这一特征的某种基因组合决定的。因此,在一开始需要实现从表现型到基因型的映射即编码工作。

由于仿照基因编码的工作很复杂,往往进行简化，如二进制编码,初代种群产生之后,按照适者生存和优胜劣汰的原理,逐代(generation)演化产生出越来越好的近似解。

在每一代,根据问题域中个体的适应度(fitness)大小选择(selection)个体，并借助于自然遗传学的遗传算子(genetic operators)进行组合交叉(crossover)和变异(mutation)，产生出代表新的解集的种群。

这个过程将导致种群像自然进化一样的后生代种群比前代更加适应于环境,末代种群中的最优个体经过解码(decoding),可以作为问题近似最优解。

来源：(42条消息) 遗传算法原理介绍_逝年！但知行好事，莫要问前程。的博客-CSDN博客_遗传算法原理

1、生物学启发

相信你还记得这句话：「细胞是所有生物的基石。」由此可知，在一个生物的任何一个细胞中，都有着相同的一套染色体。所谓染色体，就是指由 DNA 组成的聚合体。

传统上看，这些染色体可以被由数字 0 和 1 组成的字符串表达出来。

一条染色体由基因组成，这些基因其实就是组成 DNA 的基本结构，DNA 上的每个基因都编码了一个独特的性状

2、遗传算法的一般步骤

1. 评估每条染色体所对应个体的适应度。

2. 遵照适应度越高，选择概率越大的原则，从种群中选择两个个体作为父方和母方。

3. 抽取父母双方的染色体，进行交叉，产生子代。

4. 对子代的染色体进行变异。

5. 重复2，3，4步骤，直到新种群的产生。

二、图像分割

0、图像阈值

数字图像可以视为二维矩阵或者二元函数。它由称为像素的离散点组成。在彩色图像中，每个像素有三个值：红、绿和蓝。每个值在0到L-1的范围中，其中L是精度的数量。另一方面，灰度图像由像素组成，其中每个像素只有0到L-1之间的一个值，称之为灰度。对于许多图像处理问题，处理灰度图像比处理彩色图像更简单、有效。因此，在使用图像处理算法前，会将彩色图像转换为灰度图像。最常用的灰度为256（即每个像素值在0到255之间）。

图像阈值是一种用于灰度图像的图像分割方法。该想法是为了找到一个阈值，如果像素低于阈值，就认为是背景，否则认为是目标的一部分。

1、Otsu算法（最大类间方差）

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

关于Otsu算法的详解参考：otsu（大津法）_百度百科 (baidu.com)

设图像包含L个灰度级，灰度值为i的像素点个数为Ni，像素总点数为N
$$P_i=\frac{N_i}{N}$$
根据期望公式，灰度图像的均值：
$${\mu }_T=\sum _{i=0}^{L-1}i{P}_i$$
按图像的灰度特性，使用阈值T将图像分成目标c0和背景c1两类，则ω0(T)和ω1(T)分别表示阈值为T时，c0和c1发生的概率，即：
$$\begin{gathered} {\omega }_0(T)=\sum _{i=0}^T{P}_i \\ {\omega }_1(T)=1-{\omega }_0(T) \end{gathered}$$
c0和c1的均值为：
$$\begin{gathered} {\mu }_0(T)=\frac{{\sum }_{i=0}^{T}i{P}_i}{{\omega }_0(T)} \\ {\mu }_1(T)=\frac{{\sum }_{i=T+1}^{L-1}i{P}_i}{{\omega }_1(T)} \end{gathered}$$
${\sigma }^2(T)$表示直方图中阈值为T的类间方差，定义为：
$${\sigma }^2(T)={\omega }_0(T)[{\mu }_0(T)-{\mu }_T{]}^2+{\omega }_1(T)[{\mu }_1(T)-{\mu }_T{]}^2$$
最优阈值定义为类间方差最大时对应的T值，即：
σ 2 ( T ∗ ) = max ⁡ 0 ≤ T ≤ L − 1 { σ 2 ( T ) } σ^2(T^*) = \max_{0\leq T \leq L-1}\{σ^2(T)\} σ2(T∗)=0≤T≤L−1max​{σ2(T)}
将单阈值推广到多阈值：$threshold=[T_1,T_2,\cdot \cdot \cdot ,T_m]$

则最大类间方差为：
$${\sigma }_B^2(T_1^{\ast },T_2^{\ast },\cdots ,T_m^{\ast })=\underset{0\le T1\le T2\le \cdots \le L-1}{\max }{\sigma }_B^2(T_1,T_2,\cdots ,T_m)$$

2、遗传算法解OTSU

1. 初始化种群

   # 各参数含义
   #	 threshold_num = 10 # 阈值个数
   #    threshold = [] # 阈值
   #    population_size = 20 # 种群大小
   #    iter_num = 2000 # 迭代次数
   # 初始化种群
   def init_population(img,threshold_num,population_size):
       seeds = []
       for _ in range(population_size):
           chromosome = create_chromosome(threshold_num)
           seeds.append(chromosome)  #生成第一代
       return seeds
   

2. 染色体编码

   # 染色体编码（随机产生）
   # 阈值为[127,10]对应的染色体为 '0111_1111_0000_1010' 以字符串形式储存(没有下划线)
   def create_chromosome(threshold_num):
       chromosome = ''
       for _ in range(threshold_num):
           #将多个阈值连起来组成一个8*m比特的基因串
           chromosome = chromosome + "{:0>8b}".format(int(random.random()*255))
       return chromosome
   

3. 染色体解码

   #  解析多阈值基因串
   #  chromosome为m*8长度的二进制字符串
   #  '01111111_00001010' -> [127,10]
   #  返回值为十进制阈值列表
   def getThreshold(chromosome, m):
       threshold = [0, 256]
       for i in range(m):
           child_seed = chromosome[i*8:i*8+8]
           seedInt = int(child_seed,2)
           tmp = seedInt & 255
           if tmp != 0:
               threshold.append(tmp)
       threshold.sort()#将阈值从小到大排序
       return threshold
   

4. 评估每条染色体所对应个体的适应度。

   def fitness(population, p, average, threshold_num):
       Var = [0.0] * len(population)
       g_muT = 0.0
       for i in range(256):
           g_muT = g_muT + i * p[i]
       for i in range(len(population)):
           th = getThreshold(population[i], threshold_num)
           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
   

5. 遵照适应度越高，选择概率越大的原则，从筛选新的种群

   #轮盘赌法选择下一代种群
   def wheel_selection(population, Var,threshold_num):
       var = [0.0]*len(Var)
       s = 0.0
       next_seeds = ['']*len(population)
       sumV = sum(Var)
       max_var = max(Var)
   
       for i in range(len(Var)):# 最优良的性状永远被保留
           if Var[i] == max_var:
               next_seeds[0] = population[i]
       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(1,len(population)):
           s = random.random()
           for j in range(len(var)):
               if s <= var[j]:
                   next_seeds[i] = population[j]
       for i in range(len(population)):
           if not next_seeds[i] :
               next_seeds[i] = create_chromosome(threshold_num)
       return next_seeds
   

6. 抽取父母双方的染色体，进行交叉，产生子代。

   #染色体单点交叉
   def Crossover(population, threshold_num):
       for i in range(0, len(population) - 1, 2):
           if random.random() < 0.7:
               if threshold_num > 2:
                   tmp = population[i][10:]
                   population[i] = population[i][:10] + population[i+1][10:]
                   population[i+1] = population[i+1][:10] + tmp
               else:
                   tmp = population[i][6:]
                   population[i] = population[i][:6] + population[i+1][6:]
                   population[i+1] = population[i+1][:6] + tmp
       return population
   

7. 对子代的染色体进行变异。

   #基因变异
   def Mutations(population):
      for i in range(len(population)):
           if random.random() < 0.06:
               j = int(random.random()*len(population[i]))
               if len(population[i]) > 8:
                   population[i] = "{}{:>08b}".format(population[i][:-8],(int(population[i][-8:])+2)//256)
               else:
                   population[i] = "{:>08b}".format((int(population[i],2)+2)%256)
      return population
   

8. main函数

       filename = 'npu.png'
       image = cv2.imread(filename)
       gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
       threshold_num = 10 # 阈值个数
       threshold = 0 # 阈值
       population_size = 20 # 种群大小
       iter_num = 2000 # 迭代次数
       P, average = Hist(gray)  #计算直方图，P为各灰度的概率的数组，average为均值
       
       population = init_population(gray, threshold_num,population_size)
       Var = 0.0
   
       # 迭代
       for _ in range(iter_num):
           Var = fitness(population, P, average, threshold_num)
           next_population = wheel_selection(population, Var,threshold_num)
           next_population = Crossover(next_population, threshold_num)
           next_population = Mutations(next_population)
           population = next_population
   
       # 求出最大var对应的阈值
       for j in range(len(Var)):
           if Var[j] == max(Var):
               threshold = getThreshold(next_population[j], threshold_num)
   
       print(threshold)
       '''绘图'''
       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([])
       new_image = genetic_thres(gray, threshold, threshold_num)
       cv2.imwrite(f'{filename.split(".")[0]}_th1.jpg',new_image)
       plt.subplot(133), plt.imshow(new_image, "gray")
       titleName = ''
       for i in range(1, len(threshold)-1):
           titleName = titleName + str(threshold[i]) + ', '
       titleName = titleName[:len(titleName)-2]
       (plt.title(f"threshold is {titleName}"), plt.xticks([]), plt.yticks([]))
       plt.savefig(f'{filename.split(".")[0]}_result.jpg')
       plt.show()
   

9. 运行结果

   原图：

   运行图：

   OTSU单阈值分割图像：

遗传算法多阈值分割图像：

三、完整代码

import cv2
import itertools
import numpy as np
from matplotlib import pyplot as plt
import random


# 获取每个灰度值像素的个数和平均灰度值
def Hist(image):
    pixels = cv2.calcHist([image], [0], None, [256], [0, 256])
    average = np.mean(image)
    pixels = np.array(pixels)/image.size
    return pixels.T[0],average

#  解析多阈值基因串
# chromosome为m*8长度的二进制字符串
# 返回值为十进制阈值列表
def getThreshold(chromosome, m):
    threshold = [0, 256]
    for i in range(m):
        child_seed = chromosome[i*8:i*8+8]
        seedInt = int(child_seed,2)
        tmp = seedInt & 255
        if tmp != 0:
            threshold.append(tmp)
    threshold.sort()#将阈值从小到大排序
    return threshold

#适应度函数 Otsu全局算法
def fitness(population, p, average, threshold_num):
    Var = [0.0] * len(population)
    g_muT = 0.0
    for i in range(256):
        g_muT = g_muT + i * p[i]
    for i in range(len(population)):
        th = getThreshold(population[i], threshold_num)
        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(population, Var,threshold_num):
    var = [0.0]*len(Var)
    s = 0.0
    next_seeds = ['']*len(population)
    sumV = sum(Var)
    max_var = max(Var)

    for i in range(len(Var)):# 最优良的性状永远被保留
        if Var[i] == max_var:
            next_seeds[0] = population[i]

    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(1,len(population)):
        s = random.random()
        for j in range(len(var)):
            if s <= var[j]:
                next_seeds[i] = population[j]

    for i in range(len(population)):
        if not next_seeds[i] :
            next_seeds[i] = create_chromosome(threshold_num)
    return next_seeds

#染色体单点交叉
def Crossover(population, threshold_num):
    for i in range(0, len(population) - 1, 2):
        if random.random() < 0.7:
            if threshold_num > 2:
                tmp = population[i][10:]
                population[i] = population[i][:10] + population[i+1][10:]
                population[i+1] = population[i+1][:10] + tmp
            else:
                tmp = population[i][6:]
                population[i] = population[i][:6] + population[i+1][6:]
                population[i+1] = population[i+1][:6] + tmp
    return population

#基因变异
def Mutations(population):
   for i in range(len(population)):
        if random.random() < 0.06:
            j = int(random.random()*len(population[i]))
            if len(population[i]) > 8:
                population[i] = "{}{:>08b}".format(population[i][:-8],(int(population[i][-8:])+2)//256)
            else:
                population[i] = "{:>08b}".format((int(population[i],2)+2)%256)
   return population


# 初始化种群
def init_population(img,threshold_num,population_size):
    seeds = []
    for _ in range(population_size):
        chromosome = create_chromosome(threshold_num)
        seeds.append(chromosome)  #生成第一代
    return seeds
# 随机产生一个染色体
def create_chromosome(threshold_num):
    chromosome = ''
    for _ in range(threshold_num):
        chromosome = chromosome + "{:0>8b}".format(int(random.random()*255))#将阈值连起来组成一个8*m比特的基因串
    return chromosome
# 根据阈值分割图像
def genetic_thres(image, k, m):
    new_image = image
    for i, j in itertools.product(range(image.shape[0]), range(image.shape[1])):
        for t in range(1, len(k)):
            if k[t-1] <= image[i][j] < k[t]:
                new_image[i][j] = int(k[t-1])
    return new_image
def main():
    filename = 'screenshot.png'
    image = cv2.imread(filename)
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    threshold_num = 10 # 阈值个数
    threshold = 0 # 阈值
    population_size = 20 # 种群大小
    iter_num = 2000 # 迭代次数
    P, average = Hist(gray)  #计算直方图，P为各灰度的概率的数组，average为均值
    
    population = init_population(gray, threshold_num,population_size)
    Var = 0.0

    # 迭代
    for _ in range(iter_num):
        Var = fitness(population, P, average, threshold_num)
        next_population = wheel_selection(population, Var,threshold_num)
        next_population = Crossover(next_population, threshold_num)
        next_population = Mutations(next_population)
        population = next_population

    # 求出最大var对应的阈值
    for j in range(len(Var)):
        if Var[j] == max(Var):
            threshold = getThreshold(next_population[j], threshold_num)

    print(threshold)
    '''绘图'''
    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([])
    new_image = genetic_thres(gray, threshold, threshold_num)
    cv2.imwrite(f'{filename.split(".")[0]}_th1.jpg',new_image)
    plt.subplot(133), plt.imshow(new_image, "gray")
    titleName = ''
    for i in range(1, len(threshold)-1):
        titleName = titleName + str(threshold[i]) + ', '
    titleName = titleName[:len(titleName)-2]
    (plt.title(f"threshold is {titleName}"), plt.xticks([]), plt.yticks([]))
    plt.savefig(f'{filename.split(".")[0]}_result.jpg')
    plt.show()
main()