* Now at Julius-Wolff-Institut, Charité Universitätsmedizin Berlin, Augustenburger Platz 1, 13353 Berlin, Germany.

# Quantitative spectroscopic photoacoustic imaging: a review

**Ben Cox**,

**Jan G. Laufer**,

**Paul C. Beard**

University College London, Department of Medical Physics and Bioengineering, Gower Street, London WC1E 6BT, United Kingdom

**Simon R. Arridge**

University College London, Department of Computer Science, Gower Street, London WC1E 6BT, United Kingdom

*J. Biomed. Opt*. 17(6), 061202 (Jun 01, 2012). doi:10.1117/1.JBO.17.6.061202

#### Open Access

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

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

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.

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.

There are essentially two aspects to PA imaging: acoustic and optical. The focus of this paper is on the optical inversion, and the acoustic inverse problem (PA image reconstruction) is not covered here. It has been studied extensively, and we refer the interested reader to Refs. ^{4}^{5}^{6}^{7}^{8}^{9}–^{10} and to the references contained therein. This review does not discuss the recovery of properties of the tissue other than optical properties (e.g., sound speed, density, and acoustic attenuation estimation) although the accurate determination of these acoustic quantities could have an indirect effect on QPAI if they are used to improve the accuracy of the acoustic image reconstruction. A great deal of work has been done in recent years on finding and designing suitable contrast agents for PAI, but the focus here is on techniques to facilitate the estimation of their concentrations rather than on the chemistry of the substances themselves. Techniques to determine the optimal choice of wavelengths for QPAI is outside the scope of this article as it will be very dependent on the particular arrangement (the object, imaging mode, main source of contrast, other chromophores present, etc.). This question has been tackled for some specific scenarios.^{11}^{–}^{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.

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.

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.

There are several mechanisms by which light can be used to generate sound waves.^{15}^{–}^{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 $vibrational/collisional$ relaxation on a sub-nanosecond timescale. This localized injection of heat will lead to a small rise in the local temperature and a related rise in the local pressure, $p0$. 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

*initial acoustic pressure distribution*, $p0$. The absorbed energy density can be written as

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.

**F1 :**

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

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

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 $\u223c45\u2009$ *μ*m to a depth of a few millimeters, with the vertical resolution limited to half the shortest wavelength detectable ($\u223c15$ *μ*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 $\u223c10$ *μ*m should be achievable.^{22}

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.

A set of PA images $p0(x,\lambda )$ are obtained experimentally, each at a different optical wavelength $\lambda $. (The position vector $x$ will in general be in three dimensions although it can be reduced to one or two under certain conditions.) If the system is calibrated and the PA efficiency is known, then the multiwavelength set of images of absorbed energy density $H(x,\lambda )=\Gamma ^\u22121p0(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 $Ck(x)$ of $K$ chromophores with known molar absorption coefficient spectra $\alpha k(\lambda )$, given the absorbed energy images $H(x,\lambda )$ when the chromophores are linked to the absorption coefficients with the linear mapping $L\lambda $:

The light fluence, $\Phi (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

Equations (4) and (5) suggest a two-stage inversion strategy: first recover the absorption coefficients $\mu a=T\u22121(H)$ and then find the concentrations $Ck=L\lambda \u22121(\mu a)$, whereas Eq. (6) suggests a single inversion, $Ck=T\lambda \u22121(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

**F2 :**

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

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

A cursory glance at Eq. (6) might suggest that there is a linear relationship between the chromophore concentrations $Ck$ and the absorbed energy density $H$. If this were the case, then it would allow the chromophore concentrations $Ck$ 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 $\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 )\u221d\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 $\Phi (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(x)=\mu a(x)\Phi (x)$. The unknown and nonuniform light distribution $\Phi (x)$ will therefore distort what would, in its absence, have been an image proportional to $\mu a(x)$, the absorption coefficient as a function of position. This cannot be avoided by ensuring that the illumination of the target is uniform at the surface, because as soon as there is any absorption or scattering in the medium, the fluence will vary spatially. When the absorption or scattering is spatially varying (as is often the case in tissue), the unknown fluence also will be. This structural distortion and the spectral coloring are two manifestations of the same phenomenon, namely the effect of the fluence on the PA image.^{23}

The measured PA signals $p(t)$ or images $p0$ can be related to the chromophore concentrations $Ck$ 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 $Ck$. 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.

The essential difficulty with QPAI is the unknown light fluence. Inverting $T\lambda $ or $T$ to obtain $\mu a$ or $Ck$ 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}^{–}^{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.

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:

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.

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 $(s^\xb7\u2207+\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 $\u2202\Phi /\u2202z=\u2212\mu a\Phi $. If the source illuminates the surface with a fluence of $\Phi 0$ and $\mu a$ is constant, then the fluence as a function of depth may be written

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 $PN$ 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 $P1$ approximation can be written as the two equations^{30}

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

The condition $\mu s\u2032\u226b\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

*effective attenuation coefficient*. This solution will be used in some of the inversions in Sec. 5. In one dimension (1-D), the solution is $\Phi (z)=\Phi 0\u2009exp(\u2212\mu effz)$ where $\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

For a collimated incident beam, the condition for the accuracy of the DA, $\mu s\u2032\u226b\mu a$, must be supplemented with a second condition. As the collimated part of an incident beam will decay exponentially with depth $z$ as $exp(\u2212\mu tz)$, the fluence will only be dominated by diffuse photons at depths $z\u226b1/\mu t$. The optical coefficients in PA imaging applications are of the order of $\mu a=0.1\u2009\u2009mm\u22121$, $\mu s=10\u2009\u2009mm\u22121$, $g=0.9$, $\mu s\u2032=1\u2009\u2009mm\u22121$ so the condition that $\mu s\u2032\u226b\mu a$ is met, but the condition that the fluence is dominated by diffuse photons, only holds for distances $z\u226b1/\mu t\u22480.1\u2009\u2009mm$, so perhaps for $z>1\u2009\u2009mm$ 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 $PN$ 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}^{–}^{35}

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 $\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.

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.

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.

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=\u2212ln[\Phi (z)/\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.

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/[\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.

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.

Least-squares minimization is a common and well-developed framework for solving inverse problems. If the data is a measured (observed) absorbed energy distribution, $Hobs$, then the squared error between it and the output of a forward model, $Hmodel(Ck)$, is $\Vert Hmodel(Ck)\u2212Hobs\Vert 2=\u2211i[Himodel(Ck)\u2212Hiobs]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 $Ck$:

Solving Eq. (18) means finding the minimum of the scalar functional $\epsilon $. Gradient-based methods use the functional gradient (the vector of first derivatives), $gi=\u2202\epsilon /\u2202Ck,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, $Hij=\u22022\epsilon /(\u2202Ck,i\u2202Ck,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 $Jij=\u2202Himodel/\u2202Ck,j$ to estimate the Hessian matrix as $H\u2248JTJ$, 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}

Estimating the chromophore concentrations $Ck$ 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\u22121$, $T\lambda \u22121$, 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}

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 $107$ or more voxels, there may be of the order of $108$ 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 $1016$ 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 $An$, then the absorption coefficient can be written as the sum

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 $An$ 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.

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.

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

^{49}

^{–}

^{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 )|=(\Gamma ^\Phi 0\mu a/2)[(\mu ac0)2+\omega 2]\u221212$, the ratio of two measurements made at different optical wavelengths, $\lambda 1$ and $\lambda 2$, may be written as

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 ($\u223c5\u2009$ *μ*m) much smaller than a typical vessel diameter ($\u223c30$ *μ*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 $\u223c5$ *μ*m and a depth of perhaps one order of magnitude greater will not generate purely plane waves propagating in the $z$ direction, even within the acoustic focus. There will be components propagating at angles to the $z$-axis, so the detected acoustic wave will therefore not vary exponentially according to $exp(\u2212\mu ac0t)$, and Eq. (22) will not hold.