## 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 analytically^{2}3.4.^{–}^{5} or numerically,^{6}^{,}^{7} often within a minimization framework.^{8}9.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 developed^{17}18.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}^{,}^{23}24.^{–}^{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 measurements^{26}27.^{–}^{28} and in dosimetric studies for a range of light-based therapies.^{29}30.^{–}^{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 cylindrical^{34} 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)

$$\underset{{\mu}_{a}(\mathbf{x}),{\mu}_{s}(\mathbf{x})}{\mathrm{argmin}}\text{\hspace{0.17em}}\epsilon [{\mu}_{a}(\mathbf{x}),{\mu}_{s}(\mathbf{x})],$$## (2)

$$\epsilon =\frac{1}{2}{\int}_{\mathrm{\Omega}}{[{H}^{\text{meas}}(\mathbf{x})-H(\mathbf{x};{\mu}_{\mathrm{a}},{\mu}_{\mathrm{s}})]}^{2}\mathrm{d}\mathbf{x}.$$^{12}give expressions for these gradients in terms of the forward and adjoint fields, $\phi $ and ${\phi}^{*}$:

## (3)

$$\frac{\partial \epsilon}{\partial {\mu}_{\mathrm{a}}}=-\mathrm{\Phi}({H}^{\text{meas}}-H)+{\int}_{{\mathcal{S}}^{n-1}}{\phi}^{*}(\widehat{\mathbf{s}})\phi (\widehat{\mathbf{s}})\mathrm{d}\widehat{\mathbf{s}},$$## (4)

$$\frac{\partial \epsilon}{\partial {\mu}_{\mathrm{s}}}={\int}_{{\mathcal{S}}^{n-1}}{\phi}^{*}(\widehat{\mathbf{s}})\phi (\widehat{\mathbf{s}})\mathrm{d}\widehat{\mathbf{s}}-{\int}_{{\mathcal{S}}^{n-1}}{\int}_{{\mathcal{S}}^{n-1}}{\phi}^{*}(\widehat{\mathbf{s}})P(\widehat{\mathbf{s}},{\widehat{\mathbf{s}}}^{\prime})\phi ({\widehat{\mathbf{s}}}^{\prime})\mathrm{d}{\widehat{\mathbf{s}}}^{\prime}\text{\hspace{0.17em}}\mathrm{d}\widehat{\mathbf{s}}.$$## 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)

$$\mathcal{L}\phi =q\phantom{\rule{0ex}{0ex}}\equiv [\widehat{\mathbf{s}}\xb7\nabla +{\mu}_{\mathrm{a}}(\mathbf{x})+{\mu}_{\mathrm{s}}(\mathbf{x})]\phi (\mathbf{x},\widehat{\mathbf{s}})-{\mu}_{\mathrm{s}}(\mathbf{x}){\int}_{{\mathcal{S}}^{n-1}}P(\widehat{\mathbf{s}},{\widehat{\mathbf{s}}}^{\prime})\phi (\mathbf{x},{\widehat{\mathbf{s}}}^{\prime})\mathrm{d}{\widehat{\mathbf{s}}}^{\prime}=q(\mathbf{x},\widehat{\mathbf{s}}),$$^{36}The approach used here begins with launching a packet of energy, referred to herein as a “photon,” from a given position $\mathbf{x}$ in an initial direction $\widehat{\mathbf{s}}$. After traveling a distance $s=\mathcal{U}([0,1])/{\mu}_{\mathrm{s}}$ (using the convention $\mathbf{s}=|\mathbf{s}|\widehat{\mathbf{s}}$), where $\mathcal{U}([0,1])$ is a real uniform random variable on [0, 1], a fraction of the photon’s “weight” $W[1-\mathrm{exp}(-{\mu}_{\mathrm{a}}s)]$ is deposited in the current voxel, where $W$ is the current weight (or energy) of the photon packet. The photon weight is updated: $W\leftarrow W\text{\hspace{0.17em}}\mathrm{exp}(-{\mu}_{\mathrm{a}}s)$. Scattering into a new direction ${\widehat{\mathbf{s}}}^{\prime}$ in two-dimensional (2-D) involves sampling the scattering phase function, which describes the probability of a photon scattering from direction ${\widehat{\mathbf{s}}}^{\prime}$ into direction $\widehat{\mathbf{s}}$. The phase function used here was the 2-D Henyey–Greenstein phase function, commonly used in biomedical optics,

^{37}

^{,}

^{38}

## (6)

