matlab实现otsu算法

clear;
clc;
tic
rgb_image=imread（‘baitu3.jpg‘‘jpg‘）; 
a=rgb2gray（rgb_image）;
%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）*P