Feasibility and accuracy assessment of light field (plenoptic) PIV flow-measurement technique

. A light field camera can enable measurement of all the three velocity components of a flow field inside a three-dimensional volume when implemented in a PIV measurement. Due to the usage of only one camera, the measurement procedure gets greatly simplified, as well as measurement of the flows with limited visual access also becomes possible. Due to these advantages, light field cameras and their usage in PIV measurements are actively studied. The overall procedure of obtaining an instantaneous flow field consists of imaging a seeded flow at two closely separated time instants, reconstructing the two volumetric distributions of the particles using algorithms such as MART, followed by obtaining the flow velocity through cross-correlations. In this study, we examined effects of various configuration parameters of a light field camera on the in-plane and the depth resolutions, obtained near-optimal parameters in a given case, and then used it to simulate a PIV measurement scenario in order to assess the reconstruction accuracy.


Introduction
Light field cameras allow recording three dimensional features of a given object by recording information about the angle of the incoming rays on the image sensors. It was initially developed by Ives [1] and Lippman [2] about a century ago, and the underlying science was explored by Gershun [3] when he defined the concept of light field. In modern light field cameras, recording of the angular information is achieved by placing a microlens array between the main lens and the image sensor. The raw image that forms on the sensor happens to be utterly complex, which are mathematically processed to obtain a human-perceivable image. The historical foundations and working principle of the camera are well described by Levoy [4].
Since a light field camera can record three dimensional features through angular information, it can be used to identify a seeded three-dimensional fluid flow as well. Thus, it can be used to conceptualize a threedimensional PIV measurement technique to obtain flow velocity. Despite this possibility, its usage in flow measurement has started only recently [5,6]. Fundamentally, there are two ways of using the camera in PIV flow measurements. The first way is to refocus the camera at a given depth plane and obtain a twodimensional particle image there. Thereafter, crosscorrelation between two such focused particle images coming from two snapshots of a seeded flow, respectively, would yield the two-dimensional velocity field in the focused plane. After refocusing the camera in multiple planes, multiple such two-dimensional velocity fields can be obtained, effectively yielding a 3D-2C flow field. Kawaguchi et al [7] used this concept to obtain velocity fields in two parallel planes, after performing a dual-plane PIV measurement. Although conceptually easy, the results suffer from severe inaccuracies, which stem from appearance of blurred out-of-the-plane particles in the refocused plane. The second way of using a light field camera is to reconstruct the threedimensional instantaneous particle distribution of the seeded flow by using a tomographic algorithm, followed by performing a three-dimensional cross-correlation in order to obtain all the three components of the velocity in the reconstructed volume. Nowadays, this latter method is being mostly studied due to its superior properties.
Lynch [5] and Lynch and Thurow [6] described usage of light field cameras in PIV measurements in detail. They simulated PIV-like scenarios, performed a three-dimensional reconstruction of the particle intensity field, discretized two such intensity fields in interrogation volumes, and then performed threedimensional cross-correlation to obtain velocity vectors. In subsequent studies, Fahringer and Thurow [8,9] reconstructed similar intensity fields using MART algorithm [10], which resulted in more realistic particles in the reconstructed volume. Despite that, locations of the reconstructed particles in the study by Fahringer and Thurow [8] matched with those of the original particles only approximately. Fahringer and Thurow [9] estimated positional accuracy of the reconstructed particles through simulations, after considering 40 spherical particles lying on the optical axis of the camera system. The results showed that difference between the reconstructed particle positions and those of the original ones are unacceptably large in the depth direction. Hadfield and Nobes [11] also observed large inaccuracies in the position of the seeded particles in a particle tracking velocimetry (PTV) measurement. Kawaguchi et al [7] and Ogawa et al [12] discussed the MART algorithm and used it to reconstruct a three-dimensional synthetic intensity field, from where they demonstrated trade-off between the thickness of the measured volume and the spatial resolution in the depth direction.
Many of the above studies yielded velocity vector fields also, but their accuracy is highly questionable. The severity became apparent from a study by Thurow et al [13], who compared results obtained by a light field PIV measurement with the results obtained by a stereo PIV measurement. The flow measurements were performed above a wing. The comparison showed that the two results differ greatly, sometimes even in a qualitative sense. There are various sources of such large error, e.g. large elongation of the reconstructed particles in the depth direction, which leads to accordingly large uncertainty in the location of their centroid (which represents their positions), imprecise knowledge of the distances between the main lens and the microlens array, optical distortions present in the main lens and lenses in the microlens array, imperfect alignment of the microlens array and the image-sensor plane, to name some. In the literature, studies trying to remove some of these inaccuracies are reported [8,14,15]. Foy and Vlochos [16] simulated scenarios where multiple light field cameras are used, and reported that high-fidelity measurements can be performed even with only two cameras.
It is evident from the literature survey that further accuracy assessment of usage of light field cameras in PIV measurements is required. In the present study, we have assessed effects of some camera construction parameters on the in-plane and the depth resolutions, followed by assessing positional accuracy of the reconstructed particles in a three-dimensional volume. We performed this study using computer simulations.

