Particle tracking using IDL -- this PIV program written by Doug Anderson

# epiv.propivplot.procirc.promydiff3.pro

Sometimes our images are too noisy to effectively track particles. We wrote this software to do particle image velocimetry (PIV) with those types of images. An early version (dpiv.pro) was written by Doug Anderson (andersondoug(at)gmail.com), and was slightly modified by Eric Weeks (epiv.pro). circ is called by pivplot, it makes the user-defined plot symbol (plot symbol 8) to be a circle. mydiff3 is called by epiv just to make a picture while the analysis is being done.

The idea of PIV is fairly simple. Take two pictures that are slightly different from each other. Consider a small region of the first picture, say 10 x 10 pixels. Compare that region with 10 x 10 pixel regions of the second picture that are close to where the original region was taken from. The region from the second picture that most closely resembles the first picture's region then indicates that whatever you're looking at has moved over to the position of the second region. (If the region in the second picture has the same location as the original region, then nothing has moved.) Thus, this measures a displacement vector for the position of the original region. By checking regions everywhere in your picture, you get a displacement vector field.

 DISCLAIMER: Note that there's very fancy ways to implement PIV (see here for example) and we make no claims for accuracy or reliability of our code. For example, the program tries to use subpixel accuracy but often seems to return values that are stuck at the pixel location. If we improve the program, we'll post the updated version here (or if someone else improves it, please email Eric Weeks your improved version! My email address is at the bottom of this page.)

How to use the programs

You start with two images that are (hopefully) slightly different, such as these two:

-- this is image "aa"
-- this is image "bb", the big black particle has moved slightly to the right.

IDL> help,aa,bb

aa             BYTE     = Array[101, 101]
bb             BYTE     = Array[101, 101]
IDL> piv=epiv(aa,bb,6)

This runs the PIV algorithm on the two images, allowing the comparison region in the 2nd image to shift by as much as 6 pixels in any direction relative to the location of the comparison region in the 1st image. In other words, your displacement vectors will be a maximum of 6 pixels long.

The output data has the following form: IDL> help,piv

piv            FLOAT    = Array[96, 96, 4]

Essentially, four 96 x 96 arrays. The 96 x 96 is because the program crops in a little bit from each edge, to allow the comparison region in the 2nd image to be taken with a position up to 6 pixels shifted in any direction. The arrays are:

• piv(*,*,0) is the x displacement
• piv(*,*,1) is the y displacement
• piv(*,*,2) is the magnitude |dr| of the displacement
• piv(*,*,3) is the error, letting you know how closely the two images resemble each other. Smaller values are better.

There are several keywords you can use:

• spacing = 1: The program computes a displacement vector for every pixel, by default (in both x & y: thus actually only a quarter of the pixels are considered). You can change this with the spacing keyword, for example, spacing=5 is every 5th pixel.
• stretch = 1: You can rescale the displacement vectors, for example if you want to convert from a "pixels difference between two images" to "microns per second".
• /nopoly: The program, by default, tries to determine the displacement vector with subpixel accuracy. If you're not happy about this, use /nopoly.

unsupported/undocumented keywords:

• fits = 3: an undocumented keyword, related to the subpixel accuracy calculation
• denoise = 0.4: an undocumented keyword, which relates to removing any PIV displacements which seem quite different from any of their neighboring displacements. These are assumed to be errors. Use denoise=-1 if this undocumented keyword scares you.
• /morefit: an undocumented keyword which tries to improve on the subpixel accuracy. Not well-tested (by us).
• window = 4: an undocumented keyword that relates to the size of the comparison region; the default size is the same as the maximum displacement. Needs to be an even number.

If you want to display the PIV data, a simple but crude method is to use pivplot:

IDL> pivplot,piv

(Actually, I added the keyword "charsize=2" to the command to make the size of the axis labels larger for this figure.) pivplot takes several keywords:

• spacing = 5: The program by default plots a displacement vector for every 5th data point in "piv". For many images, plotting many more displacement vectors will overwhelm your image, but often you'll want to tinker with this. For example, you may want to zoom in on a region and plot more vectors.
• /nodot: will turn off the dots for each location. Often useful if you want to use a smaller value for "spacing". The dot is located at the original position of the displacement vector. In other words, if we were fancier and drew arrows, the arrowhead would go at the opposite end of the line, and the arrowtail is the location of the dot.
• calib = 1: The program assumes that the spacing between each data point in "piv" has a length of 1. If you want to calibrate your data, use this keyword.
• scale = 1: The program by default stretches the displacement vectors to be long, although still keeping them from potentially overlapping. If you want longer (or shorter) vectors, set scale to larger than 1 (or less than 1).
• Extra keywords: you can pass along extra keywords that will be used for plot. See example below.

For example:

IDL> pivplot,piv,charsize=2,spacing=2,xrange=[20,60],yrange=[20,60]