A novel method for correction of temporally- and spatially-variant optical distortion in planar particle image velocimetry

In-cylinder flow measurements are necessary to gain a fundamental understanding of swirl-supported, light-duty Diesel engine processes for high thermal efficiency and low emissions. Planar particle image velocimetry (PIV) can be used for non-intrusive, in situ measurement of swirl-plane velocity fields through a transparent piston. In order to keep the flow unchanged from all-metal engine operation, the geometry of the transparent piston must adapt the production-intent metal piston geometry. As a result, a temporally- and spatially-variant optical distortion is introduced to the particle images. To ensure reliable measurement of particle displacements, this work documents a systematic exploration of optical distortion quantification and a hybrid back-projection procedure that combines ray-tracing-based geometric and in situ manual back-projection approaches. The proposed hybrid back-projection method for the first time provides a time-efficient and robust way to process planar PIV measurements conducted in an optical research engine with temporally- and spatially-varying optical distortion. This method is based upon geometric ray tracing and serves as a universal tool for the correction of optical distortion with an arbitrary but axisymmetric piston crown window geometry. Analytical analysis demonstrates that the ignorance of optical distortion change during the PIV laser temporal interval may induce a significant error in instantaneous velocity measurements. With the proposed digital dewarping method, this piston-motion-induced error can be eliminated. Uncertainty analysis with simulated particle images provides guidance on whether to back-project particle images or back-project velocity fields in order to minimize dewarping-induced uncertainties. The optimal implementation is piston-geometry-dependent. For regions with significant change in nominal magnification factor, it is recommended to apply the proposed back-projection approach to particle images prior to PIV interrogation. For regions with significant dewarping-induced particle elongation (3$ ?>Ep>3), it is recommended to apply the proposed dewarping method to the vector fields resulting from PIV interrogation of raw particle image pairs.

