Test case #352 PIV STD project

Synthetic data were used to validate the algorithm. The synthetic data were from the PIV-STD library that was developed by the Visualization Society of Japan and served as the benchmark for many PIV/PTV research efforts (Okamoto et al 2000, Ruhnau et al 2005, Kim and Sung 2006). Data set 352 was chosen, which describes a jet flow impinging on a wall. The flow is transient and three dimensional. Data from Camera 0, frame 0–89 were adopted. At each frame, about 300 particles were in the observed region, which has a size of 256 by 256 pixels. The fluid motion was captured by three cameras. The particle spacing displacement ratio (PSDR), which is defined as the ratio of the average particle spacing to the mean particle displacement between two consecutive frames and serves as an important indicator to the difficulty of the tracking process (Malik et al 1993), is about 3.4 for thewhole image domain. However, compared to the left half of the images, the right half contained denser particles and larger particle displacement (see figures 4 and 5), which results in a lower PSDR. The PSDR is about 2.6 for the right half of the image domain and about 2.2 for the right quarter of the image domain. Such a small value makes it hard to track particles through 90 frames (Malik et al 1993).

[A multi-frame particle tracking algorithm robust against input noise, by Dongning Li, Yuanhui Zhang, Yigang Sun andWei Yan, Meas. Sci. Technol. 19 (2008) 105401 ]

3D-PTV software applied to the test case #352 (Standard PIV Images project)

Here we report the technical procedures, the results and the analysis of the case #352 processed using 3D-PTV software (ETH Zurich, Tel Aviv University).

Download

1. Download images and the data for the test case #352 from the original web site http://piv.vsj.or.jp/piv/image3d/image352.html (just in case the website/data will disappear, we'll keep a local copy on our web server)

2. Obtain 3D-PTV software from the website: http://ptv.origo.ethz.ch (for Windows platforms, if no source is needed, the ZIP file below includes the executable, so it's plug-n-play ready test case)

Pre-processing procedure

Download the ZIP file with the directory ready for use with the software. Extract the file into the experiment directory, e.g. C:/PTV/test352

verify it has the directory structure:

C:/PTV/test352

       /img

       /res

       /cal

       /parameters

read more about the software handling on http://ptv.origo.ethz.ch/wiki/handling_ptv_software

1. convert RAW images to TIFF images (no compression, 8 bit). On Windows , IrfanView (with plugins) used for batch image conversion and rename. 

image0XXX.raw -> cam1.1XXX (1XXX= 1000 to 1144). Note that we also name our cameras: 1,2,3 as the PIV STD project calls them 0,1,2

move the images into the the sub-directory /img

2. We already converted the calibration images: im0999.raw -> cam1.tif, im1999.raw -> cam2.tif, im2999.raw -> cam3.tif

and added the calibration images to the sub-directory /cal

Basically, we convert the information provided for this test case into the calibration files, compatible with the 3D-PTV software. The camera configuration is well explained on the Evaluation of standard PIV images website by Georges Quenot:

From: http://www-clips.imag.fr/mrim/georges.quenot/vsj-eval/evaluation3.html

or from the original web site:

This is precisely the same coordinate system as used by the 3D-PTV software. Therefore the conversion of the parameters is easy:

in the /cal directory you'll find the *.ORI files for each camera. These consist of:

   101.8263 -9.9019 65.1747 projective center X,Y,Z, [mm]
   0.4151383 -0.0069793 1.5073263 omega, phi, kappa [rad]

   0.0634259 -0.9979621 -0.0069792 rotation matrix (3x3)

   0.9130395 0.0608491 -0.4033067 [no unit]

   0.4029095 0.0192078 0.9150383

   -0.6139 -0.0622 xp, yp [mm]

   8.7308 principle distance [mm]

   0.0 0.0 30.0 window (glass) location [mm] 

so the three cameras are oriented approximately: 

