14th International Symposium on Particle Image Velocimetry – ISPIV 2021 August 1–5, 2021 Calibration correction of arbitrary optical distortions by non-parametric 3D disparity field for planar and volumetric PIV/LPT D. Michaelis1*, J. Wiswall2, A. Mychkovsky2, R. Prevost3, D. Neal3, B. Wieneke1 1 LaVision GmbH, 37081 Goettingen, Germany 2 Naval Nuclear Laboratory, West Mifflin, PA 15122, USA 3 LaVision Inc., Ypsilanti, MI 48197, USA *dmichaelis@lavision.de Abstract In this study, a new image calibration approach is presented that corrects arbitrary optical distortions by utilizing non-parametric, 3D disparity fields. A calibration plate with a high spatial resolution (i.e., high density of calibration marks) was used to identify optical distortions that remain after the initial calibration, which were then used to create a correction field for the pinhole or polynomial mapping functions. Results from a pipe flow experiment with four cameras using volume self-calibration (VSC) and Shake-the-Box Lagrangian particle tracking (STB LPT) are presented and the impact of the improved calibration is discussed. Using the calibration marks with the correction field, distortions of initially more than 20 pixels are reduced below 1 pixel. Using VSC with the correction field allows further reduction of average calibration disparities below 0.02 pixels (maximum 0.5 pixels), whereas without a correction field the remaining average disparity is much higher at 1 pixel (maximum 5 pixels). STB analysis of the data shows a considerable higher spatial resolution at the pipe wall and a consistent spatial distribution of the number of detected particles in the measurement volume. 1 Introduction Optical flow-measurement techniques, such as particle image velocimetry (PIV) and LPT, require precise imaging of tracer particles in the fluid to yield accurate velocity field measurements. Many practical internal flows are through conduits with curved surfaces and experimental setups often consist of water or air as the working fluid with glass or acrylic viewing windows. Small radius window curvatures, large camera viewing angles, and large refractive index differences create strong optical distortions; and optical flow-measurements in even simple tubes or pipes become challenging (see Van Doorne (2007)). Optical distortions can be reduced or prevented by refractive index matching the working fluid and viewing window materials: Bai and Katz (2014), Northrup et al. (1991) and Wright et al. (2017). Such experiments can yield excellent results if the setup in question allows for the application of this technique. Index matching, however, cannot be applied in every relevant case. The index matched fluid may not have the desired fluid properties, like viscosity or density, can be corrosive or even toxic, and can be very expensive. Also, the experiments may require delicate temperature control as the refractive index changes with temperature, which adds complexity to the setup. Finally, for all experiments with typical gases, including 14th International Symposium on Particle Image Velocimetry – ISPIV 2021 August 1–5, 2021 the most common case of air, index matching simply is not possible, as there are no solid materials which match an index of refraction close to unity. In the cases where refraction index matching is not applied, the optical distortions often lead to an inadequate image calibration. The mapping functions in common use (e.g., the pinhole: Tsai (1987) or polynomial model: Soloff et al. (1997)) are not able to correct for strong optical distortions, so that the corrected images, which should represent a regular Cartesian coordinate space, still contain considerable distortions after the image correction process. This is because the commonly used mapping functions, mentioned above, are not able to correct for arbitrary optical distortions due to their limited number of degrees of freedom. These high distortions result in large velocity measurement errors or prevent a successful particle matching process for particle tracking methods. Historically, these distortions have been mostly ignored or circumvented by restricting the measurement domain to regions with weaker, acceptable optical distortions. This, unfortunately, may in many cases exclude regions that are of special interest, like areas close to the surface where sometimes very complex, flow structure interactions occur. 2 Description of new correction method Concept The mapping functions M map 3D world coordinates xw to 2D camera sensor positions xc: 𝒙𝑐 = 𝑴(𝒙𝑀) (1) These mapping functions can be arbitrary functions and common choices are the aforementioned pinhole or polynomial models. To correct for strong optical distortions, the mapping function model needs to be extended to allow the description of arbitrary distortions. Therefore, a correction field C is added to the result of the original mapping function to provide a 2C correction vector for each position in 3D space: 𝒙𝑐 = 𝑴(𝒙𝑀) + π‘ͺ(𝒙𝑀) (2) Note that a different mapping function M and correction field C will exist for every camera in a multi- camera setup. In the current implementation C is a 3D2C vector field, which allows a simple and effective initial correction based on calibration marks. The real-world grid positions of the 2C correction vectors in the 3D correction field correspond to the real-world position and spacing of the calibration marks in x and y and to the position and shift of the calibration plate in z. Correction from calibration marks For each 2D mark position in the camera image xmc the corresponding 3D real word mark position xmw is known from the calibration process and is used to solve for the mapping function M. π’™π‘šπ‘ = 𝑴(π’™π‘šπ‘€) (3) For a perfect mapping function, the mapping of a world point xmw on the regular grid of calibration marks coincides with the position where the mark has been detected on the camera sensor xmc. In the presence of strong optical distortions, the remapped world positions of the marks in the corrected images can deviate considerably from the expected positions on the grid xw (Fig. 3, left) due to the limited number of degrees of freedom of the common mapping functions. The correction field C, for each mark position in all calibration z-planes, is determined by subtracting the mapped position M(xw) from the actual mark position in the camera image xmc for each position on the grid xw (note that the 2D correction vectors are defined in the 2D camera space): 14th International Symposium on Particle Image Velocimetry – ISPIV 2021 August 1–5, 2021 π‘ͺ(𝒙𝑀) = π’™π‘šπ‘ βˆ’π‘΄(𝒙𝑀) (4) Form the discrete 2C-3D correction field on a regular grid, the correction for any world position is calculated by linear vector interpolation. Remaining disparities from volume self-calibration VSC is a procedure for optimizing mapping functions based on the remapped triangulation error of individual tracer particles: Wieneke (2008). With this procedure, calibration inaccuracies can be removed down to and below 0.1 pixels. In experiments with strong or high-order optical distortions, often the remaining disparity cannot be lowered below 0.3 to 1.0 pixel. In such cases, the insufficient number of degrees of freedom of the mapping functions prevents further optimization. The disparity maps from volume self-calibration (Wieneke (2008)) have the same format as the correction fields C, which means that the remaining disparity fields after volume self-calibration can directly be used as additional correction fields CVSC. If VSC is performed using a calibration that already has a correction field, then the updated correction field Cupdate can simply be calculated by adding the remaining disparities CVSC to the existing correction field C: π‘ͺπ‘’π‘π‘‘π‘Žπ‘‘π‘’(𝒙𝑀) = 𝑰(π‘ͺ(𝒙𝑀)) + π‘ͺ𝑉𝑆𝐢(𝒙𝑀) (5) where I is an interpolator that converts correction vectors from the grid of calibration mark positions to the grid of the disparity maps. 3 Experimental setup Figure 1. Left: test setup used to obtain particle image data. Right: test section design shown in a CAD model rendering. Fig. 1, left shows a photograph of the experimental setup used to obtain data for the present study and Fig. 1, right shows an illustration of the test section used for flow visualization. The test section is a hexagonal block of acrylic with a 101.5 mm +/- 0.2 mm diameter hole bored along the centerline connecting two opposing faces. The tomographic PIV system was setup such that cross-sectional volumes, shown in Figure 1, of the pipe flow through the test section can be visualized. 14th International Symposium on Particle Image Velocimetry – ISPIV 2021 August 1–5, 2021 Water at room temperature was used as the fluid for the test. Steady, well-conditioned pipe flow was provided to the test section, and swirl in the flow was induced by an upstream swirl generator that injects a fraction of the total flow tangentially to the pipe. The design of the swirl generator followed Haden and Lorentz, (2018). The total flow through the test section, the sum of the axial and swirl injection flows, was 23 gallons per minute (GPM) with a standard deviation of 0.5 GPM. The fluid temperature was controlled to 74.5 ˚F Β± 3 ˚F. The fraction of flow delivered through the swirl generator was 20 % of the total flow and the remaining flow was delivered as well-conditioned pipe flow. Four high-speed, digital cameras were used to obtain images of the flow (Phantom VEO440L). The cameras were arranged at the four corners of a rectangular configuration, all viewing toward the upstream direction. The cameras were outfitted with 60 mm focal length lenses (Nikon) with Scheimpflug adapters (LaVision) and fluorescence bandpass filters that have 620 nm center wavelength and 56 nm band width (Edmund Optics #33-911). Fluorescent orange tracer particles, 50 Β΅m diameter with a peak emission wavelength of 606 nm, were used to visualize the flow (Cospheric UVPMS-BO-1.00 45-53Β΅m). The water was seeded to a target concentration of 0.02 particles per pixel as observed in the illuminated volume. The volume was illuminated using a dual cavity Nd:YLF laser (Quantronix Darwin duo). The z-axis of the illuminated volume was aligned to be in-line with the pipe axis. Volume optics (LaVision Inc.) were used in conjunction with an aperture to achieve an illuminated volume 7 mm thick in the z-direction and spanning the entire pipe cross section in the x- and y-directions. A single-plane calibration plate, mounted to a traverse, was used to provide the initial calibration for the images. The calibration plate consists of 1,909 dots, each with a 1 mm diameter and a center-to-center spacing of 2 mm (Fig. 2). The traverse has a resolution of 0.02 mm in its position and was aligned with the z-direction parallel to the pipe axis. The calibration plate was inserted into the test section through the slot (shown in Fig. 1, right) and rotated 90 degrees, similar to a butterfly valve, to provide the x y plane of for the PIV images. The plate was traversed along the slot length to obtain images at multiple z locations. Images of the calibration plate were taken at 1 mm increments and spanned the entire illuminated volume. Figure 2. 2D calibration plate with 2 mm dot spacing (sum of corrected images from all four cameras). 4 Results Calibration with correction field 14th International Symposium on Particle Image Velocimetry – ISPIV 2021 August 1–5, 2021 For the initial calibration, a z-scan of 12 views from the calibration target was recorded with an increment of 1 mm in the z-direction (-2 to +9 mm, where z = 0 mm corresponds to the downstream end of the illuminated volume). The calibration dot finding algorithm in LaVision DaVis 10.2 was used to identify the marks in each image. Only 0-3 marks at the border were not detected in each image due to strong distortions, so that nearly all of the marks were used to calculate the initial polynomial mapping functions (Fig. 3, left), as well as the correction vectors for each detected mark (Fig. 4). The original polynomial mapping function is not able to compensate the strong optical distortion near the pipe wall (Fig. 3, left), showing that this mapping is inadequate in these border regions. Applying the correction field (Eq. (4)) compensates for these optical distortions and results in a nearly perfectly corrected image (Fig. 3, right). Inspecting the correction fields in more details reveals that all cameras show a systematic distortion/correction field, which is related to the cameras position and viewing angle (Fig. 4). The correction field is similar between cameras, due to the symmetry of the camera setup. The scatter plot in Fig. 5 provides a quantitative visualization of the distortion range and magnitudes. Up to 22 pixels of distortion can be observed. Such large and uncorrected distortions prevent any further meaningful analysis in these regions, either by particle (particle tracking) or interrogation window based (planar PIV, Tomographic PIV) techniques. Volume Self-Calibration (VSC) with correction field VSC is a mandatory step for high-particle-density particle tracking or tomographic PIV (Wieneke (2008)). It removes any remaining disparities of the calibration and/or compensates for slight camera dislocations after the recording of calibration images. VSC was performed using the same particle images, from the 7 mm thick illuminated volume, that were later used for the final flow analysis. The disparity map and disparity vectors for each sub-volume were calculated using 20 x 20 x 2 sub-volumes. VSC was applied iteratively in up to 13 steps, to improve disparity vector precision, especially at the borders. Without a correction field, the polynomial mapping function is not able to adapt to the resulting complex disparity field. The initial average disparity magnitude (Fig. 6, left) is 1.8 pixel (6.4 pixel maximum) and Figure 3: Cutout of corrected calibration image from a pipe flow experiment. Left: polynomial calibration model (systematic misalignment between markers and Cartesian grid of up to 22 pixel), right: polynomial model with correction field (all calibration marks are centered on the Cartesian grid). The red arrow exemplifies a single correction vector according to Eq. (2). 14th International Symposium on Particle Image Velocimetry – ISPIV 2021 August 1–5, 2021 the average disparity remains at 1.0 pixel (4.8 pixel maximum) after the 8th iteration (Fig. 6, right). The large remaining disparity after VSC is problematic for high accuracy particle tracking and tomographic PIV and indicates that an accurate calibration over the measurement region cannot be achieved. Performing VSC in the present study using the calibration with a correction field, Eq. (5), behaves similarly to experimental cases without strong distortions. In the initial iteration, disparities of up to 1.7 pixels were detected (Fig. 7, left). The remaining disparity reduced to an average disparity of only 0.02 pixels after the 13th iteration, with very few outliers of up to 0.5 pixels (Fig. 7, right). Figure 4: Correction fields for the different cameras (one from 12 planes at z = 3 mm). Correction vectors outside the circle are extrapolated from the vectors based on calibration marks (inside the circle). Note how the vectors in the β€œFig. 3” area point from the expected position on the grid to the actual mark position. Particle Tracking with correction field STB software in the upcoming DaVis 10.2.1 version (LaVision GmbH) was used to calculate multi-pass STB, where the time series of particle images is utilized multiple times, both in the forward and backward direction, to improve the results from pass to pass. Three passes over 1112 time-steps are applied to obtain an average of 7488 matched particles per time-step, without, and 8463 matched particles per time-step, with Cam 1 Cam 2 Cam 3 Cam 4 Figure 3 14th International Symposium on Particle Image Velocimetry – ISPIV 2021 August 1–5, 2021 the correction field. Time-averaged velocities on a regular grid were calculated from particle tracks using Gaussian average binning with a bin size of 1 x 1 x 1 mmΒ³ at 75% overlap, resulting in 0.25 mm vector spacing and about 500 particles per bin (Figs. 8 and 9). In these figures, the swirl induced in the pipe flow is clearly visible, both with and without a correction field. Also, the overall flow structure and the peak velocity are similar in both cases. However, without a correction field, the strong gradient at the pipe wall cannot be resolved. The calculated w-component levels without using a correction field are unphysical resulting in 0.1 m/s (about 50% of peak velocity) within 2 mm distance from the pipe wall (Fig. 8, left). With a correction field, the gradient is well resolved and the measured out of plane components drops to a realistic 0 m/s in a ring of 0.25 mm width at the wall (Fig. 8, right). Figure 5: Scatterplot of vector components Vx and Vy in pixel from all correction vectors (Figure 5). Note the inverted y-scale: positive y components point downwards in image 4. Cam 1 Cam 2 Cam 3 Cam 4 14th International Symposium on Particle Image Velocimetry – ISPIV 2021 August 1–5, 2021 Figure 6: Scatter plot of the two components, Vx and Vy, of disparity vectors from the 1st (left) and the 8th (right) iteration of volume self-calibration without correction field applied (camera , Vx and Vy in pixels) Figure 7: Scatter plot of the two components, Vx and Vy, of disparity vectors from the 1st (left) and the 13th (right) iteration of volume self-calibration with correction field applied (camera 1). The units of Vx and Vy are in pixels. Note the different scales compared to Fig. 6. The pattern of the number of particles contributing to a single bin (or vector) is considerably different with and without a correction field (Fig. 9). For the case with a correction field, the number of particles contributing to a bin is nearly constant (about 500) and homogeneously distributed in the plane (Fig. 9, left). This flat distribution is expected as the seeding distribution in the particle images is also very homogeneous. The number only drops in a small ring at the wall, because the overlapping bins at the border reach outside the circular pipe region and therefore collect fewer particles than bins that are completely inside the pipe region. Without a correction field on the other hand (Fig. 9, right), the pattern of the number of particles is very uneven and complex. The number of particles is low in the regions of any camera view where the distortion is high. This is especially evident in a ring of about 2 mm width at the pipe wall. This can be explained by the allowed triangulation error (3 pixel in this case) not being high enough to find matching particles in these regions with stronger distortions. More particles are detected for the case with no correction field as compared with the case with a correction field. Such a result may seem surprising; 14th International Symposium on Particle Image Velocimetry – ISPIV 2021 August 1–5, 2021 however, this may be explained by a strategy used for STB, which favors the detection of particles with small triangulation errors. . Figure 8: Out of plane component w (color coded) and in plane vectors of the average flow field from central z-plane (out of 31 z-plane). Left: without correction field, right: with correction field applied, (every 8th vector displayed, original vector spacing is 0.25 mm). Figure 9: Number of matched particles per bin (color coded) and in plane vectors in central z-plane (out of 31 z-planes), every 8th vector displayed. Left: without correction field, right: with correction field. If no particle with a small triangulation can be detected, then there will be many possible particle-matches with similar high triangulation errors, from which only one is a real particle and the others are ghost particles. As STB cannot differentiate between them, all of them will end up as matching particles and in this way increase the number of particles per bin. 5 Conclusions 14th International Symposium on Particle Image Velocimetry – ISPIV 2021 August 1–5, 2021 A new method that addresses the longstanding problem of strong optical distortions for particle velocimetry has been proposed. Using correction fields, large distortions of up to 22 pixels could be compensated for and subsequently refined by volume self-calibration down to 0.02 pixel disparity. The resulting successful volume self-calibration enables precise particle tracking throughout the measurement volume, including the region close to the border of the measurement domain (in this case the inner wall of the pipe). Experimental results show that near-wall gradients could be resolved reliably once this correction field was implemented, whereas without the correction field, the velocity at the border was overestimated (being 50% of the peak velocity). The concept of the correction field has been demonstrated for a pipe flow. This approach is expected to apply in general to experiments that exhibit optical distortions from curved walls. Acknowledgements Notice: This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, nor any of their contractors, subcontractors or their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or any third party’s use or the results of such use of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof or its contractors or subcontractors. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. References Bai K, Katz J (2014) On the refractive index of sodium iodide solutions for index matching in PIV. Experiments in Fluids 55, 1704 Haden, R. E., and Lorentz, D. G. (2018) Apparatus and method for generating swirling flow, United States Patent US 9,956,532 Northrup MA, Kulp TJ, Angel SM (1991) Fluorescent particle image velocimetry: application to flow measurement in refractive index-matched porous media. Applied Optics (21):3034-40 Schanz, D., Gesemann, S., & SchrΓΆder, A. (2016). Shake-The-Box: Lagrangian particle tracking at high particle image densities. Experiments in Fluids, 57(5), 1-27. Soloff S, Adrian R and Liu Z (1997) Distortion compensation for generalized stereoscopic particle image velocimetry. Meas Sci Technol 8:1441-1454. Tsai RY (1987) A Versatile Camera Calibration Technique for High-Accuracy 3D Machine Vision Metrology Using Off-the-Shelf TV Cameras and Lenses. IEEE Journal of Robotics and Automation 4: RA-3 Van Doorne, C. W. H., & Westerweel, J. (2007). Measurement of laminar, transitional and turbulent pipe flow using stereoscopic-PIV. Experiments in Fluids, 42(2), 259-279. Wieneke B (2008) Volume self-calibration for 3D particle image velocimetry. Experiments in Fluids 45:549-556 Wright, S. F., Zadrazil, I., & Markides, C. N. (2017). A review of solid–fluid selection options for optical- based measurements in single-phase liquid, two-phase liquid–liquid and multiphase solid–liquid flows. Experiments in Fluids, 58(9), 1-39.