APPROXIMATE MARGINALIZATION OF UNKNOWN SCATTERING IN QUANTITATIVE PHOTOACOUSTIC TOMOGRAPHY

. Quantitative photoacoustic tomography is a hybrid imaging method, combining near-infrared optical and ultrasonic imaging. One of the in- terests of the method is the reconstruction of the optical absorption coeﬃcient within the target. The measurement depends also on the uninteresting but often unknown optical scattering coeﬃcient. In this work, we apply the approximation error method for handling uncertainty related to the unknown scattering and reconstruct the absorption only. This way the number of unknown parameters can be reduced in the inverse problem in comparison to the case of estimating all the unknown parameters. The approximation error approach is evaluated with data simulated using the diﬀusion approximation and Monte Carlo method. Estimates are inspected in four two-dimensional cases with biologically relevant parameter values. Estimates obtained with the approximation error approach are compared to estimates where both the absorption and scattering coeﬃcient are reconstructed, as well to estimates where the absorption is reconstructed, but the scattering is assumed (incorrect) ﬁxed value. The approximation error approach is found to give better estimates for absorption in comparison to estimates with the conventional measurement error model using ﬁxed scattering. When the true scattering contains stronger variations, improvement of the approximation error method over ﬁxed scattering assumption is more signiﬁcant.


