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

1、%两个程序，亲测可用clear alla=imread('moon.tif');figure,imshow(a) count=imhist(a);m,n=size(a);N=m*n;L=256;count=count/N;% 每一个像素的分布概率countfor i=1:Lif count(i)=0st=i-1;break;endendstfor i=L:-1:1if count(i)=0nd=i-1;break;endendndf=count(st+1:nd+1); %f 是每个灰度出现的概率 size(f)E=;for Th=st:nd-1 % 设定初始分割阈值为 Th a

2、v1=0;av2=0;Pth=sum(count(1:Th+1); %第一类的平均相对熵为 for i=0:Thav1=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);endE(Th-st+1)=av1+av2;endposition=find(E=(max(E);th=st+position-1for i=1:m for j=1:n if a(i,j)>th a(i,

3、j)=255; else a(i,j)=0; end end end figure,imshow(a);%2-d 最大熵法 (递推方法 ) clear all;clc;tica=imread('trial2_2.tiff'); figure,imshow(a);a0=double(a);m,n=size(a);h=1; a1=zeros(m,n);% 计算平均领域灰度的一维灰度直方图for i=1:mfor j=1:nfor k=-h:hfor w=-h:h;p=i+k;q=j+w; if (p<=0)|( p>m) p=i;endif (q<=0)|(q&g

4、t;n)q=j;end a1(i,j)=a0(p,q)+a1(i,j);end end a2(i,j)=uint8(1/9*a1(i,j);endendfxy=zeros(256,256);% 计算二维直方图for i=1:mfor j=1:nc1=a0(i,j);d=double(a2(i,j); fxy(c1+1,d+1)=fxy(c1+1,d+1)+1;endendPxy=fxy/m/n;% figure,% mesh(Pxy);% title(' 二维灰度直方图 ');%计算 HlHl=0;for i=1:256for j=1:256if Pxy(i,j)>0.0

5、0001Hl=Hl-Pxy(i,j)*log(Pxy(i,j);elseHl=Hl;endendend%计算 PA,HAPA=zeros(256,256);HA=zeros(256,256);PB=zeros(256,256);PA(1,1)=Pxy(1,1);if PA(1,1)<1e-4HA(1,1)=0;elseHA(1,1)=-PA(1,1)*log(PA(1,1);endfor i=2:256PA(i,1)=PA(i-1,1)+Pxy(i,1);if Pxy(i,1)>0.00001HA(i,1)=HA(i-1,1)-Pxy(i,1)*log(Pxy(i,1);elseH

6、A(i,1)=HA(i-1,1);end endfor j=2:256PA(1,j)=PA(1,j-1)+Pxy(1,j);if Pxy(1,j)>0.00001HA(1,j)=HA(1,j-1)-Pxy(1,j)*log(Pxy(1,j);elseHA(1,j)=HA(1,j-1);endendfor i=2:256for j=2:256PA(i,j)=PA(i-1,j)+PA(i,j-1)-PA(i-1,j-1)+Pxy(i,j);if Pxy(i,j)>0.00001HA(i,j)=HA(i-1,j)+HA(i,j-1)-HA(i-1,j-1)-Pxy(i,j)*log(Px

7、y(i,j); elseHA(i,j)=HA(i-1,j)+HA(i,j-1)-HA(i-1,j-1);endendend%计算最大熵PB=1-PA;h=zeros(256,256);hmax=0;for i=1:256for j=1:256if 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;endif h(i,j)>hmax hmax=h(i,j); s=i-1; t=j-1;end endz=ones(m,n);for i=1:mfor j=1:nif a0(i,j)<=s&a2(i,j)<