General Defocusing Particle Tracking:
Fundamentals and uncertainty assessment

The purpose of a 3D particle tracking system is to locate the physical coordinates

of a number $N_{\mathrm{p}}$ of tracer particles inside a measurement volume $V$ and to track their displacement in time. The tracking problem will not be discussed here; in GDPT analysis, the particle concentration is relatively low and conventional tracking algorithms can be used. The fundamental analysis in this work concerns only the problem of particle location determination in GDPT. As we consider single-camera acquisition, it is convenient to define the reference frame with two in-plane coordinates, $x$ and $y$, perpendicular to the optical axis of the camera objective, and one depth coordinate, $z$, parallel to the optical axis (Fig. 1). The measurement volume can be approximated by a rectangular cuboid with dimensions $l\times w\times h$, being $h$ the dimension in the depth direction. The maximum size in the in-plane direction ($l\times w$) is set by the field of view of the imaging system. The maximum depth $h$ that can be achieved depends on the imaging system, the particle size, and the illumination intensity, and corresponds to the region where the signal of defocused particle images is strong enough to be processed by the image analysis method.

When using single-camera systems, the tracer particles in the measurement volume are recorded on a single image. For a given optical system, we can define a particle image function $I_{\mathrm{c}}(X,Y,x,y,z)$ which provides the image of one particle depending on its position in the physical space. $I_{\mathrm{c}}$ can be seen as a point spread function defined for one particle, where $X$ and $Y$ are the coordinates in the image space, and where we impose that the center of the particle image is located at $(X,Y)=(0,0)$ for any $x$, $y$, and $z$.

We consider an idealized optical system with no distortions, and where $I_{\mathrm{c}}$ is independent of the in-plane coordinates (Fig. 2(a))

Under this approximation, we can represent a general image $I(X,Y)$, containing a total number $N_{\mathrm{p}}$ of particles, as

where $M$ is the magnification of the optical system and where $I_{0}$ is a term that accounts for background intensity and thermal noise of the camera.

With respect to the final result of a GDPT evaluation, two primary parameters must be considered:

• •

  Signal-to-noise ratio (SNR). Following a convention used in image analysis, we define the SNR as the ratio between the mean particle image signal $\mu_{\mathrm{p}}$ and the standard deviation of the noise $\sigma_{I}$

  $\displaystyle\mathrm{SNR}=\frac{\mu_{\mathrm{p}}}{\sigma_{I}}.$

  (4)

  Practically, $\sigma_{I}$ can be estimated as the standard deviation of the image intensity in a region without particles, and $\mu_{\mathrm{p}}$ can be calculated as the average intensity of the particle image minus the average background intensity. We can define the boundary of a particle image by setting a threshold, as shown in Fig. 2(b).

• •

  Particle image density $\bm{N_{\mathrm{S}}}$. This parameter is often given in particles per pixels (ppp). This approach, however, does not consider the particle image size, which is an important factor in defocusing applications with respect to overlapping particles, as shown in Fig. 2(c-d). For GDPT analysis, it is more convenient to consider the ratio between the sum of the particle image areas $A_{\mathrm{p}}^{(i)}$ and full image area $A_{I}$ kahler2012uncertainty

  $\displaystyle N_{\mathrm{S}}=\frac{1}{A_{I}}\sum_{i}A_{\mathrm{p}}^{(i)}\approx N_{\mathrm{p}}\frac{\bar{A}_{\mathrm{p}}}{A_{I}}.$

  (5)

  Here, $\bar{A}_{\mathrm{p}}$ is the average of the particle image areas, which can be estimated for instance from the calibration stack. One can easily translate the $N_{\mathrm{S}}$ value in ppp, by dividing it with $\bar{A}_{\mathrm{p}}$. It should be noted that this parameter is equivalent to the source density encountered in classical PIV literature adrian1984scattering ; westerweel2000theoretical .

In real optical systems, the approximation in Eq. (2) may not hold as a consequence of optical aberrations. The most common aberrations which are relevant for GDPT applications are:

• •

  Parallax or perspective error. The magnification is not constant across the measurement depth, i.e. objects closer to the lens appear larger on the image. This error is typically small and can be neglected for microscope objective lenses.

• •

  Field curvature. The focal plane is not flat, therefore the measured $z$ position must be corrected depending on the particle in-plane position. This effect is relevant in GDPT applications since a field curvature of a few micrometers can have a significant impact in the measurement.

• •

  Distortion. The particle images are distorted as they move away from the image center. This error is difficult to correct in GDPT analysis based on a single calibration stack but is normally not strong in conventional optical setups.

More details about practical strategies to deal with optical aberrations and bias errors are discussed in Section 4.

The aim of a GDPT processing is to determine the 3D particle positions from the defocused particle images. The in-plane position can be obtained using one of the many approaches developed for 2D PTV maas1993particle ; ohmi2000particle ; ouellette2006quantitative ; raffel2018particle . The out-of-plane component is determined using a reference set of calibration images and must contain the following elements:

1. 1.

   A discrete reference set of calibration images, referred to as the calibration stack. This can be seen as a discrete sampling at known positions of the particle image function $I_{\mathrm{c}}$

   	$\displaystyle I_{\mathrm{c}}^{(k)}(X,Y)=$	$\displaystyle I_{\mathrm{c}}(X,Y,z_{k})+I_{0}(X,Y)$
   		$\displaystyle\quad\mathrm{with}\,\,k=1,2,...,N_{\mathrm{cal}}.$		(6)

   $I_{\mathrm{c}}^{(k)}$ is typically obtained experimentally by taking subsequent images of a reference particle which is displaced at known positions, for instance using a motorized focusing stage barnkob2015general .

2. 2.

   A function or procedure to identify target particle images $I_{\mathrm{t}}^{(i)}$ inside the image $I$. This normally relies on segmentation algorithms that can be applied on raw or filtered/pre-processed images.

3. 3.

   A function or procedure to quantify the similarity between different images. This is used to rank the calibration images in $I_{\mathrm{c}}^{(k)}$ with respect to their similarity to a given target particle image $I_{\mathrm{t}}^{(i)}$.

4. 4.

   A function or procedure to estimate the final depth position $z^{\prime}_{i}$ of the target particle from the similarity values between $I_{\mathrm{t}}^{(i)}$ and $I_{\mathrm{c}}^{(k)}$. It should be noted that a simple identification of the most similar image in the stack would produce a discrete output. Interpolation schemes are needed to obtain a continuous output with a “sub-image” resolution, in analogy with what is typically done in digital PIV evaluations to obtain sub-pixel resolution.

More generally, the determination of the depth position in GDPT can be seen as a supervised learning problem, where the calibration stack is the training set used to calibrate a prediction algorithm. The scheme of the prediction algorithm remains the same but it can be trained on different setups just using different calibration stacks. We used in this paper the implementation based on the normalized cross-correlation proposed in Ref. 3, but other approaches using neural networks or more sophisticated classification schemes would be natural improvements of the method.

The result of a GDPT evaluation on a single image, containing a total number $N_{\mathrm{p}}$ of particles, is a set of measured particle coordinates $x^{\prime}_{i}$, $y^{\prime}_{i}$, $z^{\prime}_{i}$, where

where $\epsilon_{z_{i}}$ is the measurement error and the same definition applies for the other coordinates. To assess the performance of a GDPT measurement, the following parameters must be considered:

1. 1.

   The measurement uncertainty in the particle position determination.

2. 2.

   The relative number of measured particles.

3. 3.

   The depth of the measurement volume.

For a fair assessment of a GDPT measurement, all these three parameters should be considered, since they are interconnected among each other. For instance, a stricter validation criterion can reduce the error but at the expenses of a smaller number of measured particles.

• •

  Measurement uncertainty. When the true particle position is known, as in the case of synthetic images, the measurement uncertainty in the determination of one coordinate is normally given in terms of the root-mean-square of the error

  $$\sigma_{z}=\sqrt{\frac{\sum_{i}^{\prime}(z^{\prime}_{i}-z_{i})^{2}}{\sum_{i}^{\prime}1}},$$

  (8)

  where $\sum_{i}^{\prime}$ indicates summation over the measured values and $\sum_{i}^{\prime}1$ is equal to the total number of measured particles $N_{\mathrm{p}}^{\prime}$. Since the particle images have different shapes, a non-uniform uncertainty along $z$ is expected and it is useful to define a local uncertainty

  $$\sigma_{z}^{\delta}(z)=\sqrt{\frac{\sum_{i}^{\prime}(z^{\prime}_{i}-z_{i})^{2}\,\theta^{\delta}(z,z_{i})}{\sum_{i}^{\prime}\theta^{\delta}(z,z_{i})}},$$

  (9)

  with

  $\displaystyle\theta^{\delta}(z,z_{i})=\begin{cases}1&\text{for }\,\,z-\delta<z_{i}<z+\delta\\ 0&\text{otherwise}\end{cases}.$

  (10)

  Eq. (9) indicates the uncertainty associated to a bin centered at $z$ and having a width of $2\delta$. The same applies for the uncertainty of the in-plane coordinates $x$ and $y$. For particle images with rotational symmetry, we can assume $\sigma_{x}=\sigma_{y}$. Note that Eqs. (8) and (9) give the total error, since the true values $z_{i}$ are used. If the uncertainty is calculated around a mean value, as in case of experimental images, the uncertainty does not take into account bias errors. A discussion about how to apply these formulas to experimental images is presented in Section 4.

• •

  Relative number of measured particles. This is defined as the ratio between the number of measured particles and the total number of particles in one image

  $$\phi=\frac{\sum_{i}^{\prime}1}{\sum_{i}1}=\frac{N_{\mathrm{p}}^{\prime}}{N_{\mathrm{p}}}.$$

  (11)

  Also in this case it is useful to define a local $\phi$ defined as

  $$\phi^{\delta}(z)=\frac{\sum_{i}^{\prime}\theta^{\delta}(z,z_{i})}{\sum_{i}\theta^{\delta}(z,z_{i})}.$$

  (12)

  The relative number of measured particles $\phi$ tends to decrease as the particle concentration is increased, due to the more frequent occurrence of overlapping particle images. For a given GDPT implementation, there will be a critical particle image density $N_{\mathrm{S}}^{*}$ above which the number of measured particles $N^{\prime}_{\mathrm{p}}$ starts to decrease. This sets the maximum possible seeding density that should be used for that implementation. Furthermore, it should be noted that $\phi$ is not directly accessible in experimental images, where only the number of measured particles $N_{\mathrm{p}}^{\prime}$ is available.

• •

  Measurement volume. The choice of the depth $h$ of the measurement volume affects the total number of particles $N_{\mathrm{p}}$ to be measured and, consequently, the number of overlapping particles and the measurement uncertainty. For instance, a large depth $h$ can include highly-defocused particle images, which have a smaller SNR and are more difficult to detect. To fairly estimate the uncertainty of a given implementation, the measurement depth on which it is applied must be indicated. In GDPT systems, $h$ is practically set by the lowest and highest $z$ coordinate in the stack.

The calibration stack $I_{\mathrm{c}}^{(k)}$ represents a discrete sampling of the particle image function across $z$. The depth position of a particle is obtained by comparing a target particle image $I_{\mathrm{t}}^{(i)}$ with the images in the stack using a similarity function $S$, usually normalized between 0 and 1, with 1 being the perfect match (Fig. 3(a)). In an ideal case, the measurement resolution can be increased indefinitely by increasing the number of images $N_{\mathrm{cal}}$ in the stack (Fig. 3(c)). In a real case, however, beyond a certain $N_{\mathrm{cal}}$ the difference in shape between two neighbor images will be smaller than the difference induced by the image noise and in-plane image discretization (i.e. the light-intensity distribution is discretized into pixels).

