Otsu's method

1. Compute histogram and probabilities of each intensity level
2. Set up initial ${\displaystyle \omega _{i}(0)}$ and ${\displaystyle \mu _{i}(0)}$
3. Step through all possible thresholds ${\displaystyle t=1,\ldots }$ maximum intensity
   1. Update ${\displaystyle \omega _{i}}$ and ${\displaystyle \mu _{i}}$
   2. Compute ${\displaystyle \sigma _{b}^{2}(t)}$
4. Desired threshold corresponds to the maximum ${\displaystyle \sigma _{b}^{2}(t)}$

histogramCounts is a 256-element histogram of a grayscale image different gray-levels (typical for 8-bit images). level is the threshold for the image (double).

function level = otsu(histogramCounts)
total = sum(histogramCounts); % total number of pixels in the image 
%% OTSU automatic thresholding
top = 256;
sumB = 0;
wB = 0;
maximum = 0.0;
sum1 = dot(0:top-1, histogramCounts);
for ii = 1:top
    wF = total - wB;
    if wB > 0 && wF > 0
        mF = (sum1 - sumB) / wF;
        val = wB * wF * ((sumB / wB) - mF) * ((sumB / wB) - mF);
        if ( val >= maximum )
            level = ii;
            maximum = val;
        end
    end
    wB = wB + histogramCounts(ii);
    sumB = sumB + (ii-1) * histogramCounts(ii);
end
end

Matlab has built-in functions graythresh() and multithresh() in the Image Processing Toolbox which are implemented with Otsu's method and Multi Otsu's method, respectively.

This implementation requires the NumPy library.

import numpy as np

def compute_otsu_criteria(im, th):
    """Otsu's method to compute criteria."""
    # create the thresholded image
    thresholded_im = np.zeros(im.shape)
    thresholded_im[im >= th] = 1

    # compute weights
    nb_pixels = im.size
    nb_pixels1 = np.count_nonzero(thresholded_im)
    weight1 = nb_pixels1 / nb_pixels
    weight0 = 1 - weight1

    # if one of the classes is empty, eg all pixels are below or above the threshold, that threshold will not be considered
    # in the search for the best threshold
    if weight1 == 0 or weight0 == 0:
        return np.inf

    # find all pixels belonging to each class
    val_pixels1 = im[thresholded_im == 1]
    val_pixels0 = im[thresholded_im == 0]

    # compute variance of these classes
    var1 = np.var(val_pixels1) if len(val_pixels1) > 0 else 0
    var0 = np.var(val_pixels0) if len(val_pixels0) > 0 else 0

    return weight0 * var0 + weight1 * var1

im = # load your image as a numpy array.
# For testing purposes, one can use for example im = np.random.randint(0,255, size = (50,50))

# testing all thresholds from 0 to the maximum of the image
threshold_range = range(np.max(im)+1)
criterias = [compute_otsu_criteria(im, th) for th in threshold_range]

# best threshold is the one minimizing the Otsu criteria
best_threshold = threshold_range[np.argmin(criterias)]

Python libraries dedicated to image processing such as OpenCV and Scikit-image propose built-in implementations of the algorithm.

A popular local adaptation is the two-dimensional Otsu's method, which performs better for the object segmentation task in noisy images. Here, the intensity value of a given pixel is compared with the average intensity of its immediate neighborhood to improve segmentation results.[8]

At each pixel, the average gray-level value of the neighborhood is calculated. Let the gray level of the given pixel be divided into ${\displaystyle L}$ discrete values and the average gray level is also divided into the same ${\displaystyle L}$ values. Then a pair is formed: the pixel gray level and the average of the neighborhood ${\displaystyle (i,j)}$. Each pair belongs to one of the ${\displaystyle L\times L}$ possible 2-dimensional bins. The total number of occurrences (frequency), ${\displaystyle f_{ij}}$, of a pair ${\displaystyle (i,j)}$, divided by the total number of pixels in the image ${\displaystyle N}$, defines the joint probability mass function in a 2-dimensional histogram:

${\displaystyle P_{ij}={\frac {f_{ij}}{N}},\qquad \sum _{i=0}^{L-1}\sum _{j=0}^{L-1}P_{ij}=1}$

And the 2-dimensional Otsu's method is developed based on the 2-dimensional histogram as follows.

The probabilities of two classes can be denoted as:

${\displaystyle {\begin{aligned}\omega _{0}&=\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}P_{ij}\\\omega _{1}&=\sum _{i=s}^{L-1}\sum _{j=t}^{L-1}P_{ij}\end{aligned}}}$

The intensity mean value vectors of two classes and total mean vector can be expressed as follows:

${\displaystyle {\begin{aligned}\mu _{0}&=[\mu _{0i},\mu _{0j}]^{T}=\left[\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}i{\frac {P_{ij}}{\omega _{0}}},\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}j{\frac {P_{ij}}{\omega _{0}}}\right]^{T}\\\mu _{1}&=[\mu _{1i},\mu _{1j}]^{T}=\left[\sum _{i=s}^{L-1}\sum _{j=t}^{L-1}i{\frac {P_{ij}}{\omega _{1}}},\sum _{i=s}^{L-1}\sum _{j=t}^{L-1}j{\frac {P_{ij}}{\omega _{1}}}\right]^{T}\\\mu _{T}&=[\mu _{Ti},\mu _{Tj}]^{T}=\left[\sum _{i=0}^{L-1}\sum _{j=0}^{L-1}iP_{ij},\sum _{i=0}^{L-1}\sum _{j=0}^{L-1}jP_{ij}\right]^{T}\end{aligned}}}$

In most cases the probability off-diagonal will be negligible, so it is easy to verify:

${\displaystyle \omega _{0}+\omega _{1}\cong 1}$

${\displaystyle \omega _{0}\mu _{0}+\omega _{1}\mu _{1}\cong \mu _{T}}$

The inter-class discrete matrix is defined as

${\displaystyle S_{b}=\sum _{k=0}^{1}\omega _{k}[(\mu _{k}-\mu _{T})(\mu _{k}-\mu _{T})^{T}]}$

The trace of the discrete matrix can be expressed as

${\displaystyle {\begin{aligned}&\operatorname {tr} (S_{b})\\[4pt]={}&\omega _{0}[(\mu _{0i}-\mu _{Ti})^{2}+(\mu _{0j}-\mu _{Tj})^{2}]+\omega _{1}[(\mu _{1i}-\mu _{Ti})^{2}+(\mu _{1j}-\mu _{Tj})^{2}]\\[4pt]={}&{\frac {(\mu _{Ti}\omega _{0}-\mu _{i})^{2}+(\mu _{Tj}\omega _{0}-\mu _{j})^{2}}{\omega _{0}(1-\omega _{0})}}\end{aligned}}}$

where

${\displaystyle \mu _{i}=\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}iP_{ij}}$

${\displaystyle \mu _{j}=\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}jP_{ij}}$

Similar to one-dimensional Otsu's method, the optimal threshold ${\displaystyle (s,t)}$ is obtained by maximizing ${\displaystyle \operatorname {tr} (S_{b})}$.

The ${\displaystyle s}$ and ${\displaystyle t}$ is obtained iteratively which is similar with one-dimensional Otsu's method. The values of ${\displaystyle s}$ and ${\displaystyle t}$ are changed till we obtain the maximum of ${\displaystyle \operatorname {tr} (S_{b})}$, that is

max,s,t = 0;

for ss: 0 to L-1 do
    for tt: 0 to L-1 do
        evaluate tr(S_b);
        if tr(S_b) > max
            max = tr(S,b);
            s = ss;
            t = tt;
        end if
    end for
end for

return s,t;

Notice that for evaluating ${\displaystyle \operatorname {tr} (S_{b})}$, we can use a fast recursive dynamic programming algorithm to improve time performance.[12] However, even with the dynamic programming approach, 2d Otsu's method still has large time complexity. Therefore, much research has been done to reduce the computation cost.[13]

If summed area tables are used to build the 3 tables, sum over ${\displaystyle P_{ij}}$, sum over ${\displaystyle i*P_{ij}}$, and sum over ${\displaystyle j*P_{ij}}$, then the runtime complexity is the maximum of (O(N_pixels), O(N_bins*N_bins)). Note that if only coarse resolution is needed in terms of threshold, N_bins can be reduced.

function inputs and output:

hists is a ${\displaystyle 256\times 256}$ 2D-histogram of grayscale value and neighborhood average grayscale value pair.

total is the number of pairs in the given image.it is determined by the number of the bins of 2D-histogram at each direction.

threshold is the threshold obtained.

function threshold = otsu_2D(hists, total)
maximum = 0.0;
threshold = 0;
helperVec = 0:255;
mu_t0 = sum(sum(repmat(helperVec',1,256).*hists));
mu_t1 = sum(sum(repmat(helperVec,256,1).*hists));
p_0 = zeros(256);
mu_i = p_0;
mu_j = p_0;
for ii = 1:256
    for jj = 1:256
        if jj == 1
            if ii == 1
                p_0(1,1) = hists(1,1);
            else
                p_0(ii,1) = p_0(ii-1,1) + hists(ii,1);
                mu_i(ii,1) = mu_i(ii-1,1)+(ii-1)*hists(ii,1);
                mu_j(ii,1) = mu_j(ii-1,1);
            end
        else
            p_0(ii,jj) = p_0(ii,jj-1)+p_0(ii-1,jj)-p_0(ii-1,jj-1)+hists(ii,jj); % THERE IS A BUG HERE. INDICES IN MATLAB MUST BE HIGHER THAN 0. ii-1 is not valid
            mu_i(ii,jj) = mu_i(ii,jj-1)+mu_i(ii-1,jj)-mu_i(ii-1,jj-1)+(ii-1)*hists(ii,jj);
            mu_j(ii,jj) = mu_j(ii,jj-1)+mu_j(ii-1,jj)-mu_j(ii-1,jj-1)+(jj-1)*hists(ii,jj);
        end

        if (p_0(ii,jj) == 0)
            continue;
        end
        if (p_0(ii,jj) == total)
            break;
        end
        tr = ((mu_i(ii,jj)-p_0(ii,jj)*mu_t0)^2 + (mu_j(ii,jj)-p_0(ii,jj)*mu_t1)^2)/(p_0(ii,jj)*(1-p_0(ii,jj)));

        if ( tr >= maximum )
            threshold = ii;
            maximum = tr;
        end
    end
end
end

When the levels of gray of the classes of the image can be considered as Normal distributions but with unequal size and/or unequal variances, assumptions for the Otsu algorithm are not met. The Kittler-Illingworth algorithm (also known as Minimum Error thresholding)[11] is a variation of Otsu's method to handle such cases. There are several ways to mathematically describe this algorithm. One of them is to consider that for each threshold being tested, the parameters of the Normal distributions in the resulting binary image are estimated by Maximum likelihood estimation given the data.[9]

While this algorithm could seem superior to Otsu's method, it introduces new parameters to be estimated, and this can result in the algorithm being over-parametrized and thus unstable. In many cases where the assumptions from Otsu's method seem at least partially valid, it may be preferable to favor Otsu's method over the Kittler-Illingworth algorithm, following Occam's razor.[9]

Triclass thresholding tentatively divides histogram of an image into three classes, with the TBD class to be processed at next iterations.

One limitation of the Otsu’s method is that it cannot segment weak objects as the method searches for a single threshold to separate an image into two classes, namely, foreground and background, in one shot. Because the Otsu’s method looks to segment an image with one threshold, it tends to bias toward the class with the large variance.[14] Iterative triclass thresholding algorithm is a variation of the Otsu’s method to circumvent this limitation.[15] Given an image, at the first iteration, the triclass thresholding algorithm calculates a threshold ${\displaystyle \eta _{1}}$ using the Otsu’s method. Based on threshold ${\displaystyle \eta _{1}}$, the algorithm calculates mean ${\displaystyle \mu _{upper}^{[1]}}$ of pixels above ${\displaystyle \eta _{1}}$ and mean ${\displaystyle \mu _{lower}^{[1]}}$ of pixels below ${\displaystyle \eta _{1}}$. Then the algorithm tentatively separates the image into three classes (hence the name triclass), with the pixels above the upper mean ${\displaystyle \mu _{upper}^{[1]}}$ designated as the temporary foreground ${\displaystyle F}$ class and pixels below the lower mean ${\displaystyle \mu _{lower}^{[1]}}$ designated as the temporary background ${\displaystyle B}$ class. Pixels fall between ${\displaystyle [\mu _{lower}^{[1]},\mu _{upper}^{[1]}]}$ are denoted as a to-be-determined (TBD) region. This completes the first iteration of the algorithm. For the second iteration, the Otsu’s method is applied to the TBD region only to obtain a new threshold ${\displaystyle \eta _{2}}$. The algorithm then calculates the mean ${\displaystyle \mu _{upper}^{[2]}}$ of pixels in the TBD region that are above ${\displaystyle \eta _{2}}$ and the mean ${\displaystyle \mu _{lower}^{[2]}}$ of pixels in the TBD region that are below ${\displaystyle \eta _{2}}$. Pixels in the TBD region that are greater than the upper mean ${\displaystyle \mu _{upper}^{[2]}}$ are added to the temporary foreground ${\displaystyle F}$. And pixels in the TBD region that are less than the lower mean ${\displaystyle \mu _{lower}^{[2]}}$ are added to the temporary background ${\displaystyle B}$. Similarly, a new TBD region is obtained, which contains all the pixels falling between ${\displaystyle [\mu _{lower}^{[2]},\mu _{upper}^{[2]}]}$. This completes the second iteration. The algorithm then proceeds to the next iteration to process the new TBD region until it meets the stopping criterion. The criterion is that, when the difference between Otsu’s thresholds computed from two consecutive iterations is less than a small number, the iteration shall stop. For the last iteration, pixels above ${\displaystyle \eta _{n}}$ are assigned to the foreground class and pixels below the threshold are assigned to the background class. At the end, all the temporary foreground pixels are combined to constitute the final foreground. All the temporary background pixels are combined to become the final background. In implementation, the algorithm involves no parameter except for the stopping criterion in terminating the iterations. By iteratively applying the Otsu’s method and gradually shrinking the TBD region for segmentation, the algorithm can obtain a result that preserves weak objects better than the standard Otsu’s method does.