Direct Estimation of Optical Parameters From Photoacoustic Time Series in Quantitative Photoacoustic Tomography

Estimation of optical absorption and scattering of a target is an inverse problem associated with quantitative photoacoustic tomography. Conventionally, the problem is expressed as two folded. First, images of initial pressure distribution created by absorption of a light pulse are formed based on acoustic boundary measurements. Then, the optical properties are determined based on these photoacoustic images. The optical stage of the inverse problem can thus suffer from, for example, artefacts caused by the acoustic stage. These could be caused by imperfections in the acoustic measurement setting, of which an example is a limited view acoustic measurement geometry. In this work, the forward model of quantitative photoacoustic tomography is treated as a coupled acoustic and optical model and the inverse problem is solved by using a Bayesian approach. Spatial distribution of the optical properties of the imaged target are estimated directly from the photoacoustic time series in varying acoustic detection and optical illumination configurations. It is numerically demonstrated, that estimation of optical properties of the imaged target is feasible in limited view acoustic detection setting.

in the imaged target induces, via the photoacoustic effect, an excess pressure, called the initial acoustic pressure distribution. For absorber sizes of a few millimeters or less, this initial pressure distribution relaxes by radiating as a broadband pulse of sound in the ultrasonic frequency range. This sound can be observed on the boundaries of the imaged target and it forms photoacoustic time series (i.e. the measurement data). Common parameters of interest to be estimated in QPAT include: the spatial distribution of the light absorbing chromophores or the optical absorption, the optical scattering, and the Grüneisen parameter linking the initial pressure distribution to the absorbed optical energy distributions.
QPAT can be seen as the next step after photoacoustic tomography (PAT). In PAT, images of the initial pressure distribution, also called the photoacoustic images, are of principal interest. In many situations these images reflect, for example, the anatomical information of the target as produced by its internal distribution of absorbing chromophores. For a review on PAT in general, see e.g. [1]- [4]. Common methods to solve the acoustic inverse initial value problem of PAT include back projection, time reversal, and model based inversion approaches [5]- [8]. Reconstruction methods in specific acoustic measurement geometries, such as plane, cylinder, sphere, and polyhedra have been devised [5], [9]- [15], as well as methods that operate in more generic acoustic measurement setting [16], [17]. In addition to acoustically homogeneous situations, the acoustic inverse initial value problem with heterogeneous speed of sound [18]- [23] and attenuating propagation of sound [24]- [26] have been studied. In acoustically limited view setting, the acoustic inverse problem in QPAT becomes ill-posed and the estimates of the initial pressure distribution suffer from artefacts [27]. To alleviate the issue, regularization approaches using prior information on the spatial distribution of the initial pressure [28], [29] and model based inversion with Moore-Penrose pseudoinverse [30] have been used. Additionally, use of acoustic mirrors and scatterers can improve the estimation of the initial pressure distribution in acoustically limited view setting [31]- [34]. For a review on the mathematics of the acoustic inverse initial value problem in PAT, see e.g. [35].
The additional step that QPAT takes, the quantitation, is the estimation of the spatial distributions of the optical parameters that result in the formation of the initial pressure distribution [36]- [38]. Commonly, a model based inversion approach is used to solve the inverse problem in QPAT. Both, the radiative transfer equation (RTE) 0278-0062 © 2016 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. based optical propagation model [39]- [45] and its diffusion approximation (DA) [46]- [55] have been used. In addition to using photoacoustic images as the input data in the quantitation step, combined data formed by the photoacoustic images and diffuse optical tomography (DOT) measurement data have been investigated [42], [56], [57]. Use of measurements obtained with illuminations performed at multiple wavelengths (i.e. spectral measurements) in order to obtain estimates of the concentrations of specific chromophores based on their absorption spectra have been studied [58]- [63]. Within the DA framework, analytical reconstruction approaches have been devised [64], [65]. For mathematical analysis of the optical inverse problem in QPAT within the RTE and the DA frameworks, see e.g. [66], [67].
In addition to approaches where the initial pressure distribution is estimated first and the optical parameters of interest afterwards, direct estimation of the optical parameters from the photoacoustic time series has been performed recently. In [68], [69], the DA based optical framework was used to derive a linear Born approximation based model. The model was subsequently coupled with a homogeneous acoustic propagation model, and perturbations of optical absorption and diffusion parameters were estimated directly from the photoacoustic time series. In [70], the DA based Born approximation and the full non-linear DA was used to estimate the optical absorption, as well as the speed of sound of the imaged target. In [71], the absorption parameter was reconstructed within the nonlinear DA framework. In [72], the RTE based optical forward model was coupled with a homogeneous acoustic propagation model. The paper demonstrated improvement of the optical absorption estimates, when the estimates were formed directly from the photoacoustic time series. In all [70]- [72], scattering was assumed to be known. In [73], the studied photoacoustic illumination and acoustic detection setup was such that the acoustic inverse problem could not be solved in a stable fashion, and thus the inverse problem of QPAT required the use of a direct estimation approach.
In this paper, the inverse problem of estimating the optical parameters directly from the photoacoustic time series in QPAT is further investigated. Simultaneous estimation of both absorption and scattering, with parameter values corresponding to realistic values in soft biological tissues, is considered in two dimensional (2D) and three dimensional (3D) limited acoustic view situations with different optical illuminations. The inverse problem is described by using a Bayesian approach to ill-posed inverse problems [74]- [76]. In the Bayesian approach, the inverse problem is viewed in the framework of statistical inference. All parameters are treated as random variables which depend on each other through a model and information about these parameters is expressed by probability distributions. In quantitative tomography, point estimates of the posterior distributions of unknown parameters are determined based on the measurements, the model, and the prior information. Thus, instead of qualitative PAT images, in Bayesian approach to QPAT one is aiming at obtaining quantitative estimates of the parameters of (primary) interest and their uncertainties.
The forward model used in this paper is formed by combining an optical model with an acoustic propagation model. For the optical model, the DA is used. Green's function based solution to homogeneous acoustic wave equation is used for the acoustic model. The resulting forward model maps the optical parameters of the imaged medium to the photoacoustic time series on the boundary of the target.
The approach is tested with simulations. The effect of the acoustic measurement geometry, the effect of the noise level, and the effect of the optical illuminations used to generate the data are analysed.

