## 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.^{1}2.^{–}^{3} 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. 45.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.^{11}12.^{–}^{13} 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,^{14} 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.^{15}16.^{–}^{17} 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 $\text{vibrational}/\text{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

where $\hat{\mathrm{\Gamma}}$ is the PA efficiency. (This efficiency can be identified with the Grüneisen parameter $\mathrm{\Gamma}$ 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 where ${\mu}_{a}$ is the optical absorption coefficient (the total absorption due to the chromophores) and $\mathrm{\Phi}$ 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:

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.

## 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.^{18} 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.^{4}5.6.7.8.9.^{–}^{10} 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,^{19} and of detectors focused in a plane,^{20} 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\text{\hspace{0.17em}}$ *μ*m to a depth of a few millimeters, with the vertical resolution limited to half the shortest wavelength detectable ($\sim 15$ *μ*m).^{21} 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.^{22}

## 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}(\mathbf{x},\lambda )$ are obtained experimentally, each at a different optical wavelength $\lambda $. (The position vector $\mathbf{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(\mathbf{x},\lambda )={\hat{\mathrm{\Gamma}}}^{-1}{p}_{0}(\mathbf{x},\lambda )$ 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}(\mathbf{x})$ of $K$ chromophores with known molar absorption coefficient spectra ${\alpha}_{k}(\lambda )$, given the absorbed energy images $H(\mathbf{x},\lambda )$ when the chromophores are linked to the absorption coefficients with the linear mapping ${L}_{\lambda}$:

## (4)

$${\mu}_{a}(\mathbf{x},\lambda )={L}_{\lambda}({C}_{k})=\sum _{k=1}^{K}{C}_{k}(\mathbf{x}){\alpha}_{k}(\lambda ),$$## (5)

$$H(\mathbf{x},\lambda )=T({\mu}_{a},{\mu}_{s})={\mu}_{a}(\mathbf{x},\lambda )\mathrm{\Phi}(\mathbf{x},\lambda ;{\mu}_{a},{\mu}_{s}).$$The light fluence, $\mathrm{\Phi}(\mathbf{x},\lambda )$ is unknown and will depend on ${\mu}_{a}$ as well as the optical scattering coefficient ${\mu}_{s}$. A fluence model is therefore required (Sec. 2.2). Equations (4) and (5) can be combined to give

## (6)

$$H(\mathbf{x},\lambda )={T}_{\lambda}({C}_{k},{\mu}_{s})=\mathrm{\Phi}(\mathbf{x},\lambda ;{C}_{k},{\mu}_{s})\sum _{k=1}^{K}{C}_{k}(\mathbf{x}){\alpha}_{k}(\lambda ).$$Equations (4) and (5) suggest a two-stage inversion strategy: first recover the absorption coefficients ${\mu}_{a}={T}^{-1}(H)$ and then find the concentrations ${C}_{k}={L}_{\lambda}^{-1}({\mu}_{a})$, whereas Eq. (6) suggests a single inversion, ${C}_{k}={T}_{\lambda}^{-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

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 $\text{transport}/\text{absorption}$ model must be taken into account.For OR-PAM images obtained by scanning the focused illumination spot across the region of interest, the fluence distribution $\mathrm{\Phi}(\mathbf{x})$ will be different for each measurement, so Eq. (5) should be replaced by $H({\mathbf{x}}_{i})={\mu}_{a}({\mathbf{x}}_{i}){\mathrm{\Phi}}_{i}[{\mathbf{x}}_{i},{\mu}_{a}(\mathbf{x})]$ where ${\mathrm{\Phi}}_{i}(\mathbf{x})$ is the fluence distribution at point $\mathbf{x}$ when the illumination is focused at the point ${\mathbf{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(\lambda )$ straightforwardly, for example, using a linear matrix inversion. Unfortunately the unknown fluence $\mathrm{\Phi}$ depends in a nontrivial way on the chromophore concentrations because ${T}_{\lambda}$ is both nonlinear and nonlocal ($H$ at one point can be affected by ${\mu}_{a}$ some distance away via its effect on the fluence). It is not reasonable, therefore, to assume that $H(\lambda )\propto {\mu}_{a}(\lambda )$ 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(\lambda )$ 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 $\mathrm{\Phi}(\mathbf{x},\lambda )$ 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(\mathbf{x})={\mu}_{a}(\mathbf{x})\mathrm{\Phi}(\mathbf{x})$. The unknown and nonuniform light distribution $\mathrm{\Phi}(\mathbf{x})$ will therefore distort what would, in its absence, have been an image proportional to ${\mu}_{a}(\mathbf{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.^{23}

## 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}_{\lambda}$, $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}_{\lambda}$ or $T$ to obtain ${\mu}_{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.^{24}25.26.^{–}^{27} 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:

## (8)

$$\frac{1}{c}\frac{\partial \varphi}{\partial t}(\mathbf{x},\hat{\mathbf{s}},t)=q(\mathbf{x},\hat{\mathbf{s}},t)-[\hat{\mathbf{s}}\xb7\nabla +{\mu}_{a}(\mathbf{x})+{\mu}_{s}(\mathbf{x})]\varphi (\mathbf{x},\hat{\mathbf{s}},t)\phantom{\rule{0ex}{0ex}}+{\mu}_{s}\int \mathrm{\Theta}(\hat{\mathbf{s}},{\hat{\mathbf{s}}}^{\prime})\varphi (\mathbf{x},{\hat{\mathbf{s}}}^{\prime},t)\mathrm{d}{\hat{\mathbf{s}}}^{\prime},$$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.

## (9)

$$(\hat{\mathbf{s}}\xb7\nabla +{\mu}_{t})\varphi (\mathbf{x},\hat{\mathbf{s}})-{\mu}_{s}\int \mathrm{\Theta}(\hat{\mathbf{s}},{\hat{\mathbf{s}}}^{\prime})\varphi (\mathbf{x},{\hat{\mathbf{s}}}^{\prime})\mathrm{d}{\hat{\mathbf{s}}}^{\prime}=q(\mathbf{x},\hat{\mathbf{s}}),$$## (10)

$$\mathrm{\Phi}(\mathbf{x})=\int \varphi (\mathbf{x},{\hat{\mathbf{s}}}^{\prime})\mathrm{d}{\hat{\mathbf{s}}}^{\prime}.$$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.^{28}^{,}^{29}

When there is no scattering, Eq. (9) reduces to $(\hat{\mathbf{s}}\xb7\nabla +{\mu}_{a})\varphi =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 \mathrm{\Phi}/\partial z=-{\mu}_{a}\mathrm{\Phi}$. If the source illuminates the surface with a fluence of ${\mathrm{\Phi}}_{0}$ and ${\mu}_{a}$ is constant, then the fluence as a function of depth may be written

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^{30}

## (12)

$${\mu}_{a}\mathrm{\Phi}+\nabla \xb7\mathbf{F}={q}_{0},\phantom{\rule[-0.0ex]{2em}{0.0ex}}\frac{\mathbf{F}}{D}+\nabla \mathrm{\Phi}=3\text{\hspace{0.17em}}{\mathbf{q}}_{1},$$## (13)

$$\mathbf{F}(\mathbf{x})=\int \varphi (\mathbf{x},{\hat{\mathbf{s}}}^{\prime}){\hat{\mathbf{s}}}^{\prime}\mathrm{d}{\hat{\mathbf{s}}}^{\prime}.$$The isotropic and mildly directional source terms ${q}_{o}$ and ${\mathbf{q}}_{1}$ are defined analogously to Eqs. (10) and (13). When ${\mathbf{q}}_{1}=0$, Eqs. (12) and (13) reduce to the diffusion approximation (DA):

The condition ${\mu}_{s}^{\prime}\gg {\mu}_{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.^{31}^{,}^{32} For a homogeneous medium, the solution to Eq. (14) is given by

## (15)

$$\mathrm{\Phi}(\mathbf{x})=\int {G}_{0}(\mathbf{x}-{\mathbf{x}}^{\prime}){q}_{0}({\mathbf{x}}^{\prime})\mathrm{d}{\mathbf{x}}^{\prime}$$## (16)

$${G}_{0}(\mathbf{x},{\mathbf{x}}^{\prime})=\frac{\mathrm{exp}(-{\mu}_{\mathrm{eff}}|\mathbf{x}-{\mathbf{x}}^{\prime}|)}{4\pi D|\mathbf{x}-{\mathbf{x}}^{\prime}|},\phantom{\rule{0ex}{0ex}}{\mu}_{\mathrm{eff}}=\sqrt{3{\mu}_{a}({\mu}_{a}+{\mu}_{s}^{\prime})},$$*effective attenuation coefficient*. This solution will be used in some of the inversions in Sec. 5. In one dimension (1-D), the solution is $\mathrm{\Phi}(z)={\mathrm{\Phi}}_{0}\text{\hspace{0.17em}}\mathrm{exp}(-{\mu}_{\mathrm{eff}}z)$ where ${\mathrm{\Phi}}_{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

## (17)

$$\mathrm{\Phi}(z)\approx k{\mathrm{\Phi}}_{0}\text{\hspace{0.17em}}\mathrm{exp}(-{\mu}_{\mathrm{eff}}\text{\hspace{0.17em}}z)\phantom{\rule[-0.0ex]{1em}{0.0ex}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}\text{\hspace{0.17em}}z\gg \frac{1}{{\mu}_{t}},$$## 2.2.3.

#### Collimated source, $\delta $-Eddington approximation

For a collimated incident beam, the condition for the accuracy of the DA, ${\mu}_{s}^{\prime}\gg {\mu}_{a}$, must be supplemented with a second condition. As the collimated part of an incident beam will decay exponentially with depth $z$ as $\mathrm{exp}(-{\mu}_{t}z)$, the fluence will only be dominated by diffuse photons at depths $z\gg 1/{\mu}_{t}$. The optical coefficients in PA imaging applications are of the order of ${\mu}_{a}=0.1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, ${\mu}_{s}=10\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, $g=0.9$, ${\mu}_{s}^{\prime}=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ so the condition that ${\mu}_{s}^{\prime}\gg {\mu}_{a}$ is met, but the condition that the fluence is dominated by diffuse photons, only holds for distances $z\gg 1/{\mu}_{t}\approx 0.1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, so perhaps for $z>1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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 $\delta $-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.^{33}34.^{–}^{35}

## 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 $\mathrm{\Theta}$.^{36}^{,}^{37} 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.^{38} 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: ${\mu}_{a}=-\mathrm{ln}[\mathrm{\Phi}(z)/{\mathrm{\Phi}}_{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 ${\mu}_{a}$ on one side and $H/[\mathrm{\Phi}({\mu}_{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}^{\mathrm{obs}}$, then the squared error between it and the output of a forward model, ${H}^{\text{model}}({C}_{k})$, is ${\Vert {H}^{\text{model}}({C}_{k})-{H}^{\mathrm{obs}}\Vert}^{2}=\sum _{i}{[{H}_{i}^{\text{model}}({C}_{k})\phantom{\rule{0ex}{0ex}}-{H}_{i}^{\mathrm{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, $\epsilon $, by adjusting the unknowns ${C}_{k}$:

## (18)

$$\underset{{C}_{k}}{\mathrm{argmin}}\text{\hspace{0.17em}}\epsilon =\frac{1}{2}{\Vert {H}^{\text{model}}({C}_{k})-{H}^{\mathrm{obs}}\Vert}^{2}+\mathcal{P}({C}_{k}).$$Solving Eq. (18) means finding the minimum of the scalar functional $\epsilon $. Gradient-based methods use the functional gradient (the vector of first derivatives), ${g}_{i}=\partial \epsilon /\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, ${\mathcal{H}}_{ij}={\partial}^{2}\epsilon /(\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}_{ij}=\partial {H}_{i}^{\text{model}}/\partial {C}_{k,j}$ to estimate the Hessian matrix as $\mathcal{H}\approx {J}^{T}J$, which can be implemented with Krylov methods that require only matrix-vector products to be computed.^{39} 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.^{40}^{,}^{41}

## 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 $\epsilon $ 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;^{42}^{,}^{43} or by using multiple measurements made with different surface illumination patterns.^{44}^{,}^{45}

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}_{\lambda}$ 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}_{\lambda}^{-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.^{46}

## 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^{47}^{,}^{48} 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

## (19)

$${\mu}_{a}(\mathbf{x})=\sum _{n}{\mu}_{a,n}\text{\hspace{0.17em}}{S}_{n}(\mathbf{x}),\phantom{\rule[-0.0ex]{2em}{0.0ex}}{S}_{n}(\mathbf{x})=\{\begin{array}{cc}1,& \mathbf{x}\in {A}_{n}\\ 0,& \text{otherwise}\end{array}.$$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

## (20)

$$H(z)=\{\begin{array}{rr}{\mu}_{a}{\mathrm{\Phi}}_{0}\text{\hspace{0.17em}}\mathrm{exp}(-{\mu}_{a}z),& z\ge 0\\ 0,& z<0\end{array},$$## (21)

$$p(t)=\left(\frac{\hat{\mathrm{\Gamma}}{\mathrm{\Phi}}_{0}{\mu}_{a}}{2}\right)\mathrm{exp}(-{\mu}_{a}{c}_{0}t),\phantom{\rule[-0.0ex]{2em}{0.0ex}}t\ge 0,$$^{49}50.

^{–}

^{51}

More recently, a frequency domain method for quantification of ${\mu}_{a}$ has been proposed by Guo et al.^{52} 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(\omega )|=(\hat{\mathrm{\Gamma}}{\mathrm{\Phi}}_{0}{\mu}_{a}/2){[{({\mu}_{a}{c}_{0})}^{2}+{\omega}^{2}]}^{-\frac{1}{2}}$, the ratio of two measurements made at different optical wavelengths, ${\lambda}_{1}$ and ${\lambda}_{2}$, may be written as

## (22)

$$\frac{|P({\lambda}_{1};\omega )|}{|P({\lambda}_{2};\omega )|}=\frac{{\mathrm{\Phi}}_{0}({\lambda}_{1})}{{\mathrm{\Phi}}_{0}({\lambda}_{2})}\{\frac{{c}_{0}^{2}+{[\omega /{\mu}_{a}({\lambda}_{2})]}^{2}}{{c}_{0}^{2}+{[\omega /{\mu}_{a}({\lambda}_{1})]}^{2}}{\}}^{1/2},$$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\text{\hspace{0.17em}}$ *μ*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 $\mathrm{exp}(-{\mu}_{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$) was recovered to within about $\pm 4\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{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 ${\mathrm{sO}}_{2}$) in a superficial vein and artery pair selected from $1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{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 $\mathrm{\Delta}z$ will absorb energy per unit area $\mathrm{\Phi}[1-\mathrm{exp}(-{\mu}_{a}\mathrm{\Delta}z)]$, where $\mathrm{\Phi}$ 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.^{53}^{,}^{54} When the layers become so thin that the absorption coefficient becomes a continuous function of depth, the absorbed energy density may be written as

## (23)

$$H(z)={\mu}_{a}(z){\mathrm{\Phi}}_{0}\text{\hspace{0.17em}}\mathrm{exp}[-{\int}_{0}^{z}{\mu}_{a}(\zeta )\mathrm{d}\zeta ].$$^{55}

^{,}

^{56}inverted this expression to find the following formula for the absorption coefficient given the measured pressure $p(t)=\hat{\mathrm{\Gamma}}H(t={z/c}_{0})$:

## (24)

$${\mu}_{a}(z)=\frac{p(z/{c}_{0})}{{c}_{0}{\int}_{z/{c}_{0}}^{\infty}p(t)\mathrm{d}t},\phantom{\rule[-0.0ex]{2em}{0.0ex}}z\ge 0.$$This has been experimentally verified by recovering absorption profiles for depth-dependent solutions of magnetite particles in oil^{55}^{,}^{56} 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.^{57}^{,}^{58}

## 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 ${\mu}_{a}$. For a collimated beam, the fluence close to the surface will decay approximately as $\mathrm{exp}(-{\mu}_{t}z)$ and at deeper depths where the light is diffuse as $\mathrm{exp}(-{\mu}_{\mathrm{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 ${\mu}_{a}$ either. Oraevsky et al.^{59}^{,}^{60} used the approximate model for the fluence^{61} given by Eq. (17), according to which the absorbed energy profile is $H(z)\approx {\mu}_{a}k{\mathrm{\Phi}}_{0}\text{\hspace{0.17em}}\mathrm{exp}(-{\mu}_{\mathrm{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.^{61} By measuring this diffuse reflectance and fitting Eq. (17) to the exponentially decaying part of the measurement, ${\mu}_{a}$ could be inferred from the amplitude of the fitted curve extrapolated to $z=0$, that is, from ${\mu}_{a}k{\mathrm{\Phi}}_{0}$. Subsequently ${\mu}_{\mathrm{eff}}$ was measured from the slope of the curve, and ${\mu}_{s}^{\prime}$ calculated from it, by use of the known ${\mu}_{a}$. Rather than measuring the diffuse reflectance, Fainchtein et al.^{62}^{,}^{63} modeled it as approximately ${R}_{d\infty}\approx \mathrm{exp}(-7{\mu}_{a}/{\mu}_{\mathrm{eff}})$. Note that in both cases an additional optical measurement, or an additional assumption, is required to allow both ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ 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, ${\mathrm{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 ${\mathrm{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, ${\mathrm{sO}}_{2}$, is defined as the ratio

## (25)

$${\mathrm{sO}}_{2}=\frac{{C}_{{\mathrm{HbO}}_{2}}}{{C}_{{\mathrm{HbO}}_{2}}+{C}_{\mathrm{HHb}}},$$## (26)

$$\left(\begin{array}{c}{\mu}_{a,\text{blood}}({\lambda}_{1})\\ \vdots \\ {\mu}_{a,\text{blood}}({\lambda}_{N})\end{array}\right)=\left(\begin{array}{cc}{\alpha}_{{\mathrm{HbO}}_{2}}({\lambda}_{1})& {\alpha}_{\mathrm{HHb}}({\lambda}_{1})\\ \vdots & \vdots \\ {\alpha}_{{\mathrm{HbO}}_{2}}({\lambda}_{N})& {\alpha}_{\mathrm{HHb}}({\lambda}_{N})\end{array}\right)\left(\begin{array}{c}{C}_{{\mathrm{HbO}}_{2}}\\ {C}_{\mathrm{HHb}}\end{array}\right).$$^{12}To use the amplitudes of PA signals as relative measurements of ${\mu}_{a}$, it is necessary to find a scenario in which the constant $K$, which for PA will be $K=\hat{\mathrm{\Gamma}}\mathrm{\Phi}$, is independent of wavelength. This is not true in general, in fact it is rarely true, and assuming it to be so

^{64}is likely to result in significant errors in ${\mathrm{sO}}_{2}$ estimates. If PA measurements of ${\mathrm{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}_{{\mathrm{HbO}}_{2}}$ and ${C}_{\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{MHz}$), Fainchtein et al.^{62}^{,}^{63} 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.^{65} used a broadband detector ($\sim 40\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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.^{66} showed a similar result at 1064 nm with a lower bandwidth transducer. However, none of these authors provided estimates of absolute ${\mathrm{sO}}_{2}$.

Building on these preliminary results, Laufer et al.^{67} estimated blood ${\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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 $\mathrm{exp}(-{\mu}_{a}z)$ and instead was modeled by a 1-D version of the $\delta $-Eddington approximation. This model was fitted by a Nelder-Mead algorithm to both the measured PA amplitudes and, separately, to measurements of ${\mu}_{\mathrm{eff}}$ obtained from the slopes of the exponentially decaying PA signals. The ${\mathrm{sO}}_{2}$ was obtained with an accuracy of $\pm 2.5\%$ when using measurements of ${\mu}_{\mathrm{eff}}$ and $\pm 4\%$ from amplitude measurements. The smallest change in ${\mathrm{sO}}_{2}$ that could be accurately measured was $\pm 1\%$. These promising results suggest that PA measurements can be used to estimate ${\mathrm{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\ll 1/{\mu}_{t}$, where ballistic (unscattered) photons dominate over scattered photons, the fluence will decay one-dimensionally as $\mathrm{\Phi}(z)=\phantom{\rule{0ex}{0ex}}{\mathrm{\Phi}}_{0}\text{\hspace{0.17em}}\mathrm{exp}(-{\mu}_{t}z)$ and so the acoustic signal can be written as $p(t)=(\hat{\mathrm{\Gamma}}{\mathrm{\Phi}}_{0}{\mu}_{a}/2)\mathrm{exp}(-{\mu}_{t}{c}_{0}t)$ for $t\ll 1/({c}_{0}{\mu}_{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 $\mathrm{exp}(-{\mu}_{\mathrm{eff}}z)$ in the diffusive regime, $z\gg 1/{\mu}_{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$.^{25} With smaller illumination regions, the light fluence will decay faster than $\mathrm{exp}(-{\mu}_{\mathrm{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.”^{57}^{,}^{58}^{,}^{68}^{,}^{69} 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.^{70}

## 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 ${\mu}_{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.^{12}^{,}^{71} 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}/{\lambda}_{c}$, where ${\lambda}_{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 ${\mu}_{a}$: (1) for small vessels $a\ll 1/{\mu}_{a}$, and (2) for large vessels $a>{\lambda}_{c}$ when a transducer with ${\lambda}_{c}<1/{\mu}_{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 ${\mathrm{sO}}_{2}$. By making a measurement at an isosbestic point (a wavelength at which ${\alpha}_{{\mathrm{HbO}}_{2}}={\alpha}_{\mathrm{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 ${\mu}_{a}\approx 180\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$(${\mu}_{a}{\lambda}_{c}\approx 1$, ${\lambda}_{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 (${\mu}_{a}{\lambda}_{c}\approx 3$, ${\lambda}_{c}<a$) failed to recover the spectrum correctly.^{12}

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.^{72} 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, ${\mathrm{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.^{71}^{,}^{73}74.^{–}^{75} 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 ${\mathrm{sO}}_{2}$ in the skin of small animals using PA microscopy.^{73}^{,}^{74}^{,}^{76}

## 4.1.3.

#### Correction for wavelength dependence using a known contrast agent

In a similar spirit but less invasively, Rajian et al.^{77} 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,

## (28)

$${p}_{2}(\lambda )={B}^{\prime}{\mathrm{\Phi}}^{\prime}(\lambda )[\sum _{k}{C}_{k}^{\prime}{\alpha}_{k}(\lambda )+{C}_{\mathrm{CA}}{\alpha}_{\mathrm{CA}}(\lambda )],$$## (29)

$$\mathrm{\Phi}(\lambda )\approx \frac{{p}_{1}(\lambda )-{p}_{2}(\lambda )}{{\mathrm{BC}}_{\mathrm{CA}}{\alpha}_{\mathrm{CA}}(\lambda )}.$$## (30)

$$\frac{{C}_{\mathrm{CA}}{\alpha}_{\mathrm{CA}}(\lambda ){p}_{1}(\lambda )}{{p}_{1}(\lambda )-{p}_{2}(\lambda )}\approx \sum _{k}{C}_{k}{\alpha}_{k}(\lambda ),$$## 4.2.

### Model-Based Minimization

Laufer et al.^{78} 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^{67}^{,}^{78}^{,}^{79} 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}_{\lambda}$ and $A$). The light transport was modeled with a two-dimensional (2-D) finite element $\delta $-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:

## (31)

$$\text{tubes}\text{\hspace{0.17em}}{\mu}_{a,i}(\lambda )={C}_{\mathrm{HHb},i}\text{\hspace{0.17em}}{\alpha}_{\mathrm{HHb}}(\lambda )+{C}_{{\mathrm{HbO}}_{2},i}{\alpha}_{{\mathrm{HbO}}_{2}}(\lambda )\phantom{\rule{0ex}{0ex}}+{\mu}_{a,{\mathrm{H}}_{2}\mathrm{O}}(\lambda ),\phantom{\rule{0ex}{0ex}}\text{background}\text{\hspace{0.17em}}{\mu}_{a,\mathrm{bg}}(\lambda )={C}_{\mathrm{HHb},\mathrm{bg}}{\alpha}_{\mathrm{HHb}}(\lambda )+{C}_{{\mathrm{HbO}}_{2,\mathrm{bg}}}{\alpha}_{{\mathrm{HbO}}_{2}(\lambda )}\phantom{\rule{0ex}{0ex}}+{\mu}_{a,{\mathrm{H}}_{2}\mathrm{O}}(\lambda )+{\mu}_{a,\text{lipids}}(\lambda )\phantom{\rule{0ex}{0ex}}+{C}_{\text{dye}}{\alpha}_{\text{dye}}(\lambda ).$$The scattering coefficient was written as ${\mu}_{s}(\lambda )={k}_{s}{\alpha}_{\mathrm{scat}}(\lambda )$, where the wavelength dependence of the scattering ${\alpha}_{\mathrm{scat}}(\lambda )$ was considered known in advance, although the amplitude ${k}_{s}$ was not. There were therefore 11 unknown scalars to be determined: ${C}_{{\mathrm{HbO}}_{2}}$ and ${C}_{\mathrm{HHb}}$ for each tube and the background; ${C}_{\text{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, ${\mathrm{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,^{79} 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(\mathbf{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}_{\lambda}^{-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}_{\lambda}^{-1}$. Section 5.2.4 solves the combined spectral and optical inversion, ${T}_{\lambda}^{-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)]

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 ${\mu}_{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, ${\mu}_{a}(\mathbf{x})={\mu}_{a,0}+\delta {\mu}_{a}(\mathbf{x})$. The resulting fluence can be written as a sum of an unperturbed part, ${\mathrm{\Phi}}_{0}(\mathbf{x})$, corresponding to ${\mu}_{a,0}$, and a second part $\delta \mathrm{\Phi}(\mathbf{x})$, corresponding to the absorption perturbation $\delta {\mu}_{a}(\mathbf{x})$. The change in the absorbed energy distribution will (neglecting the second-order term) consist of two terms:## (33)

$$\delta H(\mathbf{x})\approx \delta {\mu}_{a}(\mathbf{x}){\mathrm{\Phi}}_{0}(\mathbf{x})+{\mu}_{a,0}\delta \mathrm{\Phi}(\mathbf{x}),$$## 5.1.1.

#### Unchanged fluence

In Eq. (33), the change in the absorbed energy, $\delta H$, is still nonlinearly related to the change in the absorption coefficient, $\delta {\mu}_{a}$, because $\delta \mathrm{\Phi}$ depends on $\delta {\mu}_{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:

This immediately shows why difference imaging will never remove the effect of the fluence even for very small absorption perturbations: if $\delta H$ is considered to be the difference between two images, $\delta H={H}_{2}-{H}_{1}$, then—while it depends on the difference in absorption coefficient, $\delta {\mu}_{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:

When the fluence due to the homogeneous part of the absorption, ${\mathrm{\Phi}}_{0}$, is known, Eq. (34) can be inverted trivially to obtain an estimate for the absorption perturbation $\delta {\mu}_{a}(\mathbf{x})=\phantom{\rule{0ex}{0ex}}\delta H(\mathbf{x})/{\mathrm{\Phi}}_{0}(\mathbf{x})$. Ripoll and Ntziachristos^{80} 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 $\mathrm{\Phi}$ calculated assuming homogeneous (and known) background optical properties.^{81}82.^{–}^{83} 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.”^{84} Earlier, Kiser et al.^{85} 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 ${\mathrm{\Phi}}_{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.^{86} 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 $\delta \mathrm{\Phi}$ due to the absorption perturbation $\delta {\mu}_{a}$ also obeys a diffusion equation:

## (35)

$$({\mu}_{a,0}-{D}_{0}{\nabla}^{2})\delta \mathrm{\Phi}(\mathbf{x})=-\delta {\mu}_{a}(\mathbf{x})\mathrm{\Phi}(\mathbf{x}),$$## (36)

$$\delta \mathrm{\Phi}(\mathbf{x})=-{\int}_{{\mathbb{R}}^{n}}{G}_{0}(\mathbf{x},{\mathbf{x}}^{\prime})\delta {\mu}_{a}({\mathbf{x}}^{\prime})\mathrm{\Phi}({\mathbf{x}}^{\prime})\mathrm{d}{\mathbf{x}}^{\prime}.$$## (37)

$$\delta H(\mathbf{x})\approx \delta {\mu}_{a}(\mathbf{x}){\mathrm{\Phi}}_{0}(\mathbf{x})\phantom{\rule{0ex}{0ex}}-{\mu}_{a,0}{\int}_{{\mathbb{R}}^{n}}{G}_{0}(\mathbf{x},{\mathbf{x}}^{\prime})\delta {\mu}_{a}({\mathbf{x}}^{\prime}){\mathrm{\Phi}}_{0}({\mathbf{x}}^{\prime})d{\mathbf{x}}^{\prime}.$$As it is linear, this equation may be discretized and written as a matrix equation relating vectors of $\delta H$ and $\delta {\mu}_{a}$, which can then be inverted with any suitable method.

Interestingly, Zemp^{86} 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.^{87}^{,}^{88} rearranged Eq. (32) to obtain a fixed-point iteration for ${\mu}_{a}(\mathbf{x})$ that makes no linearizing assumptions:

## (38)

$${\mu}_{a}^{(n+1)}(\mathbf{x})=\frac{H(\mathbf{x})}{{\mathrm{\Phi}}^{(n)}[\mathbf{x},{\mu}_{a}^{(n)}(\mathbf{x})]},$$^{87}

^{,}

^{88}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

^{89}(Sec. 5.3). Jetzfellner et al.

^{90}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.

^{91}92.

^{–}

^{93}

## 5.1.4.

#### Noniterative inversion

Banerjee et al.^{94} realized that the absorbed energy ${\mu}_{a}\mathrm{\Phi}$ appears directly in the DA model, Eq. (14), and the following equation can therefore be written for the fluence:

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

## (40)

$$\underset{{\mu}_{a}}{\mathrm{argmin}}\text{\hspace{0.17em}}\epsilon =\frac{1}{2}{\Vert {H}^{\text{model}}({\mu}_{a})-{H}^{\mathrm{obs}}\Vert}^{2}+\mathcal{P}({\mu}_{a}),$$^{95}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

^{96}(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={\mu}_{a}\mathrm{\Phi}$, if the fluence, $\Phi $, were known, then the inversion for ${\mu}_{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.^{97}^{,}^{98} 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.^{99}]

## 5.1.7.

#### Sparse basis decomposition

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

## (41)

$$\mathrm{log}[H(\mathbf{x})]=\sum _{n}{a}_{n}{\varphi}_{n}(\mathbf{x})+\sum _{m}{b}_{m}{\psi}_{m}(\mathbf{x}),$$## 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.^{101}

## 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,^{102}^{,}^{103} although they either only recover the scattering close to the tissue surface or give an averaged bulk value, so techniques to invert $T({\mu}_{a},{\mu}_{s})$ are needed in general. Under the DA the scattering is described by the reduced scattering coefficient ${\mu}_{s}^{\prime}$, or equivalently the optical diffusion coefficient $D$, and Eq. (32) becomes one or the other of

## (42)

$$H(\mathbf{x})={\mu}_{a}(\mathbf{x})\mathrm{\Phi}[\mathbf{x},{\mu}_{a}(\mathbf{x}),{\mu}_{s}^{\prime}(\mathbf{x})],\phantom{\rule[-0.0ex]{2em}{0.0ex}}H(\mathbf{x})={\mu}_{a}(\mathbf{x})\mathrm{\Phi}[\mathbf{x},{\mu}_{a}(\mathbf{x}),D(\mathbf{x})]\mathrm{.}$$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.^{45} The equation for the change in the fluence now contains an additional term for the optical diffusion perturbation (related to the scattering),

## (43)

$${(\mu}_{a,0}-\nabla \xb7{D}_{0}\nabla )\delta \mathrm{\Phi}(\mathbf{x})=[\nabla \xb7\delta D(\mathbf{x})\nabla -\delta {\mu}_{a}(\mathbf{x})]\mathrm{\Phi}(\mathbf{x}),$$## (44)

$$\mathrm{\Phi}={\mathrm{\Phi}}_{0}-{\int}_{{\mathbb{R}}^{n}}{G}_{0}(\mathbf{x},{\mathbf{x}}^{\prime})\delta {\mu}_{a}({\mathbf{x}}^{\prime}){\mathrm{\Phi}}_{0}({\mathbf{x}}^{\prime})\mathrm{d}{\mathbf{x}}^{\prime}\phantom{\rule{0ex}{0ex}}-{\int}_{{\mathbb{R}}^{n}}\delta D({\mathbf{x}}^{\prime})\nabla {G}_{0}(\mathbf{x},{\mathbf{x}}^{\prime})\xb7\nabla {\mathrm{\Phi}}_{0}({\mathbf{x}}^{\prime})\mathrm{d}{\mathbf{x}}^{\prime}.$$## (45)

$$\delta H(\mathbf{x})\approx \delta {\mu}_{a}(\mathbf{x}){\mathrm{\Phi}}_{0}(\mathbf{x})\phantom{\rule{0ex}{0ex}}-{\mu}_{a,0}{\int}_{{\mathbb{R}}^{n}}{G}_{0}(\mathbf{x},{\mathbf{x}}^{\prime})\delta {\mu}_{a}({\mathbf{x}}^{\prime}){\mathrm{\Phi}}_{0}({\mathbf{x}}^{\prime})\mathrm{d}{\mathbf{x}}^{\prime}\phantom{\rule{0ex}{0ex}}-{\mu}_{a,0}{\int}_{{\mathbb{R}}^{n}}\delta D({\mathbf{x}}^{\prime})\nabla {G}_{0}(\mathbf{x},{\mathbf{x}}^{\prime})\xb7\nabla {\mathrm{\Phi}}_{0}({\mathbf{x}}^{\prime})\mathrm{d}{\mathbf{x}}^{\prime},$$^{45}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:

## (46)

$$\frac{{H}^{(1)}}{{H}^{(2)}}=\frac{{\mathrm{\Phi}}^{(1)}}{{\mathrm{\Phi}}^{(2)}}\approx \frac{{\mathrm{\Phi}}_{0}^{(1)}-\int {G}_{0}(\mathbf{x},{\mathbf{x}}^{\prime})\delta {\mu}_{a}({\mathbf{x}}^{\prime}){\mathrm{\Phi}}_{0}^{(1)}({\mathbf{x}}^{\prime})\mathrm{d}{\mathbf{x}}^{\prime}-\int \delta D({\mathbf{x}}^{\prime})\nabla {G}_{0}(\mathbf{x},{\mathbf{x}}^{\prime})\xb7\nabla {\mathrm{\Phi}}_{0}^{(1)}({\mathbf{x}}^{\prime})\mathrm{d}{\mathbf{x}}^{\prime}}{{\mathrm{\Phi}}_{0}^{(2)}-\int {G}_{0}(\mathbf{x},{\mathbf{x}}^{\prime})\delta {\mu}_{a}({\mathbf{x}}^{\prime}){\mathrm{\Phi}}_{0}^{(2)}({\mathbf{x}}^{\prime})\mathrm{d}{\mathbf{x}}^{\prime}-\int \delta D({\mathbf{x}}^{\prime})\nabla {G}_{0}(\mathbf{x},{\mathbf{x}}^{\prime})\xb7\nabla {\mathrm{\Phi}}_{0}^{(2)}({\mathbf{x}}^{\prime})\mathrm{d}{\mathbf{x}}^{\prime}}.$$## 5.2.2.

#### Direct “zero-divergence” method

Bal et al.^{44}^{,}^{104} 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}={\mu}_{a}{\mathrm{\Phi}}_{1}$, ${H}_{2}={\mu}_{a}{\mathrm{\Phi}}_{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 ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ (or ${\mu}_{a}$ and $D$), the unknown parameters are defined as:

## (48)

$${\mu}_{1}\u2254\frac{\sqrt{D}}{{\mu}_{a}},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{\mu}_{2}\u2254-\frac{{\nabla}^{2}\left(\sqrt{D}\right)}{\sqrt{D}}-\frac{{\mu}_{a}}{D}.$$Multiplying the diffusion equation for the fluence ${\mathrm{\Phi}}_{1}$ by ${\mathrm{\Phi}}_{2}$, and vice versa, and subtracting gives ${\mathrm{\Phi}}_{2}[{\mu}_{a}{\mathrm{\Phi}}_{1}-(\nabla \xb7D\nabla ){\mathrm{\Phi}}_{1}]-\phantom{\rule{0ex}{0ex}}{\mathrm{\Phi}}_{1}[{\mu}_{a}{\mathrm{\Phi}}_{2}-(\nabla \xb7D\nabla ){\mathrm{\Phi}}_{2}]=0$, which can be rearranged to $\nabla \xb7[D{\mathrm{\Phi}}_{1}^{2}\nabla ({\mathrm{\Phi}}_{2}/{\mathrm{\Phi}}_{1})]=0$. Recalling $H={\mu}_{a}\mathrm{\Phi}$ and defining the vector field

## (49)

$$\beta ={H}_{1}\nabla {H}_{2}-{H}_{2}\nabla {H}_{1}={H}_{1}^{2}\nabla \left(\frac{{H}_{2}}{{H}_{1}}\right)={\mu}_{a}^{2}{\mathrm{\Phi}}_{1}^{2}\nabla \left(\frac{{\mathrm{\Phi}}_{2}}{{\mathrm{\Phi}}_{1}}\right),$$The procedure to reconstruct images of ${\mu}_{a}$ and $D$ is then the following: (1) Construct the vector field $\mathbf{\beta}$ from the two measured fields ${H}_{1}$ and ${H}_{2}$ by using Eq. (49). (2) Solve the transport equation, Eq. (50), for ${\mu}_{1}(\mathbf{x})$ that is, find the scalar field that makes the divergence of ${\mu}_{1}^{2}\mathbf{\beta}$ vanish. This is nontrivial and will not be possible for all fields $\mathbf{\beta}$, 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 ${\mu}_{2}=-{\nabla}^{2}({\mu}_{1}{H}_{1,2})/\phantom{\rule{0ex}{0ex}}({\mu}_{1}{H}_{1,2})$, which can be found by rearranging the definition in the second of Eqs. (48). (4) Solve $({\nabla}^{2}+{\mu}_{2})\sqrt{D}=-1/{\mu}_{1}$ for $\sqrt{D}$, which also comes from Eqs. (48). (5) The final step is to obtain ${\mu}_{a}=\sqrt{D}/{\mu}_{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^{104} discretize Eq. (50) and find ${\mu}_{1}$ by using a Bregman iteration to minimize a functional consisting of the norm of the divergence of ${\mu}_{1}^{2}\mathbf{\beta}$ and a regularizing total variation penalty term. They used a low pass filter to ameliorate the difficulties in the calculation of $\mathbf{\beta}$ caused by dividing by small values. A more visual, but less practical, method of solving Eq. (50) was proposed by Bal and Uhlmann.^{44} Using the chain rule, we may rewrite Eq. (50) as

## (52)

$${\mu}_{1}[{\mathbf{x}}^{*}(s)]={\mu}_{1}^{0}\text{\hspace{0.17em}}\mathrm{exp}(-\frac{1}{2}{\int}_{0}^{s}\gamma [\mathbf{x}({s}^{\prime})]\mathrm{d}{s}^{\prime}),$$## 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,^{96} 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, ${\mu}_{a}$, $D\in {\mathbb{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}^{\mathrm{obs}}\in {\mathbb{R}}^{M}$. The task of recovering the unknown coefficients (${\mu}_{a}$, $D$) becomes the task of adjusting the coefficients in order to minimize the functional

## (53)

$$\underset{{\mu}_{a},D}{\mathrm{argmin}}\text{\hspace{0.17em}}\epsilon =\frac{1}{2}{\Vert H({\mu}_{a},D)-{H}^{\mathrm{obs}}\Vert}^{2}+\mathcal{P}({\mu}_{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 (${\mu}_{a}^{(0)}$, ${D}^{(0)}$) by choosing new values ${\mu}_{a}^{(k+1)}={\mu}_{a}^{(k)}+\delta {\mu}_{a}^{(k)}$, ${D}^{(k+1)}={D}^{(k)}+{\delta}^{(k)}$ so that at each iteration the value of $\epsilon [{\mu}_{a}^{(k+1)},{D}^{(k+1)}]$ is smaller than $\epsilon [{\mu}_{a}^{(k)},{D}^{(k)}]$. The hope is that eventually the minimum of the multidimensional surface of $\epsilon ({\mu}_{a},D)$ will be reached, thereby giving the values of $({\mu}_{a},D)$ at which $H({\mu}_{a},D)$ is the best match to ${H}^{\mathrm{obs}}$ (given the constraints $\mathcal{P}$).

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

## (54)

$${g}_{a}=\frac{\partial \epsilon}{\partial {\mu}_{a}}={J}_{a}^{T}(H-{H}^{obs})+\frac{\partial \mathcal{P}}{\partial {\mu}_{a}},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{g}_{D}=\frac{\partial \epsilon}{\partial D}={J}_{D}^{T}(H-{H}^{obs})+\frac{\partial \mathcal{P}}{\partial D},$$## (55)

$${J}_{a}=\frac{\partial H}{\partial {\mu}_{a}}=\mathrm{\Phi}\delta +{\mu}_{a}(\frac{\partial \mathrm{\Phi}}{\partial {\mu}_{a}}),\phantom{\rule{0ex}{0ex}}{J}_{D}=\frac{\partial H}{\partial D}={\mu}_{a}\left(\frac{\partial \mathrm{\Phi}}{\partial D}\right),$$## (56)

$$\mathcal{H}=\frac{{\partial}^{2}\epsilon}{\partial {\mu}_{i}\partial {\mu}_{j}}\approx {J}^{T}J,\phantom{\rule[-0.0ex]{2em}{0.0ex}}J=[\begin{array}{cc}{J}_{a}& {J}_{D}\end{array}].$$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 \mathrm{\Phi}/\partial {\mu}_{a}$ and $\partial \mathrm{\Phi}/\partial D$ can be calculated from the sensitivity equations:

## (57)

$$({\mu}_{a}-\nabla \xb7D\nabla )\frac{\partial \mathrm{\Phi}}{\partial {\mu}_{a}}=-\mathrm{\Phi}\delta ,\phantom{\rule{0ex}{0ex}}({\mu}_{a}-\nabla \xb7D\nabla )\phantom{\rule{0ex}{0ex}}\frac{\partial \mathrm{\Phi}}{\partial D}=\nabla \xb7\delta \nabla \mathrm{\Phi}.$$## (58)

$$\left[\begin{array}{c}\delta {\mu}_{a}\\ \delta D\end{array}\right]=-\alpha {({J}^{T}J)}^{-1}\left[\begin{array}{c}{g}_{a}\\ {g}_{D}\end{array}\right],\phantom{\rule[-0.0ex]{2em}{0.0ex}}g=\left[\begin{array}{c}{g}_{a}\\ {g}_{D}\end{array}\right],$$## 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}^{\mathrm{obs}}(\mathbf{x},\lambda )$. The error functional therefore becomes

## (59)

$$\underset{{C}_{k},{\mu}_{s}^{\prime}}{\mathrm{argmin}}\text{\hspace{0.17em}}{\epsilon}_{\lambda}=\frac{1}{2}{\Vert H({C}_{k},{\mu}_{s}^{\prime})-{H}^{\mathrm{obs}}\Vert}^{2}+\mathcal{P}({C}_{k},{\mu}_{s}^{\prime}),$$## (60)

$$\frac{\partial}{\partial {C}_{k}}=\frac{\partial {\mu}_{a}}{\partial {C}_{k}}\frac{\partial}{\partial {\mu}_{a}}={\alpha}_{k}\frac{\partial}{\partial {\mu}_{a}}\mathrm{.}$$If the scattering is assumed to depend on wavelength as ${\mu}_{s}^{\prime}\approx a(\mathbf{x}){\lambda}^{-b}$, then the gradients can also be simply calculated:

## (61)

$$\frac{\partial}{\partial a}=\frac{\partial {\mu}_{s}^{\prime}}{\partial a}\frac{\partial}{\partial {\mu}_{s}^{\prime}}={\lambda}^{-b}\frac{\partial}{\partial {\mu}_{s}^{\prime}}.$$Using these expressions, Cox et al.^{42} 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, ${\alpha}_{k}(\lambda )$, 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(\mathbf{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, $2N$ times per iteration step, it is only necessary to run the forward and the adjoint models once each:^{96}

## (62)

$$({\mu}_{a}-\nabla \xb7D\nabla )\mathrm{\Phi}=q,\phantom{\rule[-0.0ex]{2em}{0.0ex}}\phantom{\rule{0ex}{0ex}}({\mu}_{a}-\nabla \xb7D\nabla ){\mathrm{\Phi}}^{*}={\mu}_{a}(H-{H}^{obs}).$$## (63)

$${g}_{a}=\mathrm{\Phi}(H-{H}^{obs})-\mathrm{\Phi}{\mathrm{\Phi}}^{*},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{g}_{D}=-\nabla \mathrm{\Phi}\xb7\nabla {\mathrm{\Phi}}^{*}.$$Using a gradient-based split Bregman algorithm, Gao et al.^{105}^{,}^{106} 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^{104} 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.^{107} 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.^{89} 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.^{108} 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 ${\mu}_{a}={T}^{-1}(H)$ then find the concentrations ${C}_{k}={L}_{\lambda}^{-1}({\mu}_{a})$, whereas Eq. (6) suggests a single inversion, ${C}_{k}={T}_{\lambda}^{-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 ${\mu}_{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.^{109} 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.^{110} 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.^{111}^{,}^{112} 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/{\mu}_{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.^{113} 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{\mathrm{\Gamma}}$ 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}^{\prime}$, pressure ${p}^{\prime}$, and density ${\rho}^{\prime}$. 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{\mathrm{\Gamma}}$ can be related to the Grüneisen parameter, $\Gamma $, which is a thermodynamic property of a material. Consider the case in which the absorbed power density, $\mathcal{H}$, the rate at which the optical energy is absorbed per unit volume per unit time, is separable so $\mathcal{H}(\mathbf{x},t)=H(\mathbf{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

## (64)

$$\frac{\partial {p}^{\prime}}{\partial t}=\left(\frac{1}{\rho {k}_{T}}\right)\frac{\partial {\rho}^{\prime}}{\partial t}+\left(\frac{\beta}{{k}_{T}}\right)\frac{\partial {T}^{\prime}}{\partial t},$$*stress confinement*. Under this isochoric condition the rate of increase in pressure can be written $\partial {p}^{\prime}/\partial t=[(\beta /{k}_{T})\partial {T}^{\prime}]/\partial t$. The rate at which the temperature increases is given by $\partial {T}^{\prime}/\partial t=\mathcal{H}/\rho {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

## (65)

$$\frac{\partial {p}^{\prime}}{\partial t}=\left(\frac{\beta}{\rho {C}_{v}{k}_{T}}\right)\mathcal{H}(\mathbf{x},\mathit{t})=\mathrm{\Gamma}H(\mathbf{x})f(\mathit{t}),$$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 {\rho}^{\prime}/\partial t\ll \partial {T}^{\prime}/\partial t$ does not hold, then $\hat{\mathrm{\Gamma}}<\mathrm{\Gamma}$. The latter will be true in pulses for which $\mathcal{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.^{114} It has been shown that this rate can be increased for gold nanoparticles by coating them with silica, thereby increasing $\hat{\mathrm{\Gamma}}$ and the magnitude of the PA signal.^{115}

## 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.,^{86}^{,}^{45} Guo et al.,^{52} and Bal et al.^{44}^{,}^{104} 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{\mathrm{\Gamma}}$ is to recover it simultaneously with the chromophore concentrations and scattering. Bal and Ren^{104} 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.^{45}

## 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}(\mathbf{x})=\hat{\mathrm{\Gamma}}H(\mathbf{x})$. However, $\hat{\mathrm{\Gamma}}$ may take different values in different tissue types: the Grüneisen parameter of blood has been measured to be^{65} $\Gamma (\text{blood})=0.14$, while in liver samples on average^{116} $\Gamma (\text{liver})=0.12$, which is greater than 15% difference. (For comparison, $\Gamma (\text{water})=0.11$.) If a contrast agent is used, it is even more likely that the PA efficiency will vary spatially, $\hat{\mathrm{\Gamma}}=\hat{\mathrm{\Gamma}}(\mathbf{x})$. Putting this spatial dependence into Eqs. (1) and (6) gives an expression for the initial pressure distribution as

## (67)

$${p}_{0}=\hat{\mathrm{\Gamma}}(\mathbf{x})\mathrm{\Phi}(\mathbf{x},\lambda )\sum _{k=1}^{K}{C}_{k}(\mathbf{x}){\alpha}_{k}(\lambda ).$$Laufer et al.^{117} suggest there may be a different PA efficiency for each chromophore which depends on that chromophore’s concentration, ${\hat{\mathrm{\Gamma}}}_{k}={\hat{\mathrm{\Gamma}}}_{k}[{C}_{k}(\mathbf{x})]$. The initial pressure distribution then becomes

## (68)

$${p}_{0}=\mathrm{\Phi}(\mathbf{x},\lambda )\sum _{k=1}^{K}{\hat{\mathrm{\Gamma}}}_{k}{C}_{k}(\mathbf{x}){\alpha}_{k}(\lambda ).$$## 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 ${\mathrm{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 ${\mathrm{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.