# Implementing Otsu's Method Using MATLAB/Python

### # Abstract

The Otsu's methodopen in new window[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 Matlabopen in new window/Pythonopen in new window.

### # Formula[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

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