erweijiaochashang.m

clc
I=imread（‘beijing000.png‘）; 
subplot（121）; 
imshow（I）; 
I1=double（I）; 
[mn]=size（I）; 
h=1; 
I2=zeros（mn）;
fxy=zeros（256256）;
% 计算平均领域灰度的一维灰度直方图 
mmimg=zeros（mn）;
for i=1:m
    for j=1:n
        rw=0; %rw:right width指从当前像素点可以最多往右推多少个像素点到边界
        while （j+rw<=n && rw<=h）
              rw=rw+1;
        end
        rw=rw-1;
        
        lw=0; %lw:left width指从当前像素点可以最多往左推多少个像素点到边界
        while （j-lw>=1&&lw<=h）
            lw=lw+1;
        end
        lw=lw-1;
        
        uw=0; %uw:up width指从当前像素点可以最多往上推多少个像素点到边界
        while （i-uw>=1 && uw<=h）
            uw=uw+1;
        end
        uw=uw-1;
        
        dw=0; %dw:down width指从当前像素点可以最多往下推多少个像素点到边界
        while （i+dw<=m && dw<=h）
            dw=dw+1;
        end
        dw=dw-1;
        
        %求当前像素点img（ij）的局部窗内容lwc它为一个矩阵
        lwc=I（i-uw:i+dwj-lw:j+rw）; %lwc:local windows content用于存储局部窗口内所有像素的灰度值
        lwc=lwc（:）;
        mmimg（ij）=mean（lwc）;
        clear lwc;
    end
end


% 计算二维直方图 
for i=1:m 
    for j=1:n 
        c1=I1（ij）; 
        d=round（mmimg（ij））; 
        fxy（c1+1d+1）=fxy（c1+1d+1）+1; 
     end 
end 
Pxy=fxy/m/n;    
%  figure 
%   mesh（Pxy）; 
%   title（‘二维灰度直方图‘）; 
P0=zeros（mn）;  %原图像在阈值点（s，t）目标的先验概率
for s=1:m
    for t=1:n
        P0（st）=P0（st）+1;
     end
end
P1=ones（mn）-P0;  %原图像在阈值点（s，t）背景的先验概率

%二维直方图总的均值
MUT0=0;
MUT1=0;
for i=1:m
    for j=1:n
        MUT0=MUT0+i*Pxy;
        MUT1=MUT1+j*Pxy;
    end
end

%计算原图像在阈值为s时均值
mu00=zeros（mn）;
MU00=zeros（mn）;   %原图像在阈值为s时目标的均值
MU10=zeros（mn）;   %原图像在阈值为s时背景的均值
u=0;
e=0;
for s=0:255
    for t=0:255
        if （s==0 && t==0）
            mu00（0+10+1）=0;
            MU00（0+10+1）=0;
            MU10（0+10+1）=MUT0/（P1+eps）;
        else if t==0
                mu00（s+11）=mu00（s1）+（s+1）*Pxy（s+11）;
                u=mu00（s+11）/（P0（s+11）+eps）;
                MU00（s+11）=u;
                MU10（s+11）=（MUT0-mu00（s+11））/（P1（s+11）+eps）;
            else if s==0
                mu00（1t+1）=mu00（1t）+0*Pxy（1t+1）;
                u=mu00（1t+1）/（P0（1t+1）+eps）;
                MU00（1t+1）=u;
                MU10（1t+1）=（MUT0-mu00（1t+1））/（P1（1t+1）+eps）;
                else
                    mu00（s