On the characterization of bias errors in defocusing-based 3D particle tracking velocimetry for microfluidics

In this work, we provide a systematic theoretical and experimental characterization of bias errors in defocusing particle tracking (DPT) methods based on a single calibration function, with respect to microfluidic applications, in which it is not possible to use calibration targets inside the measurement volume. This approach is widely used in microfluidic experiments, but bias errors are often neglected and to date only few works reported empirical procedures to compensate for that. A systematic characterization of the impact of such error in DPT measurements is still lacking. We show that the field curvature aberration and the refractive index mismatch are the main sources of bias error in these applications. We present a correction methodology for the bias error based on the determination of a reference surface, and in addition we propose a procedure based on a reference measurement of a Poiseuille flow to determine the reference surface on microfluidic channels with constant cross section. We discuss the impact of the refractive index mismatch and how to correctly compensate for it. We validated our methodology and quantified the bias errors on 10 different experimental setups, using different working fluids, materials, geometries, and microscope objective lenses ranging from 5×\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} to 40×\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} magnification. Our results indicate that the impact of this type of bias errors is in general not predictable and must be evaluated case by case. The proposed methodology allows to estimate and minimize the bias error in most microfluidic setups and is suitable for any single-camera DPT approach.


Introduction
Defocusing is a well-known concept in photography: objects lying in proximity of the focal plane of the optical system appear sharp or in focus in the image, whereas more distant objects appear blurred or out of focus. The degree of defocusing is proportional to the objects' distance from the focal plane, thus encoding information about their out-ofplane position. This concept was used in experimental fluid mechanics already in the early 1990s to obtain the threedimensional (3D) coordinates of tracer particles from single-camera recordings (Willert and Gharib 1992;Stolz and Köhler 1994).
Since then, a large variety of approaches have been proposed to read out efficiently the defocusing information, for instance using a measurement of the particle image diameter (Stolz and Köhler 1994;Speidel et al. 2003;Afik 2015), a pinhole mask (Willert and Gharib 1992;Pereira and Gharib 2002), or a controlled astigmatic aberration induced by a cylindrical lens (Kao and Verkman 1994;Cierpka et al. 2010). More recently, pattern recognition algorithms based on cross-correlation (Barnkob et al. 2015) or neural networks (König et al. 2020;Lammertse et al. 2022) have been used as well. We will generally refer to all these methods as defocusing particle tracking (DPT).
DPT methods found a large field of application in microfluidics where normally only a single optical access through the microscope objective is available and the multi-camera approaches commonly used in turbulence and large-scale facilities are not applicable (Schanz et al. 2016;Raffel et al. 2018). For instance, they have been used for velocimetry measurements inside evaporating droplets (Rossi et al. 2019), acoustic streaming flows (Qiu et al. 2019;Sachs et al. 2022), electrolyte convection (Tschulik et al. 2011), and for tracking bacteria (Taute et al. 2015).
DPT methods rely on a calibration procedure to map the defocusing degree of a particle image to the respective out-of-plane coordinate. In microfluidic applications, the calibration is typically performed on fixed particles (e.g., stuck to a wall) changing the focus of the microscope. Bias errors can occur if the mapping obtained from the calibration is different from the actual mapping in the experiments, for instance due to unwanted aberrations, misalignment, or modified optical conditions. It is worth mentioning that some of these bias errors, such for instance the ones due to refractive index mismatch, are not present if a calibration target could be inserted and displaced inside the measurement volume (Hain and Kähler 2006), however, this is often not possible in microfluidic experiments due to the limited accessibility and the small scale environment.
Many DPT measurement approaches currently used and developed in the literature rely on a single defocusing calibration function for the entire field of view (FOV), but this assumption is not always valid, especially if a cylindrical lens is used or the measurements are performed across a large field of view (FOV) ). An empirical procedure to overcome this problem in wall-bounded volumes has been proposed by Fuchs et al. (2016) considering reference particle images moving in proximity of the wall. A similar approach was followed by Rossi and Barnkob (2020) considering tracer particles sedimented on the bottom wall. However, these approaches do not consider the impact of a misalignment of the reference wall or the effect of refraction across different media. In general, a systematic characterization of bias errors in this type of DPT measurements and their impact on the final uncertainty is still missing.
In this work, we intend to fill this gap by a thorough theoretical and experimental characterization of such bias errors. The main sources of bias errors are discussed in Sect. 2, and in Sect. 3 we propose a theoretical framework to obtain a correction method for biased systems. The impact of the bias errors in DPT measurements is characterized experimentally for the different setups in Sect. 4. The results are presented and discussed Sect. 5, whereas the conclusions are given in Sect. 6.