This concept is illustrated in Fig. 3(b), where the self-similarity, i.e. the average $S$ between particle images at the same height but different in-plane positions, is shown as a function of $z$. Even in the case with no noise (blue line), $S$ is not unity due to tiny differences induced by sub-pixel displacements of the particle images. By adding noise (red line, $\sigma_{I}=25$), $S$ decreases significantly, especially in regions where the SNR is lower. A way to check whether a calibration stack is over-sampling $I_{\mathrm{c}}$, is to calculate the similarity between neighbor images in the stack, i.e., $S(I_{\mathrm{c}}^{(k)},I_{\mathrm{c}}^{(k+1)})$. In the example in Fig. 3(d), the stack with $N_{\mathrm{cal}}=500$ (orange line) is clearly over-sampled, since it coincides with the self-similarity curve in Fig. 3(b). Further refining the stack will not bring any advantage. For $N_{\mathrm{cal}}=15$ (blue line), the $S$ between neighbor images is small, showing that the stack under-samples the particle image function. The stack with $N_{\mathrm{cal}}=50$ (red line) is a better compromise, with a slightly under-sampled region for $z/h<0.5$ and over-sampled above.

In this section we evaluate the effect, on the depth coordinate determination, of the number of calibration images and sub-image interpolation along with the influence of image noise. More specifically, we analyze Dataset I (Fig. 8(a) in Appendix A) where particle images do not overlap; thus all the particles are detected, i.e., $N_{\mathrm{p}}^{\prime}=N_{\mathrm{p}}$.

To isolate the different effects, we first analyze images with no noise, i.e. $\sigma_{I}=0$. The results are shown in Fig. 4(a), where we plot the local depth coordinate uncertainty $\sigma_{z}^{\delta}(z)$ as a function of $z$ along with a selection of the corresponding particle image shapes. The results are shown when using $N_{\mathrm{cal}}=15$, 50, or 500 calibration images with (points) and without (open circles) the use of sub-image interpolation. It may be noted that the defocused particle images have a sharp ring shape below the focal position (low $z$) and blurred Gaussian-like shape above (high $z$). This is due to the spherical aberration which is present in many optical systems and included in MicroSIG rossi2019synthetic . Detecting the shape differences for blurred particle is more difficult, which explains the larger uncertainty in regions with $z/h>0.6$.

Generally, increasing $N_{\mathrm{cal}}$ leads to a decreasing local uncertainty $\sigma_{z}^{\delta}$, but e.g. for $z/h\sim 0.6$, the use of a sub-image detection scheme can lead to better results for lower $N_{\mathrm{cal}}$ which is seen by a better performance for $N_{\mathrm{cal}}=20$ (blue lines) than for $N_{\mathrm{cal}}=50$ (black lines) or $N_{\mathrm{cal}}=500$ (red lines). The use of a sub-image interpolation scheme allows for a continuous $z$ determination and does in general yield lower uncertainties up to a certain value of $N_{\mathrm{cal}}$. This is shown in Fig. 4(b), where we show the depth coordinate uncertainty $\sigma_{z}$ as a function of $N_{\mathrm{cal}}$. As $N_{\mathrm{cal}}$ approaches 200, the uncertainty, with and without the use of a sub-image scheme, converges to the same value.

We now fix the number of images in the calibration stack to $N_{\mathrm{cal}}=50$ and analyze images with increasing noise levels. The results are presented in Fig. 4(c). Not surprisingly, the local depth coordinate uncertainty $\sigma_{z}^{\delta}$ increases with a decreasing SNR. As explained before, the local uncertainty for larger values of $z$ is larger due to the blurred Gaussian-like shape of the particle images. The effect of noise on the local uncertainty naturally propagates into the total uncertainty as seen in Fig. 4(d). Clearly, as the noise level increases, a larger $N_{\mathrm{cal}}$ is needed in order to reach a converging uncertainty $\sigma_{z}$.

Lastly, to explore the reduction of noise through image pre-processing, we now fix the noise level to $\sigma_{I}=50$ and apply different median filters to the images. We illustrate this in Fig. 4(e) where we apply a median filter of $d_{\mathrm{med}}\times d_{\mathrm{med}}$ (with $d_{\mathrm{med}}$ equal to 5 or 7) to both the calibration and the measurement images. The median filter decreases the impact of noise in the local uncertainty $\sigma_{z}^{\delta}$, in some region even as much as one order of magnitude. Evidently, the filter size has an optimum between averaging out noise and actual particle image features — this is seen by the presence of some slightly lower values of $\sigma_{z}^{\delta}$ for $d_{\mathrm{med}}=5$ than for $d_{\mathrm{med}}=7$. This is confirmed in Fig. 4(f), where the same trend is seen for the average uncertainty $\sigma_{z}$.

