2 December 2016 Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance
Author Affiliations +
J. of Biomedical Optics, 21(12), 126004 (2016). doi:10.1117/1.JBO.21.12.126004
Abstract
Forward and adjoint Monte Carlo (MC) models of radiance are proposed for use in model-based quantitative photoacoustic tomography. A two-dimensional (2-D) radiance MC model using a harmonic angular basis is introduced and validated against analytic solutions for the radiance in heterogeneous media. A gradient-based optimization scheme is then used to recover 2-D absorption and scattering coefficients distributions from simulated photoacoustic measurements. It is shown that the functional gradients, which are a challenge to compute efficiently using MC models, can be calculated directly from the coefficients of the harmonic angular basis used in the forward and adjoint models. This work establishes a framework for transport-based quantitative photoacoustic tomography that can fully exploit emerging highly parallel computing architectures.
Hochuli, Powell, Arridge, and Cox: Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance

1.

Introduction

Quantitative photoacoustic tomography (PAT) is concerned with recovering quantitatively accurate estimates of chromophore concentration distributions, or related quantities such as optical coefficients or blood oxygenation, from photoacoustic images.1 The source of contrast in PAT is optical absorption, which is directly related to the tissue constituents. By obtaining PAT images at multiple optical wavelengths, it may be possible to recover chemically specific information about the tissue. However, such a spectroscopic use of PAT images must consider the effect of the spatially and spectrally varying light fluence distribution. As a photoacoustic image is the product of the optical absorption coefficient distribution, which carries information about the tissue constituents, and the optical fluence, which only acts to distort that information, the challenge in quantitative photoacoustic imaging is to remove the effect of the light fluence.

A common approach is to use a model of the unknown fluence and use it to extract the desired optical properties from the measured data. This has been done analytically23.4.5 or numerically,6,7 often within a minimization framework.89.10.11.12.13.14.15 The majority of this literature uses the diffusion approximation to the radiative transfer equation (RTE) to model the light distribution, which is accurate in highly scattering media.16 In PAT, the region of interest often lies close to the tissue surface where the diffusion approximation is not accurate. The RTE, on the other hand, is widely considered to be an accurate model of light transport so long as coherent effects are negligible, which is the case here. Finite element discretizations of the RTE have been developed1718.19.20.21 and proposed for quantitative PAT reconstructions,12,15 but due to the need to discretize in angle as well as space, they quickly become computationally intensive and their applicability is limited to small- and medium-scale problems. An alternative is Monte Carlo (MC) modeling,22,2324.25 which is a stochastic technique for modeling light transport that converges to the solution to the RTE. The significant advantage of the MC approach is that it is highly parallelizable, thus it scales well to the large-scale inversions that will be encountered in practice.

MC models of light transport are popular in biomedical optics and have predominantly been applied in the planning of the experimental measurements2627.28 and in dosimetric studies for a range of light-based therapies.2930.31 Many of the applications are summarized by Zhu et al.32 One early MC model of light transport, MCML,22 computes the fluence in three-dimensional (3-D) slab geometry. This model was later extended to simulate spherical inclusions in the tissue,33 and later to spheroidal and cylindrical34 inclusions. MC modeling in 3-D heterogeneous media has been shown both for voxelized media,23 which was later GPU-accelerated,24 and using a mesh-based geometry.25,35 Although the RTE is an equation for the radiance, which is a function of angle at every point, the quantity usually calculated by MC models is the fluence rate, which is the radiance integrated over all angles. The reasons are practical: most measurable quantities are related to the fluence rate rather than the radiance, storing just the integrated quantity saves on computational memory, and the estimates for the fluence rate will converge sooner than the underlying estimates for the radiance. In photoacoustics, the measurable signal is related to the fluence (the time-integrated fluence rate) so current MC models can be used in the simulation of photoacoustic signals. However, as will be discussed in Sec. 4, the full angle-dependent radiance is required when tackling the inverse problem of estimating the optical coefficients, specifically the optical scattering. A common approach in estimating the optical coefficients is to use a gradient-based inversion scheme in which functional gradients are used to minimize some objective function. This paper demonstrates the computation of these gradients through the use of forward and adjoint MC models of the radiance.

In this paper, Sec. 2 introduces the inverse problem of quantitative PAT. Sections 3 and 4 present forward and adjoint MC models of the radiance employing a harmonic angular basis. In Sec. 5, it is shown that this choice of basis allows the functional gradients for the inverse problem to be calculated straightforwardly. Inversions for absorption and scattering coefficient distributions are given in Sec. 6.

2.

Quantitative Photoacoustic Tomography

The inverse problem in QPAT can be stated as the minimization

(1)

argminμa(x),μs(x)ε[μa(x),μs(x)],
where the error functional is given by

(2)