Bias error in DPT
We consider DPT approaches performed on monodispersed particles, with a single calibration function obtained by scanning a reference particle with the microscope focus, using a distortion-free microscope objective, and in microfluidic devices with planar optical accesses. Under these conditions, verified in the majority of DPT applications in microfluidics, the field curvature aberration and refraction effects can be considered the main contributions to the bias error. Perspective errors usually present in volumetric particle image velocimetry measurements (Raffel et al. 2018) are negligible for the optical arrangements considered here (Rossi and Barnkob 2020).

Field curvature
When imaging through curved lenses, the field curvature causes a geometrical distortion between the object and the image plane i.e., transforms front-parallel planar surfaces into curved image surfaces, as shown in Fig. 1a (Born and Wolf 1975). The mathematical description of this problem is given by the primary Seidel aberrations, which shows that the in-focus position in the image space of a given planar object varies with the squared distance to the optical axis (i.e., z ∝ x 2 + y 2 ).
Due to the generally curved surface of microscope objective lenses, this phenomenon is commonly present in DPT measurements: Particles sitting on a flat plane perpendicular to the optical axis and expected to have the same defocusing pattern, will show instead different patterns giving biased measured positions (Fig. 1a). A practical approach to map the field curvature is to use a second-order polynomial surface to fit measured positions of particles lying on a flat plane.
Astigmatic systems, which are also often used in DPT measurements, do not feature a single focal plane and as a result the field curvature will vary with the distance to the objective lens (Born and Wolf 1975). Since the measurement depth is usually small with respect to the length scale of this variation, we assume an invariant field curvature described as (1) z = ax 2 + by 2 + cxy.
where a, b, and c are coefficients to be determined experimentally.

Refraction
In DPT measurements, the light emitted by tracer particles travels through different media: the working fluid, the wall of the microfluidics device and immersion medium of the lens. Since these components have typically different refractive indices (i.e., refractive index mismatch, RIM), the light rays are deflected at each interface according to Snell's law of refraction: with n i , n g , and n f as the refractive indices, and i , g , and f as the angles of incidence in the immersion medium of the lens, the wall of the microchannel, and working fluid, respectively.
As a consequence, the focal plane of the objective lens will be shifted from z i to z f , as shown in the schematic of Fig. 1b, and according to the following equation (Kim and Liu 2004) where the t wall is the thickness of the bottom wall of the microfluidics device and the scaling factors k 1 and k 2 are given by (2) n i sin i = n g sin g = n f sin f , Using small angle approximation ( sin ≈ 0 ) (Rossi et al. 2012), the shift of the focal plane from Eq. (3) is simplified to For DPT measurements, we need to consider the depthdistance between particles when refraction is present Δz f . Using Eq. (6) to obtain Δz f , the constant terms depending on the thickness of the wall t wall cancel out and Δz f takes the following form The above Eq. (7) expresses the scaling factor due to the RIM as the ratio between the refractive indexes of the fluid n f and immersion medium of the lens n i , henceforth k. In the present work, the validity of this assumption in DPT measurements is verified with systematic experiments (Sects. 3.2 and 5.3).

Theoretical framework
For each detected particle, a DPT evaluation returns a position vector X defined by three non-dimensional coordinates (7) Δz f ≈ n f n i Δz i . Fig. 1 a Impact of the field curvature aberration in a standard imaging system, including: object plane, lens, image plane ( z ∝ y 2 ) and camera CCD. b Effect of light refraction patterns in a standard microfluidics setup, including the ray tracing with and without refraction a b Lens Object plane Image plane CCD Optical axis (X, Y, Z). The in-plane coordinates X and Y represent the particle position projected into the two-dimensional image plane and are scaled with the side length of the camera pixel divided by the magnification. Practically, the non-dimensional values of X and Y are equivalent to values given in pixel units. The out-of-plane coordinate Z defines the depth particle position, and it will be here scaled with the total depth of the measurement volume h, thus having values between 0 and 1. It should be noted that h is the actual depth of the measurement volume which is obtained from the depth covered by the calibration images, corrected with the scaling factor due to RIM. We will come back to this point later.
To obtain physically meaningful results, we need to project the non-dimensional coordinates into the physical space using a suitable mapping function x = T(X).
For unbiased systems, a linear mapping can be used where d px is the side length of the camera pixel, and M x and M y are the magnifications in the x and y direction, respectively (with M x = M y for non-astigmatic systems). For a biased system, a generic formulation of the mapping function can be given by where M is the linear mapping in Eq. (8), C is the translation vector due to the field curvature, and R is a rotation matrix that will be discussed later. Now, to derive experimentally the expression for C, we can consider a reference plane in the measurement volume. This could be a flat wall as in Rossi and Barnkob (2020) or the symmetry plane of a Poiseuille flow as proposed in Sect. 3.2. Following the analysis in Sect. 2, such plane in the non-dimensional coordinates will appear as a second-order polynomial surface where we omitted the constant term, since it can always be eliminated by setting the origin so that F(0, 0) = 0.
Considering the field curvature described by Eq. (1), it is easy to show that the tangent plane to F(X, Y) at the point where the field curvature is null (i.e., (X, Y) = (0, 0) ) is given by and corresponds to the projection of the reference plane in the non-dimensional coordinates if no aberrations were present. Thus, the difference between F(X, Y) and H(X, Y) represents the contribution from the field curvature and can be used to define the translation vector C as Finally, from H(X, Y), properly scaled, it is straightforward to identify the normal unit vector k , which describes the tilt of the reference plane with respect to the optical axis as shown in Fig. 2. The corresponding orthogonal unit vectors i and j lying on the reference plane, can be oriented arbitrarily, for instance pointing i in the flow direction as described in Sect. 3.2. The three orthogonal orientation vectors define the rotation matrix R used in Eq. (9). Although not strictly necessary to correct the bias error, it is useful to introduce the rotation matrix R in the mapping function to align the measurements with the reference plane.

Determination of the reference surface from a Poiseuille flow
The reference surface F(X, Y) has been obtained in previous publications by looking at particles sitting in proximity of a flat wall (Fuchs et al. 2016;Rossi and Barnkob 2020). Here, we introduce an alternative approach based on a reference measurement of a Poiseuille flow. We believe that, when applicable, this approach offers a series of advantages: (i) a Fig. 2 Schematic of the methodology used to corrected the bias errors in DPT measurement, including the reference surface F(X, Y), reference plane H(X, Y), orientation vectors ( i , j , k ) used to the rotation and the optical axis Poiseuille flow can easily be established in most microfluidic experiments using a microchannel with constant cross section, (ii) a large amount of measurement points across the FOV can be collected in a short time, and (iii) F(X, Y) is calculated from particle trajectories measured in the center of the microchannel, where the uncertainty is typically smaller than for particles close to the wall (Barnkob et al. 2015). Furthermore, it can provide a direct measurement of the scaling factor k for the RIM correction, as discussed in Sect. 5.3.
To derive F(X, Y) from the measurement of a Poiseuille flow, we propose the following methodology, graphically summarized in Fig. 3a-f. First, we set the reference frame so that X and Y are aligned with the stream-wise and transverse direction of the Poiseuille flow, respectively. Then, we slice the measurement volume along the Y direction into N bins and determine for each bin the velocity profile along the Z direction from the average velocity and mid-position of each trajectory (Fig. 3a). A parabola is fitted to each individual velocity profile to extract the position of the maximum velocity in each bin ( Y N , Z N ) (Fig. 3b). Notice that in microchannels with a symmetrical cross section (e.g., rectangular or square channels), the maximum velocity lies on a flat plane, thus ( Y N , Z N ) represents a curve that can be seen as the intersection between F(X, Y) and the YZ plane. Following the discussion from Sect. 3.1, we use a parabolic equation to map it: Z(Y) = aY 2 + bY + c , see Fig. 3c. Finally, each individual particle trajectory is relocated in depth using the parabolic fit Z(Y), see Fig. 3d-e. In this way, the 3D measured particle trajectories are moulded into a single surface which we can fit with (Eq. (10)) using a robust least-squares fit to determine the reference surface F(X, Y), see Fig. 3f. It should be noted that the theory developed in Sect. (3.1) requires that (X,Y) = (0,0) coincides with the optical axis. Here, we assume that the center of the image corresponds to the optical axis. We will discuss the consequences of possible differences in the actual position of the optical axis in Sect. (5.2).

Experimental setup
The general schematic of the experimental setup is shown in Fig. 4a. A microchannel was placed on the working stage of an epi-fluorescent inverted microscope coupled with a recording camera and LED backlight. For the present measurements, we used two different configurations. The first one consisted of a custom-made microscope equipped with a high-sensitivity sCMOS camera (pco.edge 5.5, PCO GmbH). For this setup, the optical system included a M = 5× , NA = 0.13 objective lens (EC EPIPlan, Zeiss AG), with a cylindrical lens f cyl = 500 mm (LJ1144RM-A, Thorlabs Inc.) placed in front of the camera sensor. This system was used to perform measurements in a borosilicate glass capillary with rectangular cross section w ch × h ch = 2000 × 200 m 2 (3520-050, VitroCom Inc.) and nominal thickness of 140 m, see Fig. 4c(i). The FOV in the X and Y directions was 2560 and 2000 pixels, respectively. The second configuration consisted of an inverted microscope (Leica DM IL LED, Leica GmbH) equipped with a high-speed camera (HighSpeedStar 4 G, LaVision GmbH). This system was operated with two objective lenses: M = 20× , NA = 0.4 (NPlan Epi, Leica GmbH) and M = 40× , NA = 0.55 (NPlan L, Leica GmbH). For the astigmatic measurements, the same type of cylindrical lens (LJ1144RM-A, Thorlabs Inc.) was placed in front of the camera sensor. Two different microchannels with rectangular cross section ( w ch × h ch ) were observed with this configuration: 50 × 50 m 2 (CS-10000087, Darwin Microfluidics) and 800 × 20 m 2 (CS-10000109, Darwin Microfluidics), with a corresponding FOV in the measurements of 1024 × 256 pixels and 1024 × 1024 pixels, respectively. The microchannels are made of a rigid polymer (Topas COC) with bottom wall thickness of 140 m. In some measurements, a microscope glass with 1 mm thickness (D100001, Deltalab) was placed under the microchannel, see Fig. 4c (ii)-(iii).
The flow was seeded with fluorescent polystyrene tracer particles with 9.9 m or 2.5 m diameter (ex/em 530/607 nm, PS-FluoRed, MicroParticles GmbH). All configurations used a backlight illumination provided by a high-power green LED (SolisC525C, Thorlabs Inc.) and a filter cube for fluorescence imaging (Excitation: BP 525/50 nm; Dichroic: 570nm and Emission: 620/60 nm).

Calibration images
From the combination of the different setups, we obtained a total of twelve independent configurations (lens-microchannel-fluid/particles) that were tested in this work (Fig. 4c (i)-(iii)). For each configuration, the corresponding measurement volume was mapped with a calibration stack of images . The calibration stack was obtained by taking images of particles stuck to the top or bottom wall of the microchannel, and modifying their defocusing patterns by moving the microscope objective, as shown in Fig. 4b. The step size was set to 2 m for the 5× lens, and to 1 m for all the other cases.

