A survey of computational frameworks for solving the acoustic inverse problem in three-dimensional photoacoustic computed tomography

TL;DR

This survey reviews modern optimization-based computational frameworks for solving the 3D photoacoustic tomography inverse problem, highlighting models, methods, and challenges in image reconstruction.

Contribution

It provides a comprehensive overview of current optimization-based approaches and models for 3D PACT image reconstruction, including non-conventional formulations estimating medium parameters.

Findings

Analysis of continuous and discrete imaging models

Review of optimization techniques for PACT reconstruction

Discussion of joint estimation of acoustic parameters

Abstract

Photoacoustic computed tomography (PACT), also known as optoacoustic tomography, is an emerging imaging technique that holds great promise for biomedical imaging. PACT is a hybrid imaging method that can exploit the strong endogenous contrast of optical methods along with the high spatial resolution of ultrasound methods. In its canonical form that is addressed in this article, PACT seeks to estimate the photoacoustically-induced initial pressure distribution within the object. Image reconstruction methods are employed to solve the acoustic inverse problem associated with the image formation process. When an idealized imaging scenario is considered, analytic solutions to the PACT inverse problem are available; however, in practice, numerous challenges exist that are more readily addressed within an optimization-based, or iterative, image reconstruction framework. In this article, the PACT image reconstruction problem is reviewed within the context of modern optimization-based image reconstruction methodologies. Imaging models that relate the measured photoacoustic wavefields to the sought-after object function are described in their continuous and discrete forms. The basic principles of optimization-based image reconstruction from discrete PACT measurement data are presented, which includes a review of methods for modeling the PACT measurement system response and other important physical factors. Non-conventional formulations of the PACT image reconstruction problem, in which acoustic parameters of the medium are concurrently estimated along with the PACT image, are also introduced and reviewed.

Figures

Equations

$$\alpha (\mathbf{r},f)={\alpha }_0(\mathbf{r})f^y,$$

$$\frac{\partial }{\partial t}\dot{\mathbf{s}}(\mathbf{r},t)$$

$$\frac{\partial }{\partial t}\rho (\mathbf{r},t)$$

$$p(\mathbf{r},t)=c_0(\mathbf{r}{)}^2$$

$$\displaystyle p(\mathbf{r},t)\Big{|}_{t=0}$$

$$\mu(\mathbf{r})=-2\alpha_{0}(\mathbf{r})c_{0}(\mathbf{r})^{y-1},\ \ \ \eta(\mathbf{r})=2\alpha_{0}(\mathbf{r})c_{0}(\mathbf{r})^{y}\tan\Big{(}\frac{\pi y}{2}\Big{)}.$$

$$p({\mathbf{r}}^{\prime },t)={\mathcal{H}}_{wave}{p}_0(\mathbf{r}),$$

$$p({\mathbf{r}}^{\prime },t)=\frac{1}{4\pi c_0^2}{\int }_V{\mathrm{d}}^3\mathbf{r}{p}_0(\mathbf{r})\frac{d}{dt}\frac{\delta \left(t-\frac{|{\mathbf{r}}^{\prime }-\mathbf{r}|}{c_0}\right)}{|{\mathbf{r}}^{\prime }-\mathbf{r}|}\equiv {\mathcal{H}}_{int}{p}_0(\mathbf{r}),$$

$$g({\mathbf{r}}^{\prime },t)={\int }_V{\mathrm{d}}^3\mathbf{r}{p}_0(\mathbf{r})\delta (c_{0}t-|{\mathbf{r}}^{\prime }-\mathbf{r}|),$$

$$p({\mathbf{r}}^{\prime },t)=\frac{1}{4\pi c_0^2}\frac{\partial }{\partial t}\left(\frac{g({\mathbf{r}}^{\prime },t)}{t}\right).$$

$${\mathcal{H}}_{CD}\equiv \mathcal{D}{\mathcal{H}}_{CC},$$

$$[\mathbf{u}{]}_{qK+k}=p({\mathbf{r}}^{\prime },t){|}_{{\mathbf{r}}^{\prime }={\mathbf{r}}_q^{\prime },t=k{\Delta }_t},$$

$$[\mathbf{u}{]}_{qK+k}=[{\mathcal{H}}_{CD}{p}_0(\mathbf{r}){]}_{qK+k}={\int }_{-\infty }^{\infty }dt{\tau }_k(t){\int }_{{\Omega }_0}d{\Omega }_{0}p({\mathbf{r}}_0,t){\sigma }_q({\mathbf{r}}_0),$$

$$p_0^a(\mathbf{r})=\sum _{n=1}^N[\mathbf{\theta }{]}_n{\varphi }_n(\mathbf{r}),$$

$$\mathbf{u}={\mathcal{H}}_{CD}{\mathbf{p}}_0^a(\mathbf{r})=\sum _{n=1}^N[\theta {]}_n{\mathcal{H}}_{CD}\{{\varphi }_n(\mathbf{r})\}\equiv \mathbf{H}\mathbf{\theta },$$

$$\phi^{\rm int}_{n}(\mathbf{r})\equiv\left\{\begin{array}[]{ll}(1-\frac{|x-x_{n}|}{\Delta_{s}})(1-\frac{|y-y_{n}|}{\Delta_{s}})(1-\frac{|z-z_{n}|}{\Delta_{s}}),&\text{if}\,|x-x_{n}|,|y-y_{n}|,|z-z_{n}|\leq\Delta_{s}\\ 0,&\text{otherwise}\end{array}\right.,$$

$$[\mathbf{\theta }{]}_n={\int }_V{\mathrm{d}}^3\mathbf{r}\delta (\mathbf{r}-{\mathbf{r}}_n)p_0(\mathbf{r}),n=0,1,\dots ,N-1,$$

$$\mathbf{u}={\mathbf{H}}_{int}{\mathbf{\theta }}_{int},$$

$${\mathbf{H}}_{int}\equiv \mathbf{DG}.$$

$$[\mathbf{G}\mathbf{\theta }{]}_{qK+k}\equiv {\Delta }_s^2\sum _{n=0}^{N-1}[\mathbf{\theta }{]}_n\sum _{i=0}^{N_i-1}\sum _{j=0}^{N_j-1}{\varphi }_n^{int}({\mathbf{r}}_{k,i,j})\equiv [\mathbf{g}{]}_{qK+k},$$

$$[\mathbf{Dg}{]}_{qK+k}\equiv \frac{1}{8\pi c_0^2{\Delta }_t^2}\left(\frac{[\mathbf{g}{]}_{qK+k+1}}{k+1}-\frac{[\mathbf{g}{]}_{qK+k-1}}{k-1}\right).$$

$$b(x)=\left\{\begin{array}[]{lc}\left(\sqrt{1-x^{2}/a^{2}}\right)^{m}\frac{I_{m}(\gamma\sqrt{1-x^{2}/a^{2}})}{I_{m}(\gamma)},&\text{if }0\leq x\leq a\\ 0,&\text{if }x>a,\end{array}\right.$$

$${\varphi }_n^{\mathrm{KB}}(\mathbf{r})=b(|\mathbf{r}-{\mathbf{r}}_n|),$$

$$\tilde{\mathbf{u}}={\tilde{\mathbf{H}}}_{KB}{\mathbf{\theta }}_{KB},$$

$$[{\tilde{\mathbf{H}}}_{KB}{]}_{qL+l,n}\equiv {\tilde{p}}^{KB}(f){\tilde{h}}_q^{s,point}({\mathbf{r}}_n,f){|}_{f=l{\Delta }_f},\text{ for }l=0,1,\dots ,L-1,$$

$${\tilde{p}}^{KB}(f)=-\frac{j4{\pi }^{2}f{a}^3\beta }{C_p{I}_m(\gamma )}\frac{{\hat{j}}_{m+1}(\sqrt{4{\pi }^2{a}^2{f}^2/c_0^2-{\gamma }^2})}{(4{\pi }^2{a}^2{f}^2/c_0^2-{\gamma }^2{)}^{(m+1)/2}},$$

$${\tilde{h}}_q^{s,point}({\mathbf{r}}_n,f)=\frac{\exp \left(-j2\pi f\frac{|{\mathbf{r}}_q^s-{\mathbf{r}}_n|}{c_0}\right)}{2\pi |{\mathbf{r}}_q^s-{\mathbf{r}}_n|}.$$

$${\left[\begin{array}{c}{\mathbf{v}}_0,{\mathbf{v}}_1,\cdots ,{\mathbf{v}}_{K-1}\end{array}\right]}^T={\mathbf{T}}_{K-1}\cdot {\mathbf{T}}_{K-2}\cdots {\mathbf{T}}_1{\left[\begin{array}{c}{\mathbf{v}}_0,0_{7N\times 1},\cdots ,0_{7N\times 1}\end{array}\right]}^T.$$

$${\left[\begin{array}{c}{\mathbf{v}}_0,0_{7N\times 1},\cdots ,0_{7N\times 1}\end{array}\right]}^T={\mathbf{T}}_0\mathbf{\theta },$$

$$\mathbf{u}=\mathbf{M}{\left[\begin{array}{c}{\mathbf{v}}_0,{\mathbf{v}}_1,\cdots ,{\mathbf{v}}_{K-1}\end{array}\right]}^T,$$

$$\mathbf{u}=\mathbf{M}{\mathbf{T}}_{K-1}\cdots {\mathbf{T}}_1{\mathbf{T}}_0\mathbf{\theta }.$$

$${\mathbf{H}}_{PS}=\mathbf{M}{\mathbf{T}}_{K-1}\cdots {\mathbf{T}}_1{\mathbf{T}}_0.$$

Full text

††: Physics in Medicine and Biology

1 Introduction to PACT

Photoacoustic computed tomography (PACT), also known as optoacoustic tomography, is a rapidly emerging imaging technique that holds great promise for biomedical imaging (?; ?; ?; ?; ?). PACT is a hybrid technique that exploits the photoacoustic effect for signal generation. When a sufficiently short optical pulse is employed to irradiate an object such as biological tissue, the photoacoustic effect results in the generation of acoustic signals within the object. After propagating out of the object, these signals can be measured by use of wide-band ultrasonic transducers. The goal of PACT in its canonical formulation is to reconstruct an image that represents a map of the initial pressure distribution within the object from knowledge of the measured photoacoustically-induced acoustic signals. The initial pressure distribution is proportional to the absorbed optical energy distribution within the object, which can reveal diagnostically useful information based on endogenous hemoglobin contrast or exogenous contrast if molecular probes are utilized. As such, PACT can be viewed either as an ultrasound mediated optical modality or an ultrasound modality that exploits optical-enhanced image contrast (?).

Over the past few decades, there have been numerous fundamental studies of photoacoustic imaging of biological tissue (?; ?; ?; ?; ?; ?; ?), and the development of PACT continues to progress at a tremendous rate (?). Biomedical applications of PACT include small animal imaging (?; ?; ?) and human breast imaging (?; ?; ?; ?; ?), to name only a few. For additional information regarding applications of PACT, the reader is referred to the many review articles that have been published on this topic (?; ?; ?; ?; ?).

PACT is a computed imaging modality that utilizes an image reconstruction algorithm for image formation. From a physical perspective, the image reconstruction problem in PAT corresponds to an acoustic inverse source problem (?). When an idealized imaging scenario is considered, a variety of analytic solutions to the PACT inverse problem are available (?; ?; ?; ?); however, in practice, numerous challenges exist that can limit their applicability. Alternatively, optimization-based, or iterative, image reconstruction methods for PACT provide the opportunity to enhance image quality by compensating for physical factors, noise, and data-incompleteness. While such approaches are routinely employed in the broader image reconstruction community, relatively few research groups have explored such modern reconstruction methods for PACT. Although they can be computationally demanding, the advent and use of modern parallel computing technologies (?) have rendered these reconstruction algorithms feasible for many PACT applications.

In this article, the PACT image reconstruction problem is reviewed within the context of modern optimization-based image reconstruction methodologies. This review is restricted to the problem of estimating the initial pressure distribution and does not address the more complicated problem of recovering the optical properties of an object (?; ?). Imaging models that relate the measured photoacoustic wavefields to the sought-after object function are described in their continuous and discrete forms. These models will describe physical non-idealities in the data such as those introduced by acoustic inhomogeneity, attenuation, and the response of the imaging system. The basic principles of optimization-based PACT image reconstruction from discrete measurement data are presented, which includes descriptions of forward operators that accurately describe the physics of image acquisition. Furthermore, the derivation of the adjoints of the corresponding forward operators are also reviewed. The adjoint operators facilitate the application of gradient-based approaches in solving the optimization-based PACT image reconstruction problem. Non-conventional formulations of the PACT image reconstruction problem, in which acoustic parameters of the medium are concurrently estimated along with the PACT image, are also introduced and reviewed.

The article is organized as follows. In Section 2, two canonical forward models employed for PACT are reviewed that are based on the acoustic wave equation and integral geometry formulations. The explicit formulation of discrete imaging models are described in Section 3. Optimization-based image reconstruction methods that are based on the discrete imaging models are presented in Section 4. Furthermore, joint reconstruction approaches to PACT image reconstruction whereby acoustic parameters of the medium are concurrently estimated along with the initial pressure distribution are discussed in Section 6. Section 7 concludes the article by providing a brief overview of the challenges and opportunities for research related to image reconstruction for practical applications of PACT.

2 Canonical forward models in their continuous forms

A schematic of a general PACT imaging experiment is shown in Fig. 1. A sufficiently short (?) laser pulse is employed to irradiate an object at time $t=0$ and the photoacoustic effect results in the generation of internal pressure distribution inside of the object, which is denoted as $p_{0}(\mathbf{r})$, $\mathbf{r}\in\mathbb{R}^{3}$. This pressure distribution subsequently propagates outwards in the form of acoustic waves and is detected by the wide-band point-like ultrasonic transducers located on a measurement aperture $\Omega_{0}\subset\mathbb{R}^{3}$. The use of alternative measurement technologies, such as integrating line or area detectors (?; ?), have also been widely explored but will not be reviewed here. The measurement aperture $\Omega_{0}$ is a two-dimensional (2D) surface that partially or completely surrounds the object. The propagated pressure wavefield at time $t>0$ will be denoted as $p(\mathbf{r},t)$.

Below, two types of 3D canonical forward models that relate the measured propagated pressure wavefield to the sought after object function $p_{0}(\mathbf{r})$ are described in their continuous forms. These models, which form the basis for mathematical studies of the inverse problem in PACT, do not model transducer characteristics. The effects of finite sampling and other physical factors associated with transducer characteristics are addressed in the discrete versions of these models presented in Section 3.

2.1 Wave-equation based PACT forward model in its continuous form

Consider an idealized PACT experiment in which an object with known heterogeneous acoustic properties is irradiated with laser pulse to photoacoustically induce an initial pressure distribution $p_{0}(\mathbf{r})$. The initial pressure distribution subsequently generates outwardly propagating broadband acoustic waves in the surrounding medium. The acoustic properties of the surrounding medium modulate the behavior of the propagating acoustic waves based on the acoustic wave equation. Let $c_{0}(\mathbf{r})$ denote the medium’s ambient speed of sound (SOS) distribution, $\rho(\mathbf{r},t)$ and $\rho_{0}(\mathbf{r})$ denote the distributions of the medium’s acoustic and ambient densities, respectively. Let $\dot{\mathbf{s}}(\mathbf{r},t)$ denote the particle velocity in the medium. The initial pressure distribution $p_{0}(\mathbf{r})$ and all quantities that describe the properties of the medium are assumed to be represented by bounded functions that have compact supports.

For a wide range of PACT applications, acoustic absorption is not negligible (?; ?; ?; ?; ?). Hence, the strong dependence of the amplitude, spectrum and shape of the propagating broadband acoustic pressure signals on the absorption characteristics of the medium needs to be modeled. It is well known that over diagnostic ultrasound frequency ranges, the acoustic absorption in biological tissue can be described by a frequency power law of the form (?; ?)

[TABLE]

where $f$ is the temporal frequency in MHz, $\alpha_{0}(\mathbf{r})$ is the spatially varying frequency-independent absorption coefficient in dB MHz*-ycm-1*, and $y$ is the power law exponent that is typically in the range of 0.9-2.0 (?). Note that classical lossy wave equations predict an absorption that is either frequency independent or proportional to frequency squared (?).

To describe the effects of power law absorption and dispersion, formulations of the wave equation that model time-domain fractional derivative operators as convolutions have been proposed (?; ?; ?; ?). Numerical implementation of time-domain fractional derivative operators for accurately modeling 3D acoustic wave propagation poses a significant memory burden (?). To overcome this memory burden, some work has focused on devising a recursive algorithm to compute the time-domain fractional derivative (?). However, this approach is heuristic and requires a priori optimization for each value of $y$. To circumvent this, a lossy wave equation modeling power law absorption using fractional Laplacian operators has been proposed (?). The fractional Laplacian operators can be easily implemented numerically to effectively model the power law absorption in 3D lossy media. Although, the proposed fractional Laplacian derivative operator-based lossy wave equation can accurately model power law absorption, it does not exhibit the correct dispersive sound speed relation (?). In order to account for the dispersion inconsistencies, a lossy wave equation that utilizes two fractional Laplacian derivative operators was proposed (?; ?). These two fractional Laplacian derivative operators account for the required power law absorption and dispersion terms, separately.

Given the frequency-dependent power law absorption model, one can formulate the PACT forward model in a lossy, acoustically heterogeneous fluid media as (?; ?; ?):

[TABLE]

Here, the spatially varying quantities $\mu(\mathbf{r})$ and $\eta(\mathbf{r})$ describe the acoustic absorption and dispersion proportionality coefficients that are defined as

[TABLE]

The second and third terms on the right hand side of Eqn. (2c) account for the required power law absorption and dispersion terms separately through two fractional Laplacian derivative operators. The initial value problem defined in Eqn. (2) describes the propagation of photoacoustically generated pressure data $p(\mathbf{r},t)$ given the spatially varying sound speed $c_{0}(\mathbf{r})$, ambient density $\rho_{0}(\mathbf{r})$, and the acoustic absorption coefficient $\alpha_{0}(\mathbf{r})$.

Consider that the measured pressure data $p(\mathbf{r},t)|_{\mathbf{r}=\mathbf{r}^{\prime}}$ are recorded outside the support of the object for $\mathbf{r}^{\prime}\in$ $\Omega_{0}$ and $t\in[0,T]$. The continuous PACT forward model consists of the composition of the partial differential equation (PDE) described in Eqn. (2) and an observation operator that restricts the pressure recorded to the measurement surface $\Omega_{0}$. Hence, the mapping from the initial pressure distribution $p_{0}(\mathbf{r})$ to the pressure recorded on the continuous measurement aperture $\Omega_{0}$ in a acoustically heterogeneous fluid media can be expressed as:

[TABLE]

where the wave equation-based forward operator $\mathcal{H}_{wave}:\mathbb{L}_{2}(\mathbb{R}^{3})\mapsto\mathbb{L}_{2}(\Omega_{0})\times[0,T]$ describes a continuous-to-continuous (C-C) mapping between two function spaces. The reason for the use of such terminology is that the wave equation-based forward operator maps a function of continuous variable (as opposed to a discrete variable) to another function of a continuous variable. There is no implication that either the function itself is continous or even that the mapping is continuous.

2.2 Integral geometry-based forward model

In the special case of a lossless and acoustically homogeneous infinite medium, the solution to the wave equation in Eqn. (2) can be expressed as (?)

[TABLE]

where $c_{0}$ denotes the constant sound speed value, $\delta(t)$ denotes the one-dimensional (1D) Dirac delta function, and the integral geometry-based C-C forward operator is denoted as $\mathcal{H}_{int}:\mathbb{L}_{2}(\mathbb{R}^{3})\mapsto\mathbb{L}_{2}(\Omega_{0})\times[0,T]$.

Equation (5) can be conveniently expressed in terms of the spherical Radon transform (SRT) (?) as

[TABLE]

where the data function $g(\mathbf{r}^{\prime},t)$ is related to the measured pressure data $p(\mathbf{r}^{\prime},t)$ by

[TABLE]

In the special case of a 2D problem, the spherical Radon transform model reduces to a circular Radon transform (?; ?).

For the case of a known weakly heterogeneous acoustic media, Eqn. (6) can be replaced by a generalized Radon transform model that is derived from the wave equation by use of a geometrical acoustics approximation. In that case, the measured pressure signals are related to $p_{0}(\mathbf{r})$ through integration over nonspherical isochronous surfaces (?; ?; ?).

Heuristic strategies have also been proposed to mitigate image artifacts that result from employing the idealized forward model in Eqn. (6) in the presence of unknown acoustic heterogeneity. For example, half-time and partial-time image reconstruction methods have been proposed for PACT image reconstruction from temporally truncated measurements that exclude components of the measured data that are strongly aberrated (?; ?; ?).

3 Discrete forward models

3.1 Review of semidiscrete and discrete forward modeling

As with any digital imaging system, the data acquired in a PACT experiment represent a finite collection of numbers that form a vector. As such, the PACT forward operator is fundamentally a continuous-to-discrete (C-D) mapping (?) that relates $p_{0}(\mathbf{r})$ to the measurement vector. The C-D operator for PACT can be expressed as

[TABLE]

where $\mathcal{D}$ is a discretization operator that spatially and temporally samples the pressure wavefield $p(\mathbf{r}^{\prime},t)$ and $\mathcal{H}_{CC}$ represents a C-C PACT forward operator such as $\mathcal{H}_{int}$ or $\mathcal{H}_{wave}$.

Let ${\mathbf{u}}\in\mathbf{R}^{QK\times 1}$ denote a lexicographically ordered representation of the sampled pressure data and let $[{\mathbf{u}}]_{m}$ denote its $m$-$th$ component. When ideal point-like ultrasound transducers are assumed, the measured samples of the pressure wavefield are given by

[TABLE]

for $k=0,1,...,K-1$. Here, $k$ and $K$ denote the temporal sample index and total number of temporal samples, respectively, and the vectors $\mathbf{r}_{q}^{\prime}\in\mathbb{R}^{3}$, $q=0,1,2,...,Q-1$, describe the locations of the $Q$ transducers on the aperture $\Omega_{0}$.

More generally, as discussed in Section 3.4, the measured pressure data can be described as (?)

[TABLE]

where $\sigma_{q}(\mathbf{r}_{0})$ and $\tau_{k}(t)$ are functions that describe the spatial and temporal sampling apertures of the transducers, respectively.

In order to employ the algebraic or optimization-based image reconstruction methods described later, the C-D forward operator must be approximated by a discrete-to-discrete (D-D) one. To accomplish this, a finite dimensional representation of $p_{0}(\mathbf{r})$ can be employed. An $N$-dimensional representation of $p_{0}(\mathbf{r})$ can be described as

[TABLE]

where $N>0$ and the superscript $a$ indicates that $p_{0}^{a}(\mathbf{r})$ is an approximation of $p_{0}(\mathbf{r})$. The functions $\phi_{n}(\mathbf{r})$ are called expansion functions and the expansion coefficients $[\boldsymbol{\theta}]_{n}$ are elements of the $N$-dimensional vector $\boldsymbol{\theta}$. The goal of image reconstruction is to estimate $\boldsymbol{\theta}$ for a fixed choice of the expansion functions $\phi_{n}(\mathbf{r})$. As described below, different choices for the $\phi_{n}(\mathbf{r})$ and rules for determining $\boldsymbol{\theta}$ will result in different D-D forward models. The specific choice of expansion functions may be motivated by various theoretical and practical reasons including a desire to minimize representation error, incorporation of a priori information regarding the object, and efficient computation. Popular choices of expansion functions in PACT include cubic or radially symmetric expansion functions known as Kaiser-Bessel (KB) window functions (?; ?; ?; ?; ?; ?), and linear interpolation functions (?; ?; ?; ?). In general, these D-D imaging models will have distinct numerical properties that will affect the performance of iterative reconstruction algorithms (?).

D-D forward models can be established by substitution of a finite-dimensional object representation into the C-D imaging model:

[TABLE]

where the D-D operator $\mathbf{H}$ is commonly referred to as the system matrix. An element in the $n$-th row and $m$-th column of $\mathbf{H}$ will be denoted by $[\mathbf{H}]_{n,m}$.

Next, examples of D-D PACT imaging models for use with homogeneous and heterogeneous acoustic media are reviewed.

3.2 D-D forward models for use with homogeneous media

In this subsection, two popular D-D PACT imaging models for homogeneous acoustic media are reviewed that employ different choices for the expansion functions.

3.2.1 Interpolation-based D-D PACT model

In the interpolation-based D-D PACT imaging model, the the associated expansion functions can be expressed as (?)

[TABLE]

where $\Delta_{s}$ is the distance between neighboring points on an uniform and isotropic Cartesian grid. The coefficient vector $\boldsymbol{\theta}$ can be defined as (?):

[TABLE]

where $\mathbf{r}_{n}=(x_{n},y_{n},z_{n})^{T}$ denotes the location of the $n$-th Cartesian grid node. The corresponding D-D PACT imaging model based on the interpolation expansion functions can be expressed as

[TABLE]

where

[TABLE]

Here, $\mathbf{G}$ and $\mathbf{D}$ are discrete approximations of the SRT operator (Eqn. (6)) and the differential operator (Eqn. (7)), respectively.

The discrete SRT operator $\mathbf{G}$ can be implemented in a “temporal-sample-driven” manner; namely, the pressure data are computed by accumulating the contributions from voxels on a discretized spherical shell surface defined by the current data sample (?):

[TABLE]

where $N_{i}$, $N_{j}$ denotes the numbers of angular divisions over the polar and azimuth directions within the local spherical coordinate system centered at the $q$-th transducer $\mathbf{r}_{q}^{\prime}$ with a radius of $k\Delta_{t}c_{0}$, in which the center of $i$-th polar division and $j$-th azimuth division is denoted by $\mathbf{r}_{k,i,j}$  (?).

The differential operator $\mathbf{D}$ can be implemented as

[TABLE]

Due to their large sizes, the explicit storage of system matrices is typically infeasible with current computer technologies except for small scale problems. By use of Eqns. (17) and (18) the action of $\mathbf{H}_{int}$ can be computed without having to explicitly store its elements.

3.2.2 Kaiser-Bessel function-based D-D PACT model

The Kaiser-Bessel (KB) function-based D-D PACT model employs the KB functions of order $m$ as the expansion functions. These KB functions are defined as (?; ?; ?)

[TABLE]

where $x\in\mathbb{R}^{+}$, $a\in\mathbb{R}^{+}$ denotes the support radius of the KB function, $\gamma\in\mathbb{R}^{+}$ describes the smoothness of the KB function. The term $I_{m}(x)$ denotes the modified Bessel function of the first kind of order $m$. The associated expansion function can be expressed as

[TABLE]

which is simply the KB function centered at location $\mathbf{r}_{n}$.

When employing the KB expansion functions, it is convenient to formulate the corresponding D-D PACT imaging model in the temporal frequency domain (?). Let $\tilde{\mathbf{u}}$ denotes the discrete Fourier transform (DFT) of the measurement data $\mathbf{u}$. The KB function-based D-D PACT imaging model can be expressed as (?)

[TABLE]

where the elements of the system matrix $\tilde{\mathbf{H}}_{KB}$ can be written as the multiplication of two terms:

[TABLE]

where

[TABLE]

and

[TABLE]

Here, $L$ denotes the total number of frequency components and $\Delta_{f}$ is the frequency sampling interval. The first quantity $\tilde{p}^{KB}(f)$ in Eqn. (22) is the temporal Fourier transform of the acoustic pressure generated by a KB function located at origin, where $\hat{j}_{m}(x)$ is the $m$-th order spherical Bessel function of the first kind (?). The second quantity $\tilde{h}_{q}^{s,point}(\mathbf{r}_{n},f)$ is the spatial impulse response (SIR) in the temporal frequency domain. Under the assumption of an idealized point-transducer model, this quantity reduces to the Green function given in Eqn. (23), where $\mathbf{r}_{q}^{s}$ denotes the location of the $q$-th transducer and $\mathbf{r}_{n}$ denotes the location of the center of the $n$-th KB expansion function (?).

3.3 D-D forward models for use with heterogeneous acoustic media

3.3.1 Overview of approaches

The establishment of D-D imaging models for lossy, heterogeneous fluid media requires the introduction of finite-dimensional approximations of the object function $p_{0}(\mathbf{r})$ as well as the known medium acoustic parameters such as $c(\mathbf{r}),\rho_{0}(\mathbf{r})$ and $\alpha(\mathbf{r})$. In principle, these medium acoustic parameters can be estimated from adjunct ultrasound tomography data (?; ?). When such adjunct image data are not available, one can attempt to estimate the acoustic parameters concurrently with $p_{0}(\mathbf{r})$ as described in Section 6. For PACT in heterogeneous media, the choice of finite-dimensional representation of the object and medium parameters is dictated by the numerical method employed to solve the coupled first order differential equations described in Eqn. (2). Most methods for computing a numerical solution of the acoustic wave equation in lossy heterogeneous media fall into one of three major categories: finite difference (FD) methods, finite element (FE) methods, and spectral methods.

The FD and FE methods are also known as local methods because the wave propagation equations are solved at each point based only on conditions at neighboring points. In contrast, spectral methods are global, as information from the entire wavefield is employed to numerically propagate the wavefield. As the spectral methods leverage global information, they allow computation to be performed on coarser grid while maintaining high accuracy (?). As opposed to FD or FE methods that require grid spacing on the order 10 points per minimum wavelength, spectral methods, in theory, only require 2 points per minimum wavelength (?). Due to such relaxed spatial sampling constraints, spectral methods are well suited for large-scale, high frequency PA wave solvers. One of the most popular spectral methods used for solving the acoustic wave equation is the k-space pseudospectral method (?; ?; ?). Because of its popularity in the PACT community, the D-D imaging model presented below will be based upon the k-space pseudospectral time-domain (PSTD) method (?).

3.3.2 D-D forward model based on the k-space PSTD method

Below, a general formulation of a D-D forward model based on the k-space PSTD method is briefly outlined. Because of the highly technical nature of the formulation, the reader is referred to the literature for specific details (?).

The finite-dimensional approximations of the object function $p_{0}(\mathbf{r})$ as well as the medium parameters $c(\mathbf{r}),\rho_{0}(\mathbf{r})$ and $\alpha(\mathbf{r})$ are constructed by sampling these quantities on a Cartesian grid. If $\{\phi_{n}^{\rm int}(\mathbf{r})\}_{n=0}^{N-1}$ denotes the set of 3D Cartesian expansion functions defined in Eqn. (13), the initial pressure distribution $p_{0}(\mathbf{r})$ can be approximated as a finite dimensional vector $\boldsymbol{\theta}\in\mathbb{R}^{N\times 1}$. To ensure the stability of the k-space PSTD method, the material properties such as the speed of sound $c(\mathbf{r})$, the ambient density $\rho_{0}(\mathbf{r})$, the absorption coefficient $\alpha(\mathbf{r})$, and the wavefield parameters such as pressure $p(\mathbf{r},t)$, particle velocity $\dot{\mathbf{s}}(\mathbf{r},t)$, and acoustic density $\rho(\mathbf{r},t)$, are sampled at different points of staggered Cartesian grids (?). Furthermore, the wavefield parameters such as the pressure, particle velocity and acoustic density are staggered temporally. Thus, the finite-dimensional representations of the medium parameters and the wavefield parameters for use in a D-D imaging model will generally employ different expansion functions.

Let $\dot{\mathbf{s}}^{i}_{k}\in\mathbb{R}^{N\times 1}$, $\boldsymbol{\rho}^{i}_{k}\in\mathbb{R}^{N\times 1}$ represent the discrete approximations of the vector-valued particle velocity and acoustic density over the whole 3D Cartesian grid at the $k^{th}$ time step along the $i^{th}$ direction. Let $\mathbf{p}_{k}\in\mathbb{R}^{N\times 1}$ represent the acoustic pressure sampled at all spatial grid points at the $k^{th}$ time step.

Let $\mathbf{v}_{k}=(\dot{\mathbf{s}}_{k}^{1},\dot{\mathbf{s}}_{k}^{2},\dot{\mathbf{s}}_{k}^{3},\boldsymbol{\rho}_{k}^{1},\boldsymbol{\rho}_{k}^{2},\boldsymbol{\rho}_{k}^{3},\mathbf{p}_{k})^{T}$ denote a $7N\times 1$ vector containing all the wavefield variables at the time step $k\Delta t$. The image acquisition process in PACT can be mathematically modeled as the propagation of the wavefield parameters forward in time from $t=0$ to $t=(K-1)\Delta t$ as

[TABLE]

For specific details regarding the definition of the temporal matrix $\mathbf{T}_{k}(k=1,\cdots,K-1)\in\mathbb{R}^{7NK\times 7NK}$, the reader is referred to the literature (?). From the initial conditions defined in Eqn. (2), the vector $(\mathbf{v}_{0},\mathbf{0}_{7N\times 1},\cdots,\mathbf{0}_{7N\times 1})^{T}$ can be computed from the initial pressure distribution $\boldsymbol{\theta}$ as

[TABLE]

where $\mathbf{T}_{0}\in\mathbb{R}^{7N\times N}$ maps the initial pressure distribution to the initial value of the wavefield variables at time $t=0$ (?).

In general, the transducer locations $\mathbf{r}_{q}^{\prime}$ at which the pressure are recorded will not coincide with the vertices of the 3D Cartesian grid at which the propagated wavefields are computed. The pressure at the transducer locations can be related to the computed field quantities via an interpolation operation defined as

[TABLE]

where $\mathbf{M}\in\mathbb{R}^{QK\times 7NM}$. The choice of the interpolation operation is a parameter that will guide the construction of the matrix $\mathbf{M}$. Some of the most commonly employed interpolation strategies include trilinear interpolation or Delaunay triangulation-based interpolation (?).

By use of Eqns. (25), (24) and (26), one obtains

[TABLE]

The explicit form of the system matrix $\mathbf{H}_{PS}\in\mathbb{R}^{QK\times N}$, that describes the propagation of PA waves based on Eqn. (2) can be therefore expressed as

[TABLE]

Here, the subscript $PS$ stands for the fact that the D-D imaging model is computed using the k-space PSTD method.

3.4 Incorporation of transducer responses in D-D forward models

In practice, the ultrasound transducers employed in a PACT imager are imperfect and result in measurements of the PA wavefield that are averaged over finite temporal and spatial apertures as described in Eqn. (10). In the ultrasound imaging community, the effects of these sampling apertures are characterized by the transducer’s electrical impulse response (EIR) and spatial impulse response (SIR). Failure to account for these effects can result in reconstructed estimates of $p_{0}(\mathbf{r})$ that contain distortions and/or degraded spatial resolution (?).

When iterative image reconstruction methods are to-be-employed, a natural way to compensate for the EIR and SIR is to include them in the constructed D-D forward model (?; ?; ?). Below, the D-D forward models for use with homogeneous acoustic media introduced in Section 3.2 are extended to incorporate the transducer responses. The extension of the the D-D forward model for use with heterogeneous acoustic media introduced in Section 3.3 to include transducer responses is relatively straightforward and is described in the literature (?).

3.4.1 Incorporating the EIR in D-D forward models for use with homogeneous acoustic media

Because degradation of the measured data by the EIR is described by a linear time-invariant model, it can be readily incorporated into the system matrices. For example, the interpolation-based system matrix with EIR can be re-defined as (?)

[TABLE]

where $\mathbf{D},\mathbf{G}$ are the discrete approximation so the differential operator and the SRT operator and

[TABLE]

Here, $[\mathbf{h}^{e}]_{k}=\Delta_{t}h^{e}(t)|_{t=k\Delta t}$ is the EIR signal sampled in the temporal domain, $\mathcal{F}$ and $\mathcal{F}^{-1}$ represent the discrete Fourier transform and inverse discrete Fourier transform, respectively, and $\mathbf{p}_{ideal}$ denotes the temporally sampled ideal pressure signal vector that would be recorded by an idealized point transducer.

A similar approach can be adopted to incorporate EIR in the KB function-based system matrix. In this case, the convolution operation can be implemented as an element-wise multiplication in the temporal frequency domain (?; ?) and the elements of the system matrix are re-defined as

[TABLE]

where samples of $\tilde{h}^{e}(f)$ are computed as the discrete Fourier transform of $\mathbf{h}^{e}$.

3.4.2 Incorporating the SIR in D-D forward models for use with homogeneous acoustic media

The spatial impulse response (SIR) describes the spatial averaging of an acoustic signal that occurs when the signal is measured by use of a transducer possessing a non-zero active area. Specifically, consider a point acoustic source located at position $\mathbf{r}_{n}$ whose temporal response is described by $\delta(t)$. The SIR ${h}_{q}^{s}(\mathbf{r}_{n},t)$ represents the measurement of the radiated wavefield by a transducer indexed by $q$ that has an idealized EIR.

Various SIR models have been proposed (?; ?; ?) and employed in PACT (?; ?; ?; ?; ?; ?; ?; ?; ?). Because the degradation of the measured signal by the SIR is not generally described by a linear time-invariant model, it is not as straightforward to compensate for as is the EIR.

For the interpolation-based D-D imaging model in Section 3.2, Ding et al. proposed to approximately incorporate the SIR by summing the pressure signal over a collection of points on the transducer surface. However, in order to accurately model the effects of the SIR, it may be necessary to sample a large number of points on the transducer surface, which can result in a large computational burden. It remains generally difficult to accurately incorporate the SIR in an interpolation-based PACT D-D forward model and therefore certain implementations of this model have ignored the SIR (?).

In the KB function-based system matrix, the SIR can be readily incorporated as an additional element-wise multiplication step in the temporal frequency domain:

[TABLE]

Here, $\tilde{h}_{q}^{s}(\mathbf{r}_{n},f)$ denotes the temporal Fourier transform of ${h}_{q}^{s}(\mathbf{r}_{n},t)$ that can be computed by integrating the Green function in Eqn. (23) over the $q$-th transducer surface $S_{q}$:

[TABLE]

Note that the integral in Eqn. (33) resembles the Rayleigh integral (?). As a specific example, consider a transducer element that possesses a rectangular detecting surface of area $a\times b$. Under the validity of the far field assumption (?; ?):

[TABLE]

Eqn. (33) can be evaluated to determine $\tilde{h}_{q}^{s}(\mathbf{r}_{n},f)$ as (?)

[TABLE]

where $x_{nq}^{tr},y_{nq}^{tr}$ are the transverse components of $\mathbf{r}_{n}$ in a local coordinate system centered at $\mathbf{r}_{q}^{s}$. Given $\mathbf{r}_{n}=(x_{n},y_{n},z_{n})$ and the transducer location $\mathbf{r}^{s}_{q}$ specified in spherical polar coordinates $\mathbf{r}^{s}_{q}=(r_{q},\theta_{q},\phi_{q})$, the values of the transverse coordinates can be computed as (?):

[TABLE]

3.4.3 Patch-based estimation of the SIR

When the far field condition in Eqn. (34) is violated, the SIR expression in Eqn. (35) can be inaccurate and its use may lead to artifacts in the reconstructed images (?). Moreover, when the transducer surface is not flat, it may be difficult to determine an analytic expression for the SIR. These issues can be mitigated by use of “divide-and-integrate” approaches, including line-detector-based method (?) and patch-based method (?). In such approaches, the surface of the transducer is computationally decomposed into smaller elements whose SIRs can be more readily determined. In the line-detector based model, each transducer element surface is decomposed into a number of parallel straight lines, whose SIRs can be analytically expressed. In the patch-based method, each transducer element surface is divided into smaller planar patches that each satisfy the far field approximation (?); subsequently, in either case, the SIR of the transducer can be approximated by computing the average of the SIRs of the sub-elements (lines or patches):

[TABLE]

where $\tilde{h}_{q,i}^{s}(\mathbf{r}_{n},f)$ denotes the SIR corresponding to $i$-th sub-element. In the patch-based approach, $\tilde{h}_{q,i}^{s}(\mathbf{r}_{n},f)$ can be computed with the aid of Eqn. (35). It should be noted that, in the line-detector model, the approximation only accounts for the SIR effect in the direction that is parallel to the straight lines, it still assumes zero thickness (point-like) transducer in the perpendicular direction. The patch-based model employs the far-field approximation for both directions in the transducer plane, therefore providing compensation in both directions. In addition, the patch-based model can be extended beyond planar transducers by decomposing an arbitrary transducer surface into smaller patches to estimate its SIR.

A computer simulation study examining the effects of accurately modeling the SIR for flat rectangular transducers in the context of PACT image reconstruction is discussed below (?). The numerical phantom used for the study consisted of spherical objects placed within a PACT system as shown in Fig. 2(a). The simulated pressure data for each sphere were generated by numerically convolving a closed form expression of waves generated by a solid sphere (?), with a semi-analytical SIR specifically for spherical waves (?). During reconstruction, the Kaiser-Bessel function-based D-D PACT imaging model introduced in Section 3.2.2 was employed with different SIR models: Fig. 2(b) assumes a point-like transducer model, Fig. 2(c) employs the far-field-based SIR model as in Eqns. (33) and (35), and Fig. 2(d) employs the patch-based model in Eqn. (37). The reconstructed images where obtained by solving a penalized-least-squares optimization problem with quadratic smoothness penalty using the conjugate gradient method. These results show that ignoring transducer SIR in the PACT imaging model leads to significant distortion in reconstructed images. Incorporating an accurate transducer SIR model can greatly reduce such artifacts. For objects far away from the transducers, the far-field-based SIR model may suffice. However, for closer objects, patch-based or other ”divide-and-integrate” SIR models can further improve the reconstruction accuracy.

4 Image reconstruction approaches

In this section, a brief survey of conventional image reconstruction methods for PACT is provided. Subsequently, optimization-based image reconstruction methods that are based upon D-D forward models are presented in Section 5.

4.1 Brief overview of analytic reconstruction methods

A large number of analytic, or non-iterative, tomographic reconstruction methods for PACT are currently available. Several review articles provide detailed descriptions of these (?; ?), so only a quick overview is provided here. Analytic reconstruction methods are commonly based upon a C-C PACT forward model that assumes a homogeneous acoustic media, such as those described in Section 2.2. Analytic reconstruction formulae of filtered backprojection (FBP) form are available for special detector geometries (?; ?; ?; ?; ?; ?). Such reconstruction methods typically assume an acoustically homogeneous and lossless medium and full-view acoustic detection in which the pressure measurements are densely sampled on a closed surface that encloses the object. Additionally, the available analytic methods are unable to compensate for the ultrasound transducer responses (?). Other ad hoc reconstruction techniques such as inversion of the Radon transform (?), as well as delay-and-sum techniques (?) have been employed for PACT image reconstruction.

PACT reconstruction methods that are based on harmonic decomposition have also been proposed (?). For special geometries such as planar detector geometries, such reconstruction methods can be implemented efficiently by use of the fast Fourier transform (FFT) (?; ?). This approach has been extended for use with novel measurement geometries that utilize reverberant cavities (?), and for use with closed measurement surfaces for which the eigenfunction of the Dirichlet Laplacian are explicitly known (?). These approaches possess computationally efficient implementations and, in certain cases, can reconstruct 3D images thousands of times faster than backprojection-type methods (?). Similar to the FBP-type algorithms, the harmonic decomposition-based reconstruction methods typically assume an acoustically homogeneous media and do not compensate for the transducer responses.

4.2 Time-reversal reconstruction methods

As an alternative to the analytic reconstruction approaches discussed above, a variety of time-reversal (TR) based reconstruction methods have been proposed (?; ?; ?; ?; ?). Image reconstruction via TR is achieved by running a numerical model of the wave equation-based forward problem backwards in time. Namely, the measured acoustic pressure signals are re-transmitted into the domain in a time-reversed fashion. As such, they can readily account for wave propagation in heterogeneous media and can accommodate arbitrary detection geometries.

4.2.1 Formulation of TR image reconstruction for PACT

Consider a compactly supported initial pressure distribution $p_{0}(\mathbf{r})$ in an infinite 3D homogeneous lossless acoustic medium, and let the domain $\Omega$ correspond to the interior of a closed measurement surface $\Omega_{0}$ that encloses the to-be-imaged object. According to Huygen’s principle, there exists a time $T>0$ for which the radiated wavefield inside $\Omega$ vanishes $\forall t\geq T$. Because $\Omega$ contains no energy $\forall t\geq T$, one can solve the wave equation backwards in time inside $\Omega$ with zero initial conditions and the boundary condition given by the measured data on the surface $\Omega_{0}$. As described mathematically below, this process of rewinding the waves backward in time to reconstruct the $p_{0}(\mathbf{r})$ is defined as TR. When finite sampling effects are ignored, it has been demonstrated that time reversal yields a theoretically exact reconstruction of $p_{0}(\mathbf{r})$ for the case of a 3D acoustically homogeneous medium if $T$ is large enough to allow the energy to escape the domain $\Omega$ (?). In even dimensions or when the medium is acoustically heterogeneous, this result does not hold true. Nevertheless, the TR method can still be employed to obtain an estimate of $p_{0}(\mathbf{r})$. When the speed of sound distribution is non-trapping (i.e. all of the acoustic energy escapes the domain $\Omega$), the errors in the estimates can be bounded as a function of the acquisition time $T$ (?).

In addition, the same principle of replaying the wavefield back in time can also be applied in lossy fluid media. Intuitively, to compensate for the amplitude loss in the wavefields, the wavefields should be corrected by gain factors. To account for acoustic absorption, the absorption term in Eqn. (2c) must be reversed in sign when computing the time-reversed wavefields. In contrast, to account for the dispersion, the sign of the dispersive speed is left unchanged (?). Thus, the TR reconstruction method for use with lossy heterogeneous media solves the following system of equations for the time-reversed pressure wavefield $p(\mathbf{r},t)$  (?):

[TABLE]

Here, $u(\mathbf{r},t)$ denotes the measured pressure data for the case where finite sampling effects are neglected and idealized point-like ultrasound transducer are employed. The sought after estimate of $p_{0}(\mathbf{r})$ is specified by the TR method as

[TABLE]

In operator form, Eqns. (38) and (39) can be expressed as

[TABLE]

where $\mathcal{H}^{TR}:\mathbb{L}_{2}(\Omega_{0}\times[0,T])\mapsto\mathbb{L}_{2}(\Omega)$ denotes the time reversal operator described in Eqn. (38).

4.2.2 Modified TR image reconstruction based on a Neumann series method

As described above, the canonical TR method is mathematically exact when the radiated wavefield $p(\mathbf{r},t)$ decays to zero inside the domain $\Omega$ for some finite time $T$. However, in even dimensions and/or when the medium is acoustically heterogeneous and non-trapping, this condition may not be satisfied (?; ?). A reconstruction method based on a Neumann series (NS) expansion has been developed that can be employed for accurate image reconstruction in the presence of acoustically heterogeneous media and an irregular observation geometry for a finite acquisition time $T$ (?; ?). The NS-based reconstruction method is accurate when the pressure $p(\mathbf{r},t)$ is known on the whole boundary $\Omega_{0}$ and $T$ is greater than a stability threshold (?). To define the NS-based reconstruction method, one needs to first introduce the modified TR operator. The modified TR operator for lossy heterogeneous media solves Eqns. (38a), (38b) and (38c) subject to the initial and boundary value conditions:

[TABLE]

where $P_{\Omega}:\mathbb{L}_{2}(\Omega_{0})\mapsto\mathbb{L}_{2}(\Omega)$ is a harmonic extension operator. The harmonic extension operator is defined as $P_{\Omega}g(\mathbf{r},T)=\phi(\mathbf{r})$, where $\phi(\mathbf{r})$ is the solution to the elliptic boundary value problem

[TABLE]

One can summarize Eqns. (41),  (39) and (42) as

[TABLE]

where $\mathcal{H}_{mod}^{TR}:\mathbb{L}_{2}(\Omega_{0}\times[0,T])\mapsto\mathbb{L}_{2}(\Omega)$ is the modified TR reconstruction operator specified by Eqns. (41) and (42).

Given the modified TR operator, the NS-based reconstructed initial pressure distribution can be defined as

[TABLE]

where $\mathcal{H}_{wave}$ is the C-C wave equation-based forward operator in Eqn. (4) and $I$ is the identity operator. Although the NS-based reconstruction is guaranteed to be convergent when the speed of sound distribution is non-trapping, it has been successfully applied to the non-trapping case as well (?).

As the sum defined in Eqn. (44) extends from $m=0$ to $m=\infty$, it cannot exactly be implemented. It is interesting to note that the NS method can be interpreted as an iterative time reversal method (?). Let us define the $i^{th}$ estimate of the reconstructed initial pressure distribution as

[TABLE]

The partial sum in Eqn. (45) can be expressed as

[TABLE]

From Eqn. (46), the NS partial sum can be interpreted as an iterative update step, where the $i^{th}$ iteration refers to the $i^{th}$ partial sum.

5 Optimization-based image reconstruction

A general approach to PACT image reconstruction is to formulate the sought-after estimate of $p_{0}(\mathbf{r})$ as the solution of an optimization problem. In fact, most modern image reconstruction methods for computed imaging modalities including X-ray computed tomography and magnetic resonance imaging are formulated in this way (?; ?; ?). There are potential practical and conceptual advantages to optimization-based image reconstruction over analytic methods. For example, because they are based on D-D forward models, optimization-based image reconstruction methods can comprehensively compensate for the imaging physics and other physical factors such as responses of the measurement system that are not easily incorporated into an analytic method. Moreover, such methods provide a general framework for incorporating regularization, which can mitigate the effects of measurement noise and data incompleteness. Because optimization-based methods are often implemented by use of iterative algorithms, they are generally more computationally demanding than analytic methods; however the use of modern parallel computing technologies (?) can render iterative three dimensional model-based reconstruction scheme for arbitrary photoacoustic acquisition geometries feasible for many PACT applications.

Consider a D-D imaging model $\mathbf{u}=\mathbf{H}\boldsymbol{\theta}$ as described in Section 3. An optimization-based image reconstruction method seeks to determine an estimate of $\boldsymbol{\theta}$ by solving an optimization program that can be generally specified as

[TABLE]

Here, $F(\cdot)$ is the objective function to be minimized that depends on the D-D forward operator $\mathbf{H}$, the measured data $\mathbf{u}$, and the vector $\boldsymbol{\theta}$ that specifies the estimate of $p_{0}(\mathbf{r})$. The functions $\{f_{i}(\boldsymbol{\theta})\}_{i=0}^{N_{c}}$ and the closed set $\{C_{i}\}_{i=0}^{N_{c}}$ describe the set of $N_{c}$ constraints that the solution must satisfy. It is important to note that numerous choices must be made in the design of the reconstruction method. These include the specification of $\mathbf{H}$, $F(\cdot)$, and $f(\boldsymbol{\theta})$, as well as the optimization algorithm that is employed to solve Eqn. (47). It is the joint specification of these quantities that defines the sought after solution and the numerical properties of the reconstruction method (?). It should also be noted that many generic iterative methods that have been employed for PACT image reconstruction (?; ?; ?) can be interpreted as computing solutions to specific cases of Eqn. (47). Below, some commonly employed optimization programs employed by image reconstruction methods are reviewed.

5.1 Penalized least squares methods

A special case of Eqn. (47) corresponds to a penalized least squared (PLS) estimator, and is given by

[TABLE]

Since the initial pressure distribution is non-negative, the above problem can be constrained so the solution satisfies $\boldsymbol{\theta}\geq 0$. The objective function in the PLS estimator is expressed as two terms. The quantity $\frac{1}{2}||\mathbf{u}-\mathbf{H}\boldsymbol{\theta}||_{\mathbf{W}}^{2}$ is referred to as the data fidelity term that corresponds to a weighted least squares functional. The matrix $\mathbf{W}\in\mathbb{R}^{QK\times QK}$ that defines the weighted $l_{2}$ norm is symmetric and positive definite, whereas $\mathbf{u}\in\mathbb{R}^{QK\times 1}$ corresponds to the measurement vector. This data fidelity functional is convex and differentiable with respect to $\boldsymbol{\theta}$. For the case of Gaussian measurement noise, this data fidelity term can be interpreted as a negative log-likelihood function; however, its use is not restricted to that case.

The quantity $R(\boldsymbol{\theta})$ is a penalty term that can be designed to regularize the inverse problem, and $\gamma\in\mathbb{R}$ is a regularization parameter that controls the amount of regularization. A classic form of regularization is Tikhonov regularization (?), where the regularization term is specified as

[TABLE]

where $\mathbf{P}\in\mathbb{R}^{N\times N}$ is a symmetric positive definite matrix. As the Tikhonov regularization term is differentiable with respect to $\boldsymbol{\theta}$, a variety of gradient-based methods can be employed to solve Eqn. (48).

Gradient-based methods that can be employed to solve Eqn. (48) can be broadly grouped into two classes: first order methods and second order methods. As part of first order methods, the Taylor approximation used to compute the descent direction involves a linear or first order approximation of the cost function. Hence, information about the gradient of the function is employed to compute the descent direction. Some of the most commonly employed first order methods include the steepest descent method, and the conjugate gradient method. There also exists a class of first order methods that employ information about the past gradients/momentum to speed up the convergence rate of first order algorithms. Such methods are referred as Nesterov methods, whereby linear combinations of present and past gradients are used to compute the descent direction (?; ?). In addition, a family of methods called Krylov-based methods, the subset of which is are the conjugate gradient methods, are also employed to solve smooth, convex programs defined in Eqn. (48) (?; ?; ?). The design and study of variants of the Nesterov method and the Krylov-based methods to solve smooth, convex optimization programs is an active research area (?; ?; ?).

The second class of methods regularly employed to optimization programs are the second order methods. The Taylor approximation used to compute the descent direction for second order methods is a quadratic or second order approximation of the cost function. Hence, information about the inverse Hessian of the cost function is employed to compute the descent direction. Methods that explicitly compute the inverse Hessian are referred to as Newton methods. Although the newton methods have favorable converge properties, the computational burden associated with computing a Hessian and inverting it is prohibitively large. Hence, a variety of methods that approximate the inverse Hessian from first order gradient information are commonly employed. Such methods are referred to as Quasi-Newton methods. A variety of Quasi-Newton methods have been proposed and comprehensively studied (?; ?; ?; ?; ?; ?).

The set of methods discusssed above can handle cost functions that have smooth regularization and data fidelity terms that are differentiable. However, certain modern approaches to regularization employ non-smooth choices for $R(\boldsymbol{\theta})$ based on the $l_{1}$-norm:

[TABLE]

where $\boldsymbol{\Phi}\in\mathbb{R}^{N\times N}$ is a sparsifying transformation that is chosen such that $\boldsymbol{\theta}$ is sparse in its range. Popular choices for $\boldsymbol{\Phi}$ include the discrete wavelet transform or finite difference operators. Total variation (TV) regularization (?) can be achieved as a special case when the latter are employed. Such choices are based on the observation that many objects possess a sparse representation in the wavelet or gradient domain (?; ?). Although Eqn. (50) is convex, it is not differentiable with respect to $\boldsymbol{\theta}$. However, a variety of non-smooth optimization methods, such as proximal-gradient methods, can be employed to solve Eqn. (48) in this case. Proximal methods are a type of forward-backward splitting approach, which alternates between computing a gradient step on the data fidelity term and a solution to a proximal problem involving the regularization (?; ?; ?). Moreover, the proximal gradient methods can also be accelerated by utilizing momentum information from past gradients or by deploying second order information to solve the proximal problem (?; ?; ?; ?; ?; ?; ?).

Examples of PACT images of a live mouse reconstructed by use of two different reconstruction methods are displayed in Fig. 3. A description of the imaging system and experimental details have been described in previous works (?; ?). The small animal imager consisted of an arc-shaped probe of 64 transducers along the vertical axis. Thus, a tomographic view consisted of data recorded at 64 transducers along the arc. The complete tomographic dataset was acquired by rotating the object $360^{\circ}$ around the vertical axis. A 180-view dataset was used to reconstruct 3D images of the mouse trunk as shown in Fig. 3. The image shown in Fig. 3(a) was reconstructed by use of a FBP algorithm, while the image in Fig. 3(b) was reconstructed by use of the PLS estimator in Eqn. (48) that employed a TV penalty. The image corresponding to the PLS-TV estimate possesses a higher contrast than the images reconstructed by use of the FBP algorithm and reveal a much sharper body vascular tree. These observations are consistent with the fact that optimization-based image reconstruction methods can produce images that possess different physical and statistical characteristics than the ones produced by analytical method.

An important consideration in the formulation of an optimization-based image reconstruction method is the choice of the penalty function. To demonstrate how different choices can influence image quality, Fig. 4 displays images of a simple experimental phantom reconstructed by use of Eqn. (48) with two different choices for the penalty function. Figures 4 (a) and (b) correspond to use of quadratic smoothness Laplacian regularization and TV regularization, respectively. The quadratic Laplacian penalty corresponds to the $l_{2}$-norm of a discrete Laplacian operator acting on the image estimate. In this example, use of the TV penalty resulted in images with reduced artifact levels and enhanced spatial resolution as compared to use of the quadratic Laplacian penalty.

All the images shown in Figures 3 and 4 were reconstructed by use of an iterative image reconstruction model that assumed homogeneous acoustic media and thus could not account for the acoustic heterogeneity of the medium. However, knowledge about the spatial distribution of acoustic heterogeneity of the medium can be incorporated into an iterative reconstruction algorithm based on the imaging model described in Eqn. (28). In vivo experimental studies have been conducted to validate the aforementioned iterative algorithm whereby the spatial distribution of the acoustic properties of the medium are accounted for in the reconstruction algorithm. For the in vivo studies, the experimental data were acquired from a mouse trunk using a small animal imaging system that has been described in detail in earlier works (?). As opposed to the previous studies where the medium was assumed to be homogeneous, in this study the acoustic variation between the background and the bulk mouse tissue were accounted for in the reconstruction algorithm. A reconstructed initial pressure distribution image for fixed constant sound speed value is shown in Fig. 5(a). In the reconstructed image, strong surface and interior vessel structures are observed. While the interior vessel structures appear out of focus, the surface vessels appear to be in focus. When assuming a homogeneous acoustic media, the image reconstruction algorithm could not produce images where both the surface and the interior vessels could be concurrently focused with a single tuned speed of sound value. To account for this, an image reconstruction algorithm that assigned different speed of sound values to the bulk tissue and the background was employed. The PLS-TV algorithm that was based on the imaging model defined in Eqn. (28) was employed to reconstruct the initial pressure distribution. The reconstructed initial pressure distribution when a heterogeneous speed of sound distribution is used is shown in Fig. 5(b). Comparing Figures 5(a) and 5(b), the surface and interior vessels are observed to be focused concurrently in the latter image, while only the surface vessels appear in focus in the former image.

5.2 Computation of adjoint operators and objective function gradients in PACT

The reconstruction approaches described above are general and applicable to a wide range of computational inverse problems. Here, PACT-specific details regarding the implementation of optimization-based image reconstruction methods are reviewed. Specifically, methods for computing the adjoints of some of the D-D forward operators of Section 3 are presented. The explicit computation of the adjoint of the D-D forward operator facilitates the computation of gradients of data fidelity functionals.

A necessary step in solving optimization problems of the form given in Eqn. (48) is the computation of the gradient of the data fidelity term

[TABLE]

with respect to $\boldsymbol{\theta}$. Formally, this quantity is given by

[TABLE]

where $\mathbf{H}^{\dagger}$ denotes the adjoint of the system matrix $\mathbf{H}$. As such, a key step in computing the data fidelity gradient is computing the action of the adjoint operator. As mentioned previously, for many problems of interest, $\mathbf{H}$ is too large to store in memory and its action is commonly computed on-the-fly by use of an algorithm. The same is true for $\mathbf{H}^{\dagger}$. Although the weighted $l_{2}$-norm is one of the most commonly employed data fidelity terms, a variety of convex data fidelity terms can be employed for PACT image reconstruction. Some, of the alternative convex data fidelity terms include weighted $l_{1}$-norm, KL-divergence, etc. Furthermore, the weight matrix associated with the weighted $l_{2}$-norm is a design parameter that can be constructed to incorporate a priori information about the uncertainties associated with the imaging model. In cases where an arbitrary data fidelity term is employed, the gradient of the data fidelity term can be computed efficiently through use of the adjoint state method (?; ?). In the adjoint state method, the gradient of the data fidelity term with respect to the model parameters are computed through the use of adjoint state variables. For PACT applications, the adjoint state variables are solutions to the adjoint of the wave equation defined in Eqn. (2). Hence, the computational complexity associated with the computation of the gradient of an arbitrary data fidelity term through the use of the adjoint state method would at least be on the same order as computing the action of the discrete adjoint operator $\mathbf{H}^{\dagger}$.

Although not discussed below, it should be noted that another option for establishing $\mathbf{H}^{\dagger}$ is to employ an adjoint-then-discretize approach (?). In such an approach, the explicit analytical form of the adjoint of the C-C forward operator is established. Subsequently, the C-C adjoint operator is discretized to obtain an estimate $\mathbf{\hat{H}}^{\dagger}$ of the D-D adjoint operator $\mathbf{H}^{\dagger}$. From an implementation perspective, such approaches can sometimes be more convenient than computing the adjoint operator corresponding to the assumed D-D forward operator $\mathbf{H}$. However, in cases where $\mathbf{\hat{H}}^{\dagger}$ does not accurately mimic $\mathbf{H}^{\dagger}$, the behavior of iterative algorithms can be altered (?). The spectral properties of the forward-adjoint operator pair that govern convergence properties and the restrictions associated in choosing such operator pair have been studied (?; ?).

5.2.1 Adjoint for interpolation-based forward model for homogeneous medium

Here, the interpolation-based D-D forward model described in Secs. 3.2.1 and 3.4.1 is considered. The matched adjoint operator corresponding to $\mathbf{H}_{int}$ in Eqn. (29) is defined as $\mathbf{H}_{int}^{\dagger}=\mathbf{G}^{\dagger}\mathbf{D}^{\dagger}\mathbf{H}^{e\dagger}$, where $(\cdot)^{\dagger}$ denotes the transpose of the corresponding D-D operator (i.e., a matrix). These operators can be computed as (?)

[TABLE]

While the forward operator $\mathbf{H}_{int}$ can be efficiently implemented by use of parallel computing techniques (?), the adjoint operator $\mathbf{H}_{int}^{\dagger}$ is difficult to implement efficiently. This is partly due to the fact that when implementing $\mathbf{H}^{\dagger}$ in the interpolation-based model, it is difficult to satisfy the principle of “partition on target results rather than sources” (?) for safe and efficient GPU implementation. In addition, evaluating the backward operator $\mathbf{G}^{\dagger}$ relies on expensive atomic operations on the GPU (?), resulting in a 6-10 times longer runtime for $\mathbf{H}_{int}^{\dagger}$ than for $\mathbf{H}_{int}$.

To address this problem, an approximation of the adjoint operator can be employed that closely approximates the true adjoint but is computationally more efficient to compute (?). Such operators, denoted as $\mathbf{H}^{\dagger}_{unmatched}$, are referred to as unmatched adjoint operators. An unmatched adjoint operator that approximates $\mathbf{H}_{int}^{\dagger}$ by use of a simplified voxel-driven model can be defined as:

[TABLE]

where

[TABLE]

Here, $\tilde{k}\equiv(\frac{\mathbf{r}^{s}_{q}-\mathbf{r}_{n}}{c_{0}})/\Delta_{t}$ and the value of $[\tilde{\mathbf{g}}]_{\tilde{k}}$ is approximated by linearly interpolating from its two neighboring samples:

[TABLE]

with $k$ denoting the integer part of $\tilde{k}$. The operator $\mathbf{H}_{unmatched}^{\dagger}$ approximates the operation of $\mathbf{H}_{int}^{\dagger}$, but is much faster to compute. It can be seen from Eqns. (55) and (57) that the computation of the exact adjoint $\mathbf{G}^{\dagger}$ involves $QKN_{i}N_{j}\times\mathcal{O}(1)$ calculations, while the approximated adjoint $\mathbf{G}_{unmatched}^{\dagger}$ involves only $Q\times\mathcal{O}(1)$ calculations. Though it is difficult to express the computational complexity of the matched adjoint operator in strict big-O notation, because $N_{i}$ and $N_{j}$ change with $q$ and $k$, the product of these two terms can be approximately on the order of 10,000 for 3D OAT studies; in such cases, $QKN_{i}N_{j}\times\mathcal{O}(1)$ is significantly larger than $Q\times\mathcal{O}(1)$.

To demonstrate the use of the proposed unmatched adjoint operator, images were reconstructed from experimental whole-body mouse PACT data by use of the PLS estimator in Eqn. (48) with a quadratic penalty term specified by the $l_{2}$ norm of the 3D gradient of the object. The 3D PACT dataset was acquired by use of a previously reported small animal imaging system (?). The system employed an arc-shaped transducer array containing 64-elements that spanned 152 degrees over a circle of radius 65 mm. During scanning, the object was rotated over a full 360 degrees with 180 equispaced tomographic views. A Landweber iterative algorithm (?) was employed to compute two different PLS estimates; in one case, $\mathbf{H}_{int}^{\dagger}$ was employed and in the other $\mathbf{H}_{unmatched}^{\dagger}$ was employed.

The reconstructed PACT mouse images with regularization parameter $\lambda=0.001$ are shown in Fig. 6(a) and 6(b). Each figure contains the maximum intensity projection (MIP) images of the reconstructed 3D volumes. The top and bottom rows contain the results corresponding to use of $\mathbf{H}_{int}^{\dagger}$ and $\mathbf{H}_{unmatched}^{\dagger}$, respectively. The columns, from left to right, correspond to the reconstructed images after 40, 100, 200, 400, and 900 iterations, respectively. In addition, Fig. 6(c) displays the objective function values versus iteration number (top row) and computational time (bottom row). These results demonstrate that, for this example, use of $\mathbf{H}_{unmatched}^{\dagger}$ reduced computational times by approximately a factor of 6 while producing images that are visually comparable in image quality to those obtained by use of $\mathbf{H}_{int}^{\dagger}$.

However, as mentioned above for the adjoint-then-discretize approach, employing unmatched adjoint operators in iterative image reconstruction algorithms is not without risk. Improperly designed unmatched adjoint operators may lead to algorithm divergence. In addition, the converged solution obtained by use of an unmatched adjoint operator will generally be different from the true solution. That being said, a carefully-designed unmatched operator can still be of great practical value for accelerating iterative reconstruction for large scale 3D problems (?).

5.2.2 Adjoint for full-wave forward model

Here, the adjoint operator corresponding to the full-wave D-D forward model described in Section 3.3.2 for use with known heterogeneous acoustic media is considered. Given the definition of $\mathbf{H}_{PS}$ in Eqn. (28), the adjoint of $\mathbf{H}_{PS}$ can be computed as

[TABLE]

where the adjoint of the temporal propagator matrices $\mathbf{T}_{k}\in\mathbb{R}^{7NK\times 7NK}$ can be derived from Eqn. (2) (?). The adjoint of the forward operator described in Eqn. (58) can be interpreted as an explicit reversal of the computational steps of the forward scheme (?).

6 Joint reconstruction of initial pressure distribution and acoustic medium parameters

The forward models and image reconstruction methods surveyed above assume either that (1) the to-be-imaged object and surrounding medium are acoustically homogeneous or (2) the object and coupling medium are acoustically heterogeneous and the spatially variant acoustic parameters are known. However, it remains generally difficult and/or inconvenient to accurately estimate the spatially variant acoustic parameters. Accordingly, the vast majority of reported PACT studies assume that the medium is acoustically homogeneous and tune the values of the spatially invariant acoustic parameters to mitigate the impact of acoustic aberrations on the reconstructed image.

Below, a non-conventional approach to PACT image reconstruction is reviewed in which $p_{0}(\mathbf{r})$ and the spatially variant distribution of the acoustic parameters are jointly estimated. This is referred to as a joint reconstruction (JR) approach for PACT (?; ?). A physical motivation for attempting JR results from the observation that spatial variations in the acoustic parameters induce aberrations in the photoacoustic wavefields; consequently, the measured PACT data encodes some information about the acoustic parameters. Several JR methods have been proposed in which the initial pressure distribution and speed of sound (SOS) distribution are estimated concurrently from PACT data alone (?; ?; ?; ?; ?).

Motivated by the earlier numerical investigations, mathematical studies of the JR problem have been conducted (?; ?; ?; ?). These studies established that, when neglecting the discrete sampling effects, the initial pressure distribution and the SOS distribution can be uniquely determined from the measured PACT data only under certain restrictive assumptions (?). In addition, for a linearized version of the JR problem, it was demonstrated that the SOS distribution and $p_{0}(\mathbf{r})$ could not be stably recovered from PACT data alone (?). For a linearized version of the forward problem, where the linearization is represents a smoothing operator for $p_{0}(\mathbf{r})$ and $c(\mathbf{r})^{2}$, the JR problem is unstable in any scale of Sobolev spaces (?). This suggests that solving the JR problem from PACT data alone where the imaging model is described by Eqn. 4 is an extremely challenging undertaking. Various approaches have been proposed to overcome these challenges. A possible approach to mitigate the problem is to augment the PACT measurements with a sparse set of ultrasound tomography data (?).

Below, based on the D-D forward model for use with heterogeneous acoustic media, a general formulation of the JR problem in a discrete setting is provided. Subsequently, a low-dimensional parameterized version of the JR problem is reviewed that holds promise for practical applications.

6.1 JR algorithm

Here, a JR problem is considered in which the spatially variant SOS distribution of the object and coupling medium are unknown, but all other acoustic parameters are known. Let $\mathbf{c}\in\mathbb{R}^{N\times 1}$ denote a finite-dimensional representation of the SOS distribution, the explicit form of which will depend on the choice of numerical scheme that is adopted to solve the wave equation in Eqn. (2).

Consider the D-D PACT forward model in Eqn. (12) that is re-stated as

[TABLE]

where $\mathbf{H}(\mathbf{c})\in\mathbb{R}^{QK\times N}$ denotes a D-D forward model for use with heterogeneous acoustic media as described in Section 3.3.2. The notation $\mathbf{H}(\mathbf{c})$ is employed to make explicit the dependence of forward operator on the unknown discretized SOS distribution $\mathbf{c}$.

The JR problem can be formulated as

[TABLE]

where $R_{c}(\mathbf{c})$ and $R_{\theta}(\boldsymbol{\theta})$ are convex penalty functions whose relative weights are determined by the regularization parameters $\beta_{1}$ and $\beta_{2}$, respectively.

To solve the JR problem in Eqn. (60), an alternating minimization optimization can be employed as described in Algorithm 1 in which two subproblems are solved alternatively until a convergence condition is satisfied (?) . The two subproblems can be expressed as

[TABLE]

For a fixed $\mathbf{c}$, Subproblem 1 corresponds to the conventional PACT reconstruction problem stated in Eqn. (48), which is a convex optimization problem since the objective function is convex for fixed $\mathbf{c}$. In Subproblem 2, $\boldsymbol{\theta}$ is fixed and an estimate of the SOS distribution $\mathbf{c}$ is determined. This corresponds to a non-convex problem. When solving Eqn. (62) by use of gradient-based methods, the computation of the gradient of the data fidelity term with respect to the model parameter $\mathbf{c}$ can be accomplished by use of the adjoint state method (?; ?).

Since Subproblem 2 is non-convex, obtaining its accurate solution represents a challenge. In the subsequent section, some insights obtained from computer-simulation studies that shed a light on the numerical instability and the non-uniqueness properties of the JR problem will be reviewed.

6.1.1 Numerical instability of JR methods

To explore the numerical instability associated with solving Subproblem 2, computer-simulation studies were performed to provide insights into how small perturbations in the assumed $\boldsymbol{\theta}$ affects the accuracy of the reconstructed $\mathbf{c}$ (?). Figures 7(a) and 7(c) display two similar normalized numerical phantoms depicting $\boldsymbol{\theta}$. The RMSE between these phantoms is 0.004. Simulated ideal PACT measurements were computed corresponding to the cases where each of the two $\boldsymbol{\theta}$ was paired with a distinct SOS distribution $\mathbf{c}$. The utilized $\mathbf{c}$ are not shown but is nearly identical to the reconstructed image shown in Fig. 7(b). It is important to note that the $(\boldsymbol{\theta},\mathbf{c})$ pair in the top row of Fig. 7 produces nearly identical PA data to that produced by the $(\boldsymbol{\theta},\mathbf{c})$ pair in the bottom row. The simulated noiseless normalized pressure data at an arbitrary transducer location produced by the two $(\boldsymbol{\theta},\mathbf{c})$ pairs are shown in Fig. 8. The normalized pressure signals in Fig. 8 are observed to overlap almost completely and the RMSE between the two sets of PA data was $3.2e-4$. Figures 7(b) and 7(d) display the reconstructed estimates of $\mathbf{c}$ when the $\boldsymbol{\theta}$ specified in Fig. 7(a) and 7(c) was assumed, respectively. These results demonstrate that the problem of reconstructing $\mathbf{c}$ for a given $\boldsymbol{\theta}$ is ill-posed in the sense that small changes in $\boldsymbol{\theta}$ can produce significant changes in the reconstructed estimate of $\mathbf{c}$. This observation is consistent with the theoretical results (?).

Due to the non-convex nature of Subproblem 2, the optimization algorithm can return a local minimum (i.e., an inaccurate JR solution) as opposed to a global minimum (i.e., an accurate JR solution). Even though the estimate of $\mathbf{c}$ obtained at the local minimum is inaccurate, an accurate estimate of the initial pressure distribution $\boldsymbol{\theta}$ can sometimes be obtained. Various numerical studies have been conducted to show the non-uniqueness properties of the JR algorithm. Figures 9(a) and 9(b) display two numerical phantoms that describe the initial pressure distribution and the speed of sound distribution of a mouse phantom, respectively. The JR algorithm was applied to the forward PACT measurement data generated from the corresponding phantoms. The alternating minimization algorithm described in Algorithm 1 was used to solve the JR algorithm. The results from the JR algorithm are displayed in Fig. 10. From the reconstruction results, one can observe that even though the sound speed distribution is highly inaccurate, the initial pressure distribution can be effectively estimated. This demonstrates that JR algorithms can be employed to improve PACT image quality even if the produced estimates of the sound speed distribution is highly inaccurate.

6.2 Parameterized JR

To mitigate the instability of JR, a modified version of the JR problem can be formulated in which $p_{0}(\mathbf{r})$ and only a lower dimensional estimate of SOS distribution are estimated (?; ?; ?). The SOS representation $\mathbf{c}$ is typically formed by use of pixel expansion functions. To constrain the JR problem, a lower dimensional parameterized representation of the SOS distribution, $\mathbf{c}_{p}\in\mathbb{R}^{D}$, can be introduced as

[TABLE]

where $\boldsymbol{\Gamma}\in\mathbb{R}^{N\times D}$ is a binary matrix that maps $\mathbf{c}_{p}$ to $\mathbf{c}$. Let $\mathcal{I}_{j}$ denote the set of pixel indices in $\mathbf{c}$ that corresponds to the $j^{th}$-parameterized SOS value (i.e., component of $\mathbf{c}_{p}$), $\boldsymbol{\Gamma}$ can be defined as

[TABLE]

The parameterized JR problem is given by (?)

[TABLE]

Once $\mathbf{c}_{p}$ is estimated, an estimate of $\mathbf{c}$ can be obtained by use of Eqn. (63). The gradient of the cost function with respect to $\mathbf{c}_{p}$ can be related to the gradient with respect to $\mathbf{c}$ via the chain rule. When $\Gamma$ corresponds to a linear mapping, one obtains

[TABLE]

where the gradient of the cost function with respect to $\mathbf{c}$ can be computed by use of the adjoint state method (?). Subsequently, an algorithm similar to Algorithm 1 can be employed for image reconstruction, where $\mathbf{c}$ is replaced by $\mathbf{c}_{p}$. Alternatively, a weighted proximal gradient method has been proposed to solve Eqn. (65) in which the necessary gradients can be computed with only two wave solver runs (?).

To demonstrate the use of the parameterized JR algorithm, images were reconstructed from experimental mouse data using the PACT imager described earlier (?). Figure 11 displays zoomed-in regions of images of a live mouse reconstructed by use of an iterative method that employed a constant SOS value of 1.500 mm/$\mu$s (subfigure (a)) and a JR method that employed a two-parameter SOS model (subfigure (b)) (?). To establish the two-parameter SOS model, the outer boundary of the mouse was manually segmented from the estimate of $p_{0}(\mathbf{r})$ that was reconstructed using the fixed SOS value. The simulation grids consisted of $2048\times 2048$ pixels with a pixel size of 0.05 mm. In the JR method, the initial guess for the parameterized sound speed distribution was 1.48 mm/$\mu$s for the background and 1.54 mm/$\mu$s for the mouse, while the initial guess for $\boldsymbol{\theta}$ was the zero-vector. The JR result demonstrates improvement over the image obtained with the constant SOS. The largest difference can be seen in the rightmost surface vessel, which appears as an arc in the constant sound speed image and a point in the JR image. In addition, several interior vessels are better focused in the JR image.

7 Other challenges

There are numerous other challenges and opportunities for research related to image reconstruction for practical applications of PACT that are beyond the scope of this article. Some of these are outlined below.

7.1 Approximating PACT as a 2D problem

In practice, PACT imaging systems are often designed to image 2D sections through a 3D object. This can be accomplished by employing an array of elevationally focused ultrasound transducers. For example, several groups have developed systems that employ a circular ring array (or part of a ring) that surrounds the 2D section to be imaged (?; ?; ?; ?; ?). In such systems, a 3D object can be imaged on a 2D slice-by-slice basis and the resulting images stacked together along the direction perpendicular to the plane of the array to estimate a 3D reconstructed volume. There are practical advantages to such approaches. From a hardware perspective, 2D systems can be less costly to construct than full 3D systems and may not require mechanical scanning of the transducer array to image a 2D slice. When advanced electronics are employed, 2D systems can yield near real-time data-acquisition; this can result in excellent temporal resolution and spatial resolution that is not strongly compromised by physiological motion. Image reconstruction can also be performed in near real-time (?).

However, from an image reconstruction perspective, there are potential challenges associated with such 2D approaches to PACT associated with the fact that PACT is inherently a 3D technique. If (idealized) focused transducers could be employed in a PACT imager to isolate an infinitesimally thin slice of the object, the spherical Radon transform forward model described in Section 2.2 would reduce to a circular Radon transform (?). In this case, the reconstruction problem could be treated in 2D. However, the focused piezoelectric transducers employed in practice do not perfectly reject acoustic signals that arrive from outside the plane of interest. Moreover, the focusing properties of such transducers possess a temporal frequency dependence as described in Section 3.4.1. Such factors and the fact that the PA wavefield propagates according to a 3D wave equation result in a mismatch between the measured data and a canonical 2D PACT forward model. Finally, when 2D PACT images corresponding to different parallel slices of the object are stacked together to estimate a 3D reconstructed volume, the resolution is generally severely compromised in the direction perpendicular to the plane of the array. The presence of artifacts in 2D PACT images due to the presence of 3D out-of-plane scatterers can also significantly degrade image quality. Similar issues have been studied in the geophysics community and medical ultrasound computed tomography communities (?; ?; ?; ?; ?; ?; ?). Experienced practitioners can help mitigate some of these out of plane scattering artifacts by carefully designing the elevational focusing properties of the transducers (?; ?; ?).

The discrete forward models and optimization-based image reconstruction methods reviewed in the previous sections can address these issues and potentially improve the accuracy and spatial resolution of images reconstructed by use of ring-array systems. Namely, to compensate for the transducer focusing effects, the SIRs of the transducers can be incorporated into a 3D discrete forward model. In practice, however, it may be challenging to accurately model the SIR for the case where an acoustic lens is employed. By use of a 3D forward model, $p_{0}(\mathbf{r})$ within a thin 3D volume can be estimated. In such approaches, parameters such as voxel size must be designed carefully to balance the effects of data-incompleteness and model error.

7.2 Spatiotemporal image reconstruction

The majority of PACT reconstruction algorithms assume static imaging conditions; namely, the object of interest is considered to be fixed and independent of time. This assumption, however, is violated in many biomedical applications of PACT (?; ?; ?; ?; ?), where the sought-after initial pressure distribution can depend on time denoted as $p_{0}(\mathbf{r},t)$. Spatiotemporal, or dynamic, image reconstruction methods for PACT seek to reconstruct a sequence of images that correspond to a collection of temporal samples. In a dynamic PACT study, the measurement data collected at a particular time point is referred to as a data frame. It is assumed that $p_{0}(\mathbf{r})$ remains static during the acquisition of each data frame. This condition can approximately be satisfied if the temporal resolution of the imaging system is sufficiently high, such as when fixed arrays of transducers are employed.

A simple approach to dynamic reconstruction is to simply employ a (static) conventional PACT reconstruction method to reconstruct estimates of $p_{0}(\mathbf{r},t)$ corresponding to times $t$ at which the data frames were acquired. Such approaches are referred to as frame-by-frame reconstruction (FBFIR) methods (?; ?; ?; ?). A limitation of FBFIR methods is that they do not exploit statistical correlations between data frames and are therefore computationally and statistically suboptimal (?; ?; ?). Below, a few representative examples of spatiotemporal image reconstruction strategies for dynamic PACT that leverage correlations between data frames are reviewed.

Unlike the FBFIR methods, spatiotemporal image reconstruction (STIR) methods for dynamic PACT jointly reconstruct the sequence of object frame estimates by using all of the data frames. One such method that has received attention is the low rank matrix estimation (LRME)-STIR method (?; ?). The LRME-STIR method exploits the fact that, for many dynamic PACT applications, the image sequence can be approximately described as a low-rank matrix whose rank is typically much smaller than the number of temporal samples. Hence, the LRME-STIR method consists of applying a data denoising step followed by image reconstruction conducted in the domain of the singular system of the low rank data matrix. The performance of the LRME-STIR method was compared with that of conventional FBFIR method with both computer-simulated and 2D experimental data. The results of the study demonstrated that the LRME-STIR method was not only computationally more efficient but also yielded more accurate dynamic PACT images than conventional FBFIR methods (?).

Recently, to improve upon low-rank approximation-based methods, a generalized spatio-temporal modeling framework has been proposed that can encode information about a wide range of dynamics (?). In this approach, the image dynamics are described by an explicit partial differential equation (PDE) model chosen a priori such that it reflects the underlying dynamics. Subsequently, the image sequence and the corresponding motion parameters of the PDE are jointly estimated by minimizing a variational energy cost function. In particular, the motion model employed was based on the popular optical flow equation. The approach was applied successfully to a simulated data as well as to a challenging 3D scenario with experimental data. The results of the study showed that both the simulated and the experimental data-based reconstructions showed a significant improvement in image quality over FBFIR reconstructions (?). Furthermore, the reconstructed motion fields, or the parameters of the PDE provided additional information regarding the dynamics. For the case where the temporal dynamics are relatively simple, low-dimensional parametric models can also be effectively be used to constrain the dynamic image reconstruction problem (?).

In general, there remains an important need for the development of reconstruction methods for dynamic PACT that are computationally efficient and can compensate for important physical factors. For example, there have been no reported methods for dynamic PACT that can compensate for the effects of unknown heterogeneities in the acoustic properties of an object.

7.3 PACT in the presence of elastic media

Conventional PACT reconstruction methods assume that the to-be-imaged object can be described as a fluid acoustic medium. When only soft biological tissues are present, this assumption is well-accepted. However, calcified tissues or bone are more accurately described as an elastic, or solid, medium. In elastic media, the propagation of both longitudinal pressure waves as well as transverse shear waves are supported, each having a distinct speed of propagation (?). As such, to preserve image quality, PACT reconstruction methods for use in applications in which elastic solids are present need to be predicated upon the elastic wave equation as described below.

Let the photoacoustically-induced stress tensor at location $\mathbf{r}\in\mathbb{R}^{3}$ and time $t\geq 0$ be defined as

[TABLE]

where $\boldsymbol{\sigma}^{ij}(\mathbf{r},t)$ represents the stress in the $i^{\text{th}}$ direction acting on a plane perpendicular to the $j^{\text{th}}$ direction. Additionally, let $p_{0}(\mathbf{r})$ denote the photoacoustically-induced initial pressure distribution within the object, and $\dot{\mathbf{u}}(\mathbf{r},t)$ $\equiv$ ($\dot{u}^{1}(\mathbf{r},t),\dot{u}^{2}(\mathbf{r},t),\dot{u}^{3}(\mathbf{r},t)$) represent the vector-valued acoustic particle velocity. Let $\rho(\mathbf{r})$ denote medium’s density distribution and $\lambda(\mathbf{r})$, $\mu(\mathbf{r})$ represent the Lamé parameters that describe the full elastic tensor of the linear isotropic media.

The pressure and shear wave propagation speeds are given by,

[TABLE]

respectively. In a 3D heterogeneous linear isotropic elastic medium with an acoustic absorption coefficient $\alpha(\mathbf{r})$, the propagation of $\dot{\mathbf{u}}(\mathbf{r},t)$ and $\boldsymbol{\sigma}(\mathbf{r},t)$ can be modeled by the following two coupled equations (?; ?; ?; ?):

[TABLE]

Here, tr $(\cdot)$ is the operator that calculates the trace of a matrix and $\mathbf{I}\in\mathbb{R}^{3\times 3}$ is the identity matrix. In Eqn. (69c), it has been assumed that the initial pressure distribution $p_{0}(\mathbf{r})$ is compactly supported in a fluid medium where the shear modulus $\mu(\mathbf{r})=0$. In transcranial PACT (?), for example, this corresponds to the situation where the initial photoacoustic wavefield is produced within the soft tissue enclosed by the skull. Equation (69) represents a generalization (but with a diffusive instead of power law attenuation model) of the canonical PACT forward model for fluid media in Eqn. (2).

For the case of layered elastic media, an analytic image reconstruction formula has been established (?). More generally, a methodology for establishing a D-D forward model based on Eqn. (69) and the associated matched adjoint operator has been established and investigated within the context of transcranial PACT (?). The discretization of the elastic wave equation was based on the finite difference in time domain and space, which facilitated large scale parallelization of the forward and adjoint operators using multiple GPUs. By use of the D-D operators, a variety of optimization-based image reconstruction methods can be employed (?). An alternative formulation of the D-D forward and adjoint operators that was established by use of the k-space pseudo-spectral method has also been reported (?).

There remains an important need for the development of reconstruction methods for transcranial PACT that are computationally efficient and can compensate for important physical factors. For example, there have been no reported methods for transcranial PACT that can compensate for the effects of unknown acoustic heterogeneities within the skull of a subject.

8 Summary

In this topical review, a survey of the acoustic inverse problem in PACT has been presented. Specifically, the PACT image reconstruction problem of estimating the photoacoustically-induced initial pressure distribution has been reviewed within the context of modern optimization-based image reconstruction methodologies. Imaging models that relate the measured photoacoustic wavefields to the sought-after initial pressure distribution were described in their continuous and discrete forms. Subsequently, a description of how important physical factors relevant to PACT can be incorporated into the imaging models was also provided.

Details regarding the implementation of optimization-based PACT reconstruction methods that are not widely available in other review papers were presented. For example, descriptions of the implementations of different discrete adjoint operators employed in PACT were provided. Such adjoint operators are a key component of gradient-based algorithms that are employed to solve optimization-based reconstruction problems. Additionally, non-conventional formulations of the PACT image reconstruction problem were also reviewed, in which the acoustic parameters of the medium are concurrently estimated along with the initial pressure distribution.

The authors would like to thank Dr. Thomas P. Matthews, Dr. Kun Wang and Dr. Chao Huang for their insightful discussions regarding 3D PACT. The authors would also like to thank Dr. Umberto Villa and the anonymous reviewers for their careful reading of the manuscript. This work was supported in part by awards NSF 1614305 and NIH NS102213.