image thumbnail

Noncontact Speckle-Based Velocity Meter

by

 

Finding a velocity of substance by using radon transform.

Johann
function  Johann
clear all
global I1;
global I2;
figure('MenuBar','none','Name','Johann','NumberTitle','off','Position',[100, 100, 1000, 600]);
I2=imread('red-after.jpg');

I1=imread('red-before.jpg');


uicontrol('Style','PushButton','String','Basic','Position',[20,530,110,20],'CallBack',@BasicPressed);
uicontrol('Style','PushButton','String','Meshing','Position',[20,500,110,20],'CallBack',@MeshingPressed);
uicontrol('Style','PushButton','String','FFT','Position',[20,470,110,20],'CallBack',@FFTPressed);
uicontrol('Style','PushButton','String','Convolution','Position',[20,440,110,20],'CallBack',@ConvPressed);
uicontrol('Style','PushButton','String','Radon Transform','Position',[20,410,110,20],'CallBack',@RTPressed);
uicontrol('Style','PushButton','String','Correlation','Position',[20,380,110,20],'CallBack',@CorrelationPressed);
uicontrol('Style','PushButton','String','Velocity','Position',[20,350,110,20],'CallBack',@VelocityPressed);
uicontrol('Style','PushButton','String','Exit','Position',[20,320,110,20],'CallBack',@ExitPressed);
uicontrol('Style','PushButton','String','About','Position',[20,290,110,20],'CallBack',@AboutPressed);
    function BasicPressed(h,eventdata)
    figure (1)
    
    subplot(231)
     imshow(I1);
     xlabel('Original')
    subplot(232)
     imshow(I2);
     xlabel('After rotation')
    subplot(233)
     imshow(I1+I2)
     xlabel('Add')
    subplot(234)
     imshow(I1-I2)
     Xlabel('Subtract')
    subplot(235)
     imshow((I1+I2)/2)
     xlabel('Mean')
    subplot(236)
     I=rgb2gray(I1);
     imshow(I)
     xlabel('Grayscale Original')
    end

    function MeshingPressed(h,eventdata)
    figure(1)
    subplot(1,1,1)
        subplot (231)
         mesh(double(I1(:,:,1)))
         xlabel('Orinial Mesh Red')
        subplot (232)
         mesh(double(I1(:,:,2)))
         xlabel('Orinial Mesh Green')
        subplot (233)
         mesh(double(I1(:,:,3)))
         xlabel('Orinial Mesh Blue')
        subplot (234)
         mesh(double(I2(:,:,1)))
         xlabel('Rotated Mesh Red')
        subplot (235)
         mesh(double(I2(:,:,2)))
         xlabel('Rotated Mesh Green')
        subplot (236)
         mesh(double(I2(:,:,3)))
         xlabel('Rotated Mesh Blue')
    end
    function FFTPressed(h,eventdata)
      figure(1)
    subplot(1,1,1)
        fft_1=fft2(double(rgb2gray(I1)));
      fft_2=fft2(double(rgb2gray(I2)));
      subplot(221)
      imshow(abs(fftshift(fft_1)),[24 10000]), colormap gray
      xlabel('Original Image FFT2 Amplitude')
      subplot(222)
      imshow(angle(fftshift(fft_1)),[-pi pi]), colormap gray
        xlabel('Original Image FFT2 Phase')
       subplot(223)
      imshow(abs(fftshift(fft_2)),[24 10000]), colormap gray
      xlabel('Rotated Image FFT2 Amplitude')
      subplot(224)
      imshow(angle(fftshift(fft_2)),[-pi pi]), colormap gray
        xlabel('Rotated Image FFT2 Phase')
    end
    function ConvPressed(h,eventdata)
        figure(1)
    subplot(1,1,1)
        I1d=double(I1);
I2d=double(I2);
[x,y,z]=size(I1d);


A1=I1d(:,:,1);

B1=I2d(:,:,1);
global I1d;
global I2d;


f1=real(fft2(A1));
f2=real(fft2(B1));
C=conv2(f1,f2);
subplot(211)
mesh(C(1:50,1:50))
xlabel('Convolusion in Freq domain')
        C2=conv2(A1,B1);
        subplot(212)
        mesh(C2)
    xlabel('Convolusion 2d ')
   
figure(3) 
subplot(211)
mesh(C(1:50,1:50))
xlabel('Convolusion in Freq domain')
subplot(212)
        mesh(C2)
    xlabel('Conv 2d ')
    end

    function RTPressed(h,eventdata)
     figure(1)
     global f;
    subplot(1,1,1)
%        d=1;
        n=radon(I2,20);
        subplot(211)
        for q=1:360
         m=radon(I1,q);
         plot(m,'r')
         grid
         l=strcat('Radon Transform at: ',num2str(q),' degree');
         xlabel (l)
         pause(0.2)
                
        end
               f=20;
        N=radon(I2,f);
        M=radon(I1,[1:360]);
        MR=0;
        for i=1:360
            MR=[MSE1(M(:,i),N,573),MR];
        end
        subplot(212)
        plot(MR);grid
        xlabel('Degree')
        ylabel('Mean Square Error')
       
        title('Degree finder')
        
   
    end
    function CorrelationPressed(h,eventdata)
        Button = questdlg('Choose a color for correlation','Select color','Red','Green','Blue','Red'); 
        switch Button
            case 'Red'
                Z=1; p='r';
            case 'Green'
                Z=2;p='g';
            case 'Blue'
                Z=3;p='b';
        end
        
        figure(1)
        subplot(1,1,1)
        a=I1(:,:,Z);
        b=I2(:,:,Z);
        [x,y,z]=size(I1);
        for t=1:y
        c=xcorr(a(:,t),b(:,t));
        plot(c,p)
        xlabel('Correlation Coeffs')
        pause(0.3)
        end 
    end

function mse2=MSE1(x,y,n)
% Mean Square Error of two function
% x,y : signals wants to be compared
% n   : number of points to be compared
MS=(x-y).^2;
mse2=sum(MS)/n;
end
    function AboutPressed(h,eventdata)
        a='Copyright (c) 2012';
        b='Sagheb Kohpayeh Araghi ';
        b2='-----------------------------';
        c='Multimedia University';
        d='Malaysia';
        f=strvcat(a,b,b2,c,d);
        msgbox (f,'About');
    end
function ExitPressed(h,eventdata)
     button = questdlg('Do you want to quit?','Warning!');
     if strcmp(button,'Yes')
         close all
     end
     
        
end 
    function VelocityPressed(h,eventdata)
        prompt = {'Enter number of Rotation'};
dlg_title = 'Rotation';
num_lines = 1;
def = {'1200'};
        no= inputdlg(prompt,dlg_title,num_lines,def);
        n=str2double(no{1});
        %-------------------------------------------------
 I1=imread('gray.jpg');

J([381:482],[1:632],:)=I1([381:482],[1:632],:);


T1=I1([1:100],[1:632],:);


J=I1([381:482],[1:632],:);


J([381:482],[1:632],:)=I1([381:482],[1:632],:);


J(101:380,:,:)=I1(101:380,:,:);

figure(1)
subplot(121);imshow(I1)
xlabel('Original Image');
subplot(122);

imshow(J)
xlabel(' displacement')

% figure(2)
% subplot(121)
% imshow(T1)
% xlabel('1st piece of image')
% subplot(122)
% imshow(I1([381:482],[1:632],:))
% xlabel('2nd piece of image')
%-----------------------------------
[x,y,z]=size(J);
m=x;
n=pow2(nextpow2(m));

x=J(:,1,1);
y=fft(x,n);
power=y.*conj(y)/n;
f=0:n-1;
figure(3)

%plot(power(1:length(power)/2))
pr=real(fftshift(fft(power(200:350))));
pr=pr.^2;
plot(1:length(pr),pr)
xlabel('Frequency (Hz)')
ylabel('Power')
title('{\bf Periodogram}')
axis tight
grid
freq1=find(pr==max(pr));
ff=find(pr(freq1:100)==min(pr(freq1:100)));
%m2_peak=max(pr(freq1+ff:100));
freq=find(pr(freq1+ff:100)==max(pr(freq1+ff:100)));

Freq2=freq1+freq;
Df=Freq2-freq1;



        %-------------------------------------------------
              % frequency 
        w=2*pi*n/Df;
[r,y,z]=size(I1);
c=0.03;   % Radius (CM)
mo=c/m;         % Model Parameter  to convert to meter      
V=mo.*w.*r/2;

        mV=num2str(V);
        msg='Velocity is ';
        msg2=' m/s';
        msgbox([msg,mV,msg2],'V=')
    end

end

Contact us