auto_thresh.m

用Matlab编写的自动多阈值分割程序

function auto_thresh(im,L,S,H)
% im=imread('rice.tif');
% a=imread('mytu.tif');
% im=rgb2gray(a);
count=imhist(im);
u=diff(count);%  Difference and approximate derivative.求差分曲线
% for i=1:255
%     u(i)=count(i+1)-count(i); 
% end
k=1;
for i=1:254         %找出差分曲线由负到正的过零点的位置（谷底位置）
    if u(i)<0
        if u(i+1)<0      
            i=i+1;
        elseif u(i+1)>0
               d(k)=i+0.5;                
               k=k+1;
        elseif u(i+1)==0
               for j=1:50
                   if (i+j<255&u(i+j)==0&u(i+j+1)>0)                      
                       d(k)=i+0.5*j+0.5;                
                       k=k+1;
                       break;
                   else 
                       break;
                   end
               end
           end
       end
  end
  D=diff(d);
%   for i=1:k-2
%       D(i)=d(i+1)-d(i);     %求各小峰的峰宽
%   end
  dd=uint8(D);
  count1=imhist(dd);            %统计峰宽度的直方图
  DD=find(count1);            
  kd=size(DD);
  n=round(0.7*(kd(1)-1)*2);        %求出平滑窗口的尺寸
  imm=wiener2(im,[n,n]);    %平滑原始直方图
  p=imhist(imm);
  
% 由平滑直方图确定区域类数（加入判断独立峰的三条准则）
%  for i=1:255
%     p1(i)=p(i)-p(i+1);     %求差分曲线
%     p2(i)=p(i+1)-p(i);
%  end
p2=diff(p);
p1=-p2;
k1=1;
for i=1:254               %找出差分曲线由负到正的过零点的位置（峰点位置）
    if p1(i)<0
        if p1(i+1)<0      
            i=i+1;
        elseif p1(i+1)>0
               dp(k1)=i+0.5;                
               k1=k1+1;
        elseif p1(i+1)==0
               j=2;
               while (j<20&i+j<256)
                   if p1(i+j)==0
                       j=j+1;
                   elseif p1(i+j)>0              
                       dp(k1)=i+0.5*(j-1)+0.5;                
                       k1=k1+1;
                       break;
                   else
                       break;
                   end
               end
           end
       end
  end 
  
k2=1;
for i=1:254                %找出差分曲线由负到正的过零点的位置（谷点位置）
    if p2(i)<0
        if p2(i+1)<0      
            i=i+1;
        elseif p2(i+1)>0
               dg(k2)=i+0.5;                
               k2=k2+1;
        elseif p2(i+1)==0
               j=2;
               while (j<20&i+j<256)
                   if p2(i+j)==0
                       j=j+1;
                   elseif p2(i+j)>0              
                       dg(k2)=i+0.5*(j-1)+0.5;                
                       k2=k2+1;
                       break;
                   else
                       break;
                   end
               end
           end
       end
  end 
  
  dp=round(dp);
  dg=round(dg); 
  for i=1:k2-2
      dgg(i)=dg(i+1)-dg(i);     %独立峰具有一定的灰度范围
      if dgg(i)>=L
          dgl(i)=1;
      else 
          dgl(i)=0;
      end        
  end
  for i=1:k2-2
      dgs(i)=sum(p(dg(i):dg(i+1)));
      if dgs(i)>=S
          dgS(i)=1;            %独立峰具有一定的面积
      else
          dgS(i)=0;
      end
  end
  if dp(1)<dg(1)
      for i=1:k2-2
          h1(i)=p(dg(i))/p(dp(i+1));
          h2(i)=p(dg(i+1))/p(dp(i+1));
          if(h1(i)<=1.0/H||h2(i)<=1.0/H)   %独立峰具有一定的谷峰比
              h1h(i)=1;
          else
              h1h(i)=0;
          end         
      end
      js=k2-2;
  else 
      for i=1:k1-2
          h1(i)=p(dg(i))/p(dp(i));
          h2(i)=p(dg(i+1))/p(dp(i));
          if(h1(i)<=1.0/H||h2(i)<=1.0/H)
              h1h(i)=1;
          else
              h1h(i)=0;
          end
      end
      js=k1-2;
  end
  for i=1:js
      if(dgl(i)==1&dgS(i)==1&h1h(i)==1)
         fzd(i)=1;
      else
         fzd(i)=0;
     end
 end
% fzd(find(fzd))
% for i=1:js
%     f1=fzd(ffz(i));
%     f2=fzd(ffz(i)+1);
% end