II. MODEL
The photoacoustic effect can be modelled by coupling the optical and acoustic propagation models. The difference in propagation speed of light and sound means that the initial pressure distribution can be regarded as instantaneous with respect to the acoustic model. This means that the coupling of the models can be achieved by simple composition.

A. Optical Model
In this work, the propagation of light is modelled using the DA. The DA and the used boundary condition describing the flunce within the computation domain ⊂ R n are [77] −∇ · κ∇ + μ a = 0, where κ = (n(μ a +μ s )) −1 is the optical diffusion, μ a = μ a (r ) is the optical absorption parameter, μ s = μ s (r ) is the optical reduced scattering parameter, ζ n is a dimension dependent scaling factor (ζ 2 = π −1 and ζ 3 = 0.25). Further, A and s = s(r ) describe the light reflectivity and the inward light current (i.e. the light source) on the boundary ∂ . Throughout this work, value A = 1 is used for the reflectivity corresponding to matched refractive indices inside and outside . Absorption of light creates spatial absorbed energy density field H given by where the parameters are as in (1). In this work, the optical propagation model (1) is discretized using Galerkin finite element method (FEM) [78]. FEM discretization is performed in a regular grid composed of either square (in 2D) or cube (in 3D) elements. The optical properties μ a and μ s , as well as the fluence are discretized using Q piecewise bilinear (in 2D) or trilinear (in 3D) base functions defined at grid nodes r q ∈ R n , with q = 1, . . . , Q. Thus, the parameters and the fluence are approximated as where φ q is the base function corresponding to grid node r q , and μ a,q , μ s,q and q are the discrete representations of μ a , μ s and respectively. The discretized absorbed energy density field (2) can then be approximated as where q is the FEM solution of (1) with parameters μ h a , μ h s and s.

B. Acoustic Model
Propagation of sound, created by the instantaneous photoacoustic effect, in an infinite domain composed of homogeneous non-attenuating medium is described by the acoustic initial value problem [2], [5], [79] ⎧ ⎪ ⎪ ⎪ ⎨ where p is the acoustic pressure, c is the sound speed, t is the time, and p 0 is the initial pressure distribution created by the absorption of light pulse. The initial pressure distribution is given by where γ is the Grüneisen parameter, and H is the absorbed energy density as in (2). Throughout this work γ is treated as a known constant, although in general, this is not the case. In this work, the solution of (5) is approximated as a solution of the wave equation [80] 1 where g τ (t) is a (temporal) Gaussian distribution, with full width at half maximum τ , approximating the formation of the initial pressure distribution p 0 into the acoustic medium. As τ approaches zero, (7) becomes a better approximation to (5). Solution of (7) can be computed as a spatial convolution [81] where G τ is the Green's function corresponding to (7), which can be expressed as where F −1 is the (temporal) inverse Fourier transformation, j is the imaginary unit,ĝ τ is the Fourier transformation of g τ . G, in (9), is defined aŝ which is the solution of ∇ 2Ĝ + ω 2 c 2Ĝ = −δ(r ) [81], H 1 0 is Hankel function of the first kind.
In this work, p 0 (r ) in (7) is approximated using the same base functions as in the optical model as Thus, (8) is approximated as with the spatial integral approximated using Gaussian quadratures. The Fourier transforms, in (9), are approximated by using discrete Fourier transformation corresponding to the temporal discretization. In discrete form, the solution (8) at a position r can be written as where p(r ) = ( p(r, t 1 ), . . . , p(r, t T )) ∈ R T is a vector with values of the acoustic pressure at time instances t 1 ,..., t T , and w(r ) ∈ R T ×Q is a matrix describing the convolution integral in (8) as a linear operator, and p 0 = ( p 0,1 , . . . , p 0,Q ) ∈ R Q is the discrete initial pressure distribution (11) as a vector.

C. Photoacoustic Simulation Model
In this work, the simulations are carried out in a photoacoustic setting composed of a light absorbing target surrounded by varying amount, as well as arrangement, of point size acoustic detectors. The target and the detectors are assumed to be immersed in an infinite homogeneous and nonattenuating acoustic medium, and the target is assumed to have the same acoustic parameters as the medium surrounding it. The mathematical model for this system is the combination of (1) and (5).
Defining the acoustic detectors as being located at positions d k ∈ R n , with k = 1, . . . , K , the photoacoustic time series for illumination s m , with m = 1, . . . , M, can be expressed as where μ h a and μ h s are the optical parameters and are as in (13) and is a vector representation of the discretization of (6). Note that . The photoacoustic simulation model (14) can be written, for illumination s m , shortly as where ∈ R 2Q . For multiple illuminations, the forward model can be written as

III. RECONSTRUCTION METHOD
The inverse problem in QPAT, is to form estimates of the optical absorption and scattering based on the data. In this work, the inverse problem is approached in the Bayesian framework for ill-posed inverse problems.