ε=12Ω[Hmeas(x)H(x;μa,μs)]2dx.
H=μa(x,λ)Φ(x,λ;μa,μs) is the absorbed energy density and is the “data” for this problem. It is related to the photoacoustic image by the Grüneisen parameter, which here is set to 1. Additional regularization terms or terms reflecting prior knowledge may also be added to ε. Gradient-based approaches to solving this problem require estimates of the gradients of the error functional with respect to the parameters of interest. Saratoon et al.12 give expressions for these gradients in terms of the forward and adjoint fields, φ and φ*:

(3)

εμa=Φ(HmeasH)+Sn1φ*(s^)φ(s^)ds^,
and

(4)

εμs=Sn1φ*(s^)φ(s^)ds^Sn1Sn1φ*(s^)P(s^,s^)φ(s^)ds^ds^.
In the following two sections, the calculation of the forward radiance φ(s^) and adjoint radiance φ*(s^) using MC models is shown.

3.

Monte Carlo Modeling of Light Transport

In PAT, the optical and acoustic propagation times are so different that the optical propagation can be considered instantaneous and the time-dependence of the light transport can be neglected. The time-independent RTE is given by

(5)

Lφ=q[s^·+μa(x)+μs(x)]φ(x,s^)μs(x)Sn1P(s^,s^)φ(x,s^)ds^=q(x,s^),
where L is the forward operator of the RTE and φ is the radiance, μa and μs are the absorption and scattering coefficients, respectively, x is position, s^ and s^ are the original and scattered propagation directions, P(s^,s^) is the scattering phase function, q(x,s^) is a source term, and Sn1 is used to indicate integration over angle in n1 dimensions. To obtain approximations to the solutions to this equation, various flavors of MC have been proposed.36 The approach used here begins with launching a packet of energy, referred to herein as a “photon,” from a given position x in an initial direction s^. After traveling a distance s=U([0,1])/μs (using the convention s=|s|s^), where U([0,1]) is a real uniform random variable on [0, 1], a fraction of the photon’s “weight” W[1exp(μas)] is deposited in the current voxel, where W is the current weight (or energy) of the photon packet. The photon weight is updated: WWexp(μas). Scattering into a new direction s^ in two-dimensional (2-D) involves sampling the scattering phase function, which describes the probability of a photon scattering from direction s^ into direction s^. The phase function used here was the 2-D Henyey–Greenstein phase function, commonly used in biomedical optics,37,38

(6)

P(s^,s^)=12π1g2[1+g22g(s^·s^)].
The parameter, g, a property of the medium, is known as the anisotropy factor. Sampling this equation for the scattering angle, θ=arccos(s^·s^), by solving for θ in the cumulative integral over angle yields

(7)

θ=2arctan{1g1+gtan[πU([0,1])]}.
A new step length, s, is sampled and this process is repeated until the photon weight falls below some threshold value. Carrying out the above computation for many photons and adding the voxel weights will, for a sufficient number of photons, converge on a solution to the RTE.

By calculating photon paths through the medium, the MC models presented in the literature2223.24,35 do in fact simulate the radiance, but this is typically integrated over angle upon deposition of the weights in the voxels. In order to simulate the radiance, a method of depositing the weight in the voxels without losing the angular information is required.

3.1.

Monte Carlo Modeling of the Radiance

In order to compute the radiance using an MC model, angular as well as spatial discretization is required. One approach is to use discrete ordinates, whereby the unit circle is divided equally into sectors and the weight deposited in a voxel is also assigned to the relevant angular sector. The memory required will scale linearly with the number of sectors, and will slow convergence of the radiance estimate, compared with the fluence estimate, by a factor inversely related to the number of sectors. Here, a harmonic angular basis was used because a sufficiently diffuse field is dense in such a basis, meaning the field can be represented using relatively few orders. Less memory will, therefore, be required.

In 2-D, the expansion for the radiance in a Fourier basis is39

(8)

φ(x,θ)=12πa0(x)+1πn=1N=an(x)cos(nθ)+1πn=1N=bn(x)sin(nθ),
where an and bn are the coefficients associated with each harmonic and θ[π,π] is the angle of the photon direction s^ relative to the z-direction [i.e., θ=arccos(s^·z^)]. (The equivalent expansion in 3-D would be into spherical harmonics.40) For a given voxel, the weight is deposited into the relevant Fourier coefficients according to

(9)

a0=np=1NpdWnpS1δ(θθnp)dθ=np=1NpdWnp,

(10)

an=np=1NpdWnpS1δ(θθnp)cos(nθ)dθ=np=1NpdWnpcos(nθnp),

(11)

bn=np=1NpdWnpS1δ(θθnp)sin(nθ)dθ=np=1NpdWnpsin(nθnp),
where dWnp is the weight deposited by the npth photon traversing a voxel and θnp is the direction of the npth photon relative to z^. The radiance MC algorithm, or RMC, was implemented in the Julia programming language.41

3.2.

Validation of the Forward Model

Analytical solutions to the RTE are available for the fluence for a range of geometries and source types;23,42,43 however, there are few analytical solutions for the radiance, particularly in 2-D. The RMC model was compared to one such analytic solution for an infinite, homogeneous 2-D domain illuminated by an isotropic point source.4445.46 An isotropic point source was placed at the center of a domain of size 15  mm×15  mm, large compared to the transport mean free path in order to approximate an infinite domain. The absorption and scattering coefficients were 0.01 and 10  mm1, respectively, and the Henyey–Greenstein phase function47 was used with g set to 0.9. The pixel size was 0.05  mm×0.05  mm, and five Fourier harmonics were used. Figure 1 shows the good agreement between the analytical and RMC modeled radiance at radial distances of 2 and 3 mm from the source along the horizontal axis.

Fig 1

Polar plots of the angle-resolved radiance due to an isotropic point source in a homogeneous domain with μa=0.01  mm1, μs=10  mm1, and g=0.9. Results from an analytic method (infinite domain) and RMC simulations (15  mm×15  mm square domain) shown.

JBO_21_12_126004_f001.png

4.

Adjoint Monte Carlo Model

The adjoint equation to the RTE is given by

(12)

L*φ*=q*[s^·+μa(x)+μs(x)]φ*(x,s^)μs(x)Sn1P(s^,s^)φ*(x,s^)ds^=q*(x,s^),
where L* is the adjoint operator of the RTE and φ* is the adjoint radiance and q* is the adjoint source. This was implemented numerically using the same MC scheme as for the forward RMC model (Sec. 3.1). The principle difference is that the light sources q typically used in PAT are restricted to the boundary, but the adjoint source q* will not be, as a consequence of the fact that the “data” in QPAT—the photoacoustic images—is volumetric. The internally distributed sources, q*, were simulated by uniformly distributing the initial positions of each photon over a source pixel, with their direction randomly sampled from uniform angular distribution. The weight of each photon launched from a source pixel was scaled by μa(r)[Hmeas(r)H(r;μa,μs)], where r is the position corresponding to the source pixel. This of course can result in photons whose weight is negative. The adjoint RMC model treats a negative weight photon in the same way as one whose weight is positive; the photon’s weight decays to zero and rather than the deposition of energy into the computational grid, thus negative weight photons produce a reduction in the energy in each voxel it traverses.

4.1.

Validation of the Adjoint Model

The adjoint model was validated by checking that it satisfied the condition

(13)

Gu,v=u,G*v,
where G and G* are the (Green’s) operators solving the corresponding forward and adjoint RMC models, and u and v are the angle and position dependent source and detector, i.e.,

(14)

φ(r,s^)=GuLφ=u,φ*(r,s^)=G*vL*φ*=v.
Let u=ρs(r)Θs(s^), v=ρd(r)Θd(s^) and consider three different cases
Case1:(isotropic source and detector)  u1=δ(rrs)/2π,v1=δ(rrd)/2π.Case2:(isotropic source,angular detector)  u2=δ(rrs)/2π,v2=δ(rrd)Θd(s^).Case3:(spatial source,angular detector)  u3=ρs(r)Θs(s^),v3=δ(rrd)Θd(s^).

Substituting these into Eq. (13) yields

(15)

Φ1(rd)=Φ1*(rs),

(16)

2πφ2(rd,s^)Θd(s^)ds^=Φ2*(rs),

(17)

2πφ3(rd,s^)Θd(s^)dr=2πΩφ3*(r,s^)ρs(r)Θs(s^)ds^,
where φ1,2,3 and φ1,2,3* are the forward and adjoint radiances from computing Gu1,2,3 and G*v1,2,3, respectively. Φ is the fluence or angle-integrated radiance. It can be seen from Eq. (15) that case 1, where a pair of isotropic δ-functions are used for u1 and v1, that we expect the resulting fluence values at their respective positions, Φ1(rd) and Φ1*(rs), to be equal. This is an intuitive result given the reciprocity of the RTE and the angular independence of the source-detector combination.

Simulations were performed using a 40  mm×40  mm (101×101  pixel) domain, and 10 Fourier harmonics. Each source distribution emitted 106 photons. rs was set to be the center of the domain with rd moved along the x-direction across the domain. Comparisons are shown in Fig. 2 for case 1 with an isotropic source and detector, Fig. 3 for case 2 with an isotropic source and anisotropic detector with ρd=δ(rrs), Θd=1πsin2(2θ), and Fig. 4 for case 3 with the same ρd, Θd but with istropic Θs and the distributed Θs(r) shown in Fig. 4(a). Good agreement was obtained in all cases, showing the the RMC adjoint model is an accurate representation of the RTE adjoint. It can be seen in the comparison of Gu1,2,3,v1,2,3 with u1,2,3,G*v1,2,3 that some noise is present in the latter. As the overall number of photons simulated in the forward and adjoint simulations was the same, the use of spatially distributed sources in the adjoint case resulted in lower photon density in the domain, thereby reducing SNR. This issue was exacerbated in the case where an angularly dependent detector was used because a significant fraction of photons went undetected.

