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

• ⭐ Zeal's Blog
• 🛠 Zeal's 知乎专栏
• 🌀 Zeal's Github

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$。每个等级的像素出现的概率如下：

$$\begin{array}{}\text{(1)}& p_i=\frac{n_i}{N},p_i\ge 0,\sum _{i=1}^L{p}_i=1\end{array}$$

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

两个类的像素比例为：

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& 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{array}\end{array}$$

并且

$$\begin{array}{}\text{(3)}& w_0+w_1=1\end{array}$$

两个类的平均灰度值为： $$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\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{array}\end{array}$$

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

$$\begin{array}{}\text{(5)}& \mu (k)=\sum _{i=1}^{k}i\cdot p_i\end{array}$$

$$\begin{array}{}\text{(6)}& {\mu }_T=\sum _{i=1}^{L}i\cdot p_i\end{array}$$

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

$$\begin{array}{}\text{(7)}& w(k)=\sum _{i=1}^k{p}_i\end{array}$$

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

$$\begin{array}{}\text{(8)}& w_0{\mu }_0+w_1{\mu }_1={\mu }_T\end{array}$$

类内/类间方差

两个类各自的方差

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

$$\begin{array}{}\text{(9)}& {\sigma }_0^2=\sum _{i=1}^k(i-{\mu }_0{)}^2\cdot \frac{p_i}{w_0}\end{array}$$

$$\begin{array}{}\text{(10)}& {\sigma }_1^2=\sum _{i=k+1}^L(i-{\mu }_1{)}^2\cdot \frac{p_i}{w_1}\end{array}$$

总体方差

像素总体的方差为：

$$\begin{array}{}\text{(13)}& {\sigma }_T^2=\sum _{i=1}^L(i-{\mu }_T{)}^2\cdot p_i\end{array}$$

类内方差

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

$$\begin{array}{}\text{(11)}& {\sigma }_w^2=w_0{\sigma }_0^2+w_1{\sigma }_1^2\end{array}$$

类间方差

类间（between-class）方差为：

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\sigma }_b^2& =w_0({\mu }_0-{\mu }_T{)}^2+w_1({\mu }_1-{\mu }_T{)}^2\\ & =w_0{w}_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 }_{T}w(k){]}^2}{w(k)[1-w(k)]}\end{array}\end{array}$$

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

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

$$\begin{array}{}\text{(14)}& {\sigma }_T^2={\sigma }_w^2+{\sigma }_b^2\end{array}$$

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

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