%景曉軍文章中的二維直方圖
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