A parametric study of 3D PTV algorithms based on a two-view collimated imaging model

The camera model is established on the basis of a two-view collimated shadowgraph system [38]. As shown in figure 2, the system consists of two sets of conventional Z-type shadowgraph configurations, sharing an ordinary halogen tungsten lamp. Four identical parabolic mirrors, with 30 cm in diameter and 300 cm in focal length, are used to form two collimated beams which intersect in the test region. Two synchronized high-speed cameras (Photron SA-Z and Photron Mini AX100) are used to record the images.

Zoom In Zoom Out Reset image size

The optical geometry of the two-view collimated shadowgraph system is shown in figure 3. The projection model of each camera includes two parts: one part is the parallel projection in the collimated light path; the other is a pinhole model in the converging region, which corresponds to the converging beams projected onto a charge coupled device (CCD). In general, the pinhole model can be simplified as a center point and an image plane. In the following, subscripts l and r are used to denote elements related to the left and right camera, respectively. Point $({x_w},{y_w},{z_w})$ in the world coordinates is first projected on the parallel projection plane, which then forms a ray by connecting the center point $C$ and is projected on the image plane as $({u_l},{v_l})$ and $({u_r},{v_r})$. As derived in [38], the relationship between the object point $({x_w},{y_w},{z_w})$ in the world coordinate system and the corresponding points in the left and the right image coordinate systems $({u_l},{v_l})$ and $({u_r},{v_r})$ can be described as

$$\begin{align}\left\{ {\begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {{u_l}} \\ {{v_l}} \\ 1 \end{array}} \right] = \underbrace {\left[ {\begin{array}{*{20}{c}} {{\alpha _l}}&{{\gamma _l}}&{{u_{0l}}} \\ 0&{{\beta _l}}&{{v_{0l}}} \\ 0&0&1 \end{array}} \right]}_{{A_l}}\underbrace {\left[ {\begin{array}{*{20}{c}} 1&0&0&0 \\ 0&1&0&0 \\ 0&0&0&1 \end{array}} \right]}_{{E_0}}\left[ {\begin{array}{*{20}{c}} {{x_w}} \\ {{y_w}} \\ {{z_w}} \\ 1 \end{array}} \right]} \\ \!\!\!\!{\left[ {\begin{array}{*{20}{c}} {{u_r}} \\ {{v_r}} \\ 1 \end{array}} \right]\!\! = \!\! \underbrace {\left[ {\begin{array}{*{20}{c}} {{\alpha _r}}&{{\gamma _r}}&{{u_{0r}}} \\ 0&{{\beta _r}}&{{v_{0r}}} \\ 0&0&1 \end{array}} \right]}_{{A_r}}\underbrace {\left[ {\begin{array}{*{20}{c}} {r_{11}^{rl}}&{r_{12}^{rl}}&{r_{13}^{rl}}&{t_x^{rl}} \\ {r_{21}^{rl}}&{r_{22}^{rl}}&{r_{23}^{rl}}&{t_y^{rl}} \\ 0&0&0&1 \end{array}} \right]}_{{E_{rl}}}\!\!\left[ {\begin{array}{*{20}{c}} {{x_w}} \\ {{y_w}} \\ {{z_w}} \\ 1 \end{array}} \right]} \end{array}} \right.\!\!\!\!\!\!.\end{align} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

In equation (1), the intrinsic matrix of each camera (${A_l}$,${A_r}$) and the extrinsic matrix between two cameras (${E_{rl}}$) can be resolved through a camera calibration process.

Once the camera model is established, the synthetic particle images can be generated from 3D particle coordinates using the projection principle described by equation (1). In the current study, the projection is conducted using an idealized camera model, in which the camera calibration error and the particle identification error are not involved. Two sets of ideal camera calibration parameters are used in the current study, which are listed in table 1. The view angles in the two groups of camera parameters are both set at 30°.

Table 1. Set values of the camera parameters.