Fluid properties
In these experiments, two types of fluids were used. A solution of distilled (DI) water at 24,%wt. Optiprep (Opti-Prep T M Density Gradient Medium, Sigma-Aldrich), with density = 1102.1 kg/m 3 and refractive index n f = 1.351, and a solution of DI-water at 25 %wt. of urea (Thermo Scientific T M Urea, Thermo Fisher Scientific Inc.), with = 1064.6 kg/m 3 and n f = 1.352. The refractive index of the fluids was measured using a refractometer (Abbe 60, Bellingam & Stanley) working at the mean sodium D wavelength (589.6 nm). The refractive index at 293.15 K, n D t , was measured, with a reproducibility of 0.010 % and an accuracy of 0.035 %. All measurements were done in five replicate runs and buoyancy corrections were applied. Additionally, to confirm the density matching between working fluid and tracer particles, the fluid density was characterized using a vibrating tube densimeter (DSA-5000, Anton Paar). The standard uncertainty of measurement was found to be u( ) = 0.30 kg ⋅ m −3 and the expanded uncertainty at a 95 % confidence level ( k = 2 ) is U( ) = 0.46 kg ⋅ m −3 .

Poiseuille flow measurements
The Poiseuille flow was measured in total for eight different experimental configurations. For each optical arrangement the measurement depth h was larger than the depth of the microchannel h ch , hence no scanning procedures were needed. To check data reproducibility, the flow was recorded three times in a row for each configuration. In the 2000 x 20 m 2 microchannel, a flow rate of 50 l/min was imposed with a syringe pump (Legato110, KD Scientific). Each recording consisted of 250 frames taken at 25 frames per second (fps). For the other microchannels, the Poiseuille flow was imposed by hydrostatic pressure, by setting a height offset between the inlet and the collection reservoir. In each measurement 2500 images were recorded at 50 fps. The configurations for the different test cases, including the estimated coefficients ( p 20 , p 02 ) associated with the field curvature are summarized in Table 1 and will be discussed in details in Sects. 5.1 and 5.3. The experimental images were processed using the open-sourc software DefocusTracker (Barnkob and Rossi 2021), based on the general defocusing particle tracking (GDPT) method, available at https:// gitlab. com/ defoc ustra cking.

