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