# Consensus-based tracking for 3D PTV at high seeding densities

Jean Le Bris, Benjamin Leclaire, Philippe Cornic, Frédéric Champagnat, Benjamin Musci, Adam Cheminet

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jean Le Bris, Benjamin Leclaire, Philippe Cornic, Frédéric Champagnat, Benjamin Musci, et al. Consensus-based tracking for 3D PTV at high seeding densities. Experiments in Fluids, 2025, 66, pp.186. 10.1007/s00348-025-04109-7. hal-05264637

HAL Id: hal-05264637

https://hal.science/hal-05264637v1

Submitted on 17 Sep 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Consensus based tracking for 3D-PTV at high seeding densities

Jean Le Bris¹, Benjamin Leclaire², Philippe Cornic³, Frédéric Champagnat³, Benjamin Musci¹, and
Adam Cheminet¹

¹SPEC, CEA, CNRS, Université Paris-Saclay, CEA Saclay, F-91191 Gif-sur-Yvette, France

²DAAA, ONERA, Institut Polytechnique de Paris, 92190, Meudon, France

³DTIS, ONERA, Université Paris-Saclay, 91120, Palaiseau, France

Corresponding Author: Jean Le Bris

E-mail: jean.lebris@cea.fr

Abstract

A robust pairing algorithm with outlier removal is introduced
in the context of two-pulse 3D Particle Tracking Velocimetry
at high seeding densities, with high concentrations of ghost
particles. Integrating the Vector Field Consensus approach
from Ma et al. (2014), the algorithm, its underlying hypothe-
ses and its relevant input parameters are investigated in the
context of turbulent flow measurements. 2D synthetic tests
are first carried out to quantify the algorithm's performance
and derive simple guidelines for optimal parameter tuning
strategies based on experimental quantities. It is found that
2D vector fields with up to 90% outliers can be handled by
our algorithm. 3D synthetic tests are then implemented to
test the tracking strategy robustness to increasing image
densities and ghost particle concentrations. We show that
our algorithm can be used for particle pairing in particle
clouds with up to 50% of ghost particles. Results submitted
on the two-pulse dataset of the First LPT challenge, using
the associated data portal with automatic evaluation, also
showcase the overall excellent performances of the method.
Finally, the method is used successfully on experimental
data from our Giant Von Kármán setup (characterized by
up to 65% of ghost particles), as evidenced by comparisons
of its output with respect to results provided by the Shake-
The-Box algorithm, and with results provided by a pairing
approach using a 3D cross-correlation predictor.

# 1 Introduction

Turbulent flows at high Reynolds numbers can be character-
ized by a very large range of spatial scales, so that measuring
their velocity field, from object size down to their smallest
scales, is a tremendous challenge. In recent years, particle-
based volumetric velocimetry approaches have experienced
dramatic performance gains, allowing them to multiply the
maximum reachable seeding density by several orders of
magnitude. While Tomographic Particle Image Velocimetry
(TPIV) has been the first major step towards this increase in
resolution (Scarano, 2012), its associated spatial filtering due
to cross-correlation has then encouraged to revisit Particle

Tracking Velocimetry (PTV) approaches, making them more
robust to so-called *ghost particles*. As the latter are the result
of the non-uniqueness of 3D particle reconstruction from a
set of (generally four) images, they stem from epipolar geom-
etry only and do not exhibit temporal coherence, contrary
to real particles. This observation has in particular led to
the development of Lagrangian Particle Tracking, whereby
the reconstruction of particle tracks from a time-resolved
acquisition successfully eliminates a large number of ghosts
up to high seeding densities (see, e.g. Schanz et al., 2016).

In the more difficult context of two-pulse acquisition, and
despite significant recent advances (e.g. Cornic et al., 2020;
Novara et al., 2023), the maximum reachable seeding den-
sity, and thus the final spatial resolution, are logically lower.
Modern PTV algorithms usually involve a step of particle
reconstruction from the images (performed separately at the
two measurement instants), and a step of matching or track-
ing, where particles of the two instants are paired to obtain
the velocity fields. Ghost reduction or even elimination up
to very high seeding densities, necessary for turbulent flow
characterization, are still to be obtained during both these
steps, and not only the reconstruction step. Indeed, advanced
state-of-the-art particle reconstruction appraoches such as
the Iterative Particle Reconstruction (IPR) (Wieneke, 2012;
Jahn et al., 2021) still exhibit a performance drop due to
increase of ghost with density: for example, on noise-free
synthetic images, Jahn et al. (2021) reconstruct close to all
true particles at 0.14 ppp after 16 IPR iterations, with 5%
of ghost particles, while at 0.16 ppp, 80% of the true parti-
cles are reconstructed after 20 IPR iterations, with 28% of
ghosts. Very recently, improvements brought by Peak-CNN
(Godbersen et al., 2024) allowed for further improvements on
the rate of ghost reconstructed by IPR: Peak-CNN allowed
to reduce the ghost rate found by IPR to nil at 0.160 ppp
in noisy synthetic datasets. However, in real experimental
conditions, the rate of ghost particles can never be brought
to nil. Mie scattering, particles seen by only a subset of
all cameras or poor camera calibration are all factors that
can multiply the amount of ghost particles reconstructed
(see, e.g. Cheminet et al., 2013). In such conditions, in order
to keep the rate of found particles at a high level, a non
negligible fraction of ghost particles have to be kept. The task of filtering as many of these outliers as possible is then devoted to the pairing step, which is the topic of the present article.

To find the displacement field from both particle clouds, many tracking algorithms introduced in the literature have used the idea of spatial consistency and local flow regularity. Two-Pulse Shake-the-Box (2P-STB) for instance uses Particle Space Correlation introduced in Novara et al. (2016, 2023) which estimates a predictor of the flow via cross-correlation, to then obtain particle pairs by matching within a search radius around the predictor. PSC is a particle-wise variant of cross-correlation, in essence similar to other cross-correlation methods used in 3D-PIV (Discetti and Astarita, 2012; Cheminet et al., 2014) and used also for pairing purposes in Agüera et al. (2016) or Cornic et al. (2020). Other methods try using different approaches that still make use of spatial consistency. To mention but a few, we refer the reader to several founding articles such as Baek and Lee (1996), Ishikawa et al. (2000), Pereira et al. (2006), Mikheev and Zubtsov (2008), Brevis et al. (2011), Jia et al. (2013), Lei et al. (2012), Fuchs et al. (2017). In the latter, the goal is to establish the most probable displacement by considering the maximum of the histogram of all possible displacement in a displacement range, both for the particle of interest and its three closest neighbors. The eventual displacement of the particle of interest is then the closest available displacement to the established most probable displacement. More recently, Benkovic et al. (2020) proposed to integrate a relaxation algorithm based on pairing probabilities to the vision-based particle tracking algorithm imagined by Lei et al. (2012). The trend is nowadays to use machine learning and physics-informed neural network methods to tackle this issue, as in Liang et al. (2022) or Zhou et al. (2023). However, to the best of our knowledge, most of these methods were neither tested nor developed to face high noise levels, ie very high amount of ghost particles, in their data set. Recently, Mercier et al. (2024) proposed a method for robust particle pairing in noisy particle clouds. The method is based on the Coherent Point Drift (CPD) algorithm from Myronenko and Song (2010). CPD finds the optimal transformation that maps a source set of points (typically the particle cloud at $t = 0$) to an objective (particle cloud at $t = 1$) set. This allows to establish pairings between particle clouds, while discarding outlier samples from both time steps. CPD reaches satisfying performances from particle clouds with 16% of ghost particles at 0.05 ppp. At higher proportions of ghost particles, CPD was found to be too selective, discarding many correct vectors from the flow, which can be detrimental in further time-resolved applications. This however demonstrates the gains that can be obtained from algorithms developed for computer vision, as, to the best of our knowledge, there exists no proposed pairing method robust to higher proportions of ghost particles.

In this paper, we propose a new method, based on adapting and building on the Vector Field Consensus (VFC) approach from Ma et al. (2014), which was originally introduced in the computer vision domain for outlier detection and removal of noisy vector fields. VFC first exploits the fact that, as depicted in the right part of Fig. 1, while outliers hold no

physical sense and follow a random uniform distribution over the search radius, the true displacement vectors follow a smoother distribution, closer to a Gaussian around the local spatial mean. An Expectation-Maximization (E-M) algorithm is used to assess which samples obey which distributions, while computing the parameters of those two probability distributions. To do so, an estimate of the coarse-grain displacement is necessary: at each step, a regularized vector field is computed by Tikhonov Regularization, weighted by the probability found for each sample. By evaluating the distance between the current regularized vector field and the samples, the probability of being an inlier for each sample is re-estimated, which in turn will allow a better estimation of the regularized vector field. As we will also observe, the VFC approach is robust to a very high rate of outliers (up to 90% as reported by Ma et al., 2014), so that it can be considered as a sufficient option for ghost removal when integrated in a PTV algorithm. As depicted in the left part of Fig. 1, when dealing with particle clouds with a high amount of ghost particles, we choose to initially pair each particle with every possible particle in a given search radius ie neighborhood. This results in a large set of pairing vectors, containing a large amount of outliers. We will show that, by performing some adaptations and proposing several additions to VFC, together with guidelines for selecting optimal parameters based on simple flow physics considerations, the resulting algorithm, which we term "Flow VFC" (F-VFC), can be used to pair particles very efficiently in noisy 3D particle clouds typically obtained from particle reconstruction algorithms used in PTV applications, and thus has a strong potential, both in the context of 3D and of 4D PTV, as the initialization step of the latter requires vector matching.

The structure of the article is as follows. We first explain in section 2 the principles of the F-VFC approach, encapsulating and adapting the VFC algorithm for Particle Tracking purposes. Tests are then performed on 2D synthetic datasets in section 3, allowing to derive the optimal scaling and tuning of every operating parameter based on simple physical considerations, and yielding first performance assessments. Then, the algorithm's robustness is assessed on 3D synthetic cases in section 4. Finally, we assess its performances on actual experimental data acquired from the Giant Von Kármán (GVK) platform at CEA SPEC in section 5, by comparing F-VFC's results with results obtained by a matching approach using a 3D cross-correlation predictor, considering results obtained on full time sequences by the Shake-The-Box algorithm as ground truth.

## 2 The F-VFC algorithm

This section describes the proposed algorithm, which we name "Flow Vector Field Consensus" (F-VFC). Indeed, the VFC approach remains one of its major building blocks, while we also introduce modifications and additional steps, targeted towards an optimal performance in the context of 3D flow measurements. After introducing in subsection 2.1 the main notations and how the initial vector sets are obtained from particle reconstructions at the two instants
<figure><img src="image_3.png" /><figcaption>Fig. 1 Left: construction of initial pairing vectors between two particle clouds $\mathcal{X}$ (★) and $\mathcal{Y}$ (▲). We write $\Omega$ the volume within the red-dotted circle. With a large enough search radius, the correct displacements are always in the vector set, but so are a large amount of outliers. Right: obtaining the probability density functions of inliers and outliers, and their link with the regularized vector field $\bar{u}$. We write $\gamma$ the proportion of correct vectors in the set. At the scale of a few particles, outliers will exhibit a random behavior, while inliers can be expected to behave close to a Gaussian distribution of standard deviation $\sigma$. At a larger scale, these Gaussian distributions must follow the spatial structure of $\bar{u}$.</figcaption></figure>

(those reconstructions being supposed to be the input of the method), we recall in subsection 2.2 the main features of the VFC algorithm introduced by Ma et al. (2014), in particular to understand its working principle and highlight its main tuning parameters. Subsections 2.3 to 2.5 then expose the adaptations and additions that we introduce in order to address specificities of the PTV context.

## 2.1 Notations and initial vector set

We define

* $D \in \{2; 3\}$ the spatial dimension of the problem

* $\mathcal{X} = \{x_m \in \mathbb{R}^D, m = 1 \to M\}$ the set of particle positions at the first time step. $M$ is the total number of particles reconstructed at this step.

* $\mathcal{Y} = \{y_k \in \mathbb{R}^D, k = 1 \to K\}$ the set of particle positions at the second time step. $K$ is the total number of particles reconstructed at this step.

To find the pairing vectors between these particle clouds $\mathcal{X}$ and $\mathcal{Y}$ (which both include a large number of ghosts), without any prior information, we can pair every $x_m \in \mathcal{X}$ to every single $y_k \in \mathcal{Y}$ located within a search volume $\Omega$ around $x_m$ (see Fig. 1 left). We now denote by $\mathcal{V}$ the set of these vectors and have $\{(x_n, v_n), x_n \in \mathcal{X}, v_n \in \mathcal{V}, n = 1 \to N\}$ the set of particle positions and displacements, of which only an unknown fraction, denoted $\gamma$, is real.

