Computer Graphics & Geometry
Dmitry Chetverikov
Computer and Automation Research Institute
Budapest, Kende u.1317, H1111 Hungary
In this paper we investigate the applicability and efficiency of feature based tracking algorithms in flow velocity measurement. Flow visualisation and measurement of flow dynamics are important for the analysis of combustion processes, hydrodynamic and aeronautical phenomena, flame propagation, heat exchange problems, construction of artificial heart pumps, and many other practical tasks [6].
In Particle Image Velocimetry (PIV) applications, the fluid is seeded with particles that follow the flow and efficiently scatter light. Conventionally, multiple exposure cameras are used to capture images of the flow at different time instants. The images are recorded on photographic film. Optical correlation methods are then applied pointbypoint to the entire negative, and the velocity of the particles between two consecutive images is determined.
Digital PIV, or DPIV, refers to the newly emerged technique of using highperformance CCD cameras and framegrabbers to store and process the digitised PIV sequences directly by the computer. In DPIV, the conventional crosscorrelation methods are implemented using the Fast Fourier Transform [11].
At the same time, estimation and tracking of motion in image sequences is a wellestablished branch of computer vision (CV). Here, real world scenes with large, rigid or deformable moving objects are usually considered. Two main classes of motion estimation methods are traditionally distinguished: the optical flow based [2] and the local feature based [4] techniques. Efficient algorithms are available in both classes. Applying them to particle flow image sequences may lead to essential improvement in visualisation and measurement results.
Recently, Quénot et al. [9] presented an optical flow algorithm based on a orthogonal dynamic programming (ODP) technique and tested the algorithm with synthetic and real DPIV sequences, including complex flows. The results compare favourably to those obtained using the classical correlation based DPIV methods. A substantial increase in both accuracy and spatial resolution of the velocity field has been achieved. Superior robustness to artificially generated noise has also been demonstrated.
Unfortunately, the algorithm [9] is very slow, with the execution time increasing drastically with image size. This motivates our search for reasonable alternatives to both the conventional and the DP techniques. In this paper, we present the feature tracking approach to PIV, which includes novel techniques for coherence filtering and interpolation of a velocity field. Quantitative results of flow estimation are presented. The tests compare three groups of approaches to DPIV: those based on correlation, ODP, and feature tracking.
2 Feature tracking algorithms applied to PIV
Feature tracking techniques extract local regions of interest (features) and identify the corresponding features in each image of a sequence. In this section, we outline two particular algorithms, the KLT Tracker [10] and the IPAN Tracker [4], which we apply to Particle Image Velocimetry.
The KLT Tracker selects features which are optimal for tracking, and keeps track of these features. A good feature is a textured patch with high intensity variation in both x and y directions, such as a corner. The algorithm defines a measure of dissimilarity that quantifies the change in appearance of a feature between the first and the current frame, allowing for affine distortions. At the same time, a pure translational model of motion is used to estimate the position in the next frame.
Consider the intensity function g(x,y). The approach assumes small displacements and minimises the linearised dissimilarity. The displacement vector d = (d_{x},d_{y})^{T} is computed as the solution of the linear system [10] or, in matrix notation, Zd = e. (g_{x} is a partial derivative of g(x,y).) It is implied that Z and e = (g_{x},g_{y})^{T} are integrated over a feature window W. The size of W is normally set to 25×25. A patch defined by the window is accepted as a candidate feature if both eigenvalues of Z, l_{1} and l_{2}, exceed a predefined threshold l: min(l_{1},l_{2}) > l. This ensures that the matrix Z is wellconditioned and the solution of (1) is accurate. When applied to a PIV image, the KLT selects individual particles as feature centres. To cope with relatively large motion, the algorithm is implemented in a multiresolution way.
Each feature being tracked is monitored to determine if its current appearance is still similar, under affine distortion, to the initial appearance observed in the first frame. When the dissimilarity exceeds a predefined threshold, the feature is considered to have been lost and will no longer be tracked. At least two frames are needed for the operation of the algorithm.
Since the KLT algorithm incorporates an analytical solution to motion estimation, it is much faster than the methods that use explicit region matching, such as the conventional crosscorrelation and the ODP techniques. The source code of the tracker can be downloaded from the web site [3].
The IPAN Tracker is a noniterative, competitive feature point linking algorithm. Its original, applicationindependent version tracks a moving point set, tolerating point entry, exit and false and missing measurements. Position is the only data assigned to a point. The algorithm is based on a repetitive hypothesis testing procedure that switches between three consecutive image frames and maximises the smoothness of the evolving trajectories.
The applicationindependent version of the tracker is described in full detail in our recent paper [4]. When applied to PIV, the operation of the algorithm remains basically the same. The modifications necessary for the PIV application are as follows: 1. a PIVspecific feature selection mechanism is added; 2. the cost function is modified to include feature appearance. These modifications are described below.
Selection of feature points. A PIV image g(x,y) is smoothed by a 3×3 mean filter and the (realvalued) maxima of the smoothed image s(x,y) are selected as the features. Each bright particle is represented by a maximum of s(x,y). For more precise motion estimation, the position of each maximum is corrected by parabolic interpolation in x and y directions separately. Denote the corrected coordinates of a maximum (x�,y�) by (x_{f},y_{f}).
A feature point P(x_{f},y_{f}) is assigned a feature value