Fig. 2

Plot of Lu1,v1 and u1,L*v1 to validate the adjoint model. u1 and v1 were isotropic point sources with u1 at the center of the domain and v1 translated across the domain at y=23.6  mm.

JBO_21_12_126004_f002.png

Fig. 3

(a) Polar plot of source distribution for v2=δ(rrs)1πsin2(2θ); (b) plot of Lu2,v2=u2,L*v2 for validation of adjoint model. Plot was produced with u2 as an isotropic point source at the center of the domain. v2 was translated across the domain along a line at y=23.6  mm.

JBO_21_12_126004_f003.png

Fig. 4

(a) Isotropic source distribution u3=Ps(r); (b) Lu3,v3 and u3,L*v3 to validate the adjoint model. v3 was an anisotropic point source emitting light over angle following 1πsin2(2θ). v3 was translated along a line across the domain at y=23.6  mm, as shown by the gray line dashed line in (a).

JBO_21_12_126004_f004.png

5.

Functional Gradients

Both the radiance and the adjoint radiance can be expressed as Fourier series as in Eq. (8). By substituting these expressions into Eqs. (3) and (4) for the functional gradients, simple and easily computed expressions for the gradients can be obtained. The fluence is simply given by the isotropic component of the field a0. The other terms in the expressions for the functional gradients contain integrals of products of the radiance and its adjoint. If a0*, an*, and bn* are the Fourier coefficients of the adjoint radiance, then the gradient with respect to absorption can be written as

(18)

εμa=Φ(HmeasH)+S1φ(s^)φ*(s^)ds^=a0(Hmeasμaa0)+2π[14π2a0a0*+12π2a0m=1am*cos(mθ)+12π2a0m=1am*sin(mθ)+12π2a0*n=1ancos(nθ)+n=1m=1anam*cos(nθ)cos(mθ)+n=1m=1anbm*cos(nθ)sin(mθ)+12π2a0*n=1bmcos(mθ)+n=1m=1am*bnsin(nθ)cos(mθ)+nmbnbm*sin(nθ)sin(mθ)]dθ.
By orthogonality, all terms for which nm integrate to zero and Eq. (18) reduces to

(19)

εμa=a0(Hmeasμaa0)+2π[14π2a0a0*+1π2n=1anan*cos2(nθ)+1π2n=1bnbn*sin2(nθ)]dθ,

(20)

=a0(Hmeasμaa0)+12πa0a0*+1πn=1anan*+1πn=1bnbn*.
This expression for the absorption gradient is computationally straightforward to evaluate due to the fact that it requires simply summing over products of Fourier coefficients already loaded in memory.

The second term in Eq. (4) is

(21)

Sn1Sn1φ*(s^)P(s^,s^)φ(s^)ds^ds^,
which contains the phase function given in Eq. (22) and can be expanded using a Fourier series in powers of g37

(22)

P(s^·s^;g)=12π+1πl=1glcos(lΔθ),
where Δθ=arccos(s^·s^). Thus we can write

(23)

2π2πφ(s^)Pθ(s^,s^)φ*(s^)ds^ds^=2π2π[12πa0+1πn=1ancos(nθ)+1πn=1bnsin(nθ)]{12π+1πl=0glcos[l(θθ)]}[12πa0*+1πm=1am*cos(mθ)+1πm=1bm*sin(mθ)]dθdθ,
where θ and θ are the angles between the z-axis and s^ and s^, respectively. As such, the scattering angle between the previous direction s^ into the new direction s^ is given by (θθ). It is possible to expand cos[l(θθ)] as cos(lθ)cos(lθ)+sin(lθ)sin(lθ) which in turn allows us to employ orthogonality relationships to simplify the above integrals and write

(24)

S1S1φ(s^)Pθ(s^,s^)φ*(s^)ds^ds^=12πa0a0*+1πn=1anan*gn+1πn=1bnbn*gn.
Substituting this expression into Eq. (4), we can write the full expression for the functional gradient with respect to the scattering coefficient

(25)

εμs=12πa0a0*+1πn=1anan*+1πn=1bnbn*12πa0a0*+1πn=1anan*gn+1πn=1bnbn*gn

(26)

=1πn=1[anan*+bnbn*](1gn).
The ability to calculate these gradients is the first step to finding a computationally efficient way to solve the full QPAT inversion using an MC model of light transport.

6.

Inversions for Absorption and Scattering

