“This is the 17th day of my participation in the Gwen Challenge in November. Check out the details: The Last Gwen Challenge in 2021.”
preface
Hello! Friend!!!
Thank you very much for reading haihong’s article, if there are any mistakes in the article, please point out ~
Self-introduction ଘ(੭, ᵕ)੭
Nickname: Haihong
Tag: programmer monkey | C++ contestant | student
Introduction: because of C language to get acquainted with programming, then transferred to the computer major, had the honor to get some national awards, provincial awards… Has been confirmed. Currently learning C++/Linux/Python
Learning experience: solid foundation + more notes + more code + more thinking + learn English well!
concept
OTSU algorithm is an efficient algorithm for image binarization proposed by Japanese scholar OTSU in 1979. (Otsu algorithm)
Otsu principle
For image T (x,y), the segmentation threshold of foreground (i.e. target) and background is denoted as T, the proportion of foreground pixel points to the whole image is denoted as ω0, and the average gray level is μ0. The ratio of background pixel points to the whole image is ω1, and the average gray level is μ1. The average gray level of the whole image was denoted as μ, and the inter-class variance was denoted as G. Assume that the image size is M×N, the number of pixels in the image whose gray value is less than the threshold T is N0, and the number of pixels whose gray value is greater than the threshold T is N1
Use Otsu
Using Otsu algorithm, we can get a threshold value, which can be used for image binarization and other operations. Compared with fixed threshold with single threshold, OTSU algorithm performs better.
MATLAB Otsu algorithm is realized by garyThresh () function, generally used with IM2BW ()
Ex. :
t=rgb2gray(imread('a1.jpg')); x=graythresh(t); % Obtain the threshold t=im2bw(t,x) calculated by OTSU;Copy the codeGraythresh () source – MATLAB
function [level em] = graythresh(I)
%GRAYTHRESH Global image threshold using Otsu's method. % LEVEL = GRAYTHRESH(I) computes a global threshold (LEVEL) that can be % used to convert an intensity image to a binary image with IM2BW. LEVEL % is a normalized intensity value that lies in the range [0, 1]. % GRAYTHRESH uses Otsu's method, which chooses the threshold to minimize
% the intraclass variance of the thresholded black and white pixels.
%
% [LEVEL EM] = GRAYTHRESH(I) returns effectiveness metric, EM, as the
% second output argument. It indicates the effectiveness of thresholding
% of the input image and it is in the range [0.1]. The lower bound is
% attainable only by images having a single gray level, and the upper
% bound is attainable only by two-valued images.
%
% Class Support
% -------------
% The input image I can be uint8, uint16, int16, single, or double.and it
% must be nonsparse. LEVEL and EM are double scalars.
%
% Example
% -------
% I = imread('coins.png');
% level = graythresh(I);
% BW = im2bw(I,level);
% figure, imshow(BW)
%
narginchk(1.1);
validateattributes(I,{'uint8'.'uint16'.'double'.'single'.'int16'}, {'nonsparse'},... mfilename,'I'.1);
if ~isempty(I)
% Convert all N-D arrays into a single column. Convert to uint8 for
% fastest histogram computation.
I = im2uint8(I(:));
num_bins = 256;
counts = imhist(I,num_bins);
% Variables names are chosen to be similar to the formulas in
% the Otsu paper.
p = counts / sum(counts);
omega = cumsum(p);
mu = cumsum(p .* (1:num_bins)'); mu_t = mu(end); sigma_b_squared = (mu_t * omega - mu).^2 ./ (omega .* (1 - omega)); % Find the location of the maximum value of sigma_b_squared. % The maximum may extend over several bins, so average together the % locations. If maxval is NaN, meaning that sigma_b_squared is all NaN, % then return 0. maxval = max(sigma_b_squared); isfinite_maxval = isfinite(maxval); if isfinite_maxval idx = mean(find(sigma_b_squared == maxval)); % Normalize the threshold to the range [0, 1]. level = (idx - 1) / (num_bins - 1); The else level = 0.0; End else level = 0.0; isfinite_maxval = false; end % compute the effectiveness metric if nargout > 1 if isfinite_maxval em = maxval/(sum(p.*((1:num_bins).^2)') - mu_t^2);
else
em = 0;
end
end
Copy the codeOstu algorithm (personal implementation MATLAB version)
t=rgb2gray(imread('a1.jpg')); [m,n]=size(t); %counts =m*n; % count is one256The matrix in one row and one column records the number of pixels on each gray level count=imhist(t); % Ratio of pixels on each gray level to total pixels (probability) p=count/counts; %cumsum w0=cumsum(p); %u calculates the average gray level. For example, the gray level is0~125Average gray value %(1:256)Cumsum (p.*(1:256)) cumsum(p.*(1:256)) cumsum(p.*(1:256))); %u_end is the global average gray level u_end=u(end); %d is the variance d is 0256The matrix d=(w0*u_end-u).^2./(w0.*(1-w0)); % find the gray value of maximum variance in d, where y is the gray value of maximum variance [x,y]= Max (d); % for normalization with im2BW %x is the threshold x=(y- 1) /255
Copy the codeValidation of algorithm results
t=rgb2gray(imread('a1.jpg'));
[m,n]=size(t);
counts=m*n;
count=imhist(t);
p=count/counts;
w1=cumsum(p);
u=cumsum(p.*(1:256)'); u_end=u(end); %u_end is the global average grayscale d=(w1*u_end-u).^2./(w1.*(1-w1)); [x,y]=max(d); X =(y-1)/255 %x is the threshold xx= Graythresh (t)%graythresh calculated thresholdCopy the codeThe results of
% The threshold x = calculated by the code I wrote0.7569% the threshold xx = calculated using garythresh0.7569
Copy the codeConclusion This code is roughly identical to the threshold calculated by Gragthresh ()!
Suspects to reassure
In our own code, we calculate the variance using the following formula, which is different from the formula in the principle (g=ω0ω1(μ0-μ1)^) is the principle wrong? Ha ha no, in fact, it is a formula, the formula in the code is just a deformation of the principle formula, this is to reduce the variables, improve the speed of operation.
%d is the variance d is 0256The matrix d=(w1*u_end-u).^2./(w1.*(1-w1));
Copy the codeFormula of deformation
G =w0w1(u1-u0)^2; g=w0w1(u1-u0)^2; g=w0w1(u1-u0)^2;
C=imread('h.jpg'); C=rgb2gray(C); Figure,imshow(C); title('Raw grayscale image'); Count =imhist(C); % histogram statistics subplot(1.3.1); imhist(C); % draw histogram [r,t]=size(C); % image matrix size N=r*t; % Number of image pixels L=256; % set the highlight gray level as256count=count/N; % probability of occurrence of gray levelsfor i=2:L
if count(i)~=0 ;
st=i- 1;
breakEnd end % above the loop statement to find the probability of occurrence is not0Is the minimum gray value offor i=L:- 1:1
if count(i)~=0;
nd=i- 1;
breakEnd end % to find the probability of occurrence is not0F =count(st+)1:nd+1); p=st; q=nd-st; %p and q are the starting and ending values u= of gray scale respectively0;
for i=1:q
u=u+f(i)*(p+i- 1); ua(i)=u; End % calculates the average gray value of the imagefor i=1:q
w(i)=sum(f(1:i)); End % calculates the probability of A region d=(u*w-ua).^2./(w.*(1-w)); % Calculate the variance between classes with different k values [y, TP]= Max (d); % the gray value corresponding to the maximum variance th= TP + P;if th<=160
th=tp+p;
else
th=160End % Modify threshold Y1=zeros(r,t);for i=1:r
for j=1:t
X1(i,j)=double(C(i,j));
end
end
for i=1:r
for j=1:t
if (X1(i,j)>=th)
Y1(i,j)=X1(i,j);
else Y1(i,j)=0; Segment figure,imshow(Y1); title('Gray threshold segmentation image');
Copy the code