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

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 (La Rivière et al 2006, Burgholzer et al 2007, Treeby et al 2010, Deán-Ben et al 2011, Modgil et al 2012). 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 Szabo (1994, 2004)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:Powerlaw} \alpha(\mathbf{r},f) = \alpha_0(\mathbf{r})\,f^y, \nonumber \end{align} \tag{ 1 }$$

where f  is the temporal frequency in MHz, $\alpha_0(\mathbf{r})$ is the spatially varying frequency-independent absorption coefficient in dB MHz^−ycm⁻¹, and y  is the power law exponent that is typically in the range of 0.9–2.0 (Szabo 2004). Note that classical lossy wave equations predict an absorption that is either frequency independent or proportional to frequency squared (Markham et al 1951).

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 (Szabo 1994, 2004, Moshrefi-Torbati and Hammond 1998, Podlubny 1998). Numerical implementation of time-domain fractional derivative operators for accurately modeling 3D acoustic wave propagation poses a significant memory burden (Yuan et al 1999). To overcome this memory burden, some work has focused on devising a recursive algorithm to compute the time-domain fractional derivative (Liebler et al 2004). 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 (Chen and Holm 2004). 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 (Sushilov and Cobbold 2004). In order to account for the dispersion inconsistencies, a lossy wave equation that utilizes two fractional Laplacian derivative operators was proposed (Treeby and Cox 2010a, Treeby et al 2010). 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 Morse and Ingard (1961, 1968) and Treeby and Cox (2010a):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial}{\partial t} \dot{\mathbf{s}}(\mathbf{r},t) = - \frac{1}{\rho_0(\mathbf{r})} \nabla p(\mathbf{r},t) \nonumber \end{align} \tag{ 2a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial}{\partial t} \rho(\mathbf{r},t) = - \rho_0(\mathbf{r}) \nabla \cdot \dot{\mathbf{s}}(\mathbf{r},t) \nonumber \end{align} \tag{ 2b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p(\mathbf{r},t) = c_0(\mathbf{r})^2 \Big\{1 - \mu(\mathbf{r}) \frac{\partial}{\partial t} (-\nabla^2)^{\frac{y}{2} - 1} - \eta(\mathbf{r}) (-\nabla^2)^{\frac{y-1}{2}}\Big\} \rho(\mathbf{r},t),\label{eq:2a_lossiii} \nonumber \end{align} \tag{ 2c }$$

subject to the initial conditions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p(\mathbf{r},t)\Big|_{t = 0} = p_0(\mathbf{r}),\ \ \ \ \dot{\mathbf{s}}(\mathbf{r},t)\Big|_{t= 0}= 0. \nonumber \end{align} \tag{ 2d }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \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). \nonumber \end{align} \tag{ 3 }$$

The second and third terms on the right hand side of equation (2c) account for the required power law absorption and dispersion terms separately through two fractional Laplacian derivative operators. The initial value problem defined in equation (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 equation (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:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:2a_HCC} p(\mathbf{r}^{\prime},t) = \mathcal{H}_{wave} p_0(\mathbf{r}), \nonumber \end{align} \tag{ 4 }$$

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.

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p(\mathbf{r}', t) = \frac{1}{4 \pi c_0^2} \int_V {\rm d} ^3 \mathbf{r} p_0(\mathbf{r}) \frac{d}{dt} \frac{\delta \left(t - \frac{|\mathbf{r}' - \mathbf{r}|}{c_0} \right)}{|\mathbf{r}' - \mathbf{r}|} \label{eq:2b_homo_CC} \equiv \mathcal{H}_{int} p_0(\mathbf{r}), \nonumber \end{align} \tag{ 5 }$$

where c₀ 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) (Kuchment and Kunyansky 2011) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g(\mathbf{r}', t) = \int_V {\rm d} ^3 \mathbf{r} p_0(\mathbf{r}) \delta(c_0 t - |\mathbf{r}' - \mathbf{r}|), \label{eq:2b_SRT} \nonumber \end{align} \tag{ 6 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p(\mathbf{r}', t) = \frac{1}{4 \pi c_0^2} \frac{\partial}{\partial t} \left(\frac{g(\mathbf{r}', t)}{t} \right). \label{eq:2b_g_def} \nonumber \end{align} \tag{ 7 }$$

In the special case of a 2D problem, the spherical Radon transform model reduces to a circular Radon transform (Haltmeier et al 2007, Wang et al 2008).