Bias error from the field curvature
We determined the field curvature using a reference surface obtained from the measurement of a Poiseuille flow as described in Sect. 3.2. The second order fitting coefficients p 20 and p 02 from the reference surface are used as evaluation parameters. It should be noted that the error from the field curvature at a distance of X∕Y = 1000 pixels from the optical axis corresponds to p 20 ∕p 02 × 10 6 (see Eq. (10)).
The values of p 20 and p 02 for each optical configuration are reported in Table 1 and are the average from three consecutive measurements.
A major assumption on our methodology is to consider the field curvature to be the same along the optical axis (Sects. 2.1 and 3). This assumption can be verified by comparing the measured particle trajectories before and after the field curvature correction. For each measured trajectory, we can calculate an error Z according to where Z traj is the Z-coordinate at X ≈ 0 , and Z i represents the collection of N traj that compose a given trajectory. After correction, the error Z should contain only contributions due to the random error, whereas before correction the contribution due to the bias error is also included. If the field curvature is indeed constant, its contribution at the various depth positions is equal, and hence the variation of the error due to the field curvature Z fc through the measurement volume should be constant (with Z fc calculated as the difference between Z before and after correction). Figure 5a-c shows the raw dataset (i), corrected dataset (ii), and the mean absolute error Z including the contribution from the field curvature Z fc (iii), for several optical configurations. The results show that the contribution from the field curvature ( Z fc ) to the measurement error is independent of the depth position, as hypothesized. Similar results were also found in the remaining setups.
It is also interesting to notice the effect of the cylindrical lens in the field curvature, comparing the astigmatic and nonastigmatic cases shown in Table 1 ( w ch × h ch = 800 × 20 m 2 ). On one hand, for the non-astigmatic optical system we determined a field curvature with similar values for p 20 and p 02 , as a consequence of the axis-symmetry of the lens system. On the other hand, as we introduce the cylindrical Here, the cylindrical lens was rotated 90 degrees from the first to the second optical configuration which therefore results in a different direction of the maximum curvature. The field curvature aberration is also found to be different for each optical arrangement and no clear patterns between individual characteristics of the lens systems has emerged (Table 1). The impact of the field curvature on a DPT measurement depends also on the size of the FOV. Fig. 6 shows the evolution of three different curvatures from Table 1 as function of the distance Y to the optical axis at Y = 0 . A field curvature with p 02 = −0.274 × 10 −6 will give a bias error of 27 % the height of the measurement volume h at 1000 pixels distance from the optical axis (Fig. 6c), whereas if p 02 = −0.032 × 10 −6 the error reduces down to 3 % (Fig. 5 a(i)). On the other hand, if we consider a narrower FOV e.g., Y < 100 pixels, the impact of field curvature becomes negligible. In such cases, the proposed methodology is not able to interpret the small curvatures, hence the null entries in Table 1.
In sum, the present analysis shows that the impact of the field curvature is an individual characteristic of each optical system and experimental setup, and accordingly it must be evaluated case by case. a b c Fig. 5 Comparison of the field curvature for different optical setups, including the raw (i) and corrected particle cloud (ii) from the Poiseuille flow measurements, the mean absolute error Z along the depth-coordinate (iii); the contribution from the field curvature is given by Z fc . a Lens system: ×20∕0.4 ; b Lens system: ×20∕0.4 (A); c ×40∕0.55 . All measurements were performed in a microchannel with cross-sectional area of w ch × h ch = 50 × 50 m 2 a b c Fig. 6 Impact of the field curvature on the estimated depth-coordinate Z from a DPT evaluation as function of the distance Y from the optical axis. The associated error Z fc as indicated by the colormap. a Lens system ×5/0.13 (A), microchannel: w ch × h ch = 2000 × 200 m 2 . b Lens system ×20/0.4 (A), microchannel: w ch × h ch = 50 × 50 m 2 . c Lens system ×20/0.4, microchannel: w ch × h ch = 800 × 20 m 2 Finally, the estimation of the reference surface F(X, Y) from the Poiseuille flow (Fig. 7b) was compared to the surface estimated with sedimented particles (Fig. 7a). As for the latter, a set of reference surfaces F n (X, Y) was obtained from the measurement of the sedimented particles at n different focus position. The final surface F(X, Y) in Fig. 7a is determined from the average of all fits at the different depth positions. Here, we considered a single experimental configuration with astigmatic lens system, corresponding to Fig. 4c(iii). Whereas the estimated curvature along the Y-direction is similar for the two methods, a different value is estimated for the curvature along the X-direction. This discrepancy is most likely due to a different distribution of data points across the FOV. In particular, since the Poiseuille flow allows to collect a larger amount of data points more evenly spread across the FOV, it is reasonable to expect a more reliable result with this approach.

Tilt of the reference plane
We showed in Sect. 3.1 that the bias error encountered in this type of DPT measurements is essentially a combination of the field curvature aberration and the tilt of the reference plane with respect to the optical axis. If the tilt is neglected, a systematic error will also occur in the in-plane directions, proportional to the cosine of the tilt angle and the depth position of the particle. Furthermore, to calculate the tilt angle in Sect. 3.1, we made the assumption that the optical axis is located at the center of the image (see Fig. 2), but a slight misalignment can clearly occur in experimental systems. To quantify the impact of the tilt angle and misalignment of the optical axis we used here Monte Carlo simulations.
The conditions in the simulations were chosen to match experimental the test cases in Table 1 with a microchannel cross section of 800 × 20 m 2 . The imposed field curvature was defined according to the highest and lowest value detected in our experiments: p 20 , p 02 = −0.03 × 10 −6 and p 20 , p 02 = −0.28 × 10 −6 . For simplicity and without loss of generality, the tilt was only imposed on the Y-axis ( ) and varied from 0.025 to 3 degrees. The error ( OA ) in the position of the optical axis was varied from 20 to 100 pixels along the X direction.
For each case (different and OA), 19674 particle coordinates moving into a Poiseuille flow with a maximum displacement of 17.6 pixels were randomly created inside the measurement volume. The selected tilt and field curvature were applied to the simulated measurements. The coordinate positions obtained in this way were then corrected by applying the methodologies proposed in Sect. 3, with and without the tilt correction (i.e., directly subtracting F(X, Y) to the Z coordinates). The error was calculated as the difference between the true and corrected coordinates. The maximum error for each case is reported in Fig. 8.
As expected, the simulations show that the full correction with no error in the optical axis (tilt corr. and OA = 0 ) is able to eliminate the bias error for both X-and Z-directions and the remaining error is essentially due to random noise in the simulation and in the numerical procedure used to estimate the field curvature and rotation. If the tilt correction is neglected, the in-plane error max{ X } rapidly increases (Fig. 8a), however for tilt angles < 3 deg it remain below 2 pixels, which is non-negligible but still very small on the FOV of 870×820 pixels (used in the simulations). On the other hand, for a large field curvature, the tilt correction can have a detrimental effect on the accuracy of the correction, if the optical axis is not lying on the image center.
Overall, the Monte Carlo simulations show that the interplay between tilting and field curvature is not trivial. While including the tilt angle in the correction is clearly beneficial for small field curvature, it is not so straightforward for larger field curvature. Given the relatively small a b Fig. 7 Estimated field curvature F(X, Y). a sedimented particles on the bottom of the microchannel b Poiseuille flow: particle trajectories. Microchannel: w ch × h ch = 800 × 20 m 2 ; Fluid: DI-water with 25 wt.% Urea ( n f = 1.352 ); Lens system: ×20/0.4 (A). The colormap represents the value of the estimated F(X,Y) in the considered FOV error occurring in the in-plane direction, a possible strategy is to neglect the tilt angle correction and simply subtract the reference F(X, Y) to the data.

Refractive index mismatch
In general, the bias error due to refraction arises since the defocused particle images in the calibration and in the experiments originate from distinct optical paths. Calibration data are typically taken on static particles, while moving the microscope objective at different depth positions. Measurements are taken with a static objective and particles moving across a volume. As a consequence of the RIM, the calibration depth h ′ calculated as the difference between higher and lower position of the microscope objective will be different from the actual measurement depth h during experiments. In practice, the scaling factor k = h∕h � is normally approximated to k = n f ∕n i as discussed in Sect. 2.2. In this section, we verify this approximation for the scaling factor in different optical configurations and two working fluids with known refractive index. Besides, we evaluate if the calibration images can be obtained from sedimented particles or if this is an oversimplification of the problem.
We estimated the scaling factor k for the optical setups summarized in Table 1 using the ratio between the nominal h ch height of the microchannel and a measured value h ′ ch obtained with three different experimental approaches. Since we used air-immersed lenses ( n i = 1 ), the measured scaling factor could be considered as a measured effective refractive index of the fluid, henceforth n eff The first approach to obtain the experimental heights is based on particles stuck to the bottom and top wall of the microchannels. In total, we recorded two pairs of particles for each optical setup and microchannel, scanned through different focus positions. The chosen particles were adjacent and near the optical axis to minimize the contribution from the field curvature. Examples of the particle images can be seen in Fig. 9a-b (i), where P T denotes the particle on the top wall, and P B the particle on the bottom. Following the work of Rossi et al. (2012), the experimental height was determined using a focus function, specifically the Tenengrad function (Yeo et al. 1993). With this approach, we obtain for each case two curves representing the value of the focus function for the bottom-(dark-blue markers) and top-wall particle (light-blue markers) at different focus position of the microscope, as shown in Fig. 9a-b (i). The height is estimated based on distance between the two maxima i.e., the two focal positions. Note that the two local maxima observed in Fig. 9a(i) are due to the use of a cylindrical lens in the optical arrangement.
A second approach, based on the same stacks of images, consists of a direct determination of the depth coordinate of the adjacent particles with a DPT measurement and calculate the channel height as the difference between the two The height values should then remain constant across different focus positions as shown in Fig 9a-b (ii). In particular, the observed fluctuations of h ′ ch are within 2 % of measurement volume which is typically the uncertainty for this type of measurement  and no additional systematic effects can be observed. Similar results were also found for the other optical configurations and demonstrate that the scaling due to RIM is constant across the measurement volume. Finally, a third approach is based on the measurement of a Poiseuille flow. The z-coordinate of the top and bottom walls are extrapolated from the velocity profiles assuming no-slip boundary condition at the wall, and h ′ ch determined from the difference of the two coordinates. Figure 10 shows a comparison of the estimated n eff obtained with the three different methods, normalized with the actual refractive index of the working fluid n f and as a function of the lens systems. The error associated to the measurement of the microchannel height h ′ ch was estimated by taking the nominal microchannel height h ch as ground a b true value, and it was found to be between 1.2 % and 3.2 % for all considered cases.
In general, the estimated refractive index of the fluid n eff corresponds to the actual refractive index n f for the three methods within a small measurement error. Hence, on one hand, the approximation used for the scaling factor k holds for DPT measurements. On the other hand, the three methods lead to equivalent results. Besides, the results reveal that the impact of the bottom wall thickness t wall does not contribute significantly to the RIM.

Conclusion
We provided a characterization of the bias errors of DPT measurements in microfluidic experiments, in which a single calibration function obtained with a focus scanning procedure is used. We presented a theoretical model to describe the bias error based on the field curvature aberration of the lens, the tilt of the wall of the microfluidic device with respect to the optical axis, and refractive index mismatch of different media between the tracer particles and the lens. We showed that the bias error can be mapped and corrected using a reference surface, which can be obtained looking at particles close to a wall, as done previously. In addition, we proposed an approached based on a reference measurement of a Poiseuille flow, suitable for microfluidic experiments inside microchannels with constant cross section, that allows the determination of the reference surface from a large number of measurement points evenly distributed across the entire FOV.
We tested and validated the model and correction methodologies on different microfluidic setups with different optics (magnification from 5 × to 40× ), working fluids, and microchannel dimensions. The results are in good agreement with the model but the magnitude of the bias error has a large variability: While being negligible for some setups, we could observe an error up to 27 % of the measurement volume depth across the FOV. The study on refractive index mismatch confirms that the value of the scaling factor to be used to account for it corresponds with a good approximation to the refractive index of the working fluid.
Author Contributions GC, AM, AM, and MR defined the scientific objectives. GC and MR developed the theoretical framework. GC carried out the experimental work and data processing. MR supported the experimental work, data processing, and provided the DPT software. AR performed the measurement of refractive index of the fluid. GC, AM, AM, and MR discussed and interpreted the results. GC drafted the first version of the manuscript and prepared the figures. The finalized version was revised by AM, AM, and MR.
Funding Open access funding provided by Alma Mater Studiorum -Università di Bologna within the CRUI-CARE Agreement. GC acknowledges the PhD scholarship 2021.04780.BD granted by Fundaão para a Ciência e a Tecnologia (FCT). MR acknowledges the financial support by the VILLUM foundation (grant no. 00036098). GC, AM and AM acknowledge Fundaão para a Ciência e a Tecnologia (FCT) for partially financing the research trough Project Ref. PTDC/ EMETED/7801/202.

Data availability
The MATLAB code used in this article will be fully accessible in the public gitlab repository https:// gitlab. com/ defoc ustra cking/ defoc ustra ckerm atlab. The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no conflicts of interest to declare that are relevant to the content of this article.

Ethical approval Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.