Fundamentals, theory, and simulation details
Three typical modern configurations of light field cameras are schematically shown in Fig.1. These are called Plenoptic-1, Plenoptic-2 (Galilean), and Plenoptic-2 (Keplerian) configurations, respectively. As shown in the figure, a light field camera consists of a main lens, a microlens array, and an image sensor. In the Plenoptic-1 arrangement, the main lens is focused on the microlens array. In Plenoptic-2 arrangements, the microlenses are focused on the image plane. In the Galilean type arrangement, the image plane of the main lens lies behind the image sensor, whereas it lies before the sensor in the Keplerian type of the arrangement. In the present study, the attention is paid on the Keplerian arrangement, because it offers more flexibility when configured physically. A detailed diagram of the studied camera system is shown in Fig.2. As shown, there are eight configuration parameters, namely, A, B, a, b, F, f, D, and d. Simultaneously, there are three optical constrains. Two of these constrains are the lens equations 1/A + 1/B = 1/F and 1/a + 1/b = 1/f, and the third one is F/D f/d, which stems from the fact that, for optimal imaging by microlenses (i.e. maximum but non-overlapping micro images), the F-number (indicated by F#) of the main lens and the microlenses should be almost equal [7]. Furthermore, d is kept constant in this study at 0.25 mm. This leaves us with four parameters that can be chosen freely.
In the next section, we will examine the effects of the four freely chosen parameters B/A ( M), b/a ( m), F, and f on the in-plane and the depth resolutions of the camera system. The two resolutions are represented by R xy and R z , respectively. In order to obtain quantitative values of R xy and R z , we consider a spherical particle of diameter 15 ȝm present at the location where the optical axis of the camera intersects the object plane (see Fig.2).
Light rays emerging from this particle pass through the main lens and the microlens array, and then recorded on the image sensor. The pixel size of the image sensor, į×į, is set to 5.15×5.15 ȝm 2 , which is approximately same as that typically found in commercially available cameras. We ensured that the spot size yielded by the particle is smaller than the pixel size, so that the calculated R xy and R z are independent of the size of the particle and represent true resolutions (if the particle size is significantly larger and its spot size covers multiple pixels, the calculated R xy and R z would depend on the particle size). In general, the light rays are recorded on multiple pixels, qualitatively similar to that shown in Fig.2. The recorded light intensity is averaged over the pixel area, followed by reconstructing the particle using MART algorithm [10]. The reconstructed particle is elongated in the depth direction. The intensity of the reconstructed particle is maximum on its centre, and it decreases in a Gaussian fashion along any radial line that extends from the centre. The cut-off intensity value is set to 3. The cut-off results into a well-defined elongated particle with clearly recognized boundary. Its width in the Z = 0 plane would be same both along the X and the Y axes, which is taken as the quantitative value of R xy . Similarly, its length along the Z axis is taken as the value of R z .
where F# is the common F-number of the main lens and the micro lenses, and D represents the effective diameter of the main lens (see Fig.2). In order to understand the way the resolutions depend on the configuration parameters, we analysed the ray traces. We calculated equations of the paths that a ray emerging from the location (X, Y, Z) = (0, h, 0) at angle ĳ from the Z axis follows (i) between the main lens and the microlens array and (ii) between the microlens array and the image sensor. The two equations, respectively, can be written as follows: where: Y C = Y coordinate of the centre of the microlens on which the ray incidents, and -tan -1 {(D + 2h)/2A} ĳ + tan -1 {(D -2h)/2A}. The constrain on ĳ (-ĳ 1 ĳ ĳ 2 ) stems from the fact that the rays emerging from the location (0, h, 0) must fall inside the effective region of the main lens (see Fig.2).
For the rays emerging from the centre of the 15 ȝm spherical particle, h = 0. First, we used equations (7) and (8), and calculated the ray traces for a reference configuration where (M, m, F, f) = (1, 0.1, 30 mm, 1 mm), which yields A = B = 60 mm, a = 11 mm, b = 1.1 mm, D = 7.5 mm, and F# = 4. The ray traces near the microlens array and the image sensor are plotted in Fig.3. As evident, the rays emerging from different microlenses fall at different locations of the image sensor.
In order to qualitatively examine the way M affects the resolutions, we plotted the ray traces for M = 5 in Fig.4, while keeping m, F, and f same as those used in the reference configuration. It yields A = 36 mm, B = 180 mm, a = 11 mm, b = 1.1 mm, D = 7.5 mm, and F# = 4. As evident, the ray traces are different from that observed in Fig.3. After examining values of R xy and R z for multiple values of M, we found that R xy = (0.055 mm)/M and R z = (1.214 mm)/M. These values imply that the in-plane resolution is more than 20 times better than the depth resolution. Note that R xy § į/(M·m). The inversely proportional dependence of R z on M is due to the fact that, as M increases, the angle that the rays make onto the centre of the spherical particle (i.e. ĳ 1 +ĳ 2 ) also increases (see Fig.2).
We examined dependence of the resolutions on F by plotting ray traces for F = 150 mm, while keeping M, m, and f same as that in the reference configuration. It yields A = 300 mm, B = 300 mm, a = 11 mm, b = 1.1 mm, D = 37.5 mm, and F# = 4. The ray traces are shown in Fig.5. After comparing this figure with Fig.3, it is evident that the two sets of the ray traces are identical, except that the ray traces in Fig.5 are physically shifted in Z direction. Also, note that the effect of the increased F has been to increase A, B, and D linearly by a same proportion; thus, ĳ 1 +ĳ 2 remains unaltered. Consequently, the resolutions are also expected to remain independent of F. Indeed, after examining multiple values of F, we confirmed R xy = 0.055 mm and R z = 1.214 mm.  Here also, R xy § į/(M·m). R z , in general, is expected not to depend on m, as ĳ 1 +ĳ 2 also does not. In the present case, however, the effect appears perhaps due to some optical complexity arising from the small value of a.
We can observe the effect of f by comparing the ray traces for f = 3 mm with those in the reference configuration (see Fig.7 and Fig.3). For f = 3 mm, we obtain A = 60 mm, B = 60 mm, a = 33 mm, b = 3.3 mm, D = 2.5 mm, and F# = 12. It is evident from the comparison that the ray traces in Fig.7 are much flatter than those in Fig.3. After quantitatively examining the resolutions for multiple values of f, we obtained R xy = 0.055 mm and R z = 1.214f. The directly proportional dependence of R z on f results from smaller ĳ 1 +ĳ 2 , which is caused by the larger F#.  also exhibits the same qualitative dependences on M, F, and m, but increases when f increases. Although not presented, we observed similar qualitative dependences also when the camera is arranged in the Galilean configuration. Quantitatively, however, R xy and R z are little smaller than those in the Keplerian configuration.
In the next section, based on the qualitative insights obtained on the way R xy and R z depend on M, F, m, and f, we will obtain near-optimal set of the configuration parameters so that the values of R xy and R z are close to minimum. It is also important to note that a larger divergence angle at the image plane of the main lens results in a smaller R z . Therefore, we analysed the divergence angle ș (see Fig.2) for three types of commercially available lenses, whose (F, F#) = (35 mm, 1.4), (50 mm, 1.4), and (85 mm, 1.4). We also used two types of commercially available microlense arrays with (f, d) = (2.08 mm, 0.25 mm) and (0.6 mm, 0.20 mm). The results showed that ș depends weakly on the main lens, but strongly on the microlenses. For the three types of the main lenses used, ș varies in the range 4.5 o~6o when (f, d) = (2.08 mm, 0.25 mm) and in the range 12 o~1 6 o when (f, d) = (0.6 mm, 0.20 mm). For the Raytrix R5 camera, we found that ș is about 10 o . and the number of the pixels is set to 1280 1024. We image distribution of a set of spherical particles whose diameter is 50 ȝm and which are arranged in the Z = 0 plane. The particles are arranged on 1 mm spaced rectangular grids, as shown in Fig.8.

Reconstruction accuracy
We perform MART reconstruction from the imaged particle distribution and obtain particle distribution in the nominal image plane of the main lens, followed by back-projecting it across the main lens, in the nominal object plane. It allows comparison between the positions of the original particles and the position of the reconstructed particles, respectively. Details of the used MART algorithm are described elsewhere [7,12]. Next, we shift the plane of the particle distribution to Z = +1 mm and then repeat the imaging, the MART reconstruction, and the back-projection process. In total, the particles plane is shifted to 10 locations in the +Z directions and 10 location in the Z direction, allowing us to examine the position error in the three-dimensional space 10 mm X 10 mm, 8 mm Y 8 mm, 10 mm Z 10 mm. The combined size of the reconstructed volume becomes about 21.0×16.8×16.8 mm 3 , which we divided into 412×330×330 voxels. Consequently, the size of a voxel becomes about 51 51 51 ȝm 3 . The X and the Y positions of the centres of the particles are plotted in Fig.9 (see the top plot). The horizontal and the vertical axes show the positions of the original particles and the positions of the reconstructed particles, respectively. Ideally, the plotted data should lie on the dashed diagonal line. However, as evident, the data deviate from the line. The Z positions are also plotted in the same figure (see the bottom plot). The RMS deviation in the X and the Y directions is about 0.14 mm, and about 0.35 mm in the Z direction. Despite the significant deviation in the Z direction, our preliminary calculation of the displacement vectors based on two such reconstructed particle distributions show that the errors introduced in the vectors are rather small, implying that accuracy of PIV obtained velocity vectors perhaps depend only weakly on the position errors, making PIV flow measurement using such a light field camera possible. However, in order to assess the vector errors properly, a thorough investigation is required.

Conclusions
Dependence of R xy and R z on M, F, m, and f are examined. The results show that both R xy and R z are inversely proportional to M and m, but independent of F.
As for dependence on f, R xy remains constant, whereas R z is directly proportional to f. The dependences on m and f are stronger than that on M. The near-optimal configuration parameters for a commercially available microlens array with f = 0.6 mm and d = 0.2 mm can be given as F = 43.39 mm, D = 14.46 mm, A = 144.68 mm, B = 61.98 mm, a = 2.45 mm, and b = 0.79 mm. After simulating a light field PIV measurement scenario by using a light field camera configured with the optimal parameters, the difference between the positions of the MART reconstructed particles and the original particles came out to be small in the X and the Y directions, but significant in the Z direction. A preliminary calculation of the displacement vectors, however, did not exhibit large error. A more detailed investigation is required to thoroughly assess this error.