「这是我参与11月更文挑战的第17天，活动详情查看：2021最后一次更文挑战」

前言

Hello！小伙伴！

非常感谢您阅读海轰的文章，倘若文中有错误的地方，欢迎您指出～

 

自我介绍 ଘ(੭ˊᵕˋ)੭

昵称：海轰

标签：程序猿｜C++选手｜学生

简介：因C语言结识编程，随后转入计算机专业，有幸拿过一些国奖、省奖…已保研。目前正在学习C++/Linux/Python

学习经验：扎实基础 + 多做笔记 + 多敲代码 + 多思考 + 学好英语！

概念

OTSU算法是由日本学者OTSU于1979年提出的一种对图像进行二值化的高效算法。（大津算法）

Otsu原理

对于图像 t (x,y)，前景(即目标)和背景的分割阈值记作 T，属于前景的像素点数占整幅图像的比例记为 ω0，平均灰度为 μ0；背景像素点数占整幅图像的比例为 ω1，平均灰度为 μ1；整幅图像的平均灰度记为μ，类间方差记为g。 假设图像大小为M×N，图像中像素的灰度值小于阈值 T 的像素个数为 N0，像素灰度大于阈值T的像素个数为 N1 注：（7）式是将（5）式代入（6）式得到的，我们的重点放在（7）式上。

Otsu用处

利用Otsu算法，我们可以得到一个阈值，利用该阈值对图像进行二值化等操作。相比于单阈值的固定阈值，otsu算法效果更好。

MATLAB中实现Otsu算法的是 garythresh()函数，一般都与im2bw()配套使用

例：

t=rgb2gray(imread('a1.jpg'));
x=graythresh(t);%获取otsu算得的阈值
t=im2bw(t,x);

graythresh()源码–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); 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

Ostu算法（个人实现 MATLAB版）

t=rgb2gray(imread('a1.jpg'));
[m,n]=size(t);

%counts为图片总像素个数值
counts=m*n;

%count是一个256行一列的矩阵 记录了每个灰度级上像素点的个数
count=imhist(t);

%每个灰度级上像素点占总像素点数量的比例（概率）
p=count/counts;

%cumsum 计算累加概率 
w0=cumsum(p);

%u  计算的是平均灰度  比如灰度级为0~125的平均灰度值
%(1:256)' 矩阵转置 1行256列转换为256行1列 u=cumsum(p.*(1:256)');

%u_end是全局平均灰度
u_end=u(end);

%d 求方差 d是256行一列的矩阵
d=(w0*u_end-u).^2./(w0.*(1-w0));

%在d中寻找为最大方差的灰度值 这里y是最大方差的灰度值
[x,y]=max(d);

%为了和im2bw配合使用 进行归一化
%x就是所得阈值
x=(y-1)/255

算法结果验证

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是全局平均灰度 d=(w1*u_end-u).^2./(w1.*(1-w1)); [x,y]=max(d); x=(y-1)/255 %x就是所得阈值 xx=graythresh(t)%graythresh计算出的阈值 

结果　

%自己编写的代码计算出的阈值
x =

    0.7569

%利用garythresh计算出的阈值
xx =

    0.7569

结论 该代码与gragthresh（）计算出的阈值大体上完全一致！

疑点解惑

在我们自己编写的代码中，计算方差用的是下面的公式， 与原理中计算方差的公式（ g=ω0ω1(μ0-μ1)^ ）不一样 难道是原理错了？ 哈哈 不是，其实都是一个公式，我们在代码中的公式只是原理公式的变形罢了，这是为了减少变量，提高运算速度。

%d 求方差 d是256行一列的矩阵
d=(w1*u_end-u).^2./(w1.*(1-w1));

公式变形

注意：上图中的u就是验证代码中的u_end,代表的是图像全局的平均灰度 而u是代表灰度值 T前面的平均灰度。 这里要注意图中的w0u0就是验证代码中的u，对每个灰度级数乘以对应的概率并累加，这是因为u0的计算公式是每个灰度级乘以对应的数量，再除以前面一共的像素数量，这里前面一共的像素数量可以用总量乘以前面占的概率w0。哈哈，对，就是这里，就可以把w0u0中的w0约掉，所以w0*u0==u（这里的u是验证代码中的u，代表灰度值乘以概率的累加值）

其实只需要知道g=w0w1(u1-u0)^2就行了，验证代码只是为了减少变量，做了一下变形，哎，数学还是要学好啊 下面在给个版本的验证算法，原理差不多，只是个别地方做了精简（其实现在看来，这个算法反而搞复杂了）

C=imread('h.jpg');   
C=rgb2gray(C);
%读取图像

figure,imshow(C);title('原始灰度图像');%绘原图

count=imhist(C);                       %直方图统计
subplot(1,3,1); 
imhist(C);                  %绘制直方图

[r,t]=size(C); 

%图像矩阵大小 

N=r*t;                                 %图像像素个数 

L=256;                                 %制定凸显灰度级为256 
count=count/N;                        
%各级灰度出现概率

for i=2:L 
    if count(i)~=0  ;       
st=i-1;         
break
    end
end
%以上循环语句实现寻找出现概率不为0的最小灰度值
for i=L:-1:1 
    
if count(i)~=0; 
    
nd=i-1;  

break   


end 

end 

%实现找出出现概率不为0的最大灰度值
f=count(st+1:nd+1); 

p=st;q=nd-st; 
%p和q分别是灰度的起始和结束值

u=0;

for i=1:q
    
    u=u+f(i)*(p+i-1); 
    
    ua(i)=u; 
    
end 
%计算图像的平均灰度值
for i=1:q 
    
    w(i)=sum(f(1:i)); 
    
end 

%计算出选择不同k的时候，A区域的概率 

d=(u*w-ua).^2./(w.*(1-w));   %求出不同k值时类间方差 

[y,tp]=max(d);                %求出最大方差对应的灰度值 

th=tp+p;
if th<=160     
th=tp+p; 
else     
th=160 
end                            %根据具体情况适当修正门限 

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;     
end 
end 
end 
%上面一段代码实现分割 
figure,imshow(Y1);title('灰度门限分割图像');