untitled927.asv

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

%%处理裂纹
clear
I=uint8(imread('E:\工件图像\裂纹\D10.4Xd5.5X1(1)\20.bmp','bmp'));
%I=rgb2gray(I);
I=medfilt2(I);

%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))>1
         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

thres=(xmax1+xmax2)/2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=1:1:256;
figure;
plot(x,bw,'-');
ylim([-200,1200]);
xlim([0,300]);
grid on
I1=im2bw(I,xmin1/255);
I2=im2bw(I,xmin2/255);
figure,subplot(1,2,1),imshow(I1);
subplot(1,2,2),imshow(I2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=I;
width=768;
height=576;
%width=740;
%height=472;
%medfilt2(B);

for i=1:height
    for j=1:width
        if B(i,j)>xmin1
            B(i,j)=255;
        else
            B(i,j)=0;
        end
    end
end

leftx=zeros(1,height);
for i=1:height
    for j=1:width
        if B(i,j)>180
            leftx(i)=j;
        end
    end
end 

rightx=zeros(1,height);
for i=1:height
    for j=1:width
        if B(i,width+1-j)>180
            rightx(i)=width+1-j;
        end
    end
end 
xtime=0;
x=zeros(1,height);
for i=1:height
    if(leftx(i)-rightx(i))>0
        x(i)=(rightx(i)+leftx(i))/2;
        xtime=xtime+1;   
    end
end

topy=zeros(1,width);
for i=1:width
    for j=1:height
        if B(j,i)>180
            topy(i)=j;
        end
    end
end 

bottomy=zeros(1,width);
for i=1:width
    for j=1:height
        if B(height+1-j,i)>180
            bottomy(i)=height+1-j;
        end
    end
end

y=zeros(1,width);
for i=1:width
    if(topy(i)-bottomy(i))>0
        y(i)=(topy(i)+bottomy(i))/2;
    end
end

%统计x的真值
a=0;
sum=0;
for i=1:height
    if x(i)>0 
        sum=sum+x(i);
    end
end
Vx=0;
x1=sum/xtime;
for i=1:height
    if x(i)>0
        Vx=Vx+(x(i)-x1)*(x(i)-x1);
    end
end

a=sqrt(Vx/(xtime-1));

for i=1:height
    if abs(x(i)-x1)>a
        x(i)=0;
    end
end

xtime1=0;
sum1=0;
for i=1:height
    if x(i)>0
        sum1=sum1+x(i);  
        xtime1=xtime1+1;
    end
end

Vx1=0;
a1=0;
for i=1:height
    if x(i)>0
        Vx1=Vx1+(x(i)-sum1/xtime1)*(x(i)-sum1/xtime1);
    end
end
a1=sqrt(Vx1/(xtime1-1));
xtime2=0;
sum2=0;
for i=1:height
    if abs(x(i)-sum1/xtime1)>a1
        x(i)=0;
    end
    if x(i)>0
        sum2=sum2+x(i);
        xtime2=xtime2+1;
    end
end
Vx2=0;
a2=0;
for i=1:height
    if x(i)>0
        Vx2=Vx2+(x(i)-sum2/xtime2)*(x(i)-sum2/xtime2);
    end
end
a2=sqrt(Vx2/(xtime2-1));

xtime3=0;
sum3=0;
for i=1:height
    if abs(x(i)-sum2/xtime2)>a2
        x(i)=0;
    end
    if x(i)>0
       sum3=sum3+x(i);
       xtime3=xtime3+1;
   end
end

%统计y的真值
ysum=0;
ytime=0;

for i=1:width
    if y(i)>0
        ysum=ysum+y(i);
        ytime=ytime+1;
    end
end

Vy=0;
b=0;
for i=1:width
    if y(i)>0
        Vy=Vy+(y(i)-ysum/ytime)*(y(i)-ysum/ytime);
    end
end

b=sqrt(Vy/(ytime-1));

ysum1=0;
ytime1=0;
for i=1:width
    if abs(y(i)-ysum/ytime)>b
        y(i)=0;
    end
    if y(i)>0
        ysum1=ysum1+y(i);
        ytime1=ytime1+1;
    end
end

Vy1=0;
b1=0;
for i=1:width
    if y(i)>0
       Vy1=Vy1+(y(i)-ysum1/ytime1)*(y(i)-ysum1/ytime1);
   end
end

b1=sqrt(Vy1/(ytime1-1));
ysum2=0;
ytime2=0;
for i=1:width
    if abs(y(i)-ysum1/ytime1)>b1
        y(i)=0;
    end
    if y(i)>0
        ysum2=ysum2+y(i);
        ytime2=ytime2+1;
    end
end

Vy2=0;
b2=0;
for i=1:width
    if y(i)>0
        Vy2=Vy2+(y(i)-ysum2/ytime2)*(y(i)-ysum2/ytime2);
    end
end
b2=sqrt(Vy2/(ytime2-1));
ysum3=0;
ytime3=0;
for i=1:width
    if abs(y(i)-ysum2/ytime2)>b2
        y(i)=0;
    end
    if y(i)>0
        ysum3=ysum3+y(i);
        ytime3=ytime3+1;
    end
end

%计算圆心坐标
xcenter=fix(sum3/xtime3);
ycenter=fix(ysum3/ytime3);
%计算x方向的半径
xrsum=0;
xrtime=0;
xr=0;
for i=1:height
   if x(i)>0
     xrsum=xrsum+(sqrt((leftx(i)-xcenter)*(leftx(i)-xcenter)+(i-ycenter)*(i-ycenter))+sqrt((rightx(i)-xcenter)*(rightx(i)-xcenter)+(i-ycenter)*(i-ycenter)))/2;
     xrtime=xrtime+1;
   end
end
xr=fix(xrsum/xrtime);

%计算y方向的半径
yrsum=0;
yrtime=0;
yr=0;
for i=1:width
   if y(i)>0
     yrsum=yrsum+(sqrt((topy(i)-ycenter)*(topy(i)-ycenter)+(i-xcenter)*(i-xcenter))+sqrt((bottomy(i)-ycenter)*(bottomy(i)-ycenter)+(i-xcenter)*(i-xcenter)))/2;
     yrtime=yrtime+1;
     
   end
end
yr=fix(yrsum/yrtime);
%计算半径r
r=fix((yr+xr)/2);

xcenter
ycenter
r
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%截取有用数据
%C=zeros(r,r);
for i=(ycenter-r-1):(ycenter+r+1)
    for j=(xcenter-r-1):(xcenter+r+1)
        C(i+r+2-ycenter,j+r+2-xcenter)=I2(i,j);
        D(i+r+2-ycenter,j+r+2-xcenter)=I1(i,j);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%判断是圆环还是圆片
centerX=zeros(1,2*r+20);
centerY=zeros(1,2*r+20);
centerS=zeros(1,2*r+20);
for i=(ycenter-r-10):(ycenter+r+10)
    for j=(xcenter-r-10):(xcenter+r+10)
        if i==ycenter
            centerX(j-xcenter+r+10+1)=B(i,j);
        elseif j==xcenter
            centerY(i-ycenter+r+10+1)=B(i,j);
        end
    end
end
for i=1:2*r+20+1
    centerS(i)=(centerX(i)+centerY(i))/2;
end
%x=1:2*r+20+1;
%figure;
%plot(x,centerS);

center1=0;
center2=0;
for i=r+11:2*r+20+1
    if centerS(i)>0
         center1=i;
         break;
     end
 end
 for i=r+11:-1:1
     if centerS(i)>0
         center2=i;
         break;
     end
 end
 centerR=fix((center1-center2)/2);
 %%如果centerR>0，则内部为空，即圆环；否则为圆片
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %计算边缘光圈的最大的宽度 edgewidthMax
 xxedge=zeros(1,width);
 edgeN=zeros(1,width);
 edgeW=zeros(1,width);
 yy=zeros(1,width);
 kk=0;
 b=0;
 a=0;
 for k=0:0.1:6.28
     for i=1:1:(r+1)
         x=fix(xcenter+i*cos(k));
         y=fix(ycenter-i*sin(k));
         xxedge(i)=I(y,x);
     end
     
     for i=1:1:(r+1)
         if xxedge(i)<(xmax1+(xmin2-xmax1)/2)
             xxedge(i)=0;
         end
     end
     
     for i=1:1:(r+1)
         if xxedge(i)>0
            a=a+1;
        else
            if a>0
                b=b+1;
                yy(b)=a;
                a=0;
            end
         end
     end
     
     kk=fix(10*k)+1;
     edgeN(kk)=yy(1);
     if b>0
          edgeW(kk)=yy(b);
      end
 end
 
 Nedge=0;
 Nsum=0;
 Wedge=0;
 Wsum=0;
 for i=1:1:width
     if edgeN(i)>2
         Nsum=Nsum+1;
         Nedge=Nedge+edgeN(i);
     end
     if edgeW(i)>2
         Wsum=Wsum+1;
         Wedge=Wedge+edgeW(i);
     end
 end
 
 if Wsum>0 
     Wedge=fix(Wedge/Wsum)
 else
     Wedge=2
 end
 if Nsum>0
     Nedge=fix(Nedge/Nsum)
 else
     Nedge=2
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       
 if centerR>0 %圆环
     for i=1:2*(r+1)
         for j=1:2*(r+1)
            if (j-r-1)*(j-r-1)+(i-r-1)*(i-r-1)>(r-Wedge)*(r-Wedge)
                 C(i,j)=0;
            end
         end
     end
     %imagesc(C);
     for i=1:2*(r+1)
         for j=1:2*(r+1)
             if (j-r-1)*(j-r-1)+(i-r-1)*(i-r-1)<(centerR+Nedge)*(centerR+Nedge)
                  C(i,j)=0;  
             end
         end
     end
 elseif centerR==0 %圆片
     for i=1:2*(r+1)
         for j=1:2*(r+1)
            if (j-r-1)*(j-r-1)+(i-r-1)*(i-r-1)>(r-Wedge)*(r-Wedge)
                 C(i,j)=0;
            end
         end
     end
 end

E=medfilt2(medfilt2(C));
figure;
colormap(gray);
imagesc(E);

F=bwarea(E);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%计算缺陷区域面积
x_sign=0;
m_temp=0;
x_temp=0;
y_temp=0;
flag=0;

[L,num]=bwlabel(E,8);
result=zeros(1,num);
for i=1:2*(r+1)
    for j=1:2*(r+1)
        if num>0
            for k=1:num
                if L(i,j)==k
                    result(k)=result(k)+1;
                end
            end
        end
    end
end

%figure;
%colormap(gray);
%imagesc(L);
num
result
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%