MATLAB Examples

Example 08a - Treating test case A1 from 3rd PIV Challenge

This case is essentially the same as example 07: The test case A1 from 3rd PIV challenge [ref. 6] is treated. Mean and rms velocities are computed, velocity PDF is determined and wavenumber spectra is calculated; these results can be compared to results shown in reference [6], Figs. 3a, 16, 18 and Table 12.

Another shown feature is the possibility to add iterations during the processing. Initially, the data are processed in four passes (with IA size 64x64, 32x32, 16x16 and 8x8 pixels). Additional pass (with 8x8 pixels interrogation size) are then performed. Effect of number of passes on velocity PDF and spectra is shown.

Reference: [6] Stanislas, M., K. Okamoto, C. J. Kahler, J. Westerweel and F. Scarano, (2008): Main results of the third international PIV Challenge. Experiments in Fluids, vol. 45, pp. 27-71.

Instructions:

  1. Download images (case A2) from PIVchallenge site,
  2. Unzip them to folder ../Data/Test PIVChallenge3A1,
  3. Run this example.

Contents

Define image pairs to treat

Initialize the variable pivPar. Get list of image files.

% Initialize variables
clear;
pivPar = [];                                              % variable for settings
imagePath = ['../Data/Test PIVChallenge3A1'];  % folder with processed images

try
    % get list of "A" images
    aux = dir([imagePath, '/*a.tif']);
    for kk = 1:numel(aux), fileList{kk} = [imagePath, '/', aux(kk).name]; end   %#ok<SAGROW>
    im1 = sort(fileList);

    % get list of "B" images
    aux = dir([imagePath, '/*b.tif']);
    for kk = 1:numel(aux), fileList{kk} = [imagePath, '/', aux(kk).name]; end   %#ok<SAGROW>
    im2 = sort(fileList);
catch
    % if error on loading images, probably images were not downloaded. Give a message.
    error(['No images found. Please, download images (case A1) from http://www.pivchallenge.org/pub05/A/A1.zip, '...
        'unzip them and place them to folder ../Data/Test PIVChallenge3A1.']);
end

Running first 4 passesof PIV processing

First, image pairs will be treated with 3-pass PIV with decreasing size of interrogation area and spacing of interrogation grid. Set the corresponding parameters:

pivPar.iaSizeX = [64 32 16  8];           % interrogation area size for three passes
pivPar.iaStepX = [32 16  8  4];           % grid spacing for three passes
pivPar.anVelocityEst = 'none';            % iterate each image pair zero
pivPar.anOnDrive = true;                  % files with results will be stored in an output folder
pivPar.anForceProcessing = false;         % if false, only image pairs, for which no file with results is
        % available, will be processed. Processing is skipped if file with results is available. If true,
        % processing is carried out even if result file is present. (Set this parameter to true if all image
        % pairs should be reprocessed, for example because of different setting of processing parameters).
pivPar.anTargetPath = [imagePath,'/pivOut1'];   % folder for storing results of first three passes
pivPar.qvPair = {'U','clipLo',-0.20,'clipHi',0.20,...
    'invLoc'};                            % plot locations of invalid velocity vectors
figure(1);
[pivPar1] = pivParams([],pivPar,'defaultsSeq');

% Run the analysis now
[pivData1] = pivAnalyzeImageSequence(im1,im2,[],pivPar1);
Checking presence of 101 output files... Finished in 0.18 s. All required files found.
 Reading results from result file. This will take a while (in my Matlab, Ctrl+C does not work at this stage)... Finished in 11.08 s. 

Running additional (5th) pass of PIV processing

Images will be processed with two additional passes, during which the interrogation area size is kept at 16x16 pixels and grid spacing at 8x8. Results will be stored in another variable, and files will be placed in different folder.

pivPar.iaSizeX = [ 8];              % interrogation area size for fourth and fifth pass
pivPar.iaStepX = [ 4];              % grid spacing
pivPar.anVelocityEst = 'pivData';   % use velocity data stored in pivData as velocity estimate used for image
                                    % deformation. By this setting, results of previous passes are transfered
                                    % to additional passses.
pivPar.ccMethod = 'dcn';            % set cross-correlation method to direct convolution (faster than default
                                    % FFT, if displacements are small, which is the case of final passes)
pivPar.anTargetPath = [imagePath,'/pivOut2']; % set output folder different than that for first 16x16 px
pivPar.qvPair = {'U','clipLo',-0.20,'clipHi',0.20};
figure(1);
[pivPar2] = pivParams([],pivPar,'defaultsSeq');

% Run the analysis now. Note that pivData1 (results of first three passes) are an input variable.
[pivData2] = pivAnalyzeImageSequence(im1,im2,pivData1,pivPar2);
Checking presence of 101 output files... Finished in 0.19 s. All required files found.
 Reading results from result file. This will take a while (in my Matlab, Ctrl+C does not work at this stage)... Finished in 11.16 s. 

Show results

First, samples of the velocity field are shown

N1 = numel(pivData1.U);
N2 = numel(pivData2.U);
Nt = size(pivData2.U,3);
% figure(1);
for kk=1:10:Nt
    hold off;
    pivQuiver(pivData2,'timeSlice',kk,...
        'U','clipLo',-0.2,'clipHi',0.2);
    title('Background: U_{mag}. Quiver: velocity (black: valid, white: replaced)');
