Implementing Otsu's Method Using MATLAB/Python

Abstract

The Otsu's method^[1] was proposed in 1979, which segmented images by maximizing variance between classes. In other words, the Otsu's method is a nonparametric segmentation method that divides an image into different regions based on the intensity of the pixels. This article will introduce the Otsu's method in two sections. First, we will understand the principles of the Otsu's method by reviewing the mathematical formulas in the literature. Finally, we will implement the details using Matlab/Python.

[2]

Assume that $L$ is the pixel intensity levels of an image which has size of $M * N$:

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

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

where $n$ indicates the total number of image pixels, $n_i$ is the number of pixels for intensity level $i$, and probability distribution of all intensity levels is represented by $Ph_i$.

Assume there is a threshold $th$, in which th is between $0$ and $L-1$, then the image can be divided into two classes according to $th$. The first class, $C_1$, contains all pixels with pixel intensity levels between $0$ ~ $th$, while $C_2$ contains the rest of the pixels.

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

where $\omega_1(th)$ and $\omega_2(th)$ represent the cumulative probability distributions for $C_1$ and $C_2$, respectively.

$$\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$$

where $\mu_1(th)$ and $\mu_2(th)$ indicate the mean intensity levels for $C_1$ and $C_2$, respectively. $\mu_{th}$ denotes the mean intensity level from $0$ to $th$. $\mu_T$ represents the mean intensity level for the whole image. Hence, the objective function of the maximizing variance between classes can be expressed as below:

$$\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$$

Based on the above formula, we can find a $th^*$ in the range of $0$ ~ $L-1$ to maximize $σ_B^2$ to achieve a effective image segmentation.

Code

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

References

---

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). ↩︎