With a large enough search radius, we know that all the correct vectors, or inliers, are present in the dataset. These inliers are coherent with one another: we can expect that their deviation from the mean local displacement in the underlying flow, hereafter denoted by $\bar{u}$, will be small. Following Ma et al. (2014), we assume that inliers follow a gaussian distribution of parameter $\sigma$ around the mean local

displacement $\bar{u}$ (see Fig. 1 right); more precisely, and as will be detailed in the sequel, we choose to set VFC in a way that $\bar{u}$ and $\sigma$ have to identify respectively as the local mean flow and fluctuation amplitude. Wrong vectors, or outliers, on the other hand, will have a more chaotic behavior, with random orientations and magnitudes. Following Ma et al. (2014) again, one can assume that the outliers follow a uniform distribution of $\frac{1}{\Omega}$ (see Fig. 1 right).

To then identify inliers and filter outliers from the dataset, one needs to estimate this mean local displacement $\bar{u}$, then find what samples fit best a gaussian distribution around $\bar{u}$, to be considered as inliers. Estimating $\bar{u}$ is however very challenging due to the outliers. As detailed in the next subsection, the VFC algorithm was introduced as a means to tackle this issue, by iteratively finding the probability that each sample is an inlier/outlier and estimating $\bar{u}$ with this knowledge, i.e. by attributing a low weight in the estimation of $\bar{u}$ to samples with a low probability of being an inlier.

## 2.2 The VFC Algorithm

The VFC algorithm is detailed extensively in Ma et al. (2014), with explanations on the Expectation-Maximization (E-M) algorithm in Bishop (2006) or in the Appendix; we summarize the main underlying hypotheses and features here.

In VFC, the mean local displacement field $\bar{u}$ is represented on a basis of kernel functions $\Gamma$, with one component per vector present in the set, so that

$$ \bar{u}(x) = \sum_{n=1}^{N} \Gamma(x, x_n)c_n, c_n \in \mathcal{V} \quad (1) $$

Choice of a specific kernel basis should be adapted to the expected behavior of $\bar{u}$, and has a direct impact on computational performance. In practice, as will be detailed in
subsection 2.3, we will choose Gaussian functions. VFC being a probabilistic algorithm, every variable that has to be estimated should be given an a priori probability distribution (if it cannot be derived from that of other parameters), either fully known or parameterized by quantities to be estimated. Similar to Ma et al. (2014), we will consider for $\bar{u}$ the so-called “slow-and-smooth” model, assuming that it is a realization of the probability distribution

$$p(\bar{u}) \propto e^{-\frac{\lambda}{2}\phi(\bar{u})}$$

where $\phi$ is a smoothing filter and $\lambda$ a positive real number. As will be shown below, $\phi$ will boil down to a norm on coefficients $c_n$, while $\lambda$ will identify to a regularization parameter. More importantly, the remaining distributions to be prescribed pertain to vectors themselves, for:

* good vectors or inliers, chosen as a Gaussian, of mean $\bar{u}$ and standard deviation $\sigma$

* wrong vectors or outliers, chosen as uniform of value $\frac{1}{\Omega}$, where $\Omega$ is the volume of the search radius used to build the initial vector set

* the proportion $\gamma$ of inliers in the displacement set, chosen as uniform over [0, 1].

Given these distributions, VFC's objective is to iteratively and jointly estimate optimal values of unknown parameters $\Theta = (\bar{u}, \sigma, \gamma)$. Physically, this amounts to determine the probability that each vector is an inlier and outlier, together with local mean flow displacement $\bar{u}$ and associated fluctuation amplitude $\sigma$.

Let now $z_n$ be the latent variable:

$$z_n = \begin{cases} 1 & \text{if } (x_n, v_n) \text{ is an inlier} \\ 0 & \text{if } (x_n, v_n) \text{ is an outlier} \end{cases} \quad (2)$$

Seeking the optimal set $\Theta = (\bar{u}, \sigma, \gamma)$ is equivalent to minimizing the energy

$$E(\Theta) = -\ln p(\bar{u}) - \sum_{n=1}^{N} \ln \sum_{z_n} p(v_n, z_n | x_n, \Theta) \quad (3)$$

The minimization of (3) is extremely complicated, as it requires taking derivatives over unknown values. Ma et al. (2014) propose to use an EM (Expectation-Maximization) approach, where, given a current estimate $\Theta^k$ (at iteration number $k$) of the parameter $\Theta$, the complete-data log posterior $Q(\Theta, \Theta^k)$ is maximized to find a new estimate $\Theta$. This quantity is defined as:

$$\begin{aligned} Q(\Theta, \Theta^k) = & -\frac{1}{2\sigma^2} \sum_{n=1}^{N} \mathbb{P}(z_n = 1 | x_n, v_n, \Theta^k) \|v_n - \bar{u}(x_n)\|^2 \\ & -\frac{D}{2} \ln \sigma^2 \sum_{n=1}^{N} \mathbb{P}(z_n = 1 | x_n, v_n, \Theta^k) \\ & + \ln(1 - \gamma) \sum_{n=1}^{N} \mathbb{P}(z_n = 0 | x_n, v_n, \Theta^k) \\ & + \ln \gamma \sum_{n=1}^{N} \mathbb{P}(z_n = 1 | x_n, v_n, \Theta^k) - \frac{\lambda}{2} \phi(\bar{u}) \end{aligned} \quad (4)$$

**E-step:** Given the current parameter $\Theta^k$, probabilities $p(z_n = 1 | x_n, v_n, \Theta^k)$ are estimated with

$$p_n = \frac{\gamma e^{-\frac{\|v_n - \bar{u}(x_n)\|^2}{2\sigma^2}}}{\gamma e^{-\frac{\|v_n - \bar{u}(x_n)\|^2}{2\sigma^2}} + \frac{1-\gamma}{\Omega}(2\pi\sigma^2)^{D/2}} \quad (5)$$

The probabilities $p_n$ quantify the agreement between the samples and current estimated regularized displacement field $\bar{u}$.

**M-step:** With:

* $P = \text{diag}(p_1, ..., p_N)$ the diagonal matrix that contains the probabilities $p_n$

* $\tilde{P} = P \otimes I_D$ where $I_D$ is the identity matrix of size D, and $\otimes$ is the Kronecker product

* $V = (v_1^T, ..., v_N^T) \in \mathbb{R}_{N \times D}$

* $\tilde{V} = (v_1^T, ..., v_N^T)^T \in \mathbb{R}_{ND \times 1}$

* $\tilde{\tilde{U}} = (\bar{u}(x_1)^T, ..., \bar{u}(x_N)^T)^T \in \mathbb{R}_{ND \times 1}$

The new parameters $\Theta = \text{argmax}_{\Theta} Q(\Theta, \Theta^k)$ are estimated using $P$.

Considering $\frac{\partial Q}{\partial \sigma} = 0$ and $\frac{\partial Q}{\partial \gamma} = 0$, estimations for $\gamma$ and $\sigma$ are obtained:

$$\sigma^2 = \frac{(\tilde{V} - \tilde{\tilde{U}})^T \tilde{P} (\tilde{V} - \tilde{\tilde{U}})}{D \cdot \text{tr}(\tilde{P})} \quad (6)$$

$$\gamma = \frac{\text{tr}(\tilde{P})}{N} \quad (7)$$

Then, considering the terms of $Q$ that depend on $\bar{u}$, maximizing $Q$ boils down to minimizing

$$\mathcal{E}(\bar{u}) = \frac{1}{2\sigma^2} \sum_{n=1}^{N} p_n \|v_n - \bar{u}(x_n)\|^2 + \frac{\lambda}{2} \phi(\bar{u}) \quad (8)$$

To solve this minimization problem, given that $\bar{u}$ lies within the Reproducing Kernel Hilbert Space (RKHS) $\mathcal{H}$ defined by kernel functions $\Gamma$, the smoothing filter $\phi$ then has to be $\phi(\cdot) = \|\cdot\|_{\mathcal{H}}^2$, showing that $\lambda$ identifies as a regularization parameter on $\bar{u}$. Then, the minimization problem amounts to solving the linear system

$$(\tilde{\Gamma} + \lambda \sigma^2 \tilde{\tilde{P}}^{-1}) \tilde{C} = \tilde{V} \quad (9)$$

with $\tilde{C} = (c_1^T, ..., c_N^T)^T \in \mathbb{R}_{ND \times 1}$ and $\tilde{\Gamma}$ a $N \times N$ block matrix such as the block $\Gamma_{i,j} = \Gamma(x_i, x_j)$ is a DxD matrix. Solving for $\tilde{C}$, then yields an estimate of $\bar{u}$.

The VFC Algorithm can be summed up in the following formulation, whereby user parameters are highlighted:

Its process can be summarized as follows: after initializing parameter $\Theta$, iterations are performed while the quantity to maximize has not stabilized. At each step, probabilities $p_n$ of every sample are computed by penalizing deviations from the current estimate of $\bar{u}$ (E-step), then, using these probabilities as weights, variables $\bar{u}, \gamma$ and $\sigma$ are updated (M-step).
**Algorithm 1** VFC Algorithm, with user parameters (incl. initialization of optimized ones) in red

**Require:** $N \ge 0$, $\{(x_n, v_n), x \in \mathcal{X}, v \in \mathcal{V}, n = 1 \to N\}$ positions/displacements, user parameters $\lambda, \Omega, N_{lim}$ and $\Gamma$

Compute Kernel $\tilde{\Gamma}$

$\tilde{u} \leftarrow 0$

$\gamma \leftarrow \gamma_{init}, \sigma \leftarrow \sigma_{init}$

iter = 0

**while** iter $\le N_{lim}$ **do**

**E-step:**

$p_n \leftarrow \frac{\gamma e^{-\|v_n - \tilde{u}(x_n)\|^2}}{\gamma e^{-\|v_n - \tilde{u}(x_n)\|^2} + \frac{1-\gamma}{\Omega}(2\pi\sigma^2)^{D/2}}$

**M-step:**

$\tilde{C} \leftarrow (\tilde{\Gamma} + \lambda\sigma^2\tilde{P}^{-1})^{-1}\tilde{V}$

Compute $\tilde{\tilde{U}}$ from (1)

$\sigma^2 \leftarrow \frac{(\tilde{V} - \tilde{\tilde{U}})^T\tilde{P}(\tilde{V} - \tilde{\tilde{U}})}{D.tr(P)}$

$\gamma \leftarrow \sum_{n=1}^{N} p_n$

iter $\leftarrow$ iter + 1

**end while**

return $\Theta = (\tilde{u}, \gamma, \sigma)$ and $(p_1, ..., p_N)$

After the E-M, the algorithm returns the set $\Theta$ as well as the probabilities $(p_1, ..., p_N)$ of the samples. A filter $(p_1, ..., p_N) > \tau$ is then applied to select the inliers. We will however show that the choice for threshold $\tau$ holds in fact very little consequence (see section 3.2.3).

## 2.3 Adapting parameters to PTV

Processing algorithms in 3D Particle Tracking Velocimetry often involve a series of successive steps, each with a specific set of parameters. For a method to be as efficient as possible in practical situations, the global parameter set should be minimal, and the associated optimal values should be easy to select by the user. We advocate that the F-VFC algorithm, gathering seven parameters as shown in the end of the previous subsection (two of them being initializations of quantities to be optimized, see block "Algorithm 1"), is in line with this idea. We highlight in this subsection the choices for a first series of them, which are systematic whatever the flow, or which can be set without having to perform characterization tests of F-VFC. Those tests, allowing to complete the optimal choices for the remaining parameters, will be shown in section 3.

Firstly, $\Omega$ can be directly set to the value of the volume derived from the search radius used to build the initial vector set. Considering then the choice of kernel $\Gamma$, which can be critical as it directly impacts the estimated $\tilde{u}$ and can have a strong influence on the computational cost, we will choose a Gaussian kernel $\Gamma_{i,j}(x_i, x_j) = e^{-\beta\|x_i - x_j\|^2} I_D$. Such a choice renders computations more tractable: the blocks $\Gamma_{i,j}$ are diagonal, and the linear system (9) becomes, with $K_{i,j} = e^{-\beta\|x_i - x_j\|^2}$,

$$ (K + \lambda\sigma^2 P^{-1})C = V $$

The temporal complexity of the minimization becomes $O(N^2)$. Parameter $\beta$ appearing in the Gaussian kernel can be rewritten, introducing a typical length of interaction between the data points

