Otsu's method
In computer vision and image processing, Otsu's method, named after , is used to perform automatic image thresholding. 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. 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 performed on the intensity histogram. The extension to multi-level thresholding was described  in the original paper, and computationally efficient implementations have since been proposed. Otsu's method. 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 formula_2 and formula_3 are the probabilities of the two classes separated by a threshold formula_4 ,and formula_5 and formula_6 are variances of these two classes. The class probability formula_7 is computed from the formula_8 bins of the histogram: For 2 classes, minimizing the intra-class variance is equivalent to maximizing inter-class variance: which is expressed in terms of class  probabilities formula_11 and class means formula_12, where the class means formula_13, formula_14 and formula_15 are: The following relations can be easily verified: The class probabilities and class means can be computed iteratively. This idea yields an effective algorithm. MATLAB or Octave implementation. codice_1 is a 256-element histogram of a grayscale image different gray-levels (typical for 8-bit images). codice_2 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 end end Matlab has built-in functions codice_3 and codice_4 in the Image Processing Toolbox which are implemented with Otsu's method and Multi Otsu's method, respectively. Python implementation. This implementation requires the NumPy library. def compute_otsu_criteria(im, th): im = # load your image as a numpy array. threshold_range = range(np.max(im)+1) criterias = [compute_otsu_criteria(im, th) for th in threshold_range] 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. Otsu's method performs well when the histogram  has a bimodal distribution with a deep and sharp valley between the two peaks. 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. In those cases, local adaptations of the Otsu method have been developped. 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. 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) 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, such as the Kittler-Illingworth method. A variation for noisy images. Various extensions have been developed to address limitations of Otsu's method. One popular extension 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. At each pixel, the average gray-level value of the neighborhood is calculated. Let the gray level of the given pixel be divided into formula_8 discrete values and the average gray level is also divided into the same formula_8 values. Then a pair is formed: the pixel gray level and the average of the neighborhood formula_27. Each pair belongs to one of the formula_28 possible 2-dimensional bins. The total number of occurrences (frequency), formula_29, of a pair formula_27, divided by the total number of pixels in the image formula_31, 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 formula_41 is obtained  by maximizing formula_42. Algorithm. The formula_43 and formula_4 is obtained iteratively which is similar with one-dimensional Otsu's method. The values of formula_43 and formula_4 are changed till we obtain the maximum of formula_42, that is max,s,t = 0; for ss: 0 to L-1 do end for return s,t; Notice that for evaluating formula_42, we can use a fast recursive dynamic programming algorithm to improve time performance. 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. If summed  area tables are used to build the 3 tables, sum over formula_49, sum over formula_50, and sum over formula_51, 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. function inputs and output: codice_5 is a formula_52 2D-histogram of grayscale value and neighborhood average grayscale value pair. codice_6 is the number of pairs in the given image.it is determined by the number of the bins of 2D-histogram at each direction. codice_7 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 end end A variation for unbalanced images. 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) 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. 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.