The forward and inverse MC models of radiance described above were used with a gradient-descent (GD) scheme to estimate μa(x) and μs(x) from simulated PAT images by minimizing the error functional in Eq. (2). As the adjoint source, q*(x,s^)=μa(x)[Hmeas(x)H(x)] was independent of angle, photons were launched istropically with the launch position being spread out over the range of a source voxel using a randomly distributed number on the interval [0, 1]. The initial photon weight was scaled according to the source strength with normalization of the output quantity (i.e., radiance, absorbed energy density, and harmonic) being Np. The adjoint source may be negative in some places, so the initial photon weight is negative and weight deposition is also negative. The photon termination condition was, therefore, set to be the absolute value of the photon weight falling below the threshold value. The gradients were calculated using Eqs. (20) and (26). A GD scheme was chosen for the minimization because it is more robust to the MC noise in the functional gradients and error functional than techniques such as L-BFGS that use second-order information. A line-search algorithm presented by Hager and Zhang48 was used for the reconstruction of μa. A backtracking line search was implemented for the reconstruction of μs. This is described in Sec. 6.2.

The termination condition used by the optimization was

(27)

|ε(i)ε(i1)|/|ε(i)|<109,
where i is the iteration number. For all the reconstructions, it was assumed that the data Hmeas was given; no noise was added to the data, but MC noise from the forward simulation of the data was present at about 0.7% (evaluated by taking several runs of the forward model to estimate the average standard deviation over all positions across all model runs). For each inversion, the forward and adjoint RMC simulations used 108 photons and 10 Fourier harmonics, and was executed on a Dell 2U R820 32-core server.

6.1.

Inversion for Absorption Coefficient

The domain used in the estimation of the absorption coefficient consisted of a background absorption coefficient of 0.01  mm1 with a rectangular inclusion equal to 0.2  mm1 [shown in Fig. 5(a)], and a background scattering coefficient of 5  mm1 with a rectangular inclusion equal to 15  mm1 [shown in Fig. 5(b)]. The anisotropy was a homogeneously distributed value of 0.9. The measured data, Hmeas, were formed using an MC simulation illuminated by a collimated line source on the boundary at z=0  mm and on the adjacent boundary at x=4  mm, consisting of 108 photons. The inversion for the absorption coefficient only was performed under the assumption that the scattering coefficient was known and the starting estimate of the absorption coefficient was a homogeneous value of 0.01  mm1. The termination condition in Eq. (27) was satisfied after 11 iterations, having taken 4.1 h to run, and is shown in Fig. 5(b) with profiles through the true and reconstructed distributions of μa shown in Fig. 5(d).

Fig. 5

(a) True absorption coefficient, (b) true scattering coefficient, (c) reconstructed absorption coefficient after nine iterations, (d) profiles through true and reconstructed absorption coefficient at x=1.5  mm for all z.

JBO_21_12_126004_f005.png

Very good agreement between the true and reconstructed absorption coefficient is observed, with a value of the error function after 11 iterations being 2.9×109. The optimization routine was stopped as the change in the error function on the 12th iteration was below the function tolerance of 109, indicating convergence. The average error in the estimate of the absorption coefficient μaest, computed as |μatrueμaest|/μatrue, was 0.2% over the entire domain.

6.2.

Inversion for the Scattering Coefficient

Inversions for the scattering coefficient were performed using the same domain as above, shown in Figs. 6(a) and 6(b), with the illumination and number of photons in the forward simulation also being the same. Here, it was assumed that the absorption coefficient was known and the starting estimate of the scattering coefficient was equal to the background value of 5  mm1.

Fig. 6

(a) True absorption coefficient, (b) true scattering coefficient, (c) reconstructed scattering coefficient after 35 iterations, (d) profiles through true and reconstructed scattering coefficient at x=2.5  mm for all z.

JBO_21_12_126004_f006.png

Two modifications were necessary to achieve convergence in the optimization for the scattering coefficient. First, a custom GD algorithm was used in which a backtracking line search was implemented. A backtracking line search49 starts with a large candidate step length and progressively reduces the step size while checking for a sufficient decrease in the error functional. The sufficient decrease condition is expressed as

(28)

ε([μa(i),μs(i)]+α(i)p(i))ε([μa(i),μs(i)])+να(i)ε(i)Tp(i),
where ν was chosen to be 0.2 by inspection because this produced rapid convergence. In order to improve efficiency of the line search, step sizes were bounded between [105,109]; it was found that this range yielded sufficiently large steps to ensure reasonably efficient progress in the minimization. Second, the termination condition in Eq. (27) was relaxed due to the much slower convergence of the scattering coefficient, and instead required a relative change in the error functional of 105. This was satisfied after 35 iterations, and is shown in Fig. 6(c) with profiles through the true and reconstructed distributions of μs as shown in Fig. 6(d).

It can be seen from Figs. 6(c) and 6(d) that the inversion has partly reconstructed the inclusion in the scattering coefficient. The inability to reconstruct edges of the inclusion in the scattering coefficient is expected, given the diffusive nature of the scattering. However, the discrepancy in μs in the inclusion, evident from Fig. 6(d), suggests premature termination of the optimization. This is due to the fact that the gradient with respect to scattering is small and prone to noise in the functional gradients. This low SNR in the gradients has the impact that that search directions in the optimization routine are often suboptimal, which results in little or no progress of the optimization. The progressive reduction in SNR in the gradient means that nondescent steps are likely and can therefore trigger the termination condition.

