Quantitative photoacoustic tomography: modeling and inverse problems

Tanja Tarvainen; Ben Cox

doi:10.1117/1.JBO.29.S1.S11509

20 December 2023 Quantitative photoacoustic tomography: modeling and inverse problems

Tanja Tarvainen, Ben Cox

Author Affiliations +

Journal of Biomedical Optics, Vol. 29, Issue S1, S11509 (December 2023). https://doi.org/10.1117/1.JBO.29.S1.S11509

Abstract

Significance

Quantitative photoacoustic tomography (QPAT) exploits the photoacoustic effect with the aim of estimating images of clinically relevant quantities related to the tissue’s optical absorption. The technique has two aspects: an acoustic part, where the initial acoustic pressure distribution is estimated from measured photoacoustic time-series, and an optical part, where the distributions of the optical parameters are estimated from the initial pressure.

Aim

Our study is focused on the optical part. In particular, computational modeling of light propagation (forward problem) and numerical solution methodologies of the image reconstruction (inverse problem) are discussed.

Approach

The commonly used mathematical models of how light and sound propagate in biological tissue are reviewed. A short overview of how the acoustic inverse problem is usually treated is given. The optical inverse problem and methods for its solution are reviewed. In addition, some limitations of real-life measurements and their effect on the inverse problems are discussed.

Results

An overview of QPAT with a focus on the optical part was given. Computational modeling and inverse problems of QPAT were addressed, and some key challenges were discussed. Furthermore, the developments for tackling these problems were reviewed. Although modeling of light transport is well-understood and there is a well-developed framework of inverse mathematics for approaching the inverse problem of QPAT, there are still challenges in taking these methodologies to practice.

Conclusions

Modeling and inverse problems of QPAT together were discussed. The scope was limited to the optical part, and the acoustic aspects were discussed only to the extent that they relate to the optical aspect.

1. Introduction

As the name photoacoustic tomography suggests, there are two aspects to this emerging imaging modality: an optical part and an acoustic part. This short review paper is focused on the mathematics of the optical part. In particular, it surveys the current thinking regarding two related problems: (1) what is the best way to describe light propagation and its interaction with biological tissue mathematically? (2) Given photoacoustic measurements, what, in principle, can we learn about the optical properties of tissue (or indeed the related, and more clinically relevant, properties, such as blood oxygenation)? Because the ultimate aim is to obtain quantitative estimates of the tissue constituents, this topic is sometimes referred to as quantitative photoacoustic tomography (QPAT).

1.1.

Photoacoustic Imaging

There are a few closely related but different imaging modalities that come under the heading of photoacoustic imaging. All exploit the photoacoustic effect, which is when a sufficiently short pulse of light is absorbed by an elastic material and subsequently thermalized, the site of the absorption will act as a source of an acoustic pulse.¹^–³ In all variants, the light pulse is directed into the soft biological tissue under investigation, and the resulting acoustic pulse is measured at the tissue surface. From the measurements of the acoustic pulse, an image of where the light was absorbed can be formed. That is a photoacoustic image. Photoacoustic microscopy differs from photoacoustic tomography in the way the data is collected and the image is formed. In photoacoustic microscopy, either the light beam or the acoustic detector is sharply focused and raster-scanned across the tissue surface.¹^,⁴ Because of the localization caused by the focusing, an image can be formed directly from the measured acoustic time series; indeed, it is the tightness of the focusing that determines the resolution of the image. (The fact that the source or detector is commonly raster-scanned is not what makes this microscopy; an array of focused sources or detectors could just as well be used.) In photoacoustic tomography, by contrast, the light is unfocused—indeed, the illumination is arranged such that the whole region-of-interest is flooded with light—and an array of unfocused (or, at least, not tightly focused) detectors is used to record the resulting acoustic time series.¹^,² Because the photoacoustic source may be distributed throughout the tissue and because each time series could contain signals from anywhere (as the detectors are unfocused), the connection between the data and source is more complicated than for microscopy, and it is necessary to use an image reconstruction algorithm to form an image. Photoacoustic tomography, not microscopy, is the primary concern of this review, although the tissue optics described will be applicable to all photoacoustic imaging approaches in turbid media.

1.2.

Scope of the Review

Even with this restriction to photoacoustic tomography, the size of the field has grown rapidly in recent years and it is not possible to review all aspects of it in a short article. We have therefore limited the scope further and discussed the acoustic aspects of photoacoustic tomography only to the extent that they relate to the optical aspects, which is our primary concern. Furthermore, we do not describe experimental photoacoustic tomography, except to note how practical constraints impact the optical inversions (often very significantly) and make few references to specific applications. Finally, we do not provide a comprehensive literature review as that would comprise too long a list, but we hope that the articles we do reference can act as a route into the broader field for the interested reader.

1.3.

Layout of the Paper

In Sec. 2, the commonly used mathematical models of how light (Sec. 2.1) and sound (Secs. 2.2 and 2.3) propagate in biological tissue are introduced, and in Sec. 3, the corresponding inverse problems are introduced. In particular, Sec. 3.1 briefly overviews how the acoustic inverse problem is usually treated, leading up to the main topic of this review in Sec. 3.2: the optical inverse problem. The review ends with Sec. 4, which highlights some of the ways in which limitations of real-life measurements affects the inverse problems in practice and the trade-offs necessary when using experimental measurements.

2. Modeling Photoacoustic Waves

The key physical phenomena and quantities relevant to photoacoustic tomography are shown schematically in Fig. 1. The principal physical phenomenon of relevance to QPAT is light transport; the other physical phenomena are relevant to the inverse problem of obtaining an accurate estimate for the initial acoustic pressure distribution, which acts as the data for the optical part of the problem.

Fig. 1

Key stages in QPAT and the related physical phenomena and quantities.

2.1.

Light Transport

Light propagation in scattering media, such as biological tissue, can be described using transport theory.⁵^,⁶ In transport theory, light energy conservation within a small volume element of phase space is investigated. Wave phenomena, such as interference, Anderson localization, and enhanced backscattering, are assumed to be negligible and are ignored. Light transport can be modeled using deterministic and stochastic methods. In the deterministic approach, light transport is described with integro-differential equations that can be solved analytically or the solution can be numerically approximated. Generally in QPAT, the radiative transfer equation (RTE) or its approximations are used. The RTE is a “one-speed” approximation of the transport equation, which means it is assumed that the energy (or speed) of photons does not change in collisions (elastic collisions only) and that the refractive index is constant within the medium. For a discussion and studies of photon transport in a medium with spatially varying or piecewise constant refractive index, see Refs. 7 8.9.–10 and the references therein. In the stochastic approach, the photons’ absorption and scattering interactions with the medium are simulated directly. The most widely applied stochastic methods for simulating light transport in biological tissue are Monte Carlo methods, briefly described in Sec. 2.1.5 (see also reviews Refs. 11 and 12 and the references therein).

2.1.1.

Optical absorption and scattering

A photon of light is a rapid fluctuation in the electromagnetic field, which when close to a molecule can induce a dipole moment in the molecule, self-creating a mechanism through which the molecule then interacts with the photon. When the frequency of the fluctuation corresponds to the energy of an allowed transition between two energetic states of the molecule, energy transfer from the light to the molecule—absorption of the photon—is highly likely to occur. Some of the absorbed energy may then be reradiated, e.g., as fluorescence, but the part of interest in photoacoustics is the part converted to heat, typically via collisional relaxation with neighboring molecules, often water molecules. A measure of the likelihood of absorption occurring in a bulk material made up of many of the same molecules is given by the “molar absorption coefficient,” but most biological tissue consists of many different types of molecules, so it is useful to define the overall “absorption coefficient” $μ_{a} (m^{- 1}, often stated in {mm}^{- 1})$ of the tissue as a linear sum of the constituent components:¹³^,¹⁴

Eq. (1)

μ_{a} (λ) = \sum_{k = 1}^{K} α_{k} (λ) C_{k} (r),

where there are

K

components,

α_{k} (λ)

and

C_{k} (r)

are the molar absorption coefficient and concentration of the

k

’th component, respectively, and

λ

is the wavelength of the light. The concentrations

C_{k} (r)

, and therefore the absorption coefficient, will in general be spatially varying, where

r

denotes the spatial position. The absorption coefficient describes the likelihood of a photon being absorbed; specifically, the probability of a photon being absorbed while travelling a short distance

d s

is

μ_{a} d s

. Indeed, the absorption coefficient is also the rate in which the intensity of light beam will decay in a purely absorbing (nonscattering) medium, i.e., the intensity will decay as

\exp (- μ_{a} s)

, where

s

is the propagation distance.