Camera parameter	Set 1	Set 2
${A_l}$	$\left[ {\begin{array}{*{20}{c}} {100}&0&0 \\ 0&{100}&0 \\ 0&0&1 \end{array}} \right]$	$\left[ {\begin{array}{*{20}{c}} {63.662}&0&{200} \\ 0&{63.662}&{200} \\ 0&0&1 \end{array}} \right]$
${A_r}$	$\left[ {\begin{array}{*{20}{c}} {100}&0&0 \\ 0&{100}&0 \\ 0&0&1 \end{array}} \right]$	$\left[ {\begin{array}{*{20}{c}} {63.662}&0&{200} \\ 0&{63.662}&{200} \\ 0&0&1 \end{array}} \right]$
${E_{rl}}$	$\left[ {\begin{array}{*{20}{c}} 1&0&0&0 \\ 0&{0.5}&{0.866}&0 \\ 0&0&0&1 \end{array}} \right]$	$\left[ {\begin{array}{*{20}{c}} {0.9804}&{ - 0.007}&{0.1986}&{ - 3.1416} \\ { - 0.0048}&{0.9844}&{0.1761}&{ - 3.1416} \\ 0&0&0&1 \end{array}} \right]$

The synthetic particles are generated using the direct numerical simulation (DNS) results of homogeneous isotropic turbulence from the Johns Hopkins University Turbulence Databases [39]. The calculation domain has a size of ${\text{2}}\pi \; \times \;{\text{2}}\pi \; \times \;{\text{2}}\pi$, at a flow condition of ${\operatorname{Re} _\lambda } = 418$. The particles are distributed in the flow field randomly. An example is shown in figure 4, which presents particle trajectories with 2000 particles in 50 frames with a time interval of 0.005 s.

Zoom In Zoom Out Reset image size

Based on the 3D particle positions, the synthetic particle images can be generated using the camera calibration parameters shown in table 1. In the current study, three different algorithms are tested parametrically, including the DDT, 4BE and spatial-temporal four frame best estimate (ST-4BE) algorithms. The synthetic particle images for testing the DDT and 4BE algorithms are generated using Set 1 camera parameters, with an image resolution of 500 × 500 pixels. The images for ST-4BE algorithm employ Set 2 camera parameters, with an image resolution of 400 × 400 pixels. Images with different image particle densities are generated. Two examples are shown in figure 5, with 0.01 ppp and 0.025 ppp respectively. It can be seen that the particle number increases with an increase in image particle density. The dense particle conditions are more challenging for particle identification, stereo pairing and tracking in 3D space.

Zoom In Zoom Out Reset image size

Besides the image particle density, the displacement–spacing ratio $\xi$ is also an important parameter affecting the particle tracking performance. The dimensionless parameter $\xi$ is employed here, which is defined as [9]:

$$\begin{equation}\xi = \frac{{\Delta r}}{{{\Delta _0}}}\end{equation} \tag{ 2 }$$

$$\begin{equation}{\Delta _0} = {\left(\frac{{{N_{{\text{all}}}}}}{V}\right)^{ - 1/3}}.\end{equation} \tag{ 3 }$$

In equations (2) and (3), $\Delta r$ is the average particle displacement between two frames; ${\Delta _0}$ is the average spacing between particles in the same frame, which can be estimated using the total particle number ${N_{{\text{all}}}}$ and the total test volume $V$. The displacement–spacing ratio is related to both the image particle density and the time interval. With the same time interval, the displacement ratio increases at larger image particle densities. With the same image particle density, the displacement–spacing ratio $\xi$ can be changed by varying the time interval between two frames. As shown in figure 6, for 1600 particles in 50 frames, the particle trajectories show dramatic differences at the two displacement ratios 0.0297 and 0.1459. With increasing $\xi$, the lengths of the particle trajectories are increased. Different trajectories tangle together, which makes the tracking process more difficult.

Zoom In Zoom Out Reset image size

The summary of the test conditions for the three algorithms is listed in table 2. As will be shown in sections 3–5, the performances of the three algorithms deteriorate with increasing image particle density and displacement–spacing ratio. In the present study, the test is terminated when the correct ratio is below 60%. As shown in table 2, the test range of the displacement–spacing ratio does not significantly overlap for the three algorithms, which is mainly due to their performance differences.

Table 2. Summary of test conditions.

Algorithm	Image particle density (ppp)	$\xi$	Algorithm	Image particle density (ppp)	$\xi$
DDT	0.005	0.0094		0.005	0.0236
	0.01	0.0118		0.005	0.1158
	0.015	0.0136		0.0063	0.0563
	0.02	0.0149		0.01	0.0297
	0.025	0.0161		0.01	0.1459
				0.01	0.2861
			ST-4BE	0.018	0.0362
