untitled91603djf.m

这个程序能清楚显示灰度图形的灰度直方图

%迭代阈值选取法
clear
I=uint8(imread('E:\工件图像\麻坑\D19.1Xd4.1X1.0\17.bmp','bmp'));
%I=rgb2gray(I);
I=medfilt2(I);

figure,imshow(I);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bw=imhist(I);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5%
%消除直方图噪声
aw=bw;
bw=bw';
x=1:256;
for j=1:5
  for i=1:255
     if abs(bw(i)-bw(i+1))>10
         bw(i+1)=fix(bw(i)+bw(i+1))/2;
     end
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
%选取合适分割阈值
max=0;
xmax=0;
for i=1:256
    if bw(i)>0
       if max<bw(i)
           max=bw(i);
           xmax=i;
       end
    end
end
min=max;
xmin=0;
for i=xmax:256
    if bw(i)>=0
        if min>bw(i)
            min=bw(i);
            xmin=i;
        end
    end
end
max1=0;
xmax1=0;
for i=xmax+10:xmin-4
    if bw(i)>0
        if bw(i)>(bw(i-1)+bw(i-2)+bw(i-3)+bw(i-4))/4
            if bw(i)<(bw(i+1)+bw(i+2)+bw(i+3)+bw(i+4))/4
               if max1<bw(i)
                   max1=bw(i);
                   xmax1=i;
               end
           end
        end
    end
end
max2=0;
xmax2=0;
for i=xmin-5:-1:xmax1+10
    if bw(i)>0
        if bw(i)>(bw(i-1)+bw(i-2)+bw(i-3)+bw(i-4))/4
            if bw(i)>(bw(i+1)+bw(i+2)+bw(i+3)+bw(i+4))/4
            if max2<bw(i)
               max2=bw(i);
               xmax2=i;
               break;
            end
           end
        end
    end
end
xmin1=0;
min1=max1;
for i=xmax:xmax1
    if bw(i)>0
        if bw(i)<(bw(i+1)+bw(i+2)+bw(i+3)+bw(i+4))/4
            if bw(i)<(bw(i-1)+bw(i-2)+bw(i-3)+bw(i-4))/4
                if min1>bw(i)
                    min1=bw(i);
                    xmin1=i;
                end
            end
        end
    end
end
min2=max1;
xmin2=0;
for i=xmax2:-1:xmax1+10
    if bw(i)>0
        if bw(i)<(bw(i+1)+bw(i+2)+bw(i+3)+bw(i+4))/4
            if bw(i)<(bw(i-1)+bw(i-2)+bw(i-3)+bw(i-4))/4
                if min2>bw(i)
                    min2=bw(i);
                    xmin2=i;
                end
            end
        end
    end
end
%%%%%%%%%%%%%%%%%%
%max
xmax
%min1
xmin1
%max1
xmax1
%min2
xmin2
%max2
xmax2
%min
xmin
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=1:1:256;
figure;
plot(x,bw,'-');
ylim([-200,4000]);
xlim([0,300]);
grid on
I1=im2bw(I,xmin1/255);
I2=im2bw(I,(xmin2-30)/255);
figure,subplot(1,2,1),imshow(I1);
subplot(1,2,2),imshow(I2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sum=0;
rslt=0;
aw=aw';
for i=1:256
    sum=sum+aw(i);
    rslt=rslt+i*aw(i);
end

u=fix(rslt/sum)
g=zeros(256,0);

for t=1:1:256
    sum1=0;
    rslt1=0;
    sum2=0;
    rslt2=0;
    for i=1:1:t
        sum1=sum1+aw(i);
        rslt1=rslt1+i*aw(i);
    end
    
    for i=t+1:1:256
        sum2=sum2+aw(i);
        rslt2=rslt2+i*aw(i);
    end
    
    u1=fix(rslt1/sum1);
    w1=sum1/sum;
    
    u2=fix(rslt2/sum2);
    w2=sum2/sum;
    
    u=u1*w1+u2*w2;
    
    g(t)=w1*(u1-u)^2+w2*(u2-u)^2;    
    
end

n=0;
temp=0;
for i=1:1:256
    if temp<g(i)
        temp=g(i);
        n=i;
    end
end
n
        
I3=im2bw(I,n/255);
figure,imshow(I3);