A. Bayesian Approach
Solving an inverse problem of determining unknown parameter distribution x ∈ R 2Q given noisy measurements, or data, y ∈ R T K M requires definition of an observation model linking the parameters and the measurements. In this work, an observation model with additive error is used and it is written as (18), and e ∈ R T K M is a random variable denoting the additive error, or noise. The observation model (19) links the parameters of interest to the measurements.
Within the Bayesian approach [74]- [76], the forward model (19) is interpreted statistically by defining probability distributions for the unknown parameter x and the noise e. The probability distribution of x is called the prior distribution and it is intended to describe the existing (rough) knowledge of the unknown x, while the noise statistics essentially characterises the measurement setup and modelling errors. The prior and the noise probability density function are marked with π x (x) and π e (e) respectively. Assuming that x and e are uncorrelated in the additive noise model (19), the Bayesian approach results in posterior distribution π(x|y) for the unknown x conditioned by the measurements y, and it is given by In principle, the posterior distributions of the unknown parameters can be estimated using Markov chain Monte Carlo (MCMC) methods. However, these methods can be computationally prohibitively too expensive in large dimensional tomographic inverse problems. Therefore, point estimates such as the maximum a posteriori (MAP) estimate are computed. In this work, the distributions π x and π e are modelled as Gaussian distributions, and their parameters where η x and x are the mean and the covariance of the prior information discussed in Section III-B, and η e and e are the mean and the covariance characterizing the noise and uncertainties of the measurement setup. Values for η e and e , characterizing the additive error in (19) are defined in Section IV. With the Gaussian choice for the distributions, the negative logarithm of the posterior distribution (20) becomes where L e and L x are the Cholesky decompositions of the inverse covariance matrices, such that L e L e = −1 e and L x L x = −1 x . The point estimate of x used in this paper is the MAP estimate which is obtained by using a Gauss-Newton algorithm [82]. The Gauss-Newton proceeds by starting from some initial value x = x 1 and using iteration where α l is obtained by using a line search algorithm to find minimum of u(x l + α l l ). l in (23) is the Gauss-Newton iteration direction at step l computed by solving where An approximation x MAP ≈ x l is then obtained with sufficiently large l.
One of the main assets of employing the Bayesian approach is that, in addition to point estimates, one can also compute (approximations for) the posterior error estimates of the unknown parameters of interest. In Bayesian analysis, the credibility intervals (counterparts of confidence intervals in frequentist statistics) would be the standard choice for error estimates. The computation of credibility intervals would, however, call for Markov chain Monte Carlo methods which, in our case, are not computationally feasible approaches. Instead, in this paper, we compute approximations for the credibility intervals that are based on local Gaussian approximation of the posterior distribution at the MAP estimate. The computation of these approximations can be derived as follows. The forward model x → y = f (x) is approximated by the first order Taylor series at where J f (x MAP ) is the Jacobian matrix of f (x) evaluated at x = x MAP , which is obtained using the Gauss-Newton algorithm. Then, the Taylor series approximation is substituted into the observation model (19). By forming the mean and covariance matrix of the combined distribution of (x, y), and then using Schur complements to obtain the conditional distribution of x|y results in a Gaussian approximation for the posterior distribution, that is, x|y ∼ N (η, ), where the (approximate) posterior mean is η = x MAP and the posterior covariance is For full derivation, see [74], [75]. For true Gaussian distributions, 99.7% of the probability mass of each element In this paper, we refer to these intervals as posterior error estimates. Although these intervals do not exactly correspond to 99.7% of the respective probability mass, such posterior error intervals can be described as safe ones. The validity of the approximation is evaluated in the Appendix. Thus, we define the posterior error intervals at points r q as where σ a,q and σ s,q are the square roots of the respective diagonal elements of the covariance matrix .

B. Prior Model
In this work, the prior model for the unknown parameters μ a and μ s was chosen to be based on the Ornstein-Uhlenbeck process [45], [63], [83]. The prior is a Gaussian distribution with the covariance matrix defined as being proportional to matrix which has its elements defined as where i and j denote the row and column indices of the matrix, r i and r j denote the grid node coordinates, and ξ is the characteristic length scale of the prior describing the spatial distance that the parameter is expected to have (significant) spatial correlation for. The random field described by the covariance (28) is independent of the absolute position, with the spatial correlation between any two points r i and r j being only dependent on the distance between the points, and the characteristic length scale. The characteristic length scale is typically chosen to be descriptive of the expected size of inclusions in the target of interest. In this work, the mean of the priors for μ a and μ s were chosen as and the covariances of the priors were chosen as where min μ a , max μ a , min μ s , and max μ s denote the minimum and maximum values that μ a and μ s are expected to vary between, with the minimum and maximum being one standard deviation away from their midpoint. The combined prior parameters were then These choices correspond to prior assumption that μ a and μ s are mutually uncorrelated. In biological tissues, this assumption can be inaccurate in the sense that close to tissue boundaries it is reasonable to expect that both of the optical properties could change. This would imply some correlation between the optical parameters, however, the uncorrelation assumption of the optical properties promotes more generic solutions to the inverse problem and is thus used in this work.

