{"id":1949,"date":"2016-06-14T17:35:19","date_gmt":"2016-06-14T21:35:19","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=1949"},"modified":"2019-11-01T16:43:15","modified_gmt":"2019-11-01T20:43:15","slug":"image-binarization-otsus-method","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2016\/06\/14\/image-binarization-otsus-method\/","title":{"rendered":"Image binarization – Otsu’s method"},"content":{"rendered":"

In my 16-May-2016 post<\/a> about image binarization, I talked about the new binarization functions in R2016a. Today I want to switch gears and talk about Otsu's method, one of the algorithms underlying imbinarize<\/tt><\/a>.<\/p>

(A bonus feature of today's blog post is a demo of yyaxis<\/tt><\/a>, a new feature of MATLAB R2016a.)<\/p>

Otsu's method is named for Nobuyuki Otsu, who published it in IEEE Transactions on Systems, Man, and Cybernetics<\/i>, vol. SMC-9, no. 1, January 1979. At this time, researchers had already explored a variety of ways to choose a threshold automatically by examining the histogram of image pixel values. The basic idea is to look for two peaks, representing foreground and background pixel values, and pick a point in between the two peaks as the threshold value.<\/p>

Here's a simple example using the coins image.<\/p>

I = imread('coins.png'<\/span>);\r\nimshow(I)\r\n<\/pre>\"\" 
imhist(I)\r\n<\/pre>\"\" 
The function imbinarize<\/tt> calls otsuthresh<\/tt> to get a normalized threshold value.<\/p>
t = otsuthresh(histcounts(I,-0.5:255.5))\r\n<\/pre>
\r\nt =\r\n\r\n    0.4941\r\n\r\n<\/pre>
Let's see where that threshold is.<\/p>
hold on<\/span>\r\nplot(255*[t t], ylim, 'r'<\/span>, 'LineWidth'<\/span>, 5)\r\nhold off<\/span>\r\n<\/pre>\"\" 
And here is the thresholded coins image.<\/p>
imshow(imbinarize(I,t))\r\n<\/pre>\"\" 
How does this threshold selection work? It is based on entirely on the set of histogram counts. To show the computation, I'll adopt the notation from the paper. Pixels can take on the set of values $i = 1,2,\\ldots,L$. The histogram count for pixel value $i$ is $n_i$, and the associated probability is $p_i = n_i\/N$, where $N$ is the number of image pixels. (I'm using the word probability<\/i> here somewhat loosely, in the relative frequency sense.)<\/p>
The thresholding task is formulated as the problem of dividing image pixels into two classes. $C_0$ is the set of pixels with values $[1,\\ldots,k]$, and $C_1$ is the set of pixels with values in the range $[k+1,\\ldots,L]$.<\/p>
The overall class probabilities, $\\omega_0$ and $\\omega_1$, are:<\/p>
$$\\omega_0 = \\sum_{i=1}^k p_i = \\omega(k)$$<\/p>
$$\\omega_1 = \\sum_{i=k+1}^L p_i = 1 - \\omega_0(k)$$<\/p>
The class means, $\\mu_0$ and $\\mu_1$, are the mean values of the pixels in $C_0$ and $C_1$. They are given by:<\/p>
$$\\mu_0 = \\sum_{i=1}^k i p_i \/ \\omega_0 = \\mu(k)\/\\omega(k)$$<\/p>
$$\\mu_1 = \\sum_{i=k+1}^L i p_i \/ \\omega_1 = \\frac{\\mu_T - \\mu(k)}{1 -\r\n\\omega(k)}$$<\/p>
where<\/p>
$$\\mu(k) = \\sum_{i-1}^k i p_i$$<\/p>
and $\\mu_T$, the mean pixel value for the total image, is:<\/p>
$$\\mu_T = \\sum_{i=1}^L i p_i.$$<\/p>
The class variances, $\\sigma_0^2$ and $\\sigma_1^2$, are:<\/p>
$$\\sigma_0^2 = \\sum_{i = 1}^k (i - \\mu_0)^2 p_i \/ \\omega_0$$<\/p>
$$\\sigma_1^2 = \\sum_{i = k+1}^L (i - \\mu_1)^2 pi \/ \\omega_1.$$<\/p>
Otsu mentions three measures of \"good\" class separability: within-class variance ($\\lambda$), between-class variance ($\\kappa$), and total variance ($\\eta$). These are given by:<\/p>
$$\\lambda = \\sigma_B^2$$<\/p>
$$\\kappa = \\sigma_T^2\/\\sigma_W^2$$<\/p>
$$\\eta = \\sigma_B^2\/\\sigma_T^2$$<\/p>
where<\/p>
$$\\sigma_W^2 = \\omega_0 \\sigma_0^2 + \\omega_1 \\sigma_1^2$$<\/p>
$$\\sigma_B^2 = \\omega_0 (\\mu_0 - \\mu_T)^2 + \\omega_1 (\\mu_1 - \\mu_T)^2 =\r\n\\omega_0 \\omega_1 (\\mu_1 - \\mu_0)^2.$$<\/p>
He goes on to point out that maximizing any of these criteria is equivalent to maximizing the others. Further, maximizing $\\eta$ is the same as maximizing $\\sigma_B^2$, which can be rewritten in terms of the selected threshold, $k$:<\/p>
$$ \\sigma_B^2(k) = \\frac{[\\mu_T \\omega(k) - \\mu(k)]^2}{\\omega(k) [1 -\r\n\\omega(k)]}.$$<\/p>
The equation above is the heart of the algorithm. $\\sigma_B^2$ is computed for all possible threshold values, and we choose as our threshold the value that maximizes it.<\/p>
OK, that was a lot of equations, but there's really not that much involved in computing the key quantity, $\\sigma_B^2(k)$. Here's what the computation looks like for the coins image.<\/p>
counts = imhist(I);\r\nL = length(counts);\r\np = counts \/ sum(counts);\r\nomega = cumsum(p);\r\nmu = cumsum(p .* (1:L)');\r\nmu_t = mu(end);\r\n\r\nsigma_b_squared = (mu_t * omega - mu).^2 .\/ (omega .* (1 - omega));\r\n<\/pre>
Using yyaxis<\/tt>, a new R2016a feature, let's plot the histogram and $\\sigma_B^2$ together.<\/p>
close all<\/span>\r\nyyaxis left<\/span>\r\nplot(counts)\r\nylabel('Histogram'<\/span>)\r\n\r\nyyaxis right<\/span>\r\nplot(sigma_b_squared)\r\nylabel('\\sigma_B^2'<\/span>)\r\n\r\nxlim([1 256])\r\n<\/pre>\"\" 
Otsu's method chooses the place where $\\sigma_B^2$ is the highest as the threshold.<\/p>
[~,k] = max(sigma_b_squared);\r\nhold on<\/span>\r\nplot([k k],ylim,'LineWidth'<\/span>,5)\r\nhold off<\/span>\r\n<\/pre>\"\" 
Here's another example. This is a public-domain light microscope image of Lily mitosis. (The original image<\/a> is courtesy Andrew S. Bajer, University of Oregon, Eugene, OR. This version is slightly cropped.)<\/p>
url = 'https:\/\/blogs.mathworks.com\/steve\/files\/205.jpg'<\/span>;\r\nI = rgb2gray(imread(url));\r\nclf\r\nimshow(I)\r\n<\/pre>\"\" 
counts = imhist(I);\r\nL = length(counts);\r\np = counts \/ sum(counts);\r\nomega = cumsum(p);\r\nmu = cumsum(p .* (1:L)');\r\nmu_t = mu(end);\r\n\r\nsigma_b_squared = (mu_t * omega - mu).^2 .\/ (omega .* (1 - omega));\r\n\r\nclose all<\/span>\r\nyyaxis left<\/span>\r\nplot(counts)\r\nylabel('Histogram'<\/span>)\r\n\r\nyyaxis right<\/span>\r\nplot(sigma_b_squared)\r\nylabel('\\sigma_B^2'<\/span>)\r\n\r\n[~,k] = max(sigma_b_squared);\r\nhold on<\/span>\r\nplot([k k],ylim,'LineWidth'<\/span>,5)\r\nhold off<\/span>\r\n\r\nxlim([1 256])\r\n<\/pre>\"\" 
clf\r\nimshow(imbinarize(I))\r\ntitle('Thresholded cell image'<\/span>)\r\n<\/pre>\"\" 
If imbinarize<\/tt> handles this computation automatically, then why did we also provide a function called otsuthresh<\/tt><\/a>? The answer is that imbinarize<\/tt> takes an image as input, although Otsu's method does not require the original image, only the image's histogram. If you have a situation where you want to compute a threshold based only<\/b> on a histogram, then you can call otsuthresh<\/tt> directly. That's why it is there.<\/p>
To wrap up this week's discussion, I want to point out that a couple of blog readers recommended something called the Triangle method<\/i> for automatic gray-scale image thresholding. If you want to try this for yourself, there is an implementation on the File Exchange. I have not had a chance yet to experiment with it.<\/p>
Next time I'll talk about the algorithm used by imbinarize<\/tt><\/a> for locally adaptive thresholding.<\/p>