7.

Discussion and Conclusions

In this paper, an MC model of the RTE was presented. The model computes the radiance in a Fourier basis in 2-D and is straightforward to extend to 3-D using a spherical harmonics basis. The accuracy of the model was demonstrated by comparing the angle-resolved radiance at two positions in the domain to corresponding appropriate analytic solutions.

Sections 5 and 6 presented the application of the RMC algorithm to estimating the absorption and scattering coefficients from simulated PAT images. In Sec. 6.1, it was observed that the absorption coefficient was estimated with an average error of 0.2% over the domain relative to the true value, when the scattering coefficient is known, and in the presence of 0.7% average noise in the data. This is encouraging, particularly because noise is not only present in Hmeas but also in the Fourier harmonics computed using the forward and adjoint RMC simulations, which is propagated to the estimates of the functional gradients. Consequently, the search direction in the GD algorithm will always be suboptimal. Furthermore, noise in H(μa(l),μs), the estimate of the absorbed energy density at the l’th iteration of the line search, will be also be propagated to the error functional, resulting in a nonsmooth search trajectory for the line search because at every point μa(l), the error function will be corrupted by some different noise σ(l):ε(μa(l),μs)=HmeasH(μa,μs)(1+σ(l))2. In practice, this did not preclude reconstruction of the absorption coefficient since the calculated gradients remained descent directions despite the noise. Furthermore, the error functional in μa is sufficiently convex that the addition of some noise does not prevent the linsearch from yielding sufficiently large a step length to allow rapid convergence.

Reconstruction of the scattering coefficent correctly located the scattering perturbation in the simulated image; however, the peak value in the reconstruction was lower than the true value. This is a direct consequence of the fact that the scattering coefficient is related to the absorbed energy distribution only through the optical fluence distribution. Consequently, the SNR in εμs is typically much less than that for absorption. This causes termination of the algorithm before the peak magnitude of the parameter has been found in the search space.

Disclosures

No conflicts of interest, financial or otherwise, are declared by the authors.

Acknowledgments

R. Hochuli was funded by the EPSRC CDT in Medical Imaging, EP/L016478/1. S. Powell was funded by EPSRC Grant No. EP/J021318/1, and RAEng Fellowship RF1516\15\33. S. Arridge and B. Cox are employed and funded by their host institution, University College London. The authors declare they have no conflicts of interest, financial or otherwise, relevant to this manuscript. The authors acknowledge the contribution of Andre Liemert who kindly provided radiance data used in the validation of RMC in Sec. 3.2. The authors acknowledge the use of the UCL Legion High Performance Computing Facility (Legion@UCL), and associated support services, in the completion of this work. The authors also thank Tanja Tarvainen for the useful discussions regarding the forward model and inversion.

References

1. 

B. Cox et al., “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt. 17(6), 061202 (2012).JBOPFO1083-3668http://dx.doi.org/10.1117/1.JBO.17.6.061202Google Scholar

2. 

G. Bal, A. Jollivet and V. Jugnon, “Inverse transport theory of photoacoustics,” Inverse Probl. 26, 025011 (2010).INPEEY0266-5611http://dx.doi.org/10.1088/0266-5611/26/2/025011Google Scholar

3. 

G. Bal and G. Uhlmann, “Inverse diffusion theory of photoacoustics,” Inverse Probl. 26, 085010 (2010).INPEEY0266-5611http://dx.doi.org/10.1088/0266-5611/26/8/085010Google Scholar

4. 

K. Ren, H. Gao and H. Zhao, “A hybrid reconstruction method for quantitative PAT,” SIAM J. Imaging Sci. 6, 32–55 (2013).http://dx.doi.org/10.1137/120866130Google Scholar

5. 

H. Ammari et al., “Mathematical modeling in photoacoustic imaging of small absorbers,” SIAM Rev. 52(4), 677–695 (2010).SIREAD0036-1445http://dx.doi.org/10.1137/090748494Google Scholar

6. 

R. J. Zemp, “Quantitative photoacoustic tomography with multiple optical sources,” Appl. Opt. 49, 3566–3572 (2010).APOPAI0003-6935http://dx.doi.org/10.1364/AO.49.003566Google Scholar

7. 

T. Harrison, P. Shao and R. J. Zemp, “A least-squares fixed-point iterative algorithm for multiple illumination photoacoustic tomography,” Biomed. Opt. Express 4, 2224–2230 (2013).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.4.002224Google Scholar

8. 

J. Laufer et al., “Quantitative determination of chromophore concentrations from 2D photoacoustic images using a nonlinear model-based inversion scheme,” Appl. Opt. 49, 1219–1233 (2010).APOPAI0003-6935http://dx.doi.org/10.1364/AO.49.001219Google Scholar

9. 

B. T. Cox, S. R. Arridge and P. C. Beard, “Estimating chromophore distributions from multiwavelength photoacoustic images,” J. Opt. Soc. Am. A Opt. Image Sci. Vision 26, 443–455 (2009).http://dx.doi.org/10.1364/JOSAA.26.000443Google Scholar