IV. SIMULATION STUDIES
The approach was studied using 2D and 3D simulations. The inverse problem of estimating μ a and μ s using the models described in Section II and the reconstruction method described in Section III was investigated. In 2D, the acoustic sensor configuration, optical illuminations, and the noise level were varied. In 3D, a limited view case with acoustic measurement performed only from one side of the target was investigated.
For the inverse problems, the MAP estimates (22) were computed using the Gauss-Newton algorithm (23). Gauss-Newton algorithm started from initial value of x 1 = η x . Accuracy of the estimates were computed using relative errors where μ TRUE a and μ s TRUE are the true parameters used to compute the simulated measurement data and μ MAP a and μ s MAP are their respective MAP estimates, obtained with (22).

A. 2D Simulations
The 2D simulations were performed in square domain = [−5 mm, 5 mm] × [−5 mm, 5 mm]. Two different illuminations (M = 2) were used to generate the simulated measurement data. The illuminations were formed by functions with being combination of L , R , B , T corresponding to left, right, bottom and top edges of the square domain respectively. The illumination functions are listed in Table I. For the acoustic sensor locations, six different configurations were used: four densely positioned detector configurations and two sparsely positioned detector configurations. The densely positioned detector configurations were: four sided detection geometry with sensors on L , R , B , T (number of sensors is K = 128); three sided detection geometry with sensors on L , B , T (K = 96); two sided detection geometry with sensors on L , B (K = 63); one sided detection geometry with sensors on L (K = 31). The sparsely positioned detector configurations were: sensors on a circle with a radius of 7.8 mm (K = 31); and sensors on a semicircle with a radius of 7.8 mm (K = 15). The sensor configurations are shown in Figures 1 and 2 with red dots.
The data was simulated in a grid composed of Q = 1369 grid nodes formed by square elements with side length of 277.8 μm. The distribution of μ a and μ s parameters are shown in Figure 1 on the top row. The parameter values were chosen to be within biologically relevant variation ranges of fat and blood [84]. For the sound speed, value c = 1500 m/s was used, which is similar to the speed of sound in water and soft tissues. The temporal pressure waveforms were sampled at 10 MHz and discretized into T = 189 temporal points at each of the acoustic sensor location. This corresponds to a recorded temporal pressure time series duration of 18.9μs and an acoustic propagation distance of 28.4 mm. The temporally Gaussian light pulse g τ had full width at half maximum of τ = 333 ns. Normally distributed noise, with zero mean and standard deviation proportional to the peak-to-peak pressure amplitude of the detected acoustic time series, was added to the simulated measurement data. This was done for the data obtained with the two different illuminations separately. The standard deviation was thus defined as where is the noise level and m ∈ {1, 2}. Values of = 0.05, 0.01 and 0.001 were used corresponding to 5%, 1% and 0.1% noise levels, or 13, 20 and 30 dB signal-to-noise ratios (SNR) respectively. In the inverse problem, a grid composed of Q = 961 grid nodes formed by square elements with side length of 333.3μm was used. For the noise statistics, the following values were used: where σ 1 and σ 2 are the standard deviation of the added noise corresponding to the two illuminations used, and I ∈ R Q×Q is the identity matrix. These values correspond to accurately knowing the statistics of the noise. The prior was assembled according to (31), with the characteristic length scale ξ = 1 mm. For the expected range of variation, the true minimum and maximum of the estimated parameters were used resembling a good prior knowledge.

B. 3D Simulations
The  Table I with L , R , B , T , D , V corresponding to left, right, bottom, top, dorsal (rear) and ventral (front) faces of the cube. K = 1089 acoustic sensors were located on L face of the cube set in square 33 × 33 sensor array.
The 3D simulation domain was discretized into cube elements with side length 277.8μm. This resulted in a simulation grid composed of Q = 50653 grid nodes, which was used to compute the simulated measurements. The temporal data was sampled at 10 MHz and discretized into T = 128 temporal points, corresponding to a recorded temporal pressure time series duration of 12.7μs and an acoustic propagation distance of 19.1 mm. The same full width at half maximum, as in 2D, of τ = 333 ns was used for g τ . The parameters used for μ a and μ s to compute the data are shown in Figure 6. Noise was added in the same fashion as in 2D with the noise level = 0.001 corresponding to 0.1% noise level. In the inverse problem, the simulation grid was composed of Q = 29791 grid nodes. The statistics of the noise mean and covariance matrix were chosen as in (35), and the prior was assembled according to (31), with the characteristic length scale ξ = 1 mm. For the expected range of variation, the true minimum and maximum of the estimated parameters were used.

