gusucode.com > 红外图像增强及目标检测演示界面matlab源码程序 > code/WaveFisherImSg.m

    function WaveFisherImSg
%主要适用于斑点目标图像
  H=msgbox('主要适用于斑点目标图像');waitfor(H);
  I0=U_Open; 
  if (I0==0) H=errordlg('图像打开失败');waitfor(H);break; end
  label=inputdlg('若图像为暗目标图像，请输入0，否则输入1','图像类型标志',1,{'0'});
  label=str2num(char(label));
  if isrgb(I0)  I0=rgb2gray(I0);end
  [row col]=size(I0);
  if(label==1)%图像为天空背景弱小亮点目标图像，图像带边匡
     I=im2double(I0(30:row-35,45:col-44));
     I0=I;
  elseif(label==0)%图像为海空背景下的斑点目标图像，目标较背景暗
       I0=im2double(I0(10:row-20,10:col-10));
       m=max(I0(:));
       I=m-I0;
   else
      errordlg('输入图像类型错误');break;
   end
 tic
  [c1,s1]=wavedec2(I,4,'db2');
  sz=s1(1,:);%size of app. coef.(N=4)
  c1(1:sz(1)*sz(2))=0;%令低频部分为0
  f1=waverec2(c1,s1,'db2');%得到主要包含高频信息的图像
  f1=f1.*(f1>0);
  f1=medfilt2(f1);
   h=histog(f1);
   [sgf,Tf]=FisherSg(f1,h);
   [psgf,px,py]=ImPostProcess(sgf,label);
 [sgb,Tb]=BtwClsVarImSg(I);
 H=msgbox(['Fiser分割阈值',num2str(Tf)]);waitfor(H)
 H=msgbox(['内间方差分割阈值',num2str(Tb)]); waitfor(H)
figure
subplot(2,2,1);
imshow(I)
title('原图像')
d=10;
hold on;
   plot([py-d py+d],[px-d px-d],'w');
   plot([py-d py+d],[px+d px+d],'w');
   plot([py-d py-d],[px-d px+d],'w');
   plot([py+d py+d],[px-d px+d],'w');
hold off;
subplot(2,2,2);
imshow(sgf*0.9);
title('Fiser分割图像')
subplot(2,2,3);
imshow(psgf);
title('Fiser分割后处理图像')
subplot(2,2,4);
imshow(sgb);
title('内间方差分割图像')
msgbox(['所用时间',num2str(toc),'秒']);
%=============================================================
 function h=histog(f)
 %计算直方图
   m=max(f(:));
  f=round(f*255/m);
 [row col]=size(f);
 h=zeros(1,256);
 p=1/(row*col);
 for x=1:row
     for y=1:col
         i=f(x,y)+1;
         if(i>256)
             i=256;
         end
         h(i)=h(i)+p;
     end
 end
 %=============================================================
 function [sgf,T]=FisherSg(f,h)
%sgf为目标提取结果图像，亮点为目标，px,py为目标质心
%假定目标为亮点目标时的情况，暗点目标先运行:f=max(f(:))-f; 
m=max(f(:));
f=round(f*255/m);
D=zeros(256,256);
for i=1:255
  for j=(i+1):256
       m=0;
      for l=i:j
        m=m+h(l);
       end
      m=m/(j-i+1);
      for l=i:j
         %D(i,j)=D(i,j)+(h(l)-m)^2;
         D(i,j)=D(i,j)+abs(h(l)-m);
       end
      D(j,i)=D(i,j);
   end
end
code=zeros(256,256);
ep=code;
for n=3:256
      ep(n,2)=10000;code(n,2)=2;
      d=0;
       for j=2:n
           d=D(1,j-1)+D(j,n);
           if(d<ep(n,2))
               ep(n,2)=d;
               code(n,2)=j;
           end
       end
  end
for k=3:10
    for n=k+1:256                     
        ep(n,k)=100000;code(n,k)=2;
        d=0.0;
        for j=k:n
            d=ep(j-1,k-1)+D(j,n);
            if(d<ep(n,k))
               ep(n,k)=d;
               code(n,k)=j;
           end
       end
   end
end

class=2;a1=0.0;a2=.0;a=0.0;
 for i=3:8
     a1=atan(ep(256,i)-ep(256,i-1));
     a2=atan(ep(256,i+1)-ep(256,i));
     if(abs(a1-a2)>a)
         a=abs(a1-a2);class=i+1;
     end
 end
% class=class+1;
 T=code(256,class)-1;% 分割阈值
sgf=f>=T;%阈值分割——目标提取

%=========================================
function [f,px,py]=ImPostProcess(sgf,L)
if(L==1)
      for i=1:10 
        sgf=bwmorph(sgf,'dilate',1);
        sgf=bwmorph(sgf,'erode',1);
      end
     sgf=bwmorph(sgf,'erode',1);
     sgf=bwmorph(sgf,'dilate',1);
     sgf=im2double(sgf);
else
     for i=1:10 
         sgf=bwmorph(sgf,'erode',1);
         sgf=bwmorph(sgf,'dilate',1);
     end
     sgf=bwmorph(sgf,'erode',2);
     sgf=bwmorph(sgf,'dilate',1);
     sgf=im2double(sgf);
end
 [row col]=size(sgf);
 f=sgf;
 px=0;py=0;m=0;
 for i=1:row %目标质新计算
     for j=1:col
         if(sgf(i,j)~=0)
             px=px+i;
             py=py+j;
             m=m+1;
         end
     end
 end
 if(m~=0)
     px=px/m;
     py=py/m;
 end