Otsu's method

From Wikipedia, the free encyclopedia
Jump to navigation Jump to search
An example image thresholded using Otsu's algorithm
Original image

In computer vision and image processing, Otsu's method, named after Nobuyuki Otsu (大津展之, Ōtsu Nobuyuki), is used to perform automatic image thresholding.[1] In the simplest form, the algorithm returns a single intensity threshold that separate pixels into two classes, foreground and background. This threshold is determined by minimizing intra-class intensity variance, or equivalently, by maximizing inter-class variance.[2] Otsu's method is a one-dimensional discrete analog of Fisher's Discriminant Analysis, is related to Jenks optimization method, and is equivalent to a globally optimal k-means[3] performed on the intensity histogram. The extension to multi-level thresholding was described in the original paper,[2] and computationally efficient implementations have since been proposed.[4][5]

Otsu's method[edit]

Otsu's method visualization

The algorithm exhaustively searches for the threshold that minimizes the intra-class variance, defined as a weighted sum of variances of the two classes:

Weights and are the probabilities of the two classes separated by a threshold ,and and are variances of these two classes.

The class probability is computed from the bins of the histogram:

For 2 classes, minimizing the intra-class variance is equivalent to maximizing inter-class variance:[2]

which is expressed in terms of class probabilities and class means , where the class means , and are:

The following relations can be easily verified:

The class probabilities and class means can be computed iteratively. This idea yields an effective algorithm.

Algorithm[edit]

  1. Compute histogram and probabilities of each intensity level
  2. Set up initial and
  3. Step through all possible thresholds maximum intensity
    1. Update and
    2. Compute
  4. Desired threshold corresponds to the maximum

MATLAB implementation[edit]

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.

Python implementation[edit]

This implementation requires the NumPy library.

import numpy as np

def compute_otsu_criteria(im, th):
    # 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 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
    var0 = np.var(val_pixels0) if len(val_pixels0) > 0 else 0
    var1 = np.var(val_pixels1) if len(val_pixels1) > 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.

Limitations and variations[edit]

Otsu's method performs well when the histogram has a bimodal distribution with a deep and sharp valley between the two peaks.[6]

Like all other global thresholding methods, Otsu's method performs badly in case of heavy noise, small objects size, inhomogeneous lighting and larger intra-class than inter-class variance.[7] In those cases, local adaptations of the Otsu method have been developped.[8]

Moreover, the mathematical grounding of Otsu's method models the histogram of the image as a mixture of two Normal distributions with equal variance and equal size.[9] Otsu's thresholding may however yield satisfying results even when these assumptions are not met, in the same way statistical tests (to which Otsu's method is heavily connected[10]) can perform correctly even when the working assumptions are not fully satsified. However, several variations of Otsu's methods have been proposed to account for more severe deviations from these assumptions,[9] such as the Kittler-Illingworth method.[11]

A variation for noisy images[edit]

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 discrete values and the average gray level is also divided into the same values. Then a pair is formed: the pixel gray level and the average of the neighborhood . Each pair belongs to one of the possible 2-dimensional bins. The total number of occurrences (frequency), , of a pair , divided by the total number of pixels in the image , defines the joint probability mass function in a 2-dimensional histogram:

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:

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

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

The inter-class discrete matrix is defined as

The trace of the discrete matrix can be expressed as

where

Similar to one-dimensional Otsu's method, the optimal threshold is obtained by maximizing .

Algorithm[edit]

The and is obtained iteratively which is similar with one-dimensional Otsu's method. The values of and are changed till we obtain the maximum of , 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 , 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 , sum over , and sum over , 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.

MATLAB implementation[edit]

function inputs and output:

hists is a 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

A variation for unbalanced images[edit]

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]

References[edit]

  1. ^ M. Sezgin & B. Sankur (2004). "Survey over image thresholding techniques and quantitative performance evaluation". Journal of Electronic Imaging. 13 (1): 146–165. doi:10.1117/1.1631315.
  2. ^ a b c Nobuyuki Otsu (1979). "A threshold selection method from gray-level histograms". IEEE Trans. Sys. Man. Cyber. 9 (1): 62–66. doi:10.1109/TSMC.1979.4310076.
  3. ^ Liu, Dongju (2009). "Otsu method and K-means". Ninth International Conference on Hybrid Intelligent Systems IEEE. 1: 344–349.
  4. ^ Liao, Ping-Sung (2001). "A fast algorithm for multilevel thresholding" (PDF). J. Inf. Sci. Eng. 17 (5): 713–727. Archived from the original (PDF) on 2019-06-24.
  5. ^ Huang, Deng-Yuan (2009). "Optimal multi-level thresholding using a two-stage Otsu optimization approach". Pattern Recognition Letters. 30 (3): 275–284. doi:10.1016/j.patrec.2008.10.003.
  6. ^ Kittler, J.; Illingworth, J. (September 1985). "On threshold selection using clustering criteria". IEEE Transactions on Systems, Man, and Cybernetics. SMC-15 (5): 652–655. doi:10.1109/tsmc.1985.6313443. ISSN 0018-9472.
  7. ^ Lee, Sang Uk and Chung, Seok Yoon and Park, Rae Hong (1990). "A comparative performance study of several global thresholding techniques for segmentation". Computer Vision, Graphics, and Image Processing. 52 (2): 171–190. doi:10.1016/0734-189x(90)90053-x.{{cite journal}}: CS1 maint: multiple names: authors list (link)
  8. ^ a b Jianzhuang, Liu and Wenqing, Li and Yupeng, Tian (1991). "Automatic thresholding of gray-level pictures using two-dimension Otsu method". Circuits and Systems, 1991. Conference Proceedings, China., 1991 International Conference on: 325–327.{{cite journal}}: CS1 maint: multiple names: authors list (link)
  9. ^ a b c d Kurita, T.; Otsu, N.; Abdelmalek, N. (October 1992). "Maximum likelihood thresholding based on population mixture models". Pattern Recognition. 25 (10): 1231–1240. doi:10.1016/0031-3203(92)90024-d. ISSN 0031-3203.
  10. ^ Jing-Hao Xue; Titterington, D. M. (August 2011). "$t$-Tests, $F$-Tests and Otsu's Methods for Image Thresholding". IEEE Transactions on Image Processing. 20 (8): 2392–2396. doi:10.1109/tip.2011.2114358. ISSN 1057-7149.
  11. ^ a b Kittler, J.; Illingworth, J. (1986-01-01). "Minimum error thresholding". Pattern Recognition. 19 (1): 41–47. doi:10.1016/0031-3203(86)90030-0. ISSN 0031-3203.
  12. ^ Zhang, Jun & Hu, Jinglu (2008). "Image segmentation based on 2D Otsu method with histogram analysis". Computer Science and Software Engineering, 2008 International Conference on. 6: 105–108. doi:10.1109/CSSE.2008.206. ISBN 978-0-7695-3336-0.
  13. ^ Zhu, Ningbo and Wang, Gang and Yang, Gaobo and Dai, Weiming (2009). "A fast 2d otsu thresholding algorithm based on improved histogram". Pattern Recognition, 2009. CCPR 2009. Chinese Conference on: 1–5.{{cite journal}}: CS1 maint: multiple names: authors list (link)

External links[edit]