After investigating the effect of noise and $N_{\mathrm{cal}}$, we now analyze the more experimentally-realistic case where the presence of overlapping particle images is possible. In particular we use the noise-less ($\sigma_{I}=0$) Dataset II (Fig. 8(b)) to assess the GDPT uncertainties as a function of increasing particle concentration while keeping the number of calibration images fixed ($N_{\mathrm{cal}}=50$). We consider here particle image densities up to $N_{\mathrm{S}}=0.59$, corresponding to 0.001 particles per pixels. As the particle concentration increases, the number of overlapping particle images increases concurrently. This leads to lower similarities $C_{\mathrm{m}}$ between target particle images and calibration images, which furthermore leads to larger errors in the determination of the $z$-coordinates. By setting a higher threshold on the accepted values of $C_{\mathrm{m}}$, poorly-determined $z$-coordinates can be removed to increase the determination accuracy but with the cost of lowering the number of measured particles $N^{\prime}_{\mathrm{p}}$.

In Figs. 5(a), 5(c), and 5(e), we show the local uncertainties $\sigma_{z}^{\delta}$ and $\sigma_{x,y}^{\delta}$, and the local relative number of measured particles $\phi^{\delta}$ for two particle image densities of $N_{\mathrm{S}}=0.05$ (squares) and $N_{\mathrm{S}}=0.59$ (circles), respectively. As $N_{\mathrm{S}}$ is increased from 0.05 to 0.59, the local errors $\sigma_{z}^{\delta}$ and $\sigma_{x,y}^{\delta}$ are consistently increasing, while the local relative number of measured particles $\phi^{\delta}$ is consistently decreasing. As the threshold on $C_{\mathrm{m}}$ is increased, both $\sigma_{z}^{\delta}$, $\sigma_{x,y}^{\delta}$, and $\phi^{\delta}$ are decreasing as expected. The same trend is confirmed in Figs. 5(b), 5(d), and 5(f), where the uncertainties $\sigma_{z}$ and $\sigma_{x,y}$, and the relative number of measured particles $\phi$ are computed as a function of the particle image density $N_{\mathrm{S}}$ and $C_{\mathrm{m}}=0.5$, 0.8, and 0.9, respectively. In addition, we see in Fig. 5(f) that as $N_{\mathrm{S}}$ increases, the measured particle image density $N_{\mathrm{S}}^{\prime}$ (orange colors) reaches a peak, e.g., for $N_{\mathrm{S}}\approx 0.24$ when $C_{\mathrm{m}}=0.9$. In fact, as the particle concentration is increased, the number of particles in the image increases but so does the number of outliers or not-detectable particle images due to the increased number of overlapping particle images that cannot be processed. This can be seen as a critical value, denoted here as $N_{\mathrm{S}}^{*}$, that sets the maximum particle image density that can be used with those settings.

Experimental measurement images are prone to variations in image light intensity, for example, if the image illumination is inhomogeneous over the image plane or if the illumination amplitude changes over time. This can eventually affect the accuracy of defocusing-based particle detection, but depends on the type of image intensity variation and applied detection algorithm. In this section, we assess the GDPT measurement uncertainties when analyzing measurement images with an applied linear light intensity gradient of increasing steepness. More specifically, we use Dataset III (Fig. 8(c)), which is based on the images in Dataset II for $N_{\mathrm{S}}=0.30$ and to which an intensity gradient $\alpha$ has been applied.

