Spatiotemporal Image Reconstruction to Enable High-Frame Rate Dynamic Photoacoustic Tomography with Rotating-Gantry Volumetric Imagers

Photoacoustic computed tomography (PACT), also referred to as optoacoustic tomography, is an emerging and promising imaging modality with broad applications in the field of biomedical imaging [1, 2, 3]. By combining the high spatial resolution of ultrasound imaging with the high soft tissue contrast of optical imaging, PACT offers unique advantages for imaging biological structures while avoiding the use of ionizing radiation. In PACT, a fast laser pulse in the near infrared range illuminates the object. Absorption of optical energy by various molecules within the object (chromophore) induces a localized increase of acoustic pressure through the photoacoustic effect. The acoustic wavefield propagating through the object and coupling medium (water) is subsequently detected by ultrasonic transducers. The measured wavefield data can then be utilized to reconstruct an image that depicts the initial induced pressure distribution within the object.

In preclinical and clinical research, the ability to monitor dynamic physiological processes is of utmost importance for comprehending the progression of diseases and developing new treatments [4, 5, 6]. For example, tumor vascular perfusion is one such dynamic processes that is critical in the study of cancers. High vascular perfusion is indicative of angiogenesis, a well-established hallmark of cancerous growth[7, 8]. Due to its noninvasive nature and the combination of optical contrast and spatial resolution at depths beyond the optical diffusion limit, PACT represents a promising imaging modality for monitoring critical dynamic physiological processes in preclinical and clinical research [9, 10, 11, 12, 13].

Despite its considerable promise, current dynamic PACT technologies suffer from fundamental limitations. They often target two-dimensional (2D) spatial imaging due to shorter data acquisition times and computationally less demanding image reconstruction compared to 3D imaging [14, 15, 16, 17, 11, 18]. Most existing 3D PACT imagers developed to-date utilize a rotating measurement geometry in which the tomographic data are sequentially acquired [19, 2, 20, 21], as opposed to being acquired simultaneously at all views. This design is advantageous because it reduces system costs by employing a limited number of acoustic transducers and associated electronics. However, data-acquisition times for a complete tomographic scan can be tens of seconds. Due to the relatively slow rotational speed, the temporal resolution is significantly limited. While enhancing temporal resolution by using sparsely sampled tomographic data is possible, the associated dynamic image reconstruction problem becomes ill-posed and highly challenging.

Previous studies on dynamic PACT[14, 15, 16, 17, 11, 18, 10, 22], have primarily focused on scenarios where the sufficiently sampled tomographic data can be rapidly acquired. In such cases, a straightforward approach is to employ a frame-by-frame image reconstruction (FBFIR) method[14, 15, 16, 17, 11, 18, 10]. These techniques utilize conventional static image reconstruction methods to estimate a sequence of images from sufficiently sampled tomographic data. The temporal resolution is limited by the duration of the complete data acquisition process. Rapid data acquisition is feasible either with 2D PACT imaging or by leveraging a dense, albeit expensive, static transducer array in 3D imaging. However, for volumetric imagers with sequential scanning strategy, FBFIR methods are not applicable, primarily due to the extended time required to accumulate the complete set of tomographic measurements.

On the other hand, spatiotemporal image reconstruction (STIR) methods estimate a sequence of images simultaneously instead of frame-by-frame, and they have demonstrated their efficacy in accurately reconstructing dynamic objects from sparsely sampled data in various medical imaging modalities, including computed tomography[23], positron emission tomography[24], single photon emission computed tomography[25], and magnetic resonance imaging[26]. While a few STIR techniques have been proposed for PACT, some of them still presume the availability of sufficiently sampled tomographic data and strive for enhanced accuracy and/or reduced computational complexity compared to FBFIR techniques[27]. The STIR methods[28, 29, 30] considering sparsely sampled tomographic data, have relied on the principles of compressed sensing and required specific data sampling schemes that are different than the sequential sampling schemes employed by currently available volumetric imagers.

This study introduces a novel spatiotemporal image reconstruction method based on low-rank matrix estimation (LRME-STIR), which is applicable to currently available sequential 3D PACT imaging systems without requiring hardware modifications. By employing the LRME-STIR technique, it becomes possible to overcome the challenges caused by the sparsely sampled tomographic data. The proposed approach holds the potential to advance the field by enabling accurate and efficient spatiotemporal image reconstruction, thereby facilitating the monitoring of dynamic physiological changes using PACT.

The remainder of the paper is organized as follows. Section 2 presents an imaging model for sequential volumetric imagers and introduces the inverse problem formulation. The proposed LRME-STIR method is described in Section 3. The conducted numerical and experimental studies, and the results are provided in Sections 4 and 5, respectively. Finally, the paper concludes with a discussion in Section 6.

In the context of PACT, the sequential data acquisition strategy commonly involves utilizing one or more rotating or translating ultrasonic transducer arrays within single or multiple acoustic probes[19, 2, 20]. This approach facilitates data collection by rotating the probes along a fixed axis, resulting in the acquisition of a few tomographic measurements at each step, as depicted in Figure 1. These measurements are accumulated sequentially to form a complete tomographic measurement set.

During each step of sequential data acquisition, the object can be considered as static (quasi-static assumption), which is justified by the negligible data acquisition time during each step (on the order of $10^{-4}$ s), significantly shorter than the time between consecutive measurements (on the order of $10^{-1}$ s). An object frame is defined as the short period of time when the object is considered static. The sequence of object frames constitutes the dynamic object. Essentially, each data acquisition step, hereafter referred to as an imaging frame, corresponds to an object frame. Although object and imaging frames may seem interchangeable, the distinction lies in the fact that the set of imaging frames is essentially a subset of the set of object frames because the dynamic object might not be imaged throughout all object frames. Under the quasi-static assumption, the time-dependent object function at the $k$-th imaging frame, specifically the dynamic induced initial pressure distribution, can be represented as $f_{k}(\boldsymbol{r})=f(\boldsymbol{r},k\Delta t)$ for $k=1,...,K$. Here, $K$ represents the number of imaging frames, and $\Delta t$ is the time interval between consecutive imaging frames, which is equal to the laser repetition rate. Finally, the rotation speed of the gantry determines the angular spacing between imaging frames.

Under the quasi-static assumption, the data acquisition process at each imaging frame can be described by a continuous-to-discrete (C-D) imaging model as[31, 27]:

where $\tau$ denotes the fast-time (i.e., the arrival time of acoustic signals to the transducer), and $\Delta T$ corresponds to the fast-time sampling interval. The object function at the $k$-th imaging frame, $f_{k}(\boldsymbol{r})$, is assumed to be bounded and contained within the volume $V$. The scalar $c_{0}$ denotes the speed of sound, which is assumed to be constant throughout the volume $V$. The quantities $\boldsymbol{r}$ and $\boldsymbol{r}_{qk}^{\prime}$ specify the spatial coordinates within $V$ and the location of the $q$-th transducer at the $k$-th imaging frame, respectively. The vector $\boldsymbol{g}_{k}\in\mathbb{R}^{PQ}$ represents lexicographically ordered pressure traces measured by the transducers at the $k$-th imaging frame. Here, $Q$ stands for the number of transducers, and $P$ represents the number of electrical signals recorded by each transducer. The notation $[\boldsymbol{g}_{k}]_{(q-1)P+p}$ refers to the $(q-1)P+p$-th entry of the measurement vector $\boldsymbol{g}_{k}$. Here, the integer-valued indices $q$ and $p$ denote the transducer index and temporal sample, and $\Omega_{qk}$ denotes the detection area of the $q$-th transducer at the $k$-th imaging frame. When the transducer size is small and/or the object is located near the center of a relatively large measurement geometry, an idealized point-like transducer model can be assumed and the surface integral over $\Omega_{qk}$ can be neglected[27].