$$P(\widehat{\mathbf{s}},{\widehat{\mathbf{s}}}^{\prime})=\frac{1}{2\pi}\frac{1-{g}^{2}}{[1+{g}^{2}-2g(\widehat{\mathbf{s}}\xb7{\widehat{\mathbf{s}}}^{\prime})]}.$$## (7)

$$\theta =2\text{\hspace{0.17em}}\mathrm{arctan}\{\frac{1-g}{1+g}\text{\hspace{0.17em}}\mathrm{tan}[\pi \mathcal{U}([0,1])]\}.$$By calculating photon paths through the medium, the MC models presented in the literature^{22}23.^{–}^{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 is^{39}

## (8)

$$\phi (\mathbf{x},\theta )=\frac{1}{2\pi}{a}_{0}(\mathbf{x})+\frac{1}{\pi}\sum _{n=1}^{N=\infty}{a}_{n}(\mathbf{x})\mathrm{cos}(n\theta )+\frac{1}{\pi}\sum _{n=1}^{N=\infty}{b}_{n}(\mathbf{x})\mathrm{sin}(n\theta ),$$^{40}) For a given voxel, the weight is deposited into the relevant Fourier coefficients according to

## (9)

$${a}_{0}=\sum _{{n}_{p}=1}^{{N}_{p}}d{W}_{{n}_{p}}{\int}_{{\mathcal{S}}^{1}}\delta ({\theta}^{\prime}-{\theta}_{{n}_{p}})\mathrm{d}{\theta}^{\prime}=\sum _{{n}_{p}=1}^{{N}_{p}}d{W}_{{n}_{p}},$$## (10)

$${a}_{n}=\sum _{{n}_{p}=1}^{{N}_{p}}d{W}_{{n}_{p}}{\int}_{{\mathcal{S}}^{1}}\delta ({\theta}^{\prime}-{\theta}_{{n}_{p}})\mathrm{cos}(n\theta )\mathrm{d}{\theta}^{\prime}=\sum _{{n}_{p}=1}^{{N}_{p}}d{W}_{{n}_{p}}\text{\hspace{0.17em}}\mathrm{cos}(n{\theta}_{{n}_{p}}),$$## (11)