When the optical frequency is not close to the energy of an allowed transition, an oscillating dipole moment can be induced that reradiates the wave, i.e., scatters the photon in a random direction. In optical transport theory, the probability density function describing the random scattering of incident light from a bulk medium is called the “scattering phase function,” denoted here by $Θ (\hat{s}, {\hat{s}}^{'})$ , where the unit vectors $\hat{s}$ and ${\hat{s}}^{'}$ are the scattered and incident photon directions, respectively. The phase function will depend on the optical wavelength as well as the molecules doing the scattering. The probability that a photon will be scattered while traversing the distance $d s$ is $μ_{s} d s$ , where $μ_{s} (m^{- 1})$ is called the “scattering coefficient.” It is through these three quantities, the absorption and scattering coefficients and the scattering phase function, all wavelengths dependent that the optical properties of tissue are described in transport theory.

2.1.2.

Radiative transfer

The RTE can be derived through transport theory⁵ or from Maxwell’s equations.¹⁵^,¹⁶ It describes the distribution of radiance within a domain $Ω \subset R^{d}$ with boundary $\partial Ω$ as

Eq. (2)

{\begin{cases} \frac{1}{c} \frac{\partial ϕ (r, \hat{s}, t, λ)}{\partial t} + \hat{s} \cdot \nabla ϕ (r, \hat{s}, t, λ) + (μ_{s} (r, λ) + μ_{a} (r, λ)) ϕ (r, \hat{s}, t, λ) \\ = μ_{s} (r, λ) \int_{S^{d - 1}} Θ (\hat{s}, {\hat{s}}^{'}) ϕ (r, {\hat{s}}^{'}, t, λ) d {\hat{s}}^{'} + q (r, t, λ), r \in Ω \\ ϕ (r, \hat{s}, t, λ) = {\begin{cases} ϕ_{0} (r, \hat{s}, t, λ) + R ϕ (r, M \hat{s}, t, λ), & r \in r_{s}, \hat{s} \cdot \hat{n} < 0 \\ R ϕ (r, M \hat{s}, t, λ), & r \in \partial Ω ∖ r_{s}, \hat{s} \cdot \hat{n} < 0 \end{cases} \end{cases},

where

d

is the (spatial) dimension of the domain (

d = 2

or 3),

\hat{s} \in S^{d - 1}

is a unit vector in the direction of interest,

ϕ (r, \hat{s}, t, λ) [W / (m^{2} sr)]

is the radiance at point

r

, direction

\hat{s}

, and time instance

t

,

c

is the speed of light in the medium,

q (r, t, λ)

is an internal light source, which can also be placed on the boundary,

ϕ_{0} (r, \hat{s}, t, λ)

is a boundary light source at the source position

r_{s} \subset \partial Ω

, and

\hat{n}

is an outward unit normal (see Fig. 2).

Fig. 2

Illustration of the terms in the RTE considered for one direction of propagation $\hat{s}$ . The rate of change of the radiance $ϕ (\hat{s})$ within the control volume will depend on the net amount of light travelling in direction $\hat{s}$ entering and leaving the volume ( $\hat{s} \cdot \nabla ϕ$ ), the photon generated or absorbed within the control volume [ $q$ and $μ_{a} ϕ (\hat{s})$ , respectively], the photons scattered out of direction $\hat{s}$ into another direction [ $μ_{s} ϕ (\hat{s})$ ], and finally the photons scattered from any other direction into direction $\hat{s}$ (given by the integral term).

Further, $R$ is the Fresnel reflection coefficient and the mapping $M$ gives the change in light direction due to reflection at the boundary,⁸ thus in the case of matched refractive indices between the medium and surrounding medium, $R = 0$ . The radiance can be defined such that the amount of power transfer in the infinitesimal angle $d \hat{s}$ in direction $\hat{s}$ at time $t$ through an infinitesimal area $d S$ is given by

ϕ (r, \hat{s}, t, λ) \hat{s} \cdot \hat{n} d S d \hat{s},

where

\hat{n}

is the unit normal to the surface

d S

.⁵ The scattering phase function

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

describes the probability that a photon with an initial direction

{\hat{s}}^{'}

will have a direction

\hat{s}

after a scattering event. It is often assumed that this depends only on the angle

θ

between the incoming and outcoming directions, i.e.,

Θ (\hat{s}, {\hat{s}}^{'}) \approx Θ (\hat{s} \cdot {\hat{s}}^{'})

where

\hat{s} \cdot {\hat{s}}^{'} = \cos θ

. In biomedical optical imaging, the most commonly applied phase function is the Henyey–Greenstein phase function,¹⁷ which is of the form:

Eq. (3)

Θ_{HG} (\hat{s} \cdot {\hat{s}}^{'}) = {\begin{matrix} \frac{1}{2 π} \frac{1 - g^{2}}{(1 + g^{2} - 2 g (\hat{s} \cdot {\hat{s}}^{'}))}, & d = 2 \\ \frac{1}{4 π} \frac{1 - g^{2}}{{(1 + g^{2} - 2 g (\hat{s} \cdot {\hat{s}}^{'}))}^{3 / 2}}, & d = 3 \end{matrix},

where

g (dimension less)

is the scattering anisotropy parameter that defines the shape of the probability density. It takes values between

- 1 < g < 1

, such that, if

g = 0

, the scattering probability density is a uniform distribution,

g > 0

for forward dominated scattering, and

g < 0

for backward dominated scattering. Experimental measurements made for determining the anisotropy parameter

g

show that it is typically highly forward scattering in biological tissues.¹⁴ The Henyey–Greenstein phase function has the useful property that its Legendre expansion is given by powers of the parameter

g

:

Eq. (4)

\int_{S^{d - 1}} P_{l} (τ) Θ_{HG} (τ) d {\hat{s}}^{'} = g^{l},

and we have an expression:

Eq. (5)

Θ_{HG} (\hat{s} \cdot {\hat{s}}^{'}) = {\begin{matrix} \frac{1}{2 π} + \frac{1}{π} \sum_{n = 1}^{\infty} g^{n} \cos (n (\hat{s} \cdot {\hat{s}}^{'})), & d = 2 \\ \sum_{l}^{\infty} \sum_{m = - l}^{l} g^{l} {\bar{Y}}_{l, m} ({\hat{s}}^{'}) Y_{l, m} (\hat{s}), & d = 3 \end{matrix},

where

d = 2

is a Fourier series representation and

d = 3

is spherical harmonics

Y_{l, m}

expansion of the phase function.¹⁸^,¹⁹ Although widely applied in biomedical optics, the Henyey–Greenstein phase function may show differences to measured signals in some situations, such as spatially resolved reflectance measurements, where other phase functions could be applied.²⁰

As light propagates within the tissue, it is absorbed leading to a localized increase in pressure and generation of a pressure wave. In QPAT, propagation of the acoustic wave occurs on a microsecond time scale, orders of magnitude slower than the optical timescale that includes the optical pulse length as well as the optical propagation, absorption, and the decay of the absorbed optical energy to heat. Therefore, only the total absorbed optical energy density is of interest and not the rate of the absorption. Thus, in QPAT, light propagation can be modeled using a time-independent model for light transport. The time-independent RTE is

Eq. (6)

{\begin{cases} \hat{s} \cdot \nabla ϕ (r, \hat{s}, λ) + (μ_{s} (r, λ) + μ_{a} (r, λ)) ϕ (r, \hat{s}, λ) \\ = μ_{s} (r, λ) \int_{S^{d - 1}} Θ (\hat{s} \cdot {\hat{s}}^{'}) ϕ (r, {\hat{s}}^{'}) d {\hat{s}}^{'} + q (r, λ), r \in Ω \\ ϕ (r, \hat{s}, λ) = {\begin{cases} ϕ_{0} (r, \hat{s}, λ) + R ϕ (r, M \hat{s}, λ), & r \in r_{s}, \hat{s} \cdot \hat{n} < 0 \\ R ϕ (r, M \hat{s}, λ), & r \in \partial Ω ∖ r_{s}, \hat{s} \cdot \hat{n} < 0 \end{cases} \end{cases},

where

ϕ (r, \hat{s}, λ) = \int_{- \infty}^{\infty} ϕ (r, \hat{s}, t, λ) d t [J / (m^{2} sr)]

is the time-independent radiance, and

q (r, λ)

and

ϕ_{0} (r, \hat{s}, λ)

are time-independent light sources.

2.1.3.

Approximations to the radiative transfer equation

In biomedical optics, popular approximations of the RTE are the diffusion approximation (DA), which is valid when the field is diffuse, and the Beer–Lambert law, which is valid in nonscattering media.

Diffusion approximation

The typical approach to derive the DA is to expand the radiance, source term, and phase function into series using the spherical harmonics and truncate the series.¹⁸ The first-order spherical harmonics approximation is referred to as the $P_{1}$ approximation and the DA can be regarded as a special case of that. In the framework of the DA, the radiance is approximated by

Eq. (7)

ϕ (r, \hat{s}, t, λ) \approx \frac{1}{| S^{n - 1} |} Φ (r, t, λ) + \frac{n}{| S^{n - 1} |} \hat{s} \cdot J (r, t, λ),

where

Φ (r, t, λ) (W / m^{2})

and

J (r, t, λ) (W / m^{2})

are the photon fluence rate and photon current (flux):

Eq. (8)

Φ (r, t, λ) = \int_{S^{d - 1}} ϕ (r, \hat{s}, t, λ) d \hat{s},

Eq. (9)

J (r, t, λ) = \int_{S^{d - 1}} \hat{s} ϕ (r, \hat{s}, t, λ) d \hat{s} .

By inserting the approximation (7) and similar approximations written for the source term and phase function into Eq. (2) and following the derivation in Refs. 5 and 21, the DA can be derived

Eq. (10)

{\begin{cases} \frac{1}{c} \frac{\partial Φ (r, t, λ)}{\partial t} - \nabla \cdot κ (r, λ) \nabla Φ (r, t, λ) + μ_{a} (r, λ) Φ (r, t, λ) = q_{0} (r, λ), r \in Ω \\ Φ (r, t, λ) + \frac{1}{2 γ_{d}} κ (r, λ) A \frac{\partial Φ (r, t, λ)}{\partial \hat{n}} = {\begin{cases} \frac{I_{s} (r, t, λ)}{γ_{d}}, & r \in r_{s} \\ 0, & r \in \partial Ω ∖ r_{s} \end{cases} \end{cases},

where

κ (r, λ) = {(d (μ_{a} (r, λ) + μ_{s}^{'} (r, λ)))}^{- 1} (m^{- 1})

is the diffusion coefficient and where

μ_{s}^{'} (r, λ) = (1 - g_{1}) μ_{s} (r, λ) (m^{- 1})

is the reduced scattering coefficient,

d

is the dimension (

d = 2

or 3), and

g_{1}

is the mean of the cosine of the scattering angle that in the case of the Henyey–Greenstein phase function (3) is

g_{1} = g

. It is clear from this that, in highly scattering media, the scattering anisotropy parameter and scattering coefficient merge into a single parameter, the reduced scattering coefficient. Furthermore,

γ_{d}

is a dimension-dependent constant that takes values

γ_{2} = 1 / π

and

γ_{3} = 1 / 4

. Parameter

A = (1 + R) (1 - R)^{- 1}

governs reflection on the boundary, with

A = 1

in the case of no boundary reflection.

q_{0} (r, λ)

is an internal light source that is utilized in biomedical optics in fluorescence diffuse optical tomography (DOT) and in numerical simulations of light transport if an internal point source model is used, and

I_{s} (r, t, λ)

is a diffuse boundary source. It should be noted that, in biomedical optics, the boundary condition of the DA is sometimes numerically implemented by applying a Diriclet boundary condition on a “virtual boundary” outside the domain (so-called extrapolated boundary condition).²²

As in the case of the RTE, the time-independent DA can be derived. It takes the form:

Eq. (11)

{\begin{cases} - \nabla \cdot κ (r, λ) \nabla Φ (r, λ) + μ_{a} (r, λ) Φ (r, λ) = q_{0} (r, λ), r \in Ω \\ Φ (r, λ) + \frac{1}{2 γ_{d}} κ (r, λ) A \frac{\partial Φ (r, λ)}{\partial \hat{n}} = {\begin{cases} \frac{I_{s} (r, λ)}{γ_{d}}, & r \in r_{s} \\ 0, & r \in \partial Ω ∖ r_{s} \end{cases} \end{cases},

where

Φ (r, λ) = \int_{- \infty}^{\infty} Φ (r, t, λ) d t (J / m^{2})

is the (time-independent) fluence and

q_{0} (r, λ)

and

I_{s} (r, λ)

are time-independent light sources.

The DA is valid when the radiance is almost a uniform distribution [see approximation of radiance Eq. (7)]. In other words, when the radiance is almost independent of direction; the light reaching any point is coming from all directions. In practice, this is achieved in a scattering-dominated medium further than a few scattering lengths from light sources.⁵ In many QPAT experiments, however, the imaging depth can be small compared to the average scattering length, and thus the DA is not always a valid approximation.

Beer–Lambert law

If the medium is nonscattering (purely absorbing) and if the radiance is collimated (propagating in one direction), the RTE becomes

Eq. (12)

(\hat{s} \cdot \nabla) ϕ (r, \hat{s}, λ) = - μ_{a} (r, λ) ϕ (r, \hat{s}, λ),

from which it is clear that the decrease in intensity through a purely absorbing medium in direction

z

can be written as

Eq. (13)

d Φ (z, λ) = - μ_{a} (z, λ) Φ (z, λ) d z .

If the source is a monodirectional flux

Φ_{0} (r, λ)

incident in the positive

z

-direction, Eq. (13) has the solution:

Eq. (14)

Φ (r, λ) = Φ_{0} (r, λ) \exp (- \int_{z_{0}}^{z} μ_{a} (z^{'}, λ) d z^{'}),

which is known as the Beer–Lambert law.

A second scenario, in which exponential decay of the light flux occurs, is when a plane wave is incident on a homogeneous scattering medium for which the DA holds. In this case, Eq. (11) becomes

Eq. (15)

\frac{d^{2} Φ (z, λ)}{d z^{2}} = (\frac{μ_{a}}{κ}) Φ (z, λ),

which has solution

Eq. (16)

Φ (z, λ) = Φ_{0} (z, λ) \exp (- μ_{eff} z), μ_{eff} = \sqrt{μ_{a} / κ},

where

z

is the coordinate perpendicular to the propagation front of the plane wave. This approximation and similar approximations exploiting the parameter

μ_{eff}

have been used in a number of QPAT studies.²³^–²⁸

Other approximations

In addition to DA, other orders of the spherical harmonics expansion $P_{N}$ and approximations made for them, such as simplified spherical harmonics,²⁹ have been utilized in biomedical optics. Furthermore, hybrid models coupling nonscattering and/or low-scattering medium with highly scattering medium have been developed.³⁰^–³⁴ In QPAT, utilizing simplified spherical harmonics approximation has been utilized.³⁵^,³⁶ Numerical approximation of the RTE with Henyey–Greenstein phase function can be challenging if scattering is highly forward dominated. Therefore, approximations to the RTE that takes into account forward-peaked scattering analytically have been proposed. These include the delta-Eddington approximation, the Fokker–Planck approximation, the Fokker–Planck–Eddington approximation and the generalized Fokker–Planck–Eddington approximation.³²^,³⁷^,³⁸

2.1.4.

Numerical approximations

The analytical solutions of the RTE and the DA are typically limited to specific geometries, and therefore their utilization in biomedical optical imaging has been limited. Therefore, the typical approach has been to numerically approximate their solutions. For the numerical approximation of the DA, the most widely applied approach has been to use the finite-element method,²²^,³⁹ with free software also available.⁴⁰^,⁴¹ The numerical approaches for the RTE include finite difference,⁴²^,⁴³ finite element,³³^,⁴⁴^,⁴⁵ and finite volume⁴⁶ methods for spatial discretization. The numerical stability of these has been advanced using special basis functions, such as the streamline diffusion modification.³³^,⁴⁷ In addition, numerical approximation using the discontinuous Galerkin method has been implemented.⁴⁸ For angular discretization, discrete ordinates,⁴²^,⁴³^,⁴⁶ finite elements,³³ and spherical harmonics⁴⁴^,⁴⁵ have been used. Recently, the pseudospectral method for the numerical approximation of the RTE was proposed.⁴⁹ Recently, many computational challenges related to the numerical approximation of the RTE have been overcome through the development of computing resources and numerical methods such as parallel computing and preconditioning.⁵⁰^,⁵¹

2.1.5.

Monte Carlo method for light transport

The Monte Carlo method for light transport can be used to simulate the propagation of photons in a scattering medium in which they undergo random absorption and scattering events.⁵²^,⁵³ The methodology has been utilized in a variety of applications in biomedical optics, and different implementations of the Monte Carlo method exist¹¹^,¹²^,⁵⁴^,⁵⁵ together with various free and/or open access software available.⁵³^–⁶⁰ In addition to depositing the energy (or photon fluence), it is also possible to record the direction of an absorbed photon (packet). This is sometimes called radiance Monte Carlo.⁶¹^–⁶³

Monte Carlo simulation in a biological tissue obeys the following principles.¹¹^,⁵²^,⁵³ The scattering length follows an exponential probability distribution function:

Eq. (17)

f (l) = μ_{s} (l) \exp (- \int_{0}^{l} μ_{s} (l^{'}) d l^{'}) .

Then in the case of a scattering event, the scattering angle follows a probability distribution for scattering direction that in many studies is the Henyey–Greenstein phase function (3). For absorption, a general approach in biomedical optics has been to use a so-called photon packet method.⁵² In this approach, a photon packet with an initial weight

w_{0}

is simulated. As the photon packet propagates, its weight along trajectory is reduced due to absorption, and the photon weight can be described as

Eq. (18)

w (s) = w_{0} \exp (- \int_{0}^{s} μ_{a} (s^{'}) d s^{'}) .

where

μ_{a} (s^{'})

is the absorption coefficient along the trajectory of the photon packet. This is continued until the photon packet exits the computation domain, or its weight becomes negligible. Sampling scattering lengths from Eq. (17) with a weight factor assigned to those paths according to Eq. (18) is a form of importance sampling.⁶⁴^,⁶⁵

In QPAT, the absorbed optical energy density $H_{j}$ in a discretization element $j$ of the computation domain can be computed as

Eq. (19)

H_{j} = - \frac{1}{A_{j}} \int_{0}^{s_{t}} χ_{j} (s) \frac{d w}{d s} (s) d s,

where

A_{j}

is the area or volume of the discretization element

j

, the integral is understood as being carried from source position where the photon packet was created (

s = 0

) until where it is terminated (

s = s_{t}

),

χ_{j}

is the characteristic function of

j

’th element, and

- \frac{d w}{d s} (s)

is the energy absorbed by the medium during the photon packet propagation.⁶⁴^,⁶⁶

2.2.

Linking Light and Sound: the Photoacoustic Efficiency

Ultimately, the photons that have been multiply scattered around in the tissue either leave the tissue or are absorbed by chromophores (light-absorbing molecules). The resulting absorbed optical energy density $H (r, λ) (J / m^{3})$ can be written as

Eq. (20)

H (r, λ) = μ_{a} (r, λ) Φ (r, λ),

where

Φ (r, λ)

is the photon fluence as described by Eqs. (6), (8), and (11). Assuming there is no reradiation such as fluorescence (often a good assumption for tissue), all the energy thermalizes into heat. For efficient acoustic generation, the optical pulse must be short enough that the tissue does not have time to deform significantly (an isochoric condition). This results in a localized pressure increase, often termed the “initial acoustic pressure distribution”

p_{0} (r) (Pa)

as it is this that leads to the acoustic pressure pulse that travels to the tissue surface. The absorbed optical energy density is linked to the initial acoustic pressure distribution through the photoacoustic efficiency, which can be identified with the Grüneisen parameter

G (r) (dimensionless)

for a pure optically absorbing fluid:

Eq. (21)

p_{0} (r) = G (r) H (r, λ) = \frac{β v^{2}}{C_{p}} H (r),

where

β

(

K^{- 1}

) is the volume thermal expansivity of the fluid,

v

(m/s) is the speed of sound, and

C_{p}

(J/kg/K) is the specific heat capacity at constant pressure.¹^–³ It is worth noting that the Grüneisen parameter can depend on the wavelength,⁶⁷ although it is usually treated as constant and can also depend of the concentration of the absorber⁶⁸^–⁷⁰ and the temperature.⁷¹^,⁷²

2.3.

Acoustic Initial Value Problem

This paper, as mentioned in Sec. 1, is primarily concerned with optical aspects of photoacoustic tomography, and so this section on the acoustic propagation will be a brief summary. In photoacoustic tomography, the acoustic amplitudes are sufficiently low that the propagation of the photoacoustic wave can be modeled using equations of linear acoustics. For soft biological tissue, it is generally assumed that the medium is acoustically isotropic and quiescent. Shear waves can be neglected because they are not strongly generated by the photoacoustic effect; furthermore, the speed of shear wave propagation is orders of magnitude lower than the compressional wave speed. In a lossless medium, the linear acoustic wave equation can be written as

Eq. (22)

\frac{1}{v (r)^{2}} \frac{\partial^{2} p (r, t)}{\partial t^{2}} + \frac{1}{ρ_{0} (r)} \nabla ρ_{0} (r) \cdot \nabla p (r, t) - \nabla^{2} p (r, t) = 0,

where

p (r, t) (Pa)

is the acoustic pressure,

v (r)

is the speed of sound, and

ρ_{0} (r)

is the ambient mass density. (In the case of constant ambient density and constant speed of sound, it reduces to the canonical wave equation.) Combining the wave equation with the initial conditions relevant to photoacoustic wave generation:

Eq. (23)

p (r, t = 0) = p_{0} (r), \frac{\partial p (r, t = 0)}{\partial t} = 0,

gives an acoustic initial value problem. A spatially varying sound speed is often taken into account in photoacoustic models,⁷³^–⁷⁷ although few consider spatially varying mass density. The effect of acoustic attenuation, not modeled by the above equation, has also been considered in Refs. 78 79.80.–81. Furthermore, modeling wave propagation in elastic media, such as bones, has been considered in Refs. 82 83.84.–85. Several numerical models are freely available for simulating photoacoustic waves and for use in solving the inverse problem.⁸⁶^–⁹²

As real ultrasound detectors are, sadly, never point-like and instantaneously responsive pressure sensors, the measured data will not be the acoustic pressure but a filtered version of it (filtered in both temporal and spatial frequency space):

Eq. (24)

p_{d} (r_{d}, t) = M (p (r, t) |_{r = r_{d}}),

where

r_{d}

are the nominal positions of the detectors. The linear filter

M

can be usefully approximated in many cases as the separable product of a temporal and a spatial filter, often referred to as the frequency response and the directionality of the detector, respectively. Accounting for these filtering effects can be important when tackling the inverse problem (see Sec. 4). There is a large literature on modeling ultrasound sensor responses.⁹³^,⁹⁴

3. Photoacoustic Inverse Problems

In the inverse problem of QPAT, the aim is to estimate the spatially varying concentrations of light absorbing molecules (or related quantities) when the measured photoacoustic time series and the parameters of input light are given. As mentioned earlier, the difference in time scales of light absorption and ultrasound propagation allows the optical and acoustic parts of the problem to be decoupled. Thus the related inverse problems can be treated separately. The two inverse problems in QPAT then are: (1) acoustic inverse problem, that is, estimate the initial acoustic pressure distribution $p_{0} (r)$ from measured photoacoustic time-series and (2) optical inverse problem, that is, estimate the distributions of the optical parameters from the absorbed optical energy density (or the initial pressure).

When the two inverse problems are treated separately, and the output of the acoustic inversion becomes the data for the optical inversion, the quality of the acoustic inversion is clearly critical to the success of the optical inversion.⁹⁵

3.1.

Acoustic Inverse Problem

Before an experimental photoacoustic system is constructed, a choice must be made regarding the trade-off between the speed with which the data will be acquired and the potential accuracy of the image. (There are, of course, other considerations that affect this decision, such as cost, availability of equipment, and clinical restraints.) Without sufficient data it will never be possible to reconstruct an accurate image. However, supposing that sufficient data is available, a second and related choice must be made when trying to estimate the initial acoustic pressure distribution $p_{0} (r)$ from the measured data $p_{meas} (r_{d}, t)$ , which is regarding the trade-off of accuracy of the image and the speed with which the image is computed. There are multiple ways in which the many approaches to photoacoustic image reconstruction can be categorized, and none are perfect, but for this review, we consider three broad categories: (1) methods that are based on an analytical formula, (2) methods that rely on a numerical acoustic model, and (3) data-driven methods.

Under the first category fall all filtered backprojection-type methods,⁹⁶^–⁹⁸ which includes beamforming methods, and which relate to the inverse spherical mean Radon transform. These can be coded up to be very fast, as they require just (pre- and sometimes post-) filtering and a backprojection step, which maps the filtered data directly into an image. It also includes methods that are based on (truncated) series formulas, which can also be made fast in some cases.⁹⁹^,¹⁰⁰ However, these methods are restricted to certain measurement-surface geometries and they typically assume the speed of sound is a constant. Also incorporating a complex directional response for the detectors is nontrivial.

These problems can be overcome by methods that use numerical approximations of acoustic models. Under this category falls time-reversal methods,⁸¹^,¹⁰¹^–¹⁰⁴ regularized least squares,⁸³^,¹⁰⁵^–¹¹³ maximum entropy,¹¹⁴ and Bayesian approaches.¹¹⁵^–¹¹⁸ The models powering these approaches can be very general and can in principle incorporate the modeling of any aspect of the data acquisition process, although these may have knock-on implications for solving the inverse problem so some parsimony is recommended. Furthermore, with the latter two approaches it is possible to incorporate some types of prior knowledge of the solution to guide the inversion. Common among these are smoothness and total-variation constraints. Moreover, the Bayesian approach can even take into account the uncertainties of parameters, models and geometries.¹¹⁵^,¹¹⁹^,¹²⁰ When the number of detectors and image voxels is small, it can be possible to code the model in a matrix form which can be executed rapidly. However, to obtain accurate images for large imaging volumes with fully sampled measurement surfaces, these methods can be slow to compute as the photoacoustic wavefield within the entire domain needs to be computed at multiple iterations.

Data-driven approaches,¹²¹^–¹²³ the third category, can overcome both these limitations by learning, in advance, what kind of image is expected, leading to a fast reconstruction at run-time. (The word “data” in “data-driven” is not referring to the experimental data $p_{meas} (r_{d}, t)$ but to a training set of images or raw data, assumed to be drawn from the same distribution as the required image or data, which are available prior to the reconstruction.) In a so-called end-to-end framework, a neural network is trained to estimate the initial pressure distribution based on a set of photoacoustic data.¹²⁴^,¹²⁵ However, because this mapping is nonlocal (the acoustic waves from anywhere can reach the detectors), this will require a large interconnected network, which will require a large training set. A better approach is to incorporate learned components either before or after a backprojection step. In the preprocessing approaches, a photoacoustic dataset is first processed using a neural network to improve the signal quality and then a photoacoustic image is reconstructed using a backprojection method.¹²⁶^,¹²⁷ In the postprocessing approaches, a (rough) photoacoustic image is first reconstructed and then a neural network is used to correct artifacts or noise in the reconstructed image.¹²⁷^–¹²⁹ The trade-off to these advantages, however, is the risk that the learned components hallucinate image features, which are not consistent with the measured data. Various ways to mitigate these disadvantages, such as combining aspects of the above methods, are under investigation.¹³⁰^–¹³⁷

The short treatment of the acoustic inversion here should not be taken to suggest that it is unimportant with regard to the optical inversion. Accurate solutions to the acoustic inversion are critical to obtaining accurate solutions to the optical inversion, as they provide the data for it. With this in mind, some practical considerations that should be considered when dealing with experimental data are described below in Sec. 4. For more information on image reconstruction methods in PAT see reviews¹³⁸^–¹⁴⁰ and the references therein.

3.2.

Optical Inverse Problem

The optical inverse problem of QPAT is to estimate the spatially varying absorption coefficient, or the spatially varying chromophore concentrations, or related quantities, such as haemoglobin oxygenation state from the estimated initial pressure distribution with a given input light illumination. If the Grüneisen parameter is assumed to be known, the key task is to determine the photon fluence, $Φ (r, λ)$ [see Eq. (20)]. This, however, is complicated by the fact that fluence itself is dependent on the unknown optical properties. This inversion, unlike the acoustic inversion, is therefore nonlinear.

An example of Monte Carlo simulated data in a 2D slab of size $15 mm \times 10 mm$ with a spatially varying absorption coefficient $μ_{a} \in [0.005, 0.04] {mm}^{- 1}$ and a constant scattering $μ_{s} = 4 {mm}^{- 1}$ , and an anisotropy parameter $g = 0.9$ is shown in Fig. 3. The domain was illuminated with a planar illumination from the top of the domain with a spatially constant and angularly cosine radiance distribution. This figure shows the dependence of fluence $Φ$ on absorption distribution, and the effect of both absorption and fluence on absorbed optical energy density $H$ within the domain, demonstrating the nonlinear nature of the inverse problem of estimating the optical parameters. A similar simulation with a spatially varying scattering $μ_{s} \in [1, 8] {mm}^{- 1}$ and a constant absorption $μ_{a} = 0.02 {mm}^{- 1}$ is shown in Fig. 4. As it can be seen, also scattering affects both the fluence distribution $Φ$ and absorbed optical energy density $H$ . This change in the absorbed optical energy density distribution is, however, less clear than the change due to the absorption, which indicates that reconstruction of scattering is more ill-posed than reconstruction of absorption.

Fig. 3

Simulated photoacoustic data with a spatially varying absorption $μ_{a} \in [0.005, 0.04] {mm}^{- 1}$ and a constant scattering $μ_{s} = 4 {mm}^{- 1}$ . (a) Absorption coefficient $μ_{a} ({mm}^{- 1})$ , (b) logarithm of photon fluence $\log Φ$ (arbitrary units), and (c) absorbed optical energy density $H$ (arbitrary units). Image courtesy of N. Hänninen.

Fig. 4

Simulated photoacoustic data with a spatially varying scattering $μ_{s} \in [1, 8] {mm}^{- 1}$ and a constant absorption $μ_{a} = 0.02 {mm}^{- 1}$ . (a) Scattering coefficient $μ_{s} ({mm}^{- 1})$ , (b) logarithm of photon fluence $\log Φ$ (arbitrary units), and (c) absorbed optical energy density $H$ (arbitrary units). Image courtesy of N. Hänninen.

Different methods have been used to approach the optical inverse problem of QPAT. Here we give an overview of those. Many of the methods described here can be described as “iterative methods” since iterative optimization algorithms are utilized in their solution. In addition to these, so-called direct approaches for estimation of optical parameters have also been developed.¹⁴¹^–¹⁴³ For other reviews on the optical inverse problem of QPAT, see Refs. 2, 144, and 145.

There are two possible approaches to the optical inverse problem. Either the chromophore concentrations can be estimated directly from absorbed optical energy density data obtained using multiple wavelengths of light,⁶⁹^,¹⁴²^,¹⁴³^,¹⁴⁶^–¹⁴⁹ or they can be estimated in two stages, first by recovering the absorption coefficients at different wavelengths and then calculating the concentrations with a subsequent, linear, spectroscopic inversion.²^,¹⁴²^,¹⁴⁶^,¹⁴⁷ In both cases, the molar absorption spectra of the contributing chromophores need to be known a priori. Spectroscopic inversions are well understood mathematically, and they are used in many areas of optics, so we focus mostly in this section on the inversion for the optical parameters at a single wavelength. However, given its practical importance, a brief reminder of the spectroscopic inversion is given below.

3.2.1.

Spectroscopic inversion

A very important source of contrast in medical applications of PAT is hemoglobin. First, because it is easy to measure as it absorbs strongly in the near-infrared where the optical scattering is low and there are not too many other competing absorbers. Second, because—at least in principle—it gives a route to forming spatially resolved images of blood oxygen saturation, which is a key indicator of tissue function and pathology. The oxygen saturation is the ratio:

Eq. (25)

{sO}_{2} = C_{{HbO}_{2}} / (C_{{HbO}_{2}} + C_{Hb}),

where

C_{{HbO}_{2}}

and

C_{Hb}

are the concentrations of the chromophores oxyhemoglobin and deoxyhemoglobin. The link between these concentrations and the initial acoustic pressure distribution

p_{0} (r)

is given by Eqs. (1), (20), and (21). For the inversion of Eq. (1) to be well-conditioned, it is important to choose wavelengths in which the molar absorption coefficients form linearly independent vectors. It is usually assumed that the Grüneisen parameter is independent of wavelength in this scenario, and therefore cancels out when forming the ratio [Eq. (25)]. Concerningly, it is often assumed in experimental studies that the fluence too is independent of wavelength that the photoacoustic amplitudes are directly proportional to the absorption coefficients, but this cannot be the case. Because of the diffusive nature of the light propagation, the fluence will only be independent of wavelength if the optical properties of the tissue also remain constant with wavelength, in which case a spectroscopic inversion will not be possible. Solving the inverse problems, as described below, is therefore, critical to obtaining accurate estimates of this clinically important parameter.

3.2.2.

Uniqueness

Considering the optical inversion at a single wavelength, the simultaneous estimation of more than one optical parameter, e.g., absorption and scattering coefficients, is nonunique if only one light illumination pattern is used.¹⁴⁶^,¹⁵⁰ There are broadly two ways to overcome such ill-posedness: obtain more data to use in the inversion or make assumptions (perhaps based on adjunct data).

The simplest approach has been to assume that the scattering is known and to estimate only the absorption.¹⁵¹^–¹⁵⁶ Although widely used, the validity of this assumption is questionable in many cases of practical interest, e.g., when the scattering coefficient varies with tissue type or is not known accurately a priori. This approach can be improved by modeling the errors caused by the fixed scattering assumption by using a Bayesian approximation error modeling.¹⁵⁷ When using multiple wavelength data, assumptions can be made about the wavelength dependence¹⁴⁶^,¹⁵¹ to overcome the nonuniqueness. Several authors have reduced the parameter space by assuming that the optical properties are piecewise constant.⁷⁶^,¹⁵⁸^–¹⁶¹

In Ref. 141, it was shown that the nonuniqueness can be overcome using multiple optical illuminations to obtain more data. One way to achieve this would be by illuminating the target from different directions¹⁴¹^,¹⁴³^,¹⁵⁰^,¹⁶²^–¹⁶⁸ or using light patterns.¹⁶⁹ This approach has not yet been exploited experimentally, perhaps because experimentalists are concerned to obtain the maximum signal possible given the weakness of the photoacoustic effect, and so try to flood the tissue with light from as many directions as possible. Another approach that uses more data, but in this case also requires more hardware, is to combine QPAT with another modality to overcome the nonuniqueness, e.g., using DOT¹⁷⁰^–¹⁷³ or acousto-optic tomography¹⁷⁴ to estimate the fluence or using an additional modality to provide structural information⁷⁶^,¹⁶¹ for the light model.

As mentioned earlier, here we address estimation of optical parameters when absorbed optical energy density is given as data, which has been the approach in most research on the optical inverse problem of QPAT. This basically means that the acoustic inverse problem has been solved and the Grüneisen parameter is known. The Grüneisen parameter, however, is typically not known and to make the problem even more complicated, it is a spatially varying parameter that depends on chromophore concentration and temperature, and possibly even on wavelength as discussed in Sec. 2.2. One possibility to approach the problem of unknown Grüneisen parameter would be estimating it simultaneously with the optical parameters. In Ref. 142, it was shown that absorption coefficient, diffusion coefficient and the Grüneisen parameter can be estimated if data at multiple wavelengths are used and a prior information of the unknown parameters take specific forms. For example, the dependence of the coefficient on the spatial variable and on the wavelength variable needs to be separated for stable reconstruction. The recovery of the Grüneisen parameter simultaneously with the optical parameters in the case of multiwavelength optical inverse problem of QPAT has been addressed in Refs. 147 and 172. Further, estimation of absorption and scattering together with mapping of the temperature distribution has been studied in Ref. 175.

3.2.3.

Known scattering

When the scattering is known, a simple iteration is available for estimating the absorption coefficient. Rearranging Eq. (20) gives the fixed-point iteration:

Eq. (26)

μ_{a}^{(n + 1)} (r, λ) = H (r, λ) / Φ (r, λ; μ_{a} (r, λ)^{(n)}) .

Any suitable model can be used to calculate the fluence from the latest estimate of the absorption coefficient.¹⁵³^,¹⁵⁶ While straightforward to implement, this iteration can become sensitive to noise where the fluence is low. An improved method that finds

μ_{a}

by minimizing

‖ μ_{a} (r, λ) Φ (r, λ) - H (r, λ) ‖

is proposed in Ref. 176.

3.2.4.

Regularized least squares methods

Generally in parameter estimation problems, the unknown parameters of interest are determined by minimizing the least squares difference between the measured data and predictions of the forward model. Let us denote unknown optical parameters in a point $r$ by $x (r)$ . Further, denote the measurements by a finite dimensional vector $y \in R^{m}$ , where $m$ is the number of data points and the forward operator that maps the unknown parameters to data by $F (x (r))$ . Typically, in practical implementations, the parameters and the forward mapping are represented in discrete vector spaces $x (r) \mapsto x \in R^{n}$ , $F \mapsto f : R^{n} \mapsto R^{m}$ , where $n$ is the number of unknown parameters.

In the optical inverse problem, it has generally been assumed that the data $y$ are the absorbed optical energy density that has been obtained as the solution of the acoustic inverse problem. Most of the methods developed for regularized least squares and the Bayesian approach (described in the next section) in QPAT have used numerical approximation of the DA as the forward operator. In addition, the RTE has also been utilized.¹⁴³^,¹⁵⁶^,¹⁶⁶^,¹⁶⁷^,¹⁶⁹^,¹⁷⁷

Estimation of the unknown optical parameters can be written as a minimization problem:

Eq. (27)

\arg \min \frac{1}{2} {‖ y - f (x) ‖}^{2} .

In QPAT, estimation of absorption and scattering is generally ill-posed, and they cannot be solved using Eq. (27). The ill-posed nature of the problem can be alleviated by introducing a regularizing penalty term.¹⁷⁸ In that case, the minimization problem becomes

Eq. (28)

\arg \min \frac{1}{2} {‖ y - f (x) ‖}^{2} + B (x),

where

B (x)

is the regularizing penalty functional. The classical regularization methodology is to use Tikhonov regularization. In a widely applied generalized Tikhonov regularization, the penalty functional has the form:

Eq. (29)

B (x) = α {‖ D (x - x_{*}) ‖}_{2}^{2},

where

α

is a regularization parameter,

D

is a regularization matrix, and

x_{*}

is the prior estimate of

x

. Depending of the choice of the regularization matrix, Tikhonov regularization can enhance solutions with smaller norms by an identity matrix or smooth solutions through usage of difference matrix. In QPAT, Tikhonov regularization has been utilized in Refs. 141 142.–143 and 166.

Since many biological structures are better described as being piecewise constant, other regularizing norms (instead the second norm used in Tikhonov regularization) have been developed. An approach for supporting piecewise regular structures with sharp boundaries is to use total variation regularization that in a continuous formulation can be written as

Eq. (30)

B (x (r)) = β {‖ \nabla x (r) ‖}_{1},

that has been utilized in QPAT in Refs. 141, 164, 167, and 179. Furthermore, sparsity supporting

ℓ_{1}

regularization:

Eq. (31)

B (x (r)) = γ {‖ W x (r) ‖}_{1},

where

β

is a regularization parameter, and

W

is a sparse regularization operator has been utilized.¹⁸⁰ In addition, penalty acting as a physical limitation of the estimated parameters positivity constraint:

Eq. (32)

B (x)_{+} = Π_{k = 1}^{n} θ (x_{k}), θ (t) = {\begin{matrix} 1, & t \geq 0 \\ 0, & otherwise \end{matrix}

can be included into the minimization problem. A positivity constraint can be implemented in the minimization algorithm using projection or penalty methods.¹⁸¹

Alternatively to approaching the problem as a nonlinear optimization problem, in Refs. 182 and 183, it was proposed to extract absorption and photon fluence using a sparse signal representation. Then the resulting problem was approached as a linear problem of recovering the absorption and photon fluence without utilizing the forward model for light transport.

3.2.5.

Bayesian approach

In the framework of Bayesian inverse problems, the inverse problem is approached using statistical inference.¹¹⁵ All parameters are modeled as random variables, and information about them is expressed by probability distributions. The observation model with an additive noise model is of the form:

Eq. (33)

y = f (x) + e,

where

e \in R^{m}

denotes the noise. In the inverse problem, information about the parameters of primary interest is obtained based on the measurements, the model, and the prior information about the parameters. The solution of the inverse problem is the posterior probability distribution

π (x | y)

that according to Bayes’ theorem can be presented as a conditional probability density function of the form:

Eq. (34)

π (x | y) \propto π (y | x) π (x),

where

π (y | x)

is the likelihood density and

π (x)

is the prior density. Assuming the noise

e

and the unknown

x

mutually independent, observation model Eq. (33) leads to a likelihood density:

Eq. (35)

π (y | x) = π_{e} (y - f (x)),

where

π_{e} (\cdot)

is the probability density of the noise. Let us further model the noise

e

and unknown optical parameters

x

as Gaussian distributed, i.e.,

e \sim N (η_{e}, Γ_{e})

and

x \sim N (η_{x}, Γ_{x})

, where

η_{e}

and

Γ_{e}

are the mean and covariance of the noise, and

η_{x}

and

Γ_{x}

are the mean and covariance of the prior. In this case, the posterior probability density can be written as

Eq. (36)

π (x | y) \propto \exp {- \frac{1}{2} {‖ L_{e} (y - f (x) - η_{e}) ‖}^{2} - \frac{1}{2} {‖ L_{x} (x - η_{x}) ‖}^{2}},

where

L_{e}^{T} L_{e} = Γ_{e}^{- 1}

and

L_{x}^{T} L_{x} = Γ_{x}^{- 1}

are the square roots, such as the Cholesky decompositions of the inverse covariance matrices of the noise and prior, respectively.

In principle, the distributions of the unknown parameters can be estimated using Markov chain Monte Carlo methods.¹⁸⁴ However, these methods can be prohibitively computationally too expensive in large dimensional tomographic inverse problems. Therefore, point estimates are computed to approximate the posterior distribution. An often considered point estimate in tomographic imaging is the maximum a posteriori (MAP) estimate

Eq. (37)

x_{MAP} = \arg \min_{x} {\frac{1}{2} {‖ L_{e} (y - f (x) - η_{e}) ‖}^{2} + \frac{1}{2} {‖ L_{x} (x - η_{x}) ‖}^{2}} .

If the forward operator is linear and in the case of Gaussian noise and prior, the MAP estimate is also the mean of the posterior. However, as mentioned earlier, optical inverse problem of QPAT is nonlinear.

If the entire posterior distribution could be solved, it can be used to evaluate the reliability of the estimated parameters through examining the standard deviation. In the case of a nonlinear inverse problems, where (only) the MAP estimate has been solved, the standard deviation can be approximated through Laplace’s approximation, leading to a Gaussian approximation for the posterior distribution in the location of the MAP estimate

x | y \sim N (\hat{x}, \hat{Γ}),

where

\hat{x}

is the MAP estimate and

Eq. (38)

\hat{Γ} = {(J (\hat{x})^{T} Γ_{e}^{- 1} J (\hat{x}) + Γ_{x}^{- 1})}^{- 1}

is the covariance, where

J (\hat{x})

is the Jacobian matrix of

f (x)

.

The Bayesian approach provides a natural methodology for taking into account the uncertainties in parameters, data, and models; thus Bayesian approximation error modeling¹¹⁵ has been utilized in modeling of uncertainties in various applications. The Bayesian approach was formulated for the optical inverse problem of QPAT in Ref. 116, where also modeling of noise and errors due to the acoustic solver were studied. The approach has been further utilized in marginalization of scattering in Ref. 157 and it has been extended for spectral QPAT in Ref. 147.

3.2.6.

Learned/data-driven methods

As for the acoustic inverse problem described above, and biomedical imaging in general, data-driven learning-based approached have gained interest also for the optical inverse problem of QPAT (see Refs. 121, 122, and 185 186.–187 and the references therein). Many of the studies are numerical simulation studies using U-Nets or variants to map directly from synthetic maps of absorbed energy density to optical coefficients of interest, including blood oxygenation.¹⁸⁸^–¹⁹⁶ A hybrid approach is to use a learned segmentation algorithm to determine the tissue model and then model the fluence as above.¹⁵⁹ Although these numerical studies seem promising, the domain gap between the simulated datasets used for training and real measured data casts doubt on how well these methods will translate to experimental settings. Some attempts have been made to address this recently. One approach tackles the issue by aiming to improve the quality of synthetic datasets through the use of detailed anatomical digital phantoms.¹⁹⁷^,¹⁹⁸ Another promising approach is to make measured datasets on physical phantoms and use those for training.¹⁹¹ A third attempt tries to close the domain gap using GANs to move the simulated data distribution closer to the real data distribution.¹⁹⁹^,²⁰⁰

3.2.7.

Utilizing Monte Carlo method in the inverse problem of QPAT

Recently, utilization of Monte Carlo method for light transport in the optical inverse problem of QPAT has raised interest. The methodology has been utilized in estimating absorption while assuming scattering as known.²⁰¹^,²⁰² Utilizing Monte Carlo in the minimization approaches, such as regularized least squares or MAP estimation, for estimating both absorption and scattering in QPAT requires evaluation of a minimization direction, such as the steepest descent or Gauss–Newton direction. The gradients for a steepest descent algorithm can be obtained using a so-called adjoint Monte Carlo.⁶² The approach has been utilized in estimating absorption⁶¹^,⁶²^,²⁰³ and it has also been evaluated with experimental data.²⁰⁴ Alternatively, Jacobian matrices can be constructed by computing the derivative for the absorption coefficient directly from Eq. (19) by differentiation and computing the derivative for the scattering coefficient utilizing a so-called perturbation Monte Carlo.⁶⁴ The approach has been utilized in QPAT for estimation of both absorption and scattering simultaneously.⁶⁴

Due to the stochastic nature of Monte Carlo, the minimization directions evaluated using Monte Carlo simulations are stochastic. Multiple evaluations of Monte Carlo simulation can easily create a bottleneck for the attempts to utilize the method in the solution of the inverse problem. In order to overcome this problem, it was proposed in Ref. 205 to utilize stochastic gradient methods from the machine learning community. In the approach, the number of photon packets utilized in evaluation of the minimization direction is adaptively adjusted based on examining variations of these minimization directions. The approach was utilized in estimation of absorption coefficient using the steepest descent method in Ref. 205, and it was extended to a 2D imaging geometry with a Gauss–Newton algorithm in Ref. 66 and for estimating absorption and scattering simultaneously in Ref. 206.

3.3.

Single-Stage Approaches

In addition to the two-stage approach, where acoustic and optical inverse problems of QPAT are solved consecutively, estimation of the optical parameters directly from the photoacoustic time-series has been proposed. A single-stage approach was formulated utilizing the Born approximation in Refs. 207 and 208. Furthermore, the $ℓ_{1}$ sparsity regularization was utilized,²⁰⁹ where the minimization problem was solved with BFGS algorithm also in a limited view measurement geometry, and Tikhonov regularization with a proximal gradient algorithm was utilized in Ref. 177. A single-stage approach for estimating spatially varying speed of sound and optical parameters simultaneously was presented in Ref. 210. A Bayesian approach to determine optical parameters and evaluate the reliability of the estimates was proposed in Ref. 211 and studied also in a limited view measurement geometry and 3D imaging situation. Furthermore, in Ref. 212, stochastic search algorithms were proposed for single-stage QPAT.

4. Practical Considerations

This paper is principally concerned with the mathematics of the optical inversions in QPAT, but because inverse problems of this type are inherently motivated by experimental measurements, a brief discussion of the trade-offs and limitations that face experimentalists might help inform future research in this area. The first two subsections below are concerned with the quality and completeness of the acoustic data, which will directly affect the quality of the estimates of optical properties that can be obtained.

4.1.

Detector Response Versus Sensitivity

Photoacoustic pulses, because of their impulsive nature, are broadband and, because they emanate from a distributed source, may arrive at the detectors from any angle. It was mentioned in Sec. 2.3 that real detectors have the effect of filtering the acoustic pressure both temporally (frequency response) and spatially (directionality). Typically, they are more sensitive over a certain bandwidth (centered at one frequency) and over a range of angles (centered on one angle). There are ways in which the range of frequencies or angles can be increased, but the trade-off is often a drop in the detection sensitivity. As the magnitude of the photoacoustic signals is ultimately limited by safety factors controlling the maximum permissible fluence, signal detection cannot be improved by simply turning up the power. It is often the case, then that photoacoustic measurements are made with detectors whose bandwidths do not capture the full range of available frequencies, and whose directional response does not extend to the steepest angles. Although these responses can be modeled, to some extent, and deconvolved from the data, where the signal has fallen well into the noise it will not be recoverable. The loss of the low frequencies in the data will, without amelioration, lead to a loss of low spatial frequencies in the photoacoustic image, and the loss of the waves arriving from steeper angles can have an effect on the shape or amplitude of the image and lead to “limited-view”-type artifacts.

4.2.

Array Coverage Versus Cost and Complexity

The ideal detection array would consist of point-like detectors distributed over a surface surrounding the object of interest (spaced at half the shortest acoustic wavelength to satisfy Nyquist), simultaneously streaming data. Practical arrays are prevented from reaching this ideal due to considerations over the cost and complexity of array fabrication, the cost of high-channel count data acquisition systems, and, in some cases, the physical size of the detectors making up the array. Most arrays therefore cover just part of the $2 π$ radians necessary to ensure complete data,²¹³ and in some cases, e.g., linear arrays, quite a small fraction, which can lead to “limited-view” artifacts in the image. This hardware limitation can, in some cases, be compensated through scanning or rotating the array to different positions to take new measurements or using compressed sensing strategies, although this increases data acquisition time.

4.3.

Optical Excitation

How the object is illuminated, both in terms of intensity and coverage, is key to obtaining data with a good signal-to-noise ratio. There are safety limits on the irradiance that soft tissue can be exposed to so illumination over as much surface as is available will maximize the total optical energy entering the tissue. When choosing wavelengths is critical to ensure the well-posedness of the spectroscopic inverse problem, but there is another consideration. Because the depth to which the light will reach depends on the absorbers and scatterers within the tissue, it is not always obvious a priori what the optimum optical wavelengths will be in terms of depth penetration, so finding the optimal wavelengths may require trial-and-error. Pulse-to-pulse and wavelength-dependent variations in laser energy must be known but can usually be measured and corrected for. The pulse repetition rate of the laser used affects the data acquisition time (even for systems that use arrays so all the data for one image can be obtained from one pulse, each wavelength requires a different pulse), so the faster the pulses the smaller the risk of tissue motion artifacts in the images. Finally, (this links to the next section) it is important to know how the tissue is being illuminated so that it can be modeled accurately in solving the inverse problem, so the direction (e.g., divergence angle) and intensity of the incident light distribution must be measured or estimated.

4.4.

Knowledge of Model Parameters

Even when the data are completed, the solution of inverse problem can fail if the auxiliary parameters that are fed into the models are not accurate. Some, such as the positions of the detectors, can be determined through calibration procedures, although this is harder for others, e.g., the positions, directions, and angular spread of the optical sources (mentioned above). Others may be obtainable using another imaging modality, e.g., obtaining the speed of sound using ultrasound tomography. Yet other parameters can be measured on characteristic samples, e.g., the optical anisotropy or Grüneisen parameters, although there will clearly be uncertainties both in these measurements and in how well the samples match the object. In general, all predefined parameters can contain, more or less, uncertainties that will effect on the solution of the inverse problem. Modeling of such uncertainties is one of the research interests in the community of the inverse problems research, and some approaches for tackling the problem Bayesian approximation error modeling¹¹⁵ and utilizing machine learning,¹³⁵ have been developed. For example in Ref. 119, it was shown that uncertainties in modeling of ultrasound sensor locations can be compensated using the Bayesian approach.

4.5.

Computational Considerations

Choosing a solution methodology for an inverse problem is often a trade-off between accuracy of the solution and computational cost. The most simple and straightforward approaches, such as methodologies utilizing analytical solutions of the models, are often limited to specific geometries and/or make unrealistic assumptions of the parameters, assuming scattering as known in the optical inverse problem of QPAT. On the other hand, numerical approximations in large 3D volumes can be computationally expensive, and inverse problems solution methodologies utilizing them require advanced approaches, such as model reduction techniques, to be feasible in practical imaging scenarios. Overall, in practical imaging techniques, it would be desirable that the related modeling and numerical inverse problem methodologies could be implemented such that they could be solved in standard computers without significant time delays. Therefore, although QPAT is a promising methodology for providing high-resolution 3D images of physiologically relevant parameters, there are many computational modeling-based challenges that need to be tackled before the technique can be developed as a standard clinical or preclinical tool.

5. Summary

In this survey, the modeling of photoacoustics and inverse problems methodologies for image reconstruction in QPAT were reviewed. As mentioned, QPAT consists of two parts: optical and acoustic, and here we focused on the mathematics of the optical part. Approaches to modeling light propagation in tissue, described by transport theory, and methodologies for the solution of the optical inverse problem were described, but the acoustic aspects of QPAT were discussed only to the extent that they relate to the optical part. Although modeling of light transport is well-understood and there is a well-developed framework of inverse mathematics for approaching the inverse problem of QPAT, there are still challenges in taking these methodologies to practice. For example, although the effect of the variation of the light fluence throughout tissue on photoacoustic tomography images is well-understood, there is no consensus among photoacoustic practitioners on how to remove the effect of the fluence, and thereby facilitate accurate estimation of the tissue’s optical properties.

There are several reasons for this. First, and perhaps most importantly, QPAT is an ill-posed inverse problem. In practice, this means that even small errors in measurements or modeling can cause large errors in the solution of the inverse problem. Related to this is the fact that, while the effect of the optical scattering on the absorbed optical energy density is weaker than that of the absorption, it does not seem to be weak enough to ignore. Second, photoacoustic data are not ideal, as described in Sec. 4, which is a challenging starting point. Third, numerical methods to compute fluence and to solve the inverse problem can be computationally intensive even in simplified imaging situations. Finally, the auxiliary parameters that are needed as inputs to the models (see Fig. 1) can contain uncertainties, making the computational modeling even more difficult. In particular, the Grüneisen parameter may not be known accurately, and methods for compensating its effect on modeling and inverse problem are still being developed.

In this review, we touched these key challenges while giving an overview of computational modeling and inverse problem of QPAT. Furthermore, the developments for tackling these problems were reviewed. A comprehensive literature review was not provided, but we hope that the references can provide a route into broader literature for those who are interested.

Disclosures

The authors declare no conflicts of interest.

Code and Data Availability

The data that support the findings of this study are available upon reasonable request from the authors.

Acknowledgments

This work was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Program (Grant Agreement No. 101001417-QUANTOM) and the Research Council of Finland [Centre of Excellence in Inverse Modeling and Imaging project (No. 353086) and the Flagship Program Photonics Research and Innovation (Grant No. 320166)]. The authors would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the program Rich and Nonlinear Tomography where work on this paper was undertaken, supported by EPSRC (Grant No. EP/R014604/1).

References

1.

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

2.

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

3.

C. Li and L. V. Wang, “Photoacoustic tomography and sensing in biomedicine,” Phys. Med. Biol., 54 R59 –R97 https://doi.org/10.1088/0031-9155/54/19/R01 PHMBA7 0031-9155 (2009). Google Scholar

4.

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

5.

A. Ishimaru, Wave Propagation and Scattering in Random Media, 1 Academic Press, New York (1978). Google Scholar

6.

L.W. Wang and H. Wu, Biomedical Optics: Principles and Imaging, John Wiley & Sons( (2012). Google Scholar

7.

G. Bal, “Radiative transfer equation with varying refractive index: a mathematical perspective,” J. Opt. Soc. Am. A, 23 1639 –1644 https://doi.org/10.1364/JOSAA.23.001639 JOAOD6 0740-3232 (2006). Google Scholar

8.

O. Lehtikangas et al., “Finite element approximation of the radiative transport equation in a medium with piece-wise constant refractive index,” J. Comp. Phys., 282 345 –359 https://doi.org/10.1016/j.jcp.2014.11.025 (2015). Google Scholar

9.

M. L. Shendeleva, “Radiative transfer in a turbid medium with a varying refractive index: comment,” J. Opt. Soc. Am. A, 21 (12), 2464 –2467 https://doi.org/10.1364/JOSAA.21.002464 (2004). Google Scholar

10.

J.-M. Tualle and E. Tinet, “Derivation of the radiative transfer equation for scattering media with a spatially varying refractive index,” Opt. Commun., 228 (1–3), 33 –38 https://doi.org/10.1016/j.optcom.2003.09.076 OPCOB8 0030-4018 (2003). Google Scholar

11.

A. Sassaroli and F. Martelli, “Equivalence of four Monte Carlo methods for photon migration in turbid media,” J. Opt. Soc. Am. A, 29 (10), 2110 –2117 https://doi.org/10.1364/JOSAA.29.002110 JOAOD6 0740-3232 (2012). Google Scholar

12.

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

13.

I. J. Bigio and S. Fantini, Quantitative Biomedical Optics: Theory, Methods, and Applications, Cambridge University Press( (2016). Google Scholar

14.

S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol., 58 (11), R37 https://doi.org/10.1088/0031-9155/58/11/R37 PHMBA7 0031-9155 (2013). Google Scholar

15.

M. I. Mishchenko, “Maxwell’s equations, radiative transfer, and coherent backscattering: a general perspective,” J. Quant. Spectrosc. Radiat. Transf., 101 (3), 540 –555 https://doi.org/10.1016/j.jqsrt.2006.02.065 JQSRAE 0022-4073 (2006). Google Scholar

16.

J. Ripoll, “Derivation of the scalar radiative transfer equation from energy conservation of Maxwell’s equations in the far field,” J. Opt. Soc. Am. A, 28 (8), 1765 –1775 https://doi.org/10.1364/JOSAA.28.001765 (2011). Google Scholar

17.

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

18.

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

19.

J. Heino et al., “Anisotropic effects in highly scattering media,” Phys. Rev. E, 68 (3), 031908 https://doi.org/10.1103/PhysRevE.68.031908 (2003). Google Scholar

20.

A. Kienle, F. K. Forster and R. Hibst, “Influence of the phase function on determination of the optical properties of biological tissue by spatially resolved reflectance,” Opt. Lett., 26 (20), 1571 –1573 https://doi.org/10.1364/OL.26.001571 OPLEDP 0146-9592 (2001). Google Scholar

21.

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

22.

M. Schweiger et al., “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys., 22 (11), 1779 –1792 https://doi.org/10.1118/1.597634 MPHYA6 0094-2405 (1995). Google Scholar

23.

R. Hochuli et al., “Estimating blood oxygenation from photoacoustic images: can a simple linear spectroscopic inversion ever work?,” J. Biomed. Opt., 24 (12), 121914 https://doi.org/10.1117/1.JBO.24.12.121914 JBOPFO 1083-3668 (2019). Google Scholar

24.

G.-S. Jeng et al., “Real-time interleaved spectroscopic photoacoustic and ultrasound (PAUS) scanning with simultaneous fluence compensation and motion correction,” Nat. Commun., 12 (1), 716 https://doi.org/10.1038/s41467-021-20947-5 NCAOBW 2041-1723 (2021). Google Scholar

25.

M. W. Kim et al., “Correction of wavelength-dependent laser fluence in swept-beam spectroscopic photoacoustic imaging with a hand-held probe,” Photoacoustics, 19 100192 https://doi.org/10.1016/j.pacs.2020.100192 (2020). Google Scholar

26.

S. Park et al., “Normalization of optical fluence distribution for three-dimensional functional optoacoustic tomography of the breast,” J. Biomed. Opt., 27 (3), 036001 https://doi.org/10.1117/1.JBO.27.3.036001 JBOPFO 1083-3668 (2022). Google Scholar

27.

V. V. Perekatova et al., “Quantitative techniques for extraction of blood oxygenation from multispectral optoacoustic measurements,” Laser Phys. Lett., 16 (11), 116201 https://doi.org/10.1088/1612-202X/ab4dab 1612-2011 (2019). Google Scholar

28.

J. Zalev et al., “Opto-acoustic imaging of relative blood oxygen saturation and total hemoglobin for breast cancer diagnosis,” J. Biomed. Opt., 24 (12), 121915 https://doi.org/10.1117/1.JBO.24.12.121915 JBOPFO 1083-3668 (2019). Google Scholar

29.

A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys., 220 (1), 441 –470 https://doi.org/10.1016/j.jcp.2006.07.007 JCTPAH 0021-9991 (2006). Google Scholar

30.

S. R. Arridge et al., “The finite element model for the propagation of light in scattering media: a direct method for domains with nonscattering regions,” Med. Phys., 27 (1), 252 –264 https://doi.org/10.1118/1.598868 MPHYA6 0094-2405 (2000). Google Scholar

31.

T. Hayashi, Y. Kashio and E. Okada, “Hybrid Monte Carlo-diffusion method for light propagation in tissue with a low-scattering region,” Appl. Opt., 42 (16), 2888 –2896 https://doi.org/10.1364/AO.42.002888 APOPAI 0003-6935 (2003). Google Scholar

32.

O. Lehtikangas and T. Tarvainen, “Hybrid forward-peaked-scattering-diffusion approximations for light propagation in turbid media with low-scattering regions,” J. Quant. Spectrosc. Radiat. Transf., 116 132 –144 https://doi.org/10.1016/j.jqsrt.2012.10.017 JQSRAE 0022-4073 (2013). Google Scholar

33.

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

34.

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 (8), 1746 –1752 https://doi.org/10.1364/JOSAA.10.001746 JOAOD6 0740-3232 (1993). Google Scholar

35.

C. Frederick, K. Ren and S. Vallélian, “Image reconstruction in quantitative photoacoustic tomography with the simplified P2 approximation,” SIAM J. Imaging Sci., 11 (4), 2847 –2876 https://doi.org/10.1137/18M1195656 (2018). Google Scholar

36.

Y. Wang et al., “Nonlinear iterative perturbation scheme with simplified spherical harmonics (SP3) light propagation model for quantitative photoacoustic tomography,” J. Biophotonics, 14 (6), e202000446 https://doi.org/10.1002/jbio.202000446 (2021). Google Scholar

37.

P. González-Rodríguez and A. D. Kim, “Comparison of light scattering models for diffuse optical tomography,” Opt. Express, 17 (11), 8756 –8774 https://doi.org/10.1364/OE.17.008756 OPEXFF 1094-4087 (2009). Google Scholar

38.

T. Saratoon et al., “3D quantitative photoacoustic tomography using the

δ

-Eddington approximation,” Proc. SPIE, 8581 85810V https://doi.org/10.1117/12.2004105 PSISDG 0277-786X (2013). Google Scholar

39.

S. R. Arridge et al., “A finite element approach for modeling photon transport in tissue,” Med. Phys., 20 (2), 299 –309 https://doi.org/10.1118/1.597069 MPHYA6 0094-2405 (1993). Google Scholar

40.

H. Dehghani et al., “Near infrared optical tomography using NIRFAST: algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng., 25 (6), 711 –732 https://doi.org/10.1002/cnm.1162 CANMER 0748-8025 (2009). Google Scholar

41.

M. Schweiger and S. Arridge, “The Toast++ software suite for forward and inverse modeling in optical tomography,” J. Biomed. Opt., 19 (4), 040801 https://doi.org/10.1117/1.JBO.19.4.040801 JBOPFO 1083-3668 (2014). Google Scholar

42.

O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Probl., 14 (5), 1107 https://doi.org/10.1088/0266-5611/14/5/003 INPEEY 0266-5611 (1998). Google Scholar

43.

A. H. Hielscher, R. E. Alcouffe and R. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol., 43 (5), 1285 https://doi.org/10.1088/0031-9155/43/5/017 PHMBA7 0031-9155 (1998). Google Scholar

44.

E. D. Aydin, C. R. E. De Oliveira and A. J. H. Goddard, “A comparison between transport and diffusion calculations using a finite element-spherical harmonics radiation transport method,” Med. Phys., 29 (9), 2013 –2023 https://doi.org/10.1118/1.1500404 MPHYA6 0094-2405 (2002). Google Scholar

45.

S. Wright, M. Schweiger and S. R. Arridge, “Reconstruction in optical tomography using the PN approximations,” Meas. Sci. Technol., 18 (1), 79 –86 https://doi.org/10.1088/0957-0233/18/1/010 MSTCEP 0957-0233 (2006). Google Scholar

46.

K. Ren, G. Bal and A. H. Hielscher, “Frequency domain optical tomography based on the equation of radiative transfer,” SIAM J. Sci. Comput., 28 (4), 1463 –1489 https://doi.org/10.1137/040619193 SJOCE3 1064-8275 (2006). Google Scholar

47.

G. Kanschat, “A robust finite element discretization for radiative transfer problems with scattering,” East-West J. Numer. Math., 6 (4), 265 –272 EJNMEA (1998). Google Scholar

48.

O. Balima et al., “Optical tomography with the discontinuous Galerkin formulation of the radiative transfer equation in frequency domain,” J. Quant. Spectrosc. Radiat. Transf., 113 (10), 805 –814 https://doi.org/10.1016/j.jqsrt.2012.03.003 JQSRAE 0022-4073 (2012). Google Scholar

49.

S. Powell, B. T. Cox and S. R. Arridge, “A pseudospectral method for solution of the radiative transport equation,” J. Comput. Phys., 384 376 –382 https://doi.org/10.1016/j.jcp.2019.01.024 JCTPAH 0021-9991 (2019). Google Scholar

50.

M. A. Badri et al., “Preconditioned Krylov subspace methods for solving radiative transfer problems with scattering and reflection,” Comput. Math. Appl., 77 (6), 1453 –1465 https://doi.org/10.1016/j.camwa.2018.09.041 CMAPDK 0898-1221 (2019). Google Scholar

51.

H. Gao, L. Phan and Y. Lin, “Parallel multigrid solver of radiative transfer equation for photon transport via graphics processing unit,” J. Biomed. Opt., 17 (9), 096004 https://doi.org/10.1117/1.JBO.17.9.096004 JBOPFO 1083-3668 (2012). Google Scholar

52.

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

53.

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

54.

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

55.

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 https://doi.org/10.1364/OE.17.020178 OPEXFF 1094-4087 (2009). Google Scholar

56.

J. Cassidy et al., “High-performance, robustly verified Monte Carlo simulation with FullMonte,” J. Biomed. Opt., 23 (8), 085001 https://doi.org/10.1117/1.JBO.23.8.085001 JBOPFO 1083-3668 (2018). Google Scholar

57.

D. Côté and I.A. Vitkin, “Robust concentration determination of optically active molecules in turbid media with validated three-dimensional polarization sensitive Monte Carlo calculations,” Opt. Express, 13 (1), 148 –163 https://doi.org/10.1364/OPEX.13.000148 OPEXFF 1094-4087 (2005). Google Scholar

58.

A.A. Leino, A. Pulkkinen and T. Tarvainen, “ValoMC: a Monte Carlo software and MATLAB toolbox for simulating light transport in biological tissue,” OSA Contin., 2 957 –972 https://doi.org/10.1364/OSAC.2.000957 (2019). Google Scholar

59.

Y. Liu et al., “OptogenSIM: a 3D Monte Carlo simulation platform for light delivery design in optogenetics,” Biomed. Opt. Express, 6 (12), 4859 –4870 https://doi.org/10.1364/BOE.6.004859 BOEICL 2156-7085 (2015). Google Scholar

60.

O. Yang and B. Choi, “Accelerated rescaling of single Monte Carlo simulation runs with the Graphics Processing Unit (GPU),” Biomed. Opt. Express, 4 (11), 2667 –2672 https://doi.org/10.1364/BOE.4.002667 BOEICL 2156-7085 (2013). Google Scholar

61.

J. Buchmann et al., “Three-dimensional quantitative photoacoustic tomography using an adjoint radiance Monte Carlo model and gradient descent,” J. Biomed. Opt., 24 (6), 066001 https://doi.org/10.1117/1.JBO.24.6.066001 JBOPFO 1083-3668 (2019). Google Scholar

62.

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

63.

S. Zheng et al., “Quantitative photoacoustic tomography with light fluence compensation based on radiance Monte Carlo model,” Phys. Med. Biol., 68 (6), 065009 https://doi.org/10.1088/1361-6560/acbe90 PHMBA7 0031-9155 (2023). Google Scholar

64.

A. A. Leino et al., “Perturbation Monte Carlo method for quantitative photoacoustic tomography,” IEEE Trans. Med. Imaging, 39 (10), 2985 –2995 https://doi.org/10.1109/TMI.2020.2983129 ITMID4 0278-0062 (2020). Google Scholar

65.

I. Lux and L. Koblinger, Monte Carlo Particle Transport Methods: Neutron and Photon Calculations, CRC Press( (1991). Google Scholar

66.

N. Hänninen et al., “Adaptive stochastic Gauss–Newton method with optical Monte Carlo for quantitative photoacoustic tomography,” J. Biomed. Opt., 27 (8), 083013 https://doi.org/10.1117/1.JBO.27.8.083013 JBOPFO 1083-3668 (2022). Google Scholar

67.

T Stahl et al., “Novel method for the analysis of clathrates,” Refrigeration Science and Technology, 2020 306 –313 Institut International du Froid: IIF( (2020). Google Scholar

68.

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

69.

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

70.

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

71.

J. Xia et al., “Ultrasonic-heating-encoded photoacoustic tomography with virtually augmented detection view,” Optica, 2 (4), 307 –312 https://doi.org/10.1364/OPTICA.2.000307 (2015). Google Scholar

72.

M. Pramanik and L. V. Wang, “Thermoacoustic and photoacoustic sensing of temperature,” J. Biomed. Opt., 14 (5), 054024 https://doi.org/10.1117/1.3247155 JBOPFO 1083-3668 (2009). Google Scholar

73.

X. Jin and L. V. Wang, “Thermoacoustic tomography with correction for acoustic speed variations,” Phys. Med. Biol., 51 (24), 6437 https://doi.org/10.1088/0031-9155/51/24/010 PHMBA7 0031-9155 (2006). Google Scholar

74.

S. Manohar et al., “Concomitant speed-of-sound tomography in photoacoustic imaging,” Appl. Phys. Lett., 91 (13), 131911 https://doi.org/10.1063/1.2789689 APPLAB 0003-6951 (2007). Google Scholar

75.

E. Merčep et al., “Whole-body live mouse imaging by hybrid reflection-mode ultrasound and optoacoustic tomography,” Opt. Lett., 40 (20), 4643 –4646 https://doi.org/10.1364/OL.40.004643 OPLEDP 0146-9592 (2015). Google Scholar

76.

A. Pattyn et al., “Model-based optical and acoustical compensation for photoacoustic tomography of heterogeneous mediums,” Photoacoustics, 23 100275 https://doi.org/10.1016/j.pacs.2021.100275 (2021). Google Scholar

77.

J. Xia et al., “Enhancement of photoacoustic tomography by ultrasonic computed tomography based on optical excitation of elements of a full-ring transducer array,” Opt. Lett., 38 (16), 3140 –3143 https://doi.org/10.1364/OL.38.003140 OPLEDP 0146-9592 (2013). Google Scholar

78.

C. Huang et al., “Photoacoustic computed tomography correcting for heterogeneity and attenuation,” J. Biomed. Opt., 17 (6), 061211 https://doi.org/10.1117/1.JBO.17.6.061211 JBOPFO 1083-3668 (2012). Google Scholar

79.

R. Kowar and O. Scherzer, “Photoacoustic imaging taking into account attenuation,” (2010). https://doi.org/10.48550/arXiv.1009.4350 Google Scholar

80.

B. E. Treeby et al., “Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method,” J. Acoust. Soc. Am., 131 (6), 4324 –4336 https://doi.org/10.1121/1.4712021 JASMAN 0001-4966 (2012). Google Scholar

81.

B. E. Treeby, E. Z. Zhang and B. T. Cox, “Photoacoustic tomography in absorbing acoustic media using time reversal,” Inverse Probl., 26 115003 https://doi.org/10.1088/0266-5611/26/11/115003 INPEEY 0266-5611 (2010). Google Scholar

82.

J. J. Kaufman, G. Luo and R. S. Siffert, “Ultrasound simulation in bone,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 55 (6), 1205 –1218 https://doi.org/10.1109/TUFFC.2008.784 ITUCER 0885-3010 (2008). Google Scholar

83.

K. Mitsuhashi et al., “A forward-adjoint operator pair based on the elastic wave equation for use in transcranial photoacoustic computed tomography,” SIAM J. Imaging Sci., 10 (4), 2022 –2048 https://doi.org/10.1137/16M1107619 (2017). Google Scholar

84.

A. Pulkkinen et al., “Numerical simulations of clinical focused ultrasound functional neurosurgery,” Phys. Med. Biol., 59 (7), 1679 –1700 https://doi.org/10.1088/0031-9155/59/7/1679 PHMBA7 0031-9155 (2014). Google Scholar

85.

P. J. White, G. T. Clement and K. Hynynen, “Longitudinal and shear mode ultrasound propagation in human skull bone,” Ultrasound Med. Biol., 32 (7), 1085 –1096 https://doi.org/10.1016/j.ultrasmedbio.2006.03.015 USMBA3 0301-5629 (2006). Google Scholar

86.

T. Else et al., “PATATO: a Python photoacoustic tomography analysis toolkit,” Proc. SPIE, PC12379 PC123790T https://doi.org/10.1117/12.2648830 PSISDG 0277-786X (2023). Google Scholar

87.

C. Fadden and S.-R. Kothapalli, “A single simulation platform for hybrid photoacoustic and RF-acoustic computed tomography,” Appl. Sci., 8 (9), 1568 https://doi.org/10.3390/app8091568 (2018). Google Scholar

88.

M. T. Graham and M. A. Lediju Bell, “Phocospace: an open-source simulation package to implement photoacoustic spatial coherence theory,” in IEEE Int. Ultrason. Symp. (IUS), 1 –5 (2022). https://doi.org/10.1109/IUS54386.2022.9958309 Google Scholar

89.

J. Gröhl et al., “SIMPA: an open-source toolkit for simulation and image processing for photonics and acoustics,” J. Biomed. Opt., 27 (8), 083010 https://doi.org/10.1117/1.JBO.27.8.083010 JBOPFO 1083-3668 (2022). Google Scholar

90.

P. Omidi et al., “PATLAB: a graphical computational software package for photoacoustic computed tomography research,” Photoacoustics, 28 100404 https://doi.org/10.1016/j.pacs.2022.100404 (2022). Google Scholar

91.

D. O’Kelly et al., “A scalable open-source MATLAB toolbox for reconstruction and analysis of multispectral optoacoustic tomography data,” Sci. Rep., 11 (1), 19872 https://doi.org/10.1038/s41598-021-97726-1 SRCEC3 2045-2322 (2021). Google Scholar

92.

B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt., 15 (2), 021314 https://doi.org/10.1117/1.3360308 JBOPFO 1083-3668 (2010). Google Scholar

93.

K. A. Wear, C. Baker and P. Miloro, “Directivity and frequency-dependent effective sensitive element size of needle hydrophones: predictions from four theoretical forms compared with measurements,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 65 (10), 1781 –1788 https://doi.org/10.1109/TUFFC.2018.2855967 ITUCER 0885-3010 (2018). Google Scholar

94.

K. A. Wear, C. Baker and P. Miloro, “Directivity and frequency-dependent effective sensitive element size of membrane hydrophones: theory versus experiment,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 66 (11), 1723 –1730 https://doi.org/10.1109/TUFFC.2019.2930042 ITUCER 0885-3010 (2019). Google Scholar

95.

H. Zuo et al., “Spectral crosstalk in photoacoustic computed tomography,” Photoacoustics, 26 100356 https://doi.org/10.1016/j.pacs.2022.100356 (2022). Google Scholar

96.

D. Finch, S. K. Patch and Rakesh, “Determining a function from its mean values over a family of spheres,” SIAM J. Math. Anal., 35 (5), 1213 –1240 https://doi.org/10.1137/S0036141002417814 SJMAAH 0036-1410 (2004). Google Scholar

97.

L. A. Kunyansky, “Explicit inversion formulae for the spherical mean Radon transform,” Inverse Probl., 23 (1), 373 –383 https://doi.org/10.1088/0266-5611/23/1/021 INPEEY 0266-5611 (2007). Google Scholar

98.

M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E, 71 (1), 016706 https://doi.org/10.1103/PhysRevE.71.016706 (2005). Google Scholar

99.

M. Agranovsky and P. Kuchment, “Uniqueness of reconstruction and an inversion procedure for thermoacoustic and photoacoustic tomography with variable sound speed,” Inverse Probl., 23 (5), 2089 –2102 https://doi.org/10.1088/0266-5611/23/5/016 INPEEY 0266-5611 (2007). Google Scholar

100.

L. A. Kunyansky, “A series solution and a fast algorithm for the inversion of the spherical mean Radon transform,” Inverse Probl., 23 (6), S11 –S20 https://doi.org/10.1088/0266-5611/23/6/S02 INPEEY 0266-5611 (2007). Google Scholar

101.

P. Burgholzer et al., “Exact and approximative imaging methods for photoacoustic tomography using an arbitrary detection surface,” Phys. Rev. E, 75 (4), 046706 https://doi.org/10.1103/PhysRevE.75.046706 (2007). Google Scholar

102.

B. T. Cox and B. E. Treeby, “Artifact trapping during time reversal photoacoustic imaging for acoustically heteregeneous media,” IEEE Trans. Med. Imaging, 29 (2), 387 –396 https://doi.org/10.1109/TMI.2009.2032358 ITMID4 0278-0062 (2010). Google Scholar

103.

Y. Hristova, P. Kuchment and L. Nguyen, “Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media,” Inverse Probl., 24 055006 https://doi.org/10.1088/0266-5611/24/5/055006 INPEEY 0266-5611 (2008). Google Scholar

104.

Y. Xu and L. V. Wang, “Time reversal and its application to tomography with diffracting sources,” Phys. Rev. Lett., 92 (3), 339021 –339024 https://doi.org/10.1103/PhysRevLett.92.033902 PRLTAO 0031-9007 (2004). Google Scholar

105.

S. Arridge et al., “Accelerated high-resolution photoacoustic tomography via compressed sensing,” Phys. Med. Biol., 61 8908 –8940 https://doi.org/10.1088/1361-6560/61/24/8908 PHMBA7 0031-9155 (2016). Google Scholar

106.

S. R. Arridge et al., “On the adjoint operator in photoacoustic tomography,” Inverse Probl., 32 115012 https://doi.org/10.1088/0266-5611/32/11/115012 INPEEY 0266-5611 (2016). Google Scholar

107.

Z. Belhachmi, T. Glatz and O. Scherzer, “A direct method for photoacoustic tomography with inhomogeneous sound speed,” Inverse Probl., 32 (4), 045005 https://doi.org/10.1088/0266-5611/32/4/045005 INPEEY 0266-5611 (2016). Google Scholar

108.

M. Bergounioux et al., “An optimal control problem in photoacoustic tomography,” Math. Models Methods Appl. Sci., 24 (12), 2525 –2548 https://doi.org/10.1142/S0218202514500286 MMMSEU 0218-2025 (2014). Google Scholar

109.

X. L. Deán-Ben et al., “Accurate model-based reconstruction algorithm for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imaging, 31 (10), 1922 –1928 https://doi.org/10.1109/TMI.2012.2208471 ITMID4 0278-0062 (2012). Google Scholar

110.

M. Haltmeier, R. Kowar and L. V. Nguyen, “Iterative methods for photoacoustic tomography in attenuating acoustic media,” Inverse Probl., 33 115009 https://doi.org/10.1088/1361-6420/aa8cba INPEEY 0266-5611 (2017). Google Scholar

111.

A. Javaherian and S. Holman, “A multi-grid iterative method for photoacoustic tomography,” IEEE Trans. Med. Imaging, 36 (3), 696 –706 https://doi.org/10.1109/TMI.2016.2625272 ITMID4 0278-0062 (2017). Google Scholar

112.

J. Prakash et al., “Basis pursuit deconvolution for improving model-based reconstructed images in photoacoustic tomography,” Biomed. Opt. Express, 5 (5), 1363 –1377 https://doi.org/10.1364/BOE.5.001363 BOEICL 2156-7085 (2014). Google Scholar

113.

K. Wang et al., “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Biol., 57 5399 –5423 https://doi.org/10.1088/0031-9155/57/17/5399 PHMBA7 0031-9155 (2012). Google Scholar

114.

D. Razansky et al., “Maximum entropy based non-negative optoacoustic tomographic image reconstruction,” IEEE Trans. Biomed. Eng., 66 (9), 2604 –2616 https://doi.org/10.1109/TBME.2019.2892842 IEBEAX 0018-9294 (2019). Google Scholar

115.

J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, Springer, New York (2005). Google Scholar

116.

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

117.

J. Tick et al., “Three dimensional photoacoustic tomography in Bayesian framework,” J. Acoust. Soc. Am., 144 (4), 2061 –2071 https://doi.org/10.1121/1.5057109 JASMAN 0001-4966 (2018). Google Scholar

118.

J. Tick, A. Pulkkinen and T. Tarvainen, “Image reconstruction with uncertainty quantification in photoacoustic tomography,” J. Acoust. Soc. Am., 139 1951 –1961 https://doi.org/10.1121/1.4945990 JASMAN 0001-4966 (2016). Google Scholar

119.

T. Sahlström et al., “Modeling of errors due to uncertainties in ultrasound sensor locations in photoacoustic tomography,” IEEE Trans. Med. Imaging, 39 (6), 2140 –2150 https://doi.org/10.1109/TMI.2020.2966297 ITMID4 0278-0062 (2020). Google Scholar

120.

J. Tick, A. Pulkkinen and T. Tarvainen, “Modelling of errors due to speed of sound variations in photoacoustic tomography using a Bayesian framework,” Biomed. Phys. Eng. Express, 6 (1), 015003 https://doi.org/10.1088/2057-1976/ab57d1 (2019). Google Scholar

121.

J. Gröhl et al., “Deep learning for biomedical photoacoustic imaging: a review,” Photoacoutics, 22 100241 https://doi.org/10.1016/j.pacs.2021.100241 (2021). Google Scholar

122.

A. Hauptmann and B. Cox, “Deep learning in photoacoustic tomography: current approaches and future directions,” J. Biomed. Opt., 25 112903 https://doi.org/10.1117/1.JBO.25.11.112903 JBOPFO 1083-3668 (2020). Google Scholar

123.

C. Yang et al., “Review of deep learning for photoacoustic imaging,” Photoacoustics, 21 100215 https://doi.org/10.1016/j.pacs.2020.100215 (2020). Google Scholar

124.

H. Lan et al., “A deep learning approach to reconstruct the photoacoustic image using multi-frequency data,” in IEEE Int. Ultrason. Symp., (2019). https://doi.org/10.1109/ULTSYM.2019.8926287 Google Scholar

125.

D. Waibel et al., “Reconstruction of initial pressure from limited view photoacoustic images using deep learning,” Proc. SPIE, 10494 104942S https://doi.org/10.1117/12.2288353 PSISDG 0277-786X (2018). Google Scholar

126.

D. Allman, A. Reiter and M. A. L. Bell, “Photoacoustic source detection and reflection artifact removal enabled by deep learning,” IEEE Trans. Med. Imaging, 37 1464 –1477 https://doi.org/10.1109/TMI.2018.2829662 ITMID4 0278-0062 (2018). Google Scholar

127.

N. Davoudi, X. L. Deán-Ben and D. Razansky, “Deep learning optoacoustic tomography with sparse data,” Nat. Mach. Intell., 1 453 –460 https://doi.org/10.1038/s42256-019-0095-3 (2019). Google Scholar

128.

S. Antholzer, M. Haltmeier and J. Schwab, “Deep learning for photoacoustic tomography from sparse data,” Inverse Probl. Sci. Eng., 27 987 –1005 https://doi.org/10.1080/17415977.2018.1518444 (2018). Google Scholar

129.

S. Guan et al., “Fully dense UNet for 2-D sparse photoacoustic tomography artifact removal,” IEEE J. Biomed. Health Inf., 24 568 –576 https://doi.org/10.1109/JBHI.2019.2912935 (2020). Google Scholar

130.

S. Antholzer et al., “NETT regularization for compressed sensing photoacoustic tomography,” Proc. SPIE, 10878 108783B https://doi.org/10.1117/12.2508486 PSISDG 0277-786X (2019). Google Scholar

131.

Y. E. Boink, S. Manohar and C. Brune, “A partially-learned algorithm for joint photo-acoustic reconstruction and segmentation,” IEEE Trans. Med. Imaging, 39 129 –139 https://doi.org/10.1109/TMI.2019.2922026 ITMID4 0278-0062 (2019). Google Scholar

132.

H. Goh et al., “Solving Bayesian inverse problems via variational autoencoders,” in Proc. Mach. Learn. Res., 386 –425 (2021). Google Scholar

133.

A. Hauptmann et al., “Model based learning for accelerated, limited-view 3-D photoacoustic tomography,” IEEE Trans. Med. Imaging, 37 (6), 1382 –1393 https://doi.org/10.1109/TMI.2018.2820382 ITMID4 0278-0062 (2018). Google Scholar

134.

K.-T. Hsu, S. Guan and P. V. Chitnis, “Fast iterative reconstruction for photoacoustic tomography using learned physical model: theoretical validation,” Photoacoustics, 29 100452 https://doi.org/10.1016/j.pacs.2023.100452 (2023). Google Scholar

135.

S. Lunz et al., “On learned operator correction in inverse problems,” SIAM J. Imaging Sci., 14 92 –127 https://doi.org/10.1137/20M1338460 (2021). Google Scholar

136.

T. Sahlström and T. Tarvainen, “Utilizing variational autoencoders in the bayesian inverse problem of photoacoustic tomography,” SIAM J. Imaging Sci., 16 (1), 89 –110 https://doi.org/10.1137/22M1489897 (2023). Google Scholar

137.

H. Shan et al., “Simultaneous reconstruction of the initial pressure and sound speed in photoacoustic tomography using a deep-learning approach,” Proc. SPIE, 11105 1110504 https://doi.org/10.1117/12.2529984 PSISDG 0277-786X (2019). Google Scholar

138.

P. Kuchment and L. Kunyansky, Mathematics of Photoacoustic and Thermoacoustic Tomography, 817 –865 Springer, New York (2011). Google Scholar

139.

J. Poudel, Y. Lou and M. A. Anastasio, “A survey of computational frameworks for solving the acoustic inverse problem in three-dimensional photoacoustic computed tomography,” Phys. Med. Biol., 64 14TR01 https://doi.org/10.1088/1361-6560/ab2017 PHMBA7 0031-9155 (2019). Google Scholar

140.

K. Wang, M. A. Anastasio, “Photoacoustic and thermoacoustic tomography: image formation principles,” Handbook of Mathematical Methods in Imaging, 1081 –1116 Springer, New York (2015). Google Scholar

141.

G. Bal and K. Ren, “Multi-source quantitative photoacoustic tomography in a diffusive regime,” Inverse Probl., 27 075003 https://doi.org/10.1088/0266-5611/27/7/075003 INPEEY 0266-5611 (2011). Google Scholar

142.

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

143.

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

144.

Z. Wang, W. Tao and H. Zhao, “The optical inverse problem in quantitative photoacoustic tomography: a review,” Photonics, 10 487 https://doi.org/10.3390/photonics10050487 (2023). Google Scholar

145.

X. Zhou et al., “Evaluation of fluence correction algorithms in multispectral photoacoustic imaging,” Photoacoustics, 19 100181 https://doi.org/10.1016/j.pacs.2020.100181 (2020). Google Scholar

146.

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

147.

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

148.

D. Razansky, “Multispectral optoacoustic tomography—volumetric color hearing in real time,” IEEE J. Sel. Top. Quantum Electron., 18 (3), 1234 –1243 https://doi.org/10.1109/JSTQE.2011.2172192 IJSQEN 1077-260X (2012). Google Scholar

149.

D. Razansky, J. Baeten and V. Ntziachristos, “Sensitivity of molecular target detection by multispectral optoacoustic tomography (MSOT),” Med. Phys., 36 (3), 939 –945 https://doi.org/10.1118/1.3077120 MPHYA6 0094-2405 (2009). Google Scholar

150.

P. Shao, B. Cox and R. J. Zemp, “Estimating optical absorption, scattering and Grueneisen distributions with multiple-illumination photoacoustic tomography,” Appl. Opt., 50 (19), 3145 –3154 https://doi.org/10.1364/AO.50.003145 APOPAI 0003-6935 (2011). Google Scholar

151.

H. Ammari et al., “Reconstruction of the optical absorption coefficient of a small absorber from the absorbed energy density,” SIAM J. Appl. Math., 71 (3), 676 –693 https://doi.org/10.1137/09077905X SMJMAP 0036-1399 (2011). Google Scholar

152.

B. Banerjee et al., “Quantitative photoacoustic tomography from boundary pressure measurements: noniterative recovery of optical absorption coefficient from the reconstructed absorbed energy map,” J. Opt. Soc. Am. A, 25 (9), 2347 –2356 https://doi.org/10.1364/JOSAA.25.002347 JOAOD6 0740-3232 (2008). Google Scholar

153.

B. T. Cox et al., “Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method,” Appl. Opt., 45 1866 –1875 https://doi.org/10.1364/AO.45.001866 APOPAI 0003-6935 (2006). Google Scholar

154.

T. Jetzfellner et al., “Performance of iterative optoacoustic tomography with experimental data,” Appl. Phys. Lett., 95 013703 https://doi.org/10.1063/1.3167280 APPLAB 0003-6951 (2009). Google Scholar

155.

J. Ripoll and V. Ntziachristos, “Quantitative point source photoacoustic inversion formulas for scattering and absorbing media,” Phys. Rev. E, 71 031912 https://doi.org/10.1103/PhysRevE.71.031912 (2005). Google Scholar

156.

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

157.

A. Pulkkinen et al., “Approximate marginalization of unknown scattering in quantitative photoacoustic tomography,” Inverse Probl. Imaging, 8 (3), 811 –829 https://doi.org/10.3934/ipi.2014.8.811 (2014). Google Scholar

158.

F. M. Brochu et al., “Towards quantitative evaluation of tissue absorption coefficients using light fluence correction in optoacoustic tomography,” IEEE Trans. Med. Imaging, 36 (1), 322 –331 https://doi.org/10.1109/TMI.2016.2607199 ITMID4 0278-0062 (2016). Google Scholar

159.

Z. Liang et al., “Automatic 3-D segmentation and volumetric light fluence correction for photoacoustic tomography based on optimal 3-D graph search,” Med. Image Anal., 75 102275 https://doi.org/10.1016/j.media.2021.102275 (2022). Google Scholar

160.

W. Naetar and O. Scherzer, “Quantitative photoacoustic tomography with piecewise constant material parameters,” SIAM J. Imaging Sci., 7 (3), 1755 –1774 https://doi.org/10.1137/140959705 (2014). Google Scholar

161.

S. Zhang et al., “MRI information-based correction and restoration of photoacoustic tomography,” IEEE Trans. Med. Imaging, 41 (9), 2543 –2555 https://doi.org/10.1109/TMI.2022.3165839 ITMID4 0278-0062 (2022). Google Scholar

162.

F. Asllanaj and A. Addoum, “Simultaneous reconstruction of absorption, scattering and anisotropy factor distributions in quantitative photoacoustic tomography,” Biomed. Phys. Eng. Express, 6 (4), 045010 https://doi.org/10.1088/2057-1976/ab90a0 (2020). Google Scholar

163.

E. Bonnetier, M. Choulli and F. Triki, “Stability for quantitative photoacoustic tomography revisited,” Res. Math. Sci., 9 (2), 24 https://doi.org/10.1007/s40687-022-00322-6 (2022). Google Scholar

164.

H. Gao, H. Zhao and S. Osher, Bregman methods in quantitative photoacoustic tomography (2010). Google Scholar

165.

A. Oraevsky et al., “Full-view 3D imaging system for functional and anatomical screening of the breast,” Proc. SPIE, 10494 104942Y https://doi.org/10.1117/12.2318802 PSISDG 0277-786X (2018). Google Scholar

166.

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

167.

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

168.

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

169.

A. Pulkkinen et al., “Quantitative photoacoustic tomography using illuminations from a single direction,” J. Biomed. Opt., 20 (3), 036015 https://doi.org/10.1117/1.JBO.20.3.036015 JBOPFO 1083-3668 (2015). Google Scholar

170.

A. Q. Bauer et al., “Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography,” J. Biomed. Opt., 16 (9), 096016 https://doi.org/10.1117/1.3626212 JBOPFO 1083-3668 (2011). Google Scholar

171.

X. Li and H. Jiang, “Impact of inhomogeneous optical scattering coefficient distribution on recovery of optical absorption coefficient maps using tomographic photoacoustic data,” Phys. Med. Biol., 58 999 –1011 https://doi.org/10.1088/0031-9155/58/4/999 PHMBA7 0031-9155 (2013). Google Scholar

172.

O. Nykänen, A. Pulkkinen and T. Tarvainen, “Quantitative photoacoustic tomography augmented with surface light measurements,” Biomed. Opt. Express, 8 (10), 4380 –4395 https://doi.org/10.1364/BOE.8.004380 BOEICL 2156-7085 (2017). Google Scholar

173.

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

174.

A. Hussain et al., “Photoacoustic and acousto-optic tomography for quantitative and functional imaging,” Optica, 5 (12), 1579 –1589 https://doi.org/10.1364/OPTICA.5.001579 (2018). Google Scholar

175.

Y. Zhou et al., “Deep-tissue temperature mapping by multi-illumination photoacoustic tomography aided by a diffusion optical model: a numerical study,” J. Biomed. Opt., 23 (1), 016014 https://doi.org/10.1117/1.JBO.23.1.016014 JBOPFO 1083-3668 (2018). Google Scholar

176.

S. Zhang et al., “Pixel-wise reconstruction of tissue absorption coefficients in photoacoustic tomography using a non-segmentation iterative method,” Photoacoustics, 28 100390 https://doi.org/10.1016/j.pacs.2022.100390 (2022). Google Scholar

177.

M. Haltmeier, L. Neumann and S. Rabanser, “Single-stage reconstruction algorithm for quantitative photoacoustic tomography,” Inverse Probl., 31 065005 https://doi.org/10.1088/0266-5611/31/6/065005 INPEEY 0266-5611 (2015). Google Scholar

178.

O. Scherzer et al., Variational Methods in Imaging, 167 Springer( (2009). Google Scholar

179.

A. Hannukainen et al., “Efficient inclusion of total variation type priors in quantitative photoacoustic tomography,” SIAM J. Imaging Sci., 9 (3), 1132 –1153 https://doi.org/10.1137/15M1051737 (2016). Google Scholar

180.

X. Zhang et al., “Forward-backward splitting method for quantitative photoacoustic tomography,” Inverse Probl., 30 125012 https://doi.org/10.1088/0266-5611/30/12/125012 INPEEY 0266-5611 (2014). Google Scholar

181.

J. Nocedal and S. J. Wright, Numerical Optimization, Springer( (1999). Google Scholar

182.

G. S. Alberti and H. Ammari, “Disjoint sparsity for signal separation and applications to hybrid inverse problems in medical imaging,” Appl. Comput. Harmon Anal., 42 (2), 319 –349 https://doi.org/10.1016/j.acha.2015.08.013 (2017). Google Scholar

183.

A. Rosenthal, D. Razansky and V. Ntziachristos, “Quantitative optoacoustic signal extraction using sparse signal representation,” IEEE Trans. Med. Imaging, 28 (12), 1997 –2006 https://doi.org/10.1109/TMI.2009.2027116 ITMID4 0278-0062 (2009). Google Scholar

184.

J. M. Bardsley, Computational Uncertainty Quantification for Inverse Problems, SIAM( (2018). Google Scholar

185.

H. Deng et al., “Deep learning in photoacoustic imaging: a review,” J. Biomed. Opt., 26 (4), 040901 https://doi.org/10.1117/1.JBO.26.4.040901 JBOPFO 1083-3668 (2021). Google Scholar

186.

P. Rajendran, A. Sharma and M. Pramanik, “Photoacoustic imaging aided with deep learning: a review,” Biomed. Eng. Lett., 12 155 –173 https://doi.org/10.1007/s13534-021-00210-y (2022). Google Scholar

187.

C. Yang et al., “Review of deep learning for photoacoustic imaging,” Photoacoustics, 21 100215 https://doi.org/10.1016/j.pacs.2020.100215 (2021). Google Scholar

188.

C. Bench, A. Hauptmann and B. Cox, “Toward accurate quantitative photoacoustic imaging: learning vascular blood oxygen saturation in three dimensions,” J. Biomed. Opt., 25 (8), 085003 https://doi.org/10.1117/1.JBO.25.8.085003 JBOPFO 1083-3668 (2020). Google Scholar

189.

C. Cai et al., “End-to-end deep neural network for optical inversion in quantitative photoacoustic imaging,” Opt. Lett., 43 (12), 2752 –2755 https://doi.org/10.1364/OL.43.002752 OPLEDP 0146-9592 (2018). Google Scholar

190.

T. Chen et al., “A deep learning method based on U-Net for quantitative photoacoustic imaging,” Proc. SPIE, 11240 112403V https://doi.org/10.1117/12.2543173 PSISDG 0277-786X (2020). Google Scholar

191.

J. Gröhl et al., “Moving beyond simulation: data-driven quantitative photoacoustic imaging using tissue-mimicking phantoms,” (2023). Google Scholar

192.

J. Gröhl et al., “Learned spectral decoloring enables photoacoustic oximetry,” Sci. Rep., 11 (1), 6565 https://doi.org/10.1038/s41598-021-83405-8 SRCEC3 2045-2322 (2021). Google Scholar

193.

T. Kirchner and M. Frenz, “Multiple illumination learned spectral decoloring for quantitative optoacoustic oximetry imaging,” J. Biomed. Opt., 26 (8), 085001 https://doi.org/10.1117/1.JBO.26.8.085001 JBOPFO 1083-3668 (2021). Google Scholar

194.

G. P. Luke et al., “O-Net: a convolutional neural network for quantitative photoacoustic image segmentation and oximetry,” (2019). Google Scholar

195.

A. Madasamy et al., “Deep learning methods hold promise for light fluence compensation in three-dimensional optoacoustic imaging,” J. Biomed. Opt., 27 (10), 106004 https://doi.org/10.1117/1.JBO.27.10.106004 JBOPFO 1083-3668 (2022). Google Scholar

196.

C. Yang and F. Gao, “EDA-Net: dense aggregation of deep and shallow information achieves quantitative photoacoustic blood oxygenation imaging deep in human breast,” Lect. Notes Comput. Sci., 11764 246 –254 https://doi.org/10.1007/978-3-030-32239-7_28 LNCSD9 0302-9743 (2019). Google Scholar

197.

Y. Lou et al., “Generation of anatomically realistic numerical phantoms for photoacoustic and ultrasonic breast imaging,” J. Biomed. Opt., 22 (4), 041015 https://doi.org/10.1117/1.JBO.22.4.041015 JBOPFO 1083-3668 (2017). Google Scholar

198.

S. Park et al., “Stochastic three-dimensional numerical phantoms to enable computational studies in quantitative optoacoustic computed tomography of breast cancer,” J. Biomed. Opt., 28 (6), 066002 https://doi.org/10.1117/1.JBO.28.6.066002 JBOPFO 1083-3668 (2023). Google Scholar

199.

C. Bench, “Data-driven quantitative photoacoustic tomography,” University College London, (2022). Google Scholar

200.

J. Li et al., “Deep learning-based quantitative optoacoustic tomography of deep tissues in the absence of labeled experimental data,” Optica, 9 (1), 32 –41 https://doi.org/10.1364/OPTICA.438502 (2022). Google Scholar

201.

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

202.

S. Okawa et al., “3D quantitative photoacoustic image reconstruction using Monte Carlo method and linearization,” Proc. SPIE, 10494 104944Y https://doi.org/10.1117/12.2293229 PSISDG 0277-786X (2018). Google Scholar

203.

A. Capart et al., “Quantitative photoacoustic reconstruction of the optical properties of intervertebral discs using a gradient descent scheme,” Photonics, 9 116 https://doi.org/10.3390/photonics9020116 (2022). Google Scholar

204.

J. Buchmann et al., “Quantitative PA tomography of high resolution 3-D images: experimental validation in a tissue phantom,” Photoacoustics, 17 100157 https://doi.org/10.1016/j.pacs.2019.100157 (2020). Google Scholar

205.

C. M. Macdonald, S. Arridge and S. Powell, “Efficient inversion strategies for estimating optical properties with Monte Carlo radiative transport models,” J. Biomed. Opt., 25 (8), 085002 https://doi.org/10.1117/1.JBO.25.8.085002 JBOPFO 1083-3668 (2020). Google Scholar

206.

N. Hänninen et al., “Estimating absorption and scattering in quantitative photoacoustic tomography with an adaptive Monte Carlo method for light transport,” Inverse Probl. Imaging, (). Google Scholar

207.

P. Shao, T. Harrison and R. J. Zemp, “Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data,” Biomed. Opt. Express, 3 (12), 3240 –3248 https://doi.org/10.1364/BOE.3.003240 BOEICL 2156-7085 (2012). Google Scholar

208.

N. Song, C. Deumié and A. Da Silva, “Considering sources and detectors distributions for quantitative photoacoustic tomography,” Biomed. Opt. Express, 5 (11), 3960 –3974 https://doi.org/10.1364/BOE.5.003960 BOEICL 2156-7085 (2014). Google Scholar

209.

H. Gao, J. Feng and L. Song, “Limited-view multi-source quantitative photoacoustic tomography,” Inverse Probl., 31 065004 https://doi.org/10.1088/0266-5611/31/6/065004 INPEEY 0266-5611 (2015). Google Scholar

210.

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

211.

A. Pulkkinen et al., “Direct estimation of optical parameters from photoacoustic time series in quantitative photoacoustic tomography,” IEEE Trans. Med. Imaging, 35 (11), 2497 –2508 https://doi.org/10.1109/TMI.2016.2581211 ITMID4 0278-0062 (2016). Google Scholar

212.

M. Venugopal et al., “Quantitative photoacoustic tomography by stochastic search: direct recovery of the optical absorption field,” Opt. Lett., 41 (18), 4202 –4205 https://doi.org/10.1364/OL.41.004202 OPLEDP 0146-9592 (2016). Google Scholar

213.

Y. Xu et al., “Reconstructions in limited-view thermoacoustic tomography,” Med. Phys., 31 (4), 724 –733 https://doi.org/10.1118/1.1644531 MPHYA6 0094-2405 (2004). Google Scholar

Biographies of the authors are not available.

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

Tanja Tarvainen and Ben Cox "Quantitative photoacoustic tomography: modeling and inverse problems," Journal of Biomedical Optics 29(S1), S11509 (20 December 2023). https://doi.org/10.1117/1.JBO.29.S1.S11509

Received: 20 September 2023; Accepted: 28 November 2023; Published: 20 December 2023

Access the abstract

JOURNAL ARTICLE
25 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY