Otsu传统算法和快速算法的matlab实现

说明：每个%%都表示一个模块，每个模块可以单独执行，也可以将%%去掉整体执行。

将读入图像修改成自己机子上的，Otsu传统算法执行顺序如下：cell1-cell2-cell3, Otsu快速迭代算法执行顺序如下：cell1-cell2-cell4。

%%%%%%%%%%%%%%%%%%%%%%%%%%
%% cell1 读入图像1
I=imread('e:/photo/m/rice.png');
%% cell2 显示原图及其直方图
subplot(2,2,1);imshow(I);xlabel('a,原图');
subplot(2,2,2);imhist(I);
%% cell3 Otsu法传统算法的参数说明及步骤
%参数说明：
%M,N---图像大小
%TT,T-----阈值,使得方差为最大的阈值
%Ib,Ibav,Ibavm----背景图像,背景图像的灰度平均值,合适阈值时对应的背景图像的灰度平均值
%If,Ifav,Ifavm----前景图像,前景图像的灰度平均值,合适阈值时对应的前景图像的灰度平均值
%Iav----图像的平均灰度值（数学期望）
%Ivar,Ivarm----方差,值最大的方差
%%%%%%%%%%
%步骤：
%1.阈值T从图像的最小灰度值开始遍历，直到最大灰度值。

%2.划分前景和后景：所有I>=T的为背景Ib,得出背景所占比例pb与背景的平均灰度值Ibav;
%2.划分前景和后景：所有I<T的为前景If,得出前景所占比例pf与前景的平均灰度值Ifav;
%3.求出图像的平均灰度Iav和方差Ivar
%4.遍历结束后，Tm即为所求，然后二值化分割。

%%%%%%%%%%%%%%%%%%%%%% [M,N]=size(I);
Ivarm=0;
T=0;
Iav=mean(mean(I));
tic;
for TT=min(min(I)):max(max(I))
%划分前景背景
Ib=I(I>=TT);
pb=length(Ib)/(M*N);
Ibav=mean(Ib);
If=I(I<TT);
pf=length(If)/(M*N);
Ifav=mean(If);
%求类间方差，并找出最大方差
Ivar=pb*pf*(Ibav-Ifav)^2
if Ivar>Ivarm
Ivarm=Ivar;
T=TT;
end
end
toc;
% 输出
disp(strcat('大律法计算出的阈值: ',num2str(T)));
Iotsu=im2bw(I,double(T)/255);
subplot(2,2,2);imshow(Iotsu,[]);xlabel('b,大律法'); %%%%%%%%%%%%%%%%%%%%%%%%
%% cel4 OTSU传统算法的快速实现----迭代
%参数说明：
%M,N---图像大小
%T-----阈值(即使得方差为最大的阈值)
%gmin,gmax---图像灰度级最小值和最大值
%p---图像灰度级概率分布，一维数组长度gmin*gmax
%Ibp,Ibav,Ibavm----背景图像灰度级概率,背景图像的灰度平均值,合适阈值时对应的背景图像的灰度平均值
%Ifp,Ifav,Ifavm----前景图像灰度级概率,前景图像的灰度平均值,合适阈值时对应的前景图像的灰度平均值
%Iav----图像的平均灰度值（数学期望）
%Ivar,Ivarm----方差,值最大的方差
%%%%%%%%%%
%步骤：
%1.先求解gmin,gmax,p,Iav,并令Ibp=0,Ibav=0,Ivarm=0;
%2.i从gmin+1到gmax+1迭代搜索,(下面为方便说明，使用了序列)Ibp(i)=Ibp(i-1)+p(i),
Ibav(i)=(Ibp(i-1)*Ibav(i-1)+i*p(i))/Ibp(i),Ivar(i)=(Ibp(i)*(Iav-Ibav(i))^2)/(1-Ibp(i));
%2.判断Iavm是否大于Iav,若是,则i=i+1,重复第2步;若否,则Iavm=Iav,T=i,Ibavm=Ibav,
Ifavm=Ifav,i=i+1,重复第2步;
%3.迭代结束，T即为所求，二值化分割。

%%%%%%%%%%%%%%%%%%%%%%
[M,N]=size(I);
gmin=min(min(I));gmax=max(max(I));
[p,x]=imhist(I);
p=p/(M*N);
Iav=mean(mean(I));
Ibp=0;Ibav=0;Ivarm=0;
tic;
for i=gmin+1:gmax+1
x=Ibp*double(Ibav);
Ibp=Ibp+p(i);
xx=double(i)*p(i);
Ibav=(x+xx)/Ibp;
if Ibp<=0.001 || (1-Ibp)<=0.001%避免因分母过小出现错误
Ivar=0;
else
Ivar=(Ibp*(Iav-Ibav)^2)/(1-Ibp);%即 Ivar=Ibp*Ifp*(Ibav-Ifav)^2; end
if Ivar>Ivarm
Ivarm=Ivar;
T=i;
end
end
toc;
% 输出
disp(strcat('大律法快速算法计算出的阈值: ',num2str(T))); Iotsu=im2bw(I,double(T)/255);
subplot(2,2,4);imshow(Iotsu,[]);xlabel('c,大律法快速算法');
下面是运行结果：
>Elapsed time is 0.445446 seconds.
大律法计算出的阈值:130
>Elapsed time is 0.000050 seconds.
大律法快速算法计算出的阈值:129。