4BE	0.005	0.0108		0.025	0.0404
	0.01	0.0136		0.025	0.1980
	0.015	0.0156		0.038	0.0465
	0.02	0.0171		0.05	0.0508
	0.025	0.0185		0.05	0.2495
	0.03	0.0196		0.075	0.0582
	0.035	0.0207		0.088	0.0615
	0.04	0.0216		0.1	0.0641

For two-view systems, the particles at each time step are recorded from two different views. The 3D particle coordinates can be obtained using the imaging coordinates and camera calibration parameters. Figure 7 shows a two-view imaging example with parallel beams, while the following description is also applicable to the pinhole model. For each particle, there are two projecting lines passing through it. The intersections of the projecting lines are called candidate points; the whole set is called a candidate grid. It is noted that the candidate grid may contain additional points without real particles. It is also possible for multiple particles to lie on a projecting line. For n particles, there are at most n projecting lines for each projection direction, which corresponds to a maximum of n² grid points. With the camera calibration parameters, the candidate grid can be easily resolved from the projecting lines.

Zoom In Zoom Out Reset image size

Assume that at a specific time step t, the candidate grid ${G^{(t)}}$ contains $l(t)$ candidate points. At each candidate point $g_i^{(t)}$, a discretized variable $\xi _i^{(t)}$ is designated. If there is a particle then $\xi _i^{(i)} = 1$; otherwise $\xi _i^{(i)} = 0$. The 3D reconstruction process is transformed to resolve the solution ${\underline x ^{(t)}}: = {(\underline \xi _1^{(t)},\ldots,\underline \xi _{l(t)}^{(t)})^T} \in {\{ 0,1\} ^{l(t)}}$ constrained by the linear inequalities:

$$\begin{equation}\begin{gathered} {A^{(t)}}{x^{(t)}} \geqslant {b^{(t)}},{x^{(t)}} \in {\{ 0,1\} ^{l(t)}} \\ {A^{(t)}} \in {\{ 0,1\} ^{k(t) \times l(t)}} \\ \end{gathered} \end{equation} \tag{ 4 }$$

where ${b^{(t)}}: = {(1,\ldots,1)^{\text{T}}} \in {\left\{ 1 \right\}^{k(t)}}$ represents the imaging data from two views; $k(t)$ denotes the total number of detected particles; and ${A^{(t)}}$ represents the particle imaging effect matrix. If the candidate point $g_i^{(t)}$ has contributions to the particle image $b_j^{(t)}$, $A_{ji}^{(t)} = 1$; otherwise, $A_{ji}^{(t)} = 0$. Here, we use ${x^{(t)}}$ and ${\underline x ^{(t)}}$ to denote the variable and specific solutions, respectively.

For 3D particle tracking, it does not only require the 3D reconstruction solution based on equation (4) to be resolved, but also for the particle position ${x^{(t - 1)}}$ at time step $t - 1$ to correspond to the following position ${x^{(t)}}$ at time t. Here a vector ${w^{(t)}} = {(\omega _1^{(t - 1,t)},\ldots,\omega _{l(t)}^{(t - 1,t)})^{\text{T}}}$ is defined to denote the weights associated with each grid point. Then the tracking from time step $t - 1$ to t can be described by the following discrete optimization problem:

$$\begin{equation}\begin{gathered} {\text{minimize }}{w^{(t)}} \bullet {x^{(t)}} \\ {\text{subject to }}{A^{(t)}}{x^{(t)}} \geqslant {b^{(t)}},{x^{(t)}} \in {\{ 0,1\} ^{l(t)}} \\ \end{gathered}.\end{equation} \tag{ 5 }$$

In equation (5), ${w^{(t)}} \bullet {x^{(t)}}$ represents the scalar product between weight parameter vector ${w^{(t)}}$ and solution vector ${x^{(t)}}$.

Here the weight parameter vector ${w^{(t)}}$ is defined as

$$\begin{equation}\omega _i^{(t - 1,t)}: = \mathop {\min }\limits_{j:\xi _j^{(t - 1)} = 1} \left\{ {{\text{dist}}\left(g_i^{(t)},g_j^{(t - 1)}\right)} \right\},\end{equation} \tag{ 6 }$$

where ${\text{dist}}(g_i^{(t)},g_j^{(t - 1)})$ denotes the distance of a certain particle from the previous time step t – 1 to the next frame at time t. The initialization can be set as ${w^{(0)}}: = 0$, if there is no prior information available. When the framing rate reaches the kHz or MHz order, the distance between neighboring frames can be considered as 'short' for most test conditions. Then the distance can be expressed as the Euclidian distance

$$\begin{equation}{\text{dist}}\left(g_i^{(t)},g_j^{(t - 1)}\right): = {\left\| {g_i^{(t)} - g_j^{(t - 1)}} \right\|_2}.\end{equation} \tag{ 7 }$$

It should be noted that the solutions for a finite number of projections have the possibility of non-uniqueness. The ambiguities can be improved by introducing extra constraints, while the particle position in the space domain and the maximum velocity value are employed in the current study.

In the present study, four parameters are defined to quantitatively evaluate the algorithm performance, which are listed below:

$$\begin{equation}\begin{array}{*{20}{l}} {{F_t} = \frac{{{N_t}}}{{{N_{{\text{all}}}}}},{F_{t({\text{trac}})}} = \frac{{{N_{t({\text{trac}})}}}}{{{N_{{\text{all}}}}}}} \\ [10pt] {{F_{\text{g}}} = \frac{{{N_{\text{g}}}}}{{{N_{{\text{all}}}}}},{F_{t({\text{trac}})}} = \frac{{{N_{{\text{g}}({\text{trac}})}}}}{{{N_{{\text{all}}}}}}.} \end{array}\end{equation} \tag{ 8 }$$

${N_t}$ represents the particle number being correctly reconstructed at every frame. ${N_{\text{g}}}$ is the number of ghost particles at each frame. ${N_{t({\text{trac}})}}$ is the correctly tracked trajectories at each frame, while ${N_{{\text{g}}({\text{trac}})}}$ represents the wrong ones. ${N_{{\text{all}}}}$ is the total number of correct particles. From the definition in equation (8), ${F_t}$ and ${F_{\text{g}}}$ are used to evaluate the 3D reconstruction correctness, while ${F_{t({\text{trac}})}}$ and ${F_{{\text{g}}({\text{trac}})}}$ are employed for 3D tracking.

3.2.1. Tracking results at different frames.

The DDT algorithm is first tested using 1000 particles, with an image particle density of 0.005 ppp and a displacement–spacing ratio of 0.0094. The 3D reconstruction and tracking results in 100 frames are shown in figure 8. The curves indicate that the variations of ${F_{t({\text{trac}})}}$ and ${F_{{\text{g}}({\text{trac}})}}$ with frame number are in accordance with ${F_t}$ and ${F_{\text{g}}}$ respectively. The correctness ${F_t}$ decreases from 100% in the first frame to 83% in the tenth frame. Then it increases with the frame number, reaching 99% in the 40th frame and remaining there until the final frame. Although the first frame has the highest correctness, it also involves the most ghost particles, as high as 158%. The number of ghost particles decreases dramatically with the frame number, which is reduced to only 1% in the 40th frame. This variation is mainly due to the lack of prior information in the first frame, where it is assumed that all the candidate points have particles. With the processing of the DDT program, the ghost particles are excluded gradually, which also includes some real particles. After 40 frames, the constraint is more reliable, which helps to track the lost particles again.

Zoom In Zoom Out Reset image size

Considering that the poor performance in the first 40 frames is mainly due to the lack of constraints from the time domain, an iteration strategy is introduced to improve the DDT algorithm. The iteration is conducted from forth-to-back for odd frame numbers using time steps $t - 1$ and t, and back-to-forth for even frame numbers using $t + 1$ and t. The results after two iteration operations are shown in figure 9. It can be seen that the correctness of the DDT algorithm is greatly improved by the iteration. The correctness of ${F_t}$ and ${F_{t({\text{trac}})}}$ is higher than 99% from the first to the end frame, where the ratio of ghost particles is kept at less than 1%.

Zoom In Zoom Out Reset image size

3.2.2. Tracking results at different image particle densities.