(Communicated by John Schotland)
Abstract. Quantitative photoacoustic tomography is a hybrid imaging method, combining near-infrared optical and ultrasonic imaging. One of the interests of the method is the reconstruction of the optical absorption coefficient within the target. The measurement depends also on the uninteresting but often unknown optical scattering coefficient. In this work, we apply the approximation error method for handling uncertainty related to the unknown scattering and reconstruct the absorption only. This way the number of unknown parameters can be reduced in the inverse problem in comparison to the case of estimating all the unknown parameters. The approximation error approach is evaluated with data simulated using the diffusion approximation and Monte Carlo method. Estimates are inspected in four two-dimensional cases with biologically relevant parameter values. Estimates obtained with the approximation error approach are compared to estimates where both the absorption and scattering coefficient are reconstructed, as well to estimates where the absorption is reconstructed, but the scattering is assumed (incorrect) fixed value. The approximation error approach is found to give better estimates for absorption in comparison to estimates with the conventional measurement error model using fixed scattering. When the true scattering contains stronger variations, improvement of the approximation error method over fixed scattering assumption is more significant.
1. Introduction. Photoacoustic imaging is a technique utilizing both optics and acoustics [41,9]. As the surface of a target is illuminated with a short near-infrared light pulse, it propagates through the medium while part of its energy is being absorbed. If the light pulse is short enough, the absorbed optical energy translates into excess pressure. The excess pressure field propagates as sound in ultrasonic frequencies which can be measured on the boundary of the target.
The inverse problem associated with quantitative photoacoustic tomography (QPAT) is to determine the optical properties of the target based on the surface measurements of ultrasound. The inverse problem can be thought of as being twofold. First one has to reconstruct the initial pressure distribution based on the measurements. This alone can provide useful information about the target. To get quantitative data of the target, one has to also reconstruct the optical properties. This can be done by treating the reconstructed initial pressure distribution as data for the optical reconstruction.
Optical reconstruction methods in QPAT often utilize the diffusion approximation (DA) of the radiative transfer equation (RTE) as the light propagation model in the measurement model relating the optical properties to the reconstructed pressure data. Although the reconstruction problem has been studied using the RTE [4,59,19,54,43,51], the DA can provide faster reconstruction approaches due to its simplicity. Model based inversion is almost exclusively used. The optical inversion is ill-posed and regularization methods are often applied to overcome the ill-posedness [24,25]. Further, reduction of ill-posedness can be achieved by obtaining measurements with multiple illuminations [6,52], or at multiple wavelengths in combination with spectral models of tissue properties [17,62,40,7]. Multiple measurements are often needed as well to eliminate non-uniqueness [17,5] of the reconstruction problem. Additional improvements can be obtained by performing diffuse optical tomography measurements along with photoacoustic measurements [60,62,63,42].
In situations where the scattering coefficient is known beforehand, reconstruction methods can assume fixed scattering and reconstruct only the absorption distributions [8,15,10]. The approach leads to potentially fast reconstruction methods, especially in the case of QPAT where number of estimated parameters is typically large due to the capability of the approach to produce high resolution images of the optical properties. However, the fixed scattering assumption will result in systematic errors if the value of scattering is inexact. The work presented attempts to alleviate the issue by taking into account the uncertainties caused by the fixed scattering assumption.
In this work the approximation error method [33,34] is applied to the optical inverse problem of QPAT for approximate marginalization over the unknown scattering coefficient [37]. The approach was used in [37] with application in diffuse optical tomography. In the approach the modelling error caused by using inaccurate fixed scattering is modelled as an additive noise process in the measurement model by using the approximation error method. This noise process is called the approximation error. It depends on the unknowns and its realization is obviously unknown. However, given the prior probability density models of the unknowns, approximate marginalization of the likelihood model can be carried out by using a Gaussian approximation of the joint statistics of the unknown parameters and the approximation error [37,35]. This way the number of unknowns in the inverse problem is significantly reduced, while taking the uncertainties of the fixed scattering approximation into account. Reconstructions with the approximation error method are compared to reconstructions with the conventional measurement error model where the modelling error caused by the fixed scattering assumption is neglected. Previously the approximation error method has been applied e.g. in diffuse optical tomography to compensate for errors caused by inaccurate modelling of the boundary shape [27], to perform model reduction from dense to sparse inversion mesh [3,36], as well as to perform model reduction from RTE to DA [53]. Recently applicability of the approximation error method was studied in reducing artifacts due to scalp blood flow in brain activation imaging [28].
In this paper the DA is used in the optical inverse problem of QPAT, although the methods presented can be utilized for other models as well, such as the RTE. It is assumed that the acoustical reconstruction has been performed and hence the initial pressure field is known but may be polluted by noise. Grüneisen coefficient, which relates the initial pressure distribution and the optical absorbed energy, is assumed to be a known constant. The data is simulated using both the DA and Monte Carlo method. The latter model is a statistical method which has been widely accepted to describe accurately the propagation of light in medium with scattering particles. Monte Carlo can be regarded as a sampling approach for the distribution described by the RTE, whereas the DA is essentially a low order deterministic approximation of the same physical model, that assumes a simplification of the distribution of photon propagation directions. MC is therefore a "gold standard" for simulated data and is utilized in this work to form an idea of how good reconstructions can be obtained using an inverse approach based on the DA applied on physically more realistic measurement data.
The paper is organized as follows. In Section 2 the physics associated with optical inversion of QPAT is reviewed. In Section 3 the Bayesian approach for reconstructing the optical properties, as well as the approximate marginalization scheme, are presented. In Section 4 description of the numerical implementation is given. In Section 5 the numerical results are presented and in Section 6 the paper is finalized with a recap of the results and conclusions.
2. Simulation and observation model. The photoacoustical phenomena is generated by absorption of short light pulse into the target. The absorbed energy translates into excess pressure p 0 in the tissue according to relation where G is the Grüneisen coefficient related to thermodynamical properties of the target and H is the absorbed energy density of the light pulse, both functions of the spatial coordinates r. The excess pressure p 0 propagates through the target as ultrasonic sound wave. Measurement of the transient ultrasonic pressure signal on the surface makes it possible to reconstruct the initial pressure field p 0 . The absorbed energy density is defined as where µ a is the optical absorption coefficient and Φ is the optical fluence. In the rest of the paper the spatial dependence will be omitted in the notation. For highly scattering soft-tissue the optical fluence can be described by the DA [13,30] (3) where κ = n(µ a + µ s ) −1 is the optical diffusion with n being the number of spatial dimensions and µ s is the reduced scattering coefficient. A physical boundary conditions for the DA is where A describes the reflectivity of the boundary, γ n is related to number of spatial dimensions with γ 2 = π −1 and γ 3 = 1/4, n is the boundary normal, and g describes the spatial light source on the boundary. In this paper the acoustical reconstruction of p 0 based on the boundary measurements of the ultrasound wave is assumed given. Further, the Grüneisen coefficient is assumed to be a known constant and the focus is on the discretized observation model where z = (z 1 , ..., z N ) T ∈ N are the measurements z k = H(r k ) at points r k obtained from the acoustical reconstruction, f (µ a , µ s ) : 2×M → N is the discretized forward model (2) , µ a , µ s ∈ M , and e ∈ N is additive noise left by the acoustical reconstruction.

