概念
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
endOstu算法（个人实现 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('灰度门限分割图像');