OTSU算法在图像处理中的应用

博主近期需要处理原子力显微镜图，测量AFM相图中片晶厚度，本着一劳永逸的想法，直接代码实现这个目标。这里博主放上实现这个目标基础但很关键的一步代码，即通过大津算法（OTSU算法）求出图片前景和背景的阈值，并基于此对图像进行分割处理，已实现后面对图片中信息的处理和挖掘。

这里博主编写的是单阈值OTSU法，采用的是面向对象的写法。首先简单介绍一下OTUS算法的原理，一张图片粗略来分可以分为前景和背景，当两部分类间方差最大时，对应的阈值（threshold)即为二分的最佳阈值。具体数学实现过程如下：

首先将图像灰度化（0-255），假设图中包含m个灰度级，各个灰度级像素点数分别为$N_0,N_1,…,N_{m-1}$,那么灰度值为i(i$\epsilon$m)的占比$P_i$为:

$P_i=\frac{N_i}{N}$

图像灰度的均值为:

$Avg = \sum_{i=0}^{m-1}{i*Pi}$

这里先定义阈值为thr, 则thr将图像分为前景(foreground, fg)和背景(blackground, bg), 两者占比分别为$P_{fg}$和$P_{bg}$, 前景和背景的灰度均值为$Avg_{fg}$和$Avg_{bg}$则有:

$P_{fg}=\sum_{i=0}^{thr}Pi$ $P_{bg}=1-\sum_{i=0}^{thr}Pi$ $Avg_{fg}=\frac{\sum_{i=0}^{thr}{i*Pi}}{P_{fg}}$ $Avg_{bg}=\frac{Avg-\sum_{i=0}^{thr}{i*Pi}}{P_{bg}}$

定义前景和背景的类间方差${\sigma}^2$,则有:

${\sigma}^2=P_{fg}*(Avg-Avg_{fg})^2+P_{bg}*(Avg-Avg_{bg})^2$

${\sigma}^2$最大时对应的thr即为最优阈值$thr_{optimal}$。

from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import font_manager as fm, rcParams


class img_process(object):
    """用来处理图像以及分析图中的数据"""

    def __init__(self, img):
        img.show()
        self.img_array = np.array(img)
        self.img_array_fl = self.img_array.flatten()
        # 图像名义尺寸和实际尺寸
        self.pix_weight, self.pix_height = img.size
        self.N = self.pix_weight * self.pix_height
        self.Avg = 0
        self.P_list = []
        self.Avg_list = []
        # 将图像数组转化为灰度直方图
        self.gray_list = img.histogram()
        self.kdet = 0
        self.kt = 0

    def gray_his(self):
        """绘制图像的灰度直方图"""
        plt.rcParams['font.sans-serif'] = ['SimHei']
        plt.rcParams['axes.unicode_minus'] = False
        n, bins, patches = plt.hist(self.img_array_fl, bins=256, align='mid')
        plt.xlim(0, 256)
        plt.ylim((0, max(self.gray_list)))
        plt.xlabel('灰度值')
        plt.ylabel('像素点频数')
        plt.title('灰度分布直方图')
        # plt.savefig('C:/Users/yplit/Pictures/灰度直方图.jpg')
        plt.show()

    def otsu_process(self):

        # 大津算法求图像分割阈值

        """ 列出所有要用到的参数,总的像素点个数N,灰度值为i的像素点个数Ni,灰度值为i占比为Pi,，
         像素点总均值Avg,图像分割阈值thr,前景包含的像素点的类概率及类均值为P_fg和Avg_fg,背景包含
         的像素点的类概率及类均值为P_bg和Avg_bg,类间方差为det
        """

        for i in range(len(self.gray_list)):
            Pi = self.gray_list[i] / self.N
            self.P_list.append(Pi)
            Avg_i = i * Pi
            self.Avg_list.append(Avg_i)
            self.Avg += Avg_i

        for thr in range(len(self.gray_list)):
            P_fg = sum(self.P_list[:thr + 1])
            P_bg = 1 - P_fg
            Avg_fg = sum(self.Avg_list[:thr + 1]) / P_fg
            if P_bg != 0:
                Avg_bg = (self.Avg - Avg_fg *  P_fg) / P_bg
                det = P_fg * (Avg_fg - self.Avg) ** 2 + P_bg * (self.Avg - Avg_bg) ** 2
                if det >= self.kdet:
                    self.kdet, self.kt = det, thr
        print('图像分割阈值为%a' % self.kt)
        return self.kt

    def b_w(self, t, j = True):
        """已计算得到的阈值为依据将图像转化为二进制图,j默认为True，
        根据图片的实际情况选择"""
        if j:
            for i in range(self.pix_height):
                for j in range(self.pix_weight):
                    if self.img_array[i][j] <= t:
                        self.img_array[i][j] = 0
                    else:
                        self.img_array[i][j] = 255
        else:
            for i in range(self.pix_height):
                for j in range(self.pix_weight):
                    if self.img_array[i][j] <= t:
                        self.img_array[i][j] = 255
                    else:
                        self.img_array[i][j] = 0
        img = Image.fromarray(self.img_array)
        img.show()
        return self.img_array


if __name__ == '__main__':
    # 打开图像
    img = Image.open(
        'C:/Users/yplit/Desktop/微信图片_20200823190646.jpg')
    # RGB转灰度图 两种方法(直接调用covert转化或者利用numpy写函数)
    img = img.convert('L')
    img_obj = img_process(img)
    t = img_obj.otsu_process()
    img_distance = img_obj.b_w(t, True)
    img_obj.gray_his()

    # img_array = np.array(img)
    # gray = np.dot(img_array[..., :3], (0.299, 0.587, 0.114))
    # img = Image.fromarray(gray)
    # img.show()

事实上python有自带的otsu算法函数，这里博主只是根据算法的数学原理去实现了对最优阈值的求解，那么通过这个最优阈值对图像进行二分后就可以继续图像信息的挖掘。