UEF//eRepository
DSpace https://erepo.uef.fi
Rinnakkaistallenteet Luonnontieteiden ja metsätieteiden tiedekunta
2021
Perturbation Monte Carlo in
Quantitative Photoacoustic Tomography
Pulkkinen, Aki
SPIE
Artikkelit ja abstraktit tieteellisissä konferenssijulkaisuissa
© 2021 Society of Photo-Optical Instrumentation Engineers (SPIE) All rights reserved
http://dx.doi.org/10.1117/12.2615870
https://erepo.uef.fi/handle/123456789/26805
Downloaded from University of Eastern Finland's eRepository
Perturbation Monte Carlo in Quantitative Photoacoustic Tomography
Aki Pulkkinen∗,1, Aleksi Leino1,2, Tuomas Lunttila1, Meghdoot Mozumder1, Tanja Tarvainen1,3
1Department of Applied Physics, University of Eastern Finland, 70211 Kuopio, Finland
2Department of Physics, University of Helsinki, 00014 Helsinki, Finland
3Department of Computer Science, University College London, London WC1E 6BT, United Kingdom
∗Corresponding authors’ email address: Aki.Pulkkinen@uef.fi
Abstract: In this work, use of perturbation Monte Carlo is extended to solving inverse problem of quantitative photoacoustic tomography. The approach is demonstrated feasible for estimating optical absorption and scattering distributions. © 2021 The Author(s)
1. Introduction
In quantitative photoacoustic tomography, images of an initial pressure distribution are used as the data to form es- timates of optical parameter distributions of an imaged target [1]. The initial pressure distribution itself is obtained using image reconstruction on ultrasonic measurements following photoacoustic effect induced by absorption of an excitation light pulse. The estimation is commonly done using a model based inversion and relies on an optical propagation model, such as the radiative transfer equation [2], its diffusion approximation [3], or optical Monte Carlo [4]. In this work, the optical Monte Carlo is extended from being used as a forward model for simulating light propagation to the inverse problem of quantitative photoacoustic tomography [5]. This is achieved by using model based inversion and computing Jacobian matrices of the forward model using so called perturbation Monte Carlo approach [6,7].
2. Principles of perturbation Monte Carlo
The perturbation Monte Carlo approach is based on importance sampling. Given some probability distribution defined by its density function f(ξ)of some parametersξ defined in parameter spaceΞ, and a functiong(ξ)that computes some measurementMas an integral over the parameter space, the said integration can be approximated as a finite sum as
M= Z
Ξ
g(ξ)f(ξ)dξ = lim
N→∞
N
∑
n=1
g(ξn)≈
N
∑
n=1
g(ξn), (1)
whereξnare samples drawn from distribution defined by f(ξ), andNis some fixed large number. From importance sampling it follows that, for some perturbed distribution fδ(ξ), the perturbed measurementMδ can be computed using the same samplesξnof f(ξ)as [6]
Mδ = Z
Ξ
g(ξ)fδ(ξ)dξ= Z
Ξ
g(ξ)fδ(ξ)
f(ξ) f(ξ)dξ = lim
N→∞
N n=1
∑
g(ξn)fδ(ξn) f(ξn) ≈
N n=1
∑
g(ξn)fδ(ξn)
f(ξn), (2) provided that fraction ffδ(ξ(ξn)
n) can be evaluated. It turns out that the optical Monte Carlo can be expressed in the above form with f(ξ)describing stochastics of the light field,ξn being paths of individual computed photons (or photon packets), andMbeing e.g. the absorbed energy density or the fluence i.e. the output of the stochastic model [6,7]. After a lengthy derivation, the importance sampling together with difference quotient formed using MandMδ can be used to derive expressions for the Jacobians of the output of the optical Monte Carlo model with respect to the optical parameters defining the problem. For details on the derivation of the Jacobians, see Ref. [5].
3. Numerical simulations
A two dimensional origo-centric circular geometry (diameter 5 mm) was investigated. The true parameters of optical absorption µa and scatteringµs used to compute the synthetic data are shown in Fig.1. The data was generated using two illuminating irradiances defined using functions 1+cos(2θ)and 1+sin(2θ)withθ being the polar angle on the boundary of the circle. 107photons were used in each simulation. This was sufficient to ensure low stochastic noise in Monte Carlo simulations. The data was computed in a mesh composed of 2406
Fig. 1. From left to right: simulated optical absorptionµaand its MAP-estimateµa,MAP, simulated optical scatteringµsand its MAP-estimateµs,MAP, and above the estimates their relative error.
triangular elements and interpolated into a sparser mesh of 2292 triangular elements in which the estimates were computed. Normal distributed random noise (standard deviation 1% of peak-to-peak amplitude of noise free data) was added to mimick a real life noisy measurement.
Maximum a posteriori(MAP) estimatexMAP= (µa,MAP,µs,MAP)was computed following a Bayesian approach [8] by minimizing
xMAP=arg min
x
1
2kLe(y−f(x)−ηe)k2+1
2kLx(x−ηx)k2, (3)
whereyis the data, f(x)is the stochastic forward model, and an observation model of form y= f(x) +ehas been assumed, with the noiseeand prior information of unknownxhas been modelled as Gaussian distributeds.t.
e∼N (ηe,Γe)andx∼N(ηx,Γx)withηeandηxbeing the means andΓe= (L|eLe)−1andΓx= (L|xLx)−1being the covariances of the noise and the prior. For the noise model, accurate statistics of the noise added to the data was assumed. For the prior, an Ornstein-Uhlenbeck spatial prior of characteristic length of 1 mm was used with the mean corresponding to mean of the true parameter range and the standard deviation chosen to cover 95% of the range [9]. Minimization problem was solved using a Gauss-Newton algorithm iteratively as
xn+1=xn+αn Jf(xn)|L|eLeJf(xn) +L|xLx−1
−Jf(xn)|L|eLe(y−f(xn)−ηe) +L|xLx(xn−ηx)
, (4)
whereαn is a step-length parameter that could be obtained using a line-search algorithm but was fixed in our simulations, andJf(xn)is the Jacobian matrix of f(x)with respect toxatx=xncomputed using the perturbation Monte Carlo approach.
Fig.1shows the estimates of the optical absorptionµa,MAPand scatteringµs,MAPand their relative errors. As evident, the estimates qualitatively resemble the true parameter distributions. As it can be seen, the relative error of the absorption is small and the estimate is accurate, whereas the estimate of scattering is of lower quality with higher relative error. This is due to more ill-posed nature of estimating the optical scattering than estimating the absorption. The numerical simulations demonstrate that perturbation Monte Carlo is a feasible approach to solve the optical inverse problem of quantitative photoacoustic tomography.
Acknowledgments
This work was supported by the Academy of Finland under projects 314411, 336799, 312342 Centre of Excellence in Inverse Modelling and Imaging, 320166 Flagship Programme Photonics Research and Innovation, and Jane and Aatos Erkko Foundation.
References
1. B. Cox, J. G. Laufer, S. R. Arridge, P. C. Beard, “Quantitative spectroscopic photoacoustic imaging: A review,” J Biomed Opt17(6), 061202 (2012).
2. A. Ishimaru,Wave Propagation and Scattering in Random Media, vol. 1.(Academic, 1978).
3. S. R. Arridge, “Optical tomography in medical imaging,” Inv Probl15, R41–R93 (1999).
4. A. A. Leino, A. Pulkkinen, and T. Tarvainen, “ValoMC: a Monte Carlo software and MATLAB toolbox for simulating light transport in biological tissue,” OSA Continuum 2 (3), 957-972 (2019).
5. A. A. Leino, T. Lunttila, M. Mozumder, A. Pulkkinen, and T. Tarvainen, “Perturbation Monte Carlo Method for Quan- titative Photoacoustic Tomography,” IEEE Trans Med Imaging39(10), 2985-2995 (2020).
6. I. Lux, L. Koblinger,Monte Carlo Particle Transport Methods: Neutron Photon Calculations(CRC Press, 1991).
7. C. K. Hayakawa, J. Spanier, F. Bevilacque, A. K. Dunn, J. S. You, B. J. Tromberg, V. Venugopalan, “Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,” Opt Lett26(17), 1335 (2001).
8. J. Kaipio, E. Somersalo,Statistical and Computational Inverse Problems(Springer, 2005).
9. C. E. Rasmussen, C. K. I. Williams,Gaussian Processes for Machine Learning(MIT Press, 2006).