Figures 6(a) and 6(c) show the local depth coordinate uncertainty $\sigma_{z}^{\delta}(z)$ and the local relative number of measured particles $\phi^{\delta}(z)$. For small $\alpha$, the changes are small as expected from the fact that the normalized cross-correlation is insensitive to changes in intensity level. However, as $\alpha$ increases, the intensity gradient becomes comparable to the intensities of the particle image features. This affects the uncertainty and the number of measured particles specially in regions where the SNR is low (outside the green area in Fig. 6(a) and (c)). Overall, the effect of intensity gradient is small in the region with large SNR, as shown in Fig. 6(b) and (d). Note also the previously discussed connection between the uncertainty and number of measured particles, e.g., without considering the latter, one would erroneously observe in Fig. 6(b) a better performance for $\alpha=10$ than for $\alpha=5$.

A common way to obtain the calibration stack experimentally is to use one particle at a fixed position (e.g. stuck to a wall) and to move the objective lens at different depth positions. The reference $z$ coordinates are obtained from the reading of the scanning device. During the measurement, however, the situation is reversed: The objective lens is fixed and the particles move. This has consequences if the immersion medium of the lens is different from the medium where the particles are dispersed, and a correction coefficient must be multiplied to the measured $z$. In most applications, this coefficient is equal to the ratio between the refractive index coefficient of the fluid (typically water) and the one of the immersion medium of the lens (typically air) rossi2012effect .

In real applications, the particle image function $I_{\mathrm{c}}$ is not the same across in-plane positions due to optical aberrations. The more straightforward way to address this problem would be to take multiple calibration stacks for different in-plane positions; however, this approach is difficult to implement. First, the size of the calibration stack would increase considerably and so the time necessary for the calibration (a particle must be scanned in different positions). Second, higher complexity and computational costs would be needed to navigate the different calibration stacks and an interpolation scheme must be implemented to cover the entire area of the sensor (only a discrete set of in-plane positions can be taken). Future implementations of the GDPT method using neural network might be able to tackle this task.

In many practically applications, however, the aberrations are weak and a single calibration stack (typically obtained in the center of the image) can be used in the whole image area rossi2019interfacial ; barnkob2018acoustically ; qiu2019experimental . Additional bias errors due to field curvature or perspective errors can be corrected using a reference experiment, for instance, using a fixed array of tracer particles on a flat plane perpendicular to the optical axis and scan it at different depth positions.

In microfluidics applications, where microscope objectives are used, the perspective error is typically negligible, but the field curvature can give errors of few microns cierpka2010calibration , which might be relevant at those scales. If the measurements are performed on a straight microchannel, one practical way to correct this bias is by a direct measurement of a Poiseuille flow inside the duct. In this case, we know that the streamlines in such flow must be straight lines, and we can use this information to determine and correct the field curvature.

To estimate the measurement uncertainty of an experimental GDPT setup, we can use Eqs. (8) and (9), but we need to estimate the unknown true values. A practical approach can be to measure particle displacements in flows having zero velocity in one direction. We can estimate then the displacement uncertainty, considering 0 as the true value. Modeling the error as normally distributed, the positioning uncertainty will be equal to the displacement error divided by $\sqrt{2}$. Of course, particular care should be taken with reducing as much as possible the bias errors and eliminating false positive with an appropriate outlier rejection scheme. Together with the uncertainty, it is important to indicate the depth $h$ of the measurement volume. When it is possible to estimate $N_{\mathrm{S}}$, for instance, if the experimental particle concentration is known, the relative number of measured particles $\phi$ can be calculated; otherwise, $N_{\mathrm{S}}^{\prime}$ should be indicated.

We give an example of the uncertainty estimation on experimental measurements in Fig. 7 using data from Ref. 3. The data contain GDPT measurements of a Poiseuille flow in a microchannel with a rectangular cross-section of 380 $\upmu$m $\times$ 100 $\upmu$m. The flow is directed along the $x$ direction and the uncertainty was estimated along the perpendicular directions $y$ and $z$, i.e., the directions of zero flow velocity. For a $N_{\mathrm{S}}^{\prime}=0.03$ (very few overlapping particles) and across a height of $h=100~\upmu$m, we obtained uncertainties $\sigma_{x,y}=0.08$ px and $\sigma_{z}/h=0.01$, in good agreement with the values obtained in the analysis of the synthetic images.