Otsu's method

Let, H be the normalised histogram of the pixels in an image (s.t. it becomes the probability distribution of pixel intensities) with L bins. There are two classes of this histogram: C_0 for background pixels, and C_1 for foreground pixels. The primary disciminator of pixels (to assort them into classes) is the threshold t. C_0 includes pixels from 0 to (t-1) , and C_1 includes from t to (L-1) . The algorithm is then global search for an optimal threshold t^* such that intra-class variance (variance of pixels intensities in C_0 or C_1 ) is minimised. Let, \omega_0 denote the cumulative probability of C_0 , and \omega_1denote of C_1 . \begin{align} \omega_0(t) & =\sum_{i=0}^{t-1} P(i), \\ \omega_1(t) & =\sum_{i=t}^{L-1} P(i). \end{align} For a classes C_0 and C_1 , the conditional probability of selecting the i -th pixel in those classes is P(i | C_0) and P(i | C_1) respectively. Now, let \mu_0(t) and \mu_1(t) be the mean (pixel intensity) of C_0 and C_1 respectively. \begin{align} \mu_0(t) &= \sum^{t-1}_{i=0} iP(i | C_0) = \sum^{t-1}_{i=0} \frac{iP(i)}{\omega_0(t)} = \frac{\sum^{t-1}_{i=0} iP(i)}{\omega_0(t)}. \end{align} Similarly, \mu_1(t) = \frac{\sum^{L-1}_{i=t} iP(i)}{\omega_1(t)}. Now, let \sigma^2_0(t) and \sigma^2_1(t) be the (pixel intensity) variance of C_0 and C_1 respectively. \begin{align} \sigma_0^2(t) &= \sum^{t-1}_{i=0} (i - \mu_0)^2 P(i | C_0) = \sum^{t-1}_{i=0} \frac{(i - \mu_0)^2 P(i)}{\omega_0} = \frac{\sum^{t-1}_{i=0} (i - \mu_0)^2P(i)}{\omega_0(t)}. \end{align} Similarly, \sigma_1^2(t) = \frac{\sum^{L-1}_{i=t} (i - \mu_1)^2 P(i)}{\omega_1(t)}. Let, \sigma_b^2(t) be the inter-class (pixel intensity) variance, which is defined as the weighted sum of variances of aforementioned two classes. \begin{align} \sigma^2_b(t) &= \sigma^2_T - \left[\omega_0(t) \sigma^2_0(t) + \omega_1(t) \sigma^2_1(t)\right] \\ &= \omega_0(\mu_0 - \mu_T)^2 + w_1 (\mu_1 - \mu_T)^2 \\ &= \omega_0\omega_1(\mu_0 - \mu_1)^2. \end{align} Where, \sigma^2_T(t) variance of the total histogram. {{Math proof|proof=Considering \omega_0 + \omega_1 = 1 and \omega_0\mu_0 + \omega_1\mu_1 = \mu_T, we can prove the following. \begin{align} \sigma^2_b(t) &= \omega_0(\mu_0 - \mu_T)^2 + w_1 (\mu_1 - \mu_T)^2 \\ &= \omega_0\mu_0^2 + \omega_1\mu_1^2 - 2 \mu_T(\omega_0\mu_0 + \omega_1\mu_1) + \mu_T^2(\omega_0 + \omega_1) \\ &= \omega_0\mu_0^2 + \omega_1\mu_1^2 - 2 \mu_T^2 + \mu_T^2 \\ &= \omega_0\mu_0^2 + \omega_1\mu_1^2 - \mu_T^2 \\ &= \omega_0\mu_0^2 + \omega_1\mu_1^2 - (\omega_0\mu_0 + \omega_1\mu_1)^2 \\ &= \omega_0\mu_0^2 - \omega_0^2\mu_0^2 + \omega_1\mu_1^2 - \omega_1^2\mu_1^2 - 2\omega_0\mu_0\omega_1\mu_1 \\ &= \omega_0\mu_0^2 (1-\omega_0) + \omega_1\mu_1^2(1-\omega_1) - 2\omega_0\mu_0\omega_1\mu_1 \\ &= \omega_0\mu_0^2 (1-\omega_0) - \omega_0\mu_0\omega_1\mu_1 + \omega_1\mu_1^2(1-\omega_1) - \omega_0\mu_0\omega_1\mu_1 \\ &= \omega_0\omega_1\mu_0^2 - \omega_0\omega_1\mu_0\mu_1 + \omega_0\omega_1\mu_1^2 - \omega_0\omega_1\mu_0\mu_1 \\ &= \omega_0\omega_1(\mu_0^2 - \mu_0\mu_1) + \omega_0\omega_1(\mu_1^2 - \mu_0\mu_1) \\ &= \omega_0\omega_1(\mu_0^2 - 2\mu_0\mu_1 + \mu_1^2) \\ &= \omega_0\omega_1(\mu_0 - \mu_1)^2. \end{align}}} The algorithm is now to maximise \sigma_b^2(t), i.e. inter-class variance. This standpoint is motivated by a conjecture that well-thresholded classes would be separated in pixel intensities, and conversely a threshold t^* giving the best separation of classes in pixel intensities would be the best threshold. Formally, this problem is summarised as the following. {{Equation box 1|equation=\sigma^2_b(t^*) = \max_{0 }} Algorithm • Compute histogram and probabilities of each intensity level. • Set up initial \omega_0(0), \mu_0(0) and \omega_1(0) and \mu_1(0). • Step through all possible thresholds from t = 1 to maximum intensity. • Update \omega_0(0), \mu_0(0) and \omega_1(0) and \mu_1(0). • Compute \sigma^2_b(t). • Desired threshold t^* corresponds to the maximum \sigma^2_b(t). MATLAB implementation 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 wB = wB + histogramCounts(ii); wF = total - wB; sumB = sumB + (ii-1) * histogramCounts(ii); 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 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 This implementation requires the NumPy library. import numpy as np def otsu_intraclass_variance(image, threshold): """ Otsu's intra-class variance. If all pixels are above or below the threshold, this will throw a warning that can safely be ignored. """ return np.nansum( [ np.mean(cls) * np.var(image, where=cls) # weight · intra-class variance for cls in [image >= threshold, image Python libraries dedicated to image processing such as OpenCV and Scikit-image provide built-in implementations of the algorithm. ==Limitations and variations==