Research on PIV Algorithm for Matching of Particle Images Based on Affine Transformation


 For the problems of distortion and rotation in the matching of particle images of turbulent motion, according to the nature of affine transformation, using log-polar coordinate transformation, the matching is achieved by performing correlation calculations on the image line by line, and developed a matching algorithm (Turbulent Particle Image Matching, abbreviation: TPIM) for particle image pairs with affine transformation and rigid body transformation: by moving the interpretation window, the algorithm is no longer restricted by displacements of particles; by setting the affine lines according to the angle of the image in the log-polar coordinate system and using the affine line as the matching unit, the decoupling of different transformation factors is realized; according to the characteristic of non-uniform sampling in log-polar coordinate transformation, based on the principle of not losing image information, by reasonably setting the image mask and the rate of sampling, establishing the image pyramid and the relative coordinate system, the algorithm complexity is reduced to about 15% of the original. The experimental results of various types of particle images show that the matching accuracy of the TPIM algorithm can reach more than 99%.


Introduction
Particle Image Velocimetry (abbreviation: PIV) is an undisturbed, transient, full-field measurement method [1] , and this technology uses the interpretation window of the adjacent small area as the matching unit, and calculates the displacement of fluids by matching the interpretation window in two consecutive frames of images [2] . It should be noted that, the algorithm of image matching is the most core part of the PIV technology, and the performance of the algorithm directly determines the practicality and economy of the PIV technology. Among various algorithms of image matching, the algorithm of cross-correlation based on fast Fourier transform (abbreviation: FFT) is the most widely used in PIV technology, and the main reason is that the algorithm of butterfly in FFT can greatly reduce the amount of calculation and improve the real-time performance of PIV [3,4] . In the later period, related researchers improved the algorithm of cross-correlation, and proposed algorithm of normalized cross-correlation and algorithm of de-averaging normalized cross-correlation, which improved the algorithm's robustness to illumination and noise, and expanded the scope of application of the algorithm [5,6] . However, there are still defects in matching algorithm of particle image based on cross-correlation: the algorithm divides the particle images to be matched into interpretation windows (the position of the interpretation window does not change), then calculates the matrix of correlation coefficients by cross-correlating the gray values of the pixels in the interpretation window pair, finally, the position corresponding to the largest correlation coefficient is used as the best matching position. Since the position of the interpretation windows does not change during the matching process, in order to ensure that the two interpretation windows to be matched have sufficient similarity (due to the movement of the particles, the particle information in the interpretation windows to be matched are different), the PIV technology based on the cross-correlation algorithm stipulates the criterion named 1/4: the horizontal and vertical displacements of particles in two adjacent frames of images cannot be greater than 1/4 of the side length of the interpretation window [7] .
Under the situation that the cross-frame time of the camera is fixed, this criterion directly limits the range of velocity measurement of PIV technology, that is, the application range of PIV technology. In addition, because the PIV technology uses the average value of the motion states of the particles in a small area to represent the motion states of all the particles in the area, the PIV technology is only suitable for flow fields with small velocity gradients and small fluid deformations. As the velocity gradient in the flow field increases, the relative displacements between the particles in the interpretation window increase. At this situation, if the interpretation windows with fixed shape are still used for image matching, it will inevitably reduce the correlation coefficients and reduce the matching accuracy.
For the situation that the range of velocity measurement is limited, the matching process can be completed by setting search area in the second frame of image (the size of the search area can be determined according to the estimated velocity of the flow to be measured), and then the interpretation window in the first frame of image traverses each pixel in the search area. However, this method will increase the complexity of the algorithm while making the PIV technology lose its applicability to the turbulent velocity field with rotational transformation (by shifting the interpretation window, the image pair with rotation transformation cannot be overlapped; for the search area and the interpretation window have different information, the frequency spectrums are different, it's unable to use FFT to match the image pair with rotation transformation, too). For the situation where the velocity gradient is too large and the image is deformed, related scholars proposed a kind of matching algorithm based on image deformation: the algorithm performs multiple matching on the image pair, and uses the current data of the velocity as a guide to deform the image, and then matches the deformed images until the difference of velocity between two adjacent calculations converges to the expected range [8,9] . Although this kind of algorithm can improve the applicability of PIV technology to flow fields with large velocity gradients, it greatly increases the complexity, and the algorithm may not necessarily converge during the iterative calculation process.
For the coupling of scaling, rotation and translation in particle images, Matthew N Giarra et al.
proposed a matching algorithm which robust to rotation and scaling based on Fourier Merlin transform in 2015 [10] . the algorithm uses the properties of the Fourier transform：translation-invariant，rotationinvariant and scaling-opposite, by performing log-polar coordinate transformation on the frequency spectrum of the images (removing the translation factor in the phase), then the decoupling of the translation factor, rotation factor and zoom factor is completed, and the matching of particle images with large velocity gradients is achieved. However, the algorithm defaults that the scale transformation of the image is isotropic, which means the algorithm is only suitable for matching between images zoomed in the same proportion in each direction, and for deformed images (different zoom ratios in each direction), the algorithm fails to match. However, the zoom ratio in each direction of the actual particle image of turbulent is not the same.
Based on the above analysis, according to the characteristics of particle images of turbulence, combined with the properties of log-polar coordinate transformation and affine transformation, this paper proposes an algorithm named TPIM (Turbulent Particle Image Matching).

Proposed method
For image pair with affine transformation, the coordinate relationship between the corresponding pixels in the images can be expressed by formula (1) [11] : In formula (1), (x, y) represents the coordinates of the pixel in the first frame of image; (m, n) represents the corresponding position of the pixel in the first frame of image in the second frame of image after affine transformation; a0, a1, a3, and a4 are the factors of affine transformation to be determined; a2 and a5 are the factors of rigid body transformation to be determined.
Since the TPIM algorithm matches by shifting the interpretation window within the search area, it can be considered that the factors of rigid body transformation in formula (1) satisfy: 2 Which equivalents to： In formulas (4) and (5), R represents the distance of the pixel in the first frame of image relative to the pole in the polar coordinate system; α represents the angle between the pixels and the polar axis.
When performing logarithmic coordinate transformation on formula (5) a a a a n m a a a a It can be derived from formula (6) In formula (6) and formula (7), log r and  respectively represent the logarithmic distance and offset of angle of the pixels in the image relative to the origin of the coordinate after the log-polar coordinate transformation; log r  and   respectively represent the logarithmic distance and offset of angle between two adjacent frames of images in the log-polar coordinate system.
As shown in formula (7), the parameters of the image pair to be matched are not decoupled after the log-polar coordinate transformation, so the Fourier Merlin transform can't be used to match particle Images. In addition, it can be seen from formula (7) that the value of log r  is only related to α: when α remains unchanged, log r  also remains the same, that is, relative to the origin of the coordinate system, the pixels under the same angle have the same logarithmic offset.

Affine transformation
According to the principle of affine transformation, the following properties can be obtained [12] : a. Collinearity (the collinear point set remains collinear after affine transformation); b. Parallelism (two parallel lines remain parallel after affine transformation); c. Invariance of collinear ratio (after affine transformation of two line segments on the same straight line, the ratio of their lengths remains unchanged).
According to the above properties, combined with classical geometric theory, it can be deduced that affine transformation can map a circle to an ellipse.
(a) (b) Fig. 1 Affine transformation of the circle. a represents the concentric circles before affine transformation, b represents the concentric ellipses after affine transformation Fig. 2 The results of log-polar transformation on an image pair with affine transformation. a: performing log-polar coordinate transformation on concentric circles; b: performing log-polar coordinate transformation on concentric ellipses It can be seen from Fig. 1 and 2 that the log-polar coordinate transformation can map concentric rings to parallel straight lines, and concentric elliptical rings can be mapped to parallel curves. For image matching, in the log-polar coordinate system, only the curve corresponding to each straight line needs to be determined separately. It can be clearly seen from Fig. 2 that after affine transformation and log-polar coordinate transformation, the parallel lines relative to the rings and the parallel curves relative to the ellipse rings are no longer relevant, and traditional matching algorithms can't match this type of images.
It can be seen from formula (7) that after affine transformation and log-polar coordinate transformation, in the two frames of images to be matched, the displacements between pixels under the same angle remain the same. Assuming that the two frames of images after log-polar coordinate transformation are named f and g respectively, and the rows and columns of the images represent the logarithmic distance and offset of angle respectively. According to the above analysis, in the log-polar coordinate system, the displacements of pixels in the same row are the same. According to the above analysis, in order to achieve the matching between the images, the image can be divided into affine rays according to the offset angle (that is, each line in the image is set as an affine ray), since the zoom ratios of the line segments on the same affine ray keep the same (as shown in Fig. 3), affine rays can be used as the matching unit to cross-correlate the images to achieve the matching process [13] . Fig. 3 The log-polar coordinate transformation of affine rays. a and b are image pair to be matched, after performing log-polar coordinate transformation on a and b to form c and d.

Algorithm optimization
The complexity of the algorithm is one of the important indicators to measure the performance of the algorithm, which directly determines the practicality of the algorithm. In order to achieve the matching between image pairs with affine transformation, the TPIM algorithm uses log-polar coordinate transformation, however, due to the characteristic of non-uniform sampling of the log-polar coordinate transformation (the pixels close to the origin of the coordinate are stretched, and the pixels far away from the origin are compressed), in order to ensure that the image pair after coordinate transformation can be successfully matched, it is necessary to ensure that the information in the image will not change with the coordinate transformation (or ensure that the change of information is within an acceptable range), in turn, the size of image after coordinate transformation is greatly increased, which greatly increases the complexity of the algorithm, and reduces the practicality of the algorithm [14] . In order to improve the practicality of the algorithm, the complexity of the algorithm should be reduced.

Setting resolution of angular
For the matching between image pairs with rotation transformation, the images must be traversed in the range of 0 -360   . In addition, in order to ensure the matching accuracy, the algorithm has to increase the angular resolution. For example, to ensure that the matching accuracy of the algorithm reaches 0.1 , the resolution of angular after log-polar coordinate transformation cannot be less than 3600.
PIV technology is mainly used for measurement of fluid velocity, and its matching process is based on the assumption of fluid continuity [15] . This technology thinks that the flow patterns of fluid particles in the same small area are similar, so the average value of the velocity of the particles in the area can approximately represent the motion state of all the particles in the area. Based on the above analysis, before performing image matching, the matching algorithm of PIV will divide the particle images into small area interpretation windows, and then takes the interpretation window as the matching unit. In order to meet the assumption of similarity and ensure that the interpretation window contains enough information [16] , the range of the side length of the interpretation window is generally set to 16-64 (unit: pixel). According to the above analysis, if the matching accuracy is to be ensured, the size of the interpretation window will be expanded by nearly 56 times after log-polar coordinate transformation, which means that the algorithm complexity is increased by nearly 56 times. To ensure the practicability of the TPIM algorithm, it is necessary to reduce the resolution of angular.
Through analysis, it is found that since the purpose of using PIV technology is to reconstruct the velocity field of the flow to be measured, the relevant matching algorithm only needs to correctly calculate the displacements of the corresponding interpretation windows in two adjacent frames of images. There is no need to calculate the rotation angle of the interpretation window (in addition, since the actual rotation center of each interpretation window can't be determined, the rotation angle of the interpretation window has no practical meaning). Therefore, when performing log-polar coordinate transformation on the particle image to be matched, the setting of the resolution of angle only needs to ensure that the pixels in the original image will not be lost due to coordinate transformation.

Fig. 4 Schematic diagram of the intervals of angular between adjacent rays in the interpretation window
In the interpretation window shown in Fig.4 (assuming that the center of the interpretation window is the origin of the log-polar coordinate transformation, and the interpretation window is square), the intervals of angular between adjacent rays (green) are different, and the maximum interval is θ1, and the minimum interval is θ2. Therefore, sampling the image in the range of 0 -360   with θ2 as the unit can ensure that the information of the images to be matched will not be lost. Suppose the side length of the interpretation window is 2n+1, then: In formula (8), Δθ is the sampling rate of the angle, and RA is the resolution of angular after logpolar coordinate transformation. It is clear that RA is positively correlated with n, and the result of calculation shows that, when the side length of the interpretation window is within the range of 16-64, the size of image after log-polar coordinate transformation is expanded by nearly 5 times.

Setting image mask and resolution of radial distance
Resolution of angular is not the only factor that causes the size of image to expand. Due to the inherent characteristic of logarithmic transformation, the pixels near the origin of the coordinate are stretched, and the pixels far away from the origin of the coordinate are compressed (as shown in Fig. 5).

Fig. 5 Schematic diagram of non-uniform sampling of log-polar coordinate transformation
In order to ensure that the information of the image will not be lost after the log-polar coordinate transformation, the sampling rate of the image along radial distance must be set reasonably (suppose the resolution of angular after log-polar coordinate transformation is RA, and the resolution of radial distance is RD). For the interpretation window with a side length of 2n+1, the log-polar coordinate transformation only needs to satisfy the formula (9): log log 1 a a n n   (9) It can be derived from the formula (9): 1 n a n   (10) Then the RD can be expressed as: 1 log n n RD n   (11) Assuming that the side length of the window is 63, It can be seen from Table 1

Establishing relative coordinate system
Before performing log-polar coordinate transformation on the image, it is necessary to determine the new coordinate of each pixel in the image after the coordinate transformation. Since the particle image uses the interpretation window as the matching unit, it means that in order to complete the matching of the entire image, the coordinates of the pixels need to be calculated multiple times, this process involves a lot of calculations (for example, for a particle image with a size of 512*512, if the side length of interpretation window is set as 33 and the interval between the interpretation windows is set as 16, then the calculation of the coordinates needs to be repeated 961 times). Through analysis, it is found that, since the size and shape of the interpretation window remain unchanged during the process of image matching, the center of the interpretation window can be used as the origin to establish a relative coordinate system, and then calculates the relative coordinate of each pixel in the interpretation window, then storing the relative coordinates in order. When the pixels in each interpretation window need to be mapped, the coordinates of pixels can be directly read to complete the log-polar coordinate transformation.

Establishing image pyramid
The TPIM algorithm implements image matching by traversing all pixels in the search area, and setting an interpretation window with each pixel as the center. If the estimated displacement per unit time is x, the relevant calculation needs to be repeated   2 1 x  times to complete the matching of one interpretation window. In order to reduce the amount of calculation, the algorithm establishes the image pyramid through layered sampling of the original image (Fig. 6). Firstly, the algorithm matches the image at the lowest scale, and uses the calculated displacement as a guide to translate the interpretation window in the next layer of the image (the image at the higher scale) to reduce the size of search area and reduce the complexity of algorithm [17,18] .

Results and discussion
The TPIM algorithm is implemented based on the platform of Microsoft Visual Studio C++, and the processing results of the image pairs by the algorithm are displayed in the form of pictures consist with displacement vectors, and the length of the vector is equal to the displacement of interpretation window.
In the experiment, 4 types of particle images were used to test the algorithm. additionally, the test results The images of the first type are standard particle image pair with translational transformation provided by the official website of the International PIV Challenge (the address of website is: www.pivchallenge.org). Fig. 7 shows the matching results: The images of the third type of are particle image pair of Rankine vortex containing two vortices, and the coordinates of the vortices' centers are respectively (100,100) and (200,200) (unit: pixel). Fig. 9 shows the matching results: Fig. 9 Comparison of processing results of particle image pairs with twin vortex of Rankine by different algorithms. The images of the fourth type of are particle mage pair of standard Hammer-Osen vortex, and the coordinate of the vortex' center is (80,135) (unit: pixel). Fig. 10 shows the matching results: Fig. 10 Comparison of processing results of particle image pairs with vortex of Hamel-Oseen by different algorithms. For the four types of particle image pairs involved in Fig. 7-Fig. 10, Table 2 shows the specific information of the images and the matching results of different algorithms (for convenience, I1, I2, I3,   and I4 in the table respectively represents: image pairs with translation transformation provided by the   PIV Challenge website, image pairs with rotation transformation provided by the PIV Challenge website, image pairs with twin vortex of Rankine, image pairs with vortex of Hamel-Oseen): As can be seen from Table 2 (a) (b) Fig. 11 Working principle and main experimental equipment of PIV system. a is the principle of two-dimensional PIV experiment, and b is the Construction of equipment for gas-solid two-phase flow experiment In the experiment for gas-solid two-phase flow, the air compressor is used to blow SiO2 tracer particles into the area to be measured, at the same time, the laser and digital camera are turned on, and the laser emitted by the laser irradiates the tracer particles to form reflected light, and the optical lens of the CCD camera receives the reflected light to form a pair of particle images. In the experiment for gasliquid two-phase flow, the smoke generator is used to atomize the mixture of oil-water and spray it into the transparent glass container, at the same time, turn on the blower to drive the atomized fluid to flow. Similar to the experiment for gas-solid two-phase flow, the tracer particles reflect the laser light emitted by the laser, and the digital camera receives the emitted light to form the particle image pairs. Fig. 12 shows the scene of the experiment. Fig. 12 Scene of the gas-liquid two-phase flow experiment.
The CCD digital camera and double pulse laser are important components in the PIV experimental system, Table 3 and Table 4 respectively show the main parameters of related components in detail.

Conclusion
For the image pairs of turbulence, the matching accuracy of traditional algorithms decreases due to the coupling of rotation and deformation. In this paper, by analyzing the characteristics of particle images, combining the theory of affine transformation, and using log-polar coordinate transformation, an algorithm named TPIM (Turbulent Particle Image Matching) with high robustness to rotation and deformation is proposed. For the high complexity of algorithm introduced by coordinate transformation, by reasonably setting the resolutions of angle and radial distance, establishing the relative coordinate system and image pyramid, the complexity of the algorithm is reduced to about 15% of the original. The test results of various types of image pairs show that, the TPIM algorithm is suitable for matching particle images of turbulence, and the matching accuracy of the algorithm can reach more than 99%.