14th International Symposium on Particle Image Velocimetry – ISPIV 2021 August 1-5, 2021 Approximate Bayesian framework for 3D reconstruction in a volumetric PIV/PTV measurement Sayantan Bhattacharya*, Ilias Bilionis and Pavlos P. Vlachos School of Mechanical Engineering, Purdue University, West Lafayette, IN 47907, USA *bhattac3@purdue.edu Abstract Non-invasive flow velocity measurement techniques like volumetric Particle Image Velocimetry (PIV) (Elsinga et al., 2006; Adrian and Westerweel, 2011) and Particle Tracking Velocimetry (PTV) (Maas, Gruen and Papantoniou, 1993) use multi-camera projections of tracer particle motion to resolve three- dimensional flow structures. A key step in the measurement chain involves reconstructing the 3D intensity field (PIV) or particle positions (PTV) given the projected images and known camera correspondence. Due to limited number of camera-views the projected particle images are non-unique making the inverse problem of volumetric reconstruction underdetermined. Moreover, higher particle concentration (>0.05 ppp) increases erroneous reconstructions or β€œghost” particles and decreases reconstruction accuracy. Current reconstruction methods either use voxel-based representation for intensity reconstruction (e.g. MART (Elsinga et al., 2006)) or a particle-based approach (e.g. IPR (Wieneke, 2013)) for 3D position estimation. The former method is computationally intensive and has a lesser positional accuracy due to stretched shape of the reconstructed particle along the line of sight. The latter compromises triangulation accuracy (Maas, Gruen and Papantoniou, 1993) due to overlapping particle images for higher particle concentrations. Thus, each method has its own challenges and the error in 3D reconstruction significantly affects the accuracy of the velocity measurement. Though, other methods like maximum-a-posteriori (MAP) estimation have been previously developed (Levitan and Herman, 1987; Bouman and Sauer, 1996) for computed Tomography data, it has not been explored for PIV/ PTV 3D reconstruction. Here, we use a MAP estimation framework to model and solve the inverse problem. The cost function is optimized using a stochastic gradient ascent (SGA) algorithm. Such an optimization can converge to a better local maximum and also use smaller image patches for efficient iterations. Our method uses a uniform prior model for the unknown 3d particle positions (π‘₯) and a forward model (𝑓) of the measurement chain to estimate the posterior distribution. Here, we assume a uniform probability distribution, π‘₯𝑖 ∼ π‘ˆ([π‘₯π‘šπ‘–π‘›, π‘₯π‘šπ‘Žπ‘₯]) , for each of the 𝑁 particle positions within the measurement domain bounds, π‘₯π‘šπ‘–π‘› and π‘₯π‘šπ‘Žπ‘₯. The observations 𝑦1:𝑀 = (𝑦1, … , 𝑦𝑀) are recorded images from M cameras with image size π‘ƒπ‘Ÿ Γ— 𝑃𝑐 pixels. We define the row and column position of a pixel using π‘Ÿ and 𝑐 respectively. The camera calibration function (β„‚π‘Ÿ π‘š(π‘₯), ℂ𝑐 π‘š(π‘₯)) provides the forward mapping from 3D particle space (π‘₯1:𝑁) to the image pixel coordinates (π‘Ÿ, 𝑐) for π‘šth camera. Thus, the projected pixel intensity using the forward model is written as, π‘“π‘šπ‘Ÿπ‘(π‘₯1:𝑁) = 𝐼0 2πœ‹πœŽπ‘‘ 2βˆ‘ exp⁑ (βˆ’ (π‘Ÿβˆ’β„‚π‘Ÿ π‘š(π‘₯𝑝)) 2 +(π‘βˆ’β„‚π‘ π‘š(π‘₯𝑝)) 2 2πœŽπ‘‘ 2 )𝑁 𝑝=1 (1). In equation (1), 𝐼0, πœŽπ‘‘ are image parameters estimated from the optical transfer function (OTF) of particle images. The likelihood of all observed data conditioned on particle positions is: 𝑝(𝑦1:𝑀| π‘₯1:𝑁) = ∏ ∏ ∏ 𝑝(π‘¦π‘šπ‘Ÿπ‘|π‘“π‘šπ‘Ÿπ‘(π‘₯1:𝑁)) 𝑃𝑐 𝑐=1 π‘ƒπ‘Ÿ π‘Ÿ=1 𝑀 π‘š=1 (2), where, each pixel is considered as an independent observation. Assuming a Log-normal noise the individual pixel likelihood is defined as a normal distribution 𝒩(log π‘¦π‘šπ‘Ÿπ‘| log π‘“π‘šπ‘Ÿπ‘(π‘₯1:𝑁) , 𝜎 2), with noise variance 𝜎2.The posterior 𝑝(π‘₯1:𝑁 | 𝑦1:𝑀) is expressed as a function of the likelihood and the prior and finally, the MAP estimate is given by π‘₯𝑀𝐴𝑃 = argmax[log(𝑝(π‘₯1:𝑁 | 𝑦1:𝑀)), which we solve using stochastic gradient descent. This framework was implemented using PyTorch with tensor formulations. The ADAM optimizer was used with varying learning rates (0.001 to 0.1) and a polynomial calibration function was used to build the forward model. Synthetic PIV images of size 64x64 pixels were generated with varying particle concentrations (0.01 to 0.1 ppp) to validate the algorithm. Figure 1a) shows better convergence with increasing learning rate. However, the initial results with random initializations show lower particle yield across the range of seeding density (Figure 1b)). The RMS error for the valid measurements were close to 0.2 pixels. The SGA algorithm was also tested with different batch size, and suitable learning rates. Lower batch size showed higher oscillations (Figure 1c)) in mean square loss as smaller random patches from different camera images did not contain the projections of the same particles. Subsequently, random image patches corresponding to the same sub-volume resulted in better convergence. The reconstructed positions for a random initialization tended to get trapped in local minima and to avoid such cases perturbed stochastic gradient ascent (PSGA) (Jin et al., 2021) was also implemented. Every 100 iterations if the change in loss function was less than 1%, the gradients were perturbed by a zero mean Normal distribution with varying standard deviations. The methodology is also tested for triangulation-based initializations and for random corresponding patches for each camera image. Furthermore, the sensitivity of the optimization framework to different image noise models (e.g., Bernoulli or Poisson distributions) is explored. Finally, we compare the current estimation of the reconstruction error with existing methods for a synthetic vortex ring case. Figure 1: Subplot (a) shows algorithm convergence for different learning rate for a whole image optimization. Subplot (b) shows RMS error in reconstructed positions for valid measurements and % of valid measurements or particle yield for different particle concentrations in particles per pixel (ppp). Subplot (c) shows variation in mean square loss for the cost function for different batch size. References Adrian, R. J. and Westerweel, J. (2011) Particle Image Velocimetry, Cambridge Aerospace Series. Bouman, C. A. and Sauer, K. (1996) β€˜A unified approach to statistical tomography using coordinate descent optimization’, IEEE Transactions on Image Processing, 5(3), pp. 480–492. doi: 10.1109/83.491321. Elsinga, G. E. et al. (2006) β€˜Tomographic particle image velocimetry’, Experiments in Fluids, 41(6), pp. 933–947. doi: 10.1007/s00348-006-0212-z. Jin, C. et al. (2021) β€˜On Nonconvex Optimization for Machine Learning: Gradients, Stochasticity, and Saddle Points’, J. ACM. New York, NY, USA: Association for Computing Machinery, 68(2). doi: 10.1145/3418526. Levitan, E. and Herman, G. T. (1987) β€˜A Maximum a Posteriori Probability Expectation Maximization Algorithm for Image Reconstruction in Emission Tomography’, IEEE Transactions on Medical Imaging, 6(3), pp. 185–192. doi: 10.1109/TMI.1987.4307826. Maas, H. G., Gruen, A. and Papantoniou, D. (1993) β€˜Particle tracking velocimetry in three-dimensional flows’, Experiments in Fluids, 15(2), pp. 133–146. doi: 10.1007/BF00190953. Wieneke, B. (2013) β€˜Iterative reconstruction of volumetric particle distribution’, Meas. Sci. Technol. Meas. Sci. Technol, 24(24), pp. 24008–14. doi: 10.1088/0957-0233/24/2/024008.