$${b}_{n}=\sum _{{n}_{p}=1}^{{N}_{p}}d{W}_{{n}_{p}}{\int}_{{\mathcal{S}}^{1}}\delta ({\theta}^{\prime}-{\theta}_{{n}_{p}})\mathrm{sin}(n\theta )\mathrm{d}{\theta}^{\prime}=\sum _{{n}_{p}=1}^{{N}_{p}}d{W}_{{n}_{p}}\text{\hspace{0.17em}}\mathrm{sin}(n{\theta}_{{n}_{p}}),$$^{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.^{44}45.^{–}^{46} An isotropic point source was placed at the center of a domain of size $15\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 15\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$, respectively, and the Henyey–Greenstein phase function^{47} was used with $g$ set to 0.9. The pixel size was $0.05\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 0.05\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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.

## 4.

## Adjoint Monte Carlo Model

The adjoint equation to the RTE is given by

## (12)

$${\mathcal{L}}^{*}{\phi}^{*}={q}^{*}\equiv [-\widehat{\mathbf{s}}\xb7\nabla +{\mu}_{\mathrm{a}}(\mathbf{x})+{\mu}_{\mathrm{s}}(\mathbf{x})]{\phi}^{*}(\mathbf{x},\widehat{\mathbf{s}})-{\mu}_{\mathrm{s}}(\mathbf{x}){\int}_{{\mathcal{S}}^{n-1}}P({\widehat{\mathbf{s}}}^{\prime},\widehat{\mathbf{s}}){\phi}^{*}(\mathbf{x},{\widehat{\mathbf{s}}}^{\prime})\mathrm{d}{\widehat{\mathbf{s}}}^{\prime}={q}^{*}(\mathbf{x},\widehat{\mathbf{s}}),$$## 4.1.

### Validation of the Adjoint Model

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

## (13)

$$\u27e8\mathcal{G}\mathbf{u},\mathbf{v}\u27e9=\u27e8\mathbf{u},{\mathcal{G}}^{*}\mathbf{v}\u27e9,$$## (14)

$$\phi (\mathbf{r},\widehat{\mathbf{s}})=\mathcal{G}\mathbf{u}\iff \mathcal{L}\phi =\mathbf{u},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{\phi}^{*}(\mathbf{r},\widehat{\mathbf{s}})={\mathcal{G}}^{*}\mathbf{v}\iff {\mathcal{L}}^{*}{\phi}^{*}=\mathbf{v}.$$Substituting these into Eq. (13) yields

## (16)

$${\int}_{2\pi}{\phi}_{2}({\mathbf{r}}_{d},\widehat{\mathbf{s}}){\mathrm{\Theta}}_{\mathrm{d}}(\widehat{\mathbf{s}})\mathrm{d}\widehat{\mathbf{s}}={\mathrm{\Phi}}_{2}^{*}({\mathbf{r}}_{s}),$$## (17)

$${\int}_{2\pi}{\phi}_{3}({\mathbf{r}}_{d},\widehat{\mathbf{s}}){\mathrm{\Theta}}_{\mathrm{d}}(\widehat{\mathbf{s}})\mathrm{d}\mathbf{r}={\int}_{2\pi}{\int}_{\mathrm{\Omega}}{\phi}_{3}^{*}(\mathbf{r},\widehat{\mathbf{s}}){\rho}_{\mathrm{s}}(\mathbf{r}){\mathrm{\Theta}}_{\mathrm{s}}(\widehat{\mathbf{s}})\mathrm{d}\widehat{\mathbf{s}},$$Simulations were performed using a $40\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 40\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ ($101\times 101\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$) domain, and 10 Fourier harmonics. Each source distribution emitted ${10}^{6}$ photons. ${\mathbf{r}}_{s}$ was set to be the center of the domain with ${\mathbf{r}}_{d}$ 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 ${\rho}_{\mathrm{d}}=\delta (\mathbf{r}-{\mathbf{r}}_{s})$, ${\mathrm{\Theta}}_{\mathrm{d}}=\frac{1}{\pi}\text{\hspace{0.17em}}{\mathrm{sin}}^{2}(2\theta )$, and Fig. 4 for case 3 with the same ${\rho}_{\mathrm{d}}$, ${\mathrm{\Theta}}_{\mathrm{d}}$ but with istropic ${\mathrm{\Theta}}_{\mathrm{s}}$ and the distributed ${\mathrm{\Theta}}_{\mathrm{s}}(\mathbf{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 $\u27e8\mathcal{G}{\mathbf{u}}_{1,2,3},{\mathbf{v}}_{1,2,3}\u27e9$ with $\u27e8{\mathbf{u}}_{1,2,3},{\mathcal{G}}^{*}{\mathbf{v}}_{1,2,3}\u27e9$ 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.

## 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 ${a}_{0}$. The other terms in the expressions for the functional gradients contain integrals of products of the radiance and its adjoint. If ${a}_{0}^{*}$, ${a}_{n}^{*}$, and ${b}_{n}^{*}$ are the Fourier coefficients of the adjoint radiance, then the gradient with respect to absorption can be written as

## (18)

$$\frac{\partial \epsilon}{\partial {\mu}_{\mathrm{a}}}=-\mathrm{\Phi}({H}^{\text{meas}}-H)+{\int}_{{\mathcal{S}}^{1}}\phi (\widehat{\mathbf{s}}){\phi}^{*}(\widehat{\mathbf{s}})\mathrm{d}\widehat{\mathbf{s}}\phantom{\rule{0ex}{0ex}}=-{a}_{0}({H}^{\text{meas}}-{\mu}_{\mathrm{a}}{a}_{0})+{\int}_{2\pi}[\frac{1}{4{\pi}^{2}}{a}_{0}{a}_{0}^{*}+\frac{1}{2{\pi}^{2}}{a}_{0}\sum _{m=1}^{\infty}{a}_{m}^{*}\text{\hspace{0.17em}}\mathrm{cos}(m{\theta}^{\prime})+\frac{1}{2{\pi}^{2}}{a}_{0}\sum _{m=1}^{\infty}{a}_{m}^{*}\text{\hspace{0.17em}}\mathrm{sin}(m{\theta}^{\prime})+\frac{1}{2{\pi}^{2}}{a}_{0}^{*}\sum _{n=1}^{\infty}{a}_{n}\text{\hspace{0.17em}}\mathrm{cos}(n\theta )+\sum _{n=1}^{\infty}\sum _{m=1}^{\infty}{a}_{n}{a}_{m}^{*}\text{\hspace{0.17em}}\mathrm{cos}(n\theta )\mathrm{cos}(m\theta )+\sum _{n=1}^{\infty}\sum _{m=1}^{\infty}{a}_{n}{b}_{m}^{*}\text{\hspace{0.17em}}\mathrm{cos}(n\theta )\mathrm{sin}(m\theta )\phantom{\rule{0ex}{0ex}}+\frac{1}{2{\pi}^{2}}{a}_{0}^{*}\sum _{n=1}^{\infty}{b}_{m}\text{\hspace{0.17em}}\mathrm{cos}(m\theta )+\sum _{n=1}^{\infty}\sum _{m=1}^{\infty}{a}_{m}^{*}{b}_{n}\text{\hspace{0.17em}}\mathrm{sin}(n\theta )\mathrm{cos}(m\theta )+\sum _{n}\sum _{m}{b}_{n}{b}_{m}^{*}\text{\hspace{0.17em}}\mathrm{sin}(n\theta )\mathrm{sin}(m\theta )]\mathrm{d}\theta .$$## (19)

$$\frac{\partial \epsilon}{\partial {\mu}_{\mathrm{a}}}=-{a}_{0}({H}^{\text{meas}}-{\mu}_{\mathrm{a}}{a}_{0})+{\int}_{2\pi}[\frac{1}{4{\pi}^{2}}{a}_{0}{a}_{0}^{*}+\frac{1}{{\pi}^{2}}\sum _{n=1}^{\infty}{a}_{n}{a}_{n}^{*}\text{\hspace{0.17em}}{\mathrm{cos}}^{2}(n\theta )+\frac{1}{{\pi}^{2}}\sum _{n=1}^{\infty}{b}_{n}{b}_{n}^{*}\text{\hspace{0.17em}}{\mathrm{sin}}^{2}(n\theta )]\mathrm{d}\theta ,$$## (20)

$$=-{a}_{0}({H}^{\text{meas}}-{\mu}_{\mathrm{a}}{a}_{0})+\frac{1}{2\pi}{a}_{0}{a}_{0}^{*}+\frac{1}{\pi}\sum _{n=1}^{\infty}{a}_{n}{a}_{n}^{*}+\frac{1}{\pi}\sum _{n=1}^{\infty}{b}_{n}{b}_{n}^{*}.$$The second term in Eq. (4) is

## (21)

$${\int}_{{\mathcal{S}}^{n-1}}{\int}_{{\mathcal{S}}^{n-1}}{\phi}^{*}(\widehat{\mathbf{s}})P(\widehat{\mathbf{s}},{\widehat{\mathbf{s}}}^{\prime})\phi ({\widehat{\mathbf{s}}}^{\prime})\mathrm{d}{\widehat{\mathbf{s}}}^{\prime}\text{\hspace{0.17em}}\mathrm{d}\widehat{\mathbf{s}},$$^{37}

## (22)

$$P(\widehat{\mathbf{s}}\xb7{\widehat{\mathbf{s}}}^{\prime};g)=\frac{1}{2\pi}+\frac{1}{\pi}\sum _{l=1}^{\infty}{g}^{l}\text{\hspace{0.17em}}\mathrm{cos}(l\mathrm{\Delta}\theta ),$$## (23)

$${\int}_{2\pi}{\int}_{2\pi}\phi ({\widehat{\mathbf{s}}}^{\prime}){P}_{\theta}(\widehat{\mathbf{s}},{\widehat{\mathbf{s}}}^{\prime}){\phi}^{*}(\widehat{\mathbf{s}})\mathrm{d}\widehat{\mathbf{s}}\text{\hspace{0.17em}}\mathrm{d}{\widehat{\mathbf{s}}}^{\prime}={\int}_{2\pi}{\int}_{2\pi}[\frac{1}{2\pi}{a}_{0}+\frac{1}{\pi}\sum _{n=1}^{\infty}{a}_{n}\text{\hspace{0.17em}}\mathrm{cos}(n{\theta}^{\prime})+\frac{1}{\pi}\sum _{n=1}^{\infty}{b}_{n}\text{\hspace{0.17em}}\mathrm{sin}(n{\theta}^{\prime})]\phantom{\rule{0ex}{0ex}}\{\frac{1}{2\pi}+\frac{1}{\pi}\sum _{l=0}^{\infty}{g}^{l}\text{\hspace{0.17em}}\mathrm{cos}[l(\theta -{\theta}^{\prime})]\}\phantom{\rule{0ex}{0ex}}[\frac{1}{2\pi}{a}_{0}^{*}+\frac{1}{\pi}\sum _{m=1}^{\infty}{a}_{m}^{*}\text{\hspace{0.17em}}\mathrm{cos}(m\theta )+\frac{1}{\pi}\sum _{m=1}^{\infty}{b}_{m}^{*}\text{\hspace{0.17em}}\mathrm{sin}(m\theta )]\mathrm{d}\theta \text{\hspace{0.17em}}\mathrm{d}{\theta}^{\prime},$$## (24)

$${\int}_{{\mathcal{S}}^{1}}{\int}_{{\mathcal{S}}^{1}}\phi ({\widehat{\mathbf{s}}}^{\prime}){P}_{\theta}(\widehat{\mathbf{s}},{\widehat{\mathbf{s}}}^{\prime}){\phi}^{*}(\widehat{\mathbf{s}})\mathrm{d}\widehat{\mathbf{s}}\text{\hspace{0.17em}}\mathrm{d}{\widehat{\mathbf{s}}}^{\prime}=\frac{1}{2\pi}{a}_{0}{a}_{0}^{*}+\frac{1}{\pi}\sum _{n=1}^{\infty}{a}_{n}{a}_{n}^{*}{g}^{n}+\frac{1}{\pi}\sum _{n=1}^{\infty}{b}_{n}{b}_{n}^{*}{g}^{n}.$$## (25)

$$\frac{\partial \epsilon}{\partial {\mu}_{s}}=\frac{1}{2\pi}{a}_{0}{a}_{0}^{*}+\frac{1}{\pi}\sum _{n=1}^{\infty}{a}_{n}{a}_{n}^{*}+\frac{1}{\pi}\sum _{n=1}^{\infty}{b}_{n}{b}_{n}^{*}-\frac{1}{2\pi}{a}_{0}{a}_{0}^{*}+\frac{1}{\pi}\sum _{n=1}^{\infty}{a}_{n}{a}_{n}^{*}{g}^{n}+\frac{1}{\pi}\sum _{n=1}^{\infty}{b}_{n}{b}_{n}^{*}{g}^{n}$$## 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 ${\mu}_{\mathrm{a}}(\mathbf{x})$ and ${\mu}_{\mathrm{s}}(\mathbf{x})$ from simulated PAT images by minimizing the error functional in Eq. (2). As the adjoint source, ${q}^{*}(\mathbf{x},\widehat{\mathbf{s}})={\mu}_{\mathrm{a}}(\mathbf{x})[{H}^{\text{meas}}(\mathbf{x})-H(\mathbf{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 ${N}_{p}$. 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 Zhang^{48} was used for the reconstruction of ${\mu}_{\mathrm{a}}$. A backtracking line search was implemented for the reconstruction of ${\mu}_{\mathrm{s}}$. This is described in Sec. 6.2.

The termination condition used by the optimization was

where $i$ is the iteration number. For all the reconstructions, it was assumed that the data ${H}^{\text{meas}}$ 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 ${10}^{8}$ 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\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$ with a rectangular inclusion equal to $0.2\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$ [shown in Fig. 5(a)], and a background scattering coefficient of $5\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$ with a rectangular inclusion equal to $15\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$ [shown in Fig. 5(b)]. The anisotropy was a homogeneously distributed value of 0.9. The measured data, ${H}^{\text{meas}}$, were formed using an MC simulation illuminated by a collimated line source on the boundary at $z=0\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ and on the adjacent boundary at $x=4\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$, consisting of ${10}^{8}$ 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\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$. 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 ${\mu}_{\mathrm{a}}$ shown in Fig. 5(d).

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\times {10}^{-9}$. The optimization routine was stopped as the change in the error function on the 12th iteration was below the function tolerance of ${10}^{-9}$, indicating convergence. The average error in the estimate of the absorption coefficient ${\mu}_{\mathrm{a}}^{\mathrm{est}}$, computed as $|{\mu}_{\mathrm{a}}^{\text{true}}-{\mu}_{\mathrm{a}}^{\mathrm{est}}|/{\mu}_{\mathrm{a}}^{\text{true}}$, 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\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$.

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 search^{49} 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)

$$\epsilon ([{\mu}_{\mathrm{a}}^{(i)},{\mu}_{\mathrm{s}}^{(i)}]+{\alpha}^{(i)}{p}^{(i)})\le \epsilon ([{\mu}_{\mathrm{a}}^{(i)},{\mu}_{\mathrm{s}}^{(i)}])+\nu {\alpha}^{(i)}\nabla {\epsilon}^{(i)T}{p}^{(i)},$$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 ${\mu}_{\mathrm{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 ${H}^{\text{meas}}$ 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({\mu}_{\mathrm{a}}^{(l)},{\mu}_{\mathrm{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 ${\mu}_{\mathrm{a}}^{(l)}$, the error function will be corrupted by some different noise ${\sigma}^{(l)}:\epsilon ({\mu}_{\mathrm{a}}^{(l)},{\mu}_{\mathrm{s}})={\Vert {H}^{\text{meas}}-H({\mu}_{\mathrm{a}},{\mu}_{\mathrm{s}})(1+{\sigma}^{(l)})\Vert}^{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 ${\mu}_{\mathrm{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 $\frac{\partial \epsilon}{\partial {\mu}_{\mathrm{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.

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