clc
close all
clear all
k=imread('A7.jpg');
k=im2double(k);
s=size(k);
figure(1),imshow(k,[])
title('input image')

%step1: zeropadding
f=zeros(s(1)+90,s(2)+90);
f(46:s(1)+45,46:s(2)+45)=k;
s=size(f);

%step2: teta angle
h=fspecial('sobel');
dfx=imfilter(f,h);
dfy=imfilter(f,h');
for i=1:s(1)
    for j=1:s(2)
        if dfx(i,j)==0,dfx(i,j)=0.0001;
        end
        D=dfy(i,j)/dfx(i,j);
        teta(i,j)=(180/pi)*atan(D);
    end
end
figure(2),imshow(teta,[])
title('teta angle')

%step3: Quantizing teta
for i=1:s(1)
    for j=1:s(2)
        if ( teta(i,j)<0 )  teta(i,j)=teta(i,j)+180;               end
        if ( teta(i,j)>=0 && teta(i,j)<22.5 )      teta1(i,j)=0;    end
        if ( teta(i,j)>=22.5 && teta(i,j)<67.5 )   teta1(i,j)=45;   end
        if ( teta(i,j)>=67.5 && teta(i,j)<112.5 )  teta1(i,j)=90;   end
        if ( teta(i,j)>=112.5 && teta(i,j)<157.5)  teta1(i,j)=135; end
        if ( teta(i,j)>=157.5 && teta(i,j) <=180 ) teta1(i,j)=0;    end
    end
end
figure(3),imshow(teta1,[])
title('quantized teta angle')

%step4: edge detector
bw=edge(f,'canny',0.3);
figure(4),imshow(bw,[])
title('edge detector with canny')

%step5: calculate diameter
m=1;d=5;
for i=1:s(1)
    for j=1:s(2)
        if bw(i,j)==1
            if teta1(i,j)==0
                k=20;
                while k<45
                    if (bw(i-k,j)==1 && teta(i,j)-d<=teta(i-k,j) && teta(i-k,j)<=teta(i,j)+d)
                        q(m)=k;m=m+1;break;end
                    k=k+1;
                end
                k=20;
                while k<45
                    if (bw(i+k,j)==1 && teta(i,j)-d<=teta(i+k,j) && teta(i+k,j)<=teta(i,j)+d)
                        q(m)=k;m=m+1;break;end
                    k=k+1;
                end
            end
            if teta1(i,j)==45
                k=20;
                while k<45
                    if (bw(i-k,j-k)==1 && teta(i,j)-d<=teta(i-k,j-k) && teta(i-k,j-k)<=teta(i,j)+d)
                        q(m)=k;m=m+1;break;end
                    k=k+1;
                end
                k=20;
                while k<45
                    if (bw(i+k,j+k)==1 && teta(i,j)-d<=teta(i+k,j+k) && teta(i+k,j+k)<=teta(i,j)+d)
                        q(m)=k;m=m+1;break;end
                    k=k+1;
                end
            end
            if teta1(i,j)==90
                k=20;
                while k<45
                    if (bw(i,j-k)==1 && teta(i,j)-d<=teta(i,j-k) && teta(i,j-k)<=teta(i,j)+d)
                        q(m)=k;m=m+1;break;end
                    k=k+1;
                end
                k=20;
                while k<45
                    if (bw(i,j+k)==1 && teta(i,j)-d<=teta(i,j+k) && teta(i,j+k)<=teta(i,j)+d)
                        q(m)=k;m=m+1;break;end
                    k=k+1;
                end
            end
            if teta1(i,j)==135
                k=20;
                while k<45
                    if (bw(i-k,j+k)==1 && teta(i,j)-d<=teta(i-k,j+k) && teta(i-k,j+k)<=teta(i,j)+d)
                        q(m)=k;m=m+1;break;end
                    k=k+1;
                end
                k=20;
                while k<45
                    if (bw(i+k,j-k)==1 && teta(i,j)-d<=teta(i+k,j-k) && teta(i+k,j-k)<=teta(i,j)+d)
                        q(m)=k;m=m+1;break;end
                    k=k+1;
                end
            end
        end
    end
end

%step6: regularity
figure(5),
for i=20:44
    p=find(q==i);
    s=size(p);
    w(i-19)=s(2);
    stem(i,w(i-19)),hold on
end
title('diameter of fiber')
xlabel('diameter/pixel')
ylabel('number of each diameter')
z=max(w);
y=find(z==w);
s=0;
for i=y-3:y+3
    s=sum(w(i))+s;
end
s1=sum(w);
percent=(s/s1)*100
mean=mean(q)
variance=var(q)