$$ \beta \sim \frac{1}{2(l^*)^2} $$

This length $l^*$ represents a typical length scale at which we assume some similarity between the samples. Also, since it underlines displacement field $\tilde{u}$, it should be chosen so that inliers can effectively obey a Gaussian distribution around $\tilde{u}$ (of standard deviation $\sigma$). As will be shown in subsection 3.2.1, $l^*$ can in fact be chosen easily, as a function of the mean inter-particle distance in the flow.

As mentioned above, further parameter choices (i.e. for $\lambda, N_{lim}, \tau$ and initialization values $\gamma_{init}$ and $\sigma_{init}$) will be detailed in the next section, as their optimal value will be confirmed through synthetic tests.

## 2.4 Handling large datasets: volume sub-division and double-pass processing

The VFC algorithm as implemented in Ma et al. (2014) returns a single parameter $\Theta$ to describe the whole vector field. This means that a single value of $\sigma$ must describe deviations of inliers from the smoothed flow $\tilde{u}$. This can be problematic in the case of turbulent flows, where behaviors and scales of interest can vary widely even in a small region. Furthermore, applying VFC on a complete experimental vector field can lead to significant increase of computational cost, given the very high number of vectors to be addressed in nowadays 3D PTV experiments. To address both these issues, we propose to divide the volume in sub-regions and to apply VFC separately on each. This allows to turn the large problem into a series of smaller problems that are much faster to solve, as well as to obtain a $\Theta = (\tilde{u}, \gamma, \sigma)$ that fits best each subdomain. In that case, the overall computational burden is dominated for matrix inversion on each subdivision, that grows as $N_s^2$ for each subdivision, where $N$ is the number of particles and $N_s$ the number of subdivisions. In the present implementation, we perform these inversions sequentially, whereas a parallelized version would save 90% of the actual processing time. Our implementation still provides results in reasonable computational times: this will be discussed in section 5.

One issue to tackle in that process is that the regularization's quality will tend to decrease around the edges of each subdomain, due to the absence of interaction between adjacent sub-domains. To address this, we choose to also apply VFC on overlapping subdomains, with an overlap of 50%: each sample $(x_n, y_n)$ is then located in two different subdomains A and B, for which two different results $(\tilde{u}_A(x_n), \sigma_A, \gamma_A)$ and $(\tilde{u}_B(x_n), \sigma_B, \gamma_B)$ will be available. We then reconstruct the final $\tilde{u}$ by computing the weighted average

$$ \tilde{u}(x_n) = \mu_A \tilde{u}_A(x_n) + \mu_B \tilde{u}_B(x_n) $$

where $\mu_A$ and $\mu_B$ measure the normalized distance of sample $x_n$ to the center of sub-domains A and B. Similarly, we
estimate $\sigma_n$ and $\gamma_n$ at the location $\boldsymbol{x}_n$ and use those to estimate the final probabilities $p_n$.

For large and noisy datasets (such as those that we will encounter in our 3D test cases), having a large number of E-M iterations may be necessary to reach convergence, but will be time-consuming. At each iteration, the matrix inversion is slowed down by large amounts of spurious vectors. We propose a double-pass method to address this issue. F-VFC is first ran on the full dataset with a limited number $N_{lim}$ of E-M iterations. Although tuning this $N_{lim}$ accordingly to the dataset might seem challenging, we will show in 3.2.2 that the behavior of the complete-data log-posterior $\mathcal{Q}$ can easily help stop the E-M iterations at the right time. After this first pass, the most obvious spurious vectors are identified and removed, while outliers that are “close enough” to the correct behavior are kept. At this point, the estimated $\Theta$, although much better than its initialization, is still far from optimal: $\sigma$ is large, and $\bar{u}$ is poorly estimated. However, the corresponding dataset is already much less noisy than the initial one, and, above all, much smaller. This allows to run a much faster second pass, on the output of the first pass, to reach convergence in a much quicker way than in the single-pass mode. We will show in subsection 3.2.2 that on top of increasing processing speeds, this double-pass implementation is also beneficial to the final performance.

## 2.5 Complete algorithm

A full overview of the F-VFC algorithm, encapsulating the VFC method with adaptations and additions dedicated to PTV, is given in Fig. 2. Starting from particle reconstructions at the two instants, potentially gathering a large number of ghosts, the initial vector set is built by pairing each particle of the first time step with all possible matches in the second time step within a given search radius. If the flow domain is large, in particular with respect to the spatial scales, the vector set can then be optionally subdivided into subvolumes. The adapted VFC algorithm is then ran on each subvolume (or on the full volume), using a Gaussian kernel, whereby parameter $\beta$ is tuned according to the mean inter-particle distance. E-M iterations are stopped when the variation of the complete-data log-posterior $\mathcal{Q}(\Theta, \Theta^k)$ has reached a maximum. After this first pass, a second pass can be ran on the obtained output, if the vector field is judged to still contain a too large number of outliers; in this second pass, the number $N_{lim}$ of iterations is fixed, instead of considering $\mathcal{Q}$. The justification for choosing a different stopping criterion between these two passes will be provided in subsection 3.2.2. While the necessity of it has not been encountered yet in our tests, additional passes could be ran next. After the last pass, the vector field is scanned for duplicates, i.e. multiple vectors originating from the same particle in the first time step; for each of these situations, the vector kept among the duplicates is the one with largest probability of being an inlier. This yields the final displacement set.

# 3 Performance and operating conditions on a 2D synthetic flow

In this section, we introduce 2D synthetic tests in order to provide first performance characterizations of F-VFC, together with guidelines for optimal parameter tuning by the user. As the rule for selecting $\Omega$ is known (see subsection 2.3) these pertain to regularization $(\lambda, \beta)$, stopping criterion for the E-M iterations, outlier selection $(\tau)$, and initialization $(\sigma_{init}, \gamma_{init})$. After a description of the synthetic dataset in subsection 3.1, subsection 3.2 uses analysis of the obtained results to confirm that optimal values for each of these parameters, successively, can be chosen either on a physical basis, within a large range of values due to a weak impact on overall performance. Subsection 3.3 then analyzes the robustness of F-VFC to outliers on this 2D case.

## 3.1 Test case setting

### 3.1.1 Flow field

The synthetic flow fields used to generate the 2D cases are taken from a well resolved 3D Direct Numerical Simulation (DNS) of a forced Taylor Green Turbulent flow at $R_\lambda = 62$. A pseudo-spectral scheme on a $768^3$ periodic cubic box of length $2\pi$ was used to solve the Navier-Stokes equations. The forcing is introduced through the Taylor-Green modes. From this 3D flow, a 2D velocity field is simply obtained by performing an extraction along a chosen plane, getting rid of the velocity component normal to the plane. It should be emphasized that this does not yield a strictly planar flow; in particular, its 2D divergence will not be zero. Particles are then randomly generated in the plane and advected along the obtained 2D flow. This process is repeated along different planes of the volume, at different time steps, so as to obtain 24 different vector fields. One of these fields is represented in Fig. 3 (Left).

<figure><img src="image_4.png" /><figcaption>**Fig. 3 (Left):** Sample velocity field, featuring the way for selecting domain size $D_{Size}$, and thus the variability of the flow. (Right): inliers and outliers inside a domain of size $26\eta$, with a fraction of outliers %O= 15%</figcaption></figure>

### 3.1.2 Experimental parameter sweep

**Inter-particle distance** The flow seeding density, known to be the most influential parameter in 3D PTV, is here controlled through the mean inter-particle distance, denoted
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Overview of the F-VFC algorithm. Given an input dataset of two particle clouds, we initially pair all particles with candidates inside a given search volume. This leads to an initial vector set with an extremely large amount of outlier vectors. The user can choose to divide the domain in N subdivisions. The VFC algorithm is run in each subdivision, with a gaussian Kernel Γ. Iterations are stopped when variations of the log-posterior (the quantity the EM aims at maximizing) is maximal. Then a second pass is run on the output of the first pass: a fixed number of iterations is used. Eventually, the user gets the probability of being an inlier for each vector: vectors of low probability are discarded, and duplicates (vectors originating from/arriving at the same point) are handled by keeping only the most probable duplicate.</figcaption></figure>

by $\langle d \rangle$. We will also use its dimensionless counterpart with respect to the typical flow scales, which we define as $d^* = \langle d \rangle / \eta$. We can anticipate that variation of $d^*$ will allow characterizing the amount of information necessary for optimal performance of F-VFC; indeed, having a small $d^*$, i.e. a significant number of inliers separated by only a few $\eta$ (thus sharing very similar behavior), will make the consensus easy to find, whereas this will become more difficult with increasing $d^*$. This reasoning shows in particular that other physical quantities,

such as the flow turbulence rate $I = \frac{u_{RMS}}{\langle u \rangle}$ (for instance used in Mercier et al. (2024)), will not be of influence here.

**Domain size** As depicted in Fig. 3 (Left), we consider variable domain sizes $D_{Size}$, which are multiples of the Kolmogorov scale $\eta$, up to 52 times. Considering this specific maximal value results in the fact that $\eta$ is in effect not much smaller, instead quite comparable, to the largest flow scales,
and should be considered as such (rather as being separated by a significant range of scales from the largest scales, as in developed turbulence). Considering different and in particular increasing $D_{Size}$ will allow tackling volumes with an increasing variety of scales, leading to a performance drop of single-volume F-VFC processing (as the latter assumes a single value of $\sigma$ for the whole volume); symmetrically, this situation will help grasp the benefit brought by a subdivision approach instead.

**Fraction of outliers** In order to preserve a simple way to control the amount of outliers in these 2D datasets, while choosing them in a realistic way compared to relevant PTV situation, we here deviate from the process arising from systematic pairing described in section 2. In practice, outliers will arise both from ghost particles present in the reconstruction, as well as from incorrect pairing (be it between two true particles, two ghosts or between a true particle and a ghost). In order to keep the fraction of outliers as a single parameter, we restrict to the case of outliers resulting from ghost to ghost pairing; however, we choose them to be quite close in norm to inliers. More precisely, we introduce outliers by adding vectors of random position, orientation, and of random magnitude comprised between 0.5 and 1.5 times the mean local displacement field $||\bar{u}||$ - computed as the average displacement magnitude inside the domain. These outliers have indeed relatively similar magnitude as the inliers (but random orientations, as seen in Fig. 3 (Right)). Note that this is in fact not a simplification of the situation considered in section 2, rather a worst-case scenario: indeed, outliers arising from wrong pairing of true particles, or of pairing of a true particle with a ghost, will at worst be similar to the outliers considered here (i.e. close in norm to inliers), or otherwise have a norm (potentially very) different from that of neighboring inliers, thus making them easy to filter out.

In the sequel, we will note %O the fraction of outliers in the dataset, evaluated as the number of outliers relative to the total number of vectors. If F-VFC were to yield a perfect result on these datasets, then the final estimated value of $\gamma$ would identify to 1 - %O.

### 3.1.3 Processing and performance metrics

Unless otherwise specified, the volume will be treated by F-VFC without any subdivisions: this will allow to assess the robustness of the algorithm when handling large domains with varying flow behaviors inside. In the same way, single-pass processing will only be considered.

To assess the effectiveness of the outlier removal procedure, we will quantify the number of true positive vectors TP, of false negative vectors FN (a real vector was discarded) and false positive vectors FP (an outlier was kept). With these numbers we construct two measures of performances, similarly to what is done in Champagnat et al. (2014):

* the Recall = $\frac{TP}{TP+FN}$

* the Precision = $\frac{TP}{TP+FP}$

A Recall equal to 1 means that all the real vectors in the data set have been found and kept. Similarly, a Precision

equal to 1 indicates that all the vectors that were kept are indeed correct vectors.

Note that for several of the test cases presented in the next subsections, to increase convergence, averages will be performed on results obtained with different domain sizes $D_{Sizes}$, of course picking values of this quantity yielding similar performance; in such situations, for clarity we will omit the exact corresponding values gathered in the averaging.

## 3.2 Optimal parameter selection

### 3.2.1 Regularization

The two parameters that impact the regularization and the computation of $\bar{u}$ are $\beta$ (in the Gaussian Kernel) and $\lambda$ (in the Tikhonov regularization). To comprehensively study the robustness to each parameter would require varying each parameter over a wide range of scales, for a varying range of experimental parameters. However, we have observed on all our test cases that varying $\beta$ has very little effect on conclusions drawn for $\lambda$: we can thus study the effects of these parameters as if they were decoupled.

