mgabor.m

The matlab program is effective to be used for fringe processing. It contours strain distribution us

%upload on 26-05-09. To contact author, Dr Wang Jun, please drop you message at [e-mail] 

% Mgabor implements calculation of fringe frequency using one
% interferogram and is used to directly obtain strain from one single
% interferogram.
%
% This function returns a work space, mgabor.mat, in the current directory,
%	where imdata - matrix data of the pattern to be processed
%	(imwid, imheight) - the size of the imdata
%	ref_freq - frequency of the reference grating
%	(U1,V1) - frequency components of the Gabor channel filters
%	freq - fringe frequency components in the horizontal and vertical
%	divided by the frequency of the reference grating

%   a 5x5 filber band is generated and used to calculate the fringe frequnecy
%   in the horizontal and vertical directions.
%   Wang Jun, 4 April 2007

clear all
disp('Step 1: Get the file name and read the image into matrix');
imname=input('Enter the file name: ','s');
imdata=imread(imname);
[imheight,imwid]=size(imdata);
if isrgb(imdata)
    imdata=imdata(:,:,1);
elseif isbw(imdata)
    return;       
end
imdata=imdata(:,:,1);
[imheight,imwid]=size(imdata);
figure, imshow(imdata);title('Fringe pattern to be analysed');

disp('  Determine image domain for the given pattern by selceting one level from the histogram');
hnl=figure;imhist(imdata);
thres=ginput(1);
imdomain=~(imdata<=thres(1));
close(hnl);

ref_freq=2400; %twice the frequency of the specimen grating (lines/mm)
      
% save the data into work space mgabor.mat
save mgabor imdata;
save mgabor imwid imheight ref_freq -APPEND;

disp(' ');
disp('Step 2: Design a filter bank for the given pattern');
disp('  Press any key to contine');
pause;
sptr=fft2(double(imdata));
sptr=fftshift(sptr);
sptr_power = abs(sptr);
hnl=figure;imagesc(log(sptr_power));
colormap(gray);
hold on;

disp('  Choose the region representing the first-order spectrum.')
[x y]=ginput(2);
roi = sptr_power(round(min(y)):round(max(y)),round(min(x)):round(max(x)));
[v1 u1]=find(roi == max(max(roi)));
v1=v1+round(min(y))-1;u1=u1+round(min(x))-1;
plot(u1,v1,'r+');

disp('  Choose the region representing the zero-order spectrum.')
[x y]=ginput(2);
roi = sptr_power(round(min(y)):round(max(y)),round(min(x)):round(max(x)));
[v2 u2]=find(roi == max(max(roi)));
v2=v2+round(min(y))-1;u2=u2+round(min(x))-1;
plot(u2,v2,'r+');

u=u1-u2;v=v1-v2;
u=u/imwid;v=v/imheight; % convert cycle/image to cycle/pixel
F=sqrt(u^2+v^2);
A=atan2(v,u)*180/pi;%atan2, -pi<=A<=pi
close(hnl);

disp('  Input a number from 1 to 4 to choose frequency separations from large o small. ');
option = input('  your number (1,2,3,&4) = ');
alpha = sqrt(log(2)/2);%constant
if option ==1,
    std=3*alpha/(pi*F);
    disp1=1/(2*pi*std);
elseif option ==2,
    std=3*alpha/(pi*F);
    disp1=1/sqrt(8*pi^2*std^2);
elseif option ==3,
    std=1/F;
    disp1=1/(2*pi*std);
elseif option ==4,
    std=1/F;
    disp1=1/sqrt(8*pi^2*std^2);
else
    disply('please input a valid option.')
    return;
end

U1=zeros(5,5);V1=zeros(5,5);
for l=-2:1:2,
    U1(:,l+2+1)=F*cos(A*pi/180)+l*disp1;
end
for k=-2:1:2,
    V1(k+2+1,:)=F*sin(A*pi/180)+k*disp1;
end

figure,imagesc(log(sptr_power));
colormap(gray);
hold on;
plot(U1*imwid+fix(imwid/2),V1*imheight+fix(imheight/2),'r+');
title('Gabor channel filters');
truesize;

disp(' ');
disp('Step 3: Calculate response function for all channels.');
inc = 1;
g=gaborfilter(U1(2,3),V1(2,3));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1=abs(h);
fprintf(' Channel %d |',inc);

inc = inc+1;
g=gaborfilter(U1(1,1),V1(1,1));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_1=abs(h);
fprintf(' %d |',inc);

inc = inc+1;
g=gaborfilter(U1(2,1),V1(2,1));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_2=abs(h);
fprintf(' %d |',inc);

