otsu算法---matlab实现，和一种改进算法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   灰度图阈值分割提取前后背景
%   otsu法和一种改进的方法
%   th为最终确定的阈值
%   2015.12.04  by小海添加注释
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
I=imread('C:\Users\Administrator\Desktop\a.bmp');
%a=rgb2gray(I);
a=I;
figure;
imshow(a);
figure;
[n,xout]=hist(a(:),0:255);  %n是落在该灰度值的个数，xout是中心点的刻度
bar(xout,n);
xlim([0 255]);      %把x轴的范围设置在0-255之间
count=imhist(a);    %统计的是图像上对应像素值的个数，是1-256的列向量            
[m,n]=size(a);     %m=768， n=1024  图像的列和行
count=count/(m*n); %归一化数据，概率百分比，每一个灰度级别的概率

%%%%%%%%%%%%%%%%%%%%%%%%%%
%方法11111111111111111111  
%   otsu法来求阈值分割图像
%%%%%%%%%%%%%%%%%%%%%%%%%%%
E=[];              %E是一个行向量，存储中间变量的值
for TH=1:1:256     %阈值从1开始计算   
    av0=0;
    av1=0;
    w0=sum(count(1:TH));    %统计  1~（阈值+1）处的总的概率和w0
    w1 = 1 - w0;
    %第一类  背景区域
    for i=1:TH
        av0= av0 + i*count(i);   % 求背景区域的平均灰度值
    end 
    u0 = av0/w0;                 %平均灰度值/背景的总概率
    %第二类
    for i=TH+1:256
        av1= av1 + i*count(i);   %这是一个计算前景均值的方法
    end 
    u1 = av1/w1;      %平均灰度值/背景的总概率
    E(TH) = w0*w1*(u1-u0)*(u1-u0);   %将两者的平均值相加
end 
postion=find(E==max(E));  %找到E最大处的值的位置，这个位置就可以确定为阈值  
th=postion; 


%%%%%%%%%%%%%%%%%%%%%%%
%方法22222222222222222222
%不是otsu法，是一种新的方法
%%%%%%%%%%%%%%%%%%%%%%

% L=256;
% for i=1:L             %找到开始的位置
%     if count(i)~=0    %频数不等于0的地方开始
%         st=i-1;       %从小到大，找到第一个非零点st
%         break;
%     end
% end 
% for j=L:-1:1 %找到结束的位置
%     if count(j)~=0        %从最后开始L=256
%     nd=j-1;               %从后到前找到第一个非零点nd
%     break;
%     end 
% end 
% %曲线拟合  多项式拟合,
% cline=[st:nd];    %cline从st到nd，中间差值为1，是一个行向量，包含了非零点的数目
% p=polyfit(cline,count(st:nd)',5);    %曲线拟合函数，cline作为一个横坐标值，count是一个包含频数的矩阵，对应着频数值，一阶直线拟合，二阶抛物线拟合，并非阶次越高越好，看拟合情况而定。
%                                      %为什么对count进行一个转置，因为统计出来的频数count是一个列向量，转置成行向量
% y=polyval(p,cline);                  %多项式在每一个cline处的值y
% plot([1:256],count,'*r',cline,y);    %画出拟合曲线图
% %[value,index]=max(y)
% IndMin=find(diff(sign(diff(y)))>0)+1;   %极小值点
% IndMax=find(diff(sign(diff(y)))<0)+1;   %极大值点
% f=count(st+1:nd+1);
% 
% %%%%%%%%%%%%%%%%%
% %%%%%%%  计算阈值
% %%%%%%%%%%%%%%%%%
% E=[];              %E是一个行向量，存储中间变量的值
% for TH=st:1:nd-1   %阈值从非空的点开始计算   st=32
%     av1=0;
%     av2=0;
%     pth=sum(count(1:TH+1));    %统计  1~（阈值+1）处的总的概率和w0
%     %第一类  背景区域
%     for i=0:TH
%         av1=av1-count(i+1)/pth*log(count(i+1)/pth+0.00001);   %这是一个计算背景均值的方法   求和 pi/w0*log（pi/w0+0.0001）   每一个点的概率除以第一类的总概率
%     end 
%     %第二类
%     for i=TH+1:L-1
%         av2=av2-count(i+1)/(1-pth)*log(count(i+1)/(1-pth)+0.00001);   %这是一个计算前景均值的方法
%     end 
%     E(TH-st+1)=av1+av2;   %将两者的平均值相加
% end 
% postion=find(E==max(E));  %找到E最大处的值的位置，这个位置就可以确定为阈值   
% th=st+postion-1; 


 




%%%%%%%%%%%%%下面是根据阈值来分割图像
for i=1:m
    for j=1:n
        if a(i,j)>th
            a(i,j)=255;
        else 
            a(i,j)=0;
        end 
    end 
end 
figure;
imshow(a)