The most impactful and critical parameter is scalar $\beta$ appearing in our Gaussian Kernel: as explained in subsection 2.3, upon writing $\beta = 1/2l^{*2}$, $l^*$ appears as a typical length scale of $\bar{u}$, and therefore a length scale at which the velocity field of inliers obeys a Gaussian distribution around $\bar{u}$. Obviously, in order to find a consensus and fit the Gaussian distribution, enough inliers need to be present in the volume. The choice of $\beta$ must then be a compromise between having a large enough window of interaction that enough inliers lie in the window, and having a small enough window of interaction that the hypothesis of a Gaussian distribution holds. Physically, one could thus link $\beta$ to a typical turbulent scale in the flow, which, in our case, is the Kolmogorov scale $\eta$: we then expect that the scale at which the Gaussian distribution hypothesis holds will be a multiple $M_{\eta}$ of $\eta$. This leads to expressing $\beta$ as

$$ \beta = \frac{1}{2(M_{\eta}\eta)^2} \qquad (10) $$

Experimentally however, when parameter $\eta$ is often difficult to estimate, it may be easier to think of the typical range of interaction between samples as a multiple $M_d$ of the mean inter-particular distance $\langle d \rangle$, essentially looking at interactions at the scale of $M_d$ times the inter-particular distance. This leads to expressing $\beta$ as

$$ \beta = \frac{1}{2(M_d \langle d \rangle)^2} \qquad (11) $$

We show in Fig. 4 the F-VFC performances when varying $M_d$, for different values of $d^*$ and %O = 50% of outliers. These results have been averaged over 24 different test cases. We make two observations on the behavior of $\beta$:

* To optimize performances, the Gaussian kernel needs to be scaled with a length $10 \langle d \rangle$. For higher values, we see a dip in Precision: all samples are kept by the algorithm, which, concomitantly, naturally increases the Recall to 1.
* At increasingly high seeding densities (i.e. small values of $d^*$) the acceptable range for $\beta$ spreads: this indicates that in cases where the seeding density is high, misestimations of the parameter $\beta$ are inconsequential.

We better understand the interaction between F-VFC and the number of inliers in the domain: the best performances (Recall and Precision over 95%) are reached when dimensionless inter-particle distance $d^*$ is equal to or smaller than 3, and begin to drop for larger values. Given that $d^* = \langle d \rangle / \eta$, with $\eta$ the Kolmogorov scale, this limit value, below which we will stay in the remaining tests of this section, can seem worrying at first, as real world 3D experiments are often seeded at scales quite larger than the Kolmogorov scale. However, we recall here that our flow domain is already quite zoomed in around $\eta$, so that the latter is quite comparable to large flow scales. Also, we think that this quite low value of $d^*$ corresponding to optimal performances can be pinned on the non-incompressible character of the present 2D flow, which is in fact a cut within a 3D flow. Indeed, we will observe in the 3D tests performed in sections 4 and 5 that even much lower seeding densities are in fact still efficiently dealt with by F-VFC, up to $d^* = 9$.

<figure><img src="image_6.png" /><figcaption>Fig. 4 Recall and Precision when varying $M_d$ in $\beta$ (see equation (11)) with $\%O = 50\%$, for different inter-particle distances ($d^* = 1, \bullet$), ($d^* = 2, \diamond$), ($d^* = 3, \triangleleft$), ($d^* = 5, \star$), averaged over 24 different test cases.</figcaption></figure>

We now turn to parameter $\lambda$, which, in the minimization problem (8), indicates the compromise between closeness to the data ($\frac{1}{2\sigma^2} \sum_{n=1}^N p_n \|v_n - \bar{u}(x_n)\|^2$) and smoothness of the estimated field ($\phi(\bar{u})$). We study in Fig. 5 both the Recall and Precision when varying $\lambda$, for a fraction of outliers $\%O = 66\%$ and different values of $d^*$. As Fig. 5 shows, and this has been confirmed by other tests with different experimental parameters, we observe that there is a large range of $\lambda$ for which we recover optimal performances, typically from 0.001 to 1000 without significant impact on Recall or Precision. An exception to this is the case of low seeding densities, where the number of inliers in the volume can become insufficient to provide a consensus; in such situations, higher values of $\lambda$ (up to 100) have been

found to help improving performances. Note that this type of robustness is also in line with observations by Ma et al. (2014).

<figure><img src="image_7.png" /><figcaption>Fig. 5 Recall (●) and Precision (×) when varying $\lambda$ for $D_{Size} = 26\eta$, with $\%O = 66\%$ of outliers, averaged over 6 test cases. From top to bottom $d^* = 0.8, d^* = 1.2, d^* = 2$.</figcaption></figure>

### 3.2.2 Stopping Rule for E-M Iterations

We here provide justifications for our proposed choice for stopping the Expectation-Maximization iterations in the F-VFC algorithm, as depicted in Algorithm 1 and Fig. 2.

Iterations of the E-M algorithm are used to reach convergence of the complete-data log-posterior $Q(\Theta, \Theta^k)$. We write, as in Algorithm 1, $\delta E$ its relative variation between iterations $k$ and $k + 1$. Thus,

$$ \delta E = \frac{Q(\Theta^{k+1}, \Theta^k) - Q(\Theta^k, \Theta^{k-1})}{Q(\Theta^{k+1}, \Theta^k)} $$

Practically, one could choose to stop iterating when $\delta E$ is low enough (i.e. below a given threshold $\delta E_{max}$ set by a user), meaning that the complete-data log-posterior has converged to a maximum. Another option would simply be to stop once enough E-M iterations have been done; in that case, the user would have to guess an appropriate number $N_{lim}$ of iterations for a good enough convergence of $Q(\Theta, \Theta^k)$. While we will retain the stopping rule at a fixed values of iterations $N_{lim}$ in certain cases, we will also show that a finer examination of the behaviour of $\delta E$, concomitantly with the quality of the vector set, can help defining a new stopping criterion, more adapted to other cases.

In Fig. 6, we show the Recall and Precision obtained when varying the number $N_{lim}$ of E-M iterations. Results correspond to averaging over 16 different test cases, generated with varying seeding densities and proportions of outliers $\%O$, all with a domain of size $D_{Size} = 40\eta$.
<figure><img src="image_8.png" /><figcaption>Fig. 6 Recall (top) and Precision (center), relative to the maximum value they reach, and δE (bottom), when varying the amount of total E-M iterations Nlim, for different inter-particle distances and proportions of outliers.</figcaption></figure>

In the very first E-M iterations, VFC recovers more and more data, increasing both Recall and Precision: the most spurious outliers are eliminated, so that the estimate $\bar{u}$ is improved, becoming closer to the inliers, the estimated $\sigma$ remaining large. During these steps, the log-posterior varies widely and $\delta E$ increases. Interestingly, Recall peaks at the same number of iterations where the variation of the log-posterior is maximal. Additionally, it is also interesting to note that, with an increasing fraction of outliers %O, it takes more iterations to reach this peak. A distinct behavior is observed with the Precision, as it increases, irrespective of the initial behaviour of $\delta E$; however, its slower convergence towards a maximum is symmetrical with the decrease of $\delta E$ towards its minimum, i.e. convergence of $Q$. As this increase in Precision at large number of E-M iterations is very slow, adding too many of these iterations is finally detrimental, as this also leads the Recall to decrease, diminishing the overall quality of the result. As a more detailed analysis of the vector field shows (not presented here for conciseness), this behaviour is observed because $\sigma$ reaches values too small to correctly recover the various flow behaviors in the large domain considered. Certain correct vectors do not perfectly fit the estimated $\bar{u}$ and are thus considered outliers by VFC.

As a first conclusion based on these trends, we recommend limiting the amount of E-M iterations to $N_{lim} \sim 10$, or up to $N_{lim} \sim 15$ for the noisiest datasets. Any higher value will lead to over-fitting and to a decrease in overall performances. Furthermore, this study points to an attractive feature of

VFC when used for Particle Tracking purposes: the user can indeed tune the amount of E-M iterations to favor either Precision or Recall. In practice, this possibility directly justifies the interest of the double-pass implementation that we introduced in section 2, and will use for the 3D test cases (see 4.1.3). Indeed, within this implementation we choose to maximize Recall (i.e. keeping as many inliners as possible) during the first pass, by stopping the E-M process when $\delta E$ is maximal (which happens for a rather small number of iterations as observed in Fig. 6). Then, in the second pass, to ensure maximal Precision in the output (i.e. least possible number of outliers), we rather switch to a stopping criterion based on $N_{lim}$, at a value comprised between 10 to 15, so that $\delta E$ becomes low enough and nearly stable between successive iterations, ensuring convergence of the complete-data log-posterior.

### 3.2.3 Outlier selection

Parameter $\tau$ is only used at the very end of the F-VFC process, to select inliners based on their estimated probability $p_n$. In practice, we have observed that for a large majority of vectors in the sets, $p_n$ converges either towards 1, or towards $1e-5$ (the lower limit we imposed for probabilities). More precisely, in our test cases, after 20 E-M iterations, more than 99% of the samples had either converged to probabilities below 0.01 or above 0.99. The choice of $\tau$ will thus influence the selection of only a rather small amount of data, so that its precise value is not critical. We thus choose to set $\tau = 0.5$ to remain conservative with the samples, retaining every sample that is more probable than not.

### 3.2.4 Initialization

Parameters $\sigma$, $\gamma$ and $\bar{u}$ need to be initialized before the first Expectation step. As we do not include in F-VFC a preliminary algorithmic step aiming at finding a predictor for the flow displacement, setting $\bar{u} = 0$ is the most adapted and conservative option. In such a case, the value of $\sigma_{init}$ can be derived from (6), with $\tilde{\bar{U}} = 0$:

$$ \sigma_{\text{init}}^2 = \frac{\tilde{\mathbf{V}}^T \tilde{\mathbf{P}}_{\text{init}} \tilde{\mathbf{V}}}{D \cdot \text{tr}(\mathbf{P})} = \frac{1}{ND} \sum_n \|v_n\|^2 \quad (12) $$

where $\tilde{\mathbf{P}}_{\text{init}} = I_{ND}$. This strategy of initialization is consistent with the one chosen in Ma et al. (2014).

In practical cases, estimating a value of $\gamma$ can however be challenging. Interestingly enough, we have observed that its given initial value does not affect performances in any manner: a “bad” initialization will yield similar results, although in a few more iterations. To avoid having to run too many E-M iterations, we thus suggest setting $\gamma$ close to its real value if it can be inferred in any way, otherwise simply to 0.5. Note that this choice is in fact different from the implementation presented in Ma et al. (2014), where the authors initially assume $\gamma = 0.9$, due to their lower number of outliers in the dataset.
### 3.2.5 Summary of optimal parameter values

As a synthesis of this subsection, we present in Table 1 the full set of optimal values suggested for the F-VFC user parameters, together with options for finer tuning possibilities, if necessary and accessible. Its last line, regarding the number of volume subdivisions, is discussed in the next section.

## 3.3 Robustness to outliers and interest of volume subdivision

Now that guidelines for selecting optimal parameters are known, we study the robustness of our algorithm to high proportions of outliers, and consider increasing domain sizes in order to assess the interest of using volume subdivision within F-VFC. Fig. 7 shows performances of the algorithm when applied to two volumes sizes $D_{Size} = 26\eta$ and $52\eta$, when increasing the proportion of outliers %O. Consistently with the analysis pertaining to the choice of $\beta$ in subsection 3.2.1, we also consider different values of the seeding density, i.e. of $d^*$. As the focus of the present test is to evaluate the potential of F-VFC to filter out outliers, we here choose to consider a fixed number of E-M iterations, with the stopping criterion set to $N_{lim} = 12$ iterations. As mentioned in section 3.2.2, such a setting tends to favor Precision over Recall, i.e. preferentially filtering outliers, potentially at the cost of a small number of inliers.

Focusing first on the cases without volume subdivision, we find F-VFC to be overall very robust to outliers. At high seeding densities, here for $d^* = 1$, we find very satisfying results (typically, Precision over 90%, Recall over 80%) for up to 90% of outliers. At lower seeding densities, performances are still very satisfactory but lower, which is not surprising: the flow is not scanned at a high enough spatial frequency that all behaviors can be recovered. Over the various cases, F-VFC indeed favored the Precision over the Recall, as expected by our choice of stopping criterion.

