两个matlab实现最大熵法图像分割程序

%两个程序，亲测可用

clear all

a=imread('moon.tif');

figure,imshow(a)

count=imhist(a);

[m,n]=size(a);

N=m*n;

L=256;

count=count/N;%%每一个像素的分布概率

count

for i=1:L

if count(i)~=0

st=i-1;

break;

end

end

st

for i=L:-1:1

if count(i)~=0

nd=i-1;

break;

end

end

nd

f=count(st+1:nd+1); %f是每个灰度出现的概率

size(f)

E=[];

for Th=st:nd-1 %%%设定初始分割阈值为Th

av1=0;

av2=0;

Pth=sum(count(1:Th+1));

%%%第一类的平均相对熵为

for i=0:Th

av1=av1-count(i+1)/Pth*log(count(i+1)/Pth+0.00001);

end

%%%第二类的平均相对熵为

for i=Th+1:L-1

av2=av2-count(i+1)/(1-Pth)*log(count(i+1)/(1-Pth)+0.00001); end

E(Th-st+1)=av1+av2;

end

position=find(E==(max(E)));

th=st+position-1

for i=1:m

for j=1:n

if a(i,j)>th

a(i,j)=255;

else

a(i,j)=0;

end

end

end

figure,imshow(a);

%%%%%%%%%%%%%%%%%%%%%2-d 最大熵法(递推方法) %%%%%%%%%%% clear all;

clc;

tic

a=imread('trial2_2.tiff');

figure,imshow(a);

a0=double(a);

[m,n]=size(a);

h=1;

a1=zeros(m,n);

% 计算平均领域灰度的一维灰度直方图

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(i,j)=a0(p,q)+a1(i,j);

end

end

a2(i,j)=uint8(1/9*a1(i,j));

end

end

fxy=zeros(256,256);

% 计算二维直方图

for i=1:m

for j=1:n

c1=a0(i,j);

d=double(a2(i,j));

fxy(c1+1,d+1)=fxy(c1+1,d+1)+1;

end

end

Pxy=fxy/m/n;

% figure,

% mesh(Pxy);

% title('二维灰度直方图');

%计算Hl

Hl=0;

for i=1:256

for j=1:256

if Pxy(i,j)>0.00001

Hl=Hl-Pxy(i,j)*log(Pxy(i,j));

else

Hl=Hl;

end

end

end

%计算PA,HA

PA=zeros(256,256);

HA=zeros(256,256);

PB=zeros(256,256);

PA(1,1)=Pxy(1,1);

if PA(1,1)<1e-4

HA(1,1)=0;

else

HA(1,1)=-PA(1,1)*log(PA(1,1));

end

for i=2:256

PA(i,1)=PA(i-1,1)+Pxy(i,1);

if Pxy(i,1)>0.00001

HA(i,1)=HA(i-1,1)-Pxy(i,1)*log(Pxy(i,1));

else

HA(i,1)=HA(i-1,1);

end

end

for j=2:256

PA(1,j)=PA(1,j-1)+Pxy(1,j);

if Pxy(1,j)>0.00001

HA(1,j)=HA(1,j-1)-Pxy(1,j)*log(Pxy(1,j));

else

HA(1,j)=HA(1,j-1);

end

end

for i=2:256

for j=2:256

PA(i,j)=PA(i-1,j)+PA(i,j-1)-PA(i-1,j-1)+Pxy(i,j);

if Pxy(i,j)>0.00001

HA(i,j)=HA(i-1,j)+HA(i,j-1)-HA(i-1,j-1)-Pxy(i,j)*log(Pxy(i,j));

else

HA(i,j)=HA(i-1,j)+HA(i,j-1)-HA(i-1,j-1);

end

end

end

%计算最大熵

PB=1-PA;

h=zeros(256,256);

hmax=0;

for i=1:256

for j=1:256

if abs(PA(i,j))>0.00001&abs(PB(i,j))>0.00001

h(i,j)=log(PA(i,j)*PB(i,j))+HA(i,j)/PA(i,j)+(Hl-HA(i,j))/PB(i,j);

else

h(i,j)=0;

end

if h(i,j)>hmax

hmax=h(i,j);

s=i-1;

t=j-1;

end