For the case of a known weakly heterogeneous acoustic media, equation (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 (Miller et al 1987, Modgil et al 2010, Jose et al 2012).

Heuristic strategies have also been proposed to mitigate image artifacts that result from employing the idealized forward model in equation (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 (Anastasio et al 2005a, 2005b, Poudel et al 2017a).

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 (Barrett and Myers 2013) that relates $p_0(\mathbf r)$ to the measurement vector. The C-D operator for PACT can be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{H}_{CD} \equiv \mathcal{D}\mathcal{H}_{CC}, \nonumber \end{align} \tag{ 8 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle [{\mathbf{u}}]_{qK + k} = p(\mathbf{r}^{\prime},t)|_{\mathbf{r}^{\prime} = \mathbf{r}_{q}^{\prime}, t = k \Delta_t}, \nonumber \end{align} \tag{ 9 }$$

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 Wang and Anastasio (2011)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:CDaperture} [{\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), \nonumber \end{align} \tag{ 10 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{D_D:finite_dimensional_representation} p_0^a(\mathbf{r})=\sum_{n=1}^{N} [\boldsymbol \theta]_{n}\phi_{n}({\mathbf{r}}), \nonumber \end{align} \tag{ 11 }$$

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 (Paltauf et al 2002, Ephrat et al 2008, Wang et al 2011, 2013, Wang et al 2014a), and linear interpolation functions (Zhang et al 2009, Dean-Ben et al 2012, Wang et al 2013, Ding et al 2017). In general, these D–D imaging models will have distinct numerical properties that will affect the performance of iterative reconstruction algorithms (Wang and Anastasio 2011).

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{D_D:dis_to_dis_img_mod} {\mathbf{u}} = \mathcal{H}_{CD}\mathbf{p}_0^a(\mathbf r) =\sum_{n=1}^N [\theta]_{n} \mathcal{H}_{CD} \{\phi_{n}(\mathbf{r})\} \equiv \mathbf{H}\boldsymbol{\theta}, \nonumber \end{align} \tag{ 12 }$$

where the D–D operator $\mathbf{H}$ is commonly referred to as the system matrix. An element in the nth row and mth 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.

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 Kak and Slaney (2001)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \phi^{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}), & {{\rm if}}\, |x-x_n|, |y-y_n|, |z-z_n| \leqslant \Delta_s \nonumber \\ 0, & {{\rm otherwise}} \end{array}\right., \label{eq:3b_interpolation_expansion} \nonumber \end{align} \tag{ 13 }$$

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 Wang and Anastasio (2011):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle [\boldsymbol{\theta}]_n = \int_V {\rm d} ^3 \mathbf{r} \delta(\mathbf{r} - \mathbf{r}_n)\,p_0(\mathbf{r}), \quad n = 0,1,\dots,N-1, \label{eq:3b_int_coef} \nonumber \end{align} \tag{ 14 }$$

where $\mathbf{r}_n = (x_n, y_n, z_n){}^T$ denotes the location of the nth Cartesian grid node. The corresponding D–D PACT imaging model based on the interpolation expansion functions can be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathbf{u} = \mathbf{H}_{int} \boldsymbol{\theta}_{int}, \label{eq:3b_int_DD} \nonumber \end{align} \tag{ 15 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:3b_Hint} \mathbf{H}_{int} \equiv \mathbf{D} \mathbf{G}. \nonumber \end{align} \tag{ 16 }$$

Here, $\mathbf{G}$ and $\mathbf{D}$ are discrete approximations of the SRT operator (equation (6)) and the differential operator (equation (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 (Wang et al 2013):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:SRTd} [\mathbf{G} \boldsymbol{\theta}]_{qK+k} \equiv \Delta_s^2 \sum_{n=0}^{N-1} [\boldsymbol{\theta}]_n \sum_{i=0}^{N_i-1} \sum_{j=0}^{N_j-1} \phi_n^{int}(\mathbf{r}_{k,i,j}) \equiv [\mathbf{g}]_{qK+k}, \nonumber \end{align} \tag{ 17 }$$

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 qth transducer $\mathbf{r}_q'$ with a radius of $k \Delta_t c_0$, in which the center of ith polar division and j th azimuth division is denoted by $\mathbf{r}_{k, i, j}$ (Wang et al 2013).

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:differentiald} [\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). \nonumber \end{align} \tag{ 18 }$$

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 equations (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 Lewitt (1990), Wang et al (2014a), Schweiger and Arridge (2017)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 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)}, & {{\rm if~}} 0 \leqslant x \leqslant a \nonumber \\ 0, & {{\rm if~}} x > a, \end{array} \right. \nonumber \end{align} \tag{ 19 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \phi^{KB}_n(\mathbf{r}) = b(|\mathbf{r} - \mathbf{r}_n|), \nonumber \end{align} \tag{ 20 }$$

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 (Wang et al 2012). 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 Wang et al (2014a)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{\mathbf{u}} = \tilde{\mathbf{H}}_{KB} \boldsymbol{\theta}_{KB}, \nonumber \end{align} \tag{ 21 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle [\tilde{\mathbf{H}}_{KB}]_{qL+l, n} \equiv \tilde{p}^{KB}(\,f) \tilde{h}_q^{s, point}(\mathbf{r}_n, f) \vert_{f = l\Delta_f}, {{\rm ~for~}} l = 0, 1, \dots, L-1, \label{eq:3b_HKB} \nonumber \end{align} \tag{ 22 }$$

where

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \tilde{p}^{KB}(\,f) = - \frac{j 4\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}}, \nonumber \end{align*}$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{h}_q^{s, point}(\mathbf{r}_n, f) = \frac{\exp\left(-j 2\pi f \frac{|\mathbf{r}_q^s - \mathbf{r}_n|}{c_0}\right)}{2 \pi |\mathbf{r}_q^s - \mathbf{r}_n|}. \label{eq:3b_SIR_Green} \nonumber \end{align} \tag{ 23 }$$

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 equation (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 mth order spherical Bessel function of the first kind (Diebold 2009a). 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 equation (23), where $\mathbf{r}_q^s$ denotes the location of the qth transducer and $\mathbf{r}_n$ denotes the location of the center of the nth KB expansion function (Wang et al 2014a).

3.3.1. Overview of approaches

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 (Huang et al 2013).

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^{int}(\mathbf{r})\}_{n = 0}^{N-1}$ denotes the set of 3D Cartesian expansion functions defined in equation (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 (Treeby and Cox 2010b). 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 kth time step along the ith direction. Let $\mathbf{p}_k \in \mathbb{R}^{N \times 1}$ represent the acoustic pressure sampled at all spatial grid points at the kth 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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:3c_Prop} \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,\mathbf{0}_{7N \times 1},\cdots,\mathbf{0}_{7N \times 1} \end{array}\right]^T. \nonumber \end{align} \tag{ 24 }$$

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 (Huang et al 2013). From the initial conditions defined in equation (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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:3c_Init} \left[\begin{array}{@{}c@{}} \mathbf{v}_0, \mathbf{0}_{7N \times 1}, \cdots, \mathbf{0}_{7N \times 1} \end{array}\right]^T = \mathbf{T}_0 \boldsymbol{\theta}, \nonumber \end{align} \tag{ 25 }$$

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 (Huang et al 2013).

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:3c_Interp} \mathbf{u} = \mathbf{M} \left[\begin{array}{@{}c@{}} \mathbf{v}_0, \mathbf{v}_1, \cdots, \mathbf{v}_{K-1} \end{array}\right]^T, \nonumber \end{align} \tag{ 26 }$$

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 (Lee and Schachter 1980).

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:3c_Imageeqn} \mathbf{u} = \mathbf{M}\mathbf{T}_{K-1}\cdots \mathbf{T}_1 \mathbf{T}_0 \boldsymbol{\theta}. \nonumber \end{align} \tag{ 27 }$$

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 equation (2) can be therefore expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:3c_HPS} \mathbf{H}_{PS} = \mathbf{M}\mathbf{T}_{K-1}\cdots \mathbf{T}_1 \mathbf{T}_0. \nonumber \end{align} \tag{ 28 }$$

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

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 equation (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 (Xu and Wang 2003a).

When iterative image reconstruction methods are to-be-employed, a natural way to compensate for the EIR and SIR i s to include them in the constructed D–D forward model (Rosenthal et al 2011, Wang et al 2011, Ding et al 2017). 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 (Huang et al 2013).

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 Wang et al (2013)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathbf{H}_{int} = \mathbf{H}^e \mathbf{DG}, \label{eq:3d_int_DD} \nonumber \end{align} \tag{ 29 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathbf{H}^e \mathbf{p}_{ideal} \equiv \mathcal{F}^{-1}\{\mathcal{F}(\mathbf{h}^e) \mathcal{F}(\mathbf{p}_{ideal}) \}. \nonumber \end{align} \tag{ 30 }$$

Here, $[\mathbf{h}^e]_k = \Delta_t h^e(t) \vert_{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 (Wang et al 2012, 2013) and the elements of the system matrix are re-defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle [\mathbf{H}_{KB}]_{qL+l, n} = \tilde{p}^{KB}(\,f) \tilde{h}^e(\,f), \vert_{f = l \Delta_f},{{\rm ~for~}} l = 0, 1, \dots, L-1, \nonumber \end{align} \tag{ 31 }$$

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 (Stepanishen 1971, Lockwood and Willette 1973, Harris 1981) and employed in PACT (Ermilov et al 2009, Rosenthal et al 2011, Wang et al 2011, 2012, Wiskin et al 2012, Queirós et al 2013, Mitsuhashi et al 2014, Sandhu et al 2015, Ding et al 2017). 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 (Wang et al 2013).

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:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:3d_ii_forward} [\mathbf{H}_{KB}]_{qL+l, n} = \tilde{p}^{KB}(\,f) \tilde{h}^e(\,f) \tilde{h}_q^s(\mathbf{r}_n, f) \vert_{f = l \Delta_f}. \nonumber \end{align} \tag{ 32 }$$

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 equation (23) over the qth transducer surface S_q:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{h}_q^s(\mathbf{r}_n, f) = \int_{S_q} {\rm d} \mathbf{r}^3 \frac{\exp (-j 2\pi f \frac{|\mathbf{r}_q^s - \mathbf{r}_n|}{c_0})}{2 \pi |\mathbf{r}_q^s - \mathbf{r}_n|}. \label{eq:3d_SIR_1} \nonumber \end{align} \tag{ 33 }$$

Note that the integral in equation (33) resembles the Rayleigh integral (Kirkup 1994). 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 (Born and Wolf 2013, Mitsuhashi et al 2014):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle |\mathbf{r}_q^s - \mathbf{r}_n| \gg \frac{\max(a,b)^2}{\lambda}. \label{eq:3d_far_field} \nonumber \end{align} \tag{ 34 }$$

Equation (33) can be evaluated to determine $\tilde{h}_q^s(\mathbf{r}_n, f)$ as Stepanishen (1971)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{h}_q^s(\mathbf{r}_n, f) = \frac{ab \exp(-j 2\pi f \frac{|\mathbf{r}^s_q - \mathbf{r}_n|}{c_0})}{2\pi |\mathbf{r}^s_q - \mathbf{r}_n|} {\rm sinc} \left(\pi f \frac{a x_{nq}^{tr}}{c_0 |\mathbf{r}_q^s - \mathbf{r}_n|} \right) {\rm sinc} \left(\pi f \frac{b y_{nq}^{tr}}{c_0 |\mathbf{r}_q^s - \mathbf{r}_n|} \right), \label{eq:3d_SIR_2} \nonumber \end{align} \tag{ 35 }$$

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 Wang et al (2012):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & x_{nq}^{tr} = -x_n \cos \theta_q \cos \phi_q - y_n \cos \theta_q \sin \phi_q + z_n \sin \theta_q, \nonumber \\ & \qquad \qquad y_{nq}^{tr} = -x_n \sin \phi_q + y_n \cos \phi_q. \nonumber \end{align} \tag{ 36 }$$

3.4.3. Patch-based estimation of the SIR

When the far field condition in equation (34) is violated, the SIR expression in equation (35) can be inaccurate and its use may lead to artifacts in the reconstructed images (Mitsuhashi et al 2014). 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 (Rosenthal et al 2011) and patch-based method (Mitsuhashi et al 2014). 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 (Mitsuhashi et al 2014); 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):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{h}^s_q(\mathbf{r}_n, f) \approx \frac{1}{N_{patch}} \sum_{i=1}^{N_{patch}} \tilde{h}^s_{q, i}(\mathbf{r}_n, f), \label{eq:3d_SIR_patch} \nonumber \end{align} \tag{ 37 }$$

where $\tilde{h}_{q,i}^s(\mathbf{r}_n, f)$ denotes the SIR corresponding to ith sub-element. In the patch-based approach, $\tilde{h}_{q,i}^s(\mathbf{r}_n, f)$ can be computed with the aid of equation (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 (Mitsuhashi et al 2014). The numerical phantom used for the study consisted of spherical objects placed within a PACT system as shown in figure 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 (Diebold 2009a), with a semi-analytical SIR specifically for spherical waves (Jensen 1999). During reconstruction, the Kaiser–Bessel function-based D–D PACT imaging model introduced in section 3.2.2 was employed with different SIR models: figure 2(b) assumes a point-like transducer model, figure 2(c) employs the far-field-based SIR model as in equations (33) and (35), and figure 2(d) employs the patch-based model in equation (37). The reconstructed images were 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.

Zoom In Zoom Out Reset image size

A large number of analytic, or non-iterative, tomographic reconstruction methods for PACT are currently available. Several review articles provide detailed descriptions of these (Li and Wang 2009, Kuchment and Kunyansky 2011), 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 (Finch and Patch 2004, Xu and Wang 2005, Finch et al 2007, Kunyansky 2007a, 2011, Haltmeier 2014). 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 (Wang et al 2011). Other ad hoc reconstruction techniques such as inversion of the Radon transform (Kruger et al 1995), as well as delay-and-sum techniques (Hoelen and de Mul 2000) have been employed for PACT image reconstruction.

PACT reconstruction methods that are based on harmonic decomposition have also been proposed (Norton and Linzer 1981). For special geometries such as planar detector geometries, such reconstruction methods can be implemented efficiently by use of the fast Fourier transform (FFT) (Fawcett 1985, Köstli et al 2001). This approach has been extended for use with novel measurement geometries that utilize reverberant cavities (Cox et al 2007a), and for use with closed measurement surfaces for which the eigenfunction of the Dirichlet Laplacian are explicitly known (Kunyansky 2007b). These approaches possess computationally efficient implementations and, in certain cases, can reconstruct 3D images thousands of times faster than backprojection-type methods (Kunyansky 2007b). 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.

As an alternative to the analytic reconstruction approaches discussed above, a variety of time-reversal (TR) based reconstruction methods have been proposed (Xu and Wang 2004, Hristova et al 2008, Stefanov and Uhlmann 2009, Treeby et al 2010, Palacios 2016). 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\geqslant T$. Because $\Omega$ contains no energy $\forall t\geqslant 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$ (Agranovsky and Kuchment 2007). 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 (Hristova 2009).

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 equation (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 (Treeby et al 2010). 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)$ (Treeby et al 2010):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial}{\partial t} \dot{\mathbf{s}}(\mathbf{r},t) = - \frac{1}{\rho_0(\mathbf{r},t)} \nabla p(\mathbf{r},t), \label{eq:4b_lossTRi} \nonumber \end{align} \tag{ 38a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial}{\partial t} \rho(\mathbf{r},t) = - \rho_0(\mathbf{r}) \nabla \cdot \dot{\mathbf{s}}(\mathbf{r},t), \label{eq:4b_lossTRii} \nonumber \end{align} \tag{ 38b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p(\mathbf{r},t) = c_0(\mathbf{r})^2 \Big\{1 + \mu(\mathbf{r}) \frac{\partial}{\partial t} (-\nabla^2)^{\frac{y}{2} - 1} - \eta(\mathbf{r}) (-\nabla^2)^{\frac{y-1}{2}}\Big\} \rho(\mathbf{r},t),\label{eq:4b_lossTRiii} \nonumber \end{align} \tag{ 38c }$$

subject to the conditions:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p(\mathbf{r},t)\Big|_{t = 0} = 0,\ \ \dot{\mathbf{s}}(\mathbf{r},t)\Big|_{t= 0}= 0 \textrm{and}~ p(\mathbf{r},t)\Big|_{\mathbf{r} \in \Omega_0} = u(\mathbf{r},T -t)\Big|_{\mathbf{r} \in \Omega_0}\ . \nonumber \end{align} \tag{ 38d }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:4b_TRfin} \hat{p}_0^{TR}(\mathbf{r}) \equiv p(\mathbf{r},t= T) \textrm{,}~ \forall \textrm{}\mathbf{r} \in \Omega. \nonumber \end{align} \tag{ 39 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{p}_0^{TR}(\mathbf{r}) = \mathcal{H}^{TR}u(\mathbf{r}, t), \nonumber \end{align} \tag{ 40 }$$

where $\mathcal{H}^{TR}:\mathbb{L}_2(\Omega_0 \times [0,T]) \mapsto \mathbb{L}_2(\Omega)$ denotes the time reversal operator described in equation (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 (Hristova et al 2008, Hristova 2009). 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 (Stefanov and Uhlmann 2009, Qian et al 2011). 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 (Qian et al 2011). 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 equations (38a)–(38c) subject to the initial and boundary value conditions:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:4b_modTR} p(\mathbf{r},t)\Big|_{t = 0} = P_\Omega u(\mathbf{r},T),\ \ \dot{\mathbf{s}}(\mathbf{r},t)\Big|_{t= 0}= 0 \textrm{and}~ p(\mathbf{r},t)\Big|_{\mathbf{r} \in \Omega_0} = u(\mathbf{r},T -t)\Big|_{\mathbf{r} \in \Omega_0}, \nonumber \end{align} \tag{ 41 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:4b_Pos} \Delta \phi(\mathbf{r}) = 0, \ \ \phi(\mathbf{r})\Big|_{\mathbf{r} \in \Omega_0} = g(\mathbf{r}, T)\Big|_{\mathbf{r} \in \Omega_0}. \nonumber \end{align} \tag{ 42 }$$

One can summarize equations (41), (39) and (42) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{p}_0^{TR}(\mathbf{r}) = \mathcal{H}_{mod}^{TR} u(\mathbf{r}, t), \nonumber \end{align} \tag{ 43 }$$

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 equations (41) and (42).

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:4b_NS} \hat{p}_0^{NS}(\mathbf{r}) = \sum_{m = 0}^{\infty} (I - \mathcal{H}_{mod}^{TR} \mathcal{H}_{wave})^m \mathcal{H}_{mod}^{TR}u(\mathbf{r},t), \nonumber \end{align} \tag{ 44 }$$

where $\mathcal{H}_{wave}$ is the C–C wave equation-based forward operator in equation (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 (Qian et al 2011).

As the sum defined in equation (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 (Arridge et al 2016). Let us define the ith estimate of the reconstructed initial pressure distribution as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:4b_TrunkNS} \hat{p}_0^{NS,i}(\mathbf{r}) = \sum_{m = 0}^{i}(I - \mathcal{H}_{mod}^{TR} \mathcal{H}_{wave})^m \mathcal{H}_{mod}^{TR}u(\mathbf{r},t). \nonumber \end{align} \tag{ 45 }$$

The partial sum in equation (45) can be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:4b_Nsit} \hat{p}_0^{NS,i}(\mathbf{r}) = \hat{p}_0^{NS,i -1}(\mathbf{r}) - \mathcal{H}_{mod}^{TR}(\mathcal{H}_{wave} \hat{p}_0^{NS,i-1}(\mathbf{r}) - \mathbf{u}). \nonumber \end{align} \tag{ 46 }$$

From equation (46), the NS partial sum can be interpreted as an iterative update step, where the ith iteration refers to the ith partial sum.

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:4c_ii_PLS} \hat{\boldsymbol{\theta}} = \mathop{{\rm argmin}}\limits_{\boldsymbol{\theta}\geqslant 0}\frac{1}{2} ||\mathbf{u} - \mathbf{H}\boldsymbol{\theta}||_{\mathbf{W}}^2 + \gamma R(\boldsymbol{\theta}). \nonumber \end{align} \tag{ 48 }$$

Since the initial pressure distribution is non-negative, the above problem can be constrained so the solution satisfies $\boldsymbol{\theta}\geqslant 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₂ 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 (Golub et al 1999), where the regularization term is specified as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R(\theta) \equiv ||\boldsymbol{\theta}||_\mathbf{P}^2, \nonumber \end{align} \tag{ 49 }$$

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 equation (48).

Gradient-based methods that can be employed to solve equation (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 (Nesterov 1983, 1988). 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 equation (48) (Greenbaum 1997, Paige and Saunders 1982, Gutknecht 2007). 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 (Su et al 2014, Dax 2019, Ghai et al 2019).

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 (Dennis and Moré 1977, Nocedal 1980, Conn et al 1991, Khalfan et al 1993, Wright and Nocedal 1999, Dai 2002).

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₁-norm:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:l1} R(\theta)\equiv ||\boldsymbol{\Phi}\boldsymbol{\theta}||_1, \nonumber \end{align} \tag{ 50 }$$

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 (Sidky and Pan 2008) 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 (Chartrand 2009, Starck et al 2010). Although equation (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 equation (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 (Beck and Teboulle 2009b, Combettes and Pesquet 2011, Parikh et al 2014). 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 (Nesterov et al 2013, Beck and Teboulle 2009a, 2009b, Becker and Fadili 2012, Lee et al 2012, Pan et al 2013, Stella et al 2017).

Examples of PACT images of a live mouse reconstructed by use of two different reconstruction methods are displayed in figure 3. A description of the imaging system and experimental details have been described in previous works (Brecht et al 2009, Ermilov et al 2009). 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 figure 3. The image shown in figure 3(a) was reconstructed by use of a FBP algorithm, while the image in figure 3(b) was reconstructed by use of the PLS estimator in equation (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.

Zoom In Zoom Out Reset image size

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, figure 4 displays images of a simple experimental phantom reconstructed by use of equation (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₂-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.

Zoom In Zoom Out Reset image size

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 equation (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 (Li et al 2017). 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 figure 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 equation (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 figure 5(b). Comparing figures 5(a) and (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.

Zoom In Zoom Out Reset image size

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 equation (48) is the computation of the gradient of the data fidelity term

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi(\boldsymbol{\theta}) \equiv \frac{1}{2}||\mathbf{u} - \mathbf{H}\boldsymbol{\theta}||_{\mathbf{W}}^2 \nonumber \end{align} \tag{ 51 }$$

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

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \nabla \psi(\boldsymbol{\theta}) = \mathbf{H}^{\dagger}\mathbf{W} (\mathbf{H}\boldsymbol{\theta} - \mathbf{u}), \nonumber \end{align} \tag{ 52 }$$

where $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}^{\dagger}$. Although the weighted l₂-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₁-norm, KL-divergence, etc. Furthermore, the weight matrix associated with the weighted l₂-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 (Norton 1999, Plessix 2006). 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 equation (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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}^{\dagger}$.

Although not discussed below, it should be noted that another option for establishing $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}^{\dagger}$ is to employ an adjoint-then-discretize approach (Arridge et al 2016). 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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\hat {\mathbf H}^{\dagger}$ of the D–D adjoint operator $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\hat{\mathbf H}^{\dagger}$ does not accurately mimic $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}^{\dagger}$, the behavior of iterative algorithms can be altered (Zeng and Gullberg 2000). 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 (Zeng and Gullberg 2000, Lou et al 2019).

5.2.1. Adjoint for interpolation-based forward model for homogeneous medium

Here, the interpolation-based D–D forward model described in sections 3.2.1 and 3.4.1 is considered. The matched adjoint operator corresponding to $\mathbf{H}_{int}$ in equation (29) is defined as $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}_{int}^\dagger = \mathbf{G}^\dagger \mathbf{D}^\dagger \mathbf{H}^{e\dagger}$, where $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}(\cdot){}^\dagger$ denotes the transpose of the corresponding D–D operator (i.e. a matrix). These operators can be computed as Wang et al (2013)

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \mathbf{H}^{e\dagger} \mathbf{u} = \mathcal{F}^{-1}\{\mathcal{F}(\mathbf{h}^e) \mathcal{F}(\mathbf{u}) \} \equiv \tilde{\mathbf{p}}_{ideal}, \label{eq:4c_ii_3_He_adjoint} \nonumber \end{align} \tag{ 53 }$$

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle [\mathbf{D}^\dagger \tilde{\mathbf{p}}_{ideal}]_{qK+k} = \frac{\beta}{8 \pi C_p \Delta_t^2 k} ([\tilde{\mathbf{p}}_{ideal}]_{qK+k+1} - [\tilde{\mathbf{p}}_{ideal}]_{qK+k-1}) \equiv [\tilde{\mathbf{g}}]_{qK+k}, \label{eq:4c_ii_3_D_adjoint} \nonumber \end{align} \tag{ 54 }$$

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle [\mathbf{G}^\dagger \tilde{\mathbf{g}}]_{n} = \Delta_s^2 \sum_{q=0}^{Q-1} \sum_{k=0}^{K-1} [\tilde{\mathbf{g}}]_{qK+k} \sum_{i=0}^{N_i-1} \sum_{j=0}^{N_j-1} \phi_n(\mathbf{r}_{k, i, j}) \equiv [\tilde{\boldsymbol{\theta}}]_n. \label{eq:4c_ii_3_G_adjoint} \nonumber \end{align} \tag{ 55 }$$

While the forward operator $\mathbf{H}_{int}$ can be efficiently implemented by use of parallel computing techniques (Wang et al 2013), the adjoint operator $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}_{int}^{\dagger}$ is difficult to implement efficiently. This is partly due to the fact that when implementing $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}^\dagger$ in the interpolation-based model, it is difficult to satisfy the principle of 'partition on target results rather than sources' (Kirk and Wen-Mei 2016) for safe and efficient GPU implementation. In addition, evaluating the backward operator $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{G}^{\dagger}$ relies on expensive atomic operations on the GPU (Wang et al 2013), resulting in a 6–10 times longer runtime for $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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 (Lou et al 2019). Such operators, denoted as $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}^\dagger_{unmatched}$, are referred to as unmatched adjoint operators. An unmatched adjoint operator that approximates $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}_{int}^\dagger$ by use of a simplified voxel-driven model can be defined as:

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \mathbf{H}_{unmatched}^\dagger \equiv \mathbf{G}_{unmatched}^\dagger \mathbf{D}^\dagger \mathbf{H}^{e\dagger}, \label{eq:4c_ii_3_H_unmatched} \nonumber \end{align} \tag{ 56 }$$