end

Compute the mean and rms velocities, display then and compare them to the reference values.

% Statistics for the velocity fields after 4 passes:
meanU1 = mean(reshape(pivData1.U,N1,1));
meanV1 = mean(reshape(pivData1.V,N1,1));
stdU1 = std(reshape(pivData1.U,N1,1));
stdV1 = std(reshape(pivData1.V,N1,1));

% Statistics for the velocity fields after 5th pass:
meanU2 = mean(reshape(pivData2.U,N2,1));
meanV2 = mean(reshape(pivData2.V,N2,1));
stdU2 = std(reshape(pivData2.U,N2,1));
stdV2 = std(reshape(pivData2.V,N2,1));

% Print results:
fprintf('\nStatistics (8x8 px, 1 pass):   mean(U) = %+6.4f, mean(V) = %+6.4f, std(U) = %+6.4f, std(V) = %+6.4f\n',...
    meanU1,-meanV1,stdU1,stdV1);
fprintf('Statistics (8x8 px, 2 passes): mean(U) = %+6.4f, mean(V) = %+6.4f, std(U) = %+6.4f, std(V) = %+6.4f\n',...
    meanU2,-meanV2,stdU2,stdV2);
fprintf('Reference:             mean(U) = %+6.4f, mean(V) = %+6.4f, std(U) = %+6.4f, std(V) = %+6.4f\n',...
    0,0,0.0625,0.0626);
Statistics (8x8 px, 1 pass):   mean(U) = -0.0000, mean(V) = -0.0000, std(U) = +0.0572, std(V) = +0.0572
Statistics (8x8 px, 2 passes): mean(U) = -0.0000, mean(V) = -0.0000, std(U) = +0.0503, std(V) = +0.0503
Reference:             mean(U) = +0.0000, mean(V) = +0.0000, std(U) = +0.0625, std(V) = +0.0626

Compute histogram of u' and show it

% define bin range of histogram
binranges = (-0.4:0.01:0.4)';

% compute and normalize histogram
bincounts1 = histc(reshape(pivData1.U-meanU1,N1,1),binranges);
bincounts1 = bincounts1/(sum(bincounts1)*(binranges(2)-binranges(1)));
bincounts2 = histc(reshape(pivData2.U-meanU2,N2,1),binranges);
bincounts2 = bincounts2/(sum(bincounts2)*(binranges(2)-binranges(1)));

% display it
figure(2);
plot(binranges,bincounts1,'-r',binranges,bincounts2,'-k');
legend('8x8 px 1 pass','8x8 px, 2 passes');
xlabel('displacement U (px)');
ylabel('PDF (a.u.)');
title('Velocity PDF - compare to Fig. 11 in [6]');

Compute power spectra of u':

% Spectrum from the velocity field after 4 passes
uprime1 = pivData1.U-meanU1;     % velocity difference from the mean velocity
uprime1 = uprime1(1:2:end,:,:);  % reduce amount of velocity data
uspectra1 = zeros(size(uprime1,1)*size(uprime1,3),size(uprime1,2));
% Compute spectrum for every row of velocity data
for ky = 1:size(uprime1,1)
    for kt = 1:size(uprime1,3)
        uspectra1(kt+(ky-1)*size(uprime1,1),:) = abs(fft(uprime1(ky,:,kt))).^2;
    end
end
% compute mean spectrum
uspectra1 = mean(uspectra1,1);
uspectra1 = uspectra1(1:floor(numel(uspectra1)/2));
% Determane wavenumber corresponding to the spectra
dk = 1/(pivData1.X(1,end,1)-pivData1.X(1,1,1));
k1 = (0:(numel(uspectra1)-1))*dk;
% normalize spectrum
uspectra1 = uspectra1/sum(uspectra1) * std(reshape(uprime1,numel(uprime1),1),1)^2/(2*pi*dk) ;

% Spectrum from the velocity field after 5th pass
uprime2 = pivData2.U-meanU2;
uprime2 = uprime2(1:4:end,:,:);  % reduce amount of velocity data
uspectra2 = zeros(size(uprime2,1)*size(uprime2,3),size(uprime2,2));
for ky = 1:size(uprime2,1)
    for kt = 1:size(uprime2,3)
        uspectra2(kt+(ky-1)*size(uprime2,1),:) = abs(fft(uprime2(ky,:,kt))).^2;
    end
end
uspectra2 = mean(uspectra2,1);
uspectra2 = uspectra2(1:floor(numel(uspectra2)/2));
% Determane wavenumber corresponding to the spectra
dk = 1/(pivData2.X(1,end,1)-pivData2.X(1,1,1));
k2 = (0:(numel(uspectra2)-1))*dk;
% normalize spectrum
uspectra2 = uspectra2/sum(uspectra2) * std(reshape(uprime2,numel(uprime2),1),1)^2/(2*pi*dk) ;

% Display the spectra
figure(4);
loglog(2*pi*k1,uspectra1,'-r',2*pi*k2,uspectra2,'-k');
legend('8x8 px 1 pass','8x8 px, 2 passes');
xlabel('k_x (1/px)');
ylabel('E (a.u.)');
xlim([5e-3,1]); ylim([1e-4,10]);
title('Velocity spectra - compare to Figs. 3a and 8 in [6]');