10. 

T. Tarvainen et al., “Bayesian image reconstruction in quantitative photoacoustic tomography,” IEEE Trans. Med. Imaging 32, 2287–2298 (2013).ITMID40278-0062http://dx.doi.org/10.1109/TMI.2013.2280281Google Scholar

11. 

E. R. Malone et al., “Reconstruction-classification method for quantitative photoacoustic tomography,” arXiv (2015).Google Scholar

12. 

T. Saratoon et al., “A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation,” Inverse Probl. 29, 075006 (2013).INPEEY0266-5611http://dx.doi.org/10.1088/0266-5611/29/7/075006Google Scholar

13. 

A. Pulkkinen et al., “A Bayesian approach to spectral quantitative photoacoustic tomography,” Inverse Probl. 30, 065012 (2014).INPEEY0266-5611http://dx.doi.org/10.1088/0266-5611/30/6/065012Google Scholar

14. 

T. Ding, K. Ren and S. Vallélian, “A one-step reconstruction algorithm for quantitative photoacoustic imaging,” Inverse Probl. 31, 095005 (2015).INPEEY0266-5611http://dx.doi.org/10.1088/0266-5611/31/9/095005Google Scholar

15. 

L. Yao, Y. Sun and H. Jiang, “Quantitative photoacoustic tomography based on the radiative transfer equation,” Opt. Lett. 34, 1765–1767 (2009).OPLEDP0146-9592http://dx.doi.org/10.1364/OL.34.001765Google Scholar

16. 

S. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999).INPEEY0266-5611http://dx.doi.org/10.1088/0266-5611/15/2/022Google Scholar

17. 

T. Tarvainen et al., “Hybrid radiative-transfer–diffusion model for optical tomography,” Appl. Opt. 44(6), 876–886 (2005).APOPAI0003-6935http://dx.doi.org/10.1364/AO.44.000876Google Scholar

18. 

T. Tarvainen et al., “Coupled radiative transfer equation and diffusion approximation model for photon migration in turbid medium with low-scattering and nonscattering regions,” Phys. Med. Biol. 50, 4913–4930 (2005).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/50/20/011Google Scholar

19. 

T. Tarvainen et al., “Finite element model for the coupled radiative transfer equation and diffusion approximation,” Int. J. Numer. Method Eng. 65(3), 383–405 (2006).http://dx.doi.org/10.1002/(ISSN)1097-0207Google Scholar

20. 

T. Tarvainen, “Computational methods for light transport in optical tomography,” PhD Thesis, University of Kuopio (2006).Google Scholar

21. 

P. Surya Mohan et al., “Variable order spherical harmonic expansion scheme for the radiative transport equation using finite elements,” J. Comput. Phys. 230, 7364–7383 (2011).JCTPAH0021-9991http://dx.doi.org/10.1016/j.jcp.2011.06.004Google Scholar

22. 

L. Wang, S. Jacques and L. Zheng, “MCML–Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47, 131–146 (1995).CMPBEK0169-2607http://dx.doi.org/10.1016/0169-2607(95)01640-FGoogle Scholar

23. 

D. Boas et al., “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10, 159–170 (2002).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.10.000159Google Scholar

24. 

Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17, 20178–20190 (2009).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.17.020178Google Scholar

25. 

S. Powell and T. S. Leung, “Highly parallel Monte-Carlo simulations of the acousto-optic effect in heterogeneous turbid media,” J. Biomed. Opt. 17, 045002 (2012).JBOPFO1083-3668http://dx.doi.org/10.1117/1.JBO.17.4.045002Google Scholar

26. 

P. C. Beard and T. N. Mills, “Characterization of post mortem arterial tissue using time-resolved photoacoustic spectroscopy at 436, 461 and 532 nm,” Phys. Med. Biol. 42, 177–198 (1997).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/42/1/012Google Scholar

27. 

M.-R. Antonelli et al., “Impact of model parameters on Monte Carlo simulations of backscattering Mueller matrix images of colon tissue,” Biomed. Opt. Express 2, 1836–1851 (2011).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.2.001836Google Scholar

28. 

T. S. Leung et al., “Light propagation in a turbid medium with insonified microbubbles,” J. Biomed. Opt. 18, 015002 (2013).JBOPFO1083-3668http://dx.doi.org/10.1117/1.JBO.18.1.015002Google Scholar

29. 

T. Grosges and D. Barchiesi, “Nanoshells for photothermal therapy: a Monte-Carlo based numerical study of their design tolerance,” Biomed. Opt. 2(6), 243–247 (2011).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.2.001584Google Scholar

30. 

N. Manuchehrabadi et al., “Computational simulation of temperature elevations in tumors using Monte Carlo method and comparison to experimental measurements in laser photothermal therapy,” J. Biomech. Eng. 135, 121007 (2013).JBENDY0148-0731http://dx.doi.org/10.1115/1.4025388Google Scholar