where

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle [\mathbf{G}_{unmatched}^\dagger \tilde{\mathbf{g}}]_{n} = \Delta_s^2 \sum_{q=0}^{Q-1}[\tilde{\mathbf{g}}]_{\tilde{k}} \equiv [\tilde{\boldsymbol{\theta}}]_n. \label{eq:4c_ii_3_G_unmatched} \nonumber \end{align} \tag{ 57 }$$

Here, $\newcommand{\e}{{\rm e}} \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:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle [\tilde{\mathbf{g}}]_{\tilde{k}} \equiv (k+1-\tilde{k})[\tilde{\mathbf{g}}]_{qK+k} + (\tilde{k}-k)[\tilde{\mathbf{g}}]_{qK+k+1}, \nonumber \end{align*}$$

with k denoting the integer part of $\tilde{k}$. The operator $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}_{unmatched}^\dagger$ approximates the operation of $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}_{int}^\dagger$, but is much faster to compute. It can be seen from equations (55) and (57) that the computation of the exact adjoint $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{G}^\dagger$ involves $Q K N_i N_j \times \mathcal{O}(1)$ calculations, while the approximated adjoint $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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, $Q K N_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 equation (48) with a quadratic penalty term specified by the l₂ norm of the 3D gradient of the object. The 3D PACT dataset was acquired by use of a previously reported small animal imaging system (Ermilov et al 2009). 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 (Landweber 1951) was employed to compute two different PLS estimates; in one case, $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}_{int}^\dagger$ was employed and in the other $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}_{unmatched}^\dagger$ was employed.

