Otsu算法（大津法）-动态阈值

发表于 更新于 分类于 ， 阅读次数： Changyan：

Otsu算法简介

Otsu算法是一种图像二值化的算法，它的主要思想是将图像分成背景和前景两部分，背景和前景之间的差别越大越好，也就是说背景越黑，前景越白越好，这样的话，我们就可以通过阈值来将图像二值化，阈值的选取就是Otsu算法的核心。

原始论文1979：A threshold selection method from gray-level histograms

示例图片
示例图片-计算阈值并二值化

虽然Otsu算法是图像中确定阈值的方法，但也可以应用到其他方面。比如点云根据强度值进行分割，就可以使用Otsu算法。本质上Otsu就是通过统计学的方法将输入根据某个值分成两个类。

Otsu算法只能将输入分割成两类，其本质是确定一个阈值，使得用该阈值分割出来的两类的类间（between-class）方差最大，类内（within-class）方差最小。类间方差越大，说明两类之间的差别越大，类内方差越小，说明类内的差别越小。

数学形式

均值和概率

这里主要参考原始论文中的数学形式，以及维基百科中的数学形式。

假设图片中的所有像素灰度值可以分成\(L\)个等级，\([1,2,...,L]\)。比如通常对于8bit灰度图像，我们可以用256个等级表示灰度值。等级\(i\)的像素数为\(n_i\)，总像素数为\(N\)，\(N=n_1+n_2+\cdot\cdot\cdot+n_L\)。每个等级的像素出现的概率如下：

\[ p_i=\frac{n_i}{N}, \quad p_i\geq 0, \sum_{i=1}^{L}p_i=1\tag{1} \]

确定阈值\(k\)后可以将像素分成两类\(C_0,C_1\) \(C_0\)表示像素等级在\([1,...,k]\)的像素，\(C_1\)表示像素等级在\([k+1,...,L]\)的像素

两个类的像素比例为：

\[ \begin{aligned} &w_0=Pr(C_0)=\sum_{i=1}^{k}p_i=w(k)\\ &w_1=Pr(C_1)=\sum_{i=k+1}^{L}p_i=1-w(k)\\ \end{aligned}\tag{2} \]

并且

\[ w_0+w_1=1\tag{3} \]

两个类的平均灰度值为： \[ \begin{aligned} &\mu_0=E(C_0)=\sum_{i=1}^{k} i\cdot \frac{p_i}{w_0}=\mu(k)/w(k)\\ &\mu_1=E(C_1)=\sum_{i=k+1}^{L} i\cdot \frac{p_i}{w_1}=\frac{\mu_T - \mu(k)}{1 - w(k)}\\ \end{aligned}\tag{4} \]

其中，\(\mu(k)\)是等级为\(1-k\)像素的一阶累积矩（first-order cumulative moments），\(\mu_T\)是总的平均灰度值。

\[ \mu(k)=\sum_{i=1}^{k}i\cdot p_i\tag{5} \]

\[ \mu_T=\sum_{i=1}^{L}i\cdot p_i\tag{6} \]

\(w(k)\)是等级为\(1-k\)像素的概率和

\[ w(k)=\sum_{i=1}^{k}p_i\tag{7} \]

根据上面关系，我们还可以得到：

\[ w_0\mu_0+w_1\mu_1 = \mu_T\tag{8} \]

类内/类间方差

两个类各自的方差

根据上面的定义，我们可以得到两个类各自的方差

\[ \sigma_0^2=\sum_{i=1}^{k}(i-\mu_0)^2\cdot \frac{p_i}{w_0}\tag{9} \]

\[ \sigma_1^2=\sum_{i=k+1}^{L}(i-\mu_1)^2\cdot \frac{p_i}{w_1}\tag{10} \]

总体方差

像素总体的方差为：

\[ \sigma_T^2=\sum_{i=1}^{L}(i-\mu_T)^2\cdot p_i\tag{13} \]

类内方差

因此，我们可以得到类内（within-class）方差的公式：

\[ \sigma_w^2=w_0\sigma_0^2+w_1\sigma_1^2\tag{11} \]

类间方差

类间（between-class）方差为：