cam1.tif.ori (camera #0)

     -100.3781      0.0117    173.0085
        0.0000154  -0.5432731  -0.0000624

        0.8560213  0.0000534 -0.5169406
       -0.0000704  1.0000000 -0.0000132
        0.5169406  0.0000477  0.8560213

         0.0000   0.0000
        20.0000

          0.000000000000000    0.000000000000000    30.000000000000000 

cam2.tif.ori

         0.4331     -0.0829    199.4893
        0.0005046   0.0022530   0.0001731

        0.9999974 -0.0001731  0.0022530
        0.0001743  0.9999999 -0.0005046
       -0.0022529  0.0005050  0.9999973

         0.0000   0.0000
        20.0000

          0.000000000000000    0.000000000000000    30.000000000000000

cam3.tif.ori

      101.0432      0.0872    172.5569
       -0.0004084   0.5476053   0.0000590

        0.8537737 -0.0000503  0.5206442
       -0.0001537  0.9999999  0.0003487
       -0.5206442 -0.0003777  0.8537737

         0.0000   0.0000
        20.0000

          0.000000000000000    0.000000000000000    30.000000000000000

Calibration file:

   Calibration File (imX999.raw)
   Particles (Total 27 particles)
   x = -0.8, 0, 0.8 cm
   y = -0.8, 0, 0.8 cm
   z = -0.8, 0, 0.8 cm

is converted into a text file in the format, compatible with the 3D-PTV software - four columns of: ID, X, Y, Z, using the Matlab code (attached to the ZIP file):

   1    -8    -8    -8
   2    -8     0    -8
   3    -8     8    -8
   4     0    -8    -8
   5     0     0    -8
   6     0     8    -8
   7     8    -8    -8
   8     8     0    -8
   9     8     8    -8
  10    -8    -8     0
  11    -8     0     0
  12    -8     8     0
  13     0    -8     0
  14     0     0     0
  15     0     8     0
  16     8    -8     0
  17     8     0     0
  18     8     8     0
  19    -8    -8     8
  20    -8     0     8
  21    -8     8     8
  22     0    -8     8
  23     0     0     8
  24     0     8     8
  25     8    -8     8
  26     8     0     8
  27     8     8     8

Processing (particle locations, stereo-matching, 3D positions, tracking) is then according to the standard 3D-PTV procedure. All the parameters are given in the ZIP file and reproduced here as a snapshot of the respective menus:

Following 3D-PTV tutorials we proceed with:

Start -> Tracking -> Sequence/Tracking to complete the full cycle of processing. The results appear in two subdirectories:

2D particle locations in pixels, per camera in /img directory, named cam?.1XXX_targets

3D particle locations and tracking information in /res directory named ptv_is.1XXX

Comparison of results

It's difficult to build very generic code for comparing the results of all the possible types of software and data. We suggest to develop the attached code further, so it will work for various data sets and software packages. The basic subroutine i:

function compare_2d(dir1,dir2,type1,type2,num1,num2)

% COMPARE_2D(DIR1,DIR2,TYPE1,TYPE2,NUM1,NUM2) compares two ASCII files

% TYPE variable can be 'eth' or 'std'

%

%

% Usage:

% dir1 = 'c:\PTV\Software\sortgrid_from_initial_guess_Oct09\test\img';

% type1 = 'eth';

% dir2 = 'c:\PTV\Software\sortgrid_from_initial_guess_Oct09\ptc';

% type2 = 'std';

% num1 = 1;

% num2 = 144;

% >> pixelshift = []; % or 'shift' if 0.5 pixel shift is needed for ETH vs. Wesleyan

% >> compare_2d(fname1,fname2,pixelshift)

%

% Author: Greg Voth, Wesleyan Uni., August 2009.

% Modified: Alex Liberzon, 23.10.09

% for the PIV STD benchmark case

%

%%

% initialization

Y = 0;

output = [];

switch type1

   case{'eth'}

       %         fname1 = 'ptv_is.1%03d';

       fname1 = inline('fullfile(dir1,sprintf(cam%d.1%03d_targets,camid,i))','dir1','camid','i');

       loadfun1 = @loadeth2d;

   case{'std'}

       %         fname1 = inline('fullfile(dir1,sprintf(ptc%d%03d.dat,camid-1,i))','dir1','camid','i');

       fname1 = inline('fullfile(dir1,sprintf(ptc%03d.dat,i))','dir1','camid','i');

       loadfun1 = @loadstd2d;

   otherwise

       error('Unknown type 1');

end

switch type2

   case{'eth'}

       %         fname2 = 'ptv_is.1%03d';

       fname2 = inline('fullfile(dir1,sprintf(cam%d.1%03d_targets,camid,i))','dir1','camid','i');

       loadfun2 = @loadeth2d;

   case{'std'}

       fname2 = inline('fullfile(dir1,sprintf(ptc%03d.dat,i))','dir1','camid','i');

       loadfun2 = @loadstd2d;

   otherwise

       error('Unknown type 2');

end

for camid = 1:3

   for num = num1:num2

       file1 = fname1(dir1,camid,num);

       file2 = fname2(dir2,camid,num);

       [x1,y1] = feval(loadfun1,file1,camid);

       [x2,y2] = feval(loadfun2,file2,camid); % note there're zeros in STD data

       npart1 = length(x1);

       npart2 = length(x2);

       matchraw = zeros(npart1,1);

       bestraw = 1.0e6*ones(npart1,1);

       for i=1:npart1

           for j=1:npart2

               del=sqrt((x2(j)-x1(i))^2 + (y2(j)-y1(i))^2);

               if del < bestraw(i)

                   bestraw(i) = del ;

                   matchraw(i) = j;

               end

           end

       end

       matchbkg = zeros(npart2,1);

       bestbkg = 1.0e6*ones(npart2,1);

       for i=1:npart2

           for j=1:npart1

               del=sqrt((x1(j)-x2(i))^2 + (x1(j)-x2(i))^2);

               if del < bestbkg(i)

                   bestbkg(i) = del ;

                   matchbkg(i) = j;

               end

           end

       end

       goodmatch=zeros(npart1,1);

       for i=1:npart1

           if matchraw(i) > 0

               if matchbkg(matchraw(i)) > 0

                   if matchbkg(matchraw(i)) == i

                       goodmatch(i)=1;

                   end

               end

           end

       end

       bestraw(goodmatch ~= 1) = -1;

       % match(f) = 0;

       [y,xout] = hist(bestraw,31); % here some work is needed to update the x limits

       %  accumulate the results for all the frames

       Y = Y + y;

       fg = ((bestraw =0));

       output = cat(1,output,[std(bestraw(fg)),sum(fg),npart1,npart2]);

   end % of the loop for all frames

   %% Report

   figure

   plot(xout(2:end-1),smooth(Y(2:end-1)),'b-','LineWidth',2)

   title('Histogram of position differences');

   xlabel('position difference (pixels)');

   ylabel('counts');

   legend(sprintf('\n%f pixel error in %d positions matched out of \n%d particles in set 1 and %d in set 2\n\n',mean(output(:,1)),sum(output(:,2)),sum(output(:,3)),sum(output(:,4))));

   %std_dev_good_matches = std(best(fg))

   % diff = best(fg);

end

With two subroutines to load the data from ETH format (camX.XXX_targets) and PIV STD format (ptcXXX.dat):

loadeth2d.m

       function [x,y] = loadeth2d(fname,camid)
       % Loading ETH 2d format files
       fid = fopen(fname,'r');
       xy = textscan(fid,'%d%f%f%d%d%d%d%d','Headerlines',1);
       x = xy{2};
       y = xy{3};
       fclose(fid);

loadstd2d.m

   function [x,y] = loadstd2d(fname,camid)

       % Loading PIV STD 2d format files, ptc???.dat
       %{ 
       % from the project site
       Particle Files  (ptc???.dat)
                                #1 camera     #2 camera     #3 camera  intensity
      ID    x      y      z      X      Y      X      Y      X      Y
       4 -0.301  0.584  0.893 147.61  64.91  95.55  65.10  53.07  65.92 208.18
       6 -0.898  1.266 -0.258   0.00   0.00   0.00   0.00  46.22   0.15 208.69
       8  1.178  0.886  0.149 249.88  38.07 251.10  35.42 219.99  32.46 191.52
       9  0.720 -1.242 -0.530 183.76 252.88 201.24 254.40   0.00   0.00 211.21

      unit: x,y,z  [cm]
            X,Y    [pixel]

      The particle ID is unique number, therefore, the same particle could
      be easily tracked with checking the ID for serial particle files.

       %}
       fid = fopen(fname,'r');
       xy = textscan(fid,'%d%f%f%f%f%f%f%f%f%f%f');
       x = xy{(camid-1)*2+5};
       y = xy{(camid-1)*2+6};
       fclose(fid); 

The comparison of the 2D results is presented in the form of histograms of the pixel differences between the two files, and textual information regarding the number of particles found, related to the number of particles in 'theory'.