inc = inc+1;
g=gaborfilter(U1(3,1),V1(3,1));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_3=abs(h);
fprintf(' %d |',inc);

inc = inc+1;
g=gaborfilter(U1(4,1),V1(4,1));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_add_4=abs(h);
fprintf(' %d |',inc);

inc = inc+1;
g=gaborfilter(U1(5,1),V1(5,1));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_add_5=abs(h);
fprintf(' %d |',inc);

inc = inc+1;
g=gaborfilter(U1(1,2),V1(1,2));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_4=abs(h);
fprintf(' %d |',inc);

inc = inc+1;
g=gaborfilter(U1(2,2),V1(2,2));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_5=abs(h);
fprintf(' %d |',inc);

inc = inc+1;
g=gaborfilter(U1(3,2),V1(3,2));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_6=abs(h);
fprintf(' %d |',inc);

g=gaborfilter(U1(4,2),V1(4,2));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_add_7=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(5,2),V1(5,2));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_add_8=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(1,3),V1(1,3));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_7=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(3,3),V1(3,3));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_8=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(4,3),V1(4,3));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_add_9=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(5,3),V1(5,3));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_add_10=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(1,4),V1(1,4));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_9=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(2,4),V1(2,4));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_10=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(3,4),V1(3,4));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_11=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(4,4),V1(4,4));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_add_12=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(5,4),V1(5,4));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_add_13=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(1,5),V1(1,5));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_12=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(2,5),V1(2,5));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_13=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(3,5),V1(3,5));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_14=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(4,5),V1(4,5));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_add_15=abs(h);
inc = inc+1;
fprintf(' %d |',inc);

g=gaborfilter(U1(5,5),V1(5,5));
H=fft2(g,imheight,imwid);
H=fftshift(H);
H(1,1)=0;
H2=abs(H).*sptr;
h=ifft2(H2);
m1_add_16=abs(h);
inc = inc+1;
fprintf(' %d\n',inc);

save mgabor U1 V1 -APPEND;

disp(' ');
disp('Step 4: Determine fringe frequency according to the frequency of Gabor filters');
freq1=zeros(1,2);freq2=zeros(1,2);freq3=zeros(1,2);freq4=zeros(1,2);freq5=zeros(1,2);
Max1=0;Max2=0;Max3=0;Max4=0;Max5=0;
freq=zeros(imheight,imwid,2);
disp('  Computation will take some time; press any key to continue.');
pause;
disp('  Running...');

for l=1:imheight,
   for k=1:imwid,
      t1=[m1_1(l,k),m1_2(l,k),m1_3(l,k),m1_add_4(l,k),m1_add_5(l,k)];
      t2=[m1_4(l,k),m1_5(l,k),m1_6(l,k),m1_add_7(l,k),m1_add_8(l,k)];
      t3=[m1_7(l,k),m1(l,k),m1_8(l,k),m1_add_9(l,k),m1_add_10(l,k)];
      t4=[m1_9(l,k),m1_10(l,k),m1_11(l,k),m1_add_12(l,k),m1_add_13(l,k)];
      t5=[m1_12(l,k),m1_13(l,k),m1_14(l,k),m1_add_15(l,k),m1_add_16(l,k)];
      
      [Max1,I1]=max(t1);
      freq1=[U1(I1,1),V1(I1,1)];
      
      [Max2,I2]=max(t2);
      freq2=[U1(I2,2),V1(I2,2)];
      
      [ax3,I3]=max(t3);
      freq3=[U1(I3,3),V1(I3,3)];
      
      [Max4,I4]=max(t4);
      freq4=[U1(I4,4),V1(I4,4)];
      
      [Max5,I5]=max(t5);
      freq5=[U1(I5,5),V1(I5,5)];
    
      t=[freq1;freq2;freq3;freq4;freq5];
      [Max,I]=max([Max1,Max2,Max3,Max4,Max5]);
      freq(l,k,:)=t(I,:);
   end
   
   if l==fix(imheight*0.25),
      disp('  Finished 25%...');
   end
      
   if l==fix(imheight*0.5),
      disp('  Finished 50%...');
   end
         
   if l==fix(imheight*0.75),
      disp('  Finished 75%...');
   end
      
end
disp('  Finished 100%');

freq=freq/ref_freq;
figure,imagesc(medfilt2(double(imdomain).*freq(:,:,1),[3,3]));title('Derivative distribution in the horizontal');
colorbar;truesize;
figure,imagesc(medfilt2(double(imdomain).*freq(:,:,2),[3,3]));title('Derivative distribution in the vertical');
colorbar;truesize;
disp(' ');
disp('- End of the program -');

save mgabor freq -APPEND;