Bayesian inversion in QPAT.
Idea behind Bayesian inversion is to handle all the parameters of the observation model (5) as random variables [33]. Given the model (5) the conditional probability density function of z given µ a , µ s and e is According to Bayes' theorem where distributions π(µ a ), π(µ s ) and π(e) are prior probability density models of µ a , µ s and e which have been modelled mutually uncorrelated. The uncertainty related to the unknown but uninteresting measurement noise can be marginalized by integrating (7) over e, leading to (8) π(z, µ a , µ s ) = π(z, µ a , µ s , e) de where the sub-index e is used to denote the probability density function of e. The solution of the inverse problem is the posterior density representing the complete statistical model for the inverse problem, given the data, likelihood and prior models. Inference about the solution of the problem can be made by computing point estimates from the posterior model. One commonly used estimate is the maximum a posteriori (MAP) estimate, which is defined as For Gaussian prior models µ a ∼ N (η µa , Γ µa ), µ s ∼ N (η µ s , Γ µ s ) and noise statistics e ∼ N (η e , Γ e ), where η are the mean and Γ the covariances, it holds that with the Cholesky decompositions Estimates computed with (11) will be referred to as conventional error model (CEM) estimates.

3.2.
Approximate marginalization over µ s by the approximation error method. Consider the case that one seeks to estimate µ a only, while using a fixed realization µ * s for the scattering. In the approximation error method the observation model (5) is then written as [37] (12) is the approximation error describing the modelling error caused by using the fixed µ * s . The objective is to marginalize the likelihood model over the unknown noise processess and e.
To obtain an approximate but computationally efficient solution, we approximate the density of by a Gaussian density. Following the approach in [37], the MAP estimate of µ a corresponding to the model (12) is obtained as where η and Γ are the mean and covariance of . Estimates formed using (13) will be referred to as approximation error method (AEM) estimates, In this paper sampling approach is used to obtain estimates for the statistics of . This is achieved by fixing teaching distributions of µ a and µ s , which are used to make random draws of the two parameters. Both distributions will be modelled Gaussian with µ a ∼ N (η µa , Γ µa ) and µ s ∼ N (η µ s , Γ µ s ). Although the teaching distributions are referred to with the same symbol for the mean and covariance as the prior distributions in Section 3.1, they need not be the same. In this paper the teaching distributions have different statistical parameters than the priors and the meaning of the notation should be clear from the context. Using a set of K random draws µ (i) a and µ (i) s of µ a and µ s from the teaching distributions, samples (i) of with i = 1, ..., K are computed as (14) ( and the statistics of are estimated as 3.3. Fixed scattering assumption. In addition to estimates computed with the CEM and AEM, estimates using the fixed scattering assumption using the conventional measurement error model are computed. The estimates were computed in order to form a reference of how much the AEM approach can improve reconstructions over the approaches where the unknown scattering has been assumed a fixed value, while the modelling error caused by this assumption is neglected. The estimates are computed using an observation model (17) z = f (µ a , µ * s ) + e, where µ * s is the fixed realization of the scattering coefficient. The MAP-estimate is then found to be These estimates will be referred to as fixed scattering assumption (FSA) estimates.
4. Implementation. where Φ = (Φ 1 , ..., Φ L ) T , the matrix elements are and the right hand side vector elements are where g and A are the parameters in the boundary condition (4).

4.2.
Prior distributions of µ a , µ s and e. In this work proper, or informative, smoothness priors are used for µ a and µ s . The prior distributions are used for the computation of the approximation error, as well as for solving the MAP-estimates. The prior used is described thoroughly in [3], but reviewed here as well. The prior is based on modifying a smoothing preprior, for some discretized parameter x, into a Gaussian distribution. The preprior has the form where x is a vector, α is a regularization parameter, D is a differentiating matrix and B = D T D. Typically D T D is not invertible, and hence cannot describe a Gaussian distribution (the covariance is not defined.) The preprior can be modified into a proper distribution by reordering x into x = (x 1 , x 2 ) T such that for nodes x 2 some quantitative prior information of the form x 2 ∼ N (η x 2 , Γ x 2 ) is attached. In this paper, nodes x 2 are called the marginal nodes. Partitioning matrix B similarly the preprior can be expressed as which is equal to π(x) since all that was done above was rearrangement of the nodes. For conditional distribution of x 1 on condition x 2 it holds that and for combined distribution of x (and thus x) In two-dimensions the differentiating matrix D used in this work is defined as where D x and D y are matrices computing the x-and y-coordinate directional derivatives at the finite element mesh. The subindex x here should not be confused with the parameter x for which the prior is defined. By adjusting the parameter α it is possible to tune the relative point-wise standard deviations of x 1 with respect to x 2 . In this paper α is chosen such that the marginal variances of the prior are equalized.
In assembling the informative smoothing priors for discrete µ a and µ s , constant mean is used for the marginal nodes with covariances of the form [36] where 1 is a matrix of ones with the same size of the identity matrix I, and σ background and σ inclusion are the presumed standard deviations of background variation and variation of inclusions respectively. Equation (31) (with the associated mean η) corresponds to Γ x 2 (and η x 2 ) in (29). In this work, the marginal nodes had a correlation length of about 1.3 mm. The statistics of the additive error e is modelled to be zero-mean and uncorrelated with η e = 0 and Γ e = σ 2 e I, where σ e is the expected standard deviation of the additive error.

5.
Results. Use of the approximation error method for the approximate marginalization was investigated in a two-dimensional square domain with spatially distributed optical properties. The DA and Monte Carlo (MC) were used to simulate the data. MC-model that is used has been described in [53].
Estimates using CEM in equation (11), AEM in equation (13) and FSA in equation (18) are presented and compared both visually and quantitatively by computing the relative errors of the reconstructions with respect to the true optical parameters by using where Euclidean norm is used and x and x true are the reconstructed and true optical absorption or reduced scattering coefficients.

Simulation parameters.
Four two-dimensional simulation cases, shown in Figure 1 were investigated: I) square inclusions, II) smooth inclusions, III) nonoverlapping inclusions (i.e. case where µ a and µ s have different spatial structure) and IV) high contrast smooth inclusions. In each case the domain was [−5, 5] × [−5, 5] mm. The parameter values were chosen for cases I-III such that the background of the material parameters corresponded to adipose breast tissue (fat) [14] and the inclusions corresponded to blood with different oxygen saturation levels [44,45]. Case IV was constructed based on case II by scaling the reduced scattering coefficient of the background to a higher value. For each test case, four measurements were simulated in a dense triangular mesh, composed of 7820 triangular elements and 3747 grid nodes. The four measurements were computed with illuminations g = g 1 , ..., g 4 in (4) corresponding to incident light coming from left, bottom, right and top of the simulation domain defined as (33) For the MC-simulations Henyey-Greenstein scattering anisotropy parameter of zero, corresponding to value of the scattering coefficient of µ s = µ s was used.
The simulated data were interpolated into a sparser inversion mesh. The inversion mesh was composed of 5022 triangular elements and 2598 grid nodes. After the interpolation, zero-mean Gaussian random noise was added to the simulated data to form noisy measurements. The standard deviation σ e of the noise was defined such that 3σ e limits of the noise equaled to 1% of the peak amplitude of the simulated (noiseless) data.  Table 1. Values used for the reconstruction priors and the teaching distributions of µ a and µ s . η is the mean, σ background and σ inclusion are the standard deviations for the variation of background level and the inclusions. Shown are also the range in which 95% of the samples fall in.
Reconstruction Parameters used for the reconstruction priors and the teaching distributions are shown in Table 1. For the standard deviations of the background, value of σ background = 0.1 · σ inclusion was used. The teaching distributions were otherwise similar to reconstruction priors, except that they were narrower and had lower mean value.
The fixed scattering values µ * s used in AEM and FSA was chosen to be the mean of the teaching distribution of µ s , i.e. µ * s = η µ s . All the MAP-estimates were computed by using Gauss-Newton method with a line-search algorithm.

