Special Section on Photoacoustic Imaging and Sensing

Quantitative spectroscopic photoacoustic imaging: a review

[+] Author Affiliations
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
History: Received October 3, 2011; Revised March 13, 2012; Accepted March 21, 2012
Text Size: A A A

Open Access Open Access

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

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.

Figures in this Article

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.13 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.

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.

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.

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. 45678910 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.1113 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.

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.

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.

The photoacoustic effect

There are several mechanisms by which light can be used to generate sound waves.1517 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 Display Formula

p0=Γ^H,(1)
where Γ^ is the PA efficiency. (This efficiency can be identified with the Grüneisen parameter Γ for an absorbing fluid; Sec. 6.3). It is this increase in pressure that subsequently becomes a propagating acoustic (ultrasonic) wave, because of the elastic nature of tissue. It is therefore referred to as the initial acoustic pressure distribution, p0. The absorbed energy density can be written as Display Formula
H=μaΦ,(2)
where μa is the optical absorption coefficient (the total absorption due to the chromophores) and Φ is the light fluence (the radiance integrated over all directions and time; Sec. 2.2). The initial acoustic pressure distribution may therefore be written as the product of three quantities: Display Formula
p0=Γ^μaΦ.(3)

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.

Grahic Jump LocationF1 :

PA signal generation. Spatially varying chromophore concentrations (naturally occurring or contrast agents) give rise to optical absorption in the medium. The absorption and scattering coefficients, μa and μs, determine the fluence distribution Φ (how the light from a source becomes distributed in the tissue), and hence the absorbed energy distribution H. This energy generates a pressure distribution 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).

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.410 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.

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 45 μm to a depth of a few millimeters, with the vertical resolution limited to half the shortest wavelength detectable (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 10 μm should be achievable.22

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.

Formulation of the problem

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

μa(x,λ)=Lλ(Ck)=k=1KCk(x)αk(λ),(4)
and the absorption coefficients are linked to the absorbed energy images by the nonlinear mapping T: Display Formula
H(x,λ)=T(μa,μs)=μa(x,λ)Φ(x,λ;μa,μs).(5)

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

H(x,λ)=Tλ(Ck,μs)=Φ(x,λ;Ck,μs)k=1KCk(x)αk(λ).(6)

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

p(t)=AΓ^Tλ(Ck)=W(Ck),(7)
where A is the acoustic forward operator. Linearizations of W have been found useful in certain very specific circumstances (see Sec. 4), but in general the nonlinearity of the light transport/absorption model must be taken into account.

Grahic Jump LocationF2 :

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, A1; thermoelastic scaling, Γ^1; nonlinear optical inversion for the optical coefficients,T1; and finally ,Lλ1 a linear spectroscopic inversion of the absorption coefficient spectra to recover the chromophore concentrations. Alternatively, these last two stages may be combined into the operator Tλ1. A different approach is to attempt a one-step inversion of the whole process, W1. This may be linearized under specific and rather restricted conditions. The inverse light transport operators T1 and its multiwavelength equivalent Tλ1 are at the core of the methods described in this paper.

For OR-PAM images obtained by scanning the focused illumination spot across the region of interest, the fluence distribution Φ(x) will be different for each measurement, so Eq. (5) should be replaced by H(xi)=μa(xi)Φi[xi,μa(x)] where Φ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.

Spectral coloring and structural distortion

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(λ) straightforwardly, for example, using a linear matrix inversion. Unfortunately the unknown fluence Φ depends in a nontrivial way on the chromophore concentrations because Tλ is both nonlinear and nonlocal (H at one point can be affected by μa some distance away via its effect on the fluence). It is not reasonable, therefore, to assume that H(λ)μa(λ) because the spectrum of the absorbed energy density at a given point will be affected by the fluence, which will have been colored by its passage through the medium. Because biological tissue is highly scattering, each photon will take a long and convoluted path through the tissue. The spectra of the fluence and the absorption are therefore intertwined, and H(λ) will be a complicated combination of the absorption spectra of the constituent chromophores at each point in the irradiated volume. (This coloring of the spectrum is the same effect as when seeing the world through colored glass.)

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

Using a model to recover the concentrations

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λ, 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.

Fluence Models

The essential difficulty with QPAI is the unknown light fluence. Inverting Tλ or T to obtain μ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.2427 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.

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: Display Formula

1cϕt(x,s^,t)=q(x,s^,t)[s^·+μa(x)+μs(x)]ϕ(x,s^,t)+μsΘ(s^,s^)ϕ(x,s^,t)ds^,(8)
where ϕ(x,s^,t) is the light radiance, Θ(s^,s^) is called the scattering phase function and is the probability that a photon originally traveling in direction s^ ends up traveling in direction s^ if scattered, q(x,s^,t) is a source of photons, and c is the speed of light in the medium. The terms on the right-hand side of this equation account for the fact that the rate of change of the number of photons within a small region around the point x and traveling in direction s^ could be due to (1) sources q; (2) net outflow of photons due to the radiance gradient, s^·; (3) photons absorbed, μaϕ; (4) photons scattered into another direction, μsϕ; or (5) photons scattered into direction s^ from another direction (given by the phase integral). Wave effects, polarization, radiative processes, ionization, inelastic scattering, and reactions are all neglected in this model.

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

(s^·+μt)ϕ(x,s^)μsΘ(s^,s^)ϕ(x,s^)ds^=q(x,s^),(9)
where the total attenuation coefficient μt=μa+μs, and ϕ(x,s^) is now used to represent the time-integrated light radiance. The fluence Φ is the integral of the radiance over all angles s^: Display Formula
Φ(x)=ϕ(x,s^)ds^.(10)

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^·+μa)ϕ=q. In the case of a collimated source propagating in the z direction, the fluence, which equals the radiance as there is only one direction of propagation, satisfies the equation Φ/z=μaΦ. If the source illuminates the surface with a fluence of Φ0 and μa is constant, then the fluence as a function of depth may be written Display Formula

Φ(z)=Φ0exp(μaz),(11)
which is known as Beer’s law or the Beer-Lambert-Bouguer law.

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 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 equations30Display Formula

μaΦ+·F=q0,FD+Φ=3q1,(12)
where D=[3(μa+μs)]1 is the optical diffusion coefficient, μs=(1g)μs is the reduced scattering coefficient (with g the anisotropy factor), and the vector flux, F, is defined as Display Formula
F(x)=ϕ(x,s^)s^ds^.(13)

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): Display Formula

μaΦ·(D)Φ=q0.(14)

The condition μsμa is usually considered sufficient to ensure the accuracy of the diffusion approximation away from sources, because it ensures the scattered fluence is almost isotropic. For arbitrary heterogeneous media, Eq. (14) may be solved by standard numerical methods such as finite elements.31,32 For a homogeneous medium, the solution to Eq. (14) is given by Display Formula

Φ(x)=G0(xx)q0(x)dx(15)
with the (3-D) free space Green’s function, G0 [the solution to (μaD2)G0=δ(xx)] given by Display Formula
G0(x,x)=exp(μeff|xx|)4πD|xx|,μeff=3μa(μa+μs),(16)
where μeff is called the effective attenuation coefficient. This solution will be used in some of the inversions in Sec. 5. In one dimension (1-D), the solution is Φ(z)=Φ0exp(μeffz) where Φ0 is fluence incident on the surface, which is analogous to Eq. (11). Some care is required here because while this is the solution to Eq. (14) in 1-D, it only applies when the fluence everywhere is diffuse, which is rarely true close to the surface; where light enters the tissue, it often remains partially collimated. In practice, the fluence will often reach a peak at some distance beneath the surface, where both the incident and backscattered light are contributing to the fluence, before decaying exponentially. Introducing a scaling factor to account for this gives Display Formula
Φ(z)kΦ0exp(μeffz)forz1μt,(17)
which is used in the inversions in Sec. 3.3.

Collimated source, δ-Eddington approximation

For a collimated incident beam, the condition for the accuracy of the DA, μsμa, must be supplemented with a second condition. As the collimated part of an incident beam will decay exponentially with depth z as exp(μtz), the fluence will only be dominated by diffuse photons at depths z1/μt. The optical coefficients in PA imaging applications are of the order of μa=0.1mm1, μs=10mm1, g=0.9, μs=1mm1 so the condition that μsμa is met, but the condition that the fluence is dominated by diffuse photons, only holds for distances z1/μt0.1mm, so perhaps for z>1mm 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 δ-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.3335

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 Θ.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.

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.

Inversion techniques
Linearization

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

Direct inversion

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

Fixed-point iteration

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

Model-based minimization

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

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, Hobs, then the squared error between it and the output of a forward model, Hmodel(Ck), is Hmodel(Ck)Hobs2=i[Himodel(Ck)Hiobs]2, where the subscript i indicates the value at the i’th voxel. The aim in least-squares minimization is to minimize the error functional, ε, by adjusting the unknowns Ck: Display Formula

argminCkε=12Hmodel(Ck)Hobs2+P(Ck).(18)
The second term, P, is a penalty term that can be used to include regularization or other constraints. Even when ε has a well-defined minimum in the noise-free case, where the difference between Hmodel and Hobs is small, the noise in the data may mean there may be several solutions, Hmodel, that fit the measured data, Hobs, equally well. Regularization is a general term to describe methods for ameliorating this problem, and common approaches are to encourage the unknowns—the absorption or scattering coefficients, say—to be smooth (Tikhonov) or piecewise constant (total variation). The latter may be reasonable in some cases from physiological considerations; for example, when imaging a plane through a blood vessel, the absorption coefficient may take one value inside the vessel and one outside. Other prior knowledge of some aspect of the unknowns, or perhaps of an intermediate quantity, such as shape or an empirically determined maximum or minimum, may also be included in this way. The basic idea is to reduce the solution space by excluding functions with certain properties, such as those that are insufficiently smooth, from being considered as solutions.

Solving Eq. (18) means finding the minimum of the scalar functional ε. Gradient-based methods use the functional gradient (the vector of first derivatives), gi=ε/Ck,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=2ε/(Ck,iCk,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=Himodel/Ck,j to estimate the Hessian matrix as HJTJ, 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

Uniqueness and ill-posedness

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 ε 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λ 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, T1, Tλ1, will be to grow the high frequency components. This will have the effect of amplifying the noise in the measured data which, unchecked, may come to dominate the inversion. Alterations to the inverse operator that are designed to reduce this unwanted effect are termed regularization (as mentioned in the context of penalty functionals in Sec. 2.3.2 above). This type of phenomenon is common to many inverse problems, both linear and nonlinear, and a large number of regularization techniques are described in the literature.46

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 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 light47,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 Display Formula

μa(x)=nμa,nSn(x),Sn(x)={1,xAn0,otherwise.(19)

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.

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 Display Formula

H(z)={μaΦ0exp(μaz),z00,z<0,(20)
where Φ0 is the fluence at the illuminated surface of the tissue, which is assumed, without loss of generalization, to be at z=0. This initial pressure will give rise to two PA waves, traveling in the +z and z directions, each with half the amplitude but retaining the exponential shape. For a pressure detector at z=0 (backward mode detection), then the signal reaching the detector will be simply Display Formula
p(t)=(Γ^Φ0μa2)exp(μac0t),t0,(21)
where c0 is the sound speed. (We are ignoring acoustic reflections at the surface here, but they can be straightforwardly included.) The absorption coefficient, μa, can be estimated either from the maximum amplitude of the signal, if the PA efficiency, Γ^, and the incident fluence, Φ0, are known, or by fitting a curve to the exponentially decaying slope.4951

More recently, a frequency domain method for quantification of μ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(ω)|=(Γ^Φ0μa/2)[(μac0)2+ω2]12, the ratio of two measurements made at different optical wavelengths, λ1 and λ2, may be written as Display Formula

|P(λ1;ω)||P(λ2;ω)|=Φ0(λ1)Φ0(λ2){c02+[ω/μa(λ2)]2c02+[ω/μa(λ1)]2}1/2,(22)
where the system transfer function, the acoustic absorption term, and the PA efficiency have canceled out. The three unknown numbers μa(λ1), μa(λ2), and the surface fluence ratio Φ0(λ1)/Φ0(λ2) may be obtained by curve-fitting Eq. (22) to measured acoustic spectra.

Guo et al. used this technique for quantification by OR-PAM, which raises the question as to how a model that explicitly varies only with depth, z, can be justified for use with OR-PAM, where there is variation in lateral dimensions too. There are three conditions that must be satisfied for the 1-D planar assumption to hold and, therefore, for Eq. (22) to be true: (1) the medium properties, (2) the absorbed energy distribution, and (3) the acoustic propagation must all be planar (1-D) on the scale of interest. The first condition is satisfied as the light beam in OR-PAM is focused to a spot (5 μm) much smaller than a typical vessel diameter (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 5 μm and a depth of perhaps one order of magnitude greater will not generate purely plane waves propagating in the z direction, even within the acoustic focus. There will be components propagating at angles to the z-axis, so the detected acoustic wave will therefore not vary exponentially according to exp(μac0t), and Eq. (22) will not hold.