1 Introduction to PACT
Photoacoustic computed tomography (PACT), also known as optoacoustic tomography, is a rapidly emerging imaging technique that holds great promise for biomedical imaging [Kruger:95, Kruger:99, Xu:2002, Xu:2003, OraBook]. PACT is a hybrid technique that exploits the photoacoustic effect for signal generation. When a sufficiently short optical pulse is employed to irradiate an object such as biological tissue, the photoacoustic effect results in the generation of acoustic signals within the object. After propagating out of the object, these signals can be measured by use of wide-band ultrasonic transducers. The goal of PACT in its canonical formulation is to reconstruct an image that represents a map of the initial pressure distribution within the object from knowledge of the measured photoacoustically-induced acoustic signals. The initial pressure distribution is proportional to the absorbed optical energy distribution within the object, which can reveal diagnostically useful information based on endogenous hemoglobin contrast or exogenous contrast if molecular probes are utilized. As such, PACT can be viewed either as an ultrasound mediated optical modality or an ultrasound modality that exploits optical-enhanced image contrast [XuReview].
Over the past few decades, there have been numerous fundamental studies of photoacoustic imaging of biological tissue [Oraev:97, Wang:99, Kruger:99, paltauf:1526, Oraevsky:228, maslov:024006, Esen], and the development of PACT continues to progress at a tremendous rate [anastasio2017special]. Biomedical applications of PACT include small animal imaging [xia2014small, ma2009multispectral, brecht2009whole] and human breast imaging [andreev2000optoacoustic, manohar2005twente, lin2018clinical, becker2018multispectral, lou2017system], to name only a few. For additional information regarding applications of PACT, the reader is referred to the many review articles that have been published on this topic [Li-Review, beard2011biomedical, WangReview, wang2012photoacoustic, ntziachristos2010molecular].
PACT is a computed imaging modality that utilizes an image reconstruction algorithm for image formation. From a physical perspective, the image reconstruction problem in PAT corresponds to an acoustic inverse source problem [Anastasio:IP]. When an idealized imaging scenario is considered, a variety of analytic solutions to the PACT inverse problem are available [Li-Review, rosenthal2013acoustic, kuchment2011mathematics, agranovsky2007reconstruction]; however, in practice, numerous challenges exist that can limit their applicability. Alternatively, optimization-based, or iterative, image reconstruction methods for PACT provide the opportunity to enhance image quality by compensating for physical factors, noise, and data-incompleteness. While such approaches are routinely employed in the broader image reconstruction community, relatively few research groups have explored such modern reconstruction methods for PACT. Although they can be computationally demanding, the advent and use of modern parallel computing technologies [wang2013accelerating] have rendered these reconstruction algorithms feasible for many PACT applications.
In this article, the PACT image reconstruction problem is reviewed within the context of modern optimization-based image reconstruction methodologies. This review is restricted to the problem of estimating the initial pressure distribution and does not address the more complicated problem of recovering the optical properties of an object [saratoon2013gradient, yuan2006quantitative]. Imaging models that relate the measured photoacoustic wavefields to the sought-after object function are described in their continuous and discrete forms. These models will describe physical non-idealities in the data such as those introduced by acoustic inhomogeneity, attenuation, and the response of the imaging system. The basic principles of optimization-based PACT image reconstruction from discrete measurement data are presented, which includes descriptions of forward operators that accurately describe the physics of image acquisition. Furthermore, the derivation of the adjoints of the corresponding forward operators are also reviewed. The adjoint operators facilitate the application of gradient-based approaches in solving the optimization-based PACT image reconstruction problem. Non-conventional formulations of the PACT image reconstruction problem, in which acoustic parameters of the medium are concurrently estimated along with the PACT image, are also introduced and reviewed.
The article is organized as follows. In Section 2, two canonical forward models employed for PACT are reviewed that are based on the acoustic wave equation and integral geometry formulations. The explicit formulation of discrete imaging models are described in Section 3. Optimization-based image reconstruction methods that are based on the discrete imaging models are presented in Section 4. Furthermore, joint reconstruction approaches to PACT image reconstruction whereby acoustic parameters of the medium are concurrently estimated along with the initial pressure distribution are discussed in Section 6. Section 7 concludes the article by providing a brief overview of the challenges and opportunities for research related to image reconstruction for practical applications of PACT.
2 Canonical forward models in their continuous forms
A schematic of a general PACT imaging experiment is shown in Fig. 1. A sufficiently short [DieboldBook] laser pulse is employed to irradiate an object at time and the photoacoustic effect results in the generation of internal pressure distribution inside of the object, which is denoted as , . This pressure distribution subsequently propagates outwards in the form of acoustic waves and is detected by the wide-band point-like ultrasonic transducers located on a measurement aperture . The use of alternative measurement technologies, such as integrating line or area detectors [burgholzer2005thermoacoustic, paltauf2007photoacoustic], have also been widely explored but will not be reviewed here. The measurement aperture is a two-dimensional (2D) surface that partially or completely surrounds the object. The propagated pressure wavefield at time will be denoted as .
Below, two types of 3D canonical forward models that relate the measured propagated pressure wavefield to the sought after object function are described in their continuous forms. These models, which form the basis for mathematical studies of the inverse problem in PACT, do not model transducer characteristics. The effects of finite sampling and other physical factors associated with transducer characteristics are addressed in the discrete versions of these models presented in Section 3.
2.1 Wave-equation based PACT forward model in its continuous form
Consider an idealized PACT experiment in which an object with known heterogeneous acoustic properties is irradiated with laser pulse to photoacoustically induce an initial pressure distribution . 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 denote the medium’s ambient speed of sound (SOS) distribution, and denote the distributions of the medium’s acoustic and ambient densities, respectively. Let denote the particle velocity in the medium. The initial pressure distribution 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 [treeby2010photoacoustic, la2006image, burgholzer2007compensation, modgil2012photoacoustic, dean2011effects]. 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 [szabo1994time, szabo2004diagnostic]
where is the temporal frequency in MHz, is the spatially varying frequency-independent absorption coefficient in dB MHzcm, and is the power law exponent that is typically in the range of 0.9-2.0 [szabo2004diagnostic]. Note that classical lossy wave equations predict an absorption that is either frequency independent or proportional to frequency squared [markham1951absorption].
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 [szabo1994time, szabo2004diagnostic, podlubny1998fractional, moshrefi1998physical]. Numerical implementation of time-domain fractional derivative operators for accurately modeling 3D acoustic wave propagation poses a significant memory burden [yuan1999simulation]. To overcome this memory burden, some work has focused on devising a recursive algorithm to compute the time-domain fractional derivative [liebler2004full]
. However, this approach is heuristic and requiresa priori optimization for each value of . To circumvent this, a lossy wave equation modeling power law absorption using fractional Laplacian operators has been proposed [chen2004fractional]. 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 [sushilov2004frequency]. In order to account for the dispersion inconsistencies, a lossy wave equation that utilizes two fractional Laplacian derivative operators was proposed [treeby2010modeling, treeby2010photoacoustic]. 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 [treeby2010modeling, morse1968theoretical, morse1961linear]:
|subject to the initial conditions|
Here, the spatially varying quantities and describe the acoustic absorption and dispersion proportionality coefficients that are defined as
The second and third terms on the right hand side of Eqn. (2c) account for the required power law absorption and dispersion terms separately through two fractional Laplacian derivative operators. The initial value problem defined in Eqn. (2) describes the propagation of photoacoustically generated pressure data given the spatially varying sound speed , ambient density , and the acoustic absorption coefficient .
Consider that the measured pressure data are recorded outside the support of the object for and
. The continuous PACT forward model consists of the composition of the partial differential equation (PDE) described in Eqn. (2) and an observation operator that restricts the pressure recorded to the measurement surface . Hence, the mapping from the initial pressure distribution to the pressure recorded on the continuous measurement aperture in a acoustically heterogeneous fluid media can be expressed as:
where the wave equation-based forward operator describes a continuous-to-continuous (C-C) mapping between two function spaces. The reason for the use of such terminology is that the wave equation-based forward operator maps a function of continuous variable (as opposed to a discrete variable) to another function of a continuous variable. There is no implication that either the function itself is continous or even that the mapping is continuous.
2.2 Integral geometry-based forward model
In the special case of a lossless and acoustically homogeneous infinite medium, the solution to the wave equation in Eqn. (2) can be expressed as [xu2006photoacoustic]
where denotes the constant sound speed value, denotes the one-dimensional (1D) Dirac delta function, and the integral geometry-based C-C forward operator is denoted as .
Equation (5) can be conveniently expressed in terms of the spherical Radon transform (SRT) [kuchment2011mathematics] as
where the data function is related to the measured pressure data by
In the special case of a 2D problem, the spherical Radon transform model reduces to a circular Radon transform [haltmeier2007thermoacoustic, wang2008elucidation].
For the case of a known weakly heterogeneous acoustic media, Eqn. (6) can be replaced by a generalized Radon transform model that is derived from the wave equation by use of a geometrical acoustics approximation. In that case, the measured pressure signals are related to through integration over nonspherical isochronous surfaces [miller1987new, modgil2010image, jose2012speed].
Heuristic strategies have also been proposed to mitigate image artifacts that result from employing the idealized forward model in Eqn. (6) in the presence of unknown acoustic heterogeneity. For example, half-time and partial-time image reconstruction methods have been proposed for PACT image reconstruction from temporally truncated measurements that exclude components of the measured data that are strongly aberrated [poudel2017mitigation, anastasio2005feasibility, anastasio2005half].
3 Discrete forward models
3.1 Review of semidiscrete and discrete forward modeling
As with any digital imaging system, the data acquired in a PACT experiment represent a finite collection of numbers that form a vector. As such, the PACT forward operator is fundamentally a continuous-to-discrete (C-D) mapping[barrett2013foundations] that relates to the measurement vector. The C-D operator for PACT can be expressed as
where is a discretization operator that spatially and temporally samples the pressure wavefield and represents a C-C PACT forward operator such as or .
Let denote a lexicographically ordered representation of the sampled pressure data and let denote its - component. When ideal point-like ultrasound transducers are assumed, the measured samples of the pressure wavefield are given by
for . Here, and denote the temporal sample index and total number of temporal samples, respectively, and the vectors , , describe the locations of the transducers on the aperture .
More generally, as discussed in Section 3.4, the measured pressure data can be described as [wang2015photoacoustic]
where and 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 can be employed. An -dimensional representation of can be described as
where and the superscript indicates that is an approximation of . The functions are called expansion functions and the expansion coefficients are elements of the -dimensional vector . The goal of image reconstruction is to estimate for a fixed choice of the expansion functions . As described below, different choices for the and rules for determining 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 [ephrat2008three, paltauf2002iterative, wang2011imaging, wang2011imaging, wang2013accelerating, wang2014discrete]
, and linear interpolation functions[zhang2009effects, wang2013accelerating, dean2012accurate, ding2017efficient]. In general, these D-D imaging models will have distinct numerical properties that will affect the performance of iterative reconstruction algorithms [wang2011photoacoustic].
D-D forward models can be established by substitution of a finite-dimensional object representation into the C-D imaging model:
where the D-D operator is commonly referred to as the system matrix. An element in the -th row and -th column of will be denoted by .
Next, examples of D-D PACT imaging models for use with homogeneous and heterogeneous acoustic media are reviewed.
3.2 D-D forward models for use with homogeneous media
In this subsection, two popular D-D PACT imaging models for homogeneous acoustic media are reviewed that employ different choices for the expansion functions.
3.2.1 Interpolation-based D-D PACT model
In the interpolation-based D-D PACT imaging model, the the associated expansion functions can be expressed as [kak2001principles]
where is the distance between neighboring points on an uniform and isotropic Cartesian grid. The coefficient vector can be defined as [wang2011photoacoustic]:
where denotes the location of the -th Cartesian grid node. The corresponding D-D PACT imaging model based on the interpolation expansion functions can be expressed as
The discrete SRT operator 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 [wang2013accelerating]:
where , denotes the numbers of angular divisions over the polar and azimuth directions within the local spherical coordinate system centered at the -th transducer with a radius of , in which the center of -th polar division and -th azimuth division is denoted by [wang2013accelerating].
The differential operator can be implemented as
Due to their large sizes, the explicit storage of system matrices is typically infeasible with current computer technologies except for small scale problems. By use of Eqns. (17) and (18) the action of 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 as the expansion functions. These KB functions are defined as [lewitt1990multidimensional, schweiger2017basis, wang2014discrete]
where , denotes the support radius of the KB function, describes the smoothness of the KB function. The term denotes the modified Bessel function of the first kind of order . The associated expansion function can be expressed as
which is simply the KB function centered at location .
When employing the KB expansion functions, it is convenient to formulate the corresponding D-D PACT imaging model in the temporal frequency domain [wang2012investigation]. Let
denotes the discrete Fourier transform (DFT) of the measurement data. The KB function-based D-D PACT imaging model can be expressed as [wang2014discrete]
where the elements of the system matrix can be written as the multiplication of two terms:
Here, denotes the total number of frequency components and is the frequency sampling interval. The first quantity in Eqn. (22) is the temporal Fourier transform of the acoustic pressure generated by a KB function located at origin, where is the -th order spherical Bessel function of the first kind [diebold2009photoacoustic]. The second quantity is the spatial impulse response (SIR) in the temporal frequency domain. Under the assumption of an idealized point-transducer model, this quantity reduces to the Green function given in Eqn. (23), where denotes the location of the -th transducer and denotes the location of the center of the -th KB expansion function [wang2014discrete].
3.3 D-D forward models for use with heterogeneous acoustic media
3.3.1 Overview of approaches
The establishment of D-D imaging models for lossy, heterogeneous fluid media requires the introduction of finite-dimensional approximations of the object function as well as the known medium acoustic parameters such as and . In principle, these medium acoustic parameters can be estimated from adjunct ultrasound tomography data [Manohar:concomitant, Jinxing:acoustic]. When such adjunct image data are not available, one can attempt to estimate the acoustic parameters concurrently with as described in Section 6. For PACT in heterogeneous media, the choice of finite-dimensional representation of the object and medium parameters is dictated by the numerical method employed to solve the coupled first order differential equations described in Eqn. (2). Most methods for computing a numerical solution of the acoustic wave equation in lossy heterogeneous media fall into one of three major categories: finite difference (FD) methods, finite element (FE) methods, and spectral methods.
The FD and FE methods are also known as local methods because the wave propagation equations are solved at each point based only on conditions at neighboring points. In contrast, spectral methods are global, as information from the entire wavefield is employed to numerically propagate the wavefield. As the spectral methods leverage global information, they allow computation to be performed on coarser grid while maintaining high accuracy [gottlieb1977numerical]. As opposed to FD or FE methods that require grid spacing on the order 10 points per minimum wavelength, spectral methods, in theory, only require 2 points per minimum wavelength [cox2007k]. Due to such relaxed spatial sampling constraints, spectral methods are well suited for large-scale, high frequency PA wave solvers. One of the most popular spectral methods used for solving the acoustic wave equation is the k-space pseudospectral method [mast2001k, tabei2002k, cox2007k]. Because of its popularity in the PACT community, the D-D imaging model presented below will be based upon the k-space pseudospectral time-domain (PSTD) method [cox2007k].
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 [huang2013full].
The finite-dimensional approximations of the object function as well as the medium parameters and are constructed by sampling these quantities on a Cartesian grid. If denotes the set of 3D Cartesian expansion functions defined in Eqn. (13), the initial pressure distribution can be approximated as a finite dimensional vector . To ensure the stability of the k-space PSTD method, the material properties such as the speed of sound , the ambient density , the absorption coefficient , and the wavefield parameters such as pressure , particle velocity , and acoustic density , are sampled at different points of staggered Cartesian grids [treeby2010k]. 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 , represent the discrete approximations of the vector-valued particle velocity and acoustic density over the whole 3D Cartesian grid at the time step along the direction. Let represent the acoustic pressure sampled at all spatial grid points at the time step.
Let denote a vector containing all the wavefield variables at the time step . The image acquisition process in PACT can be mathematically modeled as the propagation of the wavefield parameters forward in time from to as
For specific details regarding the definition of the temporal matrix , the reader is referred to the literature [huang2013full]. From the initial conditions defined in Eqn. (2), the vector can be computed from the initial pressure distribution as
where maps the initial pressure distribution to the initial value of the wavefield variables at time [huang2013full].
In general, the transducer locations 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
where . The choice of the interpolation operation is a parameter that will guide the construction of the matrix . Some of the most commonly employed interpolation strategies include trilinear interpolation or Delaunay triangulation-based interpolation [lee1980two].
The explicit form of the system matrix , that describes the propagation of PA waves based on Eqn. (2) can be therefore expressed as
Here, the subscript stands for the fact that the D-D imaging model is computed using the k-space PSTD method.
3.4 Incorporation of transducer responses in D-D forward models
In practice, the ultrasound transducers employed in a PACT imager are imperfect and result in measurements of the PA wavefield that are averaged over finite temporal and spatial apertures as described in Eqn. (10). In the ultrasound imaging community, the effects of these sampling apertures are characterized by the transducer’s electrical impulse response (EIR) and spatial impulse response (SIR). Failure to account for these effects can result in reconstructed estimates of that contain distortions and/or degraded spatial resolution [xu2003analytic].
When iterative image reconstruction methods are to-be-employed, a natural way to compensate for the EIR and SIR is to include them in the constructed D-D forward model [wang2011imaging, rosenthal2011model, ding2017efficient]. 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 [huang2013full].
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 [wang2013accelerating]
where are the discrete approximation so the differential operator and the SRT operator and
Here, is the EIR signal sampled in the temporal domain, and represent the discrete Fourier transform and inverse discrete Fourier transform, respectively, and 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 [wang2012investigation, wang2013accelerating] and the elements of the system matrix are re-defined as
where samples of are computed as the discrete Fourier transform of .
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 whose temporal response is described by . The SIR represents the measurement of the radiated wavefield by a transducer indexed by that has an idealized EIR.
Various SIR models have been proposed [harris1981review, lockwood1973high, stepanishen1971transient] and employed in PACT [ermilov2009development, wang2011imaging, wang2012investigation, mitsuhashi2014investigation, rosenthal2011model, queiros2013modeling, ding2017efficient, wiskin2012non, sandhu2015frequency]. 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 [wang2013accelerating].
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:
Here, denotes the temporal Fourier transform of that can be computed by integrating the Green function in Eqn. (23) over the -th transducer surface :
Note that the integral in Eqn. (33) resembles the Rayleigh integral [kirkup1994computational]. As a specific example, consider a transducer element that possesses a rectangular detecting surface of area . Under the validity of the far field assumption [mitsuhashi2014investigation, born2013principles]:
Eqn. (33) can be evaluated to determine as [stepanishen1971transient]
where are the transverse components of in a local coordinate system centered at . Given and the transducer location specified in spherical polar coordinates , the values of the transverse coordinates can be computed as [wang2012investigation]:
3.4.3 Patch-based estimation of the SIR
When the far field condition in Eqn. (34) is violated, the SIR expression in Eqn. (35) can be inaccurate and its use may lead to artifacts in the reconstructed images [mitsuhashi2014investigation]. 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 [rosenthal2011model] and patch-based method [mitsuhashi2014investigation]. 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 [mitsuhashi2014investigation]; 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):
where denotes the SIR corresponding to -th sub-element. In the patch-based approach, can be computed with the aid of Eqn. (35). It should be noted that, in the line-detector model, the approximation only accounts for the SIR effect in the direction that is parallel to the straight lines, it still assumes zero thickness (point-like) transducer in the perpendicular direction. The patch-based model employs the far-field approximation for both directions in the transducer plane, therefore providing compensation in both directions. In addition, the patch-based model can be extended beyond planar transducers by decomposing an arbitrary transducer surface into smaller patches to estimate its SIR.
A computer simulation study examining the effects of accurately modeling the SIR for flat rectangular transducers in the context of PACT image reconstruction is discussed below [mitsuhashi2014investigation]. The numerical phantom used for the study consisted of spherical objects placed within a PACT system as shown in Fig. 2(a). The simulated pressure data for each sphere were generated by numerically convolving a closed form expression of waves generated by a solid sphere [diebold2009photoacoustic], with a semi-analytical SIR specifically for spherical waves [jensen1999new]. During reconstruction, the Kaiser-Bessel function-based D-D PACT imaging model introduced in Section 3.2.2 was employed with different SIR models: Fig. 2(b) assumes a point-like transducer model, Fig. 2(c) employs the far-field-based SIR model as in Eqns. (33) and (35), and Fig. 2(d) employs the patch-based model in Eqn. (37). The reconstructed images where obtained by solving a penalized-least-squares optimization problem with quadratic smoothness penalty using the conjugate gradient method. These results show that ignoring transducer SIR in the PACT imaging model leads to significant distortion in reconstructed images. Incorporating an accurate transducer SIR model can greatly reduce such artifacts. For objects far away from the transducers, the far-field-based SIR model may suffice. However, for closer objects, patch-based or other ”divide-and-integrate” SIR models can further improve the reconstruction accuracy.
4 Image reconstruction approaches
In this section, a brief survey of conventional image reconstruction methods for PACT is provided. Subsequently, optimization-based image reconstruction methods that are based upon D-D forward models are presented in Section 5.
4.1 Brief overview of analytic reconstruction methods
A large number of analytic, or non-iterative, tomographic reconstruction methods for PACT are currently available. Several review articles provide detailed descriptions of these [kuchment2011mathematics, Li-Review], 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 [finch2004determining, Xu:2005bp, kunyansky2007explicit, finch2007inversion, kunyansky2011reconstruction, haltmeier2014universal]. 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 [wang2011imaging]. Other ad hoc reconstruction techniques such as inversion of the Radon transform [Kruger:95], as well as delay-and-sum techniques [hoelen2000image] have been employed for PACT image reconstruction.
PACT reconstruction methods that are based on harmonic decomposition have also been proposed [Norton:81]. For special geometries such as planar detector geometries, such reconstruction methods can be implemented efficiently by use of the fast Fourier transform (FFT) [kostli, fawcett1985inversion]. This approach has been extended for use with novel measurement geometries that utilize reverberant cavities [cox2007photoacoustic]
, and for use with closed measurement surfaces for which the eigenfunction of the Dirichlet Laplacian are explicitly known[kunyansky2007series]. These approaches possess computationally efficient implementations and, in certain cases, can reconstruct 3D images thousands of times faster than backprojection-type methods [kunyansky2007series]. Similar to the FBP-type algorithms, the harmonic decomposition-based reconstruction methods typically assume an acoustically homogeneous media and do not compensate for the transducer responses.
4.2 Time-reversal reconstruction methods
As an alternative to the analytic reconstruction approaches discussed above, a variety of time-reversal (TR) based reconstruction methods have been proposed [xu2004time, stefanov2009thermoacoustic, hristova2008reconstruction, treeby2010photoacoustic, palacios2016reconstruction]. 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 in an infinite 3D homogeneous lossless acoustic medium, and let the domain correspond to the interior of a closed measurement surface that encloses the to-be-imaged object. According to Huygen’s principle, there exists a time for which the radiated wavefield inside vanishes . Because contains no energy , one can solve the wave equation backwards in time inside with zero initial conditions and the boundary condition given by the measured data on the surface . As described mathematically below, this process of rewinding the waves backward in time to reconstruct the is defined as TR. When finite sampling effects are ignored, it has been demonstrated that time reversal yields a theoretically exact reconstruction of for the case of a 3D acoustically homogeneous medium if is large enough to allow the energy to escape the domain [agranovsky2007uniqueness]. 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 . When the speed of sound distribution is non-trapping (i.e. all of the acoustic energy escapes the domain ), the errors in the estimates can be bounded as a function of the acquisition time [hristova2009time].
In addition, the same principle of replaying the wavefield back in time can also be applied in lossy fluid media. Intuitively, to compensate for the amplitude loss in the wavefields, the wavefields should be corrected by gain factors. To account for acoustic absorption, the absorption term in Eqn. (2c) must be reversed in sign when computing the time-reversed wavefields. In contrast, to account for the dispersion, the sign of the dispersive speed is left unchanged [treeby2010photoacoustic]. Thus, the TR reconstruction method for use with lossy heterogeneous media solves the following system of equations for the time-reversed pressure wavefield [treeby2010photoacoustic]:
|subject to the conditions:|
Here, 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 is specified by the TR method as
where denotes the time reversal operator described in Eqn. (38).
4.2.2 Modified TR image reconstruction based on a Neumann series method
As described above, the canonical TR method is mathematically exact when the radiated wavefield decays to zero inside the domain for some finite time . However, in even dimensions and/or when the medium is acoustically heterogeneous and non-trapping, this condition may not be satisfied [hristova2008reconstruction, hristova2009time]. 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 [qian2011efficient, stefanov2009thermoacoustic]. The NS-based reconstruction method is accurate when the pressure is known on the whole boundary and is greater than a stability threshold [qian2011efficient]. To define the NS-based reconstruction method, one needs to first introduce the modified TR operator. The modified TR operator for lossy heterogeneous media solves Eqns. (38a), (38b) and (38c) subject to the initial and boundary value conditions:
where is a harmonic extension operator. The harmonic extension operator is defined as , where is the solution to the elliptic boundary value problem
Given the modified TR operator, the NS-based reconstructed initial pressure distribution can be defined as
where is the C-C wave equation-based forward operator in Eqn. (4) and 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 [qian2011efficient].
As the sum defined in Eqn. (44) extends from to , it cannot exactly be implemented. It is interesting to note that the NS method can be interpreted as an iterative time reversal method [arridge2016adjoint]. Let us define the estimate of the reconstructed initial pressure distribution as
The partial sum in Eqn. (45) can be expressed as
From Eqn. (46), the NS partial sum can be interpreted as an iterative update step, where the iteration refers to the partial sum.
5 Optimization-based image reconstruction
A general approach to PACT image reconstruction is to formulate the sought-after estimate of as the solution of an optimization problem. In fact, most modern image reconstruction methods for computed imaging modalities including X-ray computed tomography and magnetic resonance imaging are formulated in this way [fessler1994penalized, wernick2004emission, pan2009commercial]. There are potential practical and conceptual advantages to optimization-based image reconstruction over analytic methods. For example, because they are based on D-D forward models, optimization-based image reconstruction methods can comprehensively compensate for the imaging physics and other physical factors such as responses of the measurement system that are not easily incorporated into an analytic method. Moreover, such methods provide a general framework for incorporating regularization, which can mitigate the effects of measurement noise and data incompleteness. Because optimization-based methods are often implemented by use of iterative algorithms, they are generally more computationally demanding than analytic methods; however the use of modern parallel computing technologies [wang2013accelerating] can render iterative three dimensional model-based reconstruction scheme for arbitrary photoacoustic acquisition geometries feasible for many PACT applications.
Consider a D-D imaging model as described in Section 3. An optimization-based image reconstruction method seeks to determine an estimate of by solving an optimization program that can be generally specified as
Here, is the objective function to be minimized that depends on the D-D forward operator , the measured data , and the vector that specifies the estimate of . The functions and the closed set describe the set of constraints that the solution must satisfy. It is important to note that numerous choices must be made in the design of the reconstruction method. These include the specification of , , and , as well as the optimization algorithm that is employed to solve Eqn. (47). It is the joint specification of these quantities that defines the sought after solution and the numerical properties of the reconstruction method [zhang2009effects]. It should also be noted that many generic iterative methods that have been employed for PACT image reconstruction [paltauf2002iterative, AnastasioTATHT, haltmeier2017analysis] can be interpreted as computing solutions to specific cases of Eqn. (47). Below, some commonly employed optimization programs employed by image reconstruction methods are reviewed.
5.1 Penalized least squares methods
A special case of Eqn. (47) corresponds to a penalized least squared (PLS) estimator, and is given by
Since the initial pressure distribution is non-negative, the above problem can be constrained so the solution satisfies . The objective function in the PLS estimator is expressed as two terms. The quantity is referred to as the data fidelity term that corresponds to a weighted least squares functional. The matrix that defines the weighted norm is symmetric and positive definite, whereas corresponds to the measurement vector. This data fidelity functional is convex and differentiable with respect to . 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 is a penalty term that can be designed to regularize the inverse problem, and is a regularization parameter that controls the amount of regularization. A classic form of regularization is Tikhonov regularization [golub1999tikhonov], where the regularization term is specified as
where is a symmetric positive definite matrix. As the Tikhonov regularization term is differentiable with respect to , a variety of gradient-based methods can be employed to solve Eqn. (48).
Gradient-based methods that can be employed to solve Eqn. (48) can be broadly grouped into two classes: first order methods and second order methods. As part of first order methods, the Taylor approximation used to compute the descent direction involves a linear or first order approximation of the cost function. Hence, information about the gradient of the function is employed to compute the descent direction. Some of the most commonly employed first order methods include the steepest descent method, and the conjugate gradient method. There also exists a class of first order methods that employ information about the past gradients/momentum to speed up the convergence rate of first order algorithms. Such methods are referred as Nesterov methods, whereby linear combinations of present and past gradients are used to compute the descent direction [nesterov1983method, nesterovapproach]. In addition, a family of methods called Krylov-based methods, the subset of which is are the conjugate gradient methods, are also employed to solve smooth, convex programs defined in Eqn. (48) [paige1982lsqr, gutknecht2007brief, greenbaum1997iterative]. 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 [ghai2019comparison, su2014differential, dax2019restarted].
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 [wright1999numerical, dennis1977quasi, nocedal1980updating, conn1991convergence, khalfan1993theoretical, dai2002convergence].
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 based on the -norm:
where is a sparsifying transformation that is chosen such that is sparse in its range. Popular choices for include the discrete wavelet transform or finite difference operators. Total variation (TV) regularization [sidky2008image] 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 [starck2010sparse, chartrand2009fast]. Although Eqn. (50) is convex, it is not differentiable with respect to . However, a variety of non-smooth optimization methods, such as proximal-gradient methods, can be employed to solve Eqn. (48) in this case. Proximal methods are a type of forward-backward splitting approach, which alternates between computing a gradient step on the data fidelity term and a solution to a proximal problem involving the regularization [parikh2014proximal, combettes2011proximal, beck2009fast]. 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 [becker2012quasi, beck2009fast, beck2009fasta, pan2013sparse, nesterov2007gradient, lee2012proximal, stella2017forward].
Examples of PACT images of a live mouse reconstructed by use of two different reconstruction methods are displayed in Fig. 3. A description of the imaging system and experimental details have been described in previous works [ermilov2009development, brecht2009whole]. 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 around the vertical axis. A 180-view dataset was used to reconstruct 3D images of the mouse trunk as shown in Fig. 3. The image shown in Fig. 3(a) was reconstructed by use of a FBP algorithm, while the image in Fig. 3(b) was reconstructed by use of the PLS estimator in Eqn. (48) that employed a TV penalty. The image corresponding to the PLS-TV estimate possesses a higher contrast than the images reconstructed by use of the FBP algorithm and reveal a much sharper body vascular tree. These observations are consistent with the fact that optimization-based image reconstruction methods can produce images that possess different physical and statistical characteristics than the ones produced by analytical method.
An important consideration in the formulation of an optimization-based image reconstruction method is the choice of the penalty function. To demonstrate how different choices can influence image quality, Fig. 4 displays images of a simple experimental phantom reconstructed by use of Eqn. (48) with two different choices for the penalty function. Figures 4 (a) and (b) correspond to use of quadratic smoothness Laplacian regularization and TV regularization, respectively. The quadratic Laplacian penalty corresponds to the -norm of a discrete Laplacian operator acting on the image estimate. In this example, use of the TV penalty resulted in images with reduced artifact levels and enhanced spatial resolution as compared to use of the quadratic Laplacian penalty.
All the images shown in Figures 3 and 4 were reconstructed by use of an iterative image reconstruction model that assumed homogeneous acoustic media and thus could not account for the acoustic heterogeneity of the medium. However, knowledge about the spatial distribution of acoustic heterogeneity of the medium can be incorporated into an iterative reconstruction algorithm based on the imaging model described in Eqn. (28). In vivo experimental studies have been conducted to validate the aforementioned iterative algorithm whereby the spatial distribution of the acoustic properties of the medium are accounted for in the reconstruction algorithm. For the in vivo studies, the experimental data were acquired from a mouse trunk using a small animal imaging system that has been described in detail in earlier works [li2017single]. As opposed to the previous studies where the medium was assumed to be homogeneous, in this study the acoustic variation between the background and the bulk mouse tissue were accounted for in the reconstruction algorithm. A reconstructed initial pressure distribution image for fixed constant sound speed value is shown in Fig. 5(a). In the reconstructed image, strong surface and interior vessel structures are observed. While the interior vessel structures appear out of focus, the surface vessels appear to be in focus. When assuming a homogeneous acoustic media, the image reconstruction algorithm could not produce images where both the surface and the interior vessels could be concurrently focused with a single tuned speed of sound value. To account for this, an image reconstruction algorithm that assigned different speed of sound values to the bulk tissue and the background was employed. The PLS-TV algorithm that was based on the imaging model defined in Eqn. (28) was employed to reconstruct the initial pressure distribution. The reconstructed initial pressure distribution when a heterogeneous speed of sound distribution is used is shown in Fig. 5(b). Comparing Figures 5(a) and 5(b), the surface and interior vessels are observed to be focused concurrently in the latter image, while only the surface vessels appear in focus in the former image.
5.2 Computation of adjoint operators and objective function gradients in PACT
The reconstruction approaches described above are general and applicable to a wide range of computational inverse problems. Here, PACT-specific details regarding the implementation of optimization-based image reconstruction methods are reviewed. Specifically, methods for computing the adjoints of some of the D-D forward operators of Section 3 are presented. The explicit computation of the adjoint of the D-D forward operator facilitates the computation of gradients of data fidelity functionals.
A necessary step in solving optimization problems of the form given in Eqn. (48) is the computation of the gradient of the data fidelity term
with respect to . Formally, this quantity is given by
where denotes the adjoint of the system matrix . 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, 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 . Although the weighted -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 -norm, KL-divergence, etc. Furthermore, the weight matrix associated with the weighted -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 [norton1999iterative, plessix2006review]. In the adjoint state method, the gradient of the data fidelity term with respect to the model parameters are computed through the use of adjoint state variables. For PACT applications, the adjoint state variables are solutions to the adjoint of the wave equation defined in Eqn. (2). Hence, the computational complexity associated with the computation of the gradient of an arbitrary data fidelity term through the use of the adjoint state method would at least be on the same order as computing the action of the discrete adjoint operator .
Although not discussed below, it should be noted that another option for establishing is to employ an adjoint-then-discretize approach [arridge2016adjoint]. 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 of the D-D adjoint operator . From an implementation perspective, such approaches can sometimes be more convenient than computing the adjoint operator corresponding to the assumed D-D forward operator . However, in cases where does not accurately mimic , the behavior of iterative algorithms can be altered [zeng2000unmatched]. 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 [zeng2000unmatched, lou2018analysis].
5.2.1 Adjoint for interpolation-based forward model for homogeneous medium
Here, the interpolation-based D-D forward model described in Secs. 3.2.1 and 3.4.1 is considered. The matched adjoint operator corresponding to in Eqn. (29) is defined as , where denotes the transpose of the corresponding D-D operator (i.e., a matrix). These operators can be computed as [wang2013accelerating]
While the forward operator can be efficiently implemented by use of parallel computing techniques [wang2013accelerating], the adjoint operator is difficult to implement efficiently. This is partly due to the fact that when implementing in the interpolation-based model, it is difficult to satisfy the principle of “partition on target results rather than sources” [kirk2016programming] for safe and efficient GPU implementation. In addition, evaluating the backward operator relies on expensive atomic operations on the GPU [wang2013accelerating], resulting in a 6-10 times longer runtime for than for .
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 [lou2018analysis]. Such operators, denoted as , are referred to as unmatched adjoint operators. An unmatched adjoint operator that approximates by use of a simplified voxel-driven model can be defined as:
Here, and the value of is approximated by linearly interpolating from its two neighboring samples:
with denoting the integer part of . The operator approximates the operation of , but is much faster to compute. It can be seen from Eqns. (55) and (57) that the computation of the exact adjoint involves