To facilitate the implementation of an iterative image reconstruction algorithm, a discrete-to-discrete (D-D) imaging model is defined as follow. The spatially continuous object functions $f_{k}$ corresponding to the $k$-th imaging frame are approximated by use of a finite linear combination of spatial expansion functions $\{\psi_{n}(\boldsymbol{r})\}_{n=1}^{N}$:

where $N$ denotes the number of spatial expansion functions. In this study, the expansion functions are piecewise trilinear Lagrangian functions defined on a uniform Cartesian grid. Their expression is given by[31]:

Here, $\boldsymbol{r}=(x,y,z)$ denotes the spatial coordinate and $\boldsymbol{r}_{n}=(x_{n},y_{n},z_{n})$ specifies the location of the $n$-th node of the uniform Cartesian grid. The parameter $\Delta s$ indicates the distance between adjacent grid points.

The coefficients $\{\alpha_{n}^{k}\}_{n=1,k=1}^{NK}$ can be organized into a matrix $\boldsymbol{F}\in\mathbb{R}^{N\times K}$ with entries $[\boldsymbol{F}]_{nk}\equiv\alpha_{n}^{k}$. The $k$-th column of $\boldsymbol{F}$ is denoted by $\boldsymbol{f}_{k}\in\mathbb{R}^{N}$ and represents the discrete approximation of the object function at the $k$-th imaging frame. That is,

where $\boldsymbol{e}_{k}\in\mathbb{R}^{K}$ represents the $k$-th column of the identity matrix in $\mathbb{R}^{K}$ and $\otimes$ denotes the vector outer product. Correspondingly, a frame-dependent D-D imaging model accounting for measurement noise can be expressed as:

where $\underline{\boldsymbol{g}}_{k}$ represents the (noisy) tomographic measurements at the $k$-th imaging frame and $\boldsymbol{\eta}_{k}$ accounts for measurement noise and modeling and discretization errors. The operator $\boldsymbol{H}_{k}\in\mathbb{R}^{QP\times N}$ stems from discretization of the C-D PACT imaging operator associated with the $k$-th imaging frame. In particular, the vector $\boldsymbol{g}_{k}\in\mathbb{R}^{QP}$ representing the action of $\boldsymbol{H}_{k}$ on $\boldsymbol{f}_{k}$ is defined as in Eqn. (1) but with the continuous in space object function $f_{k}(\boldsymbol{r})$ replaced by its finite dimensional approximation $f_{k}^{(N)}(\boldsymbol{r})$ defined in Eqn. (2).

Given the data matrix $\underline{\boldsymbol{G}}=\sum_{k=1}^{K}\underline{\boldsymbol{g}}_{k}\otimes\boldsymbol{e}_{k}\in\mathbb{R}^{PQ\times K}$ and the set of imaging operators $\boldsymbol{H_{k}}$ ($k=1,\ldots,K$) corresponding to each imaging frame, the goal of dynamic PACT image reconstruction is to find a matrix $\boldsymbol{\hat{F}}=\sum_{k=1}^{K}\boldsymbol{\hat{f}}_{k}\otimes\boldsymbol{e}_{k}\in\mathbb{R}^{N\times K}$, whose column $\boldsymbol{\hat{f}}_{k}$ represents an object estimate at the $k$-th imaging frame; for simplicity hereafter $\boldsymbol{\hat{f}}_{k}$ is referred to as the $k$-th frame of the spatiotemporal object estimate. Given the sparse tomographic measurements acquired for each imaging frame, the task of dynamic image reconstruction constitutes a significantly ill-posed inverse problem. A unique estimator, $\boldsymbol{\hat{F}}$, can be obtained by solving the following penalized least squares optimization problem:

where $\text{{R}}(\boldsymbol{F})$ is the regularization term, which is convex but possibly non-smooth. The total data fidelity term $\text{{L}}(\boldsymbol{F})=\sum_{k=1}^{K}\text{{L}}_{k}(\boldsymbol{F})$ is the sum of data fidelity terms $\text{{L}}_{k}(\boldsymbol{F})$ associated with the $k$-th imaging frame. These quantities are defined as

respectively, where $\lVert\cdot\lVert_{F}$ denotes the Frobeniuos norm.

In addition to the inherent ill-posed nature of the considered dynamic image reconstruction problem, significant implementation challenges exist. Firstly, unlike the relatively straightforward task of static image reconstruction in which a single image is estimated, spatiotemporal image reconstruction involves dealing with a considerably higher computational burden as all frames are reconstructed simultaneously, particularly when each frame is 3D in space. Secondly, as the number of frames increases, there is a growing need for memory space during computation, which necessitates effective computational strategies with minimal memory usage.

In numerous biomedical applications, the spatiotemporal object function of interest has been demonstrated to be effectively approximated by a small set of weights $\sigma_{r}$, and functions $u_{r}(\boldsymbol{r})$ and $v_{r}(t)$, depending only on the space or time variables, known as the semiseparable approximation[32, 33, 34, 35, 36, 26, 37]. Consequently, the object function can be expressed as $f(\boldsymbol{r},t)\approx\sum_{r=1}^{R}\sigma_{r}u_{r}(\boldsymbol{r})v_{r}(t)$ [32, 33, 34, 35, 36, 26, 37]. Leveraging the semiseparable approximation, the spatiotemporal reconstruction problem can be reduced to the problem of estimating $R$ weights, and $R$ spatial and temporal functions. In a discretized formulation, this is algebraically equivalent to enforcing a low-rank structure on the matrix $\boldsymbol{F}$.

A penalty scheme can then be imposed on the nuclear norm of the spatiotemporal reconstruction matrix to promote low-rankness. Specifically, the regularization term stemming from the nuclear norm of $\boldsymbol{F}$ is defined as

where $\sigma_{r}$ ($r=1,\ldots,\min(N,K)$) represent the singular values of $\boldsymbol{F}$. This approach not only effectively regularizes the ill-posed inverse problem [38] but also reduces memory demands without sacrificing accuracy. Rather than explicitly storing $\boldsymbol{F}$ in memory, its truncated singular value decomposition $\boldsymbol{U}_{R}\boldsymbol{\Sigma}_{R}\boldsymbol{V}_{R}^{T}$ is stored, resulting in decreased memory requirement. Here, $R<\min(N,K)$ denotes the truncation index, $\boldsymbol{\Sigma}_{R}$ is the diagonal matrix comprising the largest $R$ singular values, $\boldsymbol{U}_{R}$ and $\boldsymbol{V}_{R}$ are matrices with orthonormal columns collecting the left and right singular vectors, respectively, corresponding to the $R$ largest singular vaslues. This approach not only facilitates an efficient approximation of $\boldsymbol{F}$ but also enforces the space-time semiseparability. Specifically, the columns of $\boldsymbol{U}_{R}$ and $\boldsymbol{V}_{R}$ are the algebraic counterpart of the functions $u_{r}(\boldsymbol{r})$ and $v_{r}(t)$, respectively.

To further regularize the inverse problem, under the assumption that the object undergoes a smooth and slow temporal change, another regularization scheme penalizing the difference between two consecutive frames can be employed. This can be accomplished by penalizing the squared Frobenius norm of the temporal difference matrix, which is expressed through the temporal (forward) difference operator, denoted as $\boldsymbol{D}\in\mathbb{R}^{K\times(K-1)}$, applied to $\boldsymbol{F}$:

where, $\boldsymbol{D}$ is defined as:

Here, $\boldsymbol{d_{k}}$ corresponds to the $k$-th column of the temporal (forward) difference operator, $\boldsymbol{D}$. Within this column, the $k$-th and $(k+1)$-th elements possess values of $-1$ and $1$, respectively, while all other elements are set to $0$.

In this way, the sought after estimate of the dynamic reconstruction problem with consideration of a maximum rank constraint and temporal and nuclear norm penalties can be formulated as:

where $\mathbb{R}^{N\times K}_{R_{\max}}$ denotes the set of all $N$-by-$K$ matrices of rank at most $R_{\max}$ and the cost function $\text{{J}}(\boldsymbol{F})$ is comprised of the data fidelity term $\text{{L}}(\boldsymbol{F})$ in Eqn. (7), the temporal regularization term $\text{{R${}^{t}$}}(\boldsymbol{F})$ in Eqn. (9), and the nuclear norm regularization term $\text{{R${}^{nn}$}}(\boldsymbol{F})$ in Eqn. (8). Here, the parameters $\gamma\geq 0$ and $\lambda\geq 0$ control the strength of the temporal and nuclear norm regularization terms, respectively. Similarly to the data fidelity computation in Eqn. (7), the temporal regularization term $\text{{R${}^{t}$}}(\boldsymbol{F})=\sum_{k=1}^{K}\text{{R${}_{k}^{t}$}}(\boldsymbol{F})$ can written as the sum of contributions R${}_{k}^{t}$ from each imaging frame. Each term is defined as

where, for uniformity of notation, $\boldsymbol{d}_{K}\in\mathbb{R}^{K}$ is the zero vector. To highlight the contribution from each frame, the minimization problem in Eqn. (11) can then be rewritten as

In the formulated minimization problem, the data fidelity and temporal penalty terms are convex and smooth, while the nuclear norm penalty term is convex but non-smooth. Accordingly, the minimization problem can be solved by use of a proximal gradient descent (PGD) method [39, 40].

The convergence speed of the PGD method can be improved with momentum schemes [41, 42]. To further improve the convergence speed, especially in the early iterations, an ordered subsets (OS) approach[43] can be incorporated with momentum schemes[44, 45], despite the lack of theoretical guarantees of convergence. In addition to acceleration, applying the OS approach to find an approximate solution to Eqn. (13) offers several other significant benefits. Notably, it substantially reduces memory requirements by a factor proportional to the number $M$ of subset used. While the gradient with respect to all imaging frames requires $\mathcal{O}(NK)$ memory usage, the gradient corresponding to each subset only requires $\mathcal{O}(NK/M)$ storage. Furthermore, it preserve the low-rank structure of the reconstructed object function estimate at each step of the proximal gradient iteration.

For the optimization problem given by Eqn. (11), the OS-based cost function can be expressed as[44, 45]:

where, $\mathcal{K}_{j}$ represents the set of frame indices in the $j$-th ordered subset. The OS approach relies on the “subset balance” approximation[44, 45], implying that $\textbf{J}_{\mathcal{K}_{j}}(\boldsymbol{F})\approx\textbf{J}(\boldsymbol{F})$. The update procedure in PGD with OS consists of two main steps. Initially, a gradient descent step is performed moving in the negative gradient of the smooth components of the OS-based cost function, which yields

Above, $\eta$ is the step size and $\boldsymbol{F}_{j}$ denotes the spatiotemporal object estimate at the beginning of the $j$-th update. The first component of the gradient, associated with the data fidelity term for the $k$-th frame, results in an outer product that produces a rank-$1$ matrix. Similarly, the second component of the gradient, related to the temporal regularization term for the $k$-th frame, also yields a rank-$1$ matrix. Thus, the maximum rank of the gradient of the OS-based cost function $\text{{J}}_{\mathcal{K}_{j}}$ is bounded by $2|\mathcal{K}_{j}|$, where $|\mathcal{K}_{j}|$ denotes the number of imaging frames in the ordered subset.

Following the gradient descent step, the proximal step is executed, where the proximal operator is applied to account for the nonsmooth component of the objective function. The proximal step corresponds to the solution of the following minimization problem:

whose solution can be efficiently implemented via a truncated SVD factorization of $\boldsymbol{F}_{j+\frac{1}{2}}$ and the application of the soft-thresholding operator, $\mathcal{S}(.)$, to its singular values $\{\sigma_{i}\}$.[46] The solution of the Eqn. (16) is expressed as following:

where $\boldsymbol{U}_{j+\frac{1}{2}}$,$\boldsymbol{\Sigma}_{j+\frac{1}{2}}$,$\boldsymbol{V}_{j+\frac{1}{2}}$ stem from the truncated SVD of $\boldsymbol{F}_{j+\frac{1}{2}}$ with maximum rank $R_{\max}$ and $\boldsymbol{\tilde{\Sigma}}_{j+\frac{1}{2}}\coloneqq\mathcal{S}_{\eta\lambda}\left(\boldsymbol{{\Sigma}}_{j+1/2}\right)$. The soft-thresholding operator, $\mathcal{S}(.)$, is defined (component-wise) as below:

The soft-thresholding is computationally efficient and enforces a low-rank structure effectively attenuating the singular values.

Algorithm 1 summarizes the proposed accelerated proximal gradient descent algorithm, integrating both momentum and ordered subset techniques to efficiently find an approximate solution to the minimization problem in Eqn. (11). The algorithm takes the following parameters as input: the maximum allowed rank $R_{\max}$; the threshold $\epsilon$ for the stopping criterion; the regularization parameter $\gamma$ for temporal regularization; the regularization parameter $\lambda$ for nuclear-norm regularization; the step size $\eta$; and the number $M$ of subsets. At the beginning of each iteration, the sequence of frame indices undergoes random shuffling. Within the subsequent inner loop, these shuffled indices are partitioned into $M$ subsets. The gradient pertaining to both the data fidelity and temporal regularization components is computed for these subset frame indices, and the gradient descent step is executed. Then, the proximal mapping associated with nuclear norm regularization is evaluated efficiently by using randomized singular value decomposition (SVD)[47] and soft thresholding. Subsequent to this step, the FISTA [41] acceleration is deployed to enhance the convergence rate. The algorithm terminates when the squared Frobenius norm of the difference between two successive iterations $||\boldsymbol{F}^{(i)}-\boldsymbol{F}^{(i-1)}||_{F}^{2}$, normalized by its maximum $\max_{l\leq i}||\boldsymbol{F}^{(l)}-\boldsymbol{F}^{(l-1)}||_{F}^{2}\bigr{)}$ over the previous iterations, falls below the threshold $\epsilon$ defined by the user. This metric is equivalent to monitoring the norm of the gradient in smooth optimization.[40]

This study presents an accurate and computationally efficient LRME-STIR method for dynamic PACT attuned for commercially available volumetric imagers that employ a rotating measurement gantry in which the tomographic data are sequentially acquired. The implementation of the method was verified by an inverse-crime numerical validation study. The effect of varying number of tomographic measurements per imaging frame and the noise level on the method’s accuracy was investigated in an in-silico numerical study. The experimental study demonstrated the LRME-STIR method’s ability to reconstruct the flow of a contrast agent at a frame rate of 0.1 s, even when only a single tomographic measurement per imaging frame was available. Numerical and experimental studies confirm the accuracy of the proposed technique. Thus, this work will potentially have an immediate and sustained positive impact by creating a new capacity to perform dynamic 3D PACT.

The LRME-STIR method proposed in this study employs an accelerated proximal gradient descent method combined with an ordered subsets approach, distinguishing it from previously proposed LRME-STIR methods[32, 33, 26, 37]. This approach yields a rapid and computationally efficient algorithm with reduced memory demands. Nonetheless, it should be noted that the proposed approach lacks theoretical convergence guarantees. While in theory the method may exhibit instabilities, numerical and experimental studies showcase the robustness of the proposed method when the regularization strength and step size are properly chosen.

Future enhancements to the method may involve exploring different heuristics to enhance its stability[44, 45]. Subsequent investigations may focus on evaluating the method’s efficacy in imaging complex dynamic physiological processes and quantifying relevant parameters, such as wash-in and wash-out rates for tumor vascular perfusion [57, 34]. These evaluations might encompass both in-vivo and in-silico experiments, further advancing the understanding and application of the proposed method.

The author S. Ermilov discloses ownership interest in Photosound Technologies, Inc. Other authors have no relevant financial interests and no potential conflicts of interest to disclose.

This work was supported in part by the National Institutes of Health (NIH Grant Nos. OD023029 and EB031585).

Data and code are available upon request. Please contact Dr. Villa at \linkable[email protected].