Three-dimensional quantitative photoacoustic tomography using an adjoint radiance Monte Carlo model and gradient descent

Jens Buchmann; Bernhard A. Kaplan; Samuel Powell; Steffen Prohaska; Jan Laufer

doi:10.1117/1.JBO.24.6.066001

6 June 2019 Three-dimensional quantitative photoacoustic tomography using an adjoint radiance Monte Carlo model and gradient descent

Jens Buchmann, Bernhard A. Kaplan, Samuel Powell, Steffen Prohaska, Jan Laufer

Author Affiliations +

Journal of Biomedical Optics, Vol. 24, Issue 6, 066001 (June 2019). https://doi.org/10.1117/1.JBO.24.6.066001

Abstract

Quantitative photoacoustic tomography aims to recover maps of the local concentrations of tissue chromophores from multispectral images. While model-based inversion schemes are promising approaches, major challenges to their practical implementation include the unknown fluence distribution and the scale of the inverse problem. We describe an inversion scheme based on a radiance Monte Carlo model and an adjoint-assisted gradient optimization that incorporates fluence-dependent step sizes and adaptive moment estimation. The inversion is shown to recover absolute chromophore concentrations, blood oxygen saturation, and the Grüneisen parameter from in silico three-dimensional phantom images for different radiance approximations. The scattering coefficient is assumed to be homogeneous and known a priori.

1. Introduction

Biomedical photoacoustic (PA) tomography is a hybrid soft-tissue imaging modality that combines the high spatial resolution of ultrasound with the high contrast and specificity of optical imaging techniques.¹^–³ It relies on the generation of acoustic waves inside the tissue, which result from the absorption of intensity-modulated light, such as laser pulses or frequency chirps, by the tissue chromophores. From time-resolved PA signals recorded at multiple measurement points around the object, three-dimensional (3-D) data sets of the initial pressure distribution (i.e., PA images) can be calculated using acoustic reconstruction algorithms. Quantitative PA tomography (QPAT) aims to exploit the wavelength dependence of the image intensity to recover the local concentrations of endogenous tissue chromophores and exogenous contrast agents from which functional parameters, such as blood oxygen saturation, can be derived. To relate the PA image intensity to local chromophore concentrations, computational models of the physical processes during the image generation in conjunction with inversion schemes represent one approach to QPAT.⁴^,⁵ A major challenge in QPAT is the unknown light fluence in the tissue,⁵^–⁷ which is a nonlinear function of the concentrations and the scattering coefficient. Its effects on PA images have been described as spectral coloring and structural distortion.⁵ For an accurate quantification of concentrations and their ratios (e.g., blood oxygenation), the wavelength-dependent fluence distribution has to be accounted for.

Commonly used fluence correction methods include the application of empirical correction factors⁸ or simple models under the assumption of homogeneous optical properties.⁹ This is deemed sufficient to recover the absorption coefficient distribution from which maps of the chromophore concentrations can then be calculated using linear matrix inversions. The main limitation of these methods is the reliance on a priori knowledge of the distribution and wavelength dependence of the fluence, i.e., $ϕ (\vec{x}, λ)$ . For in vivo images, this assumption is often invalid and can lead to significant quantification errors especially at greater depth. Alternative approaches include data-driven methods. In the work reported by Tzoumas et al.¹⁰ the wavelength-dependent fluence is written on the basis of eigenspectra obtained using a principle component analysis of in silico training data, and Kirchner et al.¹¹ calculated the fluence maps by applying deep learning and local context encoding to a large number of training data. While these methods have the potential to offer fast inversions, they require large training sets and may thus lack generality.

Model-based inversions, incorporating light transport models to predict the fluence as a function of the spatial distribution of the absorption and scattering coefficients, remain the most promising approach to QPAT. The initial pressure distribution is obtained by multiplying the fluence with the distribution of the absorption coefficient and the Grüneisen parameter, and PA image data sets can be obtained using acoustic propagation models. The difference between the model output and the measured data, i.e., the objective function, is minimized by iteratively updating the absorption and scattering coefficients during the inversion until convergence is reached. To overcome the nonuniqueness that arises from the use of single-wavelength images,⁵ multi-illumination approaches¹²^,¹³ or multiwavelength image acquisition, in combination with a priori knowledge of the wavelength dependence of the absorption and scattering coefficients,¹⁴^,¹⁵ are employed.

A major challenge in high-resolution 3-D QPAT is the large number of variables ( $> 10^{6}$ ). While gradient-free methods can be applied to small-scale problems,¹⁶ inversions of a larger scale (tens of variables and more) quickly become computationally unfeasible. Gradient-based methods have the potential to overcome these limitations. They have been implemented using the adjoint formalism,¹⁷^–¹⁹ which is applied using a finite element model of the diffusion approximation (DA) to the inversion of measured 3-D phantom images.²⁰ While the DA is valid in the diffusive regime and can be implemented efficiently,²¹^,²² high-resolution PAT can cover depths in the ballistic and quasiballistic regimes, where the DA may not be valid.²³ Methods that aim to solve the radiative transfer equation (RTE) directly are computationally expensive and have only been demonstrated in two-dimensional (2-D) images so far.²⁴^–²⁷ Monte Carlo (MC) models have recently gained attention²⁸^–³¹ due to their highly parallelized architecture and advances in graphics-processing units and have already been applied to 2-D QPAT inversions¹⁹ and in initial studies with limited parameter space in 3-D.¹⁶^,³² In this paper, a method for inverting multiwavelength 3-D images based on an adjoint formulation of a radiance MC model is demonstrated in silico. The challenges that are addressed in this work are (1) the optimization of the objective function using inherently noisy gradients, (2) accounting for the effect of the concentration-dependent Grüneisen parameter, and (3) the representation of the radiance in terms of spherical harmonics. The capability of this approach to recover absolute chromophore concentrations and their ratios, e.g., blood oxygen saturation (blood ${sO}_{2}$ ), from high-resolution 3-D image data sets is demonstrated.

2. Methods

The forward model of the generation of the initial pressure shown in tomographic PA images is described in Sec. 2.1. The adjoint formalism²⁶ with which the gradients of the objective function are calculated is described in Sec. 2.2. The approximation of the radiance field as a finite sum of spherical harmonics³¹ within an MC light transport model is described in Secs. 2.3 and 2.4. The numerical phantom and the simulation parameters are outlined in Secs. 2.5 and 2.6, respectively. To reduce the impact of the inherent MC noise on the parameter update and to maximize the convergence speed of the gradient descent, an adaptive moment estimation (Adam) optimization algorithm³³ is employed (Sec. 2.7).

2.1.

Forward Model

Assuming that the effects of the limited detection aperture and acoustic propagation can be neglected, the image intensity represents the initial pressure distribution, $p_{0}$ , which is given as

Eq. (1)

p_{0} (\vec{r}, λ) = Γ (\vec{r}) H (\vec{r}),

where

Γ

is the Grüneisen parameter, which describes the PA efficiency;

H

is the absorbed optical energy density;

\vec{r}

is the spatial coordinate; and

λ

is the excitation wavelength. The absorbed energy density is defined as

Eq. (2)

H (\vec{r}) = μ_{a} (\vec{r}, λ) ϕ (\vec{r}, λ),

where

μ_{a}

is the absorption coefficient and

ϕ

is the light fluence, which is the radiance

ϕ

integrated over all angles:

Eq. (3)