In-cylinder flow measurements are necessary to gain a fundamental understanding of swirlsupported, light-duty Diesel engine processes for high thermal efficiency and low emissions. Planar particle image velocimetry (PIV) can be used for non-intrusive, in situ measurement of swirl-plane velocity fields through a transparent piston. In order to keep the flow unchanged from all-metal engine operation, the geometry of the transparent piston must adapt the production-intent metal piston geometry. As a result, a temporally-and spatially-variant optical distortion is introduced to the particle images. To ensure reliable measurement of particle displacements, this work documents a systematic exploration of optical distortion quantification and a hybrid back-projection procedure that combines ray-tracing-based geometric and in situ manual back-projection approaches.
The proposed hybrid back-projection method for the first time provides a time-efficient and robust way to process planar PIV measurements conducted in an optical research engine with temporally-and spatially-varying optical distortion. This method is based upon geometric ray tracing and serves as a universal tool for the correction of optical distortion with an arbitrary but axisymmetric piston crown window geometry. Analytical analysis demonstrates that the ignorance of optical distortion change during the PIV laser temporal interval may induce a significant error in instantaneous velocity measurements. With the proposed digital dewarping method, this piston-motion-induced error can be eliminated. Uncertainty analysis with simulated particle images provides guidance on whether to back-project particle images or back-project velocity fields in order to minimize dewarping-induced uncertainties. The optimal implementation is piston-geometry-dependent. For regions with significant change in nominal magnification factor, it is recommended to apply the proposed back-projection approach to particle images prior to PIV interrogation. For regions with significant dewarping-induced particle elongation ( > E 3 p ), it is recommended to apply the proposed dewarping method to the vector fields resulting from PIV interrogation of raw particle image pairs. Keywords: particle image velocimetry, optical distortion, swirl-plane, re-entrant piston geometry, ray tracing, back-projection, diesel engine (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Nomenclature α
Scaling factor of the projected object plane in the image plane (pixel mm −1 ) β Scaling factor of the 2D projected virtual surface in the image plane (mm pixel −1 ) γ 2 Parameter in Gaussian approximation d Laser temporal interval δ z Depth of field considering camera only φ Angle between particle velocity and radial direction λ Wavelength of laser

Introduction
Advanced internal combustion (IC) engines which convert chemical energy into mechanical work are considered as one of the promising solutions for propulsion purpose [1]. Directinjection, swirl-supported light-duty Diesel engines are a promising clean and efficient solution for road transportation, but represent one of the most challenging combustion environments to investigate. In-cylinder combustion is a complex, transient, 3D process that includes multi-physics interactions, multi-phase flows, coupling between turbulence and chemistry [2,3], and occurs in a matter of milliseconds. To gain fundamental understanding of IC engine processes for high thermal efficiency and low emissions, a considerable amount of research effort has been devoted to studying and characterizing in-cylinder flow [4]. In-cylinder flow plays a key role in many important processes such as fuel-air mixing [5], surface heat transfer [6], and pollutant transport [7]. Digital particle image velocimetry (PIV) has been widely used for non-intrusive in situ measurement of planar and volumetric velocity fields. Previous in-cylinder flow measurements led to observations that some flow structures are closely coupled with the piston bowl geometry [4]. In order to investigate the in-cylinder flow in a production-intent combustion device (in this case, a small-bore diesel engine), it is necessary to maintain realistic boundary conditions. Consequently, refraction at the curved piston surfaces will induce a spatially-varying optical distortion to images of horizontal laser planes taken through the optical piston. For reliable PIV results, it is crucial to improve the detection accuracy of particle displacements. To quantify true particle displacements under the presence of optical distortion, two approaches have been proposed in previous works: raytracing-based geometric back-projection [8] and exper imental in situ calibration back-projection [9,10]. The geometric back-projection approach requires an exact knowledge of the imaging parameters (such as lens focal length, the nominal magnification factor, and the actual positions of the object plane and image plane). Prasad and Adrian [8] demonstrated this approach for a twin-camera stereoscopic setup with ray tracing performed in three-dimensions based on geometric optics theory. Their validation comparison showed that the measured velocity profiles followed the same trend as the analytical solution. However, the horizontal and vertical velocities are underestimated by up to 75% compared to theor etical calculations. For the experimental in situ calibration backprojection approach, back-projection functions were approximated either as second-order polynomials [9] or cubic polynomials [10] for a 2D calibration scheme. The higher order terms in the projection functions were used to account for geometric distortions caused by imperfect imaging optics. The unknown polynomial coefficients were determined by performing a non-linear, least-squares fit with several experimentally-obtained image-object point pairs. Comparison between these two approaches by Kent et al [11] revealed that the geometric back-projection approach with ray tracing can lead to substantial errors if the geometric optical parameters are not determined accurately. On the other hand, the in situ back-projection approach yields smaller errors. Because of the robustness and straightforward implementation of the in situ back-projection approach, it has been widely used for planar PIV measurements conducted inside optical engines [7,[12][13][14][15][16]. However, the in situ back-projection approach requires a very tedious and time-consuming manual calibration process [17].
Moreover, it is challenging to exper imentally quanti fy the temporal change of optical distortion during a time interval on the order of 1-10 μs. Therefore, a simplification is usually adapted for dual-frame PIV: both frames are dewarped using the same back-projection function. The error introduced by this simplification may not be negligible.
The purpose of this work is to aid with the analysis of experimental planar PIV measurements in a real combustion device: a small-bore, light-duty diesel engine with a transparent production-intent bowl-in-piston geometry (figure 1). To achieve this purpose, a semi-automated, hybrid backprojection method-which combines geometric and in situ back-projection procedures-is developed to correct for the temporally-and spatially-varying optical distortion. Finally, an effort has been made to evaluate the error introduced by neglecting the change in optical distortion that occurs in the short laser temporal interval. An uncertainty analysis using simulation-based artificial particle images is also conducted to address the systematic error introduced by the proposed dewarping method. Figure 1 depicts a typical experimental planar PIV setup for swirl-plane studies inside a single-cylinder, light-duty optical Diesel engine [14]. The second harmonic output (λ = 532 nm) of a dual-head pulsed Nd:YAG laser was formed into a thin horizontal sheet (minimum sheet thickness is ~450 μm in the chamber center and ~800 μm near the wall) and was delivered into the combustion chamber via a side window for particle illumination. Optical access was provided by a fused-silica piston top, which is a geometrically accurate representation of a production-intent reentrant-bowl piston design. This axisymmetric re-entrant piston bowl geometry introduces significant spatially-and temporally-variant optical distortion. Double-exposure, dual-frame particle images were taken vertically through the piston while the engine was motored under non-combusting conditions without fuel injection.

Experimental setup
The aperture number ( f # ) of a monochromatic CCD camera lens (LaVision Imager Intense) was set at 11, resulting in a 2.5-pixel diffraction-limited particle image size [14]. Porous SiO 2 powder with a true density of 600 kg m −3 and particle mean diameter of 2.0 μm was used as a PIV tracer. Particle frequency response analysis at this flow condition indicates that the Stokes number is less than 0.1 throughout the full compression stroke [14]. Therefore, the particles used in this study provide acceptable flow tracing accuracy with errors below 1% [18,19]. Details of the PIV experimental setup and engine operation have been reported in a previous work [14]. Figure 2 demonstrates the complex optical distortion by imaging a calibration target (positioned in the laser plane) through piston at various crank angle degrees. Figure 2(a) is a round, 80 mm diameter calibration target filled with uniformly spaced cross symbols in a 2 mm × 2 mm pattern with concentric circles every 10 mm radii. Figure 2(b) shows the target frames imaged through the re-entrant piston, which is statically positioned for various crank angles after top dead center (°aTDC). Figure 2 shows that the optical distortion changes with piston position, and is less severe in the squish zone than in the bowl. Moreover, for all of the distorted target images shown, there exists an annular region near the bowl entrance for which the light does not pass through the piston and reach the camera. Consequently, this region will be masked and excluded from the dewarping processes.

Hybrid back-projection method for distortion correction
In principle, back-projection for optical distortion correction (dewarping) should account for the piston shape such that particle image pairs look as if no piston was present. Since the piston is axisymmetric and no azimuthal distortion occurs, the dewarping process can be simplified to cylindrical coordinates. With this simplification, a hybrid dewarping method combining ray-tracing-based geometric and manual backprojection procedures for planar PIV is proposed, as depicted in figure 3. The ray-tracing-based geometric back-projection procedure is developed to account for the transformation between the object plane and the virtual image of the object plane formed by the piston using well-characterized optical parameters (i.e. piston shape, piston location, and laser sheet position), as described by equation (1).
In equation (1), F is the back-projection function between the radial position (R (mm)) in the laser plane (object plane) and the radial position (r (mm)) in the virtual image. However, this ray-tracing-derived back-projection function does not account for the camera lens system and the position of the object with respect to the image plane. Therefore, the in situ manual backprojection approach is developed and applied to account for the impact of the dynamic optical distortion in the planar PIV imaging system as a whole. The imaged radial position on the CCD chip without the piston in place ( ′ R , undistorted domain) is related to the radial position in the projected object plane by a scaling factor (α), as described in equation (2). For a given imaging setup, α is a function of the object plane position relative to the fire deck (z).
Additionally, β accounts for the scaling of the imaged virtual surface in the distorted domain, as described in equation (3) where ′ r is the imaged radial position on the CCD chip with the piston in place. The factors α and β are calibrated with a manual back-projection approach.

Ray-tracing-based geometric back-projection
In order to quantify the changes in optical distortion caused by the piston motion, an automated ray-tracing-based geometric back-projection technique is developed ( figure 4). Using cylindrical coordinates, for a given point in the object plane located above the piston, the routine computes its virtual  image as follows. A large fan containing thousands of rays with their origins at the object-plane point is projected towards the piston at evenly spaced angular intervals. Each ray is refracted (by applying Snell's law), or absorbed (at the engine fire deck and cylinder wall) each time it interacts with a surface of the simulated combustion chamber. The placement of a stationary turning mirror (see [14], figure 1) inside the reciprocating piston assembly at 45° relative to the line-of-sight provides a view of nearly the entire combustion chamber. Only rays that pass through the bottom of the piston and through a long cylinder representing the turning mirror are assumed to reach the camera lens. These rays are divergent at the piston's bottom surface and  In this way, the radial mapping in the object plane and the radial coordinate of the virtual image can be established for any number of object plane points with any object plane location and for any piston position. Figure 5 shows the back-projection function F in the distorted domain as a function of piston position s and radial position in the object plane R. When the top of the piston is within approximately 8 mm of the object plane, the entirety of the 41 mm radii cylinder region is mapped to a 28 mm radii  cylinder region. As the piston moves away from the object plane, the virtual image of the object plane is compressed into an increasingly small space. For example, when the piston is located 66.1 mm below the fire deck (at −110°aTDC or −250°aTDC), the 41 mm radii cylinder region is mapped into a cylinder region with a radius of 15.5 mm.

In situ manual back-projection
The following paragraph presents the detailed calibration procedure for the factor α. With the optical piston removed from the engine, images are taken with the calibration target located at various vertical positions below the cylinder head. Image registration techniques (with MATLAB built-in functions) are used to determine the scaling factor α for several calibration target locations. Based on paraxial approximations [20], α can be approximated by a linear correlation with respect to the laser plane position z (figure 6).
Calibration of the scaling factor β is more complex. In figure 7, simulated rays from each radial position (red dots) in the object plane are shown with the corresponding virtual image locations (blue dots) after being refracted by the piston. Therefore, the tractability of the problem is greatly improved by the simplified assumption that the raw particle image on the CCD chip is this virtual surface (composed of blue dots) viewed through the compound lens system. Since the shape and location of the virtual surface changes with piston position due to the changing distortion pattern, the scaling factor β is spatially-and temporally-variant.
The factor β is computed as the ratio between radial coordinates in the virtual surface (r (mm)) and radial coordinates in the image plane ( ′ r (pixel)). As shown in figure 8, radial coordinates in the image plane ( ′ r ) are measured directly in pixel coordinates from the raw calibration target image. Knowledge of the coordinates of crosses in the object plane (R), combined with the ray tracing results, provides the radial locations of the virtual images of the same crosses (r). β is defined in equation (4) to be the ratio between the radial position in the virtual image and the radial position in the image plane; it is a function of ′ r at a given piston position. Figure 9 is the contour plot of calibrated β as a function of radial position in the image plane ( ′ r ) and piston position (s). The factor β depends less on radial position in the squish zone than in the bowl region, because the piston surface is much less curved in the squish region. When the piston is close to the object plane (e.g. z = 5.8 mm or CAD = −25°aTDC), β exhibits large fluctuations near the chamber center. This is because a part of object plane near the chamber center is imaged through the piston pip, which has a much smaller radius of curvature than the other portions of the piston bowl ( figure 7(b)). As the piston moves away from the object plane, the virtual surface is formed from rays which are confined to an increasingly narrow region ( figure 7(d)). This means that β becomes less sensitive to piston position as the piston moves downward (figure 9).

Evaluation of systematic error induced by piston motion with a fixed distortion correction transformation
This evaluation is devoted to justify the importance of quanti fying temporally-variant optical distortion for doubleexposure, dual-frame PIV. The following effort investigates the systematic error of radial velocity introduced in PIV measurements if the changing optical distortion due to piston motion is neglected. As shown in figure 3, the radial position in the object plane (R) is modeled as a function of radial position in the virtual surface (r) and time t, as described by equation (5) . In all figures, piston is at −45°aTDC during the compression stroke and there exists an annular blank region separating the squish and bowl regions. This blank region is masked and excluded from the velocity calculation due to its invisibility and severe optical distortion. The systematic error in → v in figure 12(b),e as a function of angle φ (angle between particle velocity and radial direction) and particle radial position (R [mm]) when piston is at −45°aTDC.
where F is the ray-tracing-derived back-projection function (figure 5) and The radial velocity v r , by definition, is computed as the total derivative of the radial coordinate: In which, . Therefore: If the changing optical distortion due to piston motion is , the radial velocity becomes: Therefore, the erroneous contribution in radial velocity (scaled by mean piston speed; here S p = 4.51 m s −1 ) will be: Figure 10(a) shows the systematic error (e) of radial velocity defined by equation (9). The radial velocity will be overestimated by up to 40.1% of the mean piston speed in the squish region and underestimated by up to −53.4% of the mean piston speed near the centrally-mounted injector exit (R = 0). e is most significant when the piston is located around −90°aTDC and decreases with reduced as the piston approaches either TDC or BDC. Figure 10(b) provides a 2D view of radial velocity error at −45°aTDC during the compression stroke. The radial velocity is biased inward near the chamber center and outward in the annular ring between radii of 21.4 mm and 41 mm. The systematic error is minimized at a radial position of 21.4 mm. This error distribution is also dependent on piston position. Since equation (9) is not a function of laser temporal interval (Δt), minimizing Δt with a fixed distortion correction transformation applied to both particle frames will not help eliminate this error in the radial velocity component. Therefore, the temporal change of optical distortion should not be neglected for dual-frame PIV. The proposed raytracing-based geometric back-projection procedure provides a time-efficient solution to quantify the change in optical distortion in microsecond-order time intervals.
The implementation of back-projection techniques is usually divided into two categories. The first approach (referred to as as Method A) is to dewarp the particle images prior to PIV interrogation [7,9,[11][12][13][15][16][17]. In the current application, once the virtual surface back-projection function F (figure 5), α (figure 6), and β (figure 9) are determined, the corresponding radial location in the object plane ( ′ R i j , ( ) (pixel)) can be computed for any radial position in the image plane ( ′ r i j , ( ) (pixel)). Therefore, the pixels of raw images can be re-located to positions ( ( ) ′ R i j , (pixel)) in the undistorted domain. Since the transverse magnification from the object plane to the virtual images is always less than unity, the field-of-view occupies a larger spatial range in the undistorted domain than in the distorted domain. Intensities of unmapped pixels are determined via bilinear interpolation. After the back-projection for each particle image is conducted, the PIV interrogation process is performed using the commercial software DaVis 8 from LaVision. Uniformlyspaced grids are used with a multi-pass interrogation strategy (window sizes: 128 × 128 pixels² and 16 × 16 pixels²). A demonstration of optical-distortion-corrected PIV results is shown in appendix A.
The second back-projection approach (referred to as Method B) is to dewarp the PIV-interrogated velocity fields into the undistorted domain [10,17]. In the current application, PIV interrogation is performed using the non-dewarped  (7). Since the piston is axisymmetric, the tangential velocity component in the undistorted domain (v t ) is calculated as follows: Finally, interrogation windows in the distorted domain can be re-located in the undistorted domain according to the back projection functions (F, α and β) evaluated at the piston position of 1st particle image for each image pair. The uncertainty in re-located interrogation window position is small when piston speed and/or laser temporal interval (Δt) is low. Each approach has its own advantages and drawbacks. Method A: dewarping the raw particle images before interrogating them may change the particle shapes and their intensity distributions, and will therefore introduce interpolation errors [21]. Method B: performing PIV interrogation in the distorted domain ignores variations in magnification factor within each interrogation window. As a consequence, results will be biased towards the motion of visually-large particles in the raw images. With the proposed hybrid backprojection method, one purpose of this paper is to evaluate the dewarping-induced uncertainties brought by each implementation.

Dewarping-induced uncertainty analysis
For both methods, it follows that the accuracy of particle displacement calculations will not only be affected by the accuracy of the ray-tracing-derived back-projection function (figure 5), and factors α (figure 6) and β (figure 9), but also by the changes in particle intensity distribution due to interpolation (Method A) or optical distortion (Method B). So, it is necessary to quantify the degree of uncertainty the above steps will bring, as a whole, to velocity results. In order to isolate the dewarping-induced uncertainty from other experimental uncertainties (e.g. particle lag error, pixelization, signal-tonoise ratio, velocity gradients, out-of-plane motion, etc [22]), simulated particle images are created and utilized. Figure 11 demonstrates an effort to evaluate the dewarpinginduced uncertainties with randomly distributed, simulated particle images generated with a MATLAB routine. To achieve some statistical significance, 50 pairs of particle images with a constant 3.68-pixel vertical shift are generated (details in appendix B) and processed using with the two back-projection methods described above. This uncertainty quantification is performed throughout the compression stroke at a range of piston positions (appendix C) and a crank angle during the compression stroke (−45°aTDC) is chosen as an example of this effort to highlight here. Assuming the laser temporal interval is 10 μs (simulated velocity is therefore → = → v j 29.8 s m s −1 ), interrogation operations produce 50 instantaneous velocity fields, one of which is shown in figure 11(d).
For Method A, the shapes of tracer particles become radially elongated after the dewarping process, especially inside the bowl region ( figure 11(c)). This is caused by the interpolation scheme that fills in unmapped data after the distortion correction has been performed. The degree of elongation can be calculated with a given ray-tracing-derived back-projection function (figure 5), and factors α (figure 6) and β (figure 9); additional details can be found in appendix D. (%) (figures 12(c) and (f)), it is evident that method B (with the exception of the region within 5 mm of the chamber center) results in less uncertainty inside the bowl region than Method A. In the squish region, both methods result in similar dewarping-induced uncertainties that are typically smaller than for the bowl region. Detailed comparison between these two methods is discussed in the next paragraph. and (e)) as a function of angle φ (angle between particle velocity and radial direction, figure 13 left) and particle radial position (R (mm)). The changes in the sign of the error for both methods cannot be explained either by the monotonic elongation (figure D1) or by monotonic back-projection ( figure 5). For method A, the error in the bowl region is considerably more complex than in the squish region: both positive and negative errors exist locally. For particles with large radial velocity components (φ is closer to 0°), the error varies from −30.0% to 36.2%. When the elongated particle has a large tangential velocity component (φ is closer to 90°), the error becomes smaller. Moreover, the error distribution exhibits local maxima and minima when φ is small. This pattern may be the result of radially-dependent elongation (appendix D) due to the change of surface curvature inside the bowl. For method B, up to 72% underestimation in velocity magnitudes occur near the injector exit region ( ⩽ R 5 mm) where the nominal magnification within each interrogation window changes significantly because of the piston pip geometry. Except injector exit region, dewarping-induced uncertainty with method B (mean error = −1.2%, root-mean-square error = 8.7%) in the bowl is much less than that with method A. In the squish region, dewarping-induced uncertainty with method B (mean error = 1.3%, root-mean-square error = 4.9%) is similar to that with method A (mean error = 1.2%, root-mean-square error = 5.0%).
Comparisons at multiple piston positions during compression stroke (figure C1, C2) indicate that dewarping-induced uncertainty with Method B is generally less severe than that with method A. The exception to this is the region near the combustion chamber center especially when piston is located farther away from the laser plane. Therefore, for the correction of temporally-and spatially-variant optical distortion in planar PIV, it is recommended to reconstruct velocity fields with the combination of both methods according to piston geometry. For regions with significant change in nominal magnification factor, less dewarping-induced uncertainty will occur if the proposed back-projection approach is applied to particle images prior to PIV interrogation. For regions with significant dewarpinginduced particle elongation ( > E 3 p ), it is recommended to apply the proposed dewarping method to velocity fields produced by PIV interrogation of distorted particle image pairs.

Summary and conclusions
The motivation of this work is to aid with PIV measurements taken in horizontal (swirl) planes in a real combustion device: a small-bore, light-duty Diesel engine with a productionintent, reentrant-bowl piston geometry. Images taken through a reciprocating optical piston with a complex geometry are distorted, which results in additional errors and uncertainties that are not present in PIV measurements with simpler optical setups. This work is devoted to the development of a semi-automated, hybrid back-projection universal approach to correct the temporally-and spatially-variant optical distortion with minimum dewarping-induced uncertainty. This work documents the following achievements: 1. A hybrid back-projection method, which combines ray-tracing-based geometric and in situ manual back-projection procedure, has been implemented. The geometric back-projection technique is based on an automated, ray-tracing based routine and is applied to quantify the optical distortion caused by refraction at the piston surfaces ( figure 4). The in situ back-projection procedure is able to account for the realistic behavior of the imaging system as a whole. This hybrid back-projection method exhibits as a universal tool for the quantification of temporally-and spatially-variant optical distortion with arbitrary, axisymmetric piston geometries. 2. Analysis of the image dewarping routine demonstrates that the change in optical distortion during the PIV laser temporal interval is piston-geometry-dependent, and therefore cannot be neglected ( figure 10). 3. The proposed hybrid back-projection method provides a time-efficient and robust way to process planar PIV measurements conducted in a temporally-and spatiallyvarying optical distortion environment. This work has demonstrated that whether the proposed back-projection method is used to dewarp particle images or dewarp velocity fields has a significant impact on the uncertainty of velocity measurements. The optimal implementation is piston-geometry-dependent. For regions with significant change in nominal magnification factor, it is recommended to apply the proposed back-projection approach to particle images prior to PIV interrogation. For regions with significant dewarping-induced particle elongation ( > E 3 p , appendix D), it is recommended to apply the proposed dewarping method after PIV interrogation is performed on distorted particle image pairs.

Appendix A. A demonstration of optical distortion correction and PIV results
With the calibration procedure detailed in the main text, once the virtual surface back-projection function (figure 5), α (figure 6) and β (figure 9) are determined, optical distortion at any swirl plane and at any crank angle position can be corrected by interpolating on the above two calibration surfaces and curve. Figure A1 demonstrates a qualitative comparison of target images before and after distortion correction at various crank angle degrees.
After the back-projection for each particle image is conducted, the PIV interrogation process is performed (LaVision DaVis 8). A multi-pass interrogation strategy with 128 × 128-pixel and 16 × 16-pixel window sizes was used. The laser temporal interval is tuned to meet the PIV quarter rule ( [22], pp 350-1). The detailed PIV interrogation settings are documented in a previous work [14]. Figure A2 shows three instantaneous velocity fields taken during the late compression stroke. It is obvious that the true swirl center (marked by the black dot) is different due to cyclic variability. All of Figure A1. A demonstration of spatial-and temporal-variant optical distortion correction with the proposed hybrid back-projection method. the velocity fields exhibit strong bore-scale tangential velocity components (swirl).

Appendix B. Artificial particle image simulation
Artificial particle images are generated to isolate the dewarping-induced uncertainties introduced to PIV interrogation from those brought on by other experimental factors, such as: in-plane and out-of-plane motion, velocity gradient, seed density, laser fluence, etc. This simulation is composed of five steps: 1. As shown in figure 11(a), simulated PIV particles are generated with a random distribution in the object plane (within ∅82 mm bore). In the object plane, the particle diameter d p follows a normal distribution with a mean diameter of 2 μm and a coefficient of variation (CV) of 1% (ideally monodisperse). In the absence of the piston, the particle intensity in the image plane can be approximated as a Gaussian distribution ( [22], pp 100-1): Considering the signal-to-noise ratio for PIV tests inside this optical engine, the peak Gaussian intensity ( ) h 0 2 is modelled to have a normal distribution with mean intensity of 128 counts and a CV of 30%, which is a typical particle image intensity distribution for PIV measurements inside the optical engine. γ 2 is 3.67 ( [22], pp 100-1). Laser fluence is assumed to be uniform in the field of view. 2. For the second exposure frame, a vertical 3.68-pixel shift is imposed ( figure 11(a)). This displacement, which is 23% of the given 16x16-pixel minimum interrogation window size, yields the maximum detection probability ( [22], p 304). 3. The effects of refraction at the piston surfaces are taken into account: particle image centroid locations in the image plane are calculated based on their locations in the object plane using the distortion correction technique developed in this work ( figure 11(b)). Assuming all of the particle virtual images are perfectly focused onto the image plane, the particle image diameter τ d is: where the diffraction-limited spot diameter is given by  9). Therefore, the particle intensity profile in the image plane follows their Gaussian profile (equation (B.1)) with changed particle image diameter τ d . To be noted, the above approximation of particle image diameter τ d is valid only when particle virtual images are perfectly focused onto the image plane. According to equation (B.3) ( [22], page 9), the depth of field for this experimental setup, δ z , is estimated to be 47.3 mm. The focal length of the lens was adjusted such that its depth of field covered the axial range of the curved virtual surface for piston positions closer than or equal to 55.6 mm (crank angles after −95°aTDC). When the piston is farther away from the object plane, equation (B.2) underestimates the particle image diameter τ d : is the camera nominal magnification (in the absence of the piston), d r is pixel pitch ( = d 6.45 r (μm/pixel)), and λ is the wavelength of the illuminating laser. 4. With the proposed hybrid back-projection procedure, optically-distorted particle images are dewarped into the undistorted domain ( figure 11(c)). 5. After the dewarping is finished, PIV interrogation (a multi-pass interrogation strategy with 128 × 128-pixel and 16 × 16-pixel window sizes) is performed on each image pair. Figure 11(d) shows an example of an instantaneous velocity field. Appendix D. Dewarping-induced particle image radial elongation

Appendix C. Uncertainty analysis (relative errors in
The particle radial elongation factor (E p ) is defined as the magnification ratio of a 1-pixel-sized region from the image plane to the object plane after dewarping: , no radial elongation is induced. If > E 1 p , particle images are elongated radially. E p is presented as a function of piston position and particle radial position in figure D1. It can be seen that radial elongation factor in the bowl region is more sensitive to particle radial position and piston location than in the squish region. Moreover, E p in the squish region stays around unity while E p in the bowl region becomes larger as piston moves away from the fire deck. 29.8 m/s Figure D1. Particle image radial elongation as a function of piston position (s (mm)) and particle radial position (R (mm)). Object plane location is z = 3 mm.