V. RESULTS
2D MAP estimates of μ a and μ s are shown in Figure 1 for the data with noise level at 1%. Illuminations s 2D LR and s 2D BT were used. Results of four simulations where the acoustic sensors were located on the sides of the square domain are shown. Visually the estimates of μ a look comparable to each other even in the case when the acoustic sensors have located on only one side of the imaged target. In the one sided sensor configuration, however, some artefacts in the estimate are evident and visually look as if streaming away from the sensors. Larger differences between the estimates of μ s obtained with different acoustic sensor configurations are evident, but even in the case of one sided acoustic sensor configuration the estimate resembles the true μ s . 2D MAP estimates of μ a and μ s are shown in Figure 2 for sensor configuration where the sensors were located on a circle or semicircle. The noise level was 1% and illuminations s 2D LR and s 2D BT were used. In the case of the circle, the number of acoustic sensors was equal to that used in the one sided configuration. Visually the estimates of μ a and μ s are closer to the true parameters and contain less artefacts in comparison to one sided acoustic configuration. The estimate of μ s is closer in amplitude to the true value as well. In the case of semicircle, degradation of the visual quality of the estimate is evident in comparison to full circle configuration. However, although the semicircle configuration has less sensors than the one sided configuration in Figure 1, the estimate is better in visual quality. Figure 3 shows 2D MAP estimates of μ a and μ s for noise levels 5%, 1% and 0.1% for two sided acoustic sensor configuration. Illuminations s 2D LR and s 2D BT were used. As it can be seen, at high noise levels, the estimates of μ s are smooth and blurry, whereas at the low noise level, the estimate becomes sharper. The effect suggests, that at high noise levels, information about the high spatial frequency components of μ s is lost. Estimates of μ a , on the other hand, appear to become polluted by spatially high frequency noise when the noise level is high. Figure 4 shows the posterior error intervals of the MAP estimates of Figure 3 on a diagonal of the domain , computed with (27). As the noise level decreases, the estimates become more accurate, the posterior error intervals decrease and, in particular, the point estimates are consistent with the error estimates. For the absorption coefficients, the posterior error intervals also reflect the intuitive hypothesis that the error intervals should be larger at points which are further away from the measurement locations. Also, as expected, the error intervals for the scattering coefficient are larger than for the absorption coefficient, and they decrease slower with decreasing noise level. Regarding the width of the error intervals, we remind that we employ ±3σ intervals which are expected to contain almost 100% of the probability mass. Figure 5 shows 2D MAP estimates of μ a and μ s when number of illuminating sides has been varied. Estimates were computed with four side (s 2D LR and s 2D BT ), three side (s 2D LB and s 2D T ), and two side (s 2D L and s 2D B ) illuminations. Noise level was 0.1%. As the total illumination surface is reduced from four sided illumination to two sided illumination, the information in the areas far from the illuminating sides becomes poorer. This is reflected in estimates of μ a as becoming more distorted in those areas, as evident in the top right corner of the estimate produced by using illumination pair s 2D L and s 2D B . Estimate of μ s becomes overall worse everywhere as the total illuminating surface is reduced, with no clearly visible spatial location where the estimate is worst.
3D MAP estimates of μ a and μ s are shown in Figure 6. Illuminations s 3D LBD and s 3D RTV were used. The figure depicts situation where the acoustic sensors are located on one side of the imaged target. The estimates of both parameters work reasonably well, although they both suffer from imaging artefacts caused by the limited view. The inclusions (volumes of higher absorption or scattering in comparison to the ambient value) are visible in both absorption and scattering estimates. Table II shows the relative errors of the MAP estimates shown in Figures 1 and 2 for the six simulated 2D acoustic sensor geometries and one 3D simulation shown in Figure 6. For the 2D simulations the relative errors are shown for noise levels of 5%, 1% and 0.1%. Illuminations s 2D LR and s 2D BT were used in 2D and s 3D LBD and s 3D RTV in 3D. Based on the relative errors it is evident that the quality of the estimates of μ a and μ s improves when the noise level is reduced, or when the sensors are located in multiple directions of target. On the smallest error level of 0.1%, the circular and semicircular detector geometries resulted in lower relative errors of the estimates than in the one sided detector geometry.  Relative errors of μ a and μ s when the illuminating functions have been varied are shown in Table III for estimates shown in Figure 5. Based on the relative errors, it is evident that the estimates become worse when the total illuminating area is reduced or when the noise level is increased. In this work, a Bayesian approach to the direct estimation of the optical parameters from photoacoustic time series was described. The forward model was based on coupling the optical and the acoustic propagation models where the DA was used for the optical propagation and the acoustic model was formed by using Green's function of the homogeneous wave equation. The MAP estimates of the optical absorption and scattering distributions were computed.
The numerical 2D simulations show, that when the acoustic view of the imaged domain becomes limited, the estimates of the optical parameters degrade in terms of relative error. This was shown with simulations where acoustic sensors were located on four, three, two and one edge of the imaged square shaped domain. Similar effect was seen when the sensors were located on either a circle or a semicircle. Qualitatively (visually) the estimates of the optical parameters were good in each sensor setting. In 2D simulations, it was found that the estimates can be improved if the noise level of the data is reduced, or if the sensors are placed such that they enclose the imaged domain at least partially. This suggests that the quality requirements on measurement systems with a limited view acoustic sensor geometry are more demanding than in settings where the acoustic measurement can be done in more enclosing fashion. In 3D, one numerical simulation was performed with the acoustic sensors located on one face of a cube shaped imaged target. Although the relative errors of the estimated optical parameters were higher, in comparison to the similar setting in 2D, the estimates were qualitatively good with the distinct shape and quantitative value of the estimated parameters being similar to the true parameters.
The work presented in this paper takes quantitative photoacoustic tomography closer to being practically realizable imaging modality. In comparison previously published similar approaches that estimate the optical properties directly form the photoacoustic time series [68]- [73], this work combines the full nonlinear photoacoustic forward model into a Bayesian inverse problem approach, providing estimates of both the optical absorption and scattering in acoustically limited view measurement setting in three spatial dimensions. The approach enables inclusion of sophisticated prior models of the parameters of interest and statistical information on the noise and uncertainties related to the imaging modality into the estimation process. In addition, the approach is able to provide posterior error intervals to aid assessing the quality of the estimates.
In this work the numerical implementation of the photoacoustic model, and the Bayesian approach used to form the estimates of the optical parameters, was done in MATLAB (R2011b, The MathWorks Inc., Natick, MA). Time to form the MAP estimates in 2D was within few minutes for a workstation equipped with Xeon E5649 (Intel Corporation, Santa Clara, CA) processor operating at 2.53 GHz, whereas in 3D the estimation took few days on a computer with Xeon E5-2630 operating at 2.40 GHz (Intel Corporation, Santa Clara, CA). Increase in computational time, when shifting from 2D to 3D, is expected as the number of unknowns grows exponentially with the number of dimensions. Part of the large difference between the estimation times is explained due code level differences: some of the optimizations done in 2D (e.g. precomputation of some of the matrices) could not be performed in 3D due to memory limitations. Part of the reason why the memory consumption was an issue in 3D is explained by the use of direct (and possibly) naive form of implementing the Gauss-Newton algorithm: this involved explicitly assembling large matrices and then using direct approaches to solve for Gauss-Newton iteration direction. Additional work needs to be done, in order to make the estimation process faster.
The current study was limited to acoustically homogeneous and non-attenuating medium, however, in practice the assumptions are not necessarily valid and imaged targets are likely to be acoustically more complex. Thus, the presented approach could be taken forward by using more complete acoustic simulation models. On the other hand, the DA was used as the optical simulation model in this work. A more realistic (as well as numerically more challenging) photoacoustic simulation model would be based on the RTE [39]- [41], [43]- [45], [67]. The use of the RTE would make it possible to investigate the parameter estimation problem of QPAT, for example, when the optical scattering has anisotropic behaviour, or when highly directed light sources are used. Recently, a direct estimation of the optical absorption from the photoacoustic time series using a Tikhonov regularized least squares scheme and the RTE as the forward model has been investigated [72].
The Bayesian approach, used in this work, makes it possible to incorporate quantitative prior information into the inverse problem. The approach also enables modelling of complex error sources, such as might be arising from physical acoustic sensors or light sources. These issues could be handled by using, for example, the approximation error method [74], [85], which has been used in QPAT in [51], [55].
In conclusion, the numerical results shown in this work demonstrate that by coupling the optical and acoustic simulation models to a single photoacoustic model, it is possible to estimate the optical parameters from the photoacoustic time series even in a acoustically limited view setting. The results thus indicate that achieving quantitative photoacoustic tomography should be feasible in practical acoustic measurement geometries.