ϕ (\vec{r}) = \int ϕ (\vec{r}, {\hat{s}}^{'}) d {\hat{s}}^{'} .

The absorption coefficient is related to the chromophore concentrations via the specific absorption coefficient, $α_{k} (λ)$ , i.e.,

Eq. (4)

μ_{a} (λ, \vec{r}) = \sum_{k}^{N_{k}} c_{k} (\vec{r}) α_{k} (λ),

where

N_{k}

is the number of chromophores and

k

indicates the chromophore type. The Grüneisen parameter

Γ

is assumed to be linearly dependent on chromophore concentrations,¹⁵^,³⁴^,³⁵ i.e.,

Eq. (5)

Γ (\vec{r}) = Γ_{water} [1 + \sum_{k} β_{k} c_{k} (\vec{r})],

where

β_{k}

is an empirical and chromophore-specific coefficient. The MC method is chosen for modeling the light fluence as it provides an accurate approximation of the radiative transport equation for superficial (1 to 2 cm), high-resolution QPAT.⁵ This involves the launch of photons (typically represented as packets of energy³⁶) according to a predefined source distribution. Their propagation within the domain is determined by the optical coefficients

μ_{a} (\vec{r})

and

μ_{s} (\vec{r})

, the scattering phase function

Θ (\hat{s}, {\hat{s}}^{'}, \vec{r})

, and the refractive index

n

. The deposition of energy when a photon traverses a voxel is determined by the absorption coefficient

μ_{a}

. The angular dependence of scattering events is described by the Henyey–Greenstein phase function.³⁷

2.2.

Adjoint-Assisted Optimization

Assuming Gaussian noise on the measured data, an estimate of the chromophore distributions is found by minimizing the objective function $ε$ , which is given as

Eq. (6)

ε = \sum_{l}^{N_{λ}} \int_{Ω} \frac{1}{2} {[p_{0}^{m} (λ_{l}, \vec{r}) - p_{0} (λ_{l}, \vec{r})]}^{2} d Ω,

where

Ω

is the imaged volume domain,

p_{0}^{m} (λ_{l}, \vec{r})

is the measured PA image at wavelength

λ_{l}

,

p_{0}

is the PA image obtained from the forward model, and

N_{λ}

is the number of excitation wavelengths.

To find the chromophore concentration maps, $c_{k} (\vec{r})$ , the derivative of $ε$ with respect to $c_{k}$ at any position ${\vec{r}}_{i}$ is required, i.e., $\frac{\partial ε}{\partial c_{i}} \equiv \frac{\partial ε}{\partial c_{k} ({\vec{r}}_{i})}$ . For the sake of brevity, only one chromophore will be considered in the remaining description, i.e., $\frac{\partial}{\partial c_{i}} \equiv \frac{\partial}{\partial c_{1} ({\vec{r}}_{i})}$ . The objective function $ε$ represents the sum of the objective functions, $ε_{λ_{l}}$ , at each excitation wavelength, which is given as $ε = \sum_{l}^{N_{λ}} ε_{λ_{l}}$ . For simplicity, only one wavelength, $λ_{l}$ , will be considered to describe the derivative:

Eq. (7)

\frac{\partial ε_{λ_{l}}}{\partial c_{i}} = - \int_{Ω} [p_{0}^{m} (\vec{r}) - p_{0} (\vec{r})] \frac{\partial p_{0}}{\partial c_{i}} d Ω .

Since $p_{0} = Γ μ_{a} ϕ$ , after applying the chain rule, $\frac{\partial p_{0}}{\partial c_{i}}$ becomes

Eq. (8)

\frac{\partial p_{0}}{\partial c_{i}} = \frac{\partial Γ}{\partial c_{i}} μ_{a} ϕ + Γ \frac{\partial μ_{a}}{\partial c_{i}} ϕ + Γ μ_{a} \frac{\partial ϕ}{\partial c_{i}}

Eq. (9)

= β_{k} δ (\vec{r} - {\vec{r}}_{i}) H (\vec{r}) + Γ (\vec{r}) α_{k} (λ_{l}) δ (\vec{r} - {\vec{r}}_{i}) ϕ (\vec{r}) + Γ (\vec{r}) μ_{a} (\vec{r}) \frac{\partial ϕ}{\partial c_{i}},

where

δ (\vec{r} - {\vec{r}}_{i})

is the Dirac delta function. The gradient of the fluence with respect to chromophore concentration at a particular position,

\frac{\partial ϕ}{\partial c_{i}}

, is generally unknown, and one can make use of the adjoint formalism.¹⁷^,²⁶ Briefly, the adjoint approach defines a source term

q^{*}

for the adjoint radiance

ϕ^{*}

Eq. (10)

q^{*} (\vec{r}, λ) = Γ (\vec{r}) μ_{a} (\vec{r}, λ) [p_{0}^{m} (\vec{r}, λ) - p_{0} (\vec{r}, λ)],

in which the adjoint radiance is expressed in terms of the difference between the measured images and the modeled images of

p_{0}

. Using this adjoint source definition enables the substitution of the term containing the unknown

\frac{\partial ϕ}{\partial c_{i}}

in Eq. (8) with a term containing the radiance

ϕ

and its adjoint counterpart

ϕ^{*}

Eq. (11)

\int_{Ω} Γ (\vec{r}) μ_{a} (\vec{r}) (p_{0}^{m} - p_{0}) \frac{\partial ϕ}{\partial c_{i}} d Ω = - α_{k} (λ_{l}) {[\int_{S^{2}} ϕ^{*} (\vec{s}) ϕ (\vec{s}) d \vec{s}]}_{\vec{r} = {\vec{r}}_{i}} V_{vox},

where

\int_{S^{2}}

is the integral over solid angle

\vec{s}

, and

S^{2}

is the unit sphere. The term

V_{vox}

is the result of the discretization of the data using a piecewise constant basis to sample (see Sec. 9). The gradient required to update the concentration of one chromophore at one wavelength is therefore given as

Eq. (12)

\frac{\partial ε_{λ_{l}}}{\partial c_{i}} = - [p_{0}^{m} ({\vec{r}}_{i}) - p_{0} ({\vec{r}}_{i})] V_{vox} [β_{k} H ({\vec{r}}_{i}) + Γ ({\vec{r}}_{i}) α_{k} ϕ ({\vec{r}}_{i})] + α_{k} {[\int_{S^{2}} ϕ^{*} (\vec{s}) ϕ (\vec{s}) d \vec{s}]}_{\vec{r} = {\vec{r}}_{i}} V_{vox} .

The adjoint model and its derivation are described in general and in more detail in Secs. 6 and 7.

2.3.

Radiance Approximation

As the radiance $ϕ (\vec{s})$ and the adjoint radiance $ϕ^{*} (\vec{s})$ are functions of solid angle and are defined on the surface of the unit sphere, both quantities can be expressed on the basis of spherical harmonic functions. The expansion of the radiance into spherical harmonics is based on previous work³¹ and is inspired by the $P_{n}$ approximations, described in Ref. 22 (similar to a Fourier expansion in 1-D or 2-D) and is outlined in detail in Sec. 8. Using a finite expansion of ${[\int_{S^{2}} ϕ^{*} (\vec{s}) ϕ (\vec{s}) d \vec{s}]}_{\vec{r} = {\vec{r}}_{i}}$ in real spherical harmonics, the gradient equation is given as

Eq. (13)

\frac{\partial ε_{λ_{l}}}{\partial c_{i}} = V_{vox} {- [p_{0}^{m} ({\vec{r}}_{i}) - p_{0} ({\vec{r}}_{i})] [β_{k} H ({\vec{r}}_{i}) + Γ ({\vec{r}}_{i}) α_{k} ϕ ({\vec{r}}_{i})] + α_{k} \sum_{l = 0}^{N_{L}} \sum_{m = - l}^{l} ψ_{l m} ({\vec{r}}_{i}) ψ_{l m}^{*} ({\vec{r}}_{i})},

where

ψ_{l m} ({\vec{r}}_{i})

is the radiance field approximated by the spherical harmonics function of degree

l

, order

m

at position

{\vec{r}}_{i}

and

ψ_{l m}^{*}

its adjoint counterpart. Equation (13) is implemented in a gradient-based optimization scheme (described in Sec. 2.7), which updates the concentrations iteratively to minimize the mismatch between measured and modeled data [Eq. (6)].

The last term of the gradient in Eq. (13) contains an expansion of the radiance and the adjoint radiance in spherical harmonics

Eq. (14)

α_{k} \sum_{l = 0}^{N_{L}} \sum_{m = - l}^{l} ψ_{l m} ({\vec{r}}_{i}) ψ_{l m}^{*} ({\vec{r}}_{i}) .

Here,

N_{L} = \infty

would give the most accurate radiance approximation, but due to constraints with respect to computation time and memory, the degree of the spherical harmonics,

N_{L}

, needs to be finite. However, it is not clear a priori up to which value of

N_{L}

the corresponding coefficients needs to be stored to approximate the radiance with sufficient accuracy. This has been investigated by evaluating the inversion scheme for three different configurations: (1)

N_{L} = 0

, i.e., using only the fluence, (2) for

N_{L} = 4

, i.e., the most accurate representation of the radiance, and (3) omitting the radiance, i.e.,

ψ_{l m} = ψ_{l m}^{*} = 0

for all

l

,

m

, thus neglecting the gradient term provided by the adjoint radiance. All other parameters remain the same during the inversions.

2.4.

Radiance Monte Carlo Simulations

Most MC simulators provide only the light fluence, which is the radiance integrated over all directions and time. To satisfy Eq. (12), a radiance MC algorithm (RMC)¹⁹^,³¹ is used. To obtain the radiance $ϕ (\vec{r})$ , the directional information of the photon traversing a voxel at position $\vec{r}$ is stored by depositing the photon weight into the relevant spherical harmonics coefficients $ψ_{l m} (\vec{r})$ . The RMC code is implemented in the Julia programming language,³⁸ which controls and dispatches the execution of kernels on both the CPU (written in Julia), and the GPU (written in NVIDIA’s CUDA language).³¹ Because the definition of the adjoint model does not change the dynamics of photon propagation, the same RMC simulation code provides the adjoint radiance $ϕ^{*} (\vec{r})$ .

2.5.

In Silico Phantom

An MC model has been used to calculate the multiwavelength 3-D PA images that represented measured data and are referred to as reference images or reference data throughout this paper. The domain of the model is divided into subvolumes (SVs) that represented simplified anatomical structures, such as a subcutaneous tumor and a number of discrete blood vessels. The depths of the structures are similar to those observed in in vivo images acquired using a Fabry–Perot scanner with a planar detection geometry.³⁹^–⁴² Absorption is assumed to originate from three chromophores, i.e., oxyhemoglobin ${(HbO}_{2})$ , deoxyhemoglobin (HHb), and methylene blue (Mb), as an exogenous contrast agent. It should be noted that the method can potentially be applied to any number of chromophores. The computational model of the phantom consists of nine different SVs, each with homogeneous optical properties (Fig. 1). This includes six tube-like structures to mimic blood vessels, a tumor consisting of an ellipsoidal rim and core SV, and the background. The tubes are positioned adjacent to the tumor at depths of 1.5 and 7.5 mm, have a circular cross section with a radius of 0.4 mm, and are filled with HHb and ${HbO}_{2}$ . The blood oxygen saturation ( ${sO}_{2}$ ) is defined as the ratio of the concentration of oxyhemoglobin and the total hemoglobin concentration

Eq. (15)

{sO}_{2} = \frac{c_{{HbO}_{2}}}{c_{{HbO}_{2}} + c_{HHb}} .

Fig. 1

A 3-D view of the in silico phantom and reference $p_{0}^{m}$ . (a) Segmented phantom with SV IDs; SV1: homogeneous background material, SV2 to SV4 and SV7 to SV9: tubes representing blood vessels, SV5 and SV6 represent an ellipsoidal tumor consisting of an inner and outer SV. (b) Initial pressure distribution $p_{0}^{m} (\vec{r})$ calculated using physiological hemoglobin concentrations and blood ${sO}_{2}$ at an excitation wavelength of $λ = 798 nm$ . The $x - y$ Gaussian profile of the excitation beam is shown below.

The tube-like structures are assumed to contain blood ${sO}_{2}$ ranging from 75% to 98% to represent the typical values found in veins and arteries.⁴³ The total hemoglobin concentration is 2.3 mM.⁴⁴ Two concentric ellipsoids represent the tumor at a depth ranging from 3.0 mm to 6.0 mm. The outer (inner) ellipsoid’s axes are $a = 4.5 (2.5) mm$ , $b = 3.5 (2.0) mm$ , and $c = 2.0 (0.8) mm$ , respectively. The tumor SVs contained 20% blood (0.46 mM).⁴⁵ The tumor shell has an ${sO}_{2}$ of 80%, whereas that of the core is 40% to mimic necrotic tissue. The tumor also contains an exogenous contrast agent, Mb, at a concentration of $10 μ M$ . The background SV contains a blood volume fraction of 1.5% with an ${sO}_{2}$ of 60%. Other parameters, such as the scattering anisotropy ( $g = 0.9$ ), the refractive index ( $n_{i} = 1.33$ inside the domain, $n_{e} = 1.5$ outside of the domain), and $μ_{s} (\vec{r}, λ)$ are held constant and uniform across the domain. The absorption spectra of HHb, ${HbO}_{2}$ , and Mb are shown in Fig. 2. The wavelength dependence of the reduced scattering coefficient $μ_{s}^{'} (λ) = μ_{s} (λ) (1 - g)$ is approximated using $μ_{s} (λ) = 6.65 \cdot 10^{3} \cdot λ_{+}^{- 1.317} {mm}^{- 1}$ with $λ_{+} = λ / nm$ , which resulted in a $μ_{s}^{'}$ of $\sim 1 {mm}^{- 1}$ at $λ = 800 nm$ . The Grüneisen parameter of water is set to 0.124. The coefficient $β_{k}$ describing the total hemoglobin concentration dependence of the Grüneisen parameter is set to $β_{Hb, HbO} = 0.02146 L / mmol$ .³⁵ It is assumed that the Mb concentration distribution does not change the Grüneisen parameter ( $β_{Mb} = 0$ ). The remaining parts of the volume are assumed to be filled with materials, such as water and lipids, whose absorption is considered negligible at the selected excitation wavelengths.

Fig. 2

Spectra of HHb, ${HbO}_{2}$ , and Mb used for forward simulation and inversion. The excitation wavelengths are indicated by gray vertical lines. The absorption spectra are taken from Ref. 46.

2.6.

Monte Carlo Simulation Parameters

To obtain reference images, i.e., data sets that represent measured multiwavelength PA images, the domain was discretized into $200 \times 200 \times 100$ (i.e., $4 \times 10^{6}$ ) isotropic voxels of size $V_{vox} = 10^{- 3} {mm}^{3}$ yielding a total volume of $20 \times 20 \times 10 {mm}^{3}$ . The source profile was a 2-D Gaussian function with a width of $σ = 4 mm$ . About $2 \cdot 10^{9} photon$ packets were used in the MC simulation of the light fluence. The angle-dependent radiance was not calculated. Reference image data sets were calculated for three different excitation wavelengths that coincided with the absorption peaks of Mb (664 nm), HHb (758 nm), and the isosbestic point of hemoglobin absorption (798 nm). Gaussian noise ( $σ = 0.1 %$ of the maximum image intensity) was added to the reference data, resulting in negative image intensities in regions of low $p_{0}$ .

The domain discretization used during the inversion was identical to that used for the reference data set, which may raise the question of whether this constitutes a so-called inverse crime. In MC models, the discretization is used merely as a basis for sampling physical quantities while photon packets can propagate freely in continuous space. This is in contrast to other methods, such as finite elements, where the discretization has a direct impact on the accuracy of the solution. Taking also into account the stochastic nature of MC models, it can be concluded that using identical discretizations does not constitute an inverse crime.

During an inversion, $1 \times 10^{7} photon$ packets were used for the calculation of the radiance and fluence; $5 \times 10^{6} photons$ were used for the calculation of the corresponding adjoint quantities. The typical running time for one inversion iteration, including the adjoint model with $N_{L} = 4$ , was 84 s on a high-end consumer GPU (NVIDIA GeForce Titan X Pascal). This was reduced to 17 s when the radiance term in Eq. (13) was neglected, i.e., $N_{L} = 0$ and no adjoint MC runs. Since $k = 3$ independent chromophore concentrations were associated with each voxel, the model contained a total of 12 million variables.

2.7.

Gradient-Based Optimization

The gradient-based optimization was initialized assuming a homogeneous $c_{HHb} = c_{HbO} = 0.023 mM$ , i.e., ${sO}_{2} = 50 %$ , and $c_{Mb} = 0.0 mM$ . Owing to the stochastic nature of MC models, the gradients [Eq. (13)] were subjected to noise. To compensate for this, the Adam optimization algorithm³³ was employed. It was developed for the first-order optimization of noisy objective functions in high-dimensional parameter spaces and can be seen as an extension of the momentum algorithm.⁴⁷^,⁴⁸ The Adam algorithm calculated an exponential moving average of the gradient (first moment) and the squared gradient (second moment) using the decay rates $β_{1, Adam}$ and $β_{2, Adam}$ with an additional bias-correction step. The final update was calculated by multiplying the step-size parameter $γ_{Adam}$ with the first-order moment divided by the square root of the second-order moment. The detailed description can be found in Ref. 33. This algorithm was found to dramatically increase the convergence speed, compared to standard gradient descent. Its efficiency depended on the values of a set of parameters, including the step-size $γ_{Adam}$ , the decay rates $β_{1, Adam}$ and $β_{2, Adam}$ , and an additional $ε_{Adam}$ , which avoided division by zero. The decay rates and $ε$ were set to recommended default values ( $β_{1, Adam} = 0.9$ , $β_{2, Adam} = 0.999$ , and $ε_{Adam} = 10^{- 8}$ ), and the step size was assigned different values depending on the chromophore according to Eqs. (16) and (17).

To ensure fast convergence for all chromophores, the step size for Mb is set to be significantly smaller than that of $HHb / {HbO}_{2}$ since Mb concentrations are in the range of $μ M$ while those of $HHb / {HbO}_{2}$ are in the range of mM. The chromophore-dependent step size is calculated as follows:

Eq. (16)

γ_{chrom} = γ_{ref} \frac{\sum_{λ} α_{chrom, λ} \cdot c_{\max, ref}}{\sum_{λ} α_{ref, λ} \cdot c_{\max, chrom}},

where

γ_{ref}

is the step size of a reference chromophore, the value of which was set ad hoc to

γ_{ref} = γ_{{HbO}_{2}} = 100

;

α_{chrom / ref, λ}

is the specific absorption coefficient of the respective chromophore at wavelength

λ

; and

c_{\max, ref / chrom}

is the anticipated maximum concentration of the respective chromophore, which is set to physiological reasonable values (

c_{\max, HbO} = c_{\max, HHb} = 2.3 mM

,

c_{\max, Mb} = 30 μ M

).

The gradient is also expressed as a function of fluence in Eq. (13), either directly or via $H$ and $Δ p$ . This would result in slow convergence in regions of low fluence. To compensate for this, a spatially dependent step size is used, which increases the step size by normalizing it by the mean fluence over all wavelengths:

Eq. (17)

γ {(\vec{r})}_{chrom, scaled} = \frac{γ_{chrom}}{{\tilde{ϕ}}_{norm} (\vec{r}) + ε_{ϕ}} .

Note that

{\tilde{ϕ}}_{norm} (\vec{r})

is the normalized mean fluence:

Eq. (18)

{\tilde{ϕ}}_{norm} (\vec{r}) = \sum_{λ} ϕ (\vec{r}, λ) / \sum_{λ} ϕ ({\vec{r}}_{\max}, λ),

where

{\vec{r}}_{\max} = {\arg \max}_{\vec{r}} \sum_{λ} ϕ (\vec{r}, λ)

represents the location where the total fluence is at a maximum. The parameter

ε_{ϕ}

avoids division by zero and determines the maximum change in the step size in regions of low fluence. In this study,

ε_{ϕ} = 10^{- 4}

leads to a sufficient speed-up in convergence for the deepest tubes. The optimization of the step size with respect to the different chromophores and the local fluence is often referred to as preconditioning.

A single iteration of the gradient-based update consisted of the execution of the MC model and its adjoint counterpart for all three wavelengths. The inversion was run for 1500 iterations to investigate the convergence. After each iteration, the updated concentrations for HHb, ${HbO}_{2}$ , and Mb were limited to a reasonable range of values (also known as projected gradient descent) to avoid spurious overshooting or undershooting that could lead to physiologically unrealistic concentrations. To compensate for the effects of noise in low-absorbing regions (i.e., negative $p_{0}$ in reference images), negative chromophore concentrations, and hence negative $μ_{a}$ values, were allowed during the gradient descent. To ensure stability, the range of negative concentration values was limited to a tenth of the maximum positive concentrations ( $- 0.35$ to 3.5 mM for $c_{HHb}$ and $c_{{HbO}_{2}}$ , $- 0.02$ to 0.2 mM for $c_{Mb}$ ). The 3-D blood ${sO}_{2}$ maps were calculated from the recovered $c_{{HbO}_{2}}$ and $c_{HHb}$ images. The scattering distribution was assumed to be known a priori.

To verify that the inversion scheme is valid over a range of physiologically plausible parameter values, multiwavelength reference images were calculated and inverted for different ${sO}_{2}$ values in the inner tumor region and the background. Two scenarios were chosen, each comprising five different combination of ${sO}_{2}$ values. First, as the tumor core (being completely enclosed by the rim) may be assumed to be most strongly affected by spectral coloring, its ${sO}_{2}$ value was varied from 10% to 90% in 20% increments (SV 6), whereas all other parameters remained fixed. Second, the ${sO}_{2}$ value of the background (SV ID 1) was varied from 10% to 90% in 20% increments.

3. Results

The accuracy of the recovered chromophore concentration and ${sO}_{2}$ maps are reported in Secs. 3.1 and 3.2. In Sec. 3.3, the results obtained using multiple inversions of image data sets, in which the chromophore concentrations and their ratios are varied, are reported.

3.1.

Absolute Concentrations

Examples of 3-D volume-rendered images of the absolute concentrations of ${HbO}_{2}$ , HHb, and Mb recovered after 1500 iterations are shown in Figs. 3(a)–3(c). The color scales are thresholded to render the background with its comparatively low chromophore concentrations transparent. To reduce the effect of the added Gaussian noise on the rendered images (particularly in regions of low fluence), a 3-D median filter is applied (non-iterative, edge- and face-connected). The error functional as a function of the number of iterations is shown in Fig. 3(d). The value of the error functional cannot reach zero due to the inherent MC noise of the forward model used during the inversion. The dashed orange line indicates the minimum value of the error functional that is reached with $1 \times 10^{7} photons$ , and the green dotted line indicates the total noise level consisting of MC noise and Gaussian noise added to the reference data. The MC noise is obtained from forward calculations with prior knowledge of the correct chromophore distributions.

Fig. 3

Absolute chromophore concentration maps recovered using the inversion scheme and its convergence. (a)–(c) The 3-D volume-rendered images of the concentrations of ${HbO}_{2}$ , (b) HHb, and (c) Mb. (d) Values of the error functional (solid blue line), the baseline of the MC noise (dashed orange line), and added Gaussian noise (dotted green line).

In Fig. 4, cross-sectional $x z$ -images of the true and recovered concentration of the three chromophores and the absolute error at $y = 10 mm$ (center plane) are shown. Excellent agreement is found in regions of high fluence [e.g., corresponding to white, light blue, and light red pixels in Figs. 4(c), 4(f), and 4(i)]. By contrast, significant quantification errors [corresponding to yellow and green pixels in Figs. 4(c), 4(f), and 4(i)] are observed in regions of low fluence. The recovered $c_{Mb}$ appears to exhibit larger quantification errors in the background compared to $c_{{HbO}_{2}}$ and $c_{HHb}$ .

Fig. 4

The 2-D cross-sectional $x z$ -images of the true and recovered chromophore concentrations, together with the absolute error at $y = 10 mm$ . Left column: true chromophore concentrations. Center column: recovered concentrations. Left and center column share the same color bar. Right column: absolute concentration error. Voxels where the error exceeds the limits of the color scale are rendered in green and yellow. In (a), the reduced background SV 1* is illustrated as a white dashed rectangle. Because background voxels far away from the source exhibit large errors due to low SNR, this reduced SV is used for calculating average concentrations.

Table 1 contains the true and recovered concentrations averaged over all voxels of each SV defined in Sec. 2.5. The brackets indicate the standard deviation (concise notation). The recovered $c_{HHb}$ and $c_{{HbO}_{2}}$ in SV 2 to 9 are in excellent agreement with the true values. While the recovered $c_{Mb}$ are generally in good agreement with the true values, SV 2 to 4 and SV 7 to 9 exhibit negative and small concentrations with large standard deviations. The background (SV 1) also exhibits large errors and standard deviations due to the low signal-to-noise ratio (SNR). In the background region closer to the light source [ $140 \times 140 \times 70$ central voxels inside the $200 \times 200 \times 100$ image domain, illustrated in Fig. 4(a)] where the SNR is greater, the recovered concentrations are in good agreement with the true values (SV 1* in Table 1).

Table 1

True and recovered (inversion) absolute chromophore concentrations and blood sO2 in the SVs of the phantom. The concentrations represent the average value over all voxels within each SV; the values in brackets indicate the standard deviations (concise notation). SV 1—background, SV 1*—background close to the source, SV 5 and SV 6—outer and inner tumor SVs, respectively, SV 2 to 4 and SV 7 to 9—blood-filled tubes.

	HHb (mM)		HbO2 (mM)		Mb (μM)		sO2 (%)
SV	True	Inversion	True	Inversion	True	Inversion	True	Inversion
1	0.0138	0.0087(1169)	0.0207	0.0320(212)	0	0.2(58)	60	78.3
1^*	0.0138	0.0138(219)	0.0207	0.0205(384)	0	$- 0.002 (922)$	60	59.8
2	0.506	0.503(70)	1.79	1.78(11)	0	$- 0.08 (266)$	78	77.9
3	0.115	0.117(63)	2.18	2.16(11)	0	$- 0.04 (244)$	95	94.9
4	0.046	0.047(31)	2.25	2.25(6)	0	$- 0.05 (106)$	98	97.9
5	0.092	0.092(18)	0.368	0.367(31)	10	9.76(99)	80	79.9
6	0.276	0.275(18)	0.184	0.183(29)	10	9.51(100)	40	39.8
7	0.575	0.560(128)	1.73	1.71(20)	0	$- 0.35 (510)$	75	75.3
8	0.230	0.229(90)	2.07	2.04(14)	0	$- 0.41 (374)$	90	89.9
9	0.460	0.450(125)	1.84	1.82(20)	0	$- 0.37 (509)$	80	80.0

As described in Sec. 2.4, the inversion is implemented using an approximation of the radiance based on spherical harmonics of varying degree $N_{L}$ , including the omission of the adjoint term. The inversions are found to converge to the same final values while propagating along different routes. Convergence is reached irrespective of the radiance approximation (see Fig. 5).

3.2.

Blood Oxygen Saturation

Figure 6 shows the cross-sectional $x z$ -images ( $y = 10 mm$ ) of the known and recovered blood ${sO}_{2}$ together with the absolute error. The accuracy is clearly affected by noise as shown in the difference image in Fig. 6(c). While most voxels in SV 2 to 9 exhibit an error within $\pm 5 %$ ${sO}_{2}$ , it is noticeably larger for objects at greater depth. Concentrations in the background SV 1 are also affected by noise, particularly in regions further away from the source. However, the accuracy improves near the source where SNR is increased. The average values of blood ${sO}_{2}$ for each SV are also calculated and are summarized in Table 1. Blood ${sO}_{2}$ in SV 2 to 9, i.e., corresponding to blood-filled tubes and the tumor SVs, are found to lie within 0.3% of the true values (i.e., while ${sO}_{2}$ errors of individual voxels can be quite large due to noise, averaging over lots of voxels greatly improves accuracy). The average blood ${sO}_{2}$ for SV 1 (i.e., the entire background SV) is found to show significant errors (18.3% ${sO}_{2}$ ) and is attributed mainly to the adverse effects of noise in low fluence regions. By contrast, the inversion results are accurate for the reduced background SV 1* due to higher SNR.

Fig. 5

Convergence of the true and recovered concentrations for inversions incorporating the gradient term of the adjoint radiance for $N_{L} = 4$ (dotted line, ⋯), $N_{L} = 0$ (dashed-dotted line, $- \cdot -$ ), and without the adjoint term (dashed line, $- -$ ). The chromophores are shown in blue (Mb), red ( ${HbO}_{2}$ ), and black (HHb). Horizontal lines in light color depict perfect convergence, i.e., $Δ c = 0$ .

Fig. 6

The 2-D cross-sectional $y z$ -images of the true (a) and recovered ${sO}_{2}$ (b), together with the absolute error (c) at the center of the phantom ( $y = 10 mm$ ). The images on the left and in the center share the same color bar. In the right color bar, voxels exceeding (dropping below) an absolute difference between true and recovered ${sO}_{2}$ of 15% are shown in green (yellow).

3.3.

Validation over a Range of Blood Oxygen Saturation

The inversion scheme is validated on image data sets where ${sO}_{2}$ is varied over a range of physiologically plausible parameter values (Sec. 2.7). The inversions are computed without including the gradient term of the radiance and $N_{L} = 0$ in order to minimize the computation time. To obtain the final ${sO}_{2}$ value for each image data set, the inversion is run for 1500 iterations after which the average ${sO}_{2}$ is obtained from the SVs. Figure 7(a) shows the true and recovered ${sO}_{2}$ values for all SVs and all image data sets together with the line of unity (dashed line) and a $\pm 5 %$ error interval (dotted lines). All recovered ${sO}_{2}$ values are in good agreement with the known values and exhibit an average error below 0.3% ${sO}_{2}$ . Figure 7(b) shows the difference between true and recovered ${sO}_{2}$ for all SVs and image data sets sorted by SV. Only the results corresponding to the reduced background SV 1* are shown as this region exhibits sufficient SNR.

Fig. 7

Validation of the inversion scheme over a range of blood ${sO}_{2}$ . (a) Average recovered ${sO}_{2}$ as a function of the true ${sO}_{2}$ of all SVs and image data sets together with the line of unity (dashed line) and a $\pm 5 %$ error interval (dotted lines). In the five image data sets, the ${sO}_{2}$ in the inner tumor material varies from 10% to 90% in 20% increments. In another five image data sets, the background ${sO}_{2}$ ranges from 10% to 90% (20% increments). (b) Box-and-whisker-plot of the absolute difference between true and recovered ${sO}_{2}$ for each SV.

4. Discussion

The 3-D maps of absolute concentrations of ${HbO}_{2}$ , HHb and Mb, and the resulting blood ${sO}_{2}$ recovered using a gradient-based MC inversion scheme showed excellent agreement with the true values. To achieve the best possible match of noise-affected PA images and the model, the inversion scheme was implemented without a non-negativity constraint for the chromophore concentrations. Even though negative concentrations were physiologically implausible, it was found that incorporating a non-negativity constraint greatly affected the recovered average concentrations and blood ${sO}_{2}$ values. However, negative concentrations can lead to negative absorption coefficients. In the MC model, it led to the photon packets gaining weight as they traversed a voxel with negative absorption. If too many voxels exhibited negative $μ_{a}$ , unstable inversions can be observed as the photon weight diverges. While this was occasionally observed in this study, it was found that a reduction of the step size and an increase in the number of iterations remedied this problem. The inversion scheme described in this paper included an expression of the radiance and its adjoint in a basis of spherical harmonics. The influence of the adjoint formalism and the spherical harmonics approximation on the accuracy and convergence of the inversion was evaluated under the assumption that the scattering coefficient was known a priori. It was found that neither accuracy nor convergence speed were affected by the radiance term and its adjoint, i.e., the last term in Eq. (13). This was also observed when the radiance term was omitted (Fig. 5) and was confirmed by the relative magnitudes of the individual gradient terms. The radiance term (irrespective of the spherical harmonic approximations) was always significantly smaller than the remaining terms of Eq. (13). Omitting the computation of the adjoint radiance resulted in a major increase in computational speed. However, from the limited investigation presented here, it can only be concluded that the adjoint term may be neglected if the scattering coefficient was known. If the recovery of the scattering coefficient was of interest, a radiance approximation to a minimum of $N_{L} = 1 \deg$ may be necessary.³¹

The gradient-based inversion was found to benefit greatly from optimization algorithms in which parameters, such as the step sizes and the exponential decay rates of the Adam algorithm, were predefined to enable a fast and accurate convergence. A potential drawback of such methods was the need to test several sets of these parameters prior to an inversion to assess whether they have a positive impact on the convergence speed. Within the scope of this study, only minor and ad hoc parameter tuning was conducted. A more thorough investigation, including the development of automated parameter selection algorithms, may yield significantly faster convergence.

The chromophore-dependent step size and the fluence-dependent spatial step size scaling [Eqs. (16) and (17)] proved to be vital for achieving convergence. Without chromophore-dependent step sizes, Mb concentrations diverged to the upper and lower fit limits. Similarly, the fluence-dependent spatial step size scaling was crucial for achieving fast convergence in the regions of low fluence.

The selection of the excitation wavelengths could also be optimized further to improve inversion accuracy and convergence speed.⁴⁹ However, such a study would exceed the scope of this paper. Despite potentially suboptimal excitation wavelengths, the inversion is shown to recover blood ${sO}_{2}$ over a wide range [Fig. 7(a)] with high accuracy ( $< 1 %$ error in ${sO}_{2}$ ) across the domain.

Gradient-based methods do not guarantee convergence to a global minimum, especially when the inversion is adversely affected by a noisy gradient. While the Adam optimization algorithm (compared to, for example, standard or momentum gradient descent) has been shown to greatly reduce the likelihood of finishing the inversion in a local minimum or on a saddle point, such a result cannot be ruled out entirely. It should also be noted that the application of this method to measured PA images does not require their segmentation into subregions. While this makes the method generally valid, some form of image segmentation may still be advantageous as it would reduce the number of variables and the risk of convergence to a local minimum and increase convergence speed. Moreover, explicit regularization of the objective function could further improve the accuracy and speed of the convergence.

While the general methodology of a gradient-based inversion using an adjoint formulation of an MC model has been demonstrated in silico, the application of this approach to experimental 3-D PA images, especially those acquired in vivo, requires further investigation. One of the perhaps most critical points is the recovery of the scattering coefficient as it is likely to have an impact on the importance of the radiance approximation using spherical harmonics. Other issues, such as the choice of inversion parameters and the selection of optimal excitation wavelengths, are also important in the translation of QPAT methods toward applications in the medical and life sciences.

5. Conclusions

An inversion scheme for recovering absolute chromophore concentrations and their ratios, such as blood ${sO}_{2}$ from 3-D multiwavelength PA images was developed and validated in silico. The scheme was based on an adjoint formulation of an MC light-transport model and allowed an approximation of the radiance using spherical harmonics. It was found that the adjoint radiance was not required to obtain accurate inversion results, provided the scattering coefficient was constant. The speed of convergence was increased by incorporating the Adam optimization algorithm, chromophore-dependent step sizes, and fluence-dependent step-size scaling. This work represented an important step in the development of robust and generally applicable methods for quantitative functional and molecular PA imaging.

6. Appendix A: Definition of the Adjoint Model

The idea of the adjoint formalism is to define nonphysical quantities, an adjoint source $q^{*} (\vec{r}, λ)$ , adjoint radiance $ϕ^{*} (\vec{r}, \vec{s}, λ)$ , and an adjoint fluence $ϕ^{*} (\vec{r}, λ) = \int_{S^{2}} ϕ^{*} (\vec{r}, \vec{s}, λ) d \vec{s}$ that help in replacing the integral term containing the unknown $\frac{\partial ϕ}{\partial c_{i}}$ in the definition of the gradient. In our case the gradient equation is

Eq. (19)

\frac{\partial ε_{λ_{l}}}{\partial c_{i}} = - [p_{0}^{m} ({\vec{r}}_{i}) - p_{0} ({\vec{r}}_{i})] V_{vox} [β_{k} H ({\vec{r}}_{i}) + Γ ({\vec{r}}_{i}) α_{k} ϕ ({\vec{r}}_{i})] + \int_{Ω} (p_{0}^{m} - p_{0}) Γ (\vec{r}) μ_{a} (\vec{r}) \frac{\partial ϕ}{\partial c_{i}} d Ω .

The adjoint approach has been used in the context of PAT earlier, see e.g., Refs. 17 18.–19, 26, and 27.

As a first step, the adjoint source term is defined as

Eq. (20)

q^{*} (\vec{r}, λ) = [p_{0}^{m} (\vec{r}, λ) - p_{0} (\vec{r}, λ)] Γ (\vec{r}) μ_{a} (\vec{r}, λ) .

Since the approach is targeted for multispectral QPAT, each wavelength requires its own definition of an adjoint source based on the difference between modeled $p_{0} (\vec{r}, λ)$ and measured data $p_{0}^{m} (\vec{r}, λ)$ . We only denote one wavelength here and omit the dependence on $λ$ for the sake of brevity.

The adjoint source is usually defined as the “prefactor” of the unknown term that is to be replaced and contains the error between modeled and measured data. By defining the behavior of the adjoint radiance $ϕ^{*}$ and adjoint fluence $ϕ^{*}$ , a relationship between these and the desired $\frac{\partial ϕ}{\partial c_{i}}$ can be established, which is outlined in the following derivation. One important aspect of the definition of the adjoint quantity is to leave the equations underlying the development similar to the ones of their physical counterpart, which is in this context the time-independent RTE.

The RTE is given as

Eq. (21)

[\vec{s} \cdot \nabla + μ_{a} (\vec{r}) + μ_{s} (\vec{r})] ϕ (\vec{r}, \vec{s}) - μ_{s} \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) ϕ (\vec{r}, {\vec{s}}^{'}) d {\vec{s}}^{'} = q (\vec{r}, \vec{s}) .

Similarly, we define the adjoint RTE (ARTE) as

Eq. (22)

[- \vec{s} \cdot \nabla + μ_{a} (\vec{r}) + μ_{s} (\vec{r})] ϕ^{*} (\vec{r}, \vec{s}) - μ_{s} \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) ϕ^{*} (\vec{r}, {\vec{s}}^{'}) d {\vec{s}}^{'} = q^{*} (\vec{r}, \vec{s}) .

One advantage of defining the adjoint radiance in that way is that the propagation dynamics are identical to that of the normal radiance defined by the RTE since the left-hand side of both the RTE and the ARTE are practically identical, the only difference being the negative sign in the ARTE indicating a change of direction in light propagation. One way to interpret the negative sign is to follow the propagation of photons in the opposite direction, which does not affect the photon’s movement and the final results in terms of energy deposit. Thus, the mechanisms for absorption and scattering remain the same as in the RTE. Because the light transport is dominated by scattering and absorption and not whether photons move in the forward or backward direction, light propagation can be seen as reciprocal and hence the numerical framework implementing the ARTE is unaffected by the additional negative sign. Thus, the same simulation code as for the RTE can be used for the ARTE. Only the difference in the source distributions needs to be taken into account, with the adjoint source being 3-D, whereas normal source distributions are usually 2-D.

It is important to note that the adjoint fluence and adjoint radiance have different physical units as their physical counterparts. The adjoint radiance’s unit is ${J / (m}^{3} sr)$ and the adjoint fluence’s unit is ${J / (m}^{3})$ .

7. Appendix B: Adjoint-Assisted Derivation of the Gradient

Using the above definition of the adjoint source in Eq. (20), the substitution of the unknown term $\frac{\partial ϕ}{\partial c_{i}}$

Eq. (23)

\int_{Ω} (p_{0}^{m} - p_{0}) Γ (\vec{r}) μ_{a} (\vec{r}) \frac{\partial ϕ}{\partial c_{i}} d Ω = - α_{k} {[\int_{S^{2}} ϕ^{*} (\vec{s}) ϕ (\vec{s}) d \vec{s}]}_{\vec{r} = {\vec{r}}_{i}} V_{vox},

can be derived.

The derivation follows the ideas presented in previous works, in particular Ref. 27. In Ref. 27, $\frac{\partial RTE}{\partial c_{i}}$ is combined with the ARTE

Eq. (24)

ϕ^{*} \cdot \frac{\partial RTE}{\partial c_{i}} - \frac{\partial ϕ}{\partial c_{i}} \cdot ARTE .

The term $\frac{\partial RTE}{\partial c_{i}}$ is

Eq. (25)

[\vec{s} \cdot \nabla + μ_{a} (\vec{r}) + μ_{s} (\vec{r})] \frac{\partial ϕ (\vec{r}, \vec{s})}{\partial c_{i}} + \frac{\partial μ_{a}}{\partial c_{i}} ϕ (\vec{r}, \vec{s}) - μ_{s} \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) \frac{\partial ϕ (\vec{r}, {\vec{s}}^{'})}{\partial c_{i}} d {\vec{s}}^{'} = 0,

as

\frac{\partial q}{\partial c_{i}} = 0

since the external light source does not depend on chromophore concentrations. The basic idea underlying the following steps is to rearrange all the terms in Eq. (24) so that all terms on the left-hand side can be set to zero after integrating over space and angles. First, we insert

\frac{\partial μ_{a}}{\partial c_{i}} = α δ (\vec{r} - {\vec{r}}_{i})

and rearrange Eq. (25)

Eq. (26)

(\vec{s} \cdot \nabla + μ_{a} (\vec{r}) + μ_{s} (\vec{r})) \frac{\partial ϕ (\vec{r}, \vec{s})}{\partial c_{i}} - μ_{s} \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) \frac{\partial ϕ (\vec{r}, {\vec{s}}^{'})}{\partial c_{i}} d {\vec{s}}^{'} = - α (\vec{r}) ϕ (\vec{r}, \vec{s}) δ (\vec{r} - {\vec{r}}_{i}) .

The combination of the ARTE and the RTE from Eq. (24) is (for brevity we omit the dependency on $\vec{r}$ and $\vec{s}$ for the moment)

Eq. (27)

ϕ^{*} (\vec{s} \cdot \nabla + μ_{a} + μ_{s}) \frac{\partial ϕ}{\partial c_{i}} - ϕ^{*} μ_{s} \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) \frac{\partial ϕ}{\partial c_{i}} d {\vec{s}}^{'} - \frac{\partial ϕ}{\partial c_{i}} (- \vec{s} \cdot \nabla + μ_{a} + μ_{s}) ϕ^{*} + \frac{\partial ϕ}{\partial c_{i}} μ_{s} \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) ϕ^{*} d {\vec{s}}^{'} = - α ϕ^{*} ϕ δ (\vec{r} - {\vec{r}}_{i}) - \frac{\partial ϕ}{\partial c_{i}} q^{*} .

The left-hand side can be simplified and can be given as

Eq. (28)

ϕ^{*} (\vec{s} \cdot \nabla) \frac{\partial ϕ}{\partial c_{i}} + \frac{\partial ϕ}{\partial c_{i}} (\vec{s} \cdot \nabla) ϕ^{*} - ϕ^{*} μ_{s} \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) \frac{\partial ϕ}{\partial c_{i}} d {\vec{s}}^{'} + \frac{\partial ϕ}{\partial c_{i}} μ_{s} \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) ϕ^{*} d {\vec{s}}^{'} = - α ϕ^{*} ϕ δ (\vec{r} - {\vec{r}}_{i}) - \frac{\partial ϕ}{\partial c_{i}} q^{*} .

The next steps of the derivation are from now on identical to Eqs. (4.11 ff) in Ref. 27. The left-hand side of Eq. (28) equates to zero, which can be seen after integrating first over all angles $s \in S^{n - 1}$ and over the volume $Ω$ with surface $\partial Ω$ :

Eq. (29)

\int_{Ω} \int_{S^{2}} ϕ^{*} (\vec{s} \cdot \nabla) \frac{\partial ϕ}{\partial c_{i}} d \vec{s} d Ω + \int_{Ω} \int_{S^{2}} \frac{\partial ϕ}{\partial c_{i}} (\vec{s} \cdot \nabla) ϕ^{*} d \vec{s} d Ω - \int_{Ω} μ_{s} \int_{S^{2}} ϕ^{*} (\vec{s}) \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) \frac{\partial ϕ ({\vec{s}}^{'})}{\partial c_{i}} d {\vec{s}}^{'} d \vec{s} d Ω + \int_{Ω} μ_{s} \int_{S^{2}} \frac{\partial ϕ (\vec{s})}{\partial c_{i}} \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) ϕ^{*} ({\vec{s}}^{'}) d {\vec{s}}^{'} d \vec{s} d Ω = - \int_{Ω} α δ (\vec{r} - {\vec{r}}_{i}) \int_{S^{2}} ϕ^{*} ϕ d \vec{s} d Ω - \int_{Ω} q^{*} \int_{S^{2}} \frac{\partial ϕ}{\partial c_{i}} d \vec{s} d Ω .

To see that the left-hand side equates to zero the volume integral with the terms involving $\nabla$ can be transformed into a surface integral using this form of the divergence theorem:

Eq. (30)

\int_{Ω} a b \cdot \nabla c d Ω + \int_{Ω} c b \cdot \nabla a d Ω = \int_{\partial Ω} b \cdot \hat{n} a c d Ω,

along with the following substitutions:

Eq. (31)

a = ϕ^{*} (\vec{s}), b = \vec{s} and c = \frac{\partial ϕ (\vec{s})}{\partial c_{i}} .

Hence, using the divergence theorem, the first two terms in Eq. (28) are replaced by a single term

\int_{Ω} \int_{S^{2}} ϕ^{*} (\vec{s} \cdot \nabla) \frac{\partial ϕ}{\partial c_{i}} d \vec{s} d Ω + \int_{Ω} \int_{S^{2}} \frac{\partial ϕ}{\partial c_{i}} (\vec{s} \cdot \nabla) ϕ^{*} d \vec{s} d Ω = \int_{\partial Ω} \int_{S^{2}} (\vec{s} \cdot \hat{\vec{n}}) ϕ^{*} (\vec{s}) \frac{\partial ϕ (\vec{s})}{\partial c_{i}} d \vec{s} d Ω .

By definition, both $ϕ^{*} \to 0$ and $\partial ϕ / \partial c_{i} \to 0$ on the boundary of the volume $\partial Ω$ . Thus, the integrand on the right-hand side and hence the integral equate to zero.

This reduces Eq. (28) to

Eq. (32)

- \int_{Ω} μ_{s} \int_{S^{2}} ϕ^{*} (\vec{s}) \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) \frac{\partial ϕ ({\vec{s}}^{'})}{\partial c_{i}} d {\vec{s}}^{'} d \vec{s} d Ω + \int_{Ω} μ_{s} \int_{S^{2}} \frac{\partial ϕ (\vec{s})}{\partial c_{i}} \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) ϕ^{*} ({\vec{s}}^{'}) d {\vec{s}}^{'} d \vec{s} d Ω = - \int_{Ω} α δ (\vec{r} - {\vec{r}}_{i}) \int_{S^{2}} ϕ^{*} ϕ d \vec{s} d Ω - \int_{Ω} q^{*} \int_{S^{2}} \frac{\partial ϕ}{\partial c_{i}} d \vec{s} d Ω .

Because we can assume that all functions are integrable, the terms on the left-hand side of Eq. (32) can be rearranged after changing the order of integration, yielding

\int_{Ω} μ_{s} \int_{S^{2}} \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) ϕ^{*} (\vec{s}) \frac{\partial ϕ ({\vec{s}}^{'})}{\partial c_{i}} d {\vec{s}}^{'} d \vec{s} d Ω - \int_{Ω} μ_{s} \int_{S^{2}} \int_{S^{2}} Θ (\vec{s}, {\vec{s}}^{'}) ϕ^{*} ({\vec{s}}^{'}) \frac{\partial ϕ (\vec{s})}{\partial c_{i}} d {\vec{s}}^{'} d \vec{s} d Ω .

Hence, the left-hand side of Eq. (32) equates to zero, which leaves us with

Eq. (33)

\int_{Ω} q^{*} \int_{S^{2}} \frac{\partial ϕ}{\partial c_{i}} d \vec{s} d Ω = - \int_{Ω} α δ (\vec{r} - {\vec{r}}_{i}) \int_{S^{2}} ϕ^{*} ϕ d \vec{s} d Ω .

The fluence $ϕ$ is by definition the integral of the time-integrated radiance $ϕ (\vec{s})$ over all directions $\vec{s}$ , that is,

Eq. (34)

ϕ (\vec{r}) = \int_{S^{2}} ϕ (\vec{r}, \vec{s}) d \vec{s} .

Applying the derivative with respect to $c_{i}$ yields

Eq. (35)

\frac{\partial ϕ (\vec{r})}{\partial c_{i}} = \int_{S^{2}} \frac{\partial ϕ (\vec{r}, \vec{s})}{\partial c_{i}} d \vec{s},

which lets us write Eq. (33) as

Eq. (36)

\int_{Ω} q^{*} \frac{\partial ϕ}{\partial c_{i}} d Ω = - α ({\vec{r}}_{i}) {[\int_{S^{2}} ϕ^{*} (\vec{s}) ϕ (\vec{s}) d \vec{s}]}_{\vec{r} = {\vec{r}}_{i}} V_{vox} .

Because we have defined the adjoint source term as $q^{*} = (p_{0}^{m} - p_{0}) Γ μ_{a}$ the left-hand side of Eq. (36) is exactly the last term in Eq. (19)

Eq. (37)

\int_{Ω} (p_{0}^{m} - p_{0}) Γ μ_{a} \frac{\partial ϕ}{\partial c_{i}} d Ω = - α ({\vec{r}}_{i}) {[\int_{S^{2}} ϕ^{*} (\vec{s}) ϕ (\vec{s}) d \vec{s}]}_{\vec{r} = {\vec{r}}_{i}} V_{vox} .

Inserting this result into the error-gradient Eq. (19) provides us with the subgradient term obtained using one wavelength

Eq. (38)

\frac{\partial ε_{λ_{l}}}{\partial c_{i}} = - [p_{0}^{m} ({\vec{r}}_{i}) - p_{0} ({\vec{r}}_{i})] α (λ_{l}) ϕ ({\vec{r}}_{i}) V_{vox} + α ({\vec{r}}_{i}) {[\int_{S^{2}} ϕ^{*} (\vec{s}) ϕ (\vec{s}) d \vec{s}]}_{\vec{r} = {\vec{r}}_{i}} V_{vox} .

In summary, the adjoint formalism enables the update of the concentration distribution using only terms obtained from running the forward model implemented by the RMC algorithm. The following section will focus on the question as to how the term $ϕ^{*} ϕ$ can be computed and approximated.

8. Appendix C: Radiance Approximations

To simulate the radiance, some discretization over angle is required. One option is to use a set of piecewise constant basis functions. However, due to the importance of ballistic and quasi-ballistic light propagation in PAT, a high number of discretization orders would be required to capture the directionality of the radiance in regions near the source. This would result in very large memory demands in finite-element implementations (see Refs. 29, 50, and 51, for more details on different approaches for estimating the radiance, their limitations, and suggested solutions. For a summary thereof, see Ref. 52). Inspired by the $P_{n}$ approximation²² and continuing the work presented by Refs. 19 and 52 we approximate the radiance in 3-D using spherical harmonics as basic functions, as in Ref. 31. Instead of discretizing the angular domain into segments, the radiance field $ϕ$ at any position $\vec{r}$ can be expanded using a series of spherical harmonics⁵³^,⁵⁴

Eq. (39)

ϕ (\vec{r}, \vec{s}) = \sum_{l = 0}^{N_{L} = \infty} \sum_{m = - l}^{l} ψ_{l m} (\vec{r}) Y_{l m} (\vec{r}, \vec{s}),

where

ψ_{l m} (\vec{r})

are the coefficients corresponding to the real spherical harmonics

Y_{l m} (\vec{r}, \vec{s})

, expressed as

Eq. (40)

Y_{l m} = {\begin{cases} \sqrt{2} \sqrt{\frac{2 l + 1}{4 π}} \sqrt{\frac{(l - | m |)!}{(l + | m |)!}} P_{l}^{| m |} \cos (θ) \sin (| m | ϕ) & if m < 0, \\ \sqrt{\frac{2 l + 1}{4 π}} P_{l}^{m} [\cos (θ)] & if m = 0, \\ \sqrt{2} \sqrt{\frac{2 l + 1}{4 π}} \sqrt{\frac{(l - m)!}{(l + m)!}} P_{l}^{m} \cos (θ) \cos (m ϕ) & if m > 0, \end{cases}

where

l

is the degree of the spherical harmonic,

m

is the order, and

P_{l}^{m}

is the associated Legendre polynomials. The coefficient

ψ_{l m} (\vec{r})

scales the total weight deposited by all simulated photons at position (voxel)

\vec{r}

for the associated spherical harmonic.

The advantage of expressing $ϕ^{*} ϕ$ in a spherical harmonics expansion lies in the fact that the $Y_{l m}$ forms an orthonormal basis, i.e.,

Eq. (41)

\int_{S^{2}} Y_{l m} \cdot Y_{l^{'} m^{'}} d \vec{s} = δ_{l l^{'}} δ_{m m^{'}} .

Using this orthonormality condition greatly simplifies the term ${[\int_{S^{2}} ϕ^{*} (\vec{s}) ϕ (\vec{s}) d \vec{s}]}_{\vec{r} = {\vec{r}}_{i}}$ from Eq. (38), when the radiance is expressed in spherical harmonics

Eq. (42)

{[\int_{S^{2}} ϕ^{*} (\vec{s}) ϕ (\vec{s}) d \vec{s}]}_{\vec{r} = {\vec{r}}_{i}} = \int_{S^{2}} [\sum_{l = 0}^{N_{L}} \sum_{m = - l}^{l} ψ_{l m}^{*} ({\vec{r}}_{i}) Y_{l m}^{*} ({\vec{r}}_{i}, \vec{s})] [\sum_{l^{'} = 0}^{N_{L}} \sum_{m^{'} = - l^{'}}^{l^{'}} ψ_{l^{'} m^{'}} ({\vec{r}}_{i}) Y_{l^{'} m^{'}} ({\vec{r}}_{i}, \vec{s})] d \vec{s},

Eq. (43)

= [\sum_{l^{'} = 0}^{N_{L}} \sum_{m^{'} = - l^{'}}^{l^{'}} ψ_{l^{'} m^{'}} ({\vec{r}}_{i})] [\binom{\underset{⏟}{\int_{S^{2}} Y_{l m}^{*} ({\vec{r}}_{i}, \vec{s}) Y_{l^{'} m^{'}} ({\vec{r}}_{i}, \vec{s}) d \vec{s}}}{δ_{l l^{'}} δ_{m m^{'}}}]

Eq. (44)

= \sum_{l = 0}^{N_{L}} \sum_{m = - l}^{l} ψ_{l m} ({\vec{r}}_{i}) ψ_{l m}^{*} ({\vec{r}}_{i}) .

Hence, with this approximation the gradient to update the distribution of chromophore $k$ finally becomes

Eq. (45)

\frac{\partial ε_{λ_{l}}}{\partial c_{i}} = V_{vox} {- [p_{0}^{m} ({\vec{r}}_{i}) - p_{0} ({\vec{r}}_{i})] [β_{k} H ({\vec{r}}_{i}) + Γ ({\vec{r}}_{i}) α_{k} ϕ ({\vec{r}}_{i})] + α_{k} \sum_{l = 0}^{N_{L}} \sum_{m = - l}^{l} ψ_{l m} ({\vec{r}}_{i}) ψ_{l m}^{*} ({\vec{r}}_{i})} .

9. Appendix D: Discretization

To solve the given equations numerically, the bases in which the data and model output are represented must be defined. Assuming a sampling of continuous fields in a point-wise basis $Ψ_{j} (\vec{r}) = δ (\vec{r} - {\vec{r}}_{i})$ , as shown in Ref. 18. Hence, the data projected onto this basis becomes a vector of coefficients ${\vec{p}}_{0}^{m}$

Eq. (46)

h_{j} = ⟨ Ψ_{j}, p_{0}^{m} (\vec{r}) ⟩ = p_{0}^{m} ({\vec{r}}_{j}) .

Equally, all other continuous fields are discretized. Transforming an integral of any integrable function $f (\vec{r})$ over the continuous domain $Ω$ into discretized space introduces a volume element $d Ω = V_{vox}$

Eq. (47)

\int_{Ω} f (\vec{r}) d Ω = \sum_{j}^{N_{vox}} ⟨ Ψ_{j}, f (\vec{r}) ⟩ V_{vox} = \sum_{j}^{N_{vox}} f ({\vec{r}}_{j}) V_{vox} .

Disclosures

The authors declare that there are no conflicts of interest related to this article.

Acknowledgments

The authors would like to thank Simon Arridge and Ben Cox for valuable discussions. This project was funded by the Deutsche Forschungsgemeinschaft (DFG project grants LA3273/1-1, PR1226/5-1), the Royal Academy of Engineering (RAEng Fellowship RF1516\15\33), and the Engineering and Physical Sciences Research Council (EPSRC grant EP/N032055/1). We acknowledge financial support within the funding programme Open Access Publishing by the German Research Foundation (DFG).

References

1.

P. Beard, “Biomedical photoacoustic imaging,” Interface Focus, 1 (4), 602 –631 (2011). https://doi.org/10.1098/rsfs.2011.0028 Google Scholar

2.

J. Xia, J. Yao and L. Wang, “Photoacoustic tomography: principles and advances,” Electromagn. Waves, 147 1 –22 (2014). https://doi.org/10.2528/PIER14032303 Google Scholar

3.

L. Wang and J. Yao, “A practical guide to photoacoustic tomography in the life sciences,” Nat. Methods, 13 627 –638 (2016). https://doi.org/10.1038/nmeth.3925 1548-7091 Google Scholar

4.

B. Cox, J. Laufer and P. Beard, “The challenges for quantitative photoacoustic imaging,” Proc. SPIE, 7177 717713 (2009). https://doi.org/10.1117/12.806788 Google Scholar

5.

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

6.

C. Lutzweiler and D. Razansky, “Optoacoustic imaging and tomography: reconstruction approaches and outstanding challenges in image performance and quantification,” Sensors, 13 (6), 7345 –7384 (2013). https://doi.org/10.3390/s130607345 SNSRES 0746-9462 Google Scholar

7.

M. Li, Y. Tang and J. Yao, “Photoacoustic tomography of blood oxygenation: a mini review,” Photoacoustics, 10 65 –73 (2018). https://doi.org/10.1016/j.pacs.2018.05.001 Google Scholar

8.

X. Wang et al., “Noninvasive imaging of hemoglobin concentration and oxygenation in the rat brain using high-resolution photoacoustic tomography,” J. Biomed. Opt., 11 (2), 024015 (2006). https://doi.org/10.1117/1.2192804 JBOPFO 1083-3668 Google Scholar

9.

D. Razansky et al., “Multispectral opto-acoustic tomography of deep-seated fluorescent proteins in vivo,” Nat. Photonics, 3 412 –417 (2009). https://doi.org/10.1038/nphoton.2009.98 NPAHBY 1749-4885 Google Scholar

10.

S. Tzoumas et al., “Eigenspectra optoacoustic tomography achieves quantitative blood oxygenation imaging deep in tissues,” Nat. Commun., 7 12121 (2016). https://doi.org/10.1038/ncomms12121 NCAOBW 2041-1723 Google Scholar

11.

T. Kirchner, J. Gröhl and L. Maier-Hein, “Context encoding enables machine learning-based quantitative photoacoustics,” J. Biomed. Opt., 23 (5), 056008 (2018). https://doi.org/10.1117/1.JBO.23.5.056008 JBOPFO 1083-3668 Google Scholar

12.

G. Bal and K. Ren, “On multi-spectral quantitative photoacoustic tomography in diffusive regime,” Inverse Prob., 28 (2), 025010 (2012). https://doi.org/10.1088/0266-5611/28/2/025010 INPEEY 0266-5611 Google Scholar

13.

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

14.

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

15.

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

16.

B. A. Kaplan et al., “Monte-Carlo-based inversion scheme for 3D quantitative photoacoustic tomography,” Proc. SPIE, 10064 100645J (2017). https://doi.org/10.1117/12.2251945 PSISDG 0277-786X Google Scholar

17.

B. Cox, S. Arridge and P. Beard, “Gradient-based quantitative photoacoustic image reconstruction for molecular imaging,” Proc. SPIE, 6437 64371T (2007). https://doi.org/10.1117/12.700031 PSISDG 0277-786X Google Scholar

18.

E. Malone et al., “Reconstruction-classification method for quantitative photoacoustic tomography,” J. Biomed. Opt., 20 (12), 126004 (2015). https://doi.org/10.1117/1.JBO.20.12.126004 JBOPFO 1083-3668 Google Scholar

19.

R. Hochuli et al., “Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance,” J. Biomed. Opt., 21 (12), 126004 (2016). https://doi.org/10.1117/1.JBO.21.12.126004 JBOPFO 1083-3668 Google Scholar

20.

M. Fonseca et al., “Three-dimensional photoacoustic imaging and inversion for accurate quantification of chromophore distributions,” Proc. SPIE, 10064 1006415 (2017). https://doi.org/10.1117/12.2250964 PSISDG 0277-786X Google Scholar

21.

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

22.

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Prob., 15 (2), R41 (1999). https://doi.org/10.1088/0266-5611/15/2/022 INPEEY 0266-5611 Google Scholar

23.

T. Tarvainen et al., “Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography,” Inverse Prob., 28 (8), 084009 (2012). https://doi.org/10.1088/0266-5611/28/8/084009 INPEEY 0266-5611 Google Scholar

24.

A. V. Mamonov and K. Ren, “Quantitative photoacoustic imaging in the radiative transport regime,” Commun. Math. Sci., 12 (2), 201 –234 (2014). https://doi.org/10.4310/CMS.2014.v12.n2.a1 1539-6746 Google Scholar

25.

S. Rabanser, L. Neumann and M. Haltmeier, “Stochastic proximal gradient algorithms for multi-source quantitative photoacoustic tomography,” Entropy, 20 (2), 121 (2018). https://doi.org/10.3390/e20020121 ENTRFG 1099-4300 Google Scholar

26.

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

27.

T. Soonthornsaratoon, “Gradient-based methods for quantitative photoacoustic tomography,” University College London (UCL), (2014). Google Scholar

28.

S. T. Flock et al., “Monte Carlo modeling of light propagation in highly scattering tissues. I. Model predictions and comparison with diffusion theory,” IEEE Trans. Biomed. Eng., 36 (12), 1162 –1168 (1989). https://doi.org/10.1109/TBME.1989.1173624 IEBEAX 0018-9294 Google Scholar

29.

L. Wang and S. L. Jacques, “Hybrid model of Monte Carlo simulation and diffusion theory for light reflectance by turbid media,” J. Opt. Soc. Am. A, 10 1746 –1752 (1993). https://doi.org/10.1364/JOSAA.10.001746 JOAOD6 0740-3232 Google Scholar

30.

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

31.

S. Powell, R. Hochuli and S. R. Arridge, “Radiance Monte-Carlo for application of the radiative transport equation in the inverse problem of diffuse optical tomography,” Proc. SPIE, 10059 100590W (2017). https://doi.org/10.1117/12.2250005 PSISDG 0277-786X Google Scholar

32.

J. Buchmann et al., “Experimental validation of a Monte-Carlo-based inversion scheme for 3D quantitative photoacoustic tomography,” Proc. SPIE, 10064 1006416 (2017). https://doi.org/10.1117/12.2252359 PSISDG 0277-786X Google Scholar

33.

D. Kingma and J. Ba, “Adam: a method for stochastic optimization,” (2014). Google Scholar

34.

M. B. Fonseca, L. An and B. T. Cox, “Sulfates as chromophores for multiwavelength photoacoustic imaging phantoms,” J. Biomed. Opt., 22 (12), 125007 (2017). https://doi.org/10.1117/1.JBO.22.12.125007 JBOPFO 1083-3668 Google Scholar

35.

D.-K. Yao et al., “Photoacoustic measurement of the Grüneisen parameter of tissue,” J. Biomed. Opt., 19 (1), 017007 (2014). https://doi.org/10.1117/1.JBO.19.1.017007 Google Scholar

36.

S. A. Prahl, “A monte carlo model of light propagation in tissue,” Proc. SPIE, 10305 1030509 (1989). https://doi.org/10.1117/12.2283590 PSISDG 0277-786X Google Scholar

37.

L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J., 93 70 –83 (1941). https://doi.org/10.1086/144246 ASJOAB 0004-637X Google Scholar

38.

J. Bezanson et al., “Julia: a fresh approach to numerical computing,” SIAM Rev., 59 (1), 65 –98 (2017). https://doi.org/10.1137/141000671 SIREAD 0036-1445 Google Scholar

39.

J. Laufer et al., “In vivo preclinical photoacoustic imaging of tumor vasculature development and therapy,” J. Biomed. Opt., 17 (5), 056016 (2012). https://doi.org/10.1117/1.JBO.17.5.056016 JBOPFO 1083-3668 Google Scholar

40.

J. Märk et al., “Dual-wavelength 3D photoacoustic imaging of mammalian cells using a photoswitchable phytochrome reporter protein,” Commun. Phys., 1 (1), 3 (2018). https://doi.org/10.1038/s42005-017-0003-2 Google Scholar

41.

E. Zhang, J. Laufer and P. Beard, “Backward-mode multiwavelength photoacoustic scanner using a planar Fabry-Pérot polymer film ultrasound sensor for high-resolution three-dimensional imaging of biological tissues,” Appl. Opt., 47 (4), 561 –577 (2008). https://doi.org/10.1364/AO.47.000561 APOPAI 0003-6935 Google Scholar

42.

J. Buchmann et al., “Characterization and modeling of Fabry–Perot ultrasound sensors with hard dielectric mirrors for photoacoustic imaging,” Appl. Opt., 56 (17), 5039 –5046 (2017). https://doi.org/10.1364/AO.56.005039 APOPAI 0003-6935 Google Scholar

43.

B. G. Barratt-Boyes and E. H. Wood, “The oxygen saturation of blood in the venae cavae, right-heart chambers, and pulmonary vessels of healthy subjects,” J. Lab. Clin. Med., 50 (1), 93 –106 (1957). JLCMAK 0022-2143 Google Scholar

44.

E. Beutler and J. Waalen, “The definition of anemia: what is the lower limit of normal of the blood hemoglobin concentration?,” Blood, 107 (5), 1747 –1750 (2006). https://doi.org/10.1182/blood-2005-07-3046 BLOOAW 0006-4971 Google Scholar

45.

G. Serⓢa et al., “Contrast enhanced MRI assessment of tumor blood volume after application of electric pulses,” Electro- Magnetobiol., 17 (2), 299 –306 (1998). https://doi.org/10.3109/15368379809022574 Google Scholar

46.

S. Prahl, “Assorted Spectra,” (2017) https://omlc.org/spectra/index.html Google Scholar

47.

B. T. Polyak, “Some methods of speeding up the convergence of iteration methods,” USSR Comput. Math. Math. Phys., 4 (5), 1 –17 (1964). https://doi.org/10.1016/0041-5553(64)90137-5 CMMPA9 0965-5425 Google Scholar

48.

G. Goh, “Why momentum really works,” Distill, (2017). https://doi.org/10.23915/distill.00006 Google Scholar

49.

G. P. Luke, S. Y. Nam and S. Y. Emelianov, “Optical wavelength selection for improved spectroscopic photoacoustic imaging,” Photoacoustics, 1 (2), 36 –42 (2013). https://doi.org/10.1016/j.pacs.2013.08.001 Google Scholar

50.

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

51.

T. Tarvainen, “Computational methods for light transport in optical tomography,” Faculty of Natural and Environmental Sciences, Department of Physics, University of Kuopio, (2006). Google Scholar

52.

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

53.

F. Inanc and A. Rohach, “A nodal solution of the multigroup neutron transport equation using spherical harmonics,” Ann. Nucl. Energy, 15 (10–11), 501 –509 (1988). https://doi.org/10.1016/0306-4549(88)90066-7 ANENDJ 0306-4549 Google Scholar

54.

A. G. Buchan et al., “A POD reduced order model for resolving angular direction in neutron/photon transport problems,” J. Comput. Phys., 296 138 –157 (2015). https://doi.org/10.1016/j.jcp.2015.04.043 JCTPAH 0021-9991 Google Scholar

Biographies of the authors are not available.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.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 Download Citation

Jens Buchmann, Bernhard A. Kaplan, Samuel Powell, Steffen Prohaska, and Jan Laufer "Three-dimensional quantitative photoacoustic tomography using an adjoint radiance Monte Carlo model and gradient descent," Journal of Biomedical Optics 24(6), 066001 (6 June 2019). https://doi.org/10.1117/1.JBO.24.6.066001

Received: 9 January 2019; Accepted: 24 April 2019; Published: 6 June 2019

Access the abstract

JOURNAL ARTICLE
13 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY