DIPvelocimetry.m

%DIPvelocimetry (developed from 04-05-2010 until 11-05-2010)
%for Micro-Particle Image Velocimetry
%Marten Darmawan, Amir Faisal, Gea Oswah Fatah Parikesit
%Microfluidics Lab
%Department of Engineering Physics
%Gadjah Mada University Jogjakarta,Indonesia
[e-mail], [e-mail], [e-mail]
%Please visit our website @ http://sites.google.com/site/geaari/microfluidicsugm
%This program is supported by DIPimage library
%DIPimage can be downloaded freely at http://www.diplib.org/download
%This program is developed using MATLAB 2009a Student Version
%--------------------------------------------------------------------------------------------------------------------------------------------

%DIPimage initialization
diproot='C:\Program Files\MATLAB\R2009a\toolbox\DIPimage';
diplib=[diproot '\win32\lib'];
dippath=[diproot '\common\dipimage'];
dipimages=[diproot '\images'];

addpath(dippath);
dip_initialise;
dipsetpref('imagefilepath', dipimages);

dipimage;
clear all
%close all
%

%--------------------------------------------------------------------------
% Get the first Particle-Image file from the sequence to be analyzed
[filename,PathName] = uigetfile('*.tif','Select the first Particle-Image file from the sequence');
addpath(PathName);
cd(PathName);
if isequal(filename,0)|isequal(PathName,0)
    disp('File not found')
    return
end
%--------------------------------------------------------------------------
% Define the basename of the Particle-Image files, for analyzing
basename=[PathName filename];
basename=strrep(basename,'_1.tif','_');


% Determine the number of files to analyze
prompt = {'The starting Particle-Image file in the sequence','The ending Particle-Image file in the sequence'};
dlg_title = 'Number of files';
num_lines= 1;
def     = {'1','10'};
answer  = inputdlg(prompt,dlg_title,num_lines,def);
answer = char(answer);
answer = str2num(answer);
START=answer(1);
STOP=answer(2);

JUMP=1;

% Determine the number of ROI to analyze
prompt = {'Input the number of ROI along x-axis','Input the number of ROI along y-axis'};
dlg_title = 'number of ROI'; % "number of ROI" may not be too small (if the particle density is too low), but also not too large (in order to allow for fast computation)
num_lines= 1;
def     = {'4','4'};
answer_2  = inputdlg(prompt,dlg_title,num_lines,def);
answer_2 = char(answer_2);
answer_2 = str2num(answer_2);
user_entry_x=answer_2(1);
user_entry_y=answer_2(2);

% Main part of PIV analysis
Counter=0;
for L=START:JUMP:(STOP-1) % iteration for every *pair* of initial-final images
    filename1=[basename num2str(L) '.tif'];
    filename2=[basename num2str(L+1) '.tif'];
   
    a = readim(filename1); %image source
    b = readim(filename2); %image target

    Wx = round(size(a,1)/user_entry_x);
    Wy = round(size(a,2)/user_entry_y);

    Final_Number_ROI_x=round((size(a,1)-Wx)/Wx);
    Final_Number_ROI_y=round((size(a,2)-Wy)/Wy);
    Number_of_ROI=(Final_Number_ROI_x*Final_Number_ROI_y); 
    
    counter = 0;
    z= zeros(Number_of_ROI,4);
    for x=0:Wx:(size(a,1)-Wx) % iteration for every *ROI*, scanned along the x-axis...
        for y=0:Wy:(size(a,2)-Wy) % ...and also scanned along the y-axis
            c=a(x:(x+Wx)-1,y:(y+Wy)-1); % crop the initial image's ROI
            d=b(x:(x+Wx)-1,y:(y+Wy)-1); % crop the final image's ROI
            e=crosscorrelation(c,d);
            f=threshold(e,'fixed',150); % the threshold value depends on the images histogram
            g=f*e;
            msr=measure(f,g,{'Gravity','MaxVal'},[],1,0,0); % measuring the correlation peaks
            [Y,I]=max(msr.MaxVal); % identify the highest correlation peak
            h3=msr.Gravity(1,I); % displacement vector magnitude along the x-axis
            h4=msr.Gravity(2,I); % displacement vector magnitude along the y-axis
            xc=x+(Wx/2);
            yc=y+(Wy/2);
            h1=xc;
            h2=yc;
            h3=h3-((size(g,1)-1)/2); % x-axis coordinate of the displacement vector
            h4=h4-((size(g,2)-1)/2); % y-axis coordinate of the displacement vector
            counter = counter + 1;
            z(counter,:) = [h1 h2 h3 h4];
        end
    end
    if Counter==0
        z_new=z;
    elseif Counter>0
        z_new=z_new+z;
    end
Counter=Counter+1;
z_final=z_new/Counter;
end

% Displaying the vector field
New_a=imread(filename1);
imshow(New_a); figure(gcf)
hold on
quiver(z_final(:,1),z_final(:,2),z_final(:,3),z_final(:,4),1,'g' )
%colorbar
hold off