APPENDICES APPENDIX [FEASIBILITY OF THE NORMAL APPROXIMATION OF THE POSTERIOR]
When approximating a (posterior) density π(x|y) with a (normal) distribution π * (x|y), s.t. x|y ∼ N (η, ), the probabilities of events are generally never equal. Here, we define a safe approximation as one for which estimated errors based on the approximative marginal distributions are larger than those based on the actual posterior model. Thus, for example, we call the model safe if for all i = 1, . . . , 2Q, and some χ > 0, where σ i = √ (i, i ). We note that for nonlinear models, the asymptotic decay of the actual posterior is often slower than for a normal model and thus the condition (36) can fail for a large enough χ. The feasibility of the approximation depends on the underlying problem and the relevance of the choice of χ, or the choice of the error interval in general.
For normal approximations, the probabilities of parameter being within ±χσ intervals corresponding to χ = 1, 2, 3 are 68.3%, 95.5%, and 99.7% respectively. In order to compare the probabilities of the approximative posterior distribution derived error intervals to the normal distribution, a Monte Carlo approach is adopted.
Random samples of optical absorption μ a and μ s , reminiscent to the parameter distribution studied in the main text, were generated, denoted by x , = 1, . . . , 100. For each sample, the corresponding noisy data y was simulated using the observation model (19) with the same noise levels (5%, 1%, and 0.1%), illuminations (s 2D LR and s 2D BT ), and acoustic detection setup (detectors on two sides of the target), as used in results of Figure 4. The MAP estimates x MAP , and the normal distribution statistics η and of the approximative posterior distribution, were computed using the approach presented in Section III for each of the data y . From the covariance matrices , standard deviations σ were computed by taking the square root of the diagonal part of the matrices. Figure 7 shows a typical generated parameter distributions, their MAP estimates, and the ±3σ error intervals.
Based on the estimates x MAP , the true parameters x , and the standard deviations σ , sample based probability of x being within the error interval [x MAP −χσ , x MAP +χσ ] was computed for χ = 1, 2, 3. The probabilities were computed  by counting the proportion of discretization points, for which the true parameter was within χ standard deviations of the estimate, with respect to number of disretization points, for all , and then taking the average of the proportions. The probabilities were computed for μ a and μ s separately. The sample based probabilities are shown in Table IV. Based on the sample based probabilities, the normal distribution approximation of the posterior distribution, and the error intervals derived from the approximation, can be considered safe for both μ a and μ s , except when the noise is small. The reason that the approximation fails at small noise levels, is likely to be caused by not including modelling errors in the observation model (19): at small error levels, the numerical errors of the forward model should be considered as well. However, even when the noise level is as low as 0.1%, the error intervals derived from the approximate posterior statistics, are less than 1.2% away from being safe estimates of error for μ a . Thus, the normal distribution approximation of the posterior distribution, and the error intervals derived from it, can be considered safe.