5.3.
Reconstructions using DA-data. Reconstructions of µ a when using DAsimulated data are shown in Figure 2. CEM reconstruction of µ s is shown in Figure  3. Reconstructions of µ a show that the absorption coefficient is generally well reconstructed, while larger errors are clearly visible in the reconstructions of µ s . Improvement of AEM reconstruction over FSA are clearly visible in the simulation case IV, where the true µ s has a high contrast. Table 2 shows the relative errors of the estimates. The benefit of approximate marginalization is evident when comparing AEM reconstructions to FSA as the relative errors are smaller in the approximation error approach in all the simulation cases. Notice also that the relative errors of the AEM are only slightly larger than the errors of CEM reconstructions for the simulation cases I-III. For the simulation case IV there is a larger discrepancy between the approaches, with AEM approach giving significantly better estimates than FSA.  Figure 4 shows the relative error of AEM estimates when the mean η µ s of the teaching distribution of µ s has been varied for the cases II and IV. The estimates were computed with three different covariances of the teaching distributions: a wide distribution with twice the standard deviations, the same standard deviations, and a narrow distribution with half the standard deviations compared to the values in Table 1. The figure also shows the relative error of the FSA estimates with the fixed reduced scattering µ * s set to the same value as in the AEM teaching distribution,  i.e. µ * s = η µ s . As evident the relative errors of AEM estimates are similar regardless of the standard deviation of the teaching distributions, especially when η µ s is close to the range of true scattering. The level of the relative error increases clearly when the mean of the teaching distribution is taken far away from the range of true scattering values. However, improvement of the AEM over the FSA estimates is found throughout the whole range of values of η µ s in case II, whereas in case IV FSA gives slightly better estimates than the AEM at high values of η µ s . E µ a (%) Figure 4. Relative error of µ a using the AEM-approach, when the mean of the scattering η µ s of the teaching prior has been varied for the case II (top) and for the case IV (bottom.) Data is shown for a teaching distribution with a wide distribution with double standard deviations (dash dotted black line), a distribution with the same standard deviations (solid black line), and for a narrow distribution with half standard deviations as in Table 1 (dashed black  line). Thick grey line shows the relative error of FSA estimates with µ * s = η µ s . The horizontal line denotes the mean ± standard deviation limits of the true scattering coefficient and circles denote the minimum and maximum values of the true scattering coefficient.  Figure 5. Figure 6 shows the CEM reconstructions of µ s . Overall the reconstructed µ a has the correct shape but it fails to match in absolute values, which is especially evident from the brightness of the inclusions in the simulation cases I and III. This is partially explained by the differences in the physics of the DA-model used in the reconstructions and the MC-model used to generate the measurements. Issues might also be caused by the size of the simulation domain with respect to the material parameters used. As was observed in [54] it might be better to utilize the RTE as the reconstruction model in the parameter scales used in the simulations. The µ s reconstructions are visibly worse than in Section 5.3. The AEM and FSA reconstructions are able to capture the general shape of µ a in cases I-III. Whereas the AEM gives good estimates also in the simulation case IV, the FSA does not. Table 3 shows the relative errors of the reconstructions when the measurements are simulated with MC. The AEM reconstructions are better than the FSA reconstructions in the simulations cases I, II and IV. However in the simulation case III the FSA reconstruction yielded the best of reconstructions approaches. Surprisingly the AEM estimates are better than the CEM estimates in cases I-III. 6. Conclusions. In this work the approximation error method was applied to the inverse problem associated with the optical reconstruction in QPAT. The statistics of approximation error caused by using fixed incorrect reduced scattering coefficient was approximated as a normal distribution. This in turn was used to reduce the number of parameters in the inverse problem by carrying out approximate marginalization over the approximation error. The approach was applied to the DA-model of the propagation of light in scattering media with biologically relevant optical parameters.
By using two-dimensional simulated data it was shown that the approximation error method can be used to reduce the number of reconstructed parameters of the inverse problem. Improvement of estimates using the approximation error method over the estimates obtained by a fixed scattering assumption was evident based on the relative errors of the estimates. This shows that the modelling error caused by the fixed scattering assumption should be taken into account when performing reconstructions with fixed scattering. The method was also tested for Monte Carlo simulated data. Larger discrepancies were observed in the estimates, although the approximation error method was beneficial in comparison to the fixed scattering estimates in most of the simulation cases. Reasons for the discrepancies arise from the physical differences between the diffusion approximation and the Monte Carlo model. Whereas diffusion approximation is a good approximation of propagation of light in highly scattering media, the Monte Carlo works well even when the scattering is not strong. The results would indicate that some criticism should be placed on the practical usability of the diffusion approximation based inversion schemes in quantitative photoacoustic tomography.
Although this work investigates the use of approximation error method to marginalize the unknown scattering in quantitative photoacoustic tomography in two dimensions, the presented approach can be extended into three dimensions. It is expected that by taking into account the uncertainties caused by the fixed scattering assumption will result in better estimates than when the uncertainties are neglected.