Based on the iteration DDT algorithm, the performances at different image particle densities from 0.005 ppp to 0.025 ppp are tested, with a constant time interval of 0.005 s. The 3D reconstruction and tracking results are shown in figure 10. It can be seen that the variations of ${F_{t({\text{trac}})}}$ and ${F_{{\text{g}}({\text{trac}})}}$ are very similar with ${F_t}$ and ${F_{\text{g}}}$. The correctness is over 99% at 0.005 ppp, and decreases with increasing image particle density. At 0.01 ppp, the correctness is around 94%, while it drops below 90% from 0.015 ppp. At 0.025 ppp, the correctness is only about 60%. The number of ghost particles increases with the image particle density, in a mirroring curve with respect to the correctness. The variations indicate a strong relevance between ghost particles and wrong detections of candidate points. With a dramatic increase in ghost particles, the correctness decreases accordingly. The quantitative analysis indicates that it is a great challenge for the global DDT algorithm to deal with dense particles.

Zoom In Zoom Out Reset image size

3.2.3. Tracking results at different displacement–spacing ratios.

The 3D reconstruction and tracking results under different displacement–spacing ratios are shown in figure 11. The displacement–spacing ratio is in a range from 0.0094 to 0.0161. It can be seen that the correctness of ${F_t}$ and ${F_{t({\text{trac}})}}$ is around 99% when $\xi$ is 0.0094. The correctness can stay as high as 92% when $\xi$ is below 0.012. With further increasing of the displacement ratio, the correctness decreases dramatically, dropping to 60% at $\xi = 0.016$. The ghost particles and ghost velocity vectors increase with $\xi$ accordingly. The results indicate that the algorithm performance is very sensitive to the displacement–spacing ratio, which should be kept small enough for practical applications.

Zoom In Zoom Out Reset image size

3.2.4. Tracking performance summary for DDT algorithm.

As aforementioned, the multi-view system can help to restrict the candidate points to a very small region. The 3D coordinates of the particles are reconstructed first, which are then tracked in time sequences. Since the candidate points are quite limited, it is easy to exclude the ghost particles using certain constraints. In two-view systems, the candidate points can only be confined to those near an epipolar line, which makes it quite difficult to exclude the large number of ghost particles in the 3D reconstruction process. The strategy of the 4BE algorithm in the present study is shown in figure 12. The 3D reconstruction is based on a searching region near the epipolar line, where a small number of ghost particles are excluded if they are out of the measurement volume. Then most of the ghost particles are identified and removed in the tracking process in space.

Zoom In Zoom Out Reset image size

In the current 4BE algorithm, the 3D coordinates in four frames, $n,n + 1,n + 2,n + 3$, are used to track the particles. Firstly, it needs to initialize the track by matching the particle $x_i^n$ to its position in next frame $x_i^{n + 1}$. In frame n+ 1, a searching region is defined to exclude ghost particles, which is set as 1.1 times the maximum particle displacement. All the particles in the searching region are considered as candidate points. The consideration is to involve all the possible points but with affordable calculation cost.

After the initialization, the two positions at frames n and n + 1 can form a track. Then the particle position $\widetilde x_i^{\,n + 2}$ at frame n + 2 can be predicted using the following equations:

$$\begin{equation}\widetilde x^{\,n + 2}_i = x_i^{\,n + 1} + \widetilde v_i^{\,n + 1}\Delta t\end{equation} \tag{ 9 }$$

$$\begin{equation}\widetilde v_i^{\,n + {\text{1}}} = (x_i^{\,n + 1} - x_i^n)/\Delta t.\end{equation} \tag{ 10 }$$

In equations (9) and (10), $\Delta t$ is the time interval between two frames. $\widetilde v_i^{\,n + 1}$ is the particle velocity, which can be calculated using equation (10), based on a constant velocity assumption.

Then the searching area for frame n + 2 is set based on the predicted positions. The searching region is set at 1.1 times the maximum possible displacement to exclude the ghost particles. If there is no particle detected in the searching region, then this track will be discarded. If there is one or more particles detected, all the tracks are kept for the following procedures. The predictions at frame n + 3 are based on the previous three frames n, n + 1 and n + 2, using both the constant velocity and acceleration assumptions, as indicated by equations (11)–(13).

$$\begin{equation}\widetilde x_i^{\,n + 3} = x_i^{n + 2} + \widetilde v_i^{\,n + 2}\Delta t\end{equation} \tag{ 11 }$$

$$\begin{equation}\widetilde v_i^{\,n + 2} = v_i^{n + 1} + \widetilde a_i^{\,n + 1}\frac{{\Delta t}}{2}\end{equation} \tag{ 12 }$$

$$\begin{equation}\widetilde a_i^{\,n + 1} = (x_i^{n + 2} - 2x_i^{n + 1} + x_i^n)/{(\Delta t)^2}.\end{equation} \tag{ 13 }$$

In equation (13), $\widetilde a_i^{\,n + 1}$ is the acceleration at frame n + 1. After the prediction, a searching region for frame n + 3 is selected, where the candidate points are tracked. If there are no particles detected, this track is discarded. If one or more particles are found, then multi-tracks remain.

Up to now, for any particles in frame n, there might be multiple candidate tracks in the following three frames. In order to resolve the most possible track, a cost function $\phi _{ij}^n$ is introduced to optimize the tracking process, which is defined as

$$\begin{equation}\phi _{ij}^n = \left\| {x_j^{n + 3} - \widetilde x_i^{\,n + 3}} \right\|\end{equation} \tag{ 14 }$$

In equations (11) and (14), $x_j^{n + 3}$ is the detected particle position at frame n + 3; $\widetilde x_i^{\,n + 3}$ is the predicted particle position. With the cost function optimization, the minimum distance between the detected and predicted particles is selected as the best track. In the current algorithm, only the trajectories longer than ten time steps are kept, which helps to exclude the large number of ghost particles in two-view systems.

4.2.1. Tracking results at different frames.

The ratios of correctly reconstructed and tracked particles at different frames are presented in figure 13. The image particle density is 0.005 ppp, with a displacement–spacing ratio $\xi$ of 0.0108. It can be seen that at relatively small image particle densities, the correctness of both ${F_t}$ and ${F_{t({\text{trac}})}}$ is kept very close to 100%. The ratio of ghost particles ${F_{\text{g}}}$ is about 7.5% in the first frame, which then decreases gradually and reaches 3.8% in the 49th frame. Due to the lack of sufficient prior knowledge at the very beginning, some ghost particles cannot be excluded. With the progression of particle tracking, the constraint from the time domain begins to play a role, which helps to remove more ghost particles. As can be seen from figure 13(b), the trends of ${F_{{\text{g}}({\text{trac}})}}$ are very similar to ${F_{\text{g}}}$, with differences of less than 1.5%.

Zoom In Zoom Out Reset image size

4.2.2. Tracking results at different image particle densities.

The 4BE algorithm is tested with different image particle densities in the range of 0.005 ppp–0.04 ppp, with a constant time interval of 0.005 s. The results are shown in figure 14, where the DDT test results are also presented for comparison. It can be seen that the ratio of ${F_t}$ is kept at more than 99% up to 0.04 ppp. However, the ghost particles increase dramatically with an increase in image particle density, which rises from 4.5% to 35.0% in the test range. The variations of ${F_{t({\text{trac}})}}$ and ${F_{{\text{g}}({\text{trac}})}}$ are very similar to those of ${F_t}$ and ${F_{\text{g}}}$. The comparison with the DDT algorithm indicates that the 4BE algorithm has better performance when the image particle density is larger than 0.01 ppp.

Zoom In Zoom Out Reset image size

4.2.3. Tracking results at different displacement–spacing ratios.

The 4BE algorithm is tested under different displacement–spacing ratios from 0.0108 to 0.0216. The results are presented in figure 15. It can be seen that both ${F_t}$ and ${F_{t({\text{trac}})}}$ remain higher than 98% for all the test cases. However, the ratio of ghost particles increases with the displacement–spacing ratio, which increases from to 7% to 36% in the test range. Compared with the DDT algorithm, the 4BE algorithm has much higher correctness and fewer ghost particles when the displacement–spacing ratio is larger than 0.012.

Zoom In Zoom Out Reset image size

4.2.4. Tracking performance summary for 4BE algorithm.

The strategy of the ST-4BE algorithm is shown in figure 16, in which the reconstruction and tracking processes are strongly coupled. The cost function involves the constraints from both the space and time domains. The constraints from epipolar geometry, predictions in 2D images and 3D space, and reprojection from space to image coordinates are utilized in a comprehensive way. The algorithm is explained in detail below.

Zoom In Zoom Out Reset image size

For an image at time n, the epipolar constraint is first applied to find all the possible stereo pairing candidates $x_1^n$ and $x_2^n$ in two images, based on which the 3D coordinates ${x^n}$ are reconstructed using the calibration parameters. For the particle image in the next frame n + 1, specific searching regions are selected based on the positions of $x_1^n$ and $x_2^n$, where the stereo pairing candidates $x_1^{n + 1}$ and $x_2^{n + 1}$ are found by epipolar constraints. Then the 3D coordinates ${x^{n + 1}}$ at frame n + 1 can be obtained. The size of the searching region is set at 1.1 times the maximum possible displacement.

After initialization, the 3D coordinates ${\widetilde x^{\,n + 2}}$ at frame n + 2 can be predicted using equations (9) and (10), which are re-projected back to two-view images with coordinates $\widetilde x_1^{\,n + 2}$ and $\widetilde x_2^{\,n + 2}$. Then the searching region is selected based on $\widetilde x_1^{\,n + 2}$ and $\widetilde x_2^{\,n + 2}$, where the stereo pairing is conducted to tackle corresponding particle coordinates $x_1^{n + 2}$ and $x_2^{n + 2}$. The 3D coordinates at frame n + 2 can then be obtained. This procedure is repeated to resolve the 3D coordinates at frame n + 3.

Up to now, all the possible tracks are established. Then a cost function is established to resolve the best particle trajectories. Several cost functions have been attempted in [35], among which the following one is selected,

$$\begin{equation}{\phi ^n} = \sum\limits_{i = 1,2} {\left( {{{\left\| {x_i^{\,n + 3} - x{^{^{\prime}}}_i^{n + 3}} \right\|}^2} + {{\left\| {x_i^{n + 3} - \widetilde x_i^{\,n + 3}} \right\|}^2}} \right)}.\end{equation} \tag{ 15 }$$

In equation (15), $i = 1,2$ represents the two views; $x_i^{n + 3}$ is the re-projected 2D image coordinates from reconstructed particles at frame $n + 3$; $\widetilde x_i^{\,n + 3}$ is the predicted 2D image coordinates at frame $n + 3$. Only the trajectories with minimum cost function value are remained; the others are discarded.

5.2.1. Tracking results at different frames.

The 3D particle reconstruction and tracking results of ST-4BE are shown in figure 17, in which the results of the 4BE algorithm are also presented for comparison. The image particle density is 0.005 ppp, with a displacement–spacing ratio $\xi$ of 0.0236. It can be seen that the ${F_t}$ and ${F_{t({\text{trac}})}}$ of the ST-4BE algorithm keep a high correctness over 99% for all the frames. The two lines coincide with the 4BE results. However, the two algorithms resolve dramatically different ratios of ghost particles. As shown in figure 17, the resolved ghost particle ratio by the ST-4BE algorithm is always less than 1%, which is much less than that of the 4BE algorithm. The comparison indicates that the comprehensive application of multi-constraints in both time and space of the ST-4BE algorithm excludes the ghost particles effectively.

Zoom In Zoom Out Reset image size

5.2.2. Tracking results at different image particle densities.

The performance of ST-4BE is investigated at different image particle densities, ranging from 0.005 ppp to 0.1 ppp, with a constant time interval of 0.005 s. As shown in figure 18, the results are compared with the other two algorithms DDT and 4BE. The comparison shows that the ST-4BE algorithm makes significant improvements. The correctness of ${F_t}$ and ${F_{t({\text{trac}})}}$ remains higher than 99% until 0.1 ppp, while the ghost particle ratio of ${F_{\text{g}}}$ and ${F_{{\text{g}}({\text{trac}})}}$ is less than 1% in the whole test range.

Zoom In Zoom Out Reset image size

5.2.3. Tracking results at different displacement–spacing ratios.

The ST-4BE algorithm is tested at different displacement–spacing ratios, which vary from 0.0236 to 0.2861, by varying both the particle image density and time interval. It can be seen from figure 19 that the correctness is higher than 92.5% at $\xi = 0.2$, and decreases dramatically with a further increase in the displacement–spacing ratio. The correctness drops below 58% at $\xi = 0.2861$. The correctness of ST-4BE and 4BE is similar when $\xi$ is smaller than 0.03, which is much higher than the DDT algorithm. However, the ST-4BE shows far fewer ghost particles than the 4BE algorithm, remaining below 5% in the wide test range.

Zoom In Zoom Out Reset image size

5.2.4. Tracking performance summary for ST-4BE algorithm.