261220205Renyi.rar

代码 close all clear all; clc; tic a=imread('rice.png'); subplot(1,2,1); imshow(a); a0=double(a); [m,n]=size(a); h=1; a1=zeros(m,n); fxy=zeros(256,256); % 计算平均领域灰度的一维灰度直方图 mmimg=zeros(m,n); 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(i,j)的局部窗内容lwc,它为一个矩阵 lwc=a(i-uw:i+dw,j-lw:j+rw); %lwc:local windows content用于存储局部窗口内所有像素的灰度值 lwc=lwc(:); mmimg(i,j)=mean(lwc); clear lwc; end end % 计算二维直方图 for i=1:m for j=1:n c1=a0(i,j); d=round(mmimg(i,j)); fxy(c1+1,d+1)=fxy(c1+1,d+1)+1; end end Pxy=fxy/m/n; % figure, % mesh(Pxy); % title('二维灰度直方图'); P0=zeros(m,n); %原图像在阈值点（s，t）目标的先验概率 for s=0:m-1 for t=0:n-1 P0(s+1,t+1)=sum(sum(Pxy(1:s+1,1:t+1))); end end P1=ones(m,n)-P0; %原图像在阈值点（s，t）背景的先验概率 %二维直方图总的均值 MUT0=0; MUT1=0; for i=1:m for j=1:n MUT0=MUT0+i*Pxy(i,j); MUT1=MUT1+j*Pxy(i,j); end end %计算原图像在阈值为s时均值 mu00=zeros(m,n); MU00=zeros(m,n); %原图像在阈值为s时目标的均值 MU10=zeros(m,n); %原图像在阈值为s时背景的均值 u=0; e=0; for s=0:m-1 for t=0:n-1 if (s==0 & t==0) mu00(0+1,0+1)=0; MU00(0+1,0+1)=0; MU10(0+1,0+1)=MUT0/(P1(0+1,0+1)+eps); else if t==0 mu00(s+1,1)=mu00(s,1)+(s+1)*Pxy(s+1,1); u=mu00(s+1,1)/(P0(s+1,1)+eps); MU00(s+1,1)=u; MU10(s+1,1)=(MUT0-mu00(s+1,1))/(P1(s+1,1)+eps); else if s==0 mu00(1,t+1)=mu00(1,t)+0*Pxy(1,t+1); u=mu00(1,t+1)/(P0(1,t+1)+eps); MU00(1,t+1)=u; MU10(1,t+1)=(MUT0-mu00(1,t+1))/(P1(1,t+1)+eps); else mu00(s+1,t+1)=mu00(s,t+1)+mu00(s+1,t)-mu00(s,t)+(s+1)*Pxy(s+1,t+1); u=mu00(s+1,t+1)/(P0(s+1,t+1)+eps); MU00(s+1,t+1)=u; MU10(s+1,t+1)=(MUT0-mu00(s+1,t+1))/(P1(s+1,t+1)+eps); end end end end end %计算邻域平均图像在阈值为t时均值 mu01=zeros(m,n); MU01=zeros(m,n); %计算邻域平均图像在阈值为t时目标均值 MU11=zeros(m,n); %计算邻域平均图像在阈值为t时背景均值 for s=0:m-1 for t=0:n-1 if (s==0 & t==0) mu01(0+1,0+1)=0; MU01(0+1,0+1)=0; MU11(0+1,0+1)=MUT1/(P0(0+1,0+1)+eps); else if t==0 mu01(s+1,1)=mu01(s,1)+0*Pxy(s+1,1); h=mu01(s+1,1)/(P1(s+1,1)+eps); MU01(s+1,1)=h; MU11(s+1,1)=(MUT1-mu01(s+1,1))/(P0(s+1,1)+eps); else if s==0 mu01(1,t+1)=mu01(1,t)+(t+1)*Pxy(1,t+1); h=mu01(1,t+1)/(P1(1,t+1)+eps); MU01(1,t+1)=h; MU11(1,t+1)=(MUT1-mu01(1,t+1))/(P0(1,t+1)+eps); else mu01(s+1,t+1)=mu01(s,t+1)+mu01(s+1,t)-mu01(s,t)+(t+1)*Pxy(s+1,t+1); h=mu01(s+1,t+1)/(P1(s+1,t+1)+eps); MU01(s+1,t+1)=h; MU11(s+1,t+1)=(MUT1-mu01(s+1,t+1))/(P0(s+1,t+1)+eps); end end end end end %计算交叉熵熵 h=zeros(256,256); hmax=0; for i=1:256 for j=1:256 if abs(P0(i,j))>0.00001&abs(P1(i,j))>0.00001 h(i,j)=P0(i,j)*(MU00(i,j)*logMU00(i,j)+MU01(i,j)*logMU01(i,j))+P1(i,j)*(MU10(i,j)*logMU10(i,j)+MU11(i,j)*logMU11(i,j)); else h(i,j)=0; end if h(i,j)>hmax hmax=h(i,j); s=i-1; t=j-1; end end end z=ones(m,n); for i=1:m for j=1:n if a0(i,j)<=s&mmimg(i,j)<=t %if double(a(i,j))+double(a2(i,j))+a3(i,j)<=s+t+q z(i,j)=0; else z(i,j)=255; end end end hmax s t subplot(1,2,2); imshow(z); toc