OTSU算法二维matlab实现

%景晓军文章中的二维直方图
clear all;close all;
clc;
tic;
I=imread（‘C:\Users\Thomas\Desktop\left.bmp‘）;
[mnk]=size（I）;
if（k==3）
    a=rgb2gray（I）;
else
    a=I;
end
%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;
[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