31. 

J. Cassidy, V. Betz and L. Lilge, “Treatment plan evaluation for interstitial photodynamic therapy in a mouse model by Monte Carlo simulation with FullMonte,” Front. Phys. 3, 1–10 (2015).http://dx.doi.org/10.3389/fphy.2015.00006Google Scholar

32. 

C. Zhu and Q. Liu, “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt. 18, 050902 (2013).JBOPFO1083-3668http://dx.doi.org/10.1117/1.JBO.18.5.050902Google Scholar

33. 

V. Periyasamy and M. Pramanik, “Monte Carlo simulation of light transport in tissue for optimizing light delivery in photoacoustic imaging of the sentinel lymph node,” J. Biomed. Opt. 18, 106008 (2013).JBOPFO1083-3668http://dx.doi.org/10.1117/1.JBO.18.10.106008Google Scholar

34. 

V. Periyasamy and M. Pramanik, “Monte Carlo simulation of light transport in turbid medium with embedded object–spherical, cylindrical, ellipsoidal, or cuboidal objects embedded within multilayered tissues,” J. Biomed. Opt. 19, 045003 (2014).JBOPFO1083-3668http://dx.doi.org/10.1117/1.JBO.19.4.045003Google Scholar

35. 

Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express 1, 165–175 (2010).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.1.000165Google Scholar

36. 

A. Sassaroli and F. Martelli, “Equivalence of four Monte Carlo methods for photon migration in turbid media,” J. Opt. Soc. Am. A Opt. Image Sci. Vision 29, 2110–2117 (2012).http://dx.doi.org/10.1364/JOSAA.29.002110Google Scholar

37. 

J. Heino et al., “Anisotropic effects in highly scattering media,” Phys. Rev. E 68, 1–8 (2003).http://dx.doi.org/10.1103/PhysRevE.68.031908Google Scholar

38. 

W. Star, J. Marijnissen and M. Gemert, “Light dosimetry in optical phantoms and in tissues: I. Multiple flux and transport theory,” Phys. Med. Biol. 33(4), 437–454 (1988).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/33/4/004Google Scholar

39. 

R. Hochuli et al., “Forward and adjoint radiance Monte Carlo models for quantitative photoacoustic imaging,” Proc. SPIE 9323, 93231P (2015).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.2081407Google Scholar

40. 

R. Hochuli, “Monte Carlo methods in quantitative photoacoustic tomography,” PhD Thesis, University College London (2016).Google Scholar

41. 

J. Bezanson et al., “Julia: a fresh approach to numerical computing,” arXiv 1411.1607, p. 37 (2014).Google Scholar

42. 

M. J. C. Van Gemert et al., “Tissue optics for a slab geometry in the diffusion approximation,” Lasers Med. Sci. 2, 295–302 (1987).http://dx.doi.org/10.1007/BF02594174Google Scholar

43. 

A. K. Jha et al., “Three-dimensional Neumann-series approach to model light transport in nonuniform media,” J. Opt. Soc. Am. A Opt. Image Sci. Vision 29, 1885–1899 (2012).http://dx.doi.org/10.1364/JOSAA.29.001885Google Scholar

44. 

A. Liemert and A. Kienle, “Analytical approach for solving the radiative transfer equation in two-dimensional layered media,” J. Quant. Spectrosc. Radiat. Transfer 113, 559–564 (2012).JQSRAE0022-4073http://dx.doi.org/10.1016/j.jqsrt.2012.01.013Google Scholar

45. 

A. Liemert and A. Kienle, “Radiative transfer in two-dimensional infinitely extended scattering media,” J. Phys. A 44, 505206 (2011).http://dx.doi.org/10.1088/1751-8113/44/50/505206Google Scholar

46. 

A. Liemert, Personal Communication (2015).Google Scholar

47. 

L. Henyey and J. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. 93, 70–83 (1941).http://dx.doi.org/10.1086/144246Google Scholar

48. 

W. Hager and H. Zhang, “A new conjugate gradient method with guaranteed descent and an efficient line search,” SIAM J. Optim. 16(1), 170–192 (2005).SJOPE81095-7189http://dx.doi.org/10.1137/030601880Google Scholar

49. 

J. Nocedal and S. Wright, Numerical Optimization, 2nd ed., Springer, New York, New York (1999).Google Scholar

Biographies for the authors are not available.

Roman Hochuli, Samuel Powell, Simon Arridge, Ben Cox, "Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance," Journal of Biomedical Optics 21(12), 126004 (2 December 2016). http://dx.doi.org/10.1117/1.JBO.21.12.126004
Submission: Received ; Accepted
JOURNAL ARTICLE
9 PAGES


SHARE
KEYWORDS
Scattering

Monte Carlo methods

Absorption

Photoacoustic tomography

3D modeling

Acquisition tracking and pointing

Inverse optics

Back to Top