\[ \begin{aligned} \sigma_b^2&=w_0(\mu_0-\mu_T)^2+w_1(\mu_1-\mu_T)^2\\ &=w_0w_1(\mu_0-\mu_1)^2\\ &=w(k)(1-w(k))[\frac{\mu(k)}{w(k)}-\frac{\mu_T-\mu(k)}{1-w(k)}]^2\\ &=\frac{[\mu(k)(1-w(k))-(\mu_T-\mu(k))w(k)]^2}{w(k)(1-w(k))}\\ &=\frac{[\mu(k)-\mu_Tw(k)]^2}{w(k)[1-w(k)]} \end{aligned}\tag{12} \]

类内方差与类间方差的关系

根据公式（11-13），我们可以得到：

\[ \sigma_T^2=\sigma_w^2+\sigma_b^2\tag{14} \]

Otsu算法的目标

Otsu算法的目标是找到一个阈值\(k\)，使得类间方差\(\sigma_b^2\)最大，也就是说，找到一个阈值\(k\)，使得类内方差\(\sigma_w^2\)最小。

Otsu算法步骤

  1. 统计每个像素等级的像素数量，并计算每个像素等级的概率\(p_i\)，\(i\in[1,...,L]\)。
  2. 计算每个像素等级的概率和\(w(k)\)，\(k\in[1,...,L]\)。
  3. 计算每个像素等级的平均灰度值\(\mu(k)\)，\(k\in[1,...,L]\)。
  4. 计算总的平均灰度值\(\mu_T\)。
  5. 根据公式（4）可以计算出每个像素等级的平均灰度值\(\mu_0\)和\(\mu_1\)，\(k\in[1,...,L]\)。
  6. 遍历所有像素等级\([1,2,3...,L]\) 6.1 根据公式（12）计算类间方差\(\sigma_b^2\)。 6.2 如果\(\sigma_b^2\)大于之前的最大值，则更新最大值，并记录当前的像素等级\(k\)

算法实现-python

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import numpy as np
import matplotlib.pyplot as plt
from skimage import io, img_as_float
from skimage.filters import threshold_otsu


def otsu(gray_img) -> int:
    """
    Otsu算法
    :param gray_img: 灰度图像
    :return: 阈值
    """
    # 统计每个像素等级的像素数量，并计算每个像素等级的概率
    # 灰度图像的像素等级为[0,1,2,...,255]
    pixel_nums = np.zeros((256, 1))  # 每个像素等级的像素数量
    pixel_probs = np.zeros((256, 1))  # 每个像素等级的概率
    for i in range(gray_img.shape[0]):
        for j in range(gray_img.shape[1]):
            pixel_nums[gray_img[i][j]] += 1
    pixel_probs = pixel_nums / (gray_img.shape[0] * gray_img.shape[1])

    # 计算每个像素等级的概率和
    pixel_probs_sum = np.zeros((256, 1))  # 每个像素等级的概率和
    for i in range(256):
        if i == 0:
            pixel_probs_sum[i] = pixel_probs[i]
        else:
            pixel_probs_sum[i] = pixel_probs_sum[i - 1] + pixel_probs[i]

    # 计算每个像素等级的平均灰度值
    pixel_means = np.zeros((256, 1))  # 每个像素等级的平均灰度值
    for i in range(256):
        if i == 0:
            pixel_means[i] = 0 * pixel_probs[i]
        else:
            pixel_means[i] = pixel_means[i - 1] + i * pixel_probs[i]

    # 计算总的平均灰度值
    img_mean = pixel_means[-1]

    # 遍历所有像素等级
    max_sigma_b = 0  # 最大类间方差
    best_threshold = 0
    for i in range(256):
        # 根据上文公式（12）计算类间方差
        # \frac{[\mu(k)-\mu_Tw(k)]^2}{w(k)[1-w(k)]}
        sigma_b = (img_mean * pixel_probs_sum[i] - pixel_means[i]) ** 2 / (pixel_probs_sum[i] * (1 - pixel_probs_sum[i]))
        # 如果类间方差大于之前的最大值，则更新最大值，并记录当前的像素等级
        if sigma_b > max_sigma_b:
            max_sigma_b = sigma
            k = i