使用Matlab/Python实现Otsu算法

摘要

Otsu算法^[1]是在1979年提出的, 它通过最大化类间方差来对图像进行分割。换句话说, Otsu算法是一种非参数分割方法, 根据像素的强度将图像分成不同的区域。本文将分为两个部分来介绍Otsu算法, 首先通过回顾文献中数学公式理解Otsu算法的原理, 最后再使用Matlab/Python实现其细节。

[2]

假设 $L$ 是一个 $M*N$ 大小的图像的像素强度级别:

$$n = n_0 + n_1 + ... + n_{L-1}$$

$$Ph_i = \frac{n_i}{n}, \sum_{i=0}^{L-1}Ph_i=1$$

其中, $n$ 表示图像像素的总数, $n_i$ 表示强度级别 $i$ 的像素数, 所有强度级别的概率分布由 $Ph_i$ 表示。

考虑二分类问题时, 假设存在一个阈值 $th$, 其中 $th$ 介于 $0$ 和 $L-1$ 之间, 那么图像可以根据 $th$ 分为两个类别, 第一个类别 $C_1$ 包含所有像素强度级别在 $0$ ~ $th$ 之间的像素，而 $C_2$ 包含剩余的像素。

$$\omega_1(th) = \sum_{i=0}^{th}Ph_i, \omega_2(th) = \sum_{i=th+1}^{L-1}Ph_i = 1 - Ph_1(th)$$

其中, $\omega_1(th)$ 和 $\omega_2(th)$ 分别表示 $C_1$ 和 $C_2$ 的累积概率分布。

$$\mu_1(th) = \sum_{i=0}^{th}iP(i|C_1) = \sum_{i=0}^{th}\frac{iP(C_1|i)P(i)}{P(C_1)} = \frac{1}{\omega_1(th)}\sum_{i=0}^{th}iPh_i$$

$$\mu_2(th) = \sum_{i=th+1}^{L-1}iP(i|C_2) = \sum_{i=th+1}^{L-1}\frac{iP(C_2|i)P(i)}{P(C_2)} = \frac{1}{\omega_2(th)}\sum_{i=th+1}^{L-1}iPh_i$$

$$\mu_{th} = \sum_{i=0}^{th}iPh_i, \mu_T = \sum_{i=0}^{L-1}iPh_i$$

其中, $\mu_1(th)$ 和 $\mu_2(th)$ 分别表示 $C_1$ 和 $C_2$ 的平均像素强度。$\mu_{th}$ 表示从 $0$ ~ $th$ 的平均像素强度。$\mu_T$ 表示整个图像的平均像素强度。因此, 最大化类别间方差的目标函数如下所示:

$$\omega_1(th) + \omega_2(th) = 1$$

$$\omega_1(th) \cdot \mu_1(th) + \omega_2(th) \cdot \mu_2(th) = \mu_T$$

$$\delta_B^2(th^*) = max\{\delta_B^2(th)\}, 0 \leq th \leq L-1$$

根据上述公式, 我们可以在 $0$ ~ $L-1$ 中找到一个 $th^*$ 来最大化 $σ_B^2$ , 从而实现有效的图像分割。

代码实现

Matlab

function fit = otsu(N, level, x, prob)
    fit = zeros(1, N);
    for j = 1:N
        fit(j) = sum(prob(1:x(j, 1) - 1)) * (sum((1:x(j, 1) - 1) .* ... 
                 prob(1:x(j, 1) - 1) / sum(prob(1:x(j, 1) - 1))) - ...
                 sum((1:255) .* prob(1:255)) ) ^ 2;
        for jlevel = 2:level - 1
            fit(j) = fit(j) + sum(prob(x(j, jlevel - 1):x(j, jlevel) - 1)) * ...
                     (sum((x(j, jlevel - 1):x(j, jlevel) - 1) .* ...
                     prob(x(j, jlevel - 1):x(j, jlevel) - 1) / ...
                     sum(prob(x(j, jlevel - 1):x(j, jlevel) - 1))) - ...
                     sum((1:255) .* prob(1:255))) ^ 2;
        end
        fit(j) = fit(j) + sum(prob(x(j, level - 1):255)) * ...
                 (sum((x(j, level - 1):255) .* prob(x(j, level - 1):255) / ...
                 sum(prob(x(j, level - 1):255))) - sum((1:255) .* prob(1:255))) ^ 2;
        if isnan(fit(j)) || isinf(fit(j))
            fit(j) = eps;
        end
    end
    fit = fit';
end

Python

import numpy as np
from typing import List, Tuple

def otsu(N: int, level: int, x: np.ndarray, prob: np.ndarray) -> np.ndarray:
    fit = np.zeros(N)
    total_mean = np.sum(np.arange(256) * prob[:256])

    for j in range(N):
        fit[j] = 0
        thresholds = [0] + x[j, :].tolist() + [256]
        for i in range(1, len(thresholds)):
            start, end = thresholds[i-1], thresholds[i]
            local_prob = prob[start:end]
            if np.sum(local_prob) == 0:
                continue
            local_mean = np.sum(np.arange(start, end) * local_prob) / np.sum(local_prob)
            fit[j] += np.sum(local_prob) * (local_mean - total_mean) ** 2

        if np.isnan(fit[j]) or np.isinf(fit[j]):
            fit[j] = np.finfo(float).eps

    return fit

• AliPay
• WeChat

参考文献

---

1. Otsu, N. A threshold selection method from gray-level histograms. IEEE transactions on systems, 9, 62-66 (1979). ↩︎

2. Wang, Z., Mo, Y. & Cui, M. An Efficient Multilevel Threshold Image Segmentation Method for COVID-19 Imaging Using Q-Learning Based Golden Jackal Optimization. J Bionic Eng, 20, 2276–2316 (2023). ↩︎