Research Papers: Imaging

Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance

[+] Author Affiliations
Roman Hochuli, Samuel Powell, Ben Cox

University College London, Department of Medical Physics and Biomedical Engineering, Malet Place, WC1E 6BT London, United Kingdom

Simon Arridge

University College London, Department of Computer Science, Malet Place, WC1E 6BT London, United Kingdom

J. Biomed. Opt. 21(12), 126004 (Dec 02, 2016). doi:10.1117/1.JBO.21.12.126004
History: Received August 25, 2016; Accepted November 7, 2016
Text Size: A A A

Open Access Open Access

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.

Figures in this Article

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 analytically25 or numerically,6,7 often within a minimization framework.815 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 developed1721 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,2325 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 measurements2628 and in dosimetric studies for a range of light-based therapies.2931 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.

The inverse problem in QPAT can be stated as the minimization Display Formula

argminμa(x),μs(x)ε[μa(x),μs(x)],(1)
where the error functional is given by Display Formula
ε=12Ω[Hmeas(x)H(x;μa,μs)]2dx.(2)
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 φ*:Display Formula
εμa=Φ(HmeasH)+Sn1φ*(s^)φ(s^)ds^,(3)
and Display Formula
εμs=Sn1φ*(s^)φ(s^)ds^Sn1Sn1φ*(s^)P(s^,s^)φ(s^)ds^ds^.(4)
In the following two sections, the calculation of the forward radiance φ(s^) and adjoint radiance φ*(s^) using MC models is shown.

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 Display Formula

Lφ=q[s^·+μa(x)+μs(x)]φ(x,s^)μs(x)Sn1P(s^,s^)φ(x,s^)ds^=q(x,s^),(5)
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,38Display Formula
P(s^,s^)=12π1g2[1+g22g(s^·s^)].(6)
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 Display Formula
θ=2arctan{1g1+gtan[πU([0,1])]}.(7)
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 literature2224,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.

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 is39Display Formula

φ(x,θ)=12πa0(x)+1πn=1N=an(x)cos(nθ)+1πn=1N=bn(x)sin(nθ),(8)
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 Display Formula
a0=np=1NpdWnpS1δ(θθnp)dθ=np=1NpdWnp,(9)
Display Formula
an=np=1NpdWnpS1δ(θθnp)cos(nθ)dθ=np=1NpdWnpcos(nθnp),(10)
Display Formula
bn=np=1NpdWnpS1δ(θθnp)sin(nθ)dθ=np=1NpdWnpsin(nθnp),(11)
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

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.4446 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.

Graphic Jump Location
Fig 1
F1 :

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.

The adjoint equation to the RTE is given by Display Formula

L*φ*=q*[s^·+μa(x)+μs(x)]φ*(x,s^)μs(x)Sn1P(s^,s^)φ*(x,s^)ds^=q*(x,s^),(12)
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.

Validation of the Adjoint Model

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

Gu,v=u,G*v,(13)
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., Display Formula
φ(r,s^)=GuLφ=u,φ*(r,s^)=G*vL*φ*=v.(14)
Let u=ρs(r)Θs(s^), v=ρd(r)Θd(s^) and consider three different cases Display Formula
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 Display Formula

Φ1(rd)=Φ1*(rs),(15)
Display Formula
2πφ2(rd,s^)Θd(s^)ds^=Φ2*(rs),(16)
Display Formula
2πφ3(rd,s^)Θd(s^)dr=2πΩφ3*(r,s^)ρs(r)Θs(s^)ds^,(17)
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.

Graphic Jump Location
Fig. 2
F2 :

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.

Graphic Jump Location
Fig. 3
F3 :

(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.

Graphic Jump Location
Fig. 4
F4 :

(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).

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 Display Formula

εμ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θ.(18)
By orthogonality, all terms for which nm integrate to zero and Eq. (18) reduces to Display Formula
εμa=a0(Hmeasμaa0)+2π[14π2a0a0*+1π2n=1anan*cos2(nθ)+1π2n=1bnbn*sin2(nθ)]dθ,(19)
Display Formula
=a0(Hmeasμaa0)+12πa0a0*+1πn=1anan*+1πn=1bnbn*.(20)
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 Display Formula

Sn1Sn1φ*(s^)P(s^,s^)φ(s^)ds^ds^,(21)
which contains the phase function given in Eq. (22) and can be expanded using a Fourier series in powers of g37Display Formula
P(s^·s^;g)=12π+1πl=1glcos(lΔθ),(22)
where Δθ=arccos(s^·s^). Thus we can write Display Formula
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θ,(23)
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 Display Formula
S1S1φ(s^)Pθ(s^,s^)φ*(s^)ds^ds^=12πa0a0*+1πn=1anan*gn+1πn=1bnbn*gn.(24)
Substituting this expression into Eq. (4), we can write the full expression for the functional gradient with respect to the scattering coefficient Display Formula
εμs=12πa0a0*+1πn=1anan*+1πn=1bnbn*12πa0a0*+1πn=1anan*gn+1πn=1bnbn*gn(25)
Display Formula
=1πn=1[anan*+bnbn*](1gn).(26)
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.

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 Display Formula

|ε(i)ε(i1)|/|ε(i)|<109,(27)
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.

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).

Graphic Jump Location
Fig. 5
F5 :

(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.

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.

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.

Graphic Jump Location
Fig. 6
F6 :

(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.

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 Display Formula

ε([μa(i),μs(i)]+α(i)p(i))ε([μa(i),μs(i)])+να(i)ε(i)Tp(i),(28)
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.

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.

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

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.

Cox  B.  et al., “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt.. 17, (6 ), 061202  (2012). 1083-3668 CrossRef
Bal  G., , Jollivet  A., and Jugnon  V., “Inverse transport theory of photoacoustics,” Inverse Probl.. 26, , 025011  (2010). 0266-5611 CrossRef
Bal  G., and Uhlmann  G., “Inverse diffusion theory of photoacoustics,” Inverse Probl.. 26, , 085010  (2010). 0266-5611 CrossRef
Ren  K., , Gao  H., and Zhao  H., “A hybrid reconstruction method for quantitative PAT,” SIAM J. Imaging Sci.. 6, , 32 –55 (2013).CrossRef
Ammari  H.  et al., “Mathematical modeling in photoacoustic imaging of small absorbers,” SIAM Rev.. 52, (4 ), 677 –695 (2010). 0036-1445 CrossRef
Zemp  R. J., “Quantitative photoacoustic tomography with multiple optical sources,” Appl. Opt.. 49, , 3566 –3572 (2010). 0003-6935 CrossRef
Harrison  T., , Shao  P., and Zemp  R. J., “A least-squares fixed-point iterative algorithm for multiple illumination photoacoustic tomography,” Biomed. Opt. Express. 4, , 2224 –2230 (2013). 2156-7085 CrossRef
Laufer  J.  et al., “Quantitative determination of chromophore concentrations from 2D photoacoustic images using a nonlinear model-based inversion scheme,” Appl. Opt.. 49, , 1219 –1233 (2010). 0003-6935 CrossRef
Cox  B. T., , Arridge  S. R., and Beard  P. C., “Estimating chromophore distributions from multiwavelength photoacoustic images,” J. Opt. Soc. Am. A Opt. Image Sci. Vision. 26, , 443 –455 (2009).CrossRef
Tarvainen  T.  et al., “Bayesian image reconstruction in quantitative photoacoustic tomography,” IEEE Trans. Med. Imaging. 32, , 2287 –2298 (2013). 0278-0062 CrossRef
Malone  E. R.  et al., “Reconstruction-classification method for quantitative photoacoustic tomography,” arXiv (2015).
Saratoon  T.  et al., “A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation,” Inverse Probl.. 29, , 075006  (2013). 0266-5611 CrossRef
Pulkkinen  A.  et al., “A Bayesian approach to spectral quantitative photoacoustic tomography,” Inverse Probl.. 30, , 065012  (2014). 0266-5611 CrossRef
Ding  T., , Ren  K., and Vallélian  S., “A one-step reconstruction algorithm for quantitative photoacoustic imaging,” Inverse Probl.. 31, , 095005  (2015). 0266-5611 CrossRef
Yao  L., , Sun  Y., and Jiang  H., “Quantitative photoacoustic tomography based on the radiative transfer equation,” Opt. Lett.. 34, , 1765 –1767 (2009). 0146-9592 CrossRef
Arridge  S., “Optical tomography in medical imaging,” Inverse Probl.. 15, , R41 –R93 (1999). 0266-5611 CrossRef
Tarvainen  T.  et al., “Hybrid radiative-transfer–diffusion model for optical tomography,” Appl. Opt.. 44, (6 ), 876 –886 (2005). 0003-6935 CrossRef
Tarvainen  T.  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). 0031-9155 CrossRef
Tarvainen  T.  et al., “Finite element model for the coupled radiative transfer equation and diffusion approximation,” Int. J. Numer. Method Eng.. 65, (3 ), 383 –405 (2006).CrossRef
Tarvainen  T., “Computational methods for light transport in optical tomography,” PhD Thesis, University of Kuopio (2006).
Surya Mohan  P.  et al., “Variable order spherical harmonic expansion scheme for the radiative transport equation using finite elements,” J. Comput. Phys.. 230, , 7364 –7383 (2011). 0021-9991 CrossRef
Wang  L., , Jacques  S., and Zheng  L., “MCML–Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed.. 47, , 131 –146 (1995). 0169-2607 CrossRef
Boas  D.  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). 1094-4087 CrossRef
Fang  Q., and Boas  D. A., “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express. 17, , 20178 –20190 (2009). 1094-4087 CrossRef
Powell  S., and Leung  T. S., “Highly parallel Monte-Carlo simulations of the acousto-optic effect in heterogeneous turbid media,” J. Biomed. Opt.. 17, , 045002  (2012). 1083-3668 CrossRef
Beard  P. C., and Mills  T. N., “Characterization of post mortem arterial tissue using time-resolved photoacoustic spectroscopy at 436, 461 and 532 nm,” Phys. Med. Biol.. 42, , 177 –198 (1997). 0031-9155 CrossRef
Antonelli  M.-R.  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). 2156-7085 CrossRef
Leung  T. S.  et al., “Light propagation in a turbid medium with insonified microbubbles,” J. Biomed. Opt.. 18, , 015002  (2013). 1083-3668 CrossRef
Grosges  T., and Barchiesi  D., “Nanoshells for photothermal therapy: a Monte-Carlo based numerical study of their design tolerance,” Biomed. Opt.. 2, (6 ), 243 –247 (2011). 2156-7085 CrossRef
Manuchehrabadi  N.  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). 0148-0731 CrossRef
Cassidy  J., , Betz  V., and Lilge  L., “Treatment plan evaluation for interstitial photodynamic therapy in a mouse model by Monte Carlo simulation with FullMonte,” Front. Phys.. 3, , 1 –10 (2015).CrossRef
Zhu  C., and Liu  Q., “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt.. 18, , 050902  (2013). 1083-3668 CrossRef
Periyasamy  V., and Pramanik  M., “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). 1083-3668 CrossRef
Periyasamy  V., and Pramanik  M., “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). 1083-3668 CrossRef
Fang  Q., “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express. 1, , 165 –175 (2010). 2156-7085 CrossRef
Sassaroli  A., and Martelli  F., “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).CrossRef
Heino  J.  et al., “Anisotropic effects in highly scattering media,” Phys. Rev. E. 68, , 1 –8 (2003).CrossRef
Star  W., , Marijnissen  J., and Gemert  M., “Light dosimetry in optical phantoms and in tissues: I. Multiple flux and transport theory,” Phys. Med. Biol.. 33, (4 ), 437 –454 (1988). 0031-9155 CrossRef
Hochuli  R.  et al., “Forward and adjoint radiance Monte Carlo models for quantitative photoacoustic imaging,” Proc. SPIE. 9323, , 93231P  (2015). 0277-786X CrossRef
Hochuli  R., “Monte Carlo methods in quantitative photoacoustic tomography,” PhD Thesis, University College London (2016).
Bezanson  J.  et al., “Julia: a fresh approach to numerical computing,” arXiv 1411.1607, p. 37  (2014).
Van Gemert  M. J. C.  et al., “Tissue optics for a slab geometry in the diffusion approximation,” Lasers Med. Sci.. 2, , 295 –302 (1987).CrossRef
Jha  A. K.  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).CrossRef
Liemert  A., and Kienle  A., “Analytical approach for solving the radiative transfer equation in two-dimensional layered media,” J. Quant. Spectrosc. Radiat. Transfer. 113, , 559 –564 (2012). 0022-4073 CrossRef
Liemert  A., and Kienle  A., “Radiative transfer in two-dimensional infinitely extended scattering media,” J. Phys. A. 44, , 505206  (2011).CrossRef
Liemert  A., Personal Communication. (2015).
Henyey  L., and Greenstein  J., “Diffuse radiation in the galaxy,” Astrophys. J.. 93, , 70 –83 (1941).CrossRef
Hager  W., and Zhang  H., “A new conjugate gradient method with guaranteed descent and an efficient line search,” SIAM J. Optim.. 16, (1 ), 170 –192 (2005). 1095-7189 CrossRef
Nocedal  J., and Wright  S., Numerical Optimization. , 2nd ed.,  Springer ,  New York, New York  (1999 ).
© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation

Roman Hochuli ; Samuel Powell ; Simon Arridge and Ben Cox
"Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance", J. Biomed. Opt. 21(12), 126004 (Dec 02, 2016). ; http://dx.doi.org/10.1117/1.JBO.21.12.126004


Figures

Graphic Jump Location
Fig. 3
F3 :

(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.

Graphic Jump Location
Fig. 2
F2 :

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.

Graphic Jump Location
Fig 1
F1 :

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.

Graphic Jump Location
Fig. 4
F4 :

(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).

Graphic Jump Location
Fig. 5
F5 :

(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.

Graphic Jump Location
Fig. 6
F6 :

(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.

Tables

References

Cox  B.  et al., “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt.. 17, (6 ), 061202  (2012). 1083-3668 CrossRef
Bal  G., , Jollivet  A., and Jugnon  V., “Inverse transport theory of photoacoustics,” Inverse Probl.. 26, , 025011  (2010). 0266-5611 CrossRef
Bal  G., and Uhlmann  G., “Inverse diffusion theory of photoacoustics,” Inverse Probl.. 26, , 085010  (2010). 0266-5611 CrossRef
Ren  K., , Gao  H., and Zhao  H., “A hybrid reconstruction method for quantitative PAT,” SIAM J. Imaging Sci.. 6, , 32 –55 (2013).CrossRef
Ammari  H.  et al., “Mathematical modeling in photoacoustic imaging of small absorbers,” SIAM Rev.. 52, (4 ), 677 –695 (2010). 0036-1445 CrossRef
Zemp  R. J., “Quantitative photoacoustic tomography with multiple optical sources,” Appl. Opt.. 49, , 3566 –3572 (2010). 0003-6935 CrossRef
Harrison  T., , Shao  P., and Zemp  R. J., “A least-squares fixed-point iterative algorithm for multiple illumination photoacoustic tomography,” Biomed. Opt. Express. 4, , 2224 –2230 (2013). 2156-7085 CrossRef
Laufer  J.  et al., “Quantitative determination of chromophore concentrations from 2D photoacoustic images using a nonlinear model-based inversion scheme,” Appl. Opt.. 49, , 1219 –1233 (2010). 0003-6935 CrossRef
Cox  B. T., , Arridge  S. R., and Beard  P. C., “Estimating chromophore distributions from multiwavelength photoacoustic images,” J. Opt. Soc. Am. A Opt. Image Sci. Vision. 26, , 443 –455 (2009).CrossRef
Tarvainen  T.  et al., “Bayesian image reconstruction in quantitative photoacoustic tomography,” IEEE Trans. Med. Imaging. 32, , 2287 –2298 (2013). 0278-0062 CrossRef
Malone  E. R.  et al., “Reconstruction-classification method for quantitative photoacoustic tomography,” arXiv (2015).
Saratoon  T.  et al., “A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation,” Inverse Probl.. 29, , 075006  (2013). 0266-5611 CrossRef
Pulkkinen  A.  et al., “A Bayesian approach to spectral quantitative photoacoustic tomography,” Inverse Probl.. 30, , 065012  (2014). 0266-5611 CrossRef
Ding  T., , Ren  K., and Vallélian  S., “A one-step reconstruction algorithm for quantitative photoacoustic imaging,” Inverse Probl.. 31, , 095005  (2015). 0266-5611 CrossRef
Yao  L., , Sun  Y., and Jiang  H., “Quantitative photoacoustic tomography based on the radiative transfer equation,” Opt. Lett.. 34, , 1765 –1767 (2009). 0146-9592 CrossRef
Arridge  S., “Optical tomography in medical imaging,” Inverse Probl.. 15, , R41 –R93 (1999). 0266-5611 CrossRef
Tarvainen  T.  et al., “Hybrid radiative-transfer–diffusion model for optical tomography,” Appl. Opt.. 44, (6 ), 876 –886 (2005). 0003-6935 CrossRef
Tarvainen  T.  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). 0031-9155 CrossRef
Tarvainen  T.  et al., “Finite element model for the coupled radiative transfer equation and diffusion approximation,” Int. J. Numer. Method Eng.. 65, (3 ), 383 –405 (2006).CrossRef
Tarvainen  T., “Computational methods for light transport in optical tomography,” PhD Thesis, University of Kuopio (2006).
Surya Mohan  P.  et al., “Variable order spherical harmonic expansion scheme for the radiative transport equation using finite elements,” J. Comput. Phys.. 230, , 7364 –7383 (2011). 0021-9991 CrossRef
Wang  L., , Jacques  S., and Zheng  L., “MCML–Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed.. 47, , 131 –146 (1995). 0169-2607 CrossRef
Boas  D.  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). 1094-4087 CrossRef
Fang  Q., and Boas  D. A., “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express. 17, , 20178 –20190 (2009). 1094-4087 CrossRef
Powell  S., and Leung  T. S., “Highly parallel Monte-Carlo simulations of the acousto-optic effect in heterogeneous turbid media,” J. Biomed. Opt.. 17, , 045002  (2012). 1083-3668 CrossRef
Beard  P. C., and Mills  T. N., “Characterization of post mortem arterial tissue using time-resolved photoacoustic spectroscopy at 436, 461 and 532 nm,” Phys. Med. Biol.. 42, , 177 –198 (1997). 0031-9155 CrossRef
Antonelli  M.-R.  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). 2156-7085 CrossRef
Leung  T. S.  et al., “Light propagation in a turbid medium with insonified microbubbles,” J. Biomed. Opt.. 18, , 015002  (2013). 1083-3668 CrossRef
Grosges  T., and Barchiesi  D., “Nanoshells for photothermal therapy: a Monte-Carlo based numerical study of their design tolerance,” Biomed. Opt.. 2, (6 ), 243 –247 (2011). 2156-7085 CrossRef
Manuchehrabadi  N.  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). 0148-0731 CrossRef
Cassidy  J., , Betz  V., and Lilge  L., “Treatment plan evaluation for interstitial photodynamic therapy in a mouse model by Monte Carlo simulation with FullMonte,” Front. Phys.. 3, , 1 –10 (2015).CrossRef
Zhu  C., and Liu  Q., “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt.. 18, , 050902  (2013). 1083-3668 CrossRef
Periyasamy  V., and Pramanik  M., “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). 1083-3668 CrossRef
Periyasamy  V., and Pramanik  M., “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). 1083-3668 CrossRef
Fang  Q., “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express. 1, , 165 –175 (2010). 2156-7085 CrossRef
Sassaroli  A., and Martelli  F., “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).CrossRef
Heino  J.  et al., “Anisotropic effects in highly scattering media,” Phys. Rev. E. 68, , 1 –8 (2003).CrossRef
Star  W., , Marijnissen  J., and Gemert  M., “Light dosimetry in optical phantoms and in tissues: I. Multiple flux and transport theory,” Phys. Med. Biol.. 33, (4 ), 437 –454 (1988). 0031-9155 CrossRef
Hochuli  R.  et al., “Forward and adjoint radiance Monte Carlo models for quantitative photoacoustic imaging,” Proc. SPIE. 9323, , 93231P  (2015). 0277-786X CrossRef
Hochuli  R., “Monte Carlo methods in quantitative photoacoustic tomography,” PhD Thesis, University College London (2016).
Bezanson  J.  et al., “Julia: a fresh approach to numerical computing,” arXiv 1411.1607, p. 37  (2014).
Van Gemert  M. J. C.  et al., “Tissue optics for a slab geometry in the diffusion approximation,” Lasers Med. Sci.. 2, , 295 –302 (1987).CrossRef
Jha  A. K.  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).CrossRef
Liemert  A., and Kienle  A., “Analytical approach for solving the radiative transfer equation in two-dimensional layered media,” J. Quant. Spectrosc. Radiat. Transfer. 113, , 559 –564 (2012). 0022-4073 CrossRef
Liemert  A., and Kienle  A., “Radiative transfer in two-dimensional infinitely extended scattering media,” J. Phys. A. 44, , 505206  (2011).CrossRef
Liemert  A., Personal Communication. (2015).
Henyey  L., and Greenstein  J., “Diffuse radiation in the galaxy,” Astrophys. J.. 93, , 70 –83 (1941).CrossRef
Hager  W., and Zhang  H., “A new conjugate gradient method with guaranteed descent and an efficient line search,” SIAM J. Optim.. 16, (1 ), 170 –192 (2005). 1095-7189 CrossRef
Nocedal  J., and Wright  S., Numerical Optimization. , 2nd ed.,  Springer ,  New York, New York  (1999 ).

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging & repositioning the boxes below.

Related Book Chapters

Topic Collections

PubMed Articles
Advertisement
  • Don't have an account?
  • Subscribe to the SPIE Digital Library
  • Create a FREE account to sign up for Digital Library content alerts and gain access to institutional subscriptions remotely.
Access This Article
Sign in or Create a personal account to Buy this article ($20 for members, $25 for non-members).
Access This Proceeding
Sign in or Create a personal account to Buy this article ($15 for members, $18 for non-members).
Access This Chapter

Access to SPIE eBooks is limited to subscribing institutions and is not available as part of a personal subscription. Print or electronic versions of individual SPIE books may be purchased via SPIE.org.