(1) 
f(x_{f},y_{f}) characterises the dominance of a feature so that large and bright particles have high dominance. The features (particles) are ranked according to their dominance and a specified number of the most dominant features are selected for tracking. Note that the KLT and IPAN algorithms select different points.
Cost function. The IPAN Tracker processes sequentially each three consecutive frames F_{n1}, F_{n} and F_{n+1}, where n is the current frame. Consider a feature point in each of the 3 frames, A � F_{n1}, B � F_{n} and C � F_{n+1}. When applied to PIV, the tracker minimises the cost function d(A,B,C) which accounts for changes in velocity and appearance of a feature:

(2) 

F(A,B,C), the appearance term of d(A,B,C), accounts for changes in feature appearance in a 3×3 neighbourhood. It is defined as

The cost function d(A,B,C) is minimised using the repetitive hypothesis testing procedure presented in [4]. At least three frames are needed for the operation of the algorithm.
3 Postprocessing of velocity vector field
Feature trackers may occasionally yield completely wrong velocity vectors. To enhance the result of measurement, coherence based postprocessing is applied to the `raw' velocity field obtained by the trackers. The coherence filter modifies a velocity vector if it is inconsistent with the dominant surrounding vectors. The solution we use is a modified version of the vector median filter [1]. The procedure operates as follows.
Given a feature point P_{c} with the velocity vector v_{c}, consider all features P_{i}, i = 1,2,�,p, lying within a distance S from P_{c}, including P_{c} itself. Let their velocities be v_{i}. Due to coherent motion, these vectors are assumed to form a cluster in the velocity space. Introduce the mean cumulative difference between a vector v_{i} and all other vectors v_{j}, j � i:

(3) 

(4) 

(5) 
Uniform sampling of the measured velocity field is normally required. A number of techniques [5] are available for resampling the results of a measurement. However, most of them perform resampling from one regular grid to another. We use the following procedure.
Given a point, G, on the required regular grid, consider all feature points P_{i}, i = 1,2,�,p, lying within a certain distance S from G. Let their velocity vectors be v_{i}. Denote by d^{2}_{i} the squared distance from P_{i} to G and by [^(d^{2})] the mean of d^{2}_{i} over all i. Introduce


(6) 
The figures below show examples of coherence filtering and resampling.
The KLT and IPAN Trackers were run on a large number of standard synthetic and real DPIV sequences available on the Internet. The postprocessing algorithms described in section 3 were applied to the initial velocity fields obtained by the KLT and IPAN Trackers. In this section, the two algorithms complemented by the postprocessing will be referred to as KLTPIV and IPANPIV, respectively. They will be compared to the correlation (CORRPIV [11]) and the Orthogonal Dynamic Programming (ODPPIV [9]) algorithms.
Two sets of synthetic test data with the ground truth were used in our experiments. The two datasets were downloaded from the Japanese PIVSTD Project web site [7] and the Quénot's web site [8].
In the tables below, the error is the mean absolute deviation from the ground truth; the variance of the absolute deviation is also given. Typical times elapsed during execution are also shown, measured on a Pentium 333 MHz under Linux OS. The execution times are only given to indicate the order of magnitude, as they depend on implementation and hardware.
The Visualisation Society of Japan has launched a project named PIVSTD [7] whose goal is to generate and distribute standard synthetic test sequences for performance evaluation of PIV algorithms. The web site [7] offers a collection of pregenerated sequences as well as a program allowing to produce sequences with desired parameters. The basic tunable parameters are the number of particles N_{p}, the average velocity v, the average particle diameter P_{a} and its standard deviation P_{d}.
An example of PIVSTD flow is given in the above figure. Table 1 compares displacement errors of the three methods, IPANPIV, KLTPIV and ODPPIV2, for six STDPIV sequences. Results of CORRPIV were not available for this dataset. ODPPIV2 is the most precise of the three variants presented in [9].
s01  s03  s04  s05  s06  
IPANPIV  4.72 � 4.3  10.3 � 4.5  4.86 � 6.2  5.41 � 4.9  4.59 � 5.7 
KLTPIV  3.37 � 2.3  9.82 � 4.0  2.02 � 2.1  3.51 � 2.7  2.02 � 1.8 
ODPPIV2  3.52 � 2.9  9.82 � 5.1  1.97 � 2.7  2.68 � 3.1  1.50 � 1.6 
To assess the PIVSTD results, consider table 2 that describes the six sequences. In this table, `standard' means the standard values which are N_{p} = 4000 particles, P_{a} = 5.0 pixel, P_{d} = 1.4 pixel, and v = 7.39 pixel/frame. `Slow' indicates a slow flow with v = 2.24. `Same particles' means that P_{d} = 0 and the particles are indistinguishable. Other parameters have the standard, or close to the standard, values.
Sequence  N_{p}  Description 
s01  4000  standard 
s03  4000  slow 
s04  10000  dense 
s05  1000  sparse 
s06  4000  same particles 
For the STDPIV dataset, the performance of KLTPIV is very close to that of ODPPIV. The accuracy of IPANPIV is lower. The high density and the indistinguishable particles pose some problem to IPANPIV. Note that the ODPPIV algorithm is much slower than the trackers, as discussed below.
Table 3 compares IPANPIV, KLTPIV, ODPPIV and the correlation method CORRPIV for a set the synthetic flow sequences called CYLINDER [8]. Typical execution times are indicated in the bottom row of the table.
IPANPIV  KLTPIV  CORRPIV  ODPPIV1  ODPPIV2  
N0  0.42 � 0.5  0.35 � 0.6  0.55 � 1.0  0.13 � 0.1  0.07 � 0.1 
N5  0.50 � 0.7  0.35 � 0.6  0.61 � 1.2  0.21 � 0.5  0.08 � 0.1 
N10  0.55 � 0.7  0.33 � 0.5  0.77 � 1.6  0.53 � 1.4  0.11 � 0.1 
N20  0.90 � 1.2  0.44 � 0.6  3.11 � 4.1  0.88 � 1.6  0.20 � 0.1 
Time  20 sec  15 sec  10 min  20 min  200 min 
ODPPIV1 is the fastest and least precise of the three variants presented in [9]. CORRPIV is a 32 ×32 window size correlator. N5, N10 and N20 are noisy versions of the original noisefree sequence N0. The numbers indicate the degrees of noise varying from 5% to 20%. CYLINDER is a complex flow with the mean displacement 7.6 pixel/frame.
For the CYLINDER dataset, the feature trackers outperform the correlation algorithm in both accuracy and execution time. For noisy (and more realistic) data, they compete with the fastest variant of ODPPIV in accuracy and are definitely superior in processing speed. ODPPIV2 is very accurate, but its computational load is almost prohibitive. In addition, KLTPIV is more robust to noise, as its accuracy deteriorates slower than that of ODPPIV2.
We have presented a novel approach to Particle Image Velocimetry based on feature tracking. Two efficient algorithms have been customised to PIV by adding the procedures for coherence filtering and resampling. The procedures have been extended to cope with flow discontinuities. The coherence filtering improves the accuracy of velocity estimation. The resampling provides a uniform sampling at the expense of a moderate decrease of the accuracy.
The results of the tests demonstrate that the proposed approach offers a good alternative to both correlation and ODP techniques. KLTPIV and IPANPIV provide higher flow estimation accuracy and are faster than the conventional correlation techniques. For noisy images, the feature tracking PIV provides accuracy comparable with that of the precise ODP algorithms, while requiring much less computation. The processing speed of the trackers can potentially make them suitable for fast flow visualisation, qualitative estimation, and analysis of timevarying flows.
The operation of the trackers is based on detection of prominent features. If no features are found in a region of the flow, interpolation provides measurements in this region as well. Such method may fail when a large area of poor visibility appears in the image because of low contrast or blur. This problem needs further research.
The most interesting, open question is that of flow complexity. Simple, smooth flows like PIVSTD can be in most cases reliably measured by all methods considered in this paper. Future research should focus on flow complexity and involve creation of realistic but complex test sequences with the ground truth provided.
Online demonstrations of IPANPIV and KLTPIV are available on the Internet at the web site of the IPAN research group:
http://visual.ipan.sztaki.hu/demo/demo.html
A remote user can select an algorithm, set the parameters and run the algorithm on a short PIV sequence. The flow sequences used in our tests are provided together with the ground truth, when available. The quantitative results of our study can therefore be checked online. Alternatively, the user can upload and process his/her own data and compare the result to the submitted ground truth.
This work was supported by the Hungarian Scientific Research Fund under grants OTKA T026592 and M28078. The author thanks Marcell Nagy and Judit Verestóy for their contributions to the implementation and testing of the algorithms described in this paper.
Computer Graphics & Geometry