MATLAB image processing - Otsu threshold segmentation (with code)

Catalog

Otsu threshold

Algorithm flow

The flow chart shows

result

Code

---

Otsu threshold

         Otsu threshold is also called maximum variance threshold , yes 1979 Proposed by Japan's Otsu exhibition in , It is derived from the principle of discrimination and least square method , The basic idea is to maximize the variance between classes , The optimal threshold is obtained .

Algorithm flow

         Make Indicates a picture with the size of Pixels in a digital image Different gray levels , Indicates that the gray level is The number of pixels , The total number of pixels in the image is . The normalized histogram has components , There is

Now? , Suppose you choose a threshold , And use it to thresholding the input image into two types and , among The gray value in the image is in the interval All pixels in the , By the gray value in the interval All pixels in the . With this threshold , Pixels are classified into classes The probability of   Given by the following cumulative sum ：

Change the angle , It's kind of Probability of occurrence . Allied , class The probability of that happening is zero ：

Assign to class The average gray value of is ：

Allied , Assign to class The average gray value of is ：

To the first k The cumulative mean value of level （ Average grayscale ） Given by the following formula ：

The average gray level of the whole image is given by the following formula ：

To evaluate the segmentation quality of the selected threshold , Define the variance between classes σ2 by ：

The above formula can be reduced to ：

As can be seen from the above formula , Calculate different thresholds k Between class variance , You just have to calculate m and Two parameters . To get the optimal threshold , from Choose different k Calculate the variance between classes , When Obtained at maximum k Is the maximum variance threshold .

The flow chart shows

result

          As you can see from the diagram , After Otsu threshold segmentation, we can basically get two parts , You can extract parts of the moon from the picture . Can quickly and effectively find the segmentation threshold between classes , But its shortcomings are also obvious , That is, it can only be divided for a single target , Or the objects of interest belong to the same gray range .

Code

clear; clc;
I=rgb2gray(imread('moon.jpg'));
subplot(1, 2, 1)
imshow(I);
xlabel('(a)  original image ');
% level = graythresh(I);      % Use MATLAB  Function to calculate the threshold 
% BW = im2bw(I, level);     
% subplot(1, 3, 2)
% imshow(BW);
% xlabel('(b) graythresh');
% disp(['graythresh  Calculate the gray threshold :', num2str(level*255)]);
T = Otsu(double(I));     % Use Otsu method to calculate the threshold 
disp([' Otsu method to calculate the gray threshold :', num2str(T)])
BW = im2bw(I, T/255);
% Threshold segmentation 
subplot(1, 2, 2)
imshow(BW);
xlabel('(c)  Dajin law ');

function ThreshValue = Otsu(Imag)
%  Otsu method to calculate the threshold 
%  Input ：
%    Imag： Two dimensional array , Numerical values represent grayscale ;
%  Output ：
%    ThreshValue： threshold 
iMax = max(Imag(:));              %  Maximum 
iMin = min(Imag(:));               %  minimum value 
T = iMin:iMax;                        %  Gray value range 
Tval = zeros(size(T));               %  variance 
[iRow, iCol] = size(Imag);        %  Data dimension size 
imagSize = iRow*iCol;            %  Number of pixels 
%  Traversal gray value , Calculate variance 
for i = 1 : length(T)
    TK = T(i);
    iFg = 0;          %  prospects 
    iBg = 0;          %  background 
    FgSum = 0;    %  Total prospects 
    BgSum = 0;    %  Total number of backgrounds 
    for j = 1 : iRow
        for k = 1 : iCol
            temp = Imag(j, k);
            if temp > TK
                iFg = iFg + 1;      %  Foreground pixel statistics 
                FgSum = FgSum + temp;
            else
                iBg = iBg + 1;      %  Background pixel statistics 
                BgSum = BgSum + temp;
            end
        end
    end
    w0 = iFg/imagSize;      %  Foreground ratio 
    w1 = iBg/imagSize;     %  Background proportion 
    u0 = FgSum/iFg;         %  Average value of foreground gray level 
    u1 = BgSum/iBg;        %  Average value of background gray level 
    Tval(i) = w0*w1*(u0 - u1)*(u0 - u1);     %  Calculate variance 
end
[~, flag] = max(Tval);             %  The maximum subscript 
ThreshValue = T(flag);
end