Quantitative spectroscopic photoacoustic imaging: a review

Benjamin T. Cox; Jan G. Laufer; Paul C. Beard; Simon R. Arridge

doi:10.1117/1.JBO.17.6.061202

1 June 2012 Quantitative spectroscopic photoacoustic imaging: a review

Benjamin T. Cox, Jan G. Laufer, Paul C. Beard, Simon R. Arridge

Author Affiliations +

Journal of Biomedical Optics, Vol. 17, Issue 6, 061202 (June 2012). https://doi.org/10.1117/1.JBO.17.6.061202

Abstract

Obtaining absolute chromophore concentrations from photoacoustic images obtained at multiple wavelengths is a nontrivial aspect of photoacoustic imaging but is essential for accurate functional and molecular imaging. This topic, known as quantitative photoacoustic imaging, is reviewed here. The inverse problems involved are described, their nature (nonlinear and ill-posed) is discussed, proposed solution techniques and their limitations are explained, and the remaining unsolved challenges are introduced.

1. Introduction

The prospect of a soft tissue imaging modality that can achieve fine spatial resolution, high sensitivity, and good specificity has led, in the last decade, to a rapid growth of interest in photoacoustic (PA) imaging.¹^–³ PA imaging shares one of the distinctive advantages of other optical imaging techniques in that it can be used spectroscopically; measurements made at multiple optical wavelengths can be used to provide information related to molecular composition. The goal of quantitative PA imaging (QPAI) is to exploit this advantage by converting multiwavelength sets of PA images into images showing quantitatively accurate estimates of the spatially varying concentrations of chromophores (light-absorbing substances) embedded within an optically scattering medium such as biological tissue.

1.1.

Scope of the Paper

The aims in writing this paper have been twofold: to describe the challenges in QPAI, and to group the attempts that have been made to tackle QPAI together with similar methods, or with methods that make similar assumptions, in a coherent way that highlights the similarities and differences. It has neither been our intention to present a chronologically accurate history of quantitative PAT, nor to provide a reference to every paper in this area, but rather to give sufficient references to act as an introduction to the literature.

1.1.1.

What does “quantitative photoacoustic imaging” cover?

There are a number of quantities that might be determined by QPAI. The contrast in PA imaging is due to chromophores, and the most fundamental quantity to be determined is their concentrations. Typical chromophores are (1) endogenous molecules such as oxyhemoglobin or deoxyhemoglobin, melanin, lipids, and water; (2) exogenous contrast agents, sometimes targeted to a particular cell-surface receptor, biomolecule, or organ; and (3) absorbing enzymes or other proteins resulting from the expression of reporter genes linked to the expression of a gene of interest. QPAI might also be used to obtain quantitative estimates of parameters of physiological interest derived from chromophore concentrations, such as blood oxygenation, the spatially varying optical absorption coefficient (the sum of contributions from all the chromophores), and the optical scattering coefficient, which is related to the tissue microstructure. Which of these is considered the primary quantity will depend on the application.

1.1.2.

What is not covered in this paper?

There are essentially two aspects to PA imaging: acoustic and optical. The focus of this paper is on the optical inversion, and the acoustic inverse problem (PA image reconstruction) is not covered here. It has been studied extensively, and we refer the interested reader to Refs. 4 5.6.7.8.9.–10 and to the references contained therein. This review does not discuss the recovery of properties of the tissue other than optical properties (e.g., sound speed, density, and acoustic attenuation estimation) although the accurate determination of these acoustic quantities could have an indirect effect on QPAI if they are used to improve the accuracy of the acoustic image reconstruction. A great deal of work has been done in recent years on finding and designing suitable contrast agents for PAI, but the focus here is on techniques to facilitate the estimation of their concentrations rather than on the chemistry of the substances themselves. Techniques to determine the optimal choice of wavelengths for QPAI is outside the scope of this article as it will be very dependent on the particular arrangement (the object, imaging mode, main source of contrast, other chromophores present, etc.). This question has been tackled for some specific scenarios.¹¹^–¹³ Finally, photoacoustic techniques that have no imaging aspect are not covered here, so “gas phase” PA spectroscopy using a PA cell, which is a mature field in its own right,¹⁴ is not described.

1.1.3.

Layout of the paper

The layout of the paper is as follows: Sec. 1.2 gives a brief introduction to PA imaging in general; Sec. 2.1 introduces the quantitative inverse problem; and Secs. 2.2 and 2.3 describe models of light transport and, in broad terms, how they might be used in inversions for estimating optical properties. The later sections discuss the specifics of methods that have been proposed assuming, respectively, that the situation is one-dimensional (Sec. 3), that optical properties can be obtained from acoustic measurements made with a single detector rather than an array (Sec. 4), and that the situation is fully three-dimensional (Sec. 5). A discussion, including PA efficiency, and a summary conclude the paper.

1.2.

Photoacoustic Imaging

This section will start with a brief description of the photoacoustic effect: how PA waves are generated. The term photoacoustic imaging is used to describe a number of related imaging modes that exploit this effect to image objects with heterogeneous optical absorption. As these PA images are a prerequisite for QPAI, Secs. 1.2.2 and 1.2.3 give very short summaries of the main imaging modalities currently in use.

1.2.1.

The photoacoustic effect

There are several mechanisms by which light can be used to generate sound waves.¹⁵^–¹⁷ PA imaging usually uses light in the non-ionizing visible or near-infrared (NIR) parts of the spectrum because the NIR “window” (a range of wavelengths over which both water absorption and tissue scattering are low) permits deeper light penetration and so greater imaging depth. For these wavelength ranges, and for light intensities below exposure safety limits, heat deposition is the dominant mechanism for the generation of acoustic pulses. The photoacoustic effect, as the term will be used here, therefore refers to the generation of sound waves through the absorption of light and conversion to heat.

Three steps are involved in the PA effect: (1) the absorption of a photon, (2) the thermalization of the absorbed energy and a corresponding localized pressure increase, and (3) the propagation of this pressure perturbation as an acoustic wave due to the elastic nature of tissue. The absorption of a photon, which typically takes place on a femtosecond timescale, raises a molecule to an excited state. There are then a number of possible subsequent chains of events for photons in the visible and NIR. The two most important of these are radiative decay and thermalization. In PA applications it is usually assumed that thermalization (nonradiative decay) dominates and that therefore all the photon energy is converted into heat via $vibrational / collisional$ relaxation on a sub-nanosecond timescale. This localized injection of heat will lead to a small rise in the local temperature and a related rise in the local pressure, $p_{0}$ . If the heat per unit volume deposited in the tissue (the absorbed energy density) is written as $H$ and the usual assumption is made that the pressure rise is linearly related to the absorbed energy, the pressure increase may be written as

Eq. (1)

p_{0} = \hat{Γ} H,

where

\hat{Γ}

is the PA efficiency. (This efficiency can be identified with the Grüneisen parameter

Γ

for an absorbing fluid; Sec. 6.3). It is this increase in pressure that subsequently becomes a propagating acoustic (ultrasonic) wave, because of the elastic nature of tissue. It is therefore referred to as the initial acoustic pressure distribution,

p_{0}

. The absorbed energy density can be written as

Eq. (2)

H = μ_{a} Φ,

where

μ_{a}

is the optical absorption coefficient (the total absorption due to the chromophores) and

Φ

is the light fluence (the radiance integrated over all directions and time; Sec. 2.2). The initial acoustic pressure distribution may therefore be written as the product of three quantities:

Eq. (3)

p_{0} = \hat{Γ} μ_{a} Φ .

The direct problem, from absorber (endogenous chromophore or contrast agent) to measured PA signal, is presented schematically in Fig. 1. It can be separated into two parts: optical propagation leading to a fluence distribution, and acoustic propagation leading to the detected ultrasonic pulse, coupled by the thermalization of the light energy.

Fig. 1

PA signal generation. Spatially varying chromophore concentrations (naturally occurring or contrast agents) give rise to optical absorption in the medium. The absorption and scattering coefficients, $μ_{a}$ and $μ_{s}$ , determine the fluence distribution $Φ$ (how the light from a source becomes distributed in the tissue), and hence the absorbed energy distribution $H$ . This energy generates a pressure distribution $p_{0}$ via thermalization, which because of the elastic nature of tissue then propagates as an acoustic wave. The resulting pulse is detected by a sensor resulting in the measured PA time series $p (t)$ .

1.2.2.

Photoacoustic tomography

PA tomography (PAT), sometimes called PA computed tomography, uses a widefield pulse of light to illuminate the tissue so that the whole tissue is flooded with light.¹⁸ An array of (ideally) omnidirectional point ultrasound detectors (or an array synthesized from many single-point measurements) is used to record the resulting PA waves. For the highest quality images, the array should enclose the object being imaged, although in practice good images are attainable when this is not the case (when a planar detection surface is used, for instance). An image reconstruction algorithm is used to reconstruct the initial acoustic pressure distribution from the recorded pressure time series.⁴^–¹⁰ This estimate of the initial acoustic pressure distribution is called a PAT image. It is inherently three dimensional. The spatial resolution is limited by the detector aperture used to measure the data, by the directionality and spacing of the detector elements, and by the maximum acoustic frequency detected. Spatial resolution is therefore higher closer to the detector because acoustic absorption (which absorbs the higher frequencies more) removes the higher frequencies for signals arriving from further away. Images with spatial resolutions of a few tens of microns for sub-centimeter depths, extending to resolutions of a millimeter or so at several centimeters, can readily be achieved. Variations on this theme include the use of integrating detectors in place of point receivers,¹⁹ and of detectors focused in a plane,²⁰ both of which reduce the image reconstruction problem to two dimensions.

1.2.3.

Photoacoustic scanning microscopy

In photoacoustic scanning microscopy (PAM) an ultrasonic detector is scanned over the sample and the A-lines (the PA time series) recorded at each scan point, which approximate depth profiles through the initial pressure distribution, are stacked up to form a three-dimensional (3-D) image. This requires no image reconstruction algorithm, in contrast to PAT. There are essentially two modes of PA microscopy, AR-PAM and OR-PAM. In acoustic resolution PA microscopy (AR-PAM), a focused ultrasound detector is used to record the PA signal, and the axial and lateral spatial resolutions are limited by ultrasonic considerations. In a typical implementation, the transducer is focused tightly to minimize the lateral resolution, and the tissue is illuminated by a ring of light sent around the transducer and weakly focused to the same point. For a high resolution system, the lateral spatial resolution is $\sim 45$ μm to a depth of a few millimeters, with the vertical resolution limited to half the shortest wavelength detectable ( $\sim 15$ μm).²¹ To improve the lateral spatial resolution to beyond the acoustic diffraction limit, optical resolution PA microscopy (OR-PAM) uses very tightly focused light to limit the illuminated region to a very small spot. As a tight focal spot can only be formed close to where the light enters the tissue (at deeper depths the scattering makes focusing very difficult) and as the PA signal is only generated in the illuminated region, OR-PAM is typically used for visualizing chromophores within the first 1 mm or so of tissue only. Sub-micron lateral resolutions and vertical resolutions limited to $\sim 10$ μm should be achievable.²²

2. Quantitative Photoacoustic Imaging

2.1.

Optical Inverse Problem

This section introduces the inverse problem that must be solved for quantitative PAI to be achieved. It is essentially the inversion of a light transport operator, hence it is called the optical inverse problem.

2.1.1.

Formulation of the problem

A set of PA images $p_{0} (x, λ)$ are obtained experimentally, each at a different optical wavelength $λ$ . (The position vector $x$ will in general be in three dimensions although it can be reduced to one or two under certain conditions.) If the system is calibrated and the PA efficiency is known, then the multiwavelength set of images of absorbed energy density $H (x, λ) = {\hat{Γ}}^{- 1} p_{0} (x, λ)$ can be taken as the measured data. The principal optical inversion in QPAI is a distributed parameter estimation problem that can then be stated as follows: find the concentrations $C_{k} (x)$ of $K$ chromophores with known molar absorption coefficient spectra $α_{k} (λ)$ , given the absorbed energy images $H (x, λ)$ when the chromophores are linked to the absorption coefficients with the linear mapping $L_{λ}$ :

Eq. (4)

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

and the absorption coefficients are linked to the absorbed energy images by the nonlinear mapping

T

:

Eq. (5)

H (x, λ) = T (μ_{a}, μ_{s}) = μ_{a} (x, λ) Φ (x, λ; μ_{a}, μ_{s}) .

The light fluence, $Φ (x, λ)$ is unknown and will depend on $μ_{a}$ as well as the optical scattering coefficient $μ_{s}$ . A fluence model is therefore required (Sec. 2.2). Equations (4) and (5) can be combined to give

Eq. (6)

H (x, λ) = T_{λ} (C_{k}, μ_{s}) = Φ (x, λ; C_{k}, μ_{s}) \sum_{k = 1}^{K} C_{k} (x) α_{k} (λ) .

Equations (4) and (5) suggest a two-stage inversion strategy: first recover the absorption coefficients $μ_{a} = T^{- 1} (H)$ and then find the concentrations $C_{k} = L_{λ}^{- 1} (μ_{a})$ , whereas Eq. (6) suggests a single inversion, $C_{k} = T_{λ}^{- 1} (H)$ . The former approach, of finding the absorption coefficients first, has been studied more than the latter, although there are potential advantages in both cases (see Sec. 6.1). No hard separation will be made between these two approaches in this paper, and wherever the absorption coefficient is treated as the unknown, it should be taken as read that a spectroscopic inversion would follow the recovery of the absorption coefficients. These two pathways are shown schematically in Fig. 2, which also shows a third approach: to combine the acoustic and optical inversions and invert directly for the chromophore concentrations from the pressure time series, $p (t)$ , that is, directly invert the mapping $W$ given by

Eq. (7)

p (t) = A \hat{Γ} T_{λ} (C_{k}) = W (C_{k}),

where

A

is the acoustic forward operator. Linearizations of

W

have been found useful in certain very specific circumstances (see Sec. 4), but in general the nonlinearity of the light

transport / absorption

model must be taken into account.

Fig. 2

The inverse problems in quantitative PA imaging. Solid lines indicate linear operators and dot-dash lines those that are inherently nonlinear. The acoustic pressure time series $p (t)$ are the measured data, and the chromophore concentrations $C_{k}$ are the unknowns. The concentrations may be obtained step by step: linear acoustic inversion, $A^{- 1}$ ; thermoelastic scaling, ${\hat{Γ}}^{- 1}$ ; nonlinear optical inversion for the optical coefficients, $T^{- 1}$ ; and finally , $L_{λ}^{- 1}$ a linear spectroscopic inversion of the absorption coefficient spectra to recover the chromophore concentrations. Alternatively, these last two stages may be combined into the operator $T_{λ}^{- 1}$ . A different approach is to attempt a one-step inversion of the whole process, $W^{- 1}$ . This may be linearized under specific and rather restricted conditions. The inverse light transport operators $T^{- 1}$ and its multiwavelength equivalent $T_{λ}^{- 1}$ are at the core of the methods described in this paper.

For OR-PAM images obtained by scanning the focused illumination spot across the region of interest, the fluence distribution $Φ (x)$ will be different for each measurement, so Eq. (5) should be replaced by $H (x_{i}) = μ_{a} (x_{i}) Φ_{i} [x_{i}, μ_{a} (x)]$ where $Φ_{i} (x)$ is the fluence distribution at point $x$ when the illumination is focused at the point $x_{i}$ . Full field inversions (Sec. 5) are not usually applied to OR-PAM images, but if they were, this detail would need to be considered.

2.1.2.

Spectral coloring and structural distortion

A cursory glance at Eq. (6) might suggest that there is a linear relationship between the chromophore concentrations $C_{k}$ and the absorbed energy density $H$ . If this were the case, then it would allow the chromophore concentrations $C_{k}$ to be recovered from measurements of the absorbed energy density spectra $H (λ)$ straightforwardly, for example, using a linear matrix inversion. Unfortunately the unknown fluence $Φ$ depends in a nontrivial way on the chromophore concentrations because $T_{λ}$ is both nonlinear and nonlocal ( $H$ at one point can be affected by $μ_{a}$ some distance away via its effect on the fluence). It is not reasonable, therefore, to assume that $H (λ) \propto μ_{a} (λ)$ because the spectrum of the absorbed energy density at a given point will be affected by the fluence, which will have been colored by its passage through the medium. Because biological tissue is highly scattering, each photon will take a long and convoluted path through the tissue. The spectra of the fluence and the absorption are therefore intertwined, and $H (λ)$ will be a complicated combination of the absorption spectra of the constituent chromophores at each point in the irradiated volume. (This coloring of the spectrum is the same effect as when seeing the world through colored glass.)

As the fluence $Φ (x, λ)$ is a function of position as well as wavelength, it will also have another effect. As well as coloring the spectrum, it will also distort the image structure. The absorbed energy density at a single wavelength is given by Eq. (2) or (4): $H (x) = μ_{a} (x) Φ (x)$ . The unknown and nonuniform light distribution $Φ (x)$ will therefore distort what would, in its absence, have been an image proportional to $μ_{a} (x)$ , the absorption coefficient as a function of position. This cannot be avoided by ensuring that the illumination of the target is uniform at the surface, because as soon as there is any absorption or scattering in the medium, the fluence will vary spatially. When the absorption or scattering is spatially varying (as is often the case in tissue), the unknown fluence also will be. This structural distortion and the spectral coloring are two manifestations of the same phenomenon, namely the effect of the fluence on the PA image.²³

2.1.3.

Using a model to recover the concentrations

The measured PA signals $p (t)$ or images $p_{0}$ can be related to the chromophore concentrations $C_{k}$ by the use of analytical or numerical models to approximate the physics encapsulated by the operators $T_{λ}$ , $T$ , or $W$ , which can then be inverted to estimate $C_{k}$ . Potential models will vary in their accuracy (how much of the physics they model), their range of validity (over what values of the known parameters and the concentrations they are valid), their complexity (how easy they are to solve) and their invertibility (how easily they can be inverted—reverse engineered—to calculate the concentrations given the measured data). For the model to be useful, it must strike a balance between accuracy (including range of validity) and simplicity (including invertibility). The two main models that have been used to describe light are summarized below.

2.2.

Fluence Models

The essential difficulty with QPAI is the unknown light fluence. Inverting $T_{λ}$ or $T$ to obtain $μ_{a}$ or $C_{k}$ would be trivial if the fluence were known, but in general it is not and so must be modeled. In some situations it may be possible to approximate the fluence distribution with an analytical expression or a simple formulation, for example, in a homogeneous or nonscattering medium. In the more general case, a numerical model of light propagation through the tissue is required. A number of numerical and mathematical models are relevant to light propagation in scattering media, and there is a considerable literature.²⁴^–²⁷ The requirements for a model are that it is sufficiently accurate to capture the essential characteristics of the light field, and perhaps, fast enough computationally to make its use in iterative inversion methods possible, where it may need to be evaluated many times. Two widely used light models are described briefly.

2.2.1.

Radiative transfer equation

Light is an electromagnetic wave satisfying Maxwell’s equations but treating it as such for propagation in turbid (highly scattering) media quickly reaches the limits of practical computation because of the spatial scales involved. It is more common, therefore, to use particle-based methods from transport theory to model the light distribution. The radiative transfer equation (which is Boltzmann’s transport equation applied to low energy, monochromatic, photons) is an integro-differential equation expressing the conservation of energy in the following form:

Eq. (8)

\frac{1}{c} \frac{\partial ϕ}{\partial t} (x, \hat{s}, t) = q (x, \hat{s}, t) - [\hat{s} \cdot \nabla + μ_{a} (x) + μ_{s} (x)] ϕ (x, \hat{s}, t) + μ_{s} \int Θ (\hat{s}, {\hat{s}}^{'}) ϕ (x, {\hat{s}}^{'}, t) d {\hat{s}}^{'},

where

ϕ (x, \hat{s}, t)

is the light radiance,

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

is called the scattering phase function and is the probability that a photon originally traveling in direction

\hat{s}

ends up traveling in direction

{\hat{s}}^{'}

if scattered,

q (x, \hat{s}, t)

is a source of photons, and

c

is the speed of light in the medium. The terms on the right-hand side of this equation account for the fact that the rate of change of the number of photons within a small region around the point

x

and traveling in direction

\hat{s}

could be due to (1) sources

q

; (2) net outflow of photons due to the radiance gradient,

\hat{s} \cdot \nabla

; (3) photons absorbed,

μ_{a} ϕ

; (4) photons scattered into another direction,

μ_{s} ϕ

; or (5) photons scattered into direction

\hat{s}

from another direction (given by the phase integral). Wave effects, polarization, radiative processes, ionization, inelastic scattering, and reactions are all neglected in this model.

Because in pulsed PA imaging the acoustic propagation occurs on a timescale several orders of magnitude longer than the heat deposition, the time-integrated absorbed power density (i.e., the absorbed energy density) is the quantity of interest, so only the time-independent radiative transfer equation (RTE) is required.

Eq. (9)

(\hat{s} \cdot \nabla + μ_{t}) ϕ (x, \hat{s}) - μ_{s} \int Θ (\hat{s}, {\hat{s}}^{'}) ϕ (x, {\hat{s}}^{'}) d {\hat{s}}^{'} = q (x, \hat{s}),

where the total attenuation coefficient

μ_{t} = μ_{a} + μ_{s}

, and

ϕ (x, \hat{s})

is now used to represent the time-integrated light radiance. The fluence

Φ

is the integral of the radiance over all angles

\hat{s}

:

Eq. (10)

Φ (x) = \int ϕ (x, {\hat{s}}^{'}) d {\hat{s}}^{'} .

To solve Eq. (9) computationally, it may, for example, be written in a weak formulation and implemented on a discretized mesh using the finite element method.²⁸^,²⁹

When there is no scattering, Eq. (9) reduces to $(\hat{s} \cdot \nabla + μ_{a}) ϕ = q$ . In the case of a collimated source propagating in the $z$ direction, the fluence, which equals the radiance as there is only one direction of propagation, satisfies the equation $\partial Φ / \partial z = - μ_{a} Φ$ . If the source illuminates the surface with a fluence of $Φ_{0}$ and $μ_{a}$ is constant, then the fluence as a function of depth may be written

Eq. (11)

Φ (z) = Φ_{0} \exp (- μ_{a} z),

which is known as Beer’s law or the Beer-Lambert-Bouguer law.

2.2.2.

Diffusion approximation

To obtain approximations to the RTE, the radiance can be written as a sum of spherical harmonics and truncated after $N$ terms. This leads to a family of approximations known as the $P_{N}$ approximations. When the directional dependence of the light distribution is weak, as within highly scattering tissue, it is often sufficient to take $N = 1$ . The time-independent $P_{1}$ approximation can be written as the two equations³⁰

Eq. (12)

μ_{a} Φ + \nabla \cdot F = q_{0}, \frac{F}{D} + \nabla Φ = 3 q_{1},

where

D = {[3 (μ_{a} + μ_{s}^{'})]}^{- 1}

is the optical diffusion coefficient,

μ_{s}^{'} = (1 - g) μ_{s}

is the reduced scattering coefficient (with

g

the anisotropy factor), and the vector flux,

F

, is defined as

Eq. (13)

F (x) = \int ϕ (x, {\hat{s}}^{'}) {\hat{s}}^{'} d {\hat{s}}^{'} .

The isotropic and mildly directional source terms $q_{o}$ and $q_{1}$ are defined analogously to Eqs. (10) and (13). When $q_{1} = 0$ , Eqs. (12) and (13) reduce to the diffusion approximation (DA):

Eq. (14)

μ_{a} Φ - \nabla \cdot (D \nabla) Φ = q_{0} .

The condition $μ_{s}^{'} ≫ μ_{a}$ is usually considered sufficient to ensure the accuracy of the diffusion approximation away from sources, because it ensures the scattered fluence is almost isotropic. For arbitrary heterogeneous media, Eq. (14) may be solved by standard numerical methods such as finite elements.³¹^,³² For a homogeneous medium, the solution to Eq. (14) is given by

Eq. (15)

Φ (x) = \int G_{0} (x - x^{'}) q_{0} (x^{'}) d x^{'}

with the (3-D) free space Green’s function,

G_{0}

[the solution to

(μ_{a} - D \nabla^{2}) G_{0} = δ (x - x^{'})

] given by

Eq. (16)

G_{0} (x, x^{'}) = \frac{\exp (- μ_{eff} | x - x^{'} |)}{4 π D | x - x^{'} |}, μ_{eff} = \sqrt{3 μ_{a} (μ_{a} + μ_{s}^{'})},

where

μ_{eff}

is called the effective attenuation coefficient. This solution will be used in some of the inversions in Sec. 5. In one dimension (1-D), the solution is

Φ (z) = Φ_{0} \exp (- μ_{eff} z)

where

Φ_{0}

is fluence incident on the surface, which is analogous to Eq. (11). Some care is required here because while this is the solution to Eq. (14) in 1-D, it only applies when the fluence everywhere is diffuse, which is rarely true close to the surface; where light enters the tissue, it often remains partially collimated. In practice, the fluence will often reach a peak at some distance beneath the surface, where both the incident and backscattered light are contributing to the fluence, before decaying exponentially. Introducing a scaling factor to account for this gives

Eq. (17)

Φ (z) \approx k Φ_{0} \exp (- μ_{eff} z) for z ≫ \frac{1}{μ_{t}},

which is used in the inversions in Sec. 3.3.

2.2.3.

Collimated source, $δ$ -Eddington approximation

For a collimated incident beam, the condition for the accuracy of the DA, $μ_{s}^{'} ≫ μ_{a}$ , must be supplemented with a second condition. As the collimated part of an incident beam will decay exponentially with depth $z$ as $\exp (- μ_{t} z)$ , the fluence will only be dominated by diffuse photons at depths $z ≫ 1 / μ_{t}$ . The optical coefficients in PA imaging applications are of the order of $μ_{a} = 0.1 {mm}^{- 1}$ , $μ_{s} = 10 {mm}^{- 1}$ , $g = 0.9$ , $μ_{s}^{'} = 1 {mm}^{- 1}$ so the condition that $μ_{s}^{'} ≫ μ_{a}$ is met, but the condition that the fluence is dominated by diffuse photons, only holds for distances $z ≫ 1 / μ_{t} \approx 0.1 mm$ , so perhaps for $z > 1 mm$ or so. This sub-millimeter surface region is often of interest in PA imaging, so it is important that the light model is accurate there. This could be achieved by using the RTE or a higher order $P_{N}$ approximation, but the simplicity of the DA is appealing. The $δ$ -Eddington approximation attempts to provide a more accurate model of the light fluence in this surface region without losing the simple form of the DA.³³^–³⁵

2.2.4.

Monte Carlo models

As an alternative to analytical models, the Monte Carlo method is a purely numerical approach. It simulates the random walk taken by “packets of energy” as they propagate one by one through the scattering medium, losing energy as they go and being scattered according to the probability as given by the phase function $Θ$ .³⁶^,³⁷ A large number of packets are required for the absorbed energy density to converge to a continuous solution, but because the path of every packet is independent, these methods are straightforwardly parallelizable. As the availability of large computing clusters and clusters of graphical processing units (GPUs) is increasing quickly, implementations are becoming significantly faster.³⁸ Monte Carlo methods are often considered the “gold standard” for modeling light transport in turbid media and are frequently used to validate other numerical models.

2.3.

Parameter Estimation

The estimation of the values of the parameters of a model (e.g., the coefficients in an equation) from measurements assumed to correspond to its solution is a problem that occurs in many fields. It has therefore has been widely studied, both for parameters that are single scalars as well as for distributed parameters. The aim in QPAI is to recover chromophore concentrations (or absorption coefficients) as a function of position, and so it is a distributed parameter estimation problem. This section gives a brief overview of a few of these methods that have been applied to quantitative PAI.

2.3.1.

Inversion techniques

Linearization

Any linear model can be written in matrix form, which allows the full repertoire of matrix inversion routines to be applied to the problem. A popular approach with nonlinear problems such as QPAI is therefore linearization: to linearize them with respect to the unknown parameters, for example, to expand them as a Taylor series about a known state (known background absorption and scattering for instance). When the unknown parameters vary little from these chosen values, then the linearized model is often a good approximation. For example, it might be argued that a small localized change in absorption will not change the fluence. Linearizations of the QPAI problem are discussed in Secs. 4.1.1, 5.1.1, 5.1.2, and 5.2.1 but are typically applicable only in very limited circumstances, such as for small perturbations of the unknowns.

Direct inversion

For light models that are available in analytical form rather than purely numerically, it might be possible to rearrange the equations in such a way as to find a closed form expression for unknown parameters in terms of the measured data, or a simple noniterative procedure for determining them. A simple example is Beer’s law, Eq. (11), which can be inverted straightforwardly to give an expression for the absorption coefficient: $μ_{a} = - \ln [Φ (z) / Φ_{0}] / z$ . There is no systematic way for finding direct inversions such as this, as each will be specific to the problem under study. They are of great interest both because of the insight into the inverse problem they provide, and because—without the need for matrix inversions or iterations—they can be fast. A direct method proposed for QPAI is described in Sec. 5.2.2. As with any inversion technique, a direct inversion needs to be stable and robust to noise in order to be useful.

Fixed-point iteration

When the model equations can be rewritten into a form such that the unknown parameter equals a known function of itself, it may be possible to use a fixed-point iteration to find it, for example, Eq. (5) can be rearranged with $μ_{a}$ on one side and $H / [Φ (μ_{a})]$ on the other. Fixed-point iterations typically converge quickly to the correct solution if they converge at all. This is applied to QPAI in Sec. 5.1.3 in the case where the scattering is known.

Model-based minimization

A very general approach that does not require the model to be known analytically is to find the unknowns by solving the forward problem iteratively, updating the unknowns at each iteration, until the output of the solver matches the measured data in some sense. This type of scheme occurs in many guises which differ both in the way that the solver output and the measured data are compared (the error functional) and in the way in which the unknown parameters are updated (the minimization algorithm). A popular and powerful subset of these methods are least-squares approaches that have been applied to QPAI and are described below.

2.3.2.

Least-squares minimization

Least-squares minimization is a common and well-developed framework for solving inverse problems. If the data is a measured (observed) absorbed energy distribution, $H^{obs}$ , then the squared error between it and the output of a forward model, $H^{model} (C_{k})$ , is ${‖ H^{model} (C_{k}) - H^{obs} ‖}^{2} = \sum_{i} {[H_{i}^{model} (C_{k}) - H_{i}^{obs}]}^{2}$ , where the subscript $i$ indicates the value at the $i$ ’th voxel. The aim in least-squares minimization is to minimize the error functional, $ε$ , by adjusting the unknowns $C_{k}$ :

Eq. (18)

\underset{C_{k}}{argmin} ε = \frac{1}{2} {‖ H^{model} (C_{k}) - H^{obs} ‖}^{2} + P (C_{k}) .

The second term,

P

, is a penalty term that can be used to include regularization or other constraints. Even when

ε

has a well-defined minimum in the noise-free case, where the difference between

H^{model}

and

H^{obs}

is small, the noise in the data may mean there may be several solutions,

H^{model}

, that fit the measured data,

H^{obs}

, equally well. Regularization is a general term to describe methods for ameliorating this problem, and common approaches are to encourage the unknowns—the absorption or scattering coefficients, say—to be smooth (Tikhonov) or piecewise constant (total variation). The latter may be reasonable in some cases from physiological considerations; for example, when imaging a plane through a blood vessel, the absorption coefficient may take one value inside the vessel and one outside. Other prior knowledge of some aspect of the unknowns, or perhaps of an intermediate quantity, such as shape or an empirically determined maximum or minimum, may also be included in this way. The basic idea is to reduce the solution space by excluding functions with certain properties, such as those that are insufficiently smooth, from being considered as solutions.

Solving Eq. (18) means finding the minimum of the scalar functional $ε$ . Gradient-based methods use the functional gradient (the vector of first derivatives), $g_{i} = \partial ε / \partial C_{k, i}$ , such as the conjugate-gradient method, to step down the function until a minimum is reached. Hessian-based methods such as Newton’s method try to reduce the number of iterations by also using the Hessian matrix of second derivatives, $H_{i j} = \partial^{2} ε / (\partial C_{k, i} \partial C_{k, j})$ (related to the curvature of the function). Some methods approximate the Hessian matrix to ameliorate the burden of calculating it explicitly. For example, the Gauss–Newton method uses the Jacobian matrix $J_{i j} = \partial H_{i}^{model} / \partial C_{k, j}$ to estimate the Hessian matrix as $H \approx J^{T} J$ , which can be implemented with Krylov methods that require only matrix-vector products to be computed.³⁹ Quasi-Newton methods such as L-BFGS estimate the Hessian matrix at each iteration by using stored values of the gradient. Gradient-free methods, such as the Nelder-Mead simplex method, use neither gradient nor Hessian information, so they tend to be slow to converge. These approaches are applied to QPAI in Secs. 5.1.5 and 5.2.3 to 5.2.5. It should be noted that the literature on all aspects of least-squares minimization (or optimization) is substantial.⁴⁰^,⁴¹

2.3.3.

Uniqueness and ill-posedness

Estimating the chromophore concentrations $C_{k}$ by minimizing the functional in Eq. (18) is not as straightforward as it might at first seem. The principal reason is that when the scattering is unknown, and therefore needs to be estimated at the same time as the concentrations, the functional $ε$ may not have one unique minimum. Consider the absorbed energy density, $H$ , at a single wavelength. An increase in the absorption coefficient at one point will increase the number of photons that are absorbed there, which will have the effect of reducing the fluence nearby. The fluence can also be altered by changing the optical scattering, so a situation may arise in which a change in the scattering occurs such that the resulting change in the fluence counteracts exactly the effect of the absorption increase on $H$ . In other words, when the scattering and absorption coefficients are both allowed to vary spatially, $H$ may not depend uniquely on the optical parameters: two different absorption and scattering distributions could lead to the same $H$ . As far as the optical inversion is concerned, this is a severe problem because it means there is no unique solution to the question “Which optical coefficients would result in the measured $H$ ?” To solve the inversion it is essential that this nonuniqueness be removed. This can be done in a number of different ways: for example, by fixing the scattering (Sec. 5.1), although this can clearly lead to bias if the scattering is not known accurately; by using knowledge of the scattering (and chromophores’) wavelength dependence;⁴²^,⁴³ or by using multiple measurements made with different surface illumination patterns.⁴⁴^,⁴⁵

A less severe form of ill-posedness is caused by the diffusive nature of the optical propagation, which tends to smear out sharp features in the fluence. This means that high spatial frequencies in the distributions of the optical properties have limited influence on the fluence distribution. In PA, because it is the absorbed energy density $H$ rather than the fluence directly that gives rise to the PA signals, the high frequencies (e.g., sharp edges) in the absorption coefficient distribution do have a significant influence on the measured data. However, the scattering coefficient only affects $H$ through its impact on the fluence so has a second-order effect on $H$ , which decreases as the spatial frequency increases. This is equivalent to saying that the operators $T$ and $T_{λ}$ act as low pass filters to reduce the amplitudes of the high frequency components of the scattering distribution, and consequently the effect of the inverse operations, $T^{- 1}$ , $T_{λ}^{- 1}$ , will be to grow the high frequency components. This will have the effect of amplifying the noise in the measured data which, unchecked, may come to dominate the inversion. Alterations to the inverse operator that are designed to reduce this unwanted effect are termed regularization (as mentioned in the context of penalty functionals in Sec. 2.3.2 above). This type of phenomenon is common to many inverse problems, both linear and nonlinear, and a large number of regularization techniques are described in the literature.⁴⁶

2.3.4.

Large-scale inversions and domain parameterization

A practical difficulty can arise when attempting the QPAI inversion on 3-D PA images. When considering the general problem of recovering the spatially varying chromophore concentrations, the value of each chromophore in each image voxel could be treated as a separate unknown. As a 3-D PAT image may consist of $10^{7}$ or more voxels, there may be of the order of $10^{8}$ unknowns. This constitutes a large-scale inverse problem and poses some practical difficulties. For example, a Hessian-based inversion would need to compute, store, and invert a Hessian matrix with perhaps $10^{16}$ elements on each iteration, which is currently impractical. Clearly, techniques to reduce the size of the problem and inversion procedures that are computationally light⁴⁷^,⁴⁸ are required. One way to reduce the number of unknowns is to divide the domain into a few regions on which the optical coefficients are assumed constant. For example, if the regions are denoted $A_{n}$ , then the absorption coefficient can be written as the sum

Eq. (19)

μ_{a} (x) = \sum_{n} μ_{a, n} S_{n} (x), S_{n} (x) = {\begin{matrix} 1, & x \in A_{n} \\ 0, & otherwise \end{matrix} .

The concentrations and scattering coefficient can similarly be separated into piecewise constant regions. Such a parameterization has two advantages: it reduces the number of unknowns to a manageable number, and it can act to regularize the inversion. The regions $A_{n}$ may either be chosen based on some prior information about the underlying tissue structure, or a multigrid approach may be used, in which the regions are initially large and are iteratively reduced.

3. One-Dimensional Quantitative Photoacoustic Imaging

One of the areas of research that fed into and helped generate the interest in PA imaging was the work done in the late 20th century on PA depth-profiling. Much of this literature is concerned with estimating the absorption coefficient of the material under study, sometimes as a function of depth, sometimes in a scattering medium, and sometimes at multiple wavelengths, so several of the key issues that arise in quantitative PA imaging in higher dimensions appear in the 1-D context. As well as the possibility of gleaning insight into the higher-dimensional versions of these problems, another reason for reviewing attempts on the 1-D problem is that under some circumstances, real situations can be modeled as 1-D or quasi-1-D, in particular in certain PA microscopy applications.

3.1.

Homogeneous Nonscattering Medium: Beer’s Law

The simplest situation is a homogeneous, nonscattering, optically absorbing medium that is illuminated (with a short pulse) by an infinitely wide beam of light. The absorbed energy density under such conditions can, following the Beer’s law, Eq. (11), be written simply as

Eq. (20)

H (z) = {\begin{array}{r} μ_{a} Φ_{0} \exp (- μ_{a} z), & z \geq 0 \\ 0, & z < 0 \end{array},

where

Φ_{0}

is the fluence at the illuminated surface of the tissue, which is assumed, without loss of generalization, to be at

z = 0

. This initial pressure will give rise to two PA waves, traveling in the

+ z

and

- z

directions, each with half the amplitude but retaining the exponential shape. For a pressure detector at

z = 0

(backward mode detection), then the signal reaching the detector will be simply

Eq. (21)

p (t) = (\frac{\hat{Γ} Φ_{0} μ_{a}}{2}) \exp (- μ_{a} c_{0} t), t \geq 0,

where

c_{0}

is the sound speed. (We are ignoring acoustic reflections at the surface here, but they can be straightforwardly included.) The absorption coefficient,

μ_{a}

, can be estimated either from the maximum amplitude of the signal, if the PA efficiency,

\hat{Γ}

, and the incident fluence,

Φ_{0}

, are known, or by fitting a curve to the exponentially decaying slope.⁴⁹^–⁵¹

More recently, a frequency domain method for quantification of $μ_{a}$ has been proposed by Guo et al.⁵² It too assumes a homogeneous, nonscattering medium, and the starting point is that the measured data will be $p (t)$ as given in Eq. (21) convolved with a frequency-dependent transfer function that will depend on both the acoustic absorption and the measurement system frequency response. As the magnitude of the Fourier transform of $p (t)$ is $| P (ω) | = (\hat{Γ} Φ_{0} μ_{a} / 2) {[{(μ_{a} c_{0})}^{2} + ω^{2}]}^{- \frac{1}{2}}$ , the ratio of two measurements made at different optical wavelengths, $λ_{1}$ and $λ_{2}$ , may be written as

Eq. (22)

\frac{| P (λ_{1}; ω) |}{| P (λ_{2}; ω) |} = \frac{Φ_{0} (λ_{1})}{Φ_{0} (λ_{2})} {\frac{c_{0}^{2} + {[ω / μ_{a} (λ_{2})]}^{2}}{c_{0}^{2} + {[ω / μ_{a} (λ_{1})]}^{2}}}^{1 / 2},

where the system transfer function, the acoustic absorption term, and the PA efficiency have canceled out. The three unknown numbers

μ_{a} (λ_{1})

,

μ_{a} (λ_{2})

, and the surface fluence ratio

Φ_{0} (λ_{1}) / Φ_{0} (λ_{2})

may be obtained by curve-fitting Eq. (22) to measured acoustic spectra.

Guo et al. used this technique for quantification by OR-PAM, which raises the question as to how a model that explicitly varies only with depth, $z$ , can be justified for use with OR-PAM, where there is variation in lateral dimensions too. There are three conditions that must be satisfied for the 1-D planar assumption to hold and, therefore, for Eq. (22) to be true: (1) the medium properties, (2) the absorbed energy distribution, and (3) the acoustic propagation must all be planar (1-D) on the scale of interest. The first condition is satisfied as the light beam in OR-PAM is focused to a spot ( $\sim 5$ μm) much smaller than a typical vessel diameter ( $\sim 30$ μm), so the vessel can be considered to be a homogeneous half-space. The second may be satisfied in the region close to the surface where ballistic photons dominate the fluence. The third, however, requires the acoustic waves to be planar (or at least planar in a sufficient region that the detector cannot tell that they are not planar everywhere), which does not seem to be the case. The restriction of the illumination to a zone with a lateral dimension of $\sim 5$ μm and a depth of perhaps one order of magnitude greater will not generate purely plane waves propagating in the $z$ direction, even within the acoustic focus. There will be components propagating at angles to the $z$ -axis, so the detected acoustic wave will therefore not vary exponentially according to $\exp (- μ_{a} c_{0} t)$ , and Eq. (22) will not hold.

Although this method appears to be justifiable only when the illumination extends over a much larger region, Guo et al. used it to analyze phantom OR-PAM measurements and showed that for a highly absorbing homogenous ink phantom (30 to $225 {mm}^{- 1}$ ) was recovered to within about $\pm 4 {mm}^{- 1}$ . The method was also used to estimate absolute values for the absorption coefficients (and subsequently the oxyhemoglobin, deoxyhemoglobin, and total hemoglobin concentrations and ${sO}_{2}$ ) in a superficial vein and artery pair selected from $1 {mm}^{2}$ OR-PAM images of a nude mouse ear obtained at 561 and 570 nm.

3.2.

Depth-Dependent Nonscattering Medium

For multilayered media in which each layer has a different absorption coefficient, different exponentials could be fitted to the parts of the curve corresponding to the different layers. However, this becomes difficult when the layers are thin and there is only a short region of curve to fit to. On the basis that a layer of thickness $Δ z$ will absorb energy per unit area $Φ [1 - \exp (- μ_{a} Δ z)]$ , where $Φ$ is the fluence of the light entering the layer, an expression for the absorbed energy density within a stack of thin layers can be found.⁵³^,⁵⁴ When the layers become so thin that the absorption coefficient becomes a continuous function of depth, the absorbed energy density may be written as

Eq. (23)

H (z) = μ_{a} (z) Φ_{0} \exp [- \int_{0}^{z} μ_{a} (ζ) d ζ] .

The fluence here appears as a form of Beer’s law, Eq. (11), generalized to depth-dependent media. Karabutov et al.⁵⁵^,⁵⁶ inverted this expression to find the following formula for the absorption coefficient given the measured pressure

p (t) = \hat{Γ} H (t = {z / c}_{0})

:

Eq. (24)

μ_{a} (z) = \frac{p (z / c_{0})}{c_{0} \int_{z / c_{0}}^{\infty} p (t) d t}, z \geq 0 .

This has been experimentally verified by recovering absorption profiles for depth-dependent solutions of magnetite particles in oil⁵⁵^,⁵⁶ and, in a slightly modified form to account for the tail of the pressure signal $p (t)$ when it is only known up to a finite time, for dyed gelatin phantoms.⁵⁷^,⁵⁸

3.3.

Homogeneous Scattering Medium

In biological tissue the light is scattered significantly for depths greater than a few hundred microns, and this must be taken into account. In 1-D the introduction of a significant level of optical scattering has two principal effects. First, the decay rate of the fluence (beyond a few mean free paths) is no longer governed purely by the absorption coefficient, so fitting an exponential to the decaying part of the curve will not recover $μ_{a}$ . For a collimated beam, the fluence close to the surface will decay approximately as $\exp (- μ_{t} z)$ and at deeper depths where the light is diffuse as $\exp (- μ_{eff} z)$ (see Sec. 2.2.2). Second, the maximum value of the fluence may not be at the surface but some distance below it (due to backscattering), so the maximum amplitude of the signal cannot be used directly to estimate $μ_{a}$ either. Oraevsky et al.⁵⁹^,⁶⁰ used the approximate model for the fluence⁶¹ given by Eq. (17), according to which the absorbed energy profile is $H (z) \approx μ_{a} k Φ_{0} \exp (- μ_{eff} z)$ in the diffuse regime. The factor $k$ accounts for the backscattered light and the resulting increase in the absorbed energy density below the surface and is given by $k = 1 + 7.1 R_{d \infty}$ , where $R_{d \infty}$ is the diffuse reflectance from the surface.⁶¹ By measuring this diffuse reflectance and fitting Eq. (17) to the exponentially decaying part of the measurement, $μ_{a}$ could be inferred from the amplitude of the fitted curve extrapolated to $z = 0$ , that is, from $μ_{a} k Φ_{0}$ . Subsequently $μ_{eff}$ was measured from the slope of the curve, and $μ_{s}^{'}$ calculated from it, by use of the known $μ_{a}$ . Rather than measuring the diffuse reflectance, Fainchtein et al.⁶²^,⁶³ modeled it as approximately $R_{d \infty} \approx \exp (- 7 μ_{a} / μ_{eff})$ . Note that in both cases an additional optical measurement, or an additional assumption, is required to allow both $μ_{a}$ and $μ_{s}^{'}$ to be estimated from the PA signal.

3.4.

Measurements in Blood

In this section we appear to make a detour from the main theme to discuss blood. The main reason is that the hemoglobin in blood is the most important source of contrast for PA imaging (although the rationale for discussing it here is just that several researchers have made 1-D PA measurements of the properties of blood). There are several reasons why it is so prominent in PA imaging studies and so important: (1) it is naturally occurring so there is no need for an exogenous contrast agent, (2) hemoglobin is the dominant endogenous chromophore in the wavelength range that corresponds to the “near infrared window” where the deepest penetrations are possible, and (3) multiwavelength measurements of the absorption of blood have the potential to provide functional information about the tissue through the oxygen saturation, ${sO}_{2}$ , which has direct physiological relevance. One goal of QPAI is to obtain 3-D images in which the voxel values are accurate estimates of ${sO}_{2}$ . The precision and accuracy with which PA methods can be used to determine the properties of blood, such as the level of oxygenation and the total hemoglobin concentration, are therefore of great interest.

3.4.1.

Blood oxygen saturation

Blood oxygen saturation, ${sO}_{2}$ , is defined as the ratio

Eq. (25)

{sO}_{2} = \frac{C_{{HbO}_{2}}}{C_{{HbO}_{2}} + C_{HHb}},

where

C_{{HbO}_{2}}

and

C_{HHb}

are the concentrations of oxyhemoglobin and deoxyhemoglobin respectively, which are related to the optical absorption coefficient of blood via their molar absorption coefficient spectra

a (λ)

:

Eq. (26)

(\begin{matrix} μ_{a, blood} (λ_{1}) \\ ⋮ \\ μ_{a, blood} (λ_{N}) \end{matrix}) = (\begin{matrix} α_{{HbO}_{2}} (λ_{1}) & α_{HHb} (λ_{1}) \\ ⋮ & ⋮ \\ α_{{HbO}_{2}} (λ_{N}) & α_{HHb} (λ_{N}) \end{matrix}) (\begin{matrix} C_{{HbO}_{2}} \\ C_{HHb} \end{matrix}) .

This matrix of molar absorption coefficients is the linear spectroscopic mapping

L_{λ}

in Fig. 2. As

{sO}_{2}

is a ratio, the concentrations

C_{{HbO}_{2}}

and

C_{HHb}

do not need to be known absolutely but only to within a multiplicative constant. As long as the constant is the same for both, then it will cancel out. This implies that a relative measurement of the absorption coefficient,

K μ_{a, blood} (λ)

, will be sufficient to determine

{sO}_{2}

as long as the multiplicative factor

K

does not depend on wavelength.¹² To use the amplitudes of PA signals as relative measurements of

μ_{a}

, it is necessary to find a scenario in which the constant

K

, which for PA will be

K = \hat{Γ} Φ

, is independent of wavelength. This is not true in general, in fact it is rarely true, and assuming it to be so⁶⁴ is likely to result in significant errors in

{sO}_{2}

estimates. If PA measurements of

{sO}_{2}

are ever to become widely used and trusted in clinical practice, then this issue needs to be addressed properly. Of course, this is just one instance of the more general problem that, to obtain absolute estimates of concentrations such as

C_{{HbO}_{2}}

and

C_{HHb}

, the wavelength dependence of the fluence must be accounted for.

3.4.2.

Cuvette measurements

Early experiments using PA to measure the properties of blood were mostly conducted with cuvettes, where the blood is presented as a 1-D homogeneous target and the experimental parameters are well-controlled. Using a low frequency ultrasound detector ( $< 1 MHz$ ), Fainchtein et al.⁶²^,⁶³ made measurements of the wavelength dependence of the amplitude of the peak of the PA wave which showed qualitative agreement with the absorption spectrum of blood (710 to 870 nm), both for blood in a cuvette and in a canine arterial-venous shunt. The measured spectrum changed approximately as expected as the level of oxygen in the blood was varied, but no quantitative estimates were made. Savateeva et al.⁶⁵ used a broadband detector ( $\sim 40 MHz$ ) to measure the exponential slope of the PA signal in order to estimate the attenuation coefficient and showed that it varied linearly with the level of oxygen saturation at 532, 757, and 1064 nm wavelengths. Esenaliev et al.⁶⁶ showed a similar result at 1064 nm with a lower bandwidth transducer. However, none of these authors provided estimates of absolute ${sO}_{2}$ .

Building on these preliminary results, Laufer et al.⁶⁷ estimated blood ${sO}_{2}$ in a cuvette by making measurements of the slopes and peak-to-peak amplitudes of PA signals with a broadband detector, $\sim 15 MHz$ , from 740 to 1040 nm in 10 nm increments. They found that the light passing through the cuvette was sufficiently affected by scattering (due to the blood cells) that it could not be modeled as $\exp (- μ_{a} z)$ and instead was modeled by a 1-D version of the $δ$ -Eddington approximation. This model was fitted by a Nelder-Mead algorithm to both the measured PA amplitudes and, separately, to measurements of $μ_{eff}$ obtained from the slopes of the exponentially decaying PA signals. The ${sO}_{2}$ was obtained with an accuracy of $\pm 2.5 %$ when using measurements of $μ_{eff}$ and $\pm 4 %$ from amplitude measurements. The smallest change in ${sO}_{2}$ that could be accurately measured was $\pm 1 %$ . These promising results suggest that PA measurements can be used to estimate ${sO}_{2}$ sufficiently accurately to be clinically useful, but the extension of this technique from 1-D to 3-D is a nontrivial matter.

3.5.

Applicability of 1-D Methods in Practice

For all methods discussed in this section, which are based on 1-D assumptions, the accuracy of the recovered parameters will depend on how well the experiment approximates to 1-D. What must be the case for this to be a good approximation? There are three things to consider, and they are not independent: the acoustic propagation, the fluence distribution, and the medium properties. For the acoustic propagation to be considered as 1-D, the PA wave that is generated must be a plane wave perpendicular to, say, the $z$ -axis [i.e., it must be invariant in the ( $x$ , $y$ ) plane], or rather it must be planar in a sufficiently large region that it appears planar to the detector. For an omnidirectional detector, no signals will be detected from outside a radius $r = c_{0} T$ from the detector, where $T$ is the time over which the signal is measured, so the wave need only be planar within this region. This requirement can be relaxed if, for example, the detector is highly directional, as it will only be sensitive to parts of the wave arriving normally.

For the initial acoustic pressure distribution to generate only plane waves, it is necessary that the fluence distribution and medium properties are invariant in the ( $x$ , $y$ ) plane too. Two separate cases will be considered for the fluence: nonscattering and scattering media. In a nonscattering medium, the condition that the fluence varies only with depth requires a collimated source of light that does not diverge significantly and has a radial profile that is flat over the region with radius $c_{0} T$ . When these conditions hold in a region of radius $R$ , then Eq. (21) is true for $t < R / c_{0}$ . In scattering media, in the superficial region close to the illumination surface where $z ≪ 1 / μ_{t}$ , where ballistic (unscattered) photons dominate over scattered photons, the fluence will decay one-dimensionally as $Φ (z) = Φ_{0} \exp (- μ_{t} z)$ and so the acoustic signal can be written as $p (t) = (\hat{Γ} Φ_{0} μ_{a} / 2) \exp (- μ_{t} c_{0} t)$ for $t ≪ 1 / (c_{0} μ_{t})$ or $t < R / c_{0}$ , whichever is smaller. When considering deeper depths where the scattered photons dominate the field, there is a much more stringent requirement. In order to assume that the light fluence decays as $\exp (- μ_{eff} z)$ in the diffusive regime, $z ≫ 1 / μ_{t}$ , it has been shown that it is necessary to ensure that the illumination is constant over a much wider region than the sensitive zone of the detector, perhaps as large as $r = 10 c_{0} T$ .²⁵ With smaller illumination regions, the light fluence will decay faster than $\exp (- μ_{eff} z)$ because of geometrical spreading. Clearly, for the fluence to depend only on depth it is not enough to have a sufficiently broad source, it is also necessary that the optical properties depend only on depth in the region of interest too.

As a final comment on 1-D methods, it is worth recalling that diffraction cannot occur in 1-D, hence divergences from 1-D are sometimes described as “problems of diffraction.”⁵⁷^,⁵⁸^,⁶⁸^,⁶⁹ Indeed, the lack of diffraction and the positivity of the initial pressure distribution give a simple way to check experimentally that the 1-D assumption holds: a truly planar 1-D PA signal detected by a sufficiently broadband detector will not contain negative components.⁷⁰

4. Quantitative Estimates from Single-Point Measurements

One-dimensional characterizations are often inappropriate for measurements made in biological tissue, where the optical properties are heterogeneous and the geometrical spreading of the acoustic waves is important. This section describes methods that have been devised to overcome this limitation but are still based on a single time series measurement from one detector. The detector may be one in an array of detectors or, in some cases, be scanned in order to generate an image, as in PAM. The difference between this section and Sec. 5 is that here the quantitative information is extracted separately for each point, that is, for each time series, whereas in Sec. 5 all the measurements are used together in the QPAI inversion, for example, by generating a PAT image and using that as the primary data in the inversion.

4.1.

Single Cylindrical Absorber

4.1.1.

Linearity with bandlimited detection

Assume that a scattering medium contains a single absorber of known shape and position relative to a detector, so for a known illumination and detector response the relationship relating $μ_{a}$ of the absorber to the amplitude of the measured PA wave can be derived. This is one instance of the operator $W$ in Fig. 2 and, in certain cases, it may be found that this function is approximately linear, which allows $W^{- 1}$ to be obtained trivially. This approach was taken by Sivaramakrishnan et al.¹²^,⁷¹ for the case of a cylindrical absorber (simulating a blood vessel). They considered a uniformly absorbing cylinder of radius $a$ illuminated from all sides equally (i.e., within the diffusive regime), and a bandlimited detector perpendicular to the vessel axis with a central frequency of $f_{c} = c_{0} / λ_{c}$ , where $λ_{c}$ is the associated center wavelength. By calculating the signal from the detector a fixed distance from the vessel, they found two scenarios for which the peak PA amplitude varied approximately linearly with $μ_{a}$ : (1) for small vessels $a ≪ 1 / μ_{a}$ , and (2) for large vessels $a > λ_{c}$ when a transducer with $λ_{c} < 1 / μ_{a}$ is used. For a suitably isolated blood vessel, this linear model provides a way of measuring the absorption coefficient spectrum to within a constant factor, which can be used with Eqs. (25) and (26) to calculate blood ${sO}_{2}$ . By making a measurement at an isosbestic point (a wavelength at which $α_{{HbO}_{2}} = α_{HHb}$ ), a relative estimate of the total hemoglobin concentration, HbT, (i.e., to within an unknown constant factor) can also be found. Measurements on an experimental phantom consisting of an ink-filled tube of diameter 0.25 mm confirmed that the signal measured by the transducer (with a central frequency of 25 MHz and 90% bandwidth) varied linearly with the absorption coefficient up to $μ_{a} \approx 180 {cm}^{- 1}$ ( $μ_{a} λ_{c} \approx 1$ , $λ_{c} < a$ ). This was used to successfully recover the wavelength dependence of the absorption coefficient of blood between 575 and 600 nm. The same experiment with a 10 MHz transducer ( $μ_{a} λ_{c} \approx 3$ , $λ_{c} < a$ ) failed to recover the spectrum correctly.¹²

The fact that this method reduces the nonlinear inversion $W^{- 1}$ to a linear problem is appealing. However, its main limitation lies in the fact that it is based on a model of a single absorber in a purely scattering (nonabsorbing) background medium, which is not representative of real tissue that will contain additional vessels, capillaries, extravascular blood, lipids, water, and potentially other absorbing molecules such as melanin. The spectrum of the light reaching the vessel of interest will therefore be colored by its passage through the background medium because of the additional chromophores present, and the measured spectrum cannot confidently be related to the absorption in the vessel alone.

4.1.2.

Invasive correction for wavelength dependence

To translate the linear method above from phantom measurements into real tissue, a measurement or estimate of the wavelength dependence of the fluence local to the vessel is required. Wang et al.⁷² used transmission measurements of light through ex vivo skin and skull samples to arrive at a first-order correction factor for the spectrum. With this correction, ${sO}_{2}$ was estimated in rat brain vasculature. This relies on the optical properties of the ex vivo samples being representative of in vivo conditions, which is doubtful. In an attempt to overcome this limitation, Maslov et al.⁷¹^,⁷³^–⁷⁵ inserted a plain black absorbing film with a spectrally flat absorption coefficient at the depth of interest into the tissue at the level of the vessel and measured the spectrum of the PA signal generated by it. Under conditions in which there is no backscattering or the backscattered light is negligible, this method would measure the correct wavelength dependence of the fluence where the light remains collimated, but where there is significant backscattering (as there is in tissue), the part of the fluence caused by photons traveling back upwards to the vessel from below will not be accounted for. Clearly, this invasive preliminary measurement cannot be performed on targets that one would like to study noninvasively and longitudinally. The invasive procedure might be performed on one animal and the noninvasive measurement made on another, but the possibility of variation between the two introduces uncertainty in the wavelength correction. Despite these uncertainties, this approach has been applied to the estimation of ${sO}_{2}$ in the skin of small animals using PA microscopy.⁷³^,⁷⁴^,⁷⁶

4.1.3.

Correction for wavelength dependence using a known contrast agent

In a similar spirit but less invasively, Rajian et al.⁷⁷ propose the use of an exogenous contrast agent instead of a black absorbing layer to help estimate the local fluence. The principle is the following: a first measurement is made in the usual way in which the PA signal is proportional to the product of the fluence and absorption from just the endogenous chromophores,

Eq. (27)

p_{1} (λ) = B Φ (λ) \sum_{k} C_{k} α_{k} (λ),

where

B

is an unknown scaling factor (which will depend on the PA efficiency, the detector sensitivity, etc.) and

Φ (λ)

is the unknown wavelength-dependent fluence. A small quantity of contrast agent is then introduced, so the second PA measurement is given by

Eq. (28)

p_{2} (λ) = B^{'} Φ^{'} (λ) [\sum_{k} C_{k}^{'} α_{k} (λ) + C_{CA} α_{CA} (λ)],

where

C_{CA}

is the contrast agent concentration and

α_{CA} (λ)

is its known absorption spectrum. The primes indicate that the quantities may have changed from the introduction of the contrast agent. Under conditions where the fluence, the concentrations of the endogenous absorbers, and the constant

B

do not change much when the contrast agent is introduced, so that

Φ (λ) \approx Φ^{'} (λ)

,

C_{k} \approx C_{k}^{'}

, and

B \approx B^{'}

, the local fluence can be estimated as

Eq. (29)

Φ (λ) \approx \frac{p_{1} (λ) - p_{2} (λ)}{{BC}_{CA} α_{CA} (λ)} .

Substituting this into Eq. (27) leads to the linear relationship

Eq. (30)

\frac{C_{CA} α_{CA} (λ) p_{1} (λ)}{p_{1} (λ) - p_{2} (λ)} \approx \sum_{k} C_{k} α_{k} (λ),

which, so long as

C_{CA}

is known, can be inverted straightforwardly to give the concentrations of the endogenous chromophores. This approach has the advantage that the (often unknown) scaling factor

B

cancels out. However, it requires the assumptions that the fluence is the same before and after the injection of the contrast agent and that the concentration of the contrast agent is known accurately. Unfortunately, it is unlikely that the fluence is insensitive to absorption changes large enough to affect the measured PA signal.

4.2.

Model-Based Minimization

Laufer et al.⁷⁸ describe a method that can both account for absorption in the tissue surrounding the vessel of interest and can be more easily generalized to different targets and measurement configurations. These are part of a series of papers⁶⁷^,⁷⁸^,⁷⁹ describing the evolution of a minimization-based inversion scheme referred to briefly in Sec. 3.4.2. They considered the PA time series recorded at a single transducer from three tubes each lying perpendicular to the transducer, each at a different fixed distance from the transducer along the line of sight, and all three simultaneously illuminated from the opposite side from which the signals were detected. To mimic real tissue, the background medium was both scattering and absorbing. The separation of the tubes was such that the PA signals from each tube could be identified separately. This allowed spectra to be recorded from the peak-to-peak amplitudes of the PA signals: one spectra per tube and three from the extraluminal space. A two-stage forward model was used in the inversion (i.e., separate light and acoustic models, $T_{λ}$ and $A$ ). The light transport was modeled with a two-dimensional (2-D) finite element $δ$ -Eddington diffusion model, and a simplified acoustic propagation model was used to convert the absorbed energy density to a pressure time series. Six modeled spectra were calculated from this modeled time series as above, by taking the peak-to-peak values. For each of the three tubes, $i = 1, 2, 3$ , the intraluminal absorption coefficient was written as the sum of contributions from three absorbers: oxyhemoglobin, deoxyhemoglobin, and water; the absorption outside the tubes also included lipids, a contrast agent (NIR dye), as well as blood and water:

Eq. (31)

tubes μ_{a, i} (λ) = C_{HHb, i} α_{HHb} (λ) + C_{{HbO}_{2}, i} α_{{HbO}_{2}} (λ) + μ_{a, H_{2} O} (λ), background μ_{a, bg} (λ) = C_{HHb, bg} α_{HHb} (λ) + C_{{HbO}_{2, bg}} α_{{HbO}_{2} (λ)} + μ_{a, H_{2} O} (λ) + μ_{a, lipids} (λ) + C_{dye} α_{dye} (λ) .

The scattering coefficient was written as $μ_{s} (λ) = k_{s} α_{scat} (λ)$ , where the wavelength dependence of the scattering $α_{scat} (λ)$ was considered known in advance, although the amplitude $k_{s}$ was not. There were therefore 11 unknown scalars to be determined: $C_{{HbO}_{2}}$ and $C_{HHb}$ for each tube and the background; $C_{dye}$ in the background; scattering $k_{s}$ ; and an overall amplitude scaling factor, $K$ , related to the incident fluence and the PA efficiency. These were found using a least-squares minimization (Sec. 2.3.2), specifically the Nelder-Mead method. This then allowed the oxygen saturation, ${sO}_{2}$ , and the total hemoglobin concentration to be estimated. The results compared favorably with those of a laboratory CO-oximeter, with a resolution of $\pm 4 %$ and accuracy in the range $- 6$ to +7%. It should be emphasized that the inclusion of chromophores of unknown concentration in the background medium as part of the inversion precludes having to make wavelength corrections of the sort described in Secs. 4.1.2 and 4.1.3.

In a subsequent paper,⁷⁹ Laufer et al. improved and generalized this technique by improving the correspondence between the model and the actual physical situation. This was done in a number of ways: by using a full-wave 3-D acoustic propagation model rather than a simplification, by including a correction factor to account for the difference between the 2-D light model used and the 3-D nature of the fluence distribution in practice (a further improvement would be to use a full 3-D light model), and most significantly by taking the spectra from 2-D PA images obtained at multiple wavelengths, rather than from the pressure time series amplitudes. Because this approach (1) can estimate the unknown (although constant) scattering rather than assuming it is known or negligible, (2) can obtain absolute concentrations (e.g., of oxyhemoglobin and deoxyhemoglobin), (3) has been demonstrated with measured multiwavelength data from known phantoms, and (4) is applicable to full 3-D PAT images (i.e., it is not restricted to superficial or single vessels), this technique has come closest to date to a practical and generally applicable method for quantitative spectroscopic PAT.

The principal remaining restriction in this model-based approach is the prior knowledge of the target geometry that is required. In each of these inversions, to reduce the number of unknowns to a manageable number, the domain was divided into regions assumed to have constant optical properties, for example, into a few regions corresponding to dye-filled tubes, and the remainder as a “background” region (see Sec. 2.3.4). The amount known about the parameterization a priori was reduced as the sophistication of the inversion increased. In Ref. 78 it was known that the absorbers were tubes lying on axis, but their depths and diameters were estimated from the measured PA time series. In Ref. 79 all that was known beforehand was that the absorbers were tubes, and both their positions and sizes were found from PA images. In principle, these positions and sizes could even be included as unknowns in the inversion. The assumption that the domain has piecewise constant optical properties is a significant assumption that simplifies the problem in two related ways: first, it reduces the number of unknowns from potentially many tens of thousands to just a few so the posedness of the inversion improves and it is much simpler to compute the inversion; second, it reduces the effect of errors in the PAT image on the quantitative estimates. To use images in the inversion requires measurements made over an array of detectors rather than single-point detection, but this latter inversion was nevertheless included in this section because it follows logically from the preceding work. In the sense that it uses images, however, it provides a link to the next section where the inversions start with an estimate of the absorbed energy density distribution.

5. Full-Field Quantitative Photoacoustic Imaging

This section is concerned with techniques to convert tomographic images of the absorbed energy density, $H (x)$ , here considered as the measured data, into images of optical coefficients or chromophore concentrations. In other words, the acoustic inversion, $A^{- 1}$ , is assumed solved, and the focus in now wholly on the optical inversions $T^{- 1}$ or $T_{λ}^{- 1}$ without the simplifying approximations that have been made in Secs. 3 and 4. Some simplifying assumptions are still made however. In particular, most of the methods in Sec. 5.1 assume the scattering is known and just recover the absorption coefficient, and the methods in Sec. 5.2 assume the diffusion approximation holds, which will not usually be true close to the surface. It is also worth emphasizing that most of the methods in Sec. 5 recover the absorption coefficient at a single wavelength, so to obtain chromophore concentrations will require the acquisition of images at different wavelengths and a subsequent multispectral inversion, $L_{λ}^{- 1}$ . Section 5.2.4 solves the combined spectral and optical inversion, $T_{λ}^{- 1}$ .

No particular distinction is made between 2-D and 3-D in this section as none of the methods described depends critically on the difference, although most have been demonstrated with 2-D rather than 3-D data sets as the computations are more manageable.

5.1.

Inversion for Absorption Coefficient Only

When the PA efficiency is known and the image reconstruction is free of errors or artifacts, the PA image is an image of the absorbed energy density which can be written as [Eq. (2)]

Eq. (32)

H (x) = μ_{a} (x) Φ [x, μ_{a} (x)],

where the dependence of the fluence on the absorption is shown explicitly. The principal challenge in this inversion is the nonlinear dependence of

H

on

μ_{a}

. A common way to tackle such problems is to notice that the nonlinear behavior is close to linear if the changes are small. Consider the case in which the absorption is a perturbation over a known and homogeneous background,

μ_{a} (x) = μ_{a, 0} + δ μ_{a} (x)

. The resulting fluence can be written as a sum of an unperturbed part,

Φ_{0} (x)

, corresponding to

μ_{a, 0}

, and a second part

δ Φ (x)

, corresponding to the absorption perturbation

δ μ_{a} (x)

. The change in the absorbed energy distribution will (neglecting the second-order term) consist of two terms:

Eq. (33)

δ H (x) \approx δ μ_{a} (x) Φ_{0} (x) + μ_{a, 0} δ Φ (x),

where the first term is due to the local absorption coefficient perturbation, and the second to changes in the local fluence caused by the absorption perturbation.

5.1.1.

Unchanged fluence

In Eq. (33), the change in the absorbed energy, $δ H$ , is still nonlinearly related to the change in the absorption coefficient, $δ μ_{a}$ , because $δ Φ$ depends on $δ μ_{a}$ . If it is assumed that the fluence remains unaffected by the absorption perturbation (which certainly does not hold in all regimes in which PA imaging may be applied, but may be useful for very weak absorbers), then the second term in Eq. (33) can be neglected to give the linearized version:

Eq. (34)

δ H (x) \approx δ μ_{a} (x) Φ_{0} (x) .

This immediately shows why difference imaging will never remove the effect of the fluence even for very small absorption perturbations: if $δ H$ is considered to be the difference between two images, $δ H = H_{2} - H_{1}$ , then—while it depends on the difference in absorption coefficient, $δ μ_{a}$ , as expected—the effect of the fluence still remains. Interestingly, using a ratio of images under the same assumption that the fluence remains unchanged gives an image proportional to the change in absorption coefficient from which the effect of the fluence has been removed:

H_{2} / H_{1} \approx (μ_{a} + δ μ_{a}) Φ_{0} / μ_{a} Φ_{0} = 1 + (δ μ_{a} / μ_{a}) .

When the fluence due to the homogeneous part of the absorption, $Φ_{0}$ , is known, Eq. (34) can be inverted trivially to obtain an estimate for the absorption perturbation $δ μ_{a} (x) = δ H (x) / Φ_{0} (x)$ . Ripoll and Ntziachristos⁸⁰ modeled this unperturbed fluence with a Green’s function based on the assumptions that (1) the light emanates from a point source embedded in the scattering medium and (2) the fluence distribution obeys the DA (Sec. 2.2.2). To calculate this background fluence, the scattering coefficient as well as the homogeneous part of the absorption coefficient must be known. They assumed the scattering is homogeneous and state that an “underlying assumption is that of insensitivity to variations of tissue scattering properties” which is justified, they say, by the “experimental demonstration of high-quality images, obtained in vivo, even though the scattering variations were not explicitly accounted for.” However, while it is certainly true that high quality qualitative PA images can be obtained without considering the scattering perturbations, it is not obvious that quantitative estimates made with those images will be similarly unaffected. (Note that Ripoll and Ntziachristos embedded this linearized optical reconstruction within an acoustic image reconstruction algorithm.)

Rather than studying the perturbed problem, a first-order correction to the light fluence can be made by simply dividing $H$ by $Φ$ calculated assuming homogeneous (and known) background optical properties.⁸¹^–⁸³ A similarly restrictive assumption of nonchanging fluence underlies the suggestion that over a small wavelength range the fluence may be invariant to changes in the wavelength, so “measurements (images) at one wavelength can be explicitly used to normalize for photon intensity heterogeneity in tissues at the other wavelengths.”⁸⁴ Earlier, Kiser et al.⁸⁵ obtained PAT images of tumors implanted in murine mammary fat pads at 758 nm (peak Hb absorption) and 798 nm (isosbestic point), and by making the assumption that the fluence distribution does not change between the two wavelengths, they obtain images estimating the “oxygenated component” and “deoxygenated component,” although no scale (relative or absolute) is given.

There are two difficulties with this approach. The first is that the background properties of the target may not be known, so $Φ_{0}$ cannot be calculated, and the second and more important is that the assumption that the fluence remains unchanged when the optical absorption changes will only hold for very small perturbations about the background values. The latter problem is slightly alleviated with the method described in the next section.

5.1.2.

Born approximation

Assumptions similar to Sec. 5.1.1 above—known homogeneous scattering, known homogeneous background absorption, and DA for the light transport—are made by Zemp.⁸⁶ However, rather than assuming the fluence remains unchanged following an absorption perturbation, a correction to the fluence is made. When the fluence obeys the DA, the perturbed fluence $δ Φ$ due to the absorption perturbation $δ μ_{a}$ also obeys a diffusion equation:

Eq. (35)

(μ_{a, 0} - D_{0} \nabla^{2}) δ Φ (x) = - δ μ_{a} (x) Φ (x),

where

Φ

is the unknown fluence and

D_{0} = {[3 (μ_{a, 0} + μ_{s}^{'})]}^{- 1}

is the optical diffusion coefficient in the unperturbed state. Using the free-space Green’s function,

G_{0}

, Eq. (16), we may write the solution to Eq. (35) as

Eq. (36)

δ Φ (x) = - \int_{R^{n}} G_{0} (x, x^{'}) δ μ_{a} (x^{'}) Φ (x^{'}) d x^{'} .

By replacing the unknown fluence,

Φ

, with its unperturbed value,

Φ_{0}

(Born approximation), a linear expression relating the absorbed energy perturbation,

δ H

, and the absorption perturbation,

δ μ_{a}

, results in:

Eq. (37)

δ H (x) \approx δ μ_{a} (x) Φ_{0} (x) - μ_{a, 0} \int_{R^{n}} G_{0} (x, x^{'}) δ μ_{a} (x^{'}) Φ_{0} (x^{'}) d x^{'} .

As it is linear, this equation may be discretized and written as a matrix equation relating vectors of $δ H$ and $δ μ_{a}$ , which can then be inverted with any suitable method.

Interestingly, Zemp⁸⁶ uses a ratio data type constructed from two images by dividing one by the other. By using two images of the same target obtained with two different optical sources (same wavelength but illuminating from different places or with different patterns), this has the effect of removing any dependence on the PA efficiency and other instrumentation-related scaling factors. Measurements of the fluence at the boundary were also included as additional data in the reconstruction (compare Sec. 5.1.6). The extension of this idea to absorption and scattering estimation is described in Sec. 5.2.1.

Although these linearized inversions may have the computational advantages that a known solution such as a Green’s function can be used for the unperturbed fluence, and the inversion is straightforward, they are only useful when the absorption perturbations are small (even when using the Born approximation) and the background coefficients are known.

5.1.3.

Fixed-point iteration

When the absorption cannot be modeled accurately as a homogeneous background plus a small perturbation, a nonlinear inversion scheme is required. Cox et al.⁸⁷^,⁸⁸ rearranged Eq. (32) to obtain a fixed-point iteration for $μ_{a} (x)$ that makes no linearizing assumptions:

Eq. (38)

μ_{a}^{(n + 1)} (x) = \frac{H (x)}{Φ^{(n)} [x, μ_{a}^{(n)} (x)]},

where

Φ^{(n)}

, the fluence estimate at iteration

n

, is calculated by a numerical model from the latest absorption coefficient estimate,

μ_{a}^{(n)}

. If implemented crudely, problems may arise wherever there is noise in the data and the fluence takes small values. To reduce this instability, the fluence in the denominator may be replaced with the term

Φ^{(n)} + σ

, where

σ

is a regularization parameter chosen by the user and dependent on the signal-to-noise ratio. Alternatively, it may be possible to restrict the region over which the inversion is performed to regions where the fluence estimate is above some noise-dependent threshold value. Cox et al.⁸⁷^,⁸⁸ demonstrated this scheme with numerical examples in 2-D using a finite-element model of the DA, although it will also apply in 3-D and does not depend on this particular choice of light model or on any particular illumination geometry. Indeed, the same scheme has been used with the RTE⁸⁹ (Sec. 5.3). Jetzfellner et al.⁹⁰ note that when the value used for the scattering is wrong, so the fluence from the model is biased, the estimates of the absorption begin to diverge after several iterations. Nevertheless, it has been used with experimental data to estimate absorption coefficients.⁹¹^–⁹³

5.1.4.

Noniterative inversion

Banerjee et al.⁹⁴ realized that the absorbed energy $μ_{a} Φ$ appears directly in the DA model, Eq. (14), and the following equation can therefore be written for the fluence:

Eq. (39)

\nabla \cdot (D \nabla) Φ = H - q_{0} .

By assuming that the optical scattering coefficient is independent of absorption,

D \approx {(3 μ_{s}^{'})}^{- 1}

, and solving Eq. (39) numerically, the absorption coefficient may be estimated by

μ_{a} = H / Φ

. When the optical diffusion is considered to depend on the absorption, the estimate of the absorption iteration could be improved by use of iteration.

5.1.5.

Least-squares minimization

The inversion for absorption has been tackled by formulating it as a least-squares minimization (Sec. 2.3.2) of the functional

Eq. (40)

\underset{μ_{a}}{argmin} ε = \frac{1}{2} {‖ H^{model} (μ_{a}) - H^{obs} ‖}^{2} + P (μ_{a}),

where

H^{obs}

is the measured image, and the penalty functional,

P

, is typically chosen to enforce Tikhonov or Total Variation regularization. Jiang and colleagues⁹⁵ applied the Gauss–Newton method, which, as it requires the calculation and storage of a Jacobian matrix as well as the inversion of the Hessian estimate, quickly becomes computationally intensive for large or high resolution images (see Sec. 5.2.3). Greater computational efficiency could be achieved with gradient-only methods, as the functional gradient vector can be calculated quickly by an adjoint model⁹⁶ (see Sec. 5.2.5), although more iterations will usually be required for convergence.

5.1.6.

Correction based on diffuse optical tomography

As the absorbed energy density $H = μ_{a} Φ$ , if the fluence, $Φ$ , were known, then the inversion for $μ_{a}$ would become trivial. This is the view taken in this and the next section. If there were some complementary imaging modality that could provide an image of the fluence with the same spatial resolution with which the PA image provides the absorbed energy density, then the absorption coefficient could be found just by dividing the PA image by the fluence image. Unfortunately, no such option is currently available. The fluence can be estimated from optical boundary measurements by diffuse optical tomography (DOT), but the high spatial frequencies are lost because of the diffusive nature of the light propagation, and only a low frequency approximation to the fluence can be obtained this way. Despite this limitation it has been proposed as a method for obtaining a fluence estimate for quantitative PA imaging.⁹⁷^,⁹⁸ In situations where the actual fluence has no high frequency components, this method will recover accurate values for the optical coefficients, so the skill to applying this approach successfully is knowing, in advance, whether the unknown optical coefficients are likely to give rise to a fluence that is sufficiently smooth for this method to work. This may be the case when only small weak absorbers are present, or when the absorption and scattering coefficients themselves are smooth, although in this case it may be that the estimates can be obtained by DOT directly. [A procedure in which the localization is achieved through “acoustic tagging” of light (i.e., acousto-optics) has also been proposed.⁹⁹]

5.1.7.

Sparse basis decomposition

By taking the log of Eq. (32), the absorption and fluence can be separated: $\log (H) = \log (μ_{a}) + \log (Φ)$ . Rosenthal et al.¹⁰⁰ propose using two different sparse representations for $\log (μ_{a})$ and $\log (Φ)$ :

Eq. (41)

\log [H (x)] = \sum_{n} a_{n} ϕ_{n} (x) + \sum_{m} b_{m} ψ_{m} (x),

where {

ϕ_{n}

} and {

ψ_{m}

} are the two sets of basis functions, and

a_{n}

and

b_{m}

are the coefficients. If {

ϕ_{n}

} and {

ψ_{m}

} can be chosen so that the absorption is sparse in one basis (only a few coefficients are large) and the fluence is sparse in the other, then this method may be able to recover both the absorption and fluence. If suitable basis functions are known a priori, then the aim is to find the smallest number of nonzero coefficients,

a_{n}

,

b_{m}

that satisfy

{‖ H - H^{obs} ‖}_{2} < ϵ

, where

ϵ

is small, so that the absorption can be reconstructed from

μ_{a} (x) \approx \exp [\sum a_{n} ϕ_{n} (x)]

where only the first few terms are included. This method applies to those cases where the fluence and the absorption coefficient distributions are sufficiently different types of function that they cannot both be represented sparsely in the same basis but can each be represented sparsely in different bases. In addition, it is necessary that the basis sets that can distinguish between the two are known. [They might be chosen based on prior knowledge of the absorption distribution. For example, Rosenthal et al. assume the fluence can be represented sparsely in a Fourier basis (smooth), and the absorption, in a Haar wavelet basis (stepped).] Even when these conditions are satisfied, a further limitation of this method is that it cannot distinguish between the components of absorption and fluence that are uniform across the domain, so it can only recover

μ_{a}

to within an additive constant. Nevertheless, its splitting of the absorption and fluence, its lack of dependence on knowledge of the scattering, and its independence from a light model are attractive.

5.1.8.

Fluence-dependent chromophores

A quite different approach would be to measure the fluence with a specially designed PA contrast agent as a “probe” that can indicate the value of the fluence at a point within the tissue. By using a contrast agent that suddenly becomes more (or less) absorbing at a particular fluence threshold, it would be possible to estimate the fluence by noticing where in the image the contrast agent stops absorbing as the illumination intensity is increased. Gold nanorods that undergo a shape change at a known threshold fluence have been suggested as a possible contrast agent.¹⁰¹

5.2.

Inversions for Absorption and Scattering: Diffusion Approximation

Most of the methods in Sec. 5.1 assumed that the optical scattering (or optical diffusion) coefficient distribution was known. However, this will not always be the case. Methods to measure the scattering noninvasively have been proposed,¹⁰²^,¹⁰³ although they either only recover the scattering close to the tissue surface or give an averaged bulk value, so techniques to invert $T (μ_{a}, μ_{s})$ are needed in general. Under the DA the scattering is described by the reduced scattering coefficient $μ_{s}^{'}$ , or equivalently the optical diffusion coefficient $D$ , and Eq. (32) becomes one or the other of

Eq. (42)

H (x) = μ_{a} (x) Φ [x, μ_{a} (x), μ_{s}^{'} (x)], H (x) = μ_{a} (x) Φ [x, μ_{a} (x), D (x)] .

As discussed in Sec. 2.3.3, one difficulty with estimating both absorption and scattering from one PA image is that there will, in general, be more than one possible solution, only one of which is the correct solution. Any method that aims to recover both coefficients must therefore include some prior information about their distributions, or the light fluence perhaps, which can restrict the solution space sufficiently to overcome this nonuniqueness.

5.2.1.

Born approximation

The linearized approach of Sec. 5.1.2 can be extended to recovering scattering and absorption perturbations simultaneously.⁴⁵ The equation for the change in the fluence now contains an additional term for the optical diffusion perturbation (related to the scattering),

Eq. (43)

{(μ}_{a, 0} - \nabla \cdot D_{0} \nabla) δ Φ (x) = [\nabla \cdot δ D (x) \nabla - δ μ_{a} (x)] Φ (x),

so the total fluence can be written as

Eq. (44)

Φ = Φ_{0} - \int_{R^{n}} G_{0} (x, x^{'}) δ μ_{a} (x^{'}) Φ_{0} (x^{'}) d x^{'} - \int_{R^{n}} δ D (x^{'}) \nabla G_{0} (x, x^{'}) \cdot \nabla Φ_{0} (x^{'}) d x^{'} .

The expression for the change in the absorbed energy therefore contains an additional term

Eq. (45)

δ H (x) \approx δ μ_{a} (x) Φ_{0} (x) - μ_{a, 0} \int_{R^{n}} G_{0} (x, x^{'}) δ μ_{a} (x^{'}) Φ_{0} (x^{'}) d x^{'} - μ_{a, 0} \int_{R^{n}} δ D (x^{'}) \nabla G_{0} (x, x^{'}) \cdot \nabla Φ_{0} (x^{'}) d x^{'},

where the unknown fluence inside the integrals has been replaced with the fluence for the homogeneous case,

Φ_{0}

, and the identity

\int_{Ω} G_{0} \nabla \cdot (δ D \nabla Φ) d Ω = \int_{\partial Ω} δ D G_{0} \nabla Φ d (\partial Ω) - \int_{Ω} δ D (\nabla G_{0} \cdot \nabla Φ) d Ω

has been used assuming the surface integral tends to zero as the radius of the surface

\partial Ω

of the domain

Ω

tends to infinity. In fact, Shao et al.⁴⁵ do not use this approximation of the absorbed energy perturbation directly but, rather, use the ratio of two images,

H^{(1)}

and

H^{(2)}

, obtained with different illumination patterns, that is, different positions of the optical source. As neither the absorption coefficient nor the PA efficiency depend on the illumination, they cancel out leaving a ratio of fluences, each of which can be written with the Born approximation above:

Eq. (46)

\frac{H^{(1)}}{H^{(2)}} = \frac{Φ^{(1)}}{Φ^{(2)}} \approx \frac{Φ_{0}^{(1)} - \int G_{0} (x, x^{'}) δ μ_{a} (x^{'}) Φ_{0}^{(1)} (x^{'}) d x^{'} - \int δ D (x^{'}) \nabla G_{0} (x, x^{'}) \cdot \nabla Φ_{0}^{(1)} (x^{'}) d x^{'}}{Φ_{0}^{(2)} - \int G_{0} (x, x^{'}) δ μ_{a} (x^{'}) Φ_{0}^{(2)} (x^{'}) d x^{'} - \int δ D (x^{'}) \nabla G_{0} (x, x^{'}) \cdot \nabla Φ_{0}^{(2)} (x^{'}) d x^{'}} .

By discretizing the integrals and rearranging, this expression can be rearranged into a matrix form

Eq. (47)

A [\begin{matrix} δ μ_{a} \\ δ D \end{matrix}] = b,

where the vector

b

consists of a combination of measured images and modeled fields for the homogeneous case,

p_{0}^{(1)} Φ_{0}^{(2)} - p_{0}^{(2)} Φ_{0}^{(1)}

, and the elements of matrix

A

contain products of the images and the (discretized) integrals. Matrix

A

can be inverted with any appropriate matrix inversion algorithm, with regularization as required to mitigate the effect of noise in the images. (Even though two images are being used with different illuminations to overcome the absorption-scattering nonuniqueness, the fact that the light propagation is diffusive in nature will result in some ill-posedness in the inversion.) Note that although the quotient of two images was used in the derivation above, it is never actually necessary to calculate that quotient, so the issue of dividing by small numbers is avoided. Interestingly, although the PA efficiency is canceled out when forming the ratio in Eq. (46), it can be recovered afterwards from the image of the initial pressure,

p_{0}

, by calculating the true fluence for one of the illumination patterns with the (now) known absorption and scattering. The drawbacks of the linearization remain, however: the Born approximation is only accurate for small perturbations in the absorption and scatter, and the values for the background coefficients need to be known in advance.

5.2.2.

Direct “zero-divergence” method

Bal et al.⁴⁴^,¹⁰⁴ have devised a method for solving the optical inversion which requires neither a linearization nor the iterative updating of the fluence. It is a direct rearrangement of the DA assuming that two (or more) PA images, $H_{1} = μ_{a} Φ_{1}$ , $H_{2} = μ_{a} Φ_{2}$ , obtained with two (or more) different illumination patterns are available. The problem is restated in an equivalent but slightly different form from most of this review, so rather than trying to estimate the optical coefficients $μ_{a}$ and $μ_{s}^{'}$ (or $μ_{a}$ and $D$ ), the unknown parameters are defined as:

Eq. (48)

μ_{1} ≔ \frac{\sqrt{D}}{μ_{a}}, μ_{2} ≔ - \frac{\nabla^{2} (\sqrt{D})}{\sqrt{D}} - \frac{μ_{a}}{D} .

Multiplying the diffusion equation for the fluence $Φ_{1}$ by $Φ_{2}$ , and vice versa, and subtracting gives $Φ_{2} [μ_{a} Φ_{1} - (\nabla \cdot D \nabla) Φ_{1}] - Φ_{1} [μ_{a} Φ_{2} - (\nabla \cdot D \nabla) Φ_{2}] = 0$ , which can be rearranged to $\nabla \cdot [D Φ_{1}^{2} \nabla (Φ_{2} / Φ_{1})] = 0$ . Recalling $H = μ_{a} Φ$ and defining the vector field

Eq. (49)

β = H_{1} \nabla H_{2} - H_{2} \nabla H_{1} = H_{1}^{2} \nabla (\frac{H_{2}}{H_{1}}) = μ_{a}^{2} Φ_{1}^{2} \nabla (\frac{Φ_{2}}{Φ_{1}}),

leads to the following concise relationship:

Eq. (50)

\nabla \cdot (μ_{1}^{2} β) = 0.

The procedure to reconstruct images of $μ_{a}$ and $D$ is then the following: (1) Construct the vector field $β$ from the two measured fields $H_{1}$ and $H_{2}$ by using Eq. (49). (2) Solve the transport equation, Eq. (50), for $μ_{1} (x)$ that is, find the scalar field that makes the divergence of $μ_{1}^{2} β$ vanish. This is nontrivial and will not be possible for all fields $β$ , as discussed below. (3) Find the second unknown scalar field through use of one of the images $H_{1}$ or $H_{2}$ and the equation $μ_{2} = - \nabla^{2} (μ_{1} H_{1, 2}) / (μ_{1} H_{1, 2})$ , which can be found by rearranging the definition in the second of Eqs. (48). (4) Solve $(\nabla^{2} + μ_{2}) \sqrt{D} = - 1 / μ_{1}$ for $\sqrt{D}$ , which also comes from Eqs. (48). (5) The final step is to obtain $μ_{a} = \sqrt{D} / μ_{1}$ from the first of Eqs. (48).

The key to this method is part (2), finding a solution to Eq. (50). This may be done in a number of ways. Bal and Ren¹⁰⁴ discretize Eq. (50) and find $μ_{1}$ by using a Bregman iteration to minimize a functional consisting of the norm of the divergence of $μ_{1}^{2} β$ and a regularizing total variation penalty term. They used a low pass filter to ameliorate the difficulties in the calculation of $β$ caused by dividing by small values. A more visual, but less practical, method of solving Eq. (50) was proposed by Bal and Uhlmann.⁴⁴ Using the chain rule, we may rewrite Eq. (50) as

Eq. (51)

2 β \cdot \nabla μ_{1} + γ μ_{1} = 0,

where

γ = \nabla \cdot β = H_{1} \nabla^{2} H_{2} - H_{2} \nabla^{2} H_{1}

. If

s

is a parametric coordinate along a field line (integral curve) of the vector field

β

that passes through position

x

then

\partial x / \partial s = β [x (s)]

. Using the chain rule shows that

β \cdot \nabla μ_{1} = \partial μ_{1} / \partial s

and so

\partial μ_{1} / \partial s = - γ μ_{1} / 2

. For a given point,

x^{*}

,

μ_{1} (x^{*})

can be calculated by integrating

γ (x)

along the field line of

β

that passes through

x^{*}

, starting from where that particular field line reaches the edge of the domain:

Eq. (52)

μ_{1} [x^{*} (s)] = μ_{1}^{0} \exp (- \frac{1}{2} \int_{0}^{s} γ [x (s^{'})] d s^{'}),

where

μ_{1}^{0}

, the value of

μ_{1} [x (s = 0)]

on the boundary of the domain, is assumed known. As

μ_{1} = \sqrt{D} / μ_{a}

, both

D

and

μ_{a}

must be known on the boundary. While not necessarily leading to a practical numerical approach, this makes clearer what condition

β

must satisfy in order that the coefficients can be recovered: every point in the domain must be connected to the edge of the domain by a field line of

β

. As

β

depends only on the measured images, and therefore on the illumination patterns, for this approach to work it is necessary to choose the illumination patterns that result in a

β

that has field lines connecting every point in the domain to the edge (or equivalently of a form such that multiplying by it by a scalar field can make its divergence zero). While there are no exact rules to guide the choice of illuminations, the examples in Ref. 104 suggest that in practice this is not a significant limitation and sensibly chosen illuminations (e.g., illumination on the top and bottom surfaces of the sample) will suffice.

5.2.3.

Least-squares minimization

The least-squares approach, introduced in Sec. 5.1.5 for recovering the absorption alone, can be extended to the case when both the scattering and absorption coefficients are unknown and must be recovered from the data,⁹⁶ provided the nonuniqueness is overcome by using multiple wavelengths or multiple illumination patterns. To simplify the notation we assume that the domain of interest has been discretized such that the unknown absorption and diffusion coefficients can be represented by $N \times 1$ column vectors, $μ_{a}$ , $D \in R^{N}$ , for example, representing the values on the elements of a mesh. The measured image is stored in an $M \times 1$ vector, for example, nodal values of some mesh, so $H^{obs} \in R^{M}$ . The task of recovering the unknown coefficients ( $μ_{a}$ , $D$ ) becomes the task of adjusting the coefficients in order to minimize the functional

Eq. (53)

\underset{μ_{a}, D}{argmin} ε = \frac{1}{2} {‖ H (μ_{a}, D) - H^{obs} ‖}^{2} + P (μ_{a}, D) .

In this section, the DA will be used as the light model, but the general approach can be applied to other light models. The idea is to iteratively improve the estimates of the unknowns starting from initial guesses ( $μ_{a}^{(0)}$ , $D^{(0)}$ ) by choosing new values $μ_{a}^{(k + 1)} = μ_{a}^{(k)} + δ μ_{a}^{(k)}$ , $D^{(k + 1)} = D^{(k)} + δ^{(k)}$ so that at each iteration the value of $ε [μ_{a}^{(k + 1)}, D^{(k + 1)}]$ is smaller than $ε [μ_{a}^{(k)}, D^{(k)}]$ . The hope is that eventually the minimum of the multidimensional surface of $ε (μ_{a}, D)$ will be reached, thereby giving the values of $(μ_{a}, D)$ at which $H (μ_{a}, D)$ is the best match to $H^{obs}$ (given the constraints $P$ ).

To help find the minimum more quickly, the gradients and second derivatives of $ε$ (related to its curvature) are usually used. The gradients of $ε$ with respect to the optical coefficients, $g_{a}$ , $g_{D} \in R^{N}$ can be written as

Eq. (54)

g_{a} = \frac{\partial ε}{\partial μ_{a}} = J_{a}^{T} (H - H^{o b s}) + \frac{\partial P}{\partial μ_{a}}, g_{D} = \frac{\partial ε}{\partial D} = J_{D}^{T} (H - H^{o b s}) + \frac{\partial P}{\partial D},

where the Jacobian matrices,

J_{a}

,

J_{D} \in R^{M \times N}

are given by

Eq. (55)

J_{a} = \frac{\partial H}{\partial μ_{a}} = Φ δ + μ_{a} (\frac{\partial Φ}{\partial μ_{a}}), J_{D} = \frac{\partial H}{\partial D} = μ_{a} (\frac{\partial Φ}{\partial D}),

where

δ

is the Dirac delta function

δ (x - x^{'})

in the continuous case, and a Kronecker delta

δ_{i j}

in the discrete case and comes from differentiating

\partial μ_{a} (x) / \partial μ_{a} (x^{'})

. (Similarly with

\partial D (x) / \partial D (x^{'})

below.) The Hessian matrix of second-order derivatives of

ε

,

H \in R^{2 N \times 2 N}

, may be written as

Eq. (56)

H = \frac{\partial^{2} ε}{\partial μ_{i} \partial μ_{j}} \approx J^{T} J, J = [\begin{matrix} J_{a} & J_{D} \end{matrix}] .

So far no explicit model for the light fluence has been introduced, but in order to evaluate the Jacobians, a model is required. If the DA is used, then $\partial Φ / \partial μ_{a}$ and $\partial Φ / \partial D$ can be calculated from the sensitivity equations:

Eq. (57)

(μ_{a} - \nabla \cdot D \nabla) \frac{\partial Φ}{\partial μ_{a}} = - Φ δ, (μ_{a} - \nabla \cdot D \nabla) \frac{\partial Φ}{\partial D} = \nabla \cdot δ \nabla Φ .

When using a particular discretization scheme, such as finite elements, the terms on the right-hand side may be calculated directly by taking the derivative of the discretized forward operator. In the Gauss–Newton method, the update to the unknown parameters is calculated by

Eq. (58)

[\begin{matrix} δ μ_{a} \\ δ D \end{matrix}] = - α {(J^{T} J)}^{- 1} [\begin{matrix} g_{a} \\ g_{D} \end{matrix}], g = [\begin{matrix} g_{a} \\ g_{D} \end{matrix}],

where the scalar

α

is chosen by linesearch to minimize

ε

at each step. One advantage of the Gauss–Newton approach is that there are many well-developed minimization algorithms available and typically only a few iterations are required.

5.2.4.

Multiwavelength minimizations

In the multiwavelength case, the unknowns are now the chromophore concentration distributions, $C_{k}$ , and the scattering and the data are multiwavelength images $H^{obs} (x, λ)$ . The error functional therefore becomes

Eq. (59)

\underset{C_{k}, μ_{s}^{'}}{argmin} ε_{λ} = \frac{1}{2} {‖ H (C_{k}, μ_{s}^{'}) - H^{obs} ‖}^{2} + P (C_{k}, μ_{s}^{'}),

where

‖ \cdot ‖

is now the norm over position and wavelength. The gradients of

ε_{λ}

and

H

with respect to the chromophores are simply related to those for the absorption coefficients. Using the chain rule and Eq. (4) gives

Eq. (60)

\frac{\partial}{\partial C_{k}} = \frac{\partial μ_{a}}{\partial C_{k}} \frac{\partial}{\partial μ_{a}} = α_{k} \frac{\partial}{\partial μ_{a}} .

If the scattering is assumed to depend on wavelength as $μ_{s}^{'} \approx a (x) λ^{- b}$ , then the gradients can also be simply calculated:

Eq. (61)

\frac{\partial}{\partial a} = \frac{\partial μ_{s}^{'}}{\partial a} \frac{\partial}{\partial μ_{s}^{'}} = λ^{- b} \frac{\partial}{\partial μ_{s}^{'}} .

Using these expressions, Cox et al.⁴² inverted for the concentrations and scattering distributions with the Gauss–Newton method. Their numerical study used multiwavelength data in 2-D inversions with the assumptions that (1) the absorption spectra of the constituent chromophores, $α_{k} (λ)$ , were known and (2) the scattering wavelength dependence $b$ was known. This proved sufficient to remove the nonuniqueness (without the use of multiple illuminations) and allow the concentration and scattering amplitude $a (x)$ to be recovered.

5.2.5.

Gradient-based minimizations

The functional gradients in Eq. (54) were calculated from the Jacobian matrices, which may be large. However, in the DA case there is a fast and memory-efficient way to calculate the gradients with an adjoint model. Rather than having to solve both Eq. (57) once per column of the Jacobian matrix, that is, $2 N$ times per iteration step, it is only necessary to run the forward and the adjoint models once each:⁹⁶

Eq. (62)

(μ_{a} - \nabla \cdot D \nabla) Φ = q, (μ_{a} - \nabla \cdot D \nabla) Φ^{*} = μ_{a} (H - H^{o b s}) .

From these two fields the gradients can be calculated directly:

Eq. (63)

g_{a} = Φ (H - H^{o b s}) - Φ Φ^{*}, g_{D} = - \nabla Φ \cdot \nabla Φ^{*} .

Using a gradient-based split Bregman algorithm, Gao et al.¹⁰⁵^,¹⁰⁶ recently demonstrated that a gradient-based approach combined with multisource data can successfully recover both absorption and scattering. They recovered piecewise constant absorption and scattering distributions in 2-D and 3-D by using a total variation regularization term in the inversion. Bal and Ren¹⁰⁴ performed similar inversions with the BFGS algorithm with multisource and multiwavelength data. A large number of variations on the basic theme of using only gradients are possible, but gradient-based approaches such as this certainly look promising as candidates for practical applications. The main drawback remaining is the reliance (in these examples) on the DA, which is rarely true for the whole region of interest in PA applications.

5.3.

Inversions for Absorption and Scattering: RTE

Most of the inversion methods described so far have relied on the DA, which is only accurate at distance further than a few mean free paths away from the surface being illuminated. In many PAT applications this surface region is of interest. For these cases it is necessary to look to more accurate models of light transport, such as the RTE (Sec. 2.2.1), Monte Carlo models (Sec. 2.2.4), or the delta-Eddington approximation (Sec. 2.2.3). Several authors have looked at using the RTE for qPAI. Bal et al.¹⁰⁷ show that in the transport regime, if the unscattered and singly scattered components of the fluence can be separated from the multiply scattered components, then not only the absorption and scattering but also the anisotropy factor, $g$ , can be recovered. Unfortunately, this decomposition is impractical for most biological media of interest. Yao et al.⁸⁹ used the RTE with the fixed-point iteration of Sec. 5.1.3, which assumes that the scattering is known, and compared it with the same inversion using the DA. They showed that the RTE gave more accurate results for both simulations and a simple experimental phantom study. Cox et al.¹⁰⁸ showed that the Gauss–Newton method (Sec. 5.2.3) with multiple illumination data can be used with the RTE to recover good estimates of both absorption and scattering in a simple numerical 2-D phantom.

6. Discussion

6.1.

Multiwavelength Versus Single-Wavelength Inversions

As mentioned in Sec. 2.1.1, Eqs. (4) and (5) suggest a two-stage inversion strategy: first recover the absorption coefficients $μ_{a} = T^{- 1} (H)$ then find the concentrations $C_{k} = L_{λ}^{- 1} (μ_{a})$ , whereas Eq. (6) suggests a single inversion, $C_{k} = T_{λ}^{- 1} (H)$ . There may be advantages to the latter strategy of recovering chromophore concentrations directly. First, when there are more wavelengths (more images) than chromophores, the number of unknown concentrations will be fewer than the number of unknown absorption coefficients (one per wavelength). As the problem can grow very large when considering 3-D images (Sec. 2.3.4), this may be an important consideration. Second, prior knowledge of spectral information can be incorporated, which can remove the absorption-scattering nonuniqueness. Third, it may require less measured data: the most widely studied way to remove the absorption-scattering nonuniqueness when inverting for absorption coefficient at a single wavelength is to use images recorded with different illumination directions, so two or more images per wavelength are required in this approach. If the nonuniqueness can be removed through the use of spectral priors, then only one image per wavelength is necessary, that is, half the amount of measured data. A possible advantage of the two-stage inversion strategy is that the single-wavelength inversions for $μ_{a}$ are independent and so can be performed straightforwardly in parallel.

6.2.

Nonlinearity

The nonlinearity of the optical inversion caused by the dependence of the fluence on absorption in Eq. (2) should not be confused with nonlinearities that appear in at least two other settings related to PA imaging. Wang et al.¹⁰⁹ observed the nonlinear dependence of a PA signal amplitude on absorption coefficient because of the bandlimiting effect of the detector, suggesting that broadband detectors should be used for quantitative PA imaging. Paltauf et al.¹¹⁰ point out that the ultrasonic attenuation will also have a bandlimiting effect that could result in a similar nonlinear dependence. Despite using a broadband detector, Karabutov et al.¹¹¹^,¹¹² noticed the nonlinear dependence of the PA signal on absorption coefficient for a beam incident on a nonscattering medium. As the absorption coefficient increases, the characteristic length of the absorbed energy distribution, $1 / μ_{a}$ reduces until it is less than the distance the PA wave can travel during the laser pulse, and the stress confinement condition no longer holds. Danielli et al.¹¹³ used the nonlinear PA response from an absorber owing to optical saturation to determine relaxation times on the picosecond timescale.

6.3.

PA Generation Efficiency

The PA generation efficiency is the thermodynamic constant connecting the acoustic and the optical inversions, labeled $\hat{Γ}$ in Figs. 1 and 2. It plays a critical role in the generation of PA signals and is therefore relevant to QPAI. Following the absorption of a photon, the energy in the excited absorbing molecule will be redistributed between the various vibrational and other modes of that molecule and neighboring molecules in a variety of processes grouped under the term thermalization. The localized increase of heat will lead to perturbations in the local temperature $T^{'}$ , pressure $p^{'}$ , and density $ρ^{'}$ . The efficiency of the PA effect indicates the size of the pressure increase that accompanies thermalization, which will depend on the local microenvironment of the absorbing molecule.

6.3.1.

Grüneisen Parameter

In an absorbing fluid, the PA efficiency $\hat{Γ}$ can be related to the Grüneisen parameter, $Γ$ , which is a thermodynamic property of a material. Consider the case in which the absorbed power density, $H$ , the rate at which the optical energy is absorbed per unit volume per unit time, is separable so $H (x, t) = H (x) f (t)$ , where $H$ is the absorbed energy density and $f (t)$ is a temporal shape function with unit integral. Some of the absorbed energy will be used doing work (changing the mass density), and the remainder will become heat and lead to a change in temperature. The balance between these two will determine the size of the PA effect: the less energy used to do work, the greater will be the increase in temperature and corresponding increase in pressure.

In fluid with homogeneous thermodynamic properties, the rates of change of the local temperature pressure, and density may be related by

Eq. (64)

\frac{\partial p^{'}}{\partial t} = (\frac{1}{ρ k_{T}}) \frac{\partial ρ^{'}}{\partial t} + (\frac{β}{k_{T}}) \frac{\partial T^{'}}{\partial t},

where

ρ

is the ambient mass density,

k_{T}

is the isothermal compressibility, and

β

is the volume thermal expansivity of the fluid. For a sufficiently short pulse,

f (t)

, the density will not have time to decrease and

| \partial ρ^{'} / \partial t | ≪ | \partial T / \partial t |

, which is sometimes called stress confinement. Under this isochoric condition the rate of increase in pressure can be written

\partial p^{'} / \partial t = [(β / k_{T}) \partial T^{'}] / \partial t

. The rate at which the temperature increases is given by

\partial T^{'} / \partial t = H / ρ C_{v}

, where

C_{v}

is the specific heat capacity at constant volume. The rate at which the local pressure rises can now be related to the absorbed optical power density by

Eq. (65)

\frac{\partial p^{'}}{\partial t} = (\frac{β}{ρ C_{v} k_{T}}) H (x, t) = Γ H (x) f (t),

where

Γ = β / (ρ C_{v} k_{T})

is the Grüneisen parameter. (An equivalent expression is

Γ = β c_{0}^{2} / C_{p}

, where

C_{p}

is the specific heat capacity at constant pressure,

c_{0} = \sqrt{B_{s} / ρ}

is the sound speed, and

B_{s} = (C_{p} / C_{v}) / k_{T}

is the isentropic bulk modulus.) The total increase in pressure (referred to as the initial acoustic pressure distribution) can be obtained by integrating over the duration of the optical pulse to get

Eq. (66)

p_{0} = Γ H .

Equation (66) is a special case of Eq. (1) that only has meaning when the absorber and the propagation medium have the same thermodynamic properties. Even with this proviso there will be conditions under which the PA efficiency will not equal the Grüneisen parameter. If there is significant radiative decay (e.g., fluorescence) following the absorption, or if the isochoric assumption $\partial ρ^{'} / \partial t ≪ \partial T^{'} / \partial t$ does not hold, then $\hat{Γ} < Γ$ . The latter will be true in pulses for which $H$ is small (energy is delivered to the tissue slowly), allowing the density to change on a timescale similar to the temperature.

6.3.2.

PA efficiency for particle suspensions

When the absorption is due to solid inclusions such as nanoparticles, with different material properties to the surrounding fluid, the PA efficiency cannot in general be simply related to the Grüneisen parameter of either material. The PA signal will arise from heating of the fluid surrounding the inclusions by heat conduction from the inclusions. In other words, rather than heating the fluid directly through the absorption of light, there is an intermediate stage: the light is absorbed by the solid inclusion and the heat is then transferred to the fluid.

The expansion of the particles themselves will be a second order effect because although the temperature rise of the inclusions themselves can be greater than for the fluid case (as the specific heat capacity tends to be lower in solids than fluids), the stress generated and expansion of the particles is small. This is because not only is the thermal expansivity is typically smaller in solids than fluids, but the density and bulk modulus are typically larger.

Whereas for an absorbing fluid the rate of heat deposition is the key factor determining the efficiency of the PA effect, in this case it is the rate at which the heat conducts into the fluid and not the rate at which it is deposited in the solid particles that is significant. The thermal resistance, rather than Grüneisen parameter, becomes the main physical parameter controlling the efficiency of the PA generation.¹¹⁴ It has been shown that this rate can be increased for gold nanoparticles by coating them with silica, thereby increasing $\hat{Γ}$ and the magnitude of the PA signal.¹¹⁵

6.3.3.

PA efficiency in quantitative PA imaging

Many of the inversion procedures in Secs. 3 through 5 assume that the PA efficiency is known and therefore take the absorbed energy density $H$ as a starting point, but there are exceptions. Zemp et al.,⁸⁶^,⁴⁵ Guo et al.,⁵² and Bal et al.⁴⁴^,¹⁰⁴ use ratios of two measurements, one result of which is to cancel out the PA efficiency. Another way around the problem of not knowing $\hat{Γ}$ is to recover it simultaneously with the chromophore concentrations and scattering. Bal and Ren¹⁰⁴ have shown in the single-wavelength case that when the light propagation obeys the DA, two measurements are insufficient to allow the absorption and scattering coefficients as well as the PA efficiency to be recovered, although Shao et al. do this in a linearized version of the problem.⁴⁵

6.3.4.

Tissue- or chromophore-dependent PA efficiency

Sometimes the PA efficiency for tissue is assumed to be a constant, so that the initial pressure distribution is directly proportional to the absorbed energy density with the same constant of proportionality throughout the tissue, $p_{0} (x) = \hat{Γ} H (x)$ . However, $\hat{Γ}$ may take different values in different tissue types: the Grüneisen parameter of blood has been measured to be⁶⁵ $Γ (blood) = 0.14$ , while in liver samples on average¹¹⁶ $Γ (liver) = 0.12$ , which is greater than 15% difference. (For comparison, $Γ (water) = 0.11$ .) If a contrast agent is used, it is even more likely that the PA efficiency will vary spatially, $\hat{Γ} = \hat{Γ} (x)$ . Putting this spatial dependence into Eqs. (1) and (6) gives an expression for the initial pressure distribution as

Eq. (67)

p_{0} = \hat{Γ} (x) Φ (x, λ) \sum_{k = 1}^{K} C_{k} (x) α_{k} (λ) .

Laufer et al.¹¹⁷ suggest there may be a different PA efficiency for each chromophore which depends on that chromophore’s concentration, ${\hat{Γ}}_{k} = {\hat{Γ}}_{k} [C_{k} (x)]$ . The initial pressure distribution then becomes

Eq. (68)

p_{0} = Φ (x, λ) \sum_{k = 1}^{K} {\hat{Γ}}_{k} C_{k} (x) α_{k} (λ) .

In this case, the algorithms of Secs. 3 to 5 could be used to estimate the product

{\hat{Γ}}_{k} C_{k}

for each chromophore, rather than the concentration

C_{k}

. If

{\hat{Γ}}_{k}

is known for the chromophores of interest, then its concentration can be recovered; it is not necessary to know

{\hat{Γ}}_{k}

for every chromophore present.

7. Conclusion

Many different inversion schemes have been presented for quantitative photoacoustic imaging, each of which makes different assumptions and will therefore have a different region of applicability. The differences are in several different directions; for example, between those schemes that assume a 1-D system and those that can incorporate the full complexity of 3-D; between those that are ad hoc, and therefore quite specific to a given scenario, and those of more generally applicability; between linearized and fully nonlinear models; between those that are computationally intensive and those whose efficiency will allow them to be extended to large scale inversions for high resolution images; and between those that assume the optical scattering is known and those for which it is inverted along with the chromophore concentrations. Some methods such as the single-point estimation of ${sO}_{2}$ in an isolated blood vessel are very specific but may well find a practically useful niche. Other methods, such as the multiwavelength gradient-based inversions, are more general and therefore look like more promising candidates for solving the complete problem of separating several spatially varying concentrations and scattering distributions.

Many of the methods described, particularly those in Sec. 6, have only been tested in simulation (and that often only in 2-D), and it is not at all clear which will be the most robust to the uncertainties that will be found in practice (such as in the “known” model parameters: detector positions, illumination beam dimensions and profile, absorption spectra, as well as to image artifacts and noise) or what other hurdles will become apparent when implementing them in practice. There is therefore a pressing need for these models to be tested under realistic conditions and with real data sets.

Others schemes, such as some of the single-point techniques that have been used to measure ${sO}_{2}$ , have been used for in vivo measurement but without sufficient prior testing in phantom experiments and simulations to show where their assumptions break down. More testing and rigorous validation with realistic phantoms would give greater confidence that the in vivo results are quantitatively correct. There is an urgent need to ensure that these techniques are valid, accurate, and robust in the settings in which they will be used, especially those that have been proposed for use in clinical applications.

Acknowledgments

The authors would like to thank Guillaume Bal, Anthony Harker, Yuroslav Kurylev, Sandy MacRobert, Stephen Price, Teedah Saratoon, Tanja Tarvainen, and Roger Zemp for helpful discussions. This work was supported by King’s College London and University College London Comprehensive Cancer Imaging Centre, Cancer Research UK & Engineering and Physical Sciences Research Council, in association with the Medical Research Council and Department of Health, UK.

References

1.

Photoacoustic Imaging and Spectroscopy, CRC Press(2009). Google Scholar

2.

P. Beard, “Biomedical photoacoustic imaging,” Interf. Focus, 1 (4), 602 –631 (2011). Google Scholar

3.

L. V. Wang, “Tutorial on photoacoustic microscopy and computed tomography,” IEEE J. Sel. Top. Quant. Electron., 14 (1), 171 –179 (2008). http://dx.doi.org/10.1109/JSTQE.2007.913398 IJSQEN 1077-260X Google Scholar

4.

P. KuchmentL. Kunyansky, “Mathematics of photoacoustic and thermoacousic tomography,” Handbook of Mathematical Methods in Imaging, 817 –865 Springer(2011). Google Scholar

5.

Y. HristovaP. KuchmentL. Nguyen, “Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media,” Inv. Prob., 24 (5), 055006 (2008). http://dx.doi.org/10.1088/0266-5611/24/5/055006 INPEEY 0266-5611 Google Scholar

6.

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

7.

K. P. Köstliet al., “Temporal backward projection of optoacoustic pressure transients using Fourier transform methods,” Phys. Med. Biol., 46 (7), 1863 –1872 (2001). http://dx.doi.org/10.1088/0031-9155/46/7/309 PHMBA7 0031-9155 Google Scholar

8.

M. XuY. XuL. V. Wang, “Time-domain reconstruction algorithms and numerical simulations for thermoacoustic tomography in various geometries,” IEEE Trans. Biomed. Eng., 50 (9), 1086 –1099 (2003). http://dx.doi.org/10.1109/TBME.2003.816081 IEBEAX 0018-9294 Google Scholar

9.

M. XuL. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E, 71 016706 (2005). PLEEE8 1063-651X Google Scholar

10.

B. E. TreebyE. Z. ZhangB. T. Cox, “Photoacoustic tomography in absorbing acoustic media using time reversal,” Inv. Prob., 26 (11), 115003 (2010). http://dx.doi.org/10.1088/0266-5611/26/11/115003 INPEEY 0266-5611 Google Scholar

11.

Z. YuanH. Jiang, “Simultaneous recovery of tissue physiological and acoustic properties and the criteria for wavelength selection in multispectral photoacoustic tomography,” Opt. Lett., 34 (11), 1714 –1716 (2009). http://dx.doi.org/10.1364/OL.34.001714 OPLEDP 0146-9592 Google Scholar

12.

M. Sivaramakrishnanet al., “Limitations of quantitative photoacoustic measurements of blood oxygenation in small vessels,” Phys. Med. Biol., 52 (5), 1349 –1361 (2007). http://dx.doi.org/10.1088/0031-9155/52/5/010 PHMBA7 0031-9155 Google Scholar

13.

J. Xiaoet al., “Quantitative multispectral photoacoustic tomography and wavelength optimization,” J. X-ray Sci. Technol., 18 (4), 415 –427 (2010). JXSTE5 0895-3996 Google Scholar

14.

A. Rosencwaig, Photoacoustics and Photoacoustic Spectroscopy, Wiley(1980). Google Scholar

15.

A. C. TamH. Coufal, “Pulsed opto-acoustics: theory and applications,” J. Physique, C6 9 –21 (1983). 1155-4304 Google Scholar

16.

L. M. Lyamshev, “Optoacoustic sources of sound,” Soviet Phys. Uspekhi, 24 (12), 977 –995 (1981). http://dx.doi.org/10.1070/PU1981v024n12ABEH004757 SOPUAP 0038-5670 Google Scholar

17.

M. W. Sigrist, “Laser generation of acoustic waves in liquids and gases,” J. Appl. Phys., 60 (7), R83 –R121 (1986). http://dx.doi.org/10.1063/1.337089 JAPIAU 0021-8979 Google Scholar

18.

R. A. KrugerP. LiuC. Appledorn, “Photoacoustic ultrasound (PAUS)—reconstruction tomography,” Med. Phys., 22 (10), 1605 –1609 (1995). http://dx.doi.org/10.1118/1.597429 MPHYA6 0094-2405 Google Scholar

19.

P. Burgholzeret al., “Thermoacoustic tomography with integrating area and line detectors,” IEEE Trans. Ultrason. Ferroelectrics Freq. Contr., 52 (9), 1577 –83 (2005). Google Scholar

20.

X. Wanget al., “Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,” Nat. Biotechnol., 21 (7), 803 –806 (2003). http://dx.doi.org/10.1038/nbt839 NABIF9 1087-0156 Google Scholar

21.

C. P. Favazzaet al., “In vivo photoacoustic microscopy of human cutaneous microvasculature and a nevus,” J. Biomed. Opt., 16 (1), 016015 (2011). http://dx.doi.org/10.1117/1.3528661 JBOPFO 1083-3668 Google Scholar

22.

C. Zhanget al., “Reflection-mode submicron-resolution in vivo photoacoustic microscopy,” J. Biomed. Opt., 17 (2), 020501 (2012). http://dx.doi.org/10.1117/1.JBO.17.2.020501 JBOPFO 1083-3668 Google Scholar

23.

B. T. CoxJ. G. LauferP. C. Beard, “The challenges for quantitative photoacoustic imaging,” Proc. SPIE, 7177 717713 (2009). http://dx.doi.org/10.1117/12.806788 PSISDG 0277-786X Google Scholar

24.

A. Ishimaru, Wave Propagation and Scattering in Random Media, Vol 1: Single Scattering and Transport Theory, Academic Press(1978). Google Scholar

25.

Optical-Thermal Response of Laser-Irradiated Tissue, Plenum Press(1995). Google Scholar

26.

F. Martelliet al., Light Propagation through Biological Tissue and Other Diffusive Media, SPIE(2010). Google Scholar

27.

L. V. WangH.-I. Wu, Biomedical Optics: Principles and Imaging, Wiley(2007). Google Scholar

28.

T. Tarvainen, “Computational Methods for Light Transport in Optical Tomography,” University of Kuopio, (2006). Google Scholar

29.

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

30.

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

31.

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

32.

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

33.

J. JosephW. Wiscombe, “The delta-Eddington approximation for radiative flux transfer,” J. Atmos. Sci., 33 2452 –2459 (1976). http://dx.doi.org/10.1175/1520-0469(1976)033<2452:TDEAFR>2.0.CO;2 JAHSAK 0022-4928 Google Scholar

34.

W. Conget al., “Modeling photon propagation in biological tissues using a generalized Delta-Eddington phase function,” Phys. Rev. E, 76 (5), 1 –5 (2007). http://dx.doi.org/10.1103/PhysRevE.76.051913 PLEEE8 1063-651X Google Scholar

35.

T. SpottL. O. Svaasand, “Collimated light sources in the diffusion approximation,” Appl. Opt., 39 (34), 6453 –6465 (2000). http://dx.doi.org/10.1364/AO.39.006453 APOPAI 0003-6935 Google Scholar

36.

L. WangS. JacquesL. Zheng, “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Prog. Biomed., 47 131 –146 (1995). http://dx.doi.org/10.1016/0169-2607(95)01640-F CMPBEK 0169-2607 Google Scholar

37.

S. BanerjeeS. K. Sharma, “Use of Monte Carlo simulations for propagation of light in biomedical tissues,” Appl. Opt., 49 (22), 4152 –4159 (2010). http://dx.doi.org/10.1364/AO.49.004152 APOPAI 0003-6935 Google Scholar

38.

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

39.

M. SchweigerS. R. ArridgeI. Nissilä, “Gauss–Newton method for image reconstruction in diffuse optical tomography,” Phys. Med. Biol., 50 (10), 2365 –86 (2005). http://dx.doi.org/10.1088/0031-9155/50/10/013 PHMBA7 0031-9155 Google Scholar

40.

J. NocedalS. J. Wright, Numerical Optimization, Springer(1999). Google Scholar

41.

R. Fletcher, Practical Methods of Optimization, 2nd ed.Wiley(2000). Google Scholar

42.

B. T. CoxS. R. ArridgeP. C. Beard, “Estimating chromophore distributions from multiwavelength photoacoustic images,” J. Opt. Soc. Am. A, 26 (2), 443 –55 (2009). http://dx.doi.org/10.1364/JOSAA.26.000443 JOAOD6 0740-3232 Google Scholar

43.

G. BalK. Ren, “On multi-spectral quantitative photoacoustic tomography in diffusive regime,” Inv. Prob., 28 (2), 025010 (2012). http://dx.doi.org/10.1088/0266-5611/28/2/025010 INPEEY 0266-5611 Google Scholar

44.

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

45.

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

46.

P. C Hansen, Rank-Deficient and Discrete Ill-Posed Problems, SIAM(1998). Google Scholar

47.

E. Haber, “Quasi-Newton methods for large scale electromagnetic inverse problems,” Inv. Prob., 21 305 –317 (2005). http://dx.doi.org/10.1088/0266-5611/21/1/019 INPEEY 0266-5611 Google Scholar

48.

E. Haber, “A parallel method for large scale time domain electromagnetic inverse problems,” Appl. Numer. Math., 58 (4), 422 –434 (2008). Google Scholar

49.

E. F. CaromeN. A. ClarkC. E. Moeller, “Generation of acoustic signals in liquids by ruby laser-induced thermal stress transients,” Appl. Phys. Lett., 4 (6), 95 –97 (1964). http://dx.doi.org/10.1063/1.1753985 APPLAB 0003-6951 Google Scholar

50.

F. Crosset al., “Time resolved photoacoustic studies of vascular tissue ablation at three laser wavelengths,” Appl. Phys. Lett., 50 (15), 1019 –1021 (1987). http://dx.doi.org/10.1063/1.97994 APPLAB 0003-6951 Google Scholar

51.

F. CrossR. Al-DhahirP. E. Dyer, “Ablative and acoustic response of pulsed UV laser-irradiated vascular tissue in a liquid environment,” J. Appl. Phys., 64 (4), 2194 –2201 (1988). http://dx.doi.org/10.1063/1.341707 JAPIAU 0021-8979 Google Scholar

52.

Z. GuoS. HuL. V. Wang, “Calibration-free absolute quantification of optical absorption coefficients using acoustic spectra in 3D photoacoustic microscopy of biological tissue,” Opt. Lett., 35 (12), 2067 –2069 (2010). http://dx.doi.org/10.1364/OL.35.002067 OPLEDP 0146-9592 Google Scholar

53.

J. A. ViatorS. L. JacquesS. A. Prahl, “Generating subsurface acoustic waves in indocyanine green stained elastin biomaterial using a Q-switched laser,” Proc. SPIE, 3254 104 –111 (1998). http://dx.doi.org/10.1117/12.308154 PSISDG 0277-786X Google Scholar

54.

J. A. ViatorS. L. JacquesS. A. Prahl, “Depth profiling of absorbing soft materials using photoacoustic methods,” IEEE J. Sel. Top. Quant. Electron., 5 (4), 989 –996 (1999). http://dx.doi.org/10.1109/2944.796321 IJSQEN 1077-260X Google Scholar

55.

A. A. KarabutovN. B. PodymovaV. S. Letokhov, “Time-resolved optoacoustic measurement of absorption of light by inhomogeneous media,” Appl. Opt., 34 (9), 1484 –1487 (1995). http://dx.doi.org/10.1364/AO.34.001484 APOPAI 0003-6935 Google Scholar

56.

A. KarabutovN. PodymovaV. Letokhov, “Time-resolved laser optoacoustic tomography of inhomogeneous media,” Appl. Phys. B, 63 545 –563 (1996). APBOEM 0946-2171 Google Scholar

57.

G. PaltaufH. Schmidt-Kloiber, “Pulsed optoacoustic characterization of layered media,” J. Appl. Phys., 88 (3), 1624 –1631 (2000). http://dx.doi.org/10.1063/1.373863 JAPIAU 0021-8979 Google Scholar

58.

M. Jaegeret al., “Diffraction-free acoustic detection for optoacoustic depth profiling of tissue using an optically transparent polyvinylidene fluoride pressure transducer operated in backward and forward mode,” J. Biomed. Opt., 10 (2), 024035 (2005). http://dx.doi.org/10.1117/1.1891443 JBOPFO 1083-3668 Google Scholar

59.

A. A. OraevskyS. L. JacquesF. K. Tittel, “Determination of tissue optical properties by piezoelectric detection of laser-induced stress waves,” Proc. SPIE, 1882 86 –101 (1993). http://dx.doi.org/10.1117/12.147694 PSISDG 0277-786X Google Scholar

60.

A. Oraevskyet al., “Direct measurement of laser fluence distribution and optoacoustic imaging in heterogeneous tissues,” Proc. SPIE, 2323 37 –46 (1995). http://dx.doi.org/10.1117/12.199215 PSISDG 0277-786X Google Scholar

61.

S. L. Jacques, “Simple theory, measurements, and rules of thumb for dosimetry during photodynamic therapy,” Proc. SPIE, 1065 100 –108 (1989). PSISDG 0277-786X Google Scholar

62.

R. Fainchtein, “In-situ determination of concentration and degree of oxygenation of hemoglobin in neural tissue by pulsed photoacoustic spectroscopy,” Proc. SPIE, 2979 417 –428 (1997). http://dx.doi.org/10.1117/12.280274 PSISDG 0277-786X Google Scholar

63.

R. Fainchtein, “Local determination of hemoglobin concentration and degree of oxygenation in tissue by pulsed photoacoustic spectroscopy,” Proc. SPIE, 3916 19 –33 (2000). http://dx.doi.org/10.1117/12.386335 PSISDG 0277-786X Google Scholar

64.

H. F. Zhanget al., “Imaging of hemoglobin oxygen saturation variations in single vessels in vivo using photoacoustic microscopy,” Appl. Phys. Lett., 90 (5), 053901 (2007). http://dx.doi.org/10.1063/1.2435697 APPLAB 0003-6951 Google Scholar

65.

E. V. Savateevaet al., “Optical properties of blood at various levels of oxygenation studied by time resolved detection of laser-induced pressure profiles,” Proc. SPIE, 4618 63 –75 (2002). http://dx.doi.org/10.1117/12.469849 PSISDG 0277-786X Google Scholar

66.

R. O. Esenalievet al., “Optoacoustic technique for noninvasive monitoring of blood oxygenation: a feasibility study,” Appl. Opt., 41 (22), 4722 –4731 (2002). http://dx.doi.org/10.1364/AO.41.004722 APOPAI 0003-6935 Google Scholar

67.

J. Lauferet al., “In vitro measurements of absolute blood oxygen saturation using pulsed near-infrared photoacoustic spectroscopy: accuracy and resolution,” Phys. Med. Biol., 50 (18), 4409 –4428 (2005). http://dx.doi.org/10.1088/0031-9155/50/18/011 PHMBA7 0031-9155 Google Scholar

68.

A. A. OraevskyS. L. JacquesF. K. Tittel, “Measurement of tissue optical properties by time-resolved detection of laser-induced transient stress,” Appl. Opt., 36 (1), 402 –425 (1997). http://dx.doi.org/10.1364/AO.36.000402 APOPAI 0003-6935 Google Scholar

69.

A. OraevskyR. EsenalievA. Karabutov, “Laser optoacoustic tomography of layered tissues: signal processing,” Proc. SPIE, 2979 59 –70 (1997). http://dx.doi.org/10.1117/12.280297 PSISDG 0277-786X Google Scholar

70.

G. J. DieboldT. Sun, “Properties of photoacoustic waves in one, two and three dimensions,” Acustica, 80 339 –351 (1994). ACUSAY 0001-7884 Google Scholar

71.

K. Maslovet al., “Technical considerations in quantitative blood oxygenation measurement using photoacoustic microscopy in vivo,” Proc. SPIE, 6086 60860R (2006). http://dx.doi.org/10.1117/12.646265 PSISDG 0277-786X Google Scholar

72.

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

73.

K. MaslovH. F. ZhangL. V. Wang, “Effects of wavelength-dependent fluence attenuation on the noninvasive photoacoustic imaging of hemoglobin oxygen saturation in subcutaneous vasculature in vivo,” Inv. Prob., 23 (6), S113 –S122 (2007). http://dx.doi.org/10.1088/0266-5611/23/6/S09 INPEEY 0266-5611 Google Scholar

74.

H. F. Zhanget al., “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nat. Biotechnol., 24 (7), 848 –851 (2006). http://dx.doi.org/10.1038/nbt1220 NABIF9 1087-0156 Google Scholar

75.

H. F. ZhangK. MaslovL. V. Wang, “In vivo imaging of subcutaneous structures using functional photoacoustic microscopy,” Nat. Protocol., 2 (4), 797 –804 (2007). http://dx.doi.org/10.1038/nprot.2007.108 NPARDW 1750-2799 Google Scholar

76.

C. Kimet al., “In vivo molecular photoacoustic tomography of melanomas targeted by bioconjugated gold nanocages,” ACS Nano., 4 (8), 4559 –4564 (2010). http://dx.doi.org/10.1021/nn100736c 1936-0851 Google Scholar

77.

J. R. RajianP. L. CarsonX. Wang, “Quantitative photoacoustic measurement of tissue optical absorption spectrum aided by an optical contrast agent,” Opt. Express, 17 (6), 4879 –4889 (2009). http://dx.doi.org/10.1364/OE.17.004879 OPEXFF 1094-4087 Google Scholar

78.

J. Lauferet al., “Quantitative spatially resolved measurement of tissue chromophore concentrations using photoacoustic spectroscopy: application to the measurement of blood oxygenation and haemoglobin concentration,” Phys. Med. Biol., 52 (1), 141 –168 (2007). http://dx.doi.org/10.1088/0031-9155/52/1/010 PHMBA7 0031-9155 Google Scholar

79.

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

80.

J. RipollV. Ntziachristos, “Quantitative point source photoacoustic inversion formulas for scattering and absorbing media,” Phys. Rev. E, 71 (3), 1 –9 (2005). http://dx.doi.org/10.1103/PhysRevE.71.031912 PLEEE8 1063-651X Google Scholar

81.

D. RazanskyV. Ntziachristos, “Hybrid photoacoustic fluorescence molecular tomography using finite-element-based inversion,” Med. Phys., 34 (11), 4293 –4301 (2007). http://dx.doi.org/10.1118/1.2786866 MPHYA6 0094-2405 Google Scholar

82.

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

83.

R. Maet al., “Multispectral optoacoustic tomography (MSOT) scanner for whole-body small animal imaging,” Opt. Express, 17 (24), 21414 –21426 (2009). http://dx.doi.org/10.1364/OE.17.021414 OPEXFF 1094-4087 Google Scholar

84.

V. NtziachristosD. Razansky, “Molecular imaging by means of multispectral optoacoustic tomography (MSOT),” Chem. Rev., 110 (5), 2783 –94 (2010). http://dx.doi.org/10.1021/cr9002566 CHREAY 0009-2665 Google Scholar

85.

W. L. Kiser Jr.et al., “Thermoacoustic in vivo determination of blood oxygenation,” Proc. SPIE, 5320 1 –7 (2004). http://dx.doi.org/10.1117/12.532666 PSISDG 0277-786X Google Scholar

86.

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

87.

B. T. Coxet al., “Quantitative photoacoustic imaging: fitting a model of light transport to the initial pressure distribution,” Proc. SPIE, 5697 49 –55 (2005). http://dx.doi.org/10.1117/12.597190 PSISDG 0277-786X Google Scholar

88.

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

89.

L. YaoY. SunH. Jiang, “Quantitative photoacoustic tomography based on the radiative transfer equation,” Opt. Lett., 34 (12), 1765 –1767 (2009). http://dx.doi.org/10.1364/OL.34.001765 OPLEDP 0146-9592 Google Scholar

90.

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

91.

Z. YuanH. Jiang, “Quantitative photoacoustic tomography: recovery of optical absorption coefficient maps of heterogeneous media,” Appl. Phys. Lett., 88 (23), 231101 (2006). http://dx.doi.org/10.1063/1.2209883 APPLAB 0003-6951 Google Scholar

92.

H. Jianget al., “Quantitative photoacoustic tomography: recovery of optical absorption coefficient maps of heterogeneous media,” Proc. SPIE, 6437 (2), 64371S (2007). http://dx.doi.org/10.1117/12.697079 PSISDG 0277-786X Google Scholar

93.

Y. SunE. SobelH. Jiang, “Quantitative three-dimensional photoacoustic tomography of the finger joints: an in vivo study,” J. Biomed. Opt., 14 (6), 064002 (2009). http://dx.doi.org/10.1117/1.3257246 JBOPFO 1083-3668 Google Scholar

94.

B. Banerjeeet 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 (2008). Google Scholar

95.

Z. YuanQ. WangH. Jiang, “Reconstruction of optical absorption coefficient maps of heterogeneous media by photoacoustic tomography coupled with diffusion equation based regularized Newton method,” Opt. Express, 15 (26), 18076 –18081 (2007). http://dx.doi.org/10.1364/OE.15.018076 OPEXFF 1094-4087 Google Scholar

96.

B. T. CoxS. R. ArridgeP. C. Beard, “Gradient-based quantitative photoacoustic image reconstruction for molecular imaging,” Proc. SPIE, 6437 64371T (2007). http://dx.doi.org/10.1117/12.700006 PSISDG 0277-786X Google Scholar

97.

L. Yinet al., “Tomographic imaging of absolute optical absorption coefficient in turbid media using combined photoacoustic and diffusing light measurements,” Opt. Lett., 32 (17), 2556 –2558 (2007). http://dx.doi.org/10.1364/OL.32.002556 OPLEDP 0146-9592 Google Scholar

98.

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

99.

K. DaoudiW. Steenbergen, “Towards quantitative tissue absorption imaging by combining photoacoustics and acousto-optics,” (2012). Google Scholar

100.

A. RosenthalD. RazanskyV. Ntziachristos, “Quantitative optoacoustic signal extraction using sparse signal representation,” IEEE Trans. Med. Imag., 28 (12), 1997 –2006 (2009). http://dx.doi.org/10.1109/TMI.2009.2027116 ITMID4 0278-0062 Google Scholar

101.

B. CoxJ. LauferP. Beard, “Quantitative Photoacoustic image reconstruction using fluence dependent chromophores,” Biomed. Opt. Express, 1 (1), 5830 –5835 (2010). http://dx.doi.org/10.1364/BOE.1.000201 BOEICL 2156-7085 Google Scholar

102.

R. J. Zempet al., “A new photoacoustic method for measuring optical transport Green’s functions in turbid media,” 1246 –1249 (2008). Google Scholar

103.

J. C. Ranasinghesagaraet al., “Photoacoustic technique for assessing optical scattering properties of turbid media,” J. Biomed. Opt., 14 (4), 040504 (2009). http://dx.doi.org/10.1117/1.3200922 JBOPFO 1083-3668 Google Scholar

104.

G. BalK. Ren, “Multi-source quantitative photoacoustic tomography in a diffusive regime,” Inv. Prob., 27 (7), 075003 (2011). http://dx.doi.org/10.1088/0266-5611/27/7/075003 INPEEY 0266-5611 Google Scholar

105.

H. GaoH. ZhaoS. Osher, “Bregman methods in quantitative photoacoustic tomography,” UCLA CAM Report 10 –42 (2010). Google Scholar

106.

H. GaoH. ZhaoS. Osher, “Quantitative photoacoustic tomography,” UCLA CAM Report 11 –28 (2011). Google Scholar

107.

G. BalA. JollivetV. Jugnon, “Inverse transport theory of photoacoustics,” Inv. Prob., 26 (2), 025011 (2010). http://dx.doi.org/10.1088/0266-5611/26/2/025011 INPEEY 0266-5611 Google Scholar

108.

B. T. CoxT. TarvainenS. Arridge, “Multiple illumination quantitative photoacoustic tomography using transport and diffusion models,” Tomography and Inverse Transport Theory, 559 1 –12 American Mathematial Society,2011). Google Scholar

109.

J. Wanget al., “Saturation effect in functional photoacoustic imaging,” J. Biomed. Opt., 15 (2), 021317 (2010). http://dx.doi.org/10.1117/1.3333549 JBOPFO 1083-3668 Google Scholar

110.

G. Paltaufet al., “Spectral optoacoustic imaging using a scanning transducer,” Proc. SPIE, 4434 81 –88 (2001). http://dx.doi.org/10.1117/12.446661 PSISDG 0277-786X Google Scholar

111.

A. KarabutovE. SavateevaA. Oraevsky, “Imaging of layered structures in biological tissues with opto-acoustic front surface transducer,” Proc. SPIE, 3601 284 –296 (1999). PSISDG 0277-786X Google Scholar

112.

A. A. Karabutovet al., “Backward mode detection of laser-induced wide-band ultrasonic transients with optoacoustic transducer,” J. Appl. Phys., 87 (4), 2003 –2014 (2000). http://dx.doi.org/10.1063/1.372127 JAPIAU 0021-8979 Google Scholar

113.

A. Danielliet al., “Picosecond absorption relaxation measured with nanosecond laser photoacoustics,” Appl. Phys. Lett., 97 (16), 163701 (2010). http://dx.doi.org/10.1063/1.3500820 APPLAB 0003-6951 Google Scholar

114.

F. Banfiet al., “Temperature dependence of the thermal boundary resistivity of glass-embedded metal nanoparticles,” Appl. Phys. Lett., 100 (1), 011902 (2012). http://dx.doi.org/10.1063/1.3673559 APPLAB 0003-6951 Google Scholar

115.

Y.-S. Chenet al., “Silica-coated gold nanorods as photoacoustic signal nanoamplifiers,” Nano Lett., 11 (2), 348 –354 (2011). http://dx.doi.org/10.1021/nl1042006 NALEFD 1530-6984 Google Scholar

116.

B. SoroushianW. M. WhelanM. C. Kolios, “Study of laser-induced thermoelastic deformation of native and coagulated ex-vivo bovine liver tissues for estimating their optical and thermomechanical properties,” J. Biomed. Opt., 15 (6), 065002 (2010). http://dx.doi.org/10.1117/1.3517455 JBOPFO 1083-3668 Google Scholar

117.

J. LauferE. ZhangP. Beard, “Evaluation of absorbing chromophores used in tissue phantoms for quantitative photoacoustic spectroscopy and imaging,” IEEE J. Sel. Top. Quant. Electron., 16 (3), 600 –607 (2010). http://dx.doi.org/10.1109/JSTQE.2009.2032513 IJSQEN 1077-260X Google Scholar

Citation Download Citation

Benjamin T. Cox, Jan G. Laufer, Paul C. Beard, and Simon R. Arridge "Quantitative spectroscopic photoacoustic imaging: a review," Journal of Biomedical Optics 17(6), 061202 (1 June 2012). https://doi.org/10.1117/1.JBO.17.6.061202

Published: 1 June 2012

Access the abstract

JOURNAL ARTICLE
23 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 583 scholarly publications and 6 patents.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Absorption

Scattering

Chromophores

Acoustics

Photoacoustic spectroscopy

Tissue optics

Photoacoustic imaging

1.

Introduction

1.1.

Scope of the Paper

1.1.1.

What does “quantitative photoacoustic imaging” cover?

1.1.2.

What is not covered in this paper?

1.1.3.

Layout of the paper

1.2.

Photoacoustic Imaging

1.2.1.

The photoacoustic effect

Eq. (1)

Eq. (2)

Eq. (3)

Fig. 1

1.2.2.

Photoacoustic tomography

1.2.3.

Photoacoustic scanning microscopy

2.

Quantitative Photoacoustic Imaging

2.1.

Optical Inverse Problem

2.1.1.

Formulation of the problem

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Fig. 2

2.1.2.

Spectral coloring and structural distortion

2.1.3.

Using a model to recover the concentrations

2.2.

Fluence Models

2.2.1.

Radiative transfer equation

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

2.2.2.

Diffusion approximation

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

Eq. (16)

Eq. (17)

2.2.3.

Collimated source, δ-Eddington approximation

2.2.4.

Monte Carlo models

2.3.

Parameter Estimation

2.3.1.

Inversion techniques

Linearization

Direct inversion

Fixed-point iteration

Model-based minimization

2.3.2.

Least-squares minimization

Eq. (18)

2.3.3.

Uniqueness and ill-posedness

2.3.4.

Large-scale inversions and domain parameterization

Eq. (19)

3.

One-Dimensional Quantitative Photoacoustic Imaging

3.1.

Homogeneous Nonscattering Medium: Beer’s Law

Eq. (20)

Eq. (21)

Eq. (22)

Collimated source, $δ$ -Eddington approximation