The reconstructed PACT mouse images with regularization parameter $\lambda = 0.001$ are shown in figures 6(a) and (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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}_{int}^\dagger$ and $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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, figure 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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\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 $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mathbf{H}_{int}^\dagger$.

Zoom In Zoom Out Reset image size

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 (Lou et al 2019).

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 equation (28), the adjoint of $\mathbf{H}_{PS}$ can be computed as

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:4c_ii_4_adj_PS} \mathbf{H}_{PS}^{\dagger} = \mathbf{T}_{0}^{\dagger} \mathbf{T}_1^{\dagger} \cdots \mathbf{T}_{K - 1}^{\dagger} \mathbf{M}^{\dagger}, \nonumber \end{align} \tag{ 58 }$$

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

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 equation (2).

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathbf{u} = \mathbf{H}(\mathbf{c}) \boldsymbol{\theta}, \nonumber \end{align} \tag{ 59 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:5_JR} \hat{\boldsymbol{\theta}}, \hat{\mathbf{c}} = \mathop{{\rm argmin}}\limits_{\boldsymbol{\theta} \geqslant 0, \mathbf{c} > 0} \frac{1}{2} ||\mathbf{u} - \mathbf{H}(\mathbf{c})\boldsymbol{\theta}||_{\mathbf{W}}^2 + \beta_1 R_{c}(\mathbf{c}) + \beta_2 R_{\theta}(\boldsymbol{\theta}), \nonumber \end{align} \tag{ 60 }$$

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 equation (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 (Huang et al 2016). The two subproblems can be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \textrm{Subproblem 1:}~\hat{\boldsymbol{\theta}} = \underset{\boldsymbol{\theta} \geqslant 0}{\textrm{argmin}}\ \frac{1}{2} ||\mathbf{u} - \mathbf{H}(\mathbf{c})\boldsymbol{\theta}||_{\mathbf{W}}^2 + + \beta_2 R_{\theta}(\boldsymbol{\theta}) \nonumber \end{align} \tag{ 61 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \textrm{Subproblem 2:}~ \hat{\mathbf{c}} = \underset{\mathbf{c} > 0}{\textrm{argmin}}\ \frac{1}{2} ||\mathbf{u} - \mathbf{H}(\mathbf{c})\boldsymbol{\theta}||_{\mathbf{W}}^2 + \beta_1 R_{c}(\mathbf{c}).\label{eq:5_SP2} \nonumber \end{align} \tag{ 62 }$$

Algorithm 1. Alternating optimization approach to JR of $\boldsymbol{\theta}$ and $\mathbf{c}$.

Input: $\newcommand{\e}{{\rm e}} \boldsymbol{\theta}^0, \mathbf{c}^0, \epsilon_{\theta}, \epsilon_{c}, \beta_1, \beta_2$.

Output: $\hat{\boldsymbol{\theta}}, \hat{\mathbf{c}}$

1: k  =  0, k is the iteration number.

2: while $\newcommand{\e}{{\rm e}} \epsilon_{\theta} < {\epsilon^{\theta}_F}$ and $\newcommand{\e}{{\rm e}} \epsilon_{c} < \epsilon^{c}_F$ do

3:   ik  =  0.

4:   $\boldsymbol{\theta}^{k + 1} = \underset{\boldsymbol{\theta^k} \geqslant 0}{\textrm{argmin}}\ \frac{1}{2} ||\mathbf{u} - \mathbf{H}(\mathbf{c}^k)\boldsymbol{\theta^k}||_{\mathbf{W}}^2 + + \beta_2 R_{\theta}(\boldsymbol{\theta}^k)$.

5:   $\mathbf{c}^{k + 1} = \underset{\mathbf{c}^k > 0}{\textrm{argmin}}\ \frac{1}{2} ||\mathbf{u} - \mathbf{H}(\mathbf{c}^k)\boldsymbol{\theta}^{k+1}||_{\mathbf{W}}^2 + \beta_1 R_{c}(\mathbf{c}^k)$.

6:   $\newcommand{\e}{{\rm e}} \epsilon^{\theta}_F = \frac{|\boldsymbol{\theta}^{k+1} - \boldsymbol{\theta}^k|}{\boldsymbol{\theta}^{k+1}}$

7:   $\newcommand{\e}{{\rm e}} \epsilon^{c}_F = \frac{|\boldsymbol{c}^{k+1} - \boldsymbol{c}^k|}{\boldsymbol{c}^{k+1}}$

8:   k  =  k  +  1,

9: end while

10: $\hat{\boldsymbol{\theta}}\leftarrow \boldsymbol{\theta}^{(k)}$, $\hat{\mathbf{c}} \leftarrow \mathbf{c}^k$

For a fixed $\mathbf{c}$, subproblem 1 corresponds to the conventional PACT reconstruction problem stated in equation (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 equation (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 (Norton 1999, Plessix 2006).

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}$ (Huang et al 2016). Figures 7(a) and (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 figures 7(b) and (d). It is important to note that the $(\boldsymbol{\theta},\mathbf{c})$ pair in the top row of figure 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 figure 8. The normalized pressure signals in figure 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 (d) display the reconstructed estimates of $\mathbf{c}$ when the $\boldsymbol{\theta}$ specified in figures 7(a) and (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 (Stefanov and Uhlmann 2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 (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 figure 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 (Zhang et al 2008, Matthews et al 2018, Poudel et al 2018).

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:5_Parameter} \mathbf{c} = \boldsymbol{\Gamma} \mathbf{c}_p, \nonumber \end{align} \tag{ 63 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:5_gamma} [\boldsymbol{\Gamma}]_{i,j} = \left\{\begin{array}{@{}ll@{}} 1,\ \ \ &i \in \mathcal{I}_j \nonumber \\ 0,\ \ \ &{\rm otherwise} \end{array}\right.. \nonumber \end{align} \tag{ 64 }$$

The parameterized JR problem is given by Matthews et al (2018)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:parameterizedJR} \hat{\boldsymbol{\theta}}, \hat{\mathbf{c}_p} = \underset{\boldsymbol{\theta} \geqslant 0, \mathbf{c}_p > 0}{\textrm{argmin}}\ \frac{1}{2} ||\mathbf{u} - \mathbf{H}(\boldsymbol{\Gamma}\mathbf{c}_p)\boldsymbol{\theta}||_{\mathbf{W}}^2 + \beta_2 R_{\theta}(\boldsymbol{\theta}). \nonumber \end{align} \tag{ 65 }$$

Once $\mathbf{c}_p$ is estimated, an estimate of $\mathbf{c}$ can be obtained by use of equation (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

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \nabla_{\mathbf{c}_p}\Big(\frac{1}{2} ||\mathbf{u} - \mathbf{H}(\boldsymbol{\Gamma}\mathbf{c}_p)\boldsymbol{\theta}||_{\mathbf{W}}^2 \Big) = \Gamma^{\dagger}\nabla_{\mathbf{c}}\Big(\frac{1}{2} ||\mathbf{u} - \mathbf{H}(\mathbf{c})\boldsymbol{\theta}||_{\mathbf{W}}^2\Big), \nonumber \end{align} \tag{ 66 }$$

where the gradient of the cost function with respect to $\mathbf{c}$ can be computed by use of the adjoint state method (Plessix 2006). 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 equation (65) in which the necessary gradients can be computed with only two wave solver runs (Matthews et al 2018).

To demonstrate the use of the parameterized JR algorithm, images were reconstructed from experimental mouse data using the PACT imager described earlier (Li et al 2017). 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)) (Matthews et al 2018). 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\times2048$ 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.

Zoom In Zoom Out Reset image size

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 (Gamelin et al 2009, Ma et al 2009, Xia and Wang 2014, Li et al 2017, Lin et al 2018b). 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 (Buehler et al 2010).

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 (Wang et al 2008). 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 (Huthwaite and Simonetti 2011, Bleistein 2012, Yedlin et al 2012, Auer et al 2013, Duric et al 2014, Gemmeke et al 2017, Wiskin et al 2017). Experienced practitioners can help mitigate some of these out of plane scattering artifacts by carefully designing the elevational focusing properties of the transducers (Duric et al 2014, Li et al 2017, Lin et al 2018a).

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.

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 (Gamelin et al 2009, Buehler et al 2010, Li et al 2010, Xiang et al 2013, Lou et al 2016), 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 (Buehler et al 2010, Li et al 2010, Chatni et al 2012, Xiang et al 2013). A limitation of FBFIR methods is that they do not exploit statistical correlations between data frames and are therefore computationally and statistically suboptimal (Wernick et al 1999, Tsao et al 2003, Lingala et al 2011). 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 (Trémoulhéac et al 2014, Wang et al 2014b). 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 (Wang et al 2014b).

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 (Lucka et al 2018). 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 (Lucka et al 2018). 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 (Chung and Nguyen 2017).

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.

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 (Fry and Barger 1978). 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 \geqslant 0$ be defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{\sigma}(\mathbf{r},t)\equiv \left[\begin{array}{@{}ccc@{}} \boldsymbol{\sigma}^{11}(\mathbf{r},t) &\boldsymbol{\sigma}^{12}(\mathbf{r},t) &\boldsymbol{\sigma}^{13}(\mathbf{r},t) \nonumber \\ \boldsymbol{\sigma}^{21}(\mathbf{r},t) &\boldsymbol{\sigma}^{22}(\mathbf{r},t) &\boldsymbol{\sigma}^{23}(\mathbf{r},t) \nonumber \\ \boldsymbol{\sigma}^{31}(\mathbf{r},t) &\boldsymbol{\sigma}^{32}(\mathbf{r},t) &\boldsymbol{\sigma}^{33}(\mathbf{r},t) \end{array}\right], \nonumber \end{align} \tag{ 67 }$$

where $\boldsymbol{\sigma}^{ij}(\mathbf{r},t)$ represents the stress in the ith direction acting on a plane perpendicular to the j 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)$ $\newcommand{\e}{{\rm e}} \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,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c_l(\mathbf{r}) = \sqrt{\frac{\lambda(\mathbf{r}) + 2 \mu(\mathbf{r})}{\rho(\mathbf{r})}} {{\rm ~and~}} c_s(\mathbf{r}) = \sqrt{\frac{\mu(\mathbf{r})}{\rho\left(\mathbf{r}\right)}}, \nonumber \end{align} \tag{ 68a }$$

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 (Alterman and Karal 1968, Boore 1972, Virieux 1986, Madariaga et al 1998):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \partial_t \dot{\mathbf{u}}\left(\mathbf{r},t\right) + \alpha\left(\mathbf{r}\right) \dot{\mathbf{u}}\left(\mathbf{r},t\right)= \frac{1}{\rho\left(\mathbf{r}\right)}\Big(\nabla \cdot \boldsymbol{\sigma}\left(\mathbf{r},t\right)\Big) \label{eq:subeq1} \nonumber \end{align} \tag{ 69a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \partial_t \boldsymbol{\sigma}\left(\mathbf{r},t\right) = \lambda(\mathbf{r}) {\bf tr}(\nabla\dot{\mathbf{u}}\left(\mathbf{r},t\right)) \mathbf{I} + \mu(\mathbf{r})(\nabla \dot{\mathbf{u}} \left(\mathbf{r},t\right) + \nabla \dot{\mathbf{u}}\left(\mathbf{r},t\right) ^ {T}), \label{eq:subeq2} \nonumber \end{align} \tag{ 69b }$$

subject to the initial conditions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:Ini} \boldsymbol{\sigma}_0(\mathbf{r}) \equiv \boldsymbol{\sigma}(\mathbf{r},t)|_{t = 0} = -\frac{1}{3} p_0(\mathbf{r})\mathbf{I}, \ \ \ \ \ \ \ \dot{\mathbf{u}}\left(\mathbf{r},t\right)|_{t=0} = 0. \nonumber \end{align} \tag{ 69c }$$

Here, ${\bf 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 equation (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 (Huang et al 2012), 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 equation (2).

For the case of layered elastic media, an analytic image reconstruction formula has been established (Schoonover and Anastasio 2011). More generally, a methodology for establishing a D–D forward model based on equation (69) and the associated matched adjoint operator has been established and investigated within the context of transcranial PACT (Mitsuhashi et al 2017). 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 (Poudel et al 2017b). 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 (Javaherian and Holman 2018).

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.