tidu.m

matlab模拟产生中心应力集中光学条纹图

function y=TiDu(m,n,Isub);
% Parameters of the Gaussian filter:
%n1=5;sigma1=1.414;n2=10;sigma2=1.414;theta1=0;

%Isub=imread('circular.tif');

Isub=double(Isub);
[M N]=size(Isub);

H1=[-1 0 1
    -2 0 2
    -1 0 1];
H2=[-1 -2 -1
    0   0  0
    1   2  1];
fy1=conv2(Isub,H1,'same');
fy2=conv2(Isub,H2,'same');
%Idiv=-atan(fy2./fy1);
%Idiff=255.0.*(Idiv-min(min(Idiv)))/(max(max(Idiv))-min(min(Idiv)));
%figure, imshow(uint8(Idiff));
%figure, imshow(uint8(fy1));
%figure, imshow(uint8(fy2));
%  fyI=sqrt(fy1.*fy1+fy2.*fy2);   % intensity;
%  fyd=atan(fy2/fy1);       % direction;

%calculate the direction of every pixel
Imult=2.0*fy1.*fy2;      % sum of 2.*fy1*fy2
Iadd=fy2.*fy2-fy1.*fy1;        % sum of fy1*fy1+fy2*fy2
[M N]=size(fy1);
for i=1:M
    for j=1:N
              Im=0.0;      % sum of 2.*fy1*fy2
              Ia=0.0;          % sum of fy1*fy1+fy2*fy2
              
        for k=-m:m
            for l=-n:n
                if ((i+k)<1 || (j+l)<1 || (i+k)>M || (j+l)>N)
                    Im=Im;
                    Ia=Ia;
                else
                    Im=Im+Imult(i+k,j+l);
                    Ia=Ia+Iadd(i+k,j+l);
                end
            end
        end
        
        Imave(i,j)=Im;
        Iaave(i,j)=Ia;
    end
end


for i=1:M
    for j=1:N
        if (abs(Iaave(i,j))<1e-4)
            if(Iaave(i,j)*Imave(i,j)>=0.0)
                Wph(i,j)=pi/4;
            else 
                Wph(i,j)=-pi/4;
            end
        else
               Wph(i,j)=0.5*abs(atan(Imave(i,j)/Iaave(i,j)));

        end
    end
end

for i=1:M
    for j=1:N
        if (Imave(i,j)<0 ) 
            if (Iaave(i,j)<0)
                Wph(i,j)=pi/2+Wph(i,j);
            else Wph(i,j)=pi-Wph(i,j);
            end
        else
            if(Iaave(i,j)<0)
                Wph(i,j)=pi/2-Wph(i,j);
            end
        end
    end
end

y=Wph;