At a given seeding density, we find that the larger the domain (compare the leftmost and rightmost columns of Fig. 7), the worse the performances. This is expected, as we use a single parameter $\theta = (\bar{u}, \sigma, \gamma)$ to describe the whole flow field, whereas it exhibits smaller scales of spatial variation than the domain size. In such a situation, a single value of $\theta$ cannot describe the behavior of all the inliers. We show in the central column of Fig. 7 the performances of F-VFC when the larger domain, of size $D_{Size} = 52\eta$, is divided in sub-windows of size $D_{Size} = 26\eta$, with 50% overlap (see section 2.4). As expected, performances become comparable to those obtained for a single smaller domain of $D_{Size} = 26\eta$: this fully justifies the relevance and interest of this subdivision method for complex and turbulent flows. This also underlines the choice we will be doing in the following sections. We will indeed look to have a number of subdivisions $N_S$ in the volume as high as possible, to both improve on performance and limit computational costs. It should however be noted that a limit should be respected, as it remains important to keep enough inliers in each subdivision that it is possible to find a consensus. We found that maintaining around 100 inliers per subdivision serves as a good rule-of-thumb.

In practice, such an information can be obtained from the estimated seeding density in the images.

# 4 Performance assessment on 3D synthetic cases

The F-VFC Algorithm was found to robustly filter out outliers from 2D vector sets with large amounts of outliers. In this section, we extend its use to 3D PTV datasets. As we saw in subsection 3.2.2 that the stopping criterion on the E-M iterations of VFC may be used as a way to select the performance metrics to be favored, between Recall and Precision, we will here choose to favor Recall, i.e. to maximize it as much as possible, but not at any cost, i.e. also maintaining a (very) good Precision. This choice, which can seem surprising, e.g. compared to a compromise between both quantities, is in fact motivated by the fact that, as shown by Cornic et al. (2021), the final position and thus velocity error of a 3D PTV algorithm can be more heavily influenced by the missed particles/vectors (FN), than by ghosts/outliers. In the future perspective of integrating F-VFC within an iterative 3D PTV algorithm, we thus believe that this choice of favoring Recall is more promising.

In the sequel of this section, we consider two different test cases. The first one, presented in subsection 4.1, was especially designed with in-house tools to allow a large and detailed range of seeding densities. In order to facilitate future comparison with other methods, we also provide results on a second one, which is publicly available: the two-pulse dataset of the First LPT Challenge (Leclaire et al., 2021; Sciacchitano et al., 2021). This is presented in section 4.2.

## 4.1 First test case

### 4.1.1 Setup

A flow domain of extensions $80 \times 50 \times 6 \text{ mm}^3$, in the $X$, $Y$ and $Z$ directions, is filled with randomly positioned particles, which are then advected using the forced homogeneous isotropic turbulence dataset from the John Hopkins Turbulent Database (Li et al., 2008; Perlman et al., 2007) (Reynolds number $Re = 157000$, Taylor-Scale Reynolds number $Re_\lambda = 418$). In order to converge statistics, we use several different instants of the database to generate distinct synthetic datasets; the velocity magnitude in the middle plane of one of them is represented in Fig. 8.
<table><caption>**Table 1** Suggested optimal values for F-VFC's user parameters.</caption><thead><tr><th></th><th>Parameter</th><th>Optimal value</th><th>Possible finer tuning</th></tr></thead><tbody><tr><td rowspan="2">Regularization</td><td>λ</td><td>1</td><td>Increase if low seeding density</td></tr><tr><td>β</td><td>β ~ <sup>1</sup>&frasl;<sub>2(10⟨d⟩)2</sub></sub></td><td></td></tr><tr><td>Stopping criterion of iterations</td><td>N<sub>lim</sub>, δE</td><td>N<sub>lim</sub> = 10 - 15<br>Stop when δE is maximal</td><td></td></tr><tr><td rowspan="2">Outlier selection</td><td>Ω</td><td>πR<sub>S</sub><sup>2</sup> (D = 2) or 4/3πR<sub>S</sub><sup>3</sup> (D = 3)<br>(with R<sub>S</sub> search radius of initial pairing)</td><td></td></tr><tr><td>τ</td><td>0.5</td><td>Increase/decrease to be more/less selective (but weak impact in practice)</td></tr><tr><td rowspan="3">Initialization</td><td>u&#773;</td><td>0</td><td></td></tr><tr><td>σ</td><td>from equation (12)</td><td></td></tr><tr><td>γ</td><td>0.5</td><td>Closer to real value if known</td></tr><tr><td>Number of subdivisions</td><td>N<sub>S</sub></td><td>Such that each subdomain should contain a minimum of 100 inliers (possible estimation from seeding density)</td><td></td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 8** Observed volume of 80 x 50 x 6 mm³ and isocontour of the velocity magnitude in the middle plane of the volume.</figcaption></figure>

A virtual illumination and imaging setup is considered, with a uniform illumination over the 6mm thickness, and a four-camera system placed along the arc of a circle, at Y = 0 and centered at the volume centre, of radius R = 0.9m. The two external cameras are positioned at an angle of ±30° relative to the Z axis, and the internal ones at ±10°. Each of them has a 1920 x 1200 pixels sensor with a 10μm pitch and is equipped with an f = 200mm lens. No Scheimpflug adapter is considered for simplicity, and calibration is supposed to be error-free. Using this setup, synthetic images are generated by considering uniformly polydispersed particle intensities, particle image sizes of 2.4 pixel, as well as typical levels of thermal noise. The inter-frame time was chosen such that the average displacement is around 5 pixels, with maximum values ranging between 14 and 20 pixels. Seeding densities ranging from 0.014 to 0.143 particles per pixel (ppp) will be considered.
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Recall(●) and Precision(×) when varying %O for different sets of extrinsic parameters, and for processings with volume subdivision (central column, highlighted in red) and without subdivision (left and right columns). We vary the domain size and the normalized mean inter-particular distance $d^*$. (Left): Domain of size $D_{Size} = 26\eta$; (Center): Domain of size $D_{Size} = 52\eta$ divided in sub-windows of size $D_{Size} = 26\eta$ with 50% overlap; (Right): Domain of size $D_{Size} = 52\eta$.</figcaption></figure>

### 4.1.2 Input particle clouds and initial pairing

To obtain the input particle reconstructions at the two in-
stants, we use the Particle Volume Reconstruction (PVR)
approach introduced by Champagnat et al. (2014). We first
discretize the volume into a voxel grid with voxel-to-pixel
ratio of 0.5, and threshold the images at a fixed level to de-
crease the impact of noise. The tomographic reconstruction
is initialized using the Multiplicative Line-Of-Sight approach
(Atkinson and Soria, 2009). Then, 25 iterations of SMART
(Simultaneous Multiplicative Algebraic Reconstruction Tech-
nique) are run, using the PVR weighting matrix modelling
the cameras' point spread function (i.e. accounting for the
size of the particle images). To finally obtain candidate 3D
particle positions out of the result, which, after the SMART
iterations, still consists of equivalent grey level on a voxel
grid, we perform a three-point 3D Gaussian fit over the local
intensity maxima; candidates for these maxima are selected
along voxels whose intensity is above the value used for initial
image thresholding. The fact that we use in this process
voxels which are half of the size of a back-projected pixel
is linked to the optimal working conditions of the PVR ap-
proach, as described in Champagnat et al. (2014). In the
sequel, we will also refer to voxels as a characteristic scale;
note that these will then be in a more traditional sense, i.e.
with size of a back-projected pixel.

From these particle reconstructions, the initial pairing
providing the input vector set to F-VFC is done, as described
in subsection 2.1, by considering all matches possible inside
a search radius around each particle of the first time step

(see Fig. 1). The search radius is scaled on the maximum
expected displacement in the dataset. This yields a vector
set with a large number of outliers, but with every correct
vector contained in the input particle clouds.

Similarly to section 3, we will adopt the Recall and Precision performance metrics. For these 3D tests however, we will have to distinguish whether these quantities pertain to the particle clouds or to the vector sets. Considering the input clouds, a given particle will be considered a TP if a ground truth (GT) particle is found in a neighborhood of 1 voxel, and otherwise an FP. If no particle exists in the input within a 1 voxel neighborhood of a GT particle, then an FN will be counted. For the vector sets, a vector will be considered a TP if formed from the matching of two TPs, these being the same GT particle. Particle Precision, as well as Vector Recall and Precision, will then be computed from these quantities. In the following, as the purpose is the performance assessment of F-VFC, we will not provide values of the Particle Recall for the different cases. Here, it is simply sufficient to mention that, consistently with the SMART approach, this quantity is systematically higher than 0.95, whatever the seeding density; this indicates that the amount of inliers in the vector sets is consistent with typical values expected at the corresponding seeding densities, in the case of a high quality reconstruction. For the same reason, Vector Recall will be understood relative to the input particle reconstructions, i.e. the counting of potential FN vectors will be performed with respect to true particles of the reconstructions, not to that of the ground truth.

We write in Table 2 the Particle Precision and Vector
<table><caption>**Table 2** Precision of the input particle reconstructions and of the initial vector set at varying seeding densities in particles per pixel (ppp).</caption><thead><tr><td>ppp</td><td>Reconstructed Particles</td><td>Initial Vector Set</td></tr></thead><tbody><tr><td>0.014</td><td>0.990</td><td>0.454</td></tr><tr><td>0.030</td><td>0.968</td><td>0.273</td></tr><tr><td>0.046</td><td>0.918</td><td>0.179</td></tr><tr><td>0.062</td><td>0.848</td><td>0.119</td></tr><tr><td>0.077</td><td>0.769</td><td>0.082</td></tr><tr><td>0.094</td><td>0.687</td><td>0.057</td></tr><tr><td>0.108</td><td>0.615</td><td>0.039</td></tr><tr><td>0.126</td><td>0.558</td><td>0.029</td></tr><tr><td>0.143</td><td>0.510</td><td>0.022</td></tr></tbody></table>

Precision of the initial vector set, for the explored range of seeding densities. With this dataset, we are able to assess the robustness of our algorithm to various proportions of ghost particles, from 1% to 49%. As expected, the method chosen to build the initial vector set generates a large amount of outliers, their fraction ranging from %O from 54.6% to 97.8 %, since %O is equal to 1 - Precision.

### 4.1.3 Processing parameters

The presents dataset turns out to be extremely large (up to $10^7$ vectors for the densest case), with high proportions of outliers. To improve on performance and computational efficiency, we leverage on the volume subdivision and double-pass options of F-VFC described in section 2.4. The volume is divided in subdomains, of size chosen such that they each contain around 100 inlier vectors, as discussed in section 3.3. For this test case as well, maintaining small subdivisions allows us to maintain a relatively low computational cost while improving on performances. Consequently, the number of subdomains is different for the different seeding densities. Then, a double-pass VFC Algorithm is run. Consistently with the approach described in section 3.2.2, in the first pass, we try to maximize Recall, by stopping the E-M iterations when the variation $\delta E$ of complete-data log-posterior is maximal. Also, $\tau$ is set at a lower value of 0.1. This ensures that we recover a maximum number of correct vectors while “cleaning” the dataset from evidently spurious vectors created by our initial pairing method. Then, a second pass is run on the output dataset, with $N_{lim} = 10$ iterations, in the objective of increasing the Precision while maintaining the Recall at the previously obtained high level.

The set of remaining user parameters is chosen as follows:

* Regularization: $\lambda = 1$ and $\beta = \frac{1}{2 \cdot (10 \langle d \rangle)^2}$ where $\langle d \rangle$ is the estimated mean inter-particular distance in the domain, based on the image seeding density.

* Outlier selection: we set $\Omega = 4/3\pi R_S^3$ with $R_S = 16$ voxels, $\tau = 0.1$ on the first pass, and $\tau = 0.5$ on the second pass.

* Initialization: $\sigma_{init}$ initially set as explained in section 3.2.4, and $\gamma_{init} = 0.5$.

### 4.1.4 Results

<figure><img src="image_11.png" /><figcaption>**Fig. 9** Vector Recall (●) and Precision (×) for the different test cases at varying ppp. Recall is here understood relative to the input particle reconstructions. The dashed lines indicate corresponding performance metrics after the first pass of the VFC Algorithm. Solid lines show results after the second pass.</figcaption></figure>

We show in Fig. 9 the Vector Recall and Precision of F-VFC’s output. It is found that, up to a seeding density of 0.126 ppp, corresponding to reconstructed particle clouds containing around 45% of ghost particles, the Recall of F-VFC remains above 90%. Only at the higher ppp of 0.143 do we get a Recall below 90% (although still quite high, above 82%). As expected and in a way targeted, Recall is kept very high at the expense of Precision. Still, we find that we get relatively good performances in terms of Precision, staying above 80% for up to ppp 0.126. We do find that, above that density, the Precision seems to collapse. Our interpretation of this behaviour is that the corresponding remaining outliers which contribute to that drop in Precision are in fact outliers whose characteristics are similar to that of inliers (i.e. within the probability distribution of inliers), so that, by construction, VFC cannot filter them out. Note that, as seeding density increases, the occurrence of such outliers close to inliers in nature should become more and more probable. To support this hypothesis, we analyze selectively the error (in voxels) with respect to the ground truth vectors of both the TP and FP vectors in F-VFC’s output. By definition of our performance metrics, a vector must deviate from the Ground Truth by a maximum of 2 voxels to be considered a TP. Our hypothesis will thus be supported if we find a large number of FP/outliers with errors lower than this value. At ppp = 0.126, we find that F-VFC recovered around 35000 False Positives, for about 216 000 True Positives. As shown in Fig. 10, of all those FP, 65% were indeed found to exhibit an error with respect to the GT similar to that of the TP. This confirms that a very significant portion of the FP are coherent with the flow.

Overall, the present analysis shows the robustness of F-VFC to high proportions (up to 45%) of ghost particles in the particle clouds. Note that the high Recall values reached, in particular, can mainly be pinned on the double-pass implementation. As shown by comparing the curves
<figure><img src="image_12.png" /><figcaption>**Fig. 10** Probability Distribution Function of the True Positives(×) and False Positives(●) at ppp 0.126, as a function of the displacement error with respect to the GT. (Inset): cumulative histogram of the False Positives as a function of the displacement error with respect to the GT. The blue vertical line shows the maximum possible error for the TP.</figcaption></figure>

corresponding to the first and second pass outputs of F-VFC in Fig. 9, we can confirm that the choice of stopping criterion proposed in 3.2.2 indeed helps reaching the targeted behaviour. After a rather non selective first pass, we succeed at maintaining a high Recall, while having eliminated few absurd vectors: the Precision is still poor. The second pass then manages to greatly improve on Precision, while not degrading the Recall too much.

## 4.2 Second test case: First LPT Challenge

To further assess the robustness of F-VFC for Particle Tracking at high seeding densities, we challenge it with the two-pulse 3D synthetic data from the First LPT Challenge (Leclaire et al. (2021), Sciacchitano et al. (2021)). Similarly to the previous case, we first reconstruct particles using the PVR-SMART algorithm (see section 4.1.2 for more details). Performances of this reconstruction step in terms of Particle Precision and Recall, obtained by performing an upload with automatic evaluation on the dedicated portal (https://www.onera.fr/flow-benchmarks/, see Leclaire et al., 2022) are shown in Table 3. Weaker reconstruction performances are observed compared to the previous test case, visible by a higher number of ghosts at same seeding density. To us, this difference is presumably to be ascribed to the difference in relative thickness of the two volumes. While the volume of the previous test case had a 6mm thickness, to be compared with the 80mm and 50mm X and Y extensions, the volume of the First LPT Challenge case is 30mm thick, with 100mm and 50mm extensions in the X and Y directions. Lines of sight projecting to a given camera pixel thus have a larger portion within the volume in the latter case, making the formation of ghosts more probable. This difference in particle reconstruction performance explains why we restrict to the Challenge cases

at ppp ≤ 0.08: indeed, at this seeding density, the Particle Precision is comparable to its value reached for 0.14 ppp in the previous case, where we saw the performances of F-VFC starting to drop. Another reason for this choice is the (expected) limited accuracy of the three-point Gaussian fit used after PVR-SMART to obtain the particle locations; we observed that it led to quite important position errors, logically increasing with the seeding density. As this error was found of 0.4 voxel at 0.08 ppp, continuing the tests at higher densities would have led to a drop in performance linked to this error, not to that of F-VFC.

<table><caption>**Table 3** Two-pulse case of the first LPT Challenge: performance metrics for particle reconstruction using PVR-SMART, and for final vector fields obtained with F-VFC. For the latter, Recall is here defined with respect to the ground truth.</caption><thead><tr><th rowspan="2">Seeding (ppp)</th><th colspan="2">Reconstruction</th><th colspan="2">VFC Performance</th></tr><tr><th>Recall</th><th>Precision</th><th>Recall</th><th>Precision</th></tr></thead><tbody><tr><td>0.005</td><td>1</td><td>0.993</td><td>0.995</td><td>0.993</td></tr><tr><td>0.025</td><td>0.998</td><td>0.958</td><td>0.99</td><td>0.974</td></tr><tr><td>0.050</td><td>0.994</td><td>0.751</td><td>0.986</td><td>0.914</td></tr><tr><td>0.08</td><td>0.979</td><td>0.462</td><td>0.888</td><td>0.852</td></tr></tbody></table>

We also write in Table 3 the values reached for Vector Recall and Precision, using the double-pass F-VFC algorithm from these particle reconstructions at the various seeding densities tested, also obtained by automatic evaluation. Note that, contrary to the previous 3D synthetic test case, the Vector Recall is here compared to the ground truth, not to the particle reconstruction. Indeed, the automatic evaluation of the portal only provides the averaged performances (whether on the reconstruction or on the vector fields), and does not indicate which of the particles or vectors are respectively TP and FP. Values reported in Table 3 confirm the observations made with the previous 3D test cases: F-VFC remains robust to very high proportions of ghost particles in particle clouds, maintaining relatively high Recall and Precision at 55% of ghost particles in the dataset. As it is nowadays known that optimal 3D PTV performances at high seeding densities are obtained by algorithms which iterate several times over the operations (reconstruction followed by matching) which we here only perform once, the level of performance that we obtain is thus quite remarkable, and shows that the F-VFC method has a strong potential for leading to even higher performances when integrated in such an iterative approach.

While the present results allow to benchmark F-VFC, together with a particle reconstruction method based on the PVR-SMART paradigm, with state-of-the-art 3D PTV algorithms, their global nature does not allow for direct comparison of F-VFC with other pairing algorithms, as performances also heavily depend on particle reconstruction. We will propose a more direct comparison, with a matching approach using a 3D cross-correlation based predictor, in the following section, on experimental data.
# 5 Experimental test case

In this section, we assess the performances of F-VFC on a von Kármán flow experiment. Pairs of images are extracted from a full time-resolved acquisition, from which 3D particle sets are derived using the Iterative Particle Reconstruction (IPR) (Wieneke, 2012; Jahn et al., 2021). F-VFC is then applied to obtain particle pairings, and assessed with respect to results provided by Shake-the-Box (STB) on the full time-resolved dataset, the latter being considered as ground truth. As just mentioned, in order to better highlight the gain brought by F-VFC compared to state-of-the-art approaches, we will compare its performance with that obtained with a matching approach using a 3D cross-correlation based predictor, such as done for instance by (Agüera et al., 2016; Cornic et al., 2020).

## 5.1 Experimental setup

The turbulent flow is generated inside the Giant von Kármán platform at CEA SPEC (Cheminet et al., 2022; Geneste, 2023). The sheer size of the platform (Radius = 0.5 m, Height = 0.8 m) allows us to optically access the smallest scales of the flow, even reaching sub-Kolmogorov measurements. Inside the Giant von Kármán, the flow is generated between two counter rotating impellers facing each other. The rotating impellers transfer momentum to the fluid by viscous diffusion, the centrifugal force propels the fluid outwards radially. This type of flow has been extensively studied in Debue (2019); Dubrulle (2019). The fluid is water, which we maintain at 20°C with cooling circuits above and below the two turbines. This ensures that the kinematic viscosity $ν$ of the flow is kept constant at $ν = 1.0 \times 10^{-6} \text{ m}^2\text{s}^{-1}$. The rotating frequency is $f = 0.1$ Hz. The radius of the von Kármán is $R= 0.5$ m. The Reynolds number of the flow is $Re = \frac{2πR^2f}{ν} = 1.57 \times 10^5$. The flow is seeded with monodisperse polystyrene beads manufactured by Spherotech. The diameter of the particles is between 5.0 and 5.9 µm, with a mean diameter of 5.33 µm, and of density 1.05. The corresponding Stokes number is $9.1 \times 10^{-5}$, ensuring they behave as passive tracers.

As depicted in Fig. 11, illumination of a $50 \times 44 \times 6 \text{ mm}^3$ volume located at the centre (both in the horizontal and vertical directions) of a longitudinal plane of the von Kármán setup is achieved by using a 30 mJ pulse high speed Nd-YLF laser. The laser has its beam waist located 900 mm from the 8 mm diameter output. A sketch of the optical setup is shown in Fig. 11.

The illuminated volume is imaged by a system of four Phantom Miro M340 cameras, equipped with CMOS sensors of $2560 \times 1600$ pixels with $10 \mu m$ pitch. Cameras 1 and 4 (resp. 2 and 3) are placed at a $\pm 67.5^\circ$ (resp. $\pm 45^\circ$) angle relative to the laser. This ensures forward light scattering, although perceived intensity varies widely depending on the angle of observation. ND8 Filters were thus positioned in front of Cameras 2 and 3 to homogenize the average intensity on all four cameras.

Each sequence of measurement is made of 3226 consecutive time-steps acquired in a time-resolved way, with an acquisition frequency of $f_{acq} = 1.2$ kHz. The acquisition frequency

<figure><img src="image_13.png" /><figcaption>**Fig. 11 (Top)** The Giant von Kármán facility (Top left) and experimental set-up (Top right) at CEA SPEC. The four cameras observe a 50 mm x 45 mm x 6 mm volume at the center of the von Kármán. (Bottom) Sketch of the laser optics used to create the illuminated volume: focal lengths are indicated in red and distances in green, both in mm.</figcaption></figure>

was tuned so as to get an average image particle displacement of 5 pixels. Calibration was done using a LaVision 3D double-sided dual-plane calibration plate.

## 5.2 Data processing

As our purpose here is to illustrate the performance of F-VFC on experimental data and to compare it to that of a state-of-the-art pairing approach, a ground truth is necessary. It is obtained by processing three full time-resolved sequences, with seeding densities estimated at around 0.02, 0.05 and 0.08 ppp on cameras 2 and 3 (see Fig. 11), with the Shake-The-Box (STB) algorithm (LaVision DaVis 10.2 version). For this processing, after calibration and volume self-calibration, images are pre-processed with a spatial filter that removes the local sliding minimum (computed over a square of side 11 pixels), in order to reduce background noise, before running STB with standard parameters recommended in Davis User Manual.

Within each sequence / ppp, we then single out seven different pairs of neighboring instants, from which we reconstruct individual particle clouds, using the Iterative Particle Reconstruction (IPR) algorithm from Davis 10.2. These two particle clouds form the input data on which F-VFC and the 3D cross-correlation based matching are performed. As it is known that the global STB approach, in which reconstruction and matching are coupled, leads to a highly enhanced particle reconstruction by IPR than when performed in an isolated way, it is expected that the particle clouds taken as input for testing the pairing methods will be characterized by a significant amount of ghosts. This is indeed the case: taking the particles reconstructed by STB as ground truth, we find that the individual IPR reconstructions are characterized by the following values for particle Precision:
* ppp = 0.02: #STB tracks = 29007, Particle Precision = 0.781

* ppp = 0.05: #STB tracks = 79517, Particle Precision = 0.646

* ppp = 0.08: #STB tracks = 124296, Particle Precision = 0.336

The reconstructed particle clouds are indeed rather noisy, with over 65% of ghost particles on average at 0.08 ppp.

The particle clouds are then processed with F-VFC similarly to the synthetic cases. All possibilities match in a search radius $R_S$ around every detection are considered to create the initial dataset. The search radius is set larger than the maximal displacement found in STB, so that all possible correct vectors are in the initial vector sets. Then, we run the double-pass F-VFC algorithm with following parameters:

* Regularization: $\lambda = 1$ and $\beta = \frac{1}{2 \cdot (10 \langle d \rangle)^2}$ where $\langle d \rangle$ is the mean inter-particular distance in the input particle clouds.

* Stopping criterion: stop the first pass when $\delta E$ is maximal, then set $N_{lim} = 10$ in the second pass.

* Outlier selection: $\Omega = 4/3\pi R_S^3$ with $R_S = 0.45$ mm = 16 voxels, $\tau = 0.1$ on the first pass, and $\tau = 0.5$ on the second pass.

* Initialization: $\sigma_{init}$ set as obtained from equation 12, $\gamma_{init} = 0.5$.

* Number of subdivisions: $N_S$ tuned at each seeding density so as to have at least around 100 inliers per subvolume.

In the current implementation, which is non parallelized and on CPU (with matrix inversion in the E-M steps done using functions from Python Jax library), the corresponding processing times along with the amount of vectors in the initial vector set created, are as follows:

* ppp = 0.02: 95.4 s for ~ 70000 vectors

* ppp = 0.05: 287.1 s for ~ 220000 vectors

* ppp = 0.08: 1690.9 s for ~ 990000 vectors

Since the current implementation has a strong potential for parallelization (see section 2.4), such processing times thus indicate that F-VFC's performances can be obtained at reasonable cost, and that the approach is thus totally suited for practical use.

The matching approach based on 3D cross-correlation to which we compare F-VFC is that adopted within the Double-Frame Tomographic PTV algorithm, introduced by Cornic et al. (2020). Based on the input particle clouds provided by IPR, a displacement predictor is obtained by extending the particles on a coarse voxel grid (with a voxel-to-pixel ratio $v/p = 3$), and performing 3D cross-correlation using the FOLKI-3D algorithm Cheminet et al. (2014), with cubic interrogation volumes of 65 voxels side. Similar to Cornic

<figure><img src="image_14.png" /><figcaption>**Fig. 12 Solid lines:** Vector Recall (●) (relative to the reconstruction) and Vector Precision (×) obtained at each seeding density with F-VFC. **Dashed lines:** Vector Recall (■) (relative to the reconstruction) and Vector Precision (★) obtained at each seeding density with 3D correlation based matching. Results have been averaged over the seven different time steps tested.</figcaption></figure>

et al. (2020), we observed that this quite large value of $v/p = 3$, which had been singled out to maintain reasonable processing times, does not degrade performance; we indeed found similar metrics when using a finer, $v/p = 1$ grid. Pairing is then done using the obtained predictor and a search radius of 2 voxels around the predicted position.

## 5.3 Results

We show in Fig. 12 the Vector Recall (relative to the reconstructed particle clouds) and Precision, averaged over the seven snapshots, for each seeding density. We find that, whereas Recall drops steadily with correlation-based matching, F-VFC recovers practically every True Positive in the reconstructed particle cloud, even at the highest density, as the Recall remains over 99%. This is not the case however for the Precision, which is around 90% at low and median densities, and even falls below 70% for the densest cases. It is interesting to notice that F-VFC and correlation-based matching maintain similar Precision values at the different seeding densities. This again raises the question of the characteristics of remaining outliers; to get a better grasp on the matter, we show in Fig. 13 the linear interpolation of the velocity, in the center plane of the volume, obtained by STB, F-VFC and correlation-based matching. Fig. 14 also shows the linear interpolation of the error on the velocity magnitude with respect to STB results. We find that at low and median densities, all approaches reconstruct a quasi identical flow: there are close to no distinguishable differences between their results, although correlation-based matching does start exhibiting some outliers at median density - this trend increasing further at high seeding density, where a large amount of outliers is observed. At this highest seeding density, although F-VFC results have increased levels of
<figure><img src="image_15.png" /><figcaption>**Fig. 13** Linear interpolation of the velocity on the center plane of the volume. Top: STB reconstruction, center: F-VFC, bottom: matching with a 3D cross-correlation based predictor. From left to right: 0.02 ppp, 0.05 ppp, 0.08 ppp.</figcaption></figure>

noise, it can be seen that the field is much smoother than the one obtained with the correlation-based matching approach. Furthermore, the amount of noise in F-VFC results is not substantial enough to be in line with the low Precision we have found, of around 70%. A first possible explanation could be that, here again, most of these outliers would be close to the Ground Truth, i.e. fall within the probability distribution of inliers. While this might be the case of some of them, we believe that, possibly, the quite high number of remaining outliers could be simply ascribed to the method used in our performance test, more precisely due to taking STB results as the ground truth. Indeed, looking at the total number of vectors (inliers and outliers) found in each case (see Fig. 15), we observe that F-VFC is very close to

reconstructing the correct amount of vectors. Knowing that STB vector fields result from a temporal fit, which leads to particles positions slightly shifted from their raw positions in the IPR reconstructions, it could be thus probable that our actual Vector Precision is in fact higher than what the present test method returns.

# 6 Conclusion

In 3D PTV, after the reconstruction step, particle clouds at both time steps are riddled with ghost particles, which make identifying the correct displacement vectors complicated. We adapted and built on the Vector Field Consensus (VFC) (Ma et al., 2014) algorithm to this purpose, leading
<figure><img src="image_16.png" /><figcaption>**Fig. 14** Linear interpolation of the velocity error on the center plane of the volume. Top: F-VFC Results; Bottom: matching based on 3D cross-correlation. From left to right: 0.02 ppp, 0.05 ppp, 0.08 ppp.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15** Average number of pairing vectors found in the volume, for ppp tested. From left to right: 0.02 ppp, 0.05 ppp, 0.08 ppp.</figcaption></figure>

to a novel algorithm which we termed Flow Vector Field Consensus, F-VFC. Both VFC and F-VFC algorithms are controlled by a seemingly large number of parameters: we have however shown that all of their optimal values can be easily and physically predicted from the input dataset. To this extent, provided the particle concentration can be estimated in the volume, the F-VFC algorithm is not subject to mis-estimations of user parameters. Furthermore, by choosing the amount of E-M iterations, the operator can have control on the output, whether Precision or Recall needs to be favored. We proposed a double-pass implementation

that aims at maximizing Recall, while strongly increasing processing speeds: its performances can be pinned on a well chosen stopping criterion for the first pass.

The F-VFC algorithm was evaluated on synthetic 3D test cases generated from the John Hopkins Turbulent Database and from the First LPT Challenge. Our findings are that it is robust enough to deal with particle clouds with up to 45% of ghost particles, and returns correct displacement fields with a Precision over 80% and Recall over 90%. The Precision values reached can seem problematic: luckily, even for the noisiest test cases, we find that a significant portion (at least 65%) of
the False Positives are in agreement with the Ground Truth. Despite the False Positives, the flow is being reconstructed with extremely high fidelity. Also, these performances were obtained by performing a single step of matching, whereas full state-of-the-art 3D PTV algorithms also leverage on iterating several times over the reconstruction and matching phases, exploiting the residual images.

We also challenged our algorithm by comparing it to results obtained from Shake-The-Box on images acquired in our GVK Platform at CEA SPEC. We have found that F-VFC reconstructed close to every correct vector, maintaining a recall above 99% at the highest seeding densities tested, relative to vectors in the reconstructed dataset. Other reconstructed vectors, despite being identified as outliers, mostly remain coherent with the flow and therefore do not affect the global quality of the reconstructed flow field. We also found that, at the highest seeding density tested, F-VFC recovered more vector than a cross-correlation based approach, while maintaining similar if not lower levels of noise.

The high performances we obtained with the double-pass F-VFC algorithm lead us to believe that it can be a useful tool in particle tracking algorithms. The high values of Recall can be an asset when initializing particle trajectories for a time-resolved implementation. Then, the high fidelity with which the flow is being reconstructed (with the few remaining outliers being largely in agreement with the Ground Truth) ensures that F-VFC can provide a great predictor to reconstruction/tracking algorithms at high seeding densities.

## Conflicts of interest

The authors acknowledge that there are no known competing financial interests or personal relationships that have influenced the work presented in this article.

## Acknowledgments

This work was funded through ANR EXPLOIT, grant agreement no. ANR-16-CE06-0006-01 and ANR BANG grant agreement no. ANR-22-CE30-0025. The authors would like to thank Pierre Bragança, Christophe Cuvier and Jean-Marc Foucaut from Laboratoire de Mécanique des Fluides de Lille (LMFL), as well as Cécile Wiertel-Gasquet, Damien Geneste, Vincent Padilla, Antoine Barlet, Bérengère Dubrulle and François Daviaud from CEA Paris-Saclay, CNRS, for their help and expertise throughout the several months of experimental measurements on the GVK Platform. Special thanks to Jean-Philippe Laval from LMFL, for providing us with the Direct Numerical Simulation results used for the 2D test cases.

# A Expectation Maximization Algorithm

## A.1 Basic principle

Assume a statistical model $P(V|\Theta)$ that is underlied by some latent variables. These variables are unobserved but their

knowledge would simplify the estimation of $\Theta$. The EM algorithm (Dempster et al., 1977) is a general framework that copes with latent variables in order to maximize the likelihood of $\Theta$ or in a Bayesian context, the posterior likelihood $P(\Theta|V)$.

Let $z_n \in \mathcal{Z}$ be the latent variable associated to $v_n \in \mathcal{V}$. Based on a current value $\Theta^k$ of the parameter, one can define

$$p_n^k = P(z_n|\mathcal{V}, \Theta^k). \qquad (13)$$

Then, one defines the expectation of the complete log-likelihood

$$Q(\Theta, \Theta^k) = E_{p^k}(\ln P(\mathcal{V}, \mathcal{Z}, \Theta)). \qquad (14)$$

Then using Bayes rule,

$$\frac{P(\Theta|\mathcal{V})}{P(\Theta^k|\mathcal{V})} = \frac{P(\mathcal{V}, \mathcal{Z}, \Theta)}{P(\mathcal{V}, \mathcal{Z}, \Theta^k)} \frac{P(\mathcal{Z}|\mathcal{V}, \Theta^k)}{P(\mathcal{Z}|\mathcal{V}, \Theta)}. \qquad (15)$$

Taking the expectation over $P(\mathcal{Z}|\mathcal{V}, \Theta^k)$ of the log of this expression

$$\ln P(\Theta|\mathcal{V}) - \ln P(\Theta^k|\mathcal{V}) = Q(\Theta, \Theta^k) - Q(\Theta^k, \Theta^k) + KL(P(\mathcal{Z}|\mathcal{V}, \Theta^k)||P(\mathcal{Z}|\mathcal{V}, \Theta)), \qquad (16)$$

where $D_{KL}(P||Q)$ is the Kullback-Leibler divergence of P from Q: it is non negative whatever P and Q, and is equal to 0 when the two distributions are identical. $D_{KL}(P||Q)$ can therefore be interpreted as a discrepancy measure between distributions P and Q. From (16), it is clear that increasing $Q(\Theta, \Theta^k)$ increases $\ln P(\Theta|\mathcal{V})$. Therefore the E-M algorithm is defined as follows

$$\begin{array}{ll}
\textit{Expectation (E)} & \text{Compute } Q(\Theta, \Theta^k) \text{ from (13)(14)} \\
\textit{Maximization (M)} & \Theta^{k+1} = \arg \max_{\Theta} Q(\Theta, \Theta^k).
\end{array}$$

It therefore means that an E-M step increases the posterior likelihood. Under some appropriate regularity properties, the EM algorithm converges to the MAP estimator of $\Theta$. These two steps are developed hereafter in our particular setting.

## A.2 Application for the VFC Algorithm

In this paper, the observed variables are the sets $\mathcal{X}$, $\mathcal{V}$, the parameter $\Theta$ is the vector $(\bar{u}, \sigma, \gamma)$, and the latent variable set associated to $\mathcal{V}$ is $\mathcal{Z}$, where the $z_n$ are either equal to 0 or 1. Note that we do not need to define a statistical model for $\mathcal{X}$, but specification of $\mathcal{X}$ is needed to define $p(\bar{u})$, so all the quantities will be conditioned by $\mathcal{X}$ to recall this dependence, but this does not impact the logic of the E-M algorithm defined in section A.1.

We assume flat priors on $\sigma$, $\gamma$. Also, we assume $p(\bar{u}) \propto e^{-\frac{1}{2}\phi(\bar{u})}$ where $\phi$ is a smoothing filter and $\lambda$ a positive real number. Hence $P(\Theta) \propto e^{-\frac{1}{2}\phi(\bar{u})}$ where the proportionality coefficients do not depend on $\Theta$.

Looking for a MAP Estimate of $\Theta$, we will use the E-M Algorithm of section A.1. Taking a current estimate $\Theta^k$:
### A.2.1 E-step

From (13), the update of distribution $p$ of the $z_n$ writes

$$ \begin{cases} p_n^k = \mathbb{P}(z_n = 1|x_n, v_n, \Theta^k) \\ 1 - p_n^k = \mathbb{P}(z_n = 0|x_n, v_n, \Theta^k) \end{cases} $$

To derive these quantities, we use Bayes formula, writing

$$ \mathbb{P}(z_n = 1|x_n, v_n, \Theta^k) = \frac{\mathbb{P}(v_n, z_n = 1|x_n, \Theta^k)}{\mathbb{P}(v_n|x_n, \Theta^k)} $$

and the law of total probability, writing

$$ \mathbb{P}(v_n|x_n, \Theta^k) = \mathbb{P}(v_n, z_n = 1|x_n, \Theta^k) + \mathbb{P}(v_n, z_n = 0|x_n, \Theta^k) $$

Finally, we obtain (5) from

$$ \begin{cases} \mathbb{P}(v_n, z_n = 1|x_n, \Theta) = \frac{\gamma}{(2\pi\sigma^2)^{D/2}} e^{-\frac{\|v_n - \bar{u}(x_n)\|^2}{2\sigma^2}} \\ \mathbb{P}(v_n, z_n = 0|x_n, \Theta) = \frac{1-\gamma}{\Omega} \end{cases} $$

### A.2.2 M-step

From (14), we find that the maximization problem to update $\Theta$ is

$$ \Theta = \arg\max_{\Theta} [\mathbb{E}_{p^k}(\ln \mathbb{P}(\mathcal{Y}, \mathcal{Z}, \Theta | \mathcal{X}))] \quad (17) $$

For maximizing with respect to $\Theta$, $\ln \mathbb{P}(\Theta)$ reduces down to $\ln p(\bar{u})$. The quantity to maximize with respect to $\Theta$ becomes

$$ \sum_n p_n^k \ln \mathbb{P}(v_n, z_n = 1|x_n, \Theta) + \sum_n (1 - p_n^k) \ln \mathbb{P}(v_n, z_n = 0|x_n, \Theta) + \ln p(\bar{u}) $$

Taking the values for $p$ found during the E-step, we eventually find (4). Note that we also discarded the term $-\ln(\Omega) \sum_n \mathbb{P}(z_n = 0|x_n, v_n, \Theta)$, as it does not depend on $\Theta$.

## References

Nereida Agüera, Gioacchino Cafiero, Tommaso Astarita, and Stefano Discetti. Ensemble 3D PTV for high resolution turbulent statistics. *Measurement Science and Technology*, 27(12): 124011, oct 2016. doi: 10.1088/0957-0233/27/12/124011. URL https://dx.doi.org/10.1088/0957-0233/27/12/124011.

Callum Atkinson and Julio Soria. An efficient simultaneous reconstruction technique for tomographic particle image velocimetry. *Experiments in fluids*, 47(4):553–568, 2009.

S.J. Baek and S.J. Lee. A new two-frame particle tracking algorithm using match probability. *Experiments in Fluids*, 22: 23–32, 1996. doi: https://doi.org/10.1007/BF01893303.

Théo Benkovic, Jean-François Krawczynski, and Philippe Druault. Vision-based correspondence using relaxation algorithms for particle tracking velocimetry. *Measurement Science and Technology*, 2020. doi: 10.1088/1361-6501/abb437. URL https://hal.sorbonne-universite.fr/hal-02947162.

Christopher Bishop. *Pattern Recognition and Machine Learning*, volume 16, pages 140–155. 01 2006. doi: 10.1117/1.2819119.

W. Brevis, Y. Niño, and G.H. Jirka. Integrating cross-correlation and relaxation algorithms for particle tracking velocimetry. *Experiments in Fluids*, 50:135–147, 2011. URL https://doi.org/10.1007/s00348-010-0907-z.

Frederic Champagnat, Philippe Cornic, Adam Cheminet, Benjamin Leclaire, Guy Le Besnerais, and Aurelien Plyer. Tomographic PIV: Particles versus blobs. *Measurement Science and Technology*, 25:084002, 07 2014. doi: 10.1088/0957-0233/25/8/084002.

A. Cheminet, B. Leclaire, F. Champagnat, A. Plyer, R. Yega-vian, and G. Le Besnerais. Accuracy assessment of a Lucas-Kanade based correlation method for 3D PIV. In *17th International Symposium on Applications of Laser Techniques to Fluid Mechanics*, Lisbonne, Portugal, July 2014. URL https://hal.science/hal-01100782.

Adam Cheminet, Benjamin Leclaire, Frédéric Champagnat, Philippe Cornic, and Guy Le Besnerais. On factors affecting the quality of tomographic reconstruction. *In Proceedings of PIV13*, 2013.

Adam Cheminet, Damien Geneste, Antoine Barlet, Yasar Osto-van, Tarek Chaabo, Valentina Valori, Paul Debue, Christophe Cuvier, François Daviaud, Jean-Marc Foucaut, J. Laval, Vincent Padilla, Cécile Wiertel, and Berengere Dubrulle. Eulerian vs Lagrangian irreversibility in an experimental turbulent swirling flow. *Physical Review Letters*, 129, 09 2022. doi: 10.1103/PhysRevLett.129.124501.

Philippe Cornic, Benjamin Leclaire, Frédéric Champagnat, Guy Le Besnerais, Adam Cheminet, Cédric Illoul, and Gilles Losfeld. Double-frame tomographic PTV at high seeding densities. *Experiments in Fluids*, 61:1–24, 2020.

Philippe Cornic, Frédéric Champagnat, and Benjamin Leclaire. Lagrangian particle tracking: a link between localisation error and fraction of missed particles. In *ISPIV 2021*, 2021.

Paul Debue. *Experimental approach to the problem of the Navier-Stokes singularities*. Theses, Université Paris Saclay (CO-MUE), September 2019. URL https://theses.hal.science/tel-02420454.

A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the EM algorithm. *Journal of the Royal Statistical Society B*, 39:1–38, 1977.

Stefano Discetti and Tommaso Astarita. Fast 3d PIV with direct sparse cross-correlation. *Experiments in fluids*, 2012.

Bérengère Dubrulle. Beyond Kolmogorov cascades. *Journal of Fluid Mechanics*, 867:P1, 2019. doi: 10.1017/jfm.2019.98.

Thomas Fuchs, Rainer Hain, and Christian Kähler. Non-iterative double-frame 2D/3D particle tracking velocimetry. *Experiments in Fluids*, 58:119, 08 2017. doi: 10.1007/s00348-017-2404-0.

Damien Geneste. *Des évènements extrêmes de dissipation aux évènements extrêmes de dispersion, une empreinte des singularités : Approche expérimentale*. Theses, Université Paris-Saclay, February 2023. URL https://theses.hal.science/tel-04046970.
Philipp Godbersen, Daniel Schanz, and Andreas Schröder. Peak-CNN: improved particle image localization using single-stage CNNs. *Experiments in Fluids*, 65:1–19, 10 2024. doi: 10.1007/s00348-024-03884-z.

M. Ishikawa, Y. Murai, and A. Wada. A novel algorithm for particle tracking velocimetry using the velocity gradient tensor. *Experiments in Fluids*, 29:519–531, 2000. URL https://doi.org/10.1007/s003480000120.

Tobias Jahn, Daniel Schanz, and Andreas Schröder. Advanced iterative particle reconstruction for lagrangian particle tracking. *Experiments in Fluids*, 62, 08 2021. doi: 10.1007/s00348-021-03276-7.

P Jia, Y Wang, and Y Zhang. Improvement in the independence of relaxation method-based particle tracking velocimetry. *Measurement Science and Technology*, 24(5):055301, mar 2013. doi: 10.1088/0957-0233/24/5/055301. URL https://dx.doi.org/10.1088/0957-0233/24/5/055301.

Benjamin Leclaire, Ivan Mary, Cedric Liauzun, Stéphanie Péron, Andrea Sciacchitano, Andreas Schröder, Philippe Cornic, and Frederic Champagnat. First challenge on lagrangian particle tracking and data assimilation: datasets description and planned evolution to an open online benchmark. In *Proceedings of 14th International Symposium on Particle Image Velocimetry*, volume 1. ILLINOIS Tech/Paul V. Galvin Library., 08 2021. doi: 10.18409/ispiv.v1i1.119.

Benjamin Leclaire, Philippe Cornic, Frédéric Champagnat, Eric Fabre, and Frédéric Calmels. A web portal for automatic performance evaluation of lagrangian particle tracking and data assimilation algorithms. In *LXLaser 2022-20th International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics*, 2022.

Y. C. Lei, W. H. Tien, J. Duncan, M. Paul, N. Ponchaut, C. Mouton, D. Dabiri, T. Rösgen, and J. Hove. A vision-based hybrid particle tracking velocimetry (PTV) technique using a modified cascade correlation peak-finding method. *Experiments in Fluids*, 53(5):1251–1268, November 2012. doi: 10.1007/s00348-012-1357-6.

Yi Li, Eric Perlman, Minping Wan, Yunke Yang, Charles Menneveau, Randal Burns, Shiyi Chen, Alexander Szalay, and Gregory Eyink. A public turbulence database cluster and applications to study lagrangian evolution of velocity increments in turbulence. *Journal of Turbulence*, 9, 04 2008. doi: 10.1080/14685240802376389.

Jiaming Liang, Shengze Cai, Chao Xu, Tehuan Chen, and Jian Chu. DeepPTV: Particle tracking velocimetry for complex flow motion via deep neural networks. *IEEE Transactions on Instrumentation and Measurement*, 71:1–16, 2022. doi: 10.1109/TIM.2021.3120127.

Jiayi Ma, Ji Zhao, Jinwen Tian, Alan Yuille, and Z. Tu. Robust point matching via vector field consensus. *Image Processing, IEEE Transactions on*, 23:1706–1721, 04 2014. doi: 10.1109/TIP.2014.2307478.

Bertrand Mercier, Lionel Thomas, Benoit Tremblais, and Laurent David. A robust pairing method for two-pulse particle tracking velocimetry based on coherent point drift. *Measurement Science and Technology*, 35(6):065301, mar 2024. doi: 10.1088/1361-6501/ad2b42. URL https://dx.doi.org/10.1088/1361-6501/ad2b42.

A V Mikheev and V M Zubtsov. Enhanced particle-tracking velocimetry (EPTV) with a combined two-component pair-matching algorithm. *Measurement Science and Technology*, 19(8):085401, jun 2008. doi: 10.1088/0957-0233/19/8/085401. URL https://dx.doi.org/10.1088/0957-0233/19/8/085401.

Andriy Myronenko and Xubo Song. Point set registration: Coherent point drift. *IEEE transactions on pattern analysis and machine intelligence*, 32:2262–75, 12 2010. doi: 10.1109/TPAMI.2010.46.

Matteo Novara, Daniel Schanz, Sebastian Gesemann, K Lynch, and Andreas Schröder. Lagrangian 3d particle tracking for multi-pulse systems: performance assessment and application of shake-the-box. In *Proceedings of 18th International Symposium on applications of laser imaging techniques to fluid mechanics*, 2016.

Matteo Novara, Daniel Schanz, and Andreas Schröder. Two-pulse 3d particle tracking with shake-the-box. *Experiments in Fluids*, 64, 04 2023. doi: 10.1007/s00348-023-03634-7.

Francisco Pereira, Heinrich Stüer, Emilio C Graff, and Morteza Gharib. Two-frame 3d particle tracking. *Measurement Science and Technology*, 17(7):1680, jun 2006. doi: 10.1088/0957-0233/17/7/006. URL https://dx.doi.org/10.1088/0957-0233/17/7/006.

Eric Perlman, Randal Burns, Yi Li, and Charles Meneveau. Data exploration of turbulence simulations using a database cluster. *Proceedings of the 2007 ACM/IEEE Conference on Supercomputing, SC'07*, page 23, 11 2007. doi: 10.1145/1362622.1362654.

Fulvio Scarano. Tomographic PIV: Principles and practice. *Measurement Science and Technology*, 24:012001, 10 2012. doi: 10.1088/0957-0233/24/1/012001.

Daniel Schanz, Sebastian Gesemann, and Andreas Schröder. Shake-the-box: Lagrangian particle tracking at high particle image densities. *Experiments in Fluids*, 57, 04 2016. doi: 10.1007/s00348-016-2157-1.

Andrea Sciacchitano, Benjamin Leclaire, and Andreas Schröder. Main results of the first lagrangian particle tracking challenge. In *Proceedings of 14th International Symposium on Particle Image Velocimetry*, volume 1. ILLINOIS Tech/Paul V. Galvin Library., 08 2021. doi: 10.18409/ispiv.v1i1.197.

Bernhard Wieneke. Iterative reconstruction of volumetric particle distribution. *Measurement Science and Technology*, 24:024008, 12 2012. doi: 10.1088/0957-0233/24/2/024008.

Ke Zhou, Jiaqi Li, Jiarong Hong, and Samuel J Grauer. Stochastic particle advection velocimetry (spav): theory, simulations, and proof-of-concept experiments. *Measurement Science and Technology*, 34(6):065302, mar 2023. doi: 10.1088/1361-6501/acc049. URL https://dx.doi.org/10.1088/1361-6501/acc049.