灰度图像的二维Otsu自动阈值分割法 matlab源程序

算法速度很快哦，只要0.4秒，风格效果稳定，用了cameraman测试的

%景晓军文章中的二维直方图
clear;
clc;
tic
a=imread（‘d:\image\xx.jpg‘‘jpg‘）;       %Read a image 
%a=imread（‘d:\image\clock.tif‘‘tif‘）;
%a=imread（‘d:\image\zhuangjia.tiff‘‘tiff‘）;
%a=a（::1）;
% % a=M;
% %load（‘source_N0.02.mat‘）;
% %a=X;
%a=imread（‘syn1-g2.gif‘）; 
%load（‘装甲车N_0.01.mat‘）;
%a=b;
%a=noise_h;
subplot（121）;
imshow（a）;  
[mn]=size（a）;
 
%b=imnoise（a‘salt & pepper‘0.003）;
%b=imnoise（b‘gaussian‘00.0015）;
%b = IMNOISE（a‘speckle‘0.09）;
%b=a;
a0=double（a）;
h=1;                          
a1=zeros（mn）;
% 计算平均领域灰度的一维灰度直方图
for i=1:m
    for j=1:n
        for k=-h:h
            for w=-h:h;
                p=i+k;
                q=j+w;
                if （p<=0）|（ p>m）
                    p=i;
                end
                if （q<=0）|（q>n）
                    q=j;
                end
                 a1（ij）=a0（pq）+a1（ij）;
             end
        end
        a2（ij）=uint8（1/9*a1（ij））;
    end
end
fxy=zeros（256256）;
% 计算二维直方图
for i=1:m
    for j=1:n
        c=a0（ij）;
        d=double（a2（ij））;
        fxy（c+1d+1）=fxy（c+1d+1）+1;
     end
  end
%  figure
%  mesh（fxy）;
%  title（‘二维灰度直方图‘）;
Pxy=fxy/m/n;
P0=zeros（256256）;
Ui=zeros（256256）;
Uj=zeros（256256）;
P0（11）=Pxy（11）;
for i=2:256
    P0（1i）=P0（1i-1）+Pxy（1i）;
end
for i=2:256
    P0（i1）=P0（i-11）+Pxy（i1）;
end
for i=2:256
    for j=2:256
        P0（ij）=P0（ij-1）+P0（i-1j）-P0（i-1j-1）+Pxy（ij）;
    end
end
P1=ones（256256）-P0;
Ui（11）=0;
for i=2:256
    Ui（1i）=Ui（1i-1）+（1-1）*Pxy（1i）;
end
for i=2:256
    Ui（i1）=Ui（i-11）+（i-1）*Pxy（i1）;
end
for i=2:256
    for j=2:256
        Ui（ij）=Ui（ij-1）+Ui（i-1j）-Ui（i-1j-1）+（i-1）*Pxy（ij）;
    end
end
Uj（11）=0;
for i=2:256
    Uj（1i）=Uj（1i-1）+（i-1）*Pxy（1i）;
end
for i=2:256
    Uj（i1）=Uj（i-11）+（1-1）*Pxy（i1）;
end
for i=2:256
    for j=2:256
        Uj（ij）=Uj（ij-1）+Uj（i-1j）-Uj（i-1j-1）+（j-1）*Pxy（ij）;
    end
end
uti=0;
utj=0;
for i=1:256
    for j=1:256
        uti=uti+（i-1）*Pxy（ij）; 
        utj=utj+（j-1）*Pxy（ij）;
    end
end
% for i=1:256
%     for j=1:256
%       Ui1（ij）=Ui（256256）-Ui（ij）;
%       Uj1（ij）=Uj（256256）-Uj（ij）;
%     end
% end
%计算类间方差
hmax=0;
for i=1:256
    for j=1:256
       if P0（ij）~=0&P1（ij）~=0
        h（ij）=（（uti*P0（ij）-Ui（ij））^2+（utj*P0（ij）-Uj（ij））^2）/（P1（ij）*P0（ij））;
       else
        h（ij）=0;
       end
    end
end
hmax=max（h（:））;
for i=1:256
    for j=1:256
        if h（ij）==hmax
             s=i-1;
             t=j-1;
             continue;
         end
     end
 end
hmax
s
t
z=ones（mn）; 
%c=zeros（mn）;  
for i=1:m
    for j=1:n
      if a（ij）<=s&a2（ij）<=t 
      %if double（b（ij））+double（a2（ij））<=s+t  
       z（ij）=0;      
      end
  end
end
%imwrite（z‘D:\MATLAB7\work\feng\文章涉及到的\曲线otsu\曲线otsu实验结果\baboon\wuprocess.jpg‘）; 
%imwrite（z‘D:\MATLAB7\work\feng\文章涉及到的\曲线otsu\曲线otsu实验结果\baboon\you_process.jpg‘）;
%imwrite（uint8（b）

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3208  2007-05-22 16:43  otsu_2d.m

----------- ---------  ---------- -----  ----

                 3426                    2