1.IntroductionPhotoacoustic 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–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. The 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–6 For example, tumor vascular perfusion is a dynamic process 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–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 with 3D imaging.11,14–18 Most existing 3D PACT imagers developed to date utilize a rotating measurement geometry in which the tomographic data are sequentially acquired2,19–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. Although enhancing temporal resolution using sparsely sampled tomographic data is possible, the associated dynamic image reconstruction problem becomes ill-posed and highly challenging. Previous studies on dynamic PACT10,11,14–18,22 have primarily focused on scenarios in which 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.10,11,14–18 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. This limitation arises because conventional static image reconstruction techniques require densely sampled tomographic measurements for accurate object estimates; if sparsely sampled measurement data are used, the reconstructions suffer from severe artifacts.23–25 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,26 positron emission tomography,27 single photon emission computed tomography,28 and magnetic resonance imaging.29 Although 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 with FBFIR techniques.30 The STIR methods,31–33 considering sparsely sampled tomographic data, rely on the principles of compressed sensing and require specific data sampling schemes that are different than the sequential sampling schemes employed by currently available volumetric imagers. This study introduces a novel STIR 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 STIR, 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 Sec. 3. The conducted numerical and experimental studies, as well as their results, are provided in Secs. 4 and 5, respectively. Finally, the paper concludes with a discussion in Sec. 6. 2.Imaging Model for Sequential Volumetric Imagers and Inverse Problem FormulationIn 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.2,19,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 Fig. 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 to be static (quasi-static assumption), which is justified by the negligible data acquisition time during each step (on the order of ), significantly shorter than the time between consecutive measurements (on the order of ). 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 - th imaging frame, specifically the dynamic induced initial pressure distribution, is represented as for . Here, represents the number of imaging frames, and 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 is described by a continuous-to-discrete (C-D) imaging model as30,34 where denotes the fast-time (i.e., the arrival time of acoustic signals to the transducer) and corresponds to the fast-time sampling interval. The object function at the ’th imaging frame, , is assumed to be bounded and contained within volume . The scalar denotes the speed of sound, which is assumed to be constant throughout volume . The quantities and specify the spatial coordinates within and the location of the ’th transducer at the ’th imaging frame, respectively. The vector represents lexicographically ordered pressure traces measured by the transducers at the ’th imaging frame. Here, stands for the number of transducers, and represents the number of electrical signals recorded by each transducer. The notation refers to the -th entry of the measurement vector . Here, the integer-valued indices and denote the transducer index and temporal sample, respectively, and denotes the detection area of the ’th transducer at the ’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 can be neglected.30To facilitate the implementation of an iterative image reconstruction algorithm, a discrete-to-discrete (D-D) imaging model is defined as follows. The spatially continuous object functions corresponding to the ’th imaging frame are approximated using a finite linear combination of spatial expansion functions and are given as where 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 as34 where denotes the spatial coordinate and specifies the location of the ’th node of the uniform Cartesian grid. The parameter indicates the distance between adjacent grid points.The coefficients are organized into a matrix with entries . The ’th column of is denoted by and represents the discrete approximation of the object function at the ’th imaging frame, that is, where represents the ’th column of the identity matrix in and denotes the vector outer product. Correspondingly, a frame-dependent D-D imaging model accounting for measurement noise is expressed as where represents the (noisy) tomographic measurements at the ’th imaging frame and accounts for the measurement noise and the modeling and discretization errors. The operator stems from discretization of the C-D PACT imaging operator associated with the ’th imaging frame. In particular, the vector representing the action of on is defined as in Eq. (1) but with the continuous in space object function replaced by its finite dimensional approximation defined in Eq. (2).Given the data matrix and the set of imaging operators () corresponding to each imaging frame, the goal of dynamic PACT image reconstruction is to find a matrix , with its column representing an object estimate at the ’th imaging frame; for simplicity, hereafter is referred to as the ’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, , is obtained by solving the following penalized least squares optimization problem: where is the regularization term, which is convex but possibly non-smooth. The total data fidelity term is the sum of data fidelity terms associated with the ’th imaging frame. These quantities are defined as respectively, where denotes the Frobeniuos norm.In addition to the inherent ill-posed nature of the considered dynamic image reconstruction problem, significant implementation challenges exist. First, unlike the relatively straightforward task of static image reconstruction in which a single image is estimated, STIR involves dealing with a considerably higher computational burden as all frames are reconstructed simultaneously, particularly when each frame is 3D in space. Second, 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. 3.Low-rank Matrix Estimation-based Spatiotemporal Image ReconstructionIn numerous biomedical applications, the spatiotemporal object function of interest has been demonstrated to be effectively approximated by a small set of weights and functions and , depending only on the space or time variables, known as the semiseparable approximation.29,35–40 Consequently, the object function is expressed as .29,35–40 Leveraging the semiseparable approximation, the spatiotemporal reconstruction problem can be reduced to the problem of estimating weights and spatial and temporal functions. In a discretized formulation, this is algebraically equivalent to enforcing a low-rank structure on the matrix . 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 is defined as where () represent the singular values of . This approach not only effectively regularizes the ill-posed inverse problem41 but also reduces memory demands without sacrificing accuracy. Rather than explicitly storing in memory, its truncated singular value decomposition (SVD) is stored, resulting in a decreased memory requirement. Here, denotes the truncation index, is the diagonal matrix comprising the largest singular values, and and are matrices with orthonormal columns collecting the left and right singular vectors, respectively, corresponding to the largest singular values. This approach not only facilitates an efficient approximation of but also enforces the space-time semiseparability. Specifically, the columns of and are the algebraic counterpart of functions and , 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 is accomplished by penalizing the squared Frobenius norm of the temporal difference matrix, which is expressed through the temporal (forward) difference operator, denoted as , being applied to , given as where is defined as where corresponds to the ’th column of the temporal (forward) difference operator, . Within this column, the ’th and ’th elements possess values of and 1, respectively, and 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 is formulated as where denotes the set of all -by- matrices of rank at most and the cost function is comprised of the data fidelity term in Eq. (7), the temporal regularization term in Eq. (9), and the nuclear norm regularization term in Eq. (8). Here, the parameters and control the strength of the temporal and nuclear norm regularization terms, respectively. Similar to the data fidelity computation in Eq. (7), the temporal regularization term is written as the sum of contributions from each imaging frame. Each term is defined as Where, for uniformity of notation, is the zero vector. To highlight the contribution from each frame, the minimization problem in Eq. (11) is then rewritten asIn the formulated minimization problem, the data fidelity and temporal penalty terms are convex and smooth, wheras the nuclear norm penalty term is convex but non-smooth. Accordingly, the minimization problem can be solved using a proximal gradient descent (PGD) method.42,43 The convergence speed of the PGD method can be improved with momentum schemes.44,45 To further improve the convergence speed, especially in the early iterations, an ordered subsets (OS) approach46 can be incorporated with momentum schemes,47,48 despite the lack of theoretical guarantees of convergence. In addition to acceleration, applying the OS approach to find an approximate solution to Eq. (13) offers several other significant benefits. Notably, it substantially reduces memory requirements by a factor proportional to the number of subsets used. The gradient with respect to all imaging frames requires memory usage, whereas the gradient corresponding to each subset only requires storage. Furthermore, it preserves the low-rank structure of the reconstructed object function estimate at each step of the proximal gradient iteration. For the optimization problem given by Eq. (11), the OS-based cost function is expressed as47,48 where represents the set of frame indices in the ’th ordered subset. The OS approach relies on the “subset balance” approximation,47,48 implying that . The update procedure in PGD with OS consists of two main steps. Initially, a gradient descent step is performed; it moves in the negative gradient of the smooth components of the OS-based cost function, which yieldsHere, is the step size, and denotes the spatiotemporal object estimate at the beginning of the ’th update. The first component of the gradient, associated with the data fidelity term for the ’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 ’th frame, also yields a rank-1 matrix. Thus, the maximum rank of the gradient of the OS-based cost function is bounded by , where denotes the number of imaging frames in the ordered subset. Following the gradient descent step, the proximal step is executed by applying the proximal operator to account for the nonsmooth component of the objective function. The proximal step corresponds to the solution of the following minimization problem: the solution of which can be efficiently implemented via a truncated SVD factorization of and the application of the soft-thresholding operator, , to its singular values .49 The solution of the Eq. (16) is expressed as follows: where , , and stem from the truncated SVD of with maximum rank and . The soft-thresholding operator, , is defined (component-wise) asThe soft-thresholding is computationally efficient and enforces a low-rank structure that effectively attenuates the singular values. Algorithm 1 summarizes the proposed accelerated PGD algorithm, integrating both momentum and ordered subset techniques to efficiently find an approximate solution to the minimization problem in Eq. (11). The algorithm takes the following parameters as input: the maximum allowed rank ; the threshold for the stopping criterion; the regularization parameter for temporal regularization; the regularization parameter for nuclear-norm regularization; the step size ; and the number 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 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 using randomized SVD50 and soft thresholding. Subsequent to this step, the fast iterative shrinkage-thresholding algorithm44 (FISTA) acceleration scheme is deployed to enhance the convergence rate. The algorithm terminates when the squared Frobenius norm of the difference between two successive iterations , normalized by its maximum over the previous iterations, falls below the threshold defined by the user. This metric is equivalent to monitoring the norm of the gradient in smooth optimization.43 Algorithm 1LRME-STIR.
4.Study Description4.1.3D PACT Imaging System Specifications for Experimental and Numerical StudiesThe TriTom preclinical imaging system19,51 developed by PhotoSound was employed for the experimental studies and emulated for the numerical in-silico studies. It integrates photoacoustic (PA) and fluorescence (FL) imaging modalities, harnessing their individual strengths. The system comprises a central rotary scanning stage, an optical excitation setup for PA and FL imaging, a curvilinear 96-element PA transducer array, and a fluorescence-enabled scientific complementary metal-oxide-semiconductor (sCMOS) camera. The configuration of the TriTom setup is depicted in Fig. 2(a), and the imaging chamber is shown in Fig. 2(b). During imaging, the object is immersed in water and continuously rotated while being optically stimulated by a short laser pulse with a 10 Hz repetition rate. The maximum rotation speed is 10 degrees per second; therefore, completing a full 360-deg scan requires 36 s. Four optical fiber bundles are located on the outer circumference of the cylindrical imaging chamber, perpendicular to the scanning plane of the PA array. The laser emits pulses in the 670 to 1064 nm wavelength range at a frequency of 10 Hz, with each pulse lasting 5 ns. Acoustic waves generated by the laser excitation are detected by the PA transducer array, which comprises piezoelectric transducer elements, each measuring . The center frequency of the transducer elements is (at ) with bandwidth . The array is vertically oriented and cylindrically focused, with the central element positioned 65 mm from the center of the imaging chamber. For fluorescence imaging, an sCMOS camera with a resolution and a field-of-view is placed outside the imaging chamber. In the acoustic modeling of the TriTom imaging system, the transducer array was assumed to consist of idealized point-like transducers positioned at the central locations of the transducer elements. The acoustic simulation was implemented with a GPU-accelerated D-D imaging model34 assuming an acoustically homogeneous medium. Although the TriTom system has a single transducer arc, other existing 3D PACT designs (e.g., Refs. 20 and 52) feature multiple transducer arcs; thus the numerical studies also explored scenarios in which multiple tomographic measurements were acquired per imaging frame. Specifically, measurement configurations in which two transducer arcs separated by 90 deg and four transducer arcs separated by 45 deg were considered, as illustrated in Fig. 3. The associated compression ratios for sparse sampling in these scenarios are 1/360, 1/180, and 1/90 for the measurement configuration with 1, 2, and 4 transducer arcs, respectively. The sampling rate for both the experimental and numerical studies was set to 31.25 MHz, with 2048 temporal samples collected per imaging frame. 4.2.Inverse Crime Validation StudyTo verify that Algorithm 1 was correctly implemented, an inverse crime validation study was conducted in silico in which a simple rank-4 dynamic phantom was employed. The phantom consisted of spatial voxels and 360 object frames, with a voxel size of . The induced pressure in the phantom was assumed constant along the -axis and piecewise constant within the -plane at each object frame. The phantom was structured into four distinct regions, each characterized by varying temporal activities as shown in Fig. 4. Figure 4(a) illustrates the central -slice of the phantom at the 120th object frame, and Fig. 4(b) displays the time activity curves (TACs) corresponding to the numbered regions. For each imaging frame (360 in total), four tomographic measurements separated by 45 deg were considered, as depicted in Fig. 3(b). Simulated acoustic pressure data were generated using a grid voxel size of , assuming the phantom was centered within the imaging system.34 The speed of sound was assumed constant at 1495 m/s, and no noise was added to the simulated measurement data. During image reconstruction, the same computational grid that was used for generating the measurement data was employed. In the algorithm, the maximum allowed rank, , was set to 4, and no temporal or nuclear norm penalties were applied (). The step-size, , was tuned empirically to ensure convergence. To explore the impact of the number of OS used in the randomized evaluation of the data fidelity term, three different numbers of OS were investigated: . 4.3.Numerical Phantom StudyA dynamic numerical phantom was utilized to evaluate the performance of the proposed method through in silico experiments. The phantom consisted of 360 object frames and contained four convex ellipsoidal blobs within a larger ellipsoidal blob, along with a vasculature mimicking structure, as shown in Fig. 5(d). The size of the numerical phantom was , with a voxel size of . The time activity at each voxel was designed to mimic a contrast agent’s flow along the paths from ellipsoidal blobs 1 to 2 and 3 to 4. Figure 5(c) illustrates the time activity at the center of each ellipsoidal blob, and the blobs are numbered in Fig. 5(d). The singular values of the dynamic numerical phantom are shown in Fig. 5(a), where a rapid singular value decay is observed. The rapid singular decay indicates that the phantom can be accurately approximated with a low-rank representation, as demonstrated by the mean squared error (MSE) versus rank plot depicted in Fig. 5(b). To generate the synthetic measurement data, the grid voxel size was set to , and 360 imaging frames (corresponding to a complete rotation of the system) were simulated. The speed of sound was assumed to be constant at . Zero-mean Gaussian noise with a standard deviation equivalent to a specified percentage of the maximum value of the simulated data (1%, 3%, or 5%) was added to the pressure signals. The primary objective of this numerical study was to investigate the effect of the following physical factors on the proposed method.
4.4.Numerical StudyFor image reconstruction, the maximum allowed rank in the algorithm was set to as it allows for an accurate approximation of the dynamic numerical phantom with an MSE of around , as shown in Fig. 5(d). To avoid the discretization inverse crime, a coarser grid with voxel size was used for the reconstruction. The speed of sound value was the same as that used for simulating the data. The threshold for the stopping criterion of Algorithm 1 was set to . The step-size, , was tuned empirically to ensure convergence. The number of subsets, , was set to 18. To select appropriate regularization parameters, the balancing principle53 was employed to reduce the number of tunable regularization parameters from two to one. The balancing principle rescales each regularization term by an estimate of their value at the object function . An iterative procedure was proposed in Ref. 53 to estimate such values; however, for simplicity, this work assumes direct knowledge of the actual nuclear norm and the Frobenius norm of the temporal difference of . Importantly, this assumption is specific to the numerical study, which examines the method’s performance under controlled conditions. In practical applications in which direct knowledge of the true values is unavailable, as in the experimental study (see Sec. 4.5), empirical parameter sweeps can be used to tune regularization parameters for dynamic image reconstruction. Specifically, in this study, the solution to the dynamic image reconstruction problem is defined as For each case, four different values of the parameter were explored: . The parameter yielding the best average normalized squared error (nSE) over frames was selected.54 The nSE for each frame was computed as Specifically, the regularization parameter value that yielded the optimal results was when four tomographic measurements per frame were available and for all other cases. 4.5.Experimental StudyAn open-ended dynamic flow phantom was constructed by bending a silicone tube into a U-shaped structure, as depicted in Fig. 6. Figure 6(a) provides an illustration of the phantom within the imaging chamber, and Fig. 6(b) shows an actual photograph of the physical phantom. The inner diameter of the silicone tube was 0.0635 cm, and the outer diameter was 0.1194 cm. The dynamic phantom, positioned at the center of the imaging system, was illuminated by a laser pulse with a wavelength of 770 nm and an energy of 100 mJ (before entering the fiber optic light delivery unit). PACT data were acquired with the TriTom imaging system during the injection of a photoacoustic-fluorescent contrast agent (PAtIR55) through one end of the tube. Over a span of 36 s, the dynamic phantom underwent scanning, which resulted in 360 imaging frames in total, with 10 imaging frames per second and a single tomographic measurement per imaging frame. Two-dimensional fluorescence images were concurrently gathered, which serve as reference for the time evolution of the contrast agent concentration within the tube. Based on the region illuminated by laser, the reconstruction volume was set to a region of located at the center of the imaging system. The voxel size was set to , resulting in spatial voxels. The maximum allowed rank, , was set to 40, and the number of subsets, , was set to 18. For the stopping criterion of Algorithm 1, the threshold was set to . The step-size, , was tuned empirically to ensure convergence. To calibrate the speed of sound, multiple static estimates of the dynamic object were reconstructed using the universal back-projection algorithm56 assuming speed of sound values in the range of 1480 to , with increments. The speed of sound value of resulted in the best visual appearance and was selected to perform the dynamic image reconstruction. To choose the regularization parameters, nine dynamic image reconstructions were performed for all possible combinations of and . Ultimately, the spatiotemporal object estimate that most closely captured the observed dynamic changes in the reference fluorescence images, as determined through visual examination, was selected. The corresponding regularization parameters were and . 5.Results5.1.Inverse Crime Validation Study ResultsIn this inverse crime study, the step sizes for were empirically tuned to ensure convergence, resulting in values of , respectively. The algorithm was allowed to run for 2500 iterations in each case to ensure convergence of the solution up to numerical precision. For the purpose of this validation study, Algorithm 1 was modified to re-evaluate the total data fidelity term , which accumulates the contributions from all time frames, at each iteration. The results of the validation study are depicted in Fig. 7. Figure 7(a) exhibits the data fidelity versus the iteration count, and Fig. 7(b) presents the average nSE versus the iteration count. For all values of , a significant decrease of approximately 11 and 13 orders of magnitude can be observed in the data fidelity term and average nSE, respectively. This indicates that the proposed method using momentum-acceleration in combination with OS (case ), although lacking a theoretical guarantee of convergence, can achieve (up to machine precision) to the same object estimate produced by the momentum-accelerated PGD without the subsampling method (). It is also evident that a larger number of subsets improves the convergence speed, particularly in the early iterations. This also translate into a possibly faster time to solution as the cost of each iteration is dominated by the evaluation of the imaging operator, whereas the time spent in performing the truncated SVD factorization using the randomized method is negligible. In the numerical studies presented here, the computational time required per iteration was approximately 15 min on a workstation (AMD EPYC 7702P 64-Core processor, 32 GB RAM, one Nvidia Geforce RTX 2080 graphic processing unit), independent of the number of subsets used. In summary, this validation study demonstrates the correct implementation, efficiency, and robust convergence, despite a lack of theoretical guarantees, of the proposed method. 5.2.Numerical Phantom Study Results5.2.1.Sensitivity to the number of tomographic measurements per imaging frameFigure 8 displays a selection of MIP images depicting the numerical phantom and the spatiotemporal estimates from simulated data with varying numbers of tomographic measurements per imaging frame. A video featuring the spatiotemporal evolution of the numerical phantom and its corresponding estimates is available as Video 1. Upon careful inspection, it is apparent that increasing the number of tomographic measurements per imaging frame leads to more accurate estimates of the dynamic object. For instance, object estimates reconstructed from data with two and four measurements per imaging frame correctly capture the increase in intensity of blob-4 after frame-301. However, this intensity change is less prominent in the dynamic estimate reconstructed from data with only one tomographic measurement per imaging frame, which also exhibits artifacts around blob-4. In addition, the estimate reconstructed using four measurements per imaging frame better captures the intensity difference between the two ends of blob-1 in frame-1. These observations are more evident in Fig. 9, which shows the TACs at each ellipsoidal blob’s center in both the numerical phantom and the spatiotemporal estimates from data with different numbers of tomographic measurements per imaging frame. Notably, the TACs for the reconstruction from four tomographic measurements per imaging frame closely align with those of the numerical phantom. The nSE versus the imaging frame plot for the reconstructions with varying numbers of tomographic measurements per imaging frame is depicted in Fig. 11(a). Both the TACs and nSE versus frame plots confirm the observation that a higher number of tomographic measurements per imaging frame leads to more accurate spatiotemporal reconstructions. 5.2.2.Sensitivity to measurement noiseFigure 10 displays the TACs at each ellipsoidal blob’s center in both the numerical phantom and its spatiotemporal estimates from data with varying noise levels, when the number of tomographic measurements per frame is kept at two. Notably, the estimated TACs remain similar across all noise levels, highlighting the algorithm’s robustness in capturing dynamic changes even in the presence of increased noise. Figure 11(b) shows the nSE versus imaging frame plot for estimates reconstructed from data with varying noise levels. As anticipated, the nSE exhibits an upward trend with increasing noise levels. The numerical phantom study results underscore the effectiveness of the proposed method in accurately estimating dynamic changes, when limited tomographic measurements per frame are available. The study reveals that the ill-posed nature of the dynamic reconstruction problem diminishes significantly as the number of tomographic measurements per imaging frame increases. Notably, the method’s robustness in faithfully capturing the temporal dynamics of the object remains evident even in the presence of increased noise levels. These findings collectively underscore the algorithm’s reliability and potential significance in addressing challenges associated with dynamic imaging scenarios. 5.3.Experimental Study ResultsFigure 12 displays dynamic PACT images reconstructed from experimental data (bottom row), along with the corresponding reference fluorescence images (top row). The fluorescence images were processed to suppress the background, and the contrast agent is highlighted in yellow. This enhancement was accomplished through manual segmentation of the tube and contrast agent from the raw fluorescence images. The spatiotemporal PACT object estimates were visualized using ParaView.57 To ensure a qualitatively close alignment of the field-of-view between the ParaView visualization of the spatiotemporal PACT object estimate and the 2D fluorescence images, the view angle in ParaView was adjusted manually at each frame; however, a slight misalignment remains. A video featuring the raw images, visually enhanced images through the segmentation of the tube and contrast agent, and ParaView visualizations of the spatiotemporal PACT object estimate is provided in Video 2. In particular, when examining frames 30 to 120, one can readily observe the accurate recovery of dynamic contrast flow within the tube. A similar trend is evident in frames 210 to 300, with the lower right section of frame-240 showing the presence of the contrast agent, albeit with reduced contrast. This reduced contrast is also noticeable in specific frames, such as 330 and 360, possibly due to light obstruction by other tube segments. Nevertheless, the overall effective recovery of dynamic flow remains apparent, as convincingly demonstrated in Video 2. This highlights the effectiveness and practicality of the proposed method for experimental data beyond simulated measurements, even when only one tomographic measurement per imaging frame is available. Furthermore, it is worth noting that the spatiotemporal reconstruction offers a frame rate that is equal to the laser pulse rate [10 frames per second (FPS) for the instrument used in the phantom studies]. By contrast, an FBFIR technique would provide a frame rate of 1/36 FPS, which is equal to the reciprocal of the time required for a complete tomographic measurement. This underscores the reconstruction method’s potential for monitoring dynamic physiological processes that demand an enhanced frame rate. This can enhance the value of 3D PACT systems that involve rotating measurement gantries for preclinical research, enabling STIR from limited tomographic measurements per imaging frame. 6.Conclusion and DiscussionVolumetric PACT imagers commonly adopt a sequential data acquisition strategy utilizing rotating measurement gantries because of their cost-effectiveness and simplified hardware.2,19–21 However, the relatively slow data acquisition times associated with this approach pose significant challenges for dynamic imaging. Consequently, previous studies on dynamic PACT have predominantly focused on 2D dynamic imaging scenarios, leveraging the rapid data acquisition and computationally less demanding image reconstruction calculations of 2D systems. Recognizing that many commercially available volumetric imagers employ sequential data acquisition strategies with rotating gantries, it becomes imperative to advance dynamic image reconstruction techniques tailored to these widely adopted configurations. Addressing this need will enhance the versatility and practical utility of volumetric PACT imagers, allowing for more effective imaging of dynamic processes. This study presented 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 numerical studies demonstrated that the proposed method is robust against increasing noise levels. They also demonstrated that, as expected, increasing the number of tomographic measurements per frame improves the reconstruction accuracy. Therefore, it is desirable and beneficial for dynamic imaging to design acquisition systems that, like the University of Twente breast imager,52 utilize multiple acoustic probes. The experimental study demonstrated the LRME-STIR method’s ability to reconstruct the flow of a contrast agent at a frame rate of 10 FPS, 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 employed an accelerated PGD method combined with an OS approach, distinguishing it from previously proposed LRME-STIR methods.29,35,36,40 It is important to note that previous studies have demonstrated that the combination of the OS approach with momentum schemes can lead to convergence instability, particularly as the number of subsets increases.47,48,58,59 This instability arises due to error accumulation in the momentum term, as the “subset balance” approximation does not perfectly hold.47 The “subset balance” approximation implies that the cost function, , defined in Eq. (11) can be approximated with the OS-based cost function, , defined in Eq. (14), i.e., . When a small number of subsets is employed, the “subset balance” approximation is more accurate; thus, the method exhibits more stable convergence characteristics; however, the acceleration from the ordered subset approach is limited in such cases.47 To enhance stability, various strategies have been introduced in the literature. For instance, the relaxation of the momentum term has been proposed to mitigate convergence issues.47 Another approach involves a monotonic adaptation of the FISTA acceleration, which utilizes objective function values to adaptively determine the subsequent update.48 Additionally, an adaptive restarting method that resets the momentum terms based on a stability metric was introduced, further promoting convergence reliability.59 Although combining OS and momentum schemes does not provide theoretical guarantees of convergence, in the numerical and experimental studies, the proposed approach demonstrated stable convergence behavior, given the proper selection of the step size and number of subsets. Future refinements of the proposed method may explore the different strategies from the literature to obtain a provably convergent method.48,58,59 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.37,60 These evaluations might encompass both in-vivo and in-silico experiments, further advancing the understanding and application of the proposed method. DisclosuresThe authors S. Ermilov, W. Thompson, and M. Anastasio disclose ownership interests in Photosound Technologies, Inc. The other authors have no relevant financial interests and no potential conflicts of interest to disclose. Code and Data AvailabilityData and code are available upon request. Please contact Dr. Villa at uvilla@oden.utexas.edu. AcknowledgmentsThe research reported in this publication was supported in part by the Office of The Director, National Institutes of Health (NIH; grant nos. R44OD023029 and R01EB031585). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. ReferencesX. Wang et al.,
“Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,”
Nat. Biotechnol., 21
(7), 803
–806 https://doi.org/10.1038/nbt839 NABIF9 1087-0156
(2003).
Google Scholar
H.-P. Brecht et al.,
“Whole-body three-dimensional optoacoustic tomography system for small animals,”
J. Biomed. Opt., 14
(6), 064007 https://doi.org/10.1117/1.3259361 JBOPFO 1083-3668
(2009).
Google Scholar
S. Yang et al.,
“Functional imaging of cerebrovascular activities in small animals using high-resolution photoacoustic tomography,”
Med. Phys., 34
(8), 3294
–3301 https://doi.org/10.1118/1.2757088 MPHYA6 0094-2405
(2007).
Google Scholar
G. C. Kagadis et al.,
“In vivo small animal imaging: current status and future prospects,”
Med. Phys., 37
(12), 6421
–6442 https://doi.org/10.1118/1.3515456 MPHYA6 0094-2405
(2010).
Google Scholar
G. Loudos, G. C. Kagadis and D. Psimadas,
“Current status and future perspectives of in vivo small animal imaging using radiolabeled nanoparticles,”
Eur. J. Radiol., 78
(2), 287
–295 https://doi.org/10.1016/j.ejrad.2010.06.025 EJRADR 0720-048X
(2011).
Google Scholar
B. L. Franc et al.,
“Small-animal SPECT and SPECT/CT: important tools for preclinical investigation,”
J. Nucl. Med., 49
(10), 1651
–1663 https://doi.org/10.2967/jnumed.108.055442 JNMEAQ 0161-5505
(2008).
Google Scholar
P. Carmeliet,
“Angiogenesis in life, disease and medicine,”
Nature, 438
(7070), 932
–936 https://doi.org/10.1038/nature04478
(2005).
Google Scholar
P. Carmeliet and R. K. Jain,
“Angiogenesis in cancer and other diseases,”
Nature, 407
(6801), 249
–257 https://doi.org/10.1038/35025220
(2000).
Google Scholar
J. Xia and L. V. Wang,
“Small-animal whole-body photoacoustic tomography: a review,”
IEEE Trans. Biomed. Eng., 61
(5), 1380
–1389 https://doi.org/10.1109/TBME.2013.2283507 IEBEAX 0018-9294
(2014).
Google Scholar
L. Li et al.,
“Single-impulse panoramic photoacoustic computed tomography of small-animal whole-body dynamics at high spatiotemporal resolution,”
Nat. Biomed. Eng., 1
(5), 0071 https://doi.org/10.1038/s41551-017-0071
(2017).
Google Scholar
P. K. Upputuri and M. Pramanik,
“Dynamic in vivo imaging of small animal brain using pulsed laser diode-based photoacoustic tomography system,”
J. Biomed. Opt., 22
(9), 090501 https://doi.org/10.1117/1.JBO.22.9.090501 JBOPFO 1083-3668
(2017).
Google Scholar
S. Shi et al.,
“Thermosensitive biodegradable copper sulfide nanoparticles for real-time multispectral optoacoustic tomography,”
ACS Appl. Biomater., 2
(8), 3203
–3211 https://doi.org/10.1021/acsabm.9b00133
(2019).
Google Scholar
S. Manohar and M. Dantuma,
“Current and future trends in photoacoustic breast imaging,”
Photoacoustics, 16 100134 https://doi.org/10.1016/j.pacs.2019.04.004
(2019).
Google Scholar
L. Lin et al.,
“Single-breath-hold photoacoustic computed tomography of the breast,”
Nat. Commun., 9
(1), 2352 https://doi.org/10.1038/s41467-018-04576-z NCAOBW 2041-1723
(2018).
Google Scholar
J. Gamelin et al.,
“A real-time photoacoustic tomography system for small animals,”
Opt. Express, 17
(13), 10489
–10498 https://doi.org/10.1364/OE.17.010489 OPEXFF 1094-4087
(2009).
Google Scholar
W. He et al.,
“Plasmonic titanium nitride nanoparticles for in vivo photoacoustic tomography imaging and photothermal cancer therapy,”
Biomaterials, 132 37
–47 https://doi.org/10.1016/j.biomaterials.2017.04.007 BIMADU 0142-9612
(2017).
Google Scholar
P. K. Upputuri et al.,
“A high-performance compact photoacoustic tomography system for in vivo small-animal brain imaging,”
J. Visual. Exp.,
(124), e55811 https://doi.org/10.3791/55811
(2017).
Google Scholar
T. Shan et al.,
“In-vivo hemodynamic imaging of acute prenatal ethanol exposure in fetal brain by photoacoustic tomography,”
J. Biophotonics, 13
(5), e201960161 https://doi.org/10.1002/jbio.201960161
(2020).
Google Scholar
H. P. Brecht et al.,
“A 3D imaging system integrating photoacoustic and fluorescence orthogonal projections for anatomical, functional and molecular assessment of rodent models,”
Proc. SPIE, 10494 1049411 https://doi.org/10.1117/12.2290880 PSISDG 0277-786X
(2018).
Google Scholar
L. Lin et al.,
“High-speed three-dimensional photoacoustic computed tomography for preclinical research and clinical translation,”
Nat. Commun., 12
(1), 882 https://doi.org/10.1038/s41467-021-21232-1 NCAOBW 2041-1723
(2021).
Google Scholar
S. A. Ermilov et al.,
“Laser optoacoustic imaging system for detection of breast cancer,”
J. Biomed. Opt., 14
(2), 024007 https://doi.org/10.1117/1.3086616 JBOPFO 1083-3668
(2009).
Google Scholar
L. Xiang et al.,
“4-D photoacoustic tomography,”
Sci. Rep., 3
(1), 1
–8 https://doi.org/10.1038/srep01113 SRCEC3 2045-2322
(2013).
Google Scholar
A. Hauptmann et al.,
“Model-based learning for accelerated, limited-view 3-D photoacoustic tomography,”
IEEE Trans. Med. Imaging, 37
(6), 1382
–1393 https://doi.org/10.1109/TMI.2018.2820382 ITMID4 0278-0062
(2018).
Google Scholar
J. Schwab, S. Antholzer and M. Haltmeier,
“Learned backprojection for sparse and limited view photoacoustic tomography,”
Proc. SPIE, 10878 1087837 https://doi.org/10.1117/12.2508438 PSISDG 0277-786X
(2019).
Google Scholar
J. Frikel and E. T. Quinto,
“Artifacts in incomplete data tomography with applications to photoacoustic tomography and sonar,”
SIAM J. Appl. Math., 75
(2), 703
–725 https://doi.org/10.1137/140977709 SMJMAP 0036-1399
(2015).
Google Scholar
J.-F. Cai et al.,
“Cine cone beam CT reconstruction using low-rank matrix factorization: algorithm and a proof-of-principle study,”
IEEE Trans. Med. Imaging, 33
(8), 1581
–1591 https://doi.org/10.1109/TMI.2014.2319055 ITMID4 0278-0062
(2014).
Google Scholar
M. Wernick, E. Infusino and M. Milosevic,
“Fast spatio-temporal image reconstruction for dynamic PET,”
IEEE Trans. Med. Imaging, 18
(3), 185
–195 https://doi.org/10.1109/42.764885 ITMID4 0278-0062
(1999).
Google Scholar
Q. Ding, M. Burger and X. Zhang,
“Dynamic SPECT reconstruction with temporal edge correlation,”
Inverse Probl., 34
(1), 014005 https://doi.org/10.1088/1361-6420/aa9a94 INPEEY 0266-5611
(2017).
Google Scholar
J. P. Haldar and Z.-P. Liang,
“Spatiotemporal imaging with partially separable functions: a matrix recovery approach,”
in IEEE Int. Symp. Biomed. Imaging: From Nano to Macro,
716
–719
(2010). https://doi.org/10.1109/ISBI.2010.5490076 Google Scholar
K. Wang et al.,
“Fast spatiotemporal image reconstruction based on low-rank matrix estimation for dynamic photoacoustic computed tomography,”
J. Biomed. Opt., 19
(5), 056007 https://doi.org/10.1117/1.JBO.19.5.056007 JBOPFO 1083-3668
(2014).
Google Scholar
S. Arridge et al.,
“Accelerated high-resolution photoacoustic tomography via compressed sensing,”
Phys. Med. Biol., 61
(24), 8908 https://doi.org/10.1088/1361-6560/61/24/8908 PHMBA7 0031-9155
(2016).
Google Scholar
A. Özbek, X. L. Deán-Ben and D. Razansky,
“Compressed optoacoustic sensing of volumetric cardiac motion,”
IEEE Trans. Med. Imaging, 39
(10), 3250
–3255 https://doi.org/10.1109/TMI.2020.2985134 ITMID4 0278-0062
(2020).
Google Scholar
F. Lucka et al.,
“Enhancing compressed sensing 4D photoacoustic tomography by simultaneous motion estimation,”
SIAM J. Imaging Sci., 11
(4), 2224
–2253 https://doi.org/10.1137/18M1170066
(2018).
Google Scholar
K. Wang et al.,
“Accelerating image reconstruction in three-dimensional optoacoustic tomography on graphics processing units,”
Med. Phys., 40
(2), 023301 https://doi.org/10.1118/1.4774361 MPHYA6 0094-2405
(2013).
Google Scholar
Z.-P. Liang,
“Spatiotemporal imaging with partially separable functions,”
in 4th IEEE Int. Symp. Biomed. Imaging: From Nano to Macro,
988
–991
(2007). https://doi.org/10.1109/ISBI.2007.357020 Google Scholar
C. Brinegar et al.,
“Real-time cardiac MRI using prior spatial-spectral information,”
in Annu. Int. Conf. IEEE Eng. in Med. and Biol. Soc.,
4383
–4386
(2009). https://doi.org/10.1109/IEMBS.2009.5333482 Google Scholar
R. M. Cam et al.,
“Dynamic image reconstruction to monitor tumor vascular perfusion in small animals using 3D photoacoustic computed-tomography imagers with rotating gantries,”
Proc. SPIE, 12379 123790F https://doi.org/10.1117/12.2647663 PSISDG 0277-786X
(2023).
Google Scholar
R. M. Cam et al.,
“Low-rank matrix estimation-based spatiotemporal image reconstruction from few tomographic measurements per frame for dynamic photoacoustic computed tomography,”
Proc. SPIE, 12463 124630R https://doi.org/10.1117/1.JBO.19.5.056007 PSISDG 0277-786X
(2023).
Google Scholar
L. Lozenski, M. A. Anastasio and U. Villa,
“A memory-efficient self-supervised dynamic image reconstruction method using neural fields,”
IEEE Trans. Comput. Imaging, 8 879
–892 https://doi.org/10.1109/TCI.2022.3208511
(2022).
Google Scholar
S. G. Lingala et al.,
“Accelerated dynamic MRI exploiting sparsity and low-rank structure: kt SLR,”
IEEE Trans. Med. Imaging, 30
(5), 1042
–1054 https://doi.org/10.1109/TMI.2010.2100850 ITMID4 0278-0062
(2011).
Google Scholar
S. Gu et al.,
“Weighted nuclear norm minimization with application to image denoising,”
in Proc. IEEE Conf. Comput. Vis. and Pattern Recognit.,
2862
–2869
(2014). https://doi.org/10.1109/CVPR.2014.366 Google Scholar
R. Tibshirani et al.,
“Proximal gradient descent and acceleration,”
(2010). http://www.stat.cmu.edu/ryantibs/convexopt-S15/lectures/08-prox-grad.pdf Google Scholar
P. L. Combettes, J.-C. Pesquet,
“Proximal splitting methods in signal processing,”
Fixed-Point Algorithms for Inverse Problems in Science and Engineering, 185
–212 Springer(
(2011). Google Scholar
A. Beck and M. Teboulle,
“A fast iterative shrinkage-thresholding algorithm for linear inverse problems,”
SIAM J. Imaging Sci., 2
(1), 183
–202 https://doi.org/10.1137/080716542
(2009).
Google Scholar
Y. Nesterov,
“Gradient methods for minimizing composite functions,”
Math. Program., 140
(1), 125
–161 https://doi.org/10.1007/s10107-012-0629-5 MHPGA4 1436-4646
(2013).
Google Scholar
H. M. Hudson and R. S. Larkin,
“Accelerated image reconstruction using ordered subsets of projection data,”
IEEE Trans. Med. Imaging, 13
(4), 601
–609 https://doi.org/10.1109/42.363108 ITMID4 0278-0062
(1994).
Google Scholar
D. Kim, S. Ramani and J. A. Fessler,
“Combining ordered subsets and momentum for accelerated X-ray CT image reconstruction,”
IEEE Trans. Med. Imaging, 34
(1), 167
–178 https://doi.org/10.1109/TMI.2014.2350962 ITMID4 0278-0062
(2014).
Google Scholar
V. Haase et al.,
“An improved combination of ordered subsets and momentum for fast model-based iterative CT reconstruction,”
Proc. SPIE, 11312 113121N https://doi.org/10.1117/12.2549733 PSISDG 0277-786X
(2020).
Google Scholar
J.-F. Cai, E. J. Candès and Z. Shen,
“A singular value thresholding algorithm for matrix completion,”
SIAM J. Optim., 20
(4), 1956
–1982 https://doi.org/10.1137/080738970 SJOPE8 1095-7189
(2010).
Google Scholar
N. Halko, P.-G. Martinsson and J. A. Tropp,
“Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions,”
SIAM Rev., 53
(2), 217
–288 https://doi.org/10.1137/090771806 SIREAD 0036-1445
(2011).).
Google Scholar
W. R. Thompson et al.,
“Characterizing a photoacoustic and fluorescence imaging platform for preclinical murine longitudinal studies,”
J. Biomed. Opt., 28
(3), 036001 https://doi.org/10.1117/1.JBO.28.3.036001 JBOPFO 1083-3668
(2023).
Google Scholar
S. M. Schoustra et al.,
“Twente Photoacoustic Mammoscope 2: system overview and three-dimensional vascular network images in healthy breasts,”
J. Biomed. Opt., 24
(12), 121909 https://doi.org/10.1117/1.JBO.24.12.121909 JBOPFO 1083-3668
(2019).
Google Scholar
K. Ito, B. Jin and T. Takeuchi,
“Multi-parameter Tikhonov regularization—an augmented approach,”
Chin. Ann. Math. Ser. B, 35
(3), 383
–398 https://doi.org/10.1007/s11401-014-0835-y
(2014).
Google Scholar
S. Pereverzev and E. Schock,
“On the adaptive selection of the parameter in regularization of ill-posed problems,”
SIAM J. Numer. Anal., 43
(5), 2060
–2076 https://doi.org/10.1137/S0036142903433819
(2005).
Google Scholar
W. Thompson et al.,
“A preclinical small animal imaging platform combining multi-angle photoacoustic and fluorescence projections into co-registered 3D maps,”
Proc. SPIE, 11240 112400L https://doi.org/10.1117/12.2549088 PSISDG 0277-786X
(2020).
Google Scholar
M. Xu and L. V. Wang,
“Universal back-projection algorithm for photoacoustic computed tomography,”
Phys. Rev. E, 71
(1), 016706 https://doi.org/10.1103/PhysRevE.71.016706
(2005).
Google Scholar
J. Ahrens, B. Geveci and C. Law, ParaView: An End-User Tool for Large Data Visualization, 717
–731 Elsevier Butterworth-Heinemann, Burlington, MA, USA
(2005). Google Scholar
D. Kim and J. A. Fessler,
“Ordered subsets acceleration using relaxed momentum for x-ray CT image reconstruction,”
in IEEE Nucl. Sci. Symp. and Med. Imaging Conf. (2013 NSS/MIC),
1
–5
(2013). https://doi.org/10.1109/NSSMIC.2013.6829298 Google Scholar
A. Sisniega et al.,
“Accelerated 3D image reconstruction with a morphological pyramid and noise-power convergence criterion,”
Phys. Med. Biol., 66
(5), 055012 https://doi.org/10.1088/1361-6560/abde97 PHMBA7 0031-9155
(2021).
Google Scholar
C. W. Hupple et al.,
“A light-fluence-independent method for the quantitative analysis of dynamic contrast-enhanced multispectral optoacoustic tomography (DCE MSOT),”
Photoacoustics, 10 54
–64 https://doi.org/10.1016/j.pacs.2018.04.003
(2018).
Google Scholar
BiographyRefik Mert Cam is a PhD candidate in electrical and computer engineering at the University of Illinois Urbana-Champaign. He obtained his BS degree in electrical and electronics engineering from the Middle East Technical University, Turkey, in 2020. His research interests include inverse problems, medical image reconstruction, and deep learning for medical imaging. He is a student member of SPIE and 2023–2024 Mavis Future Faculty Fellow at Grainger College of Engineering of the University of Illinois Urbana-Champaign. Chao Wang received her BS degree from Jilin University in China. She obtained her PhD in computational mathematics from Peking University in 2019. She worked as a postdoc at the University of Illinois at Urbana-Champaign in 2020 and as a research fellow at the National University of Singapore from 2021 to 2023. Her research involves several aspects of computational mathematics, including PDE-constrained inverse problems, medical imaging, generative neural networks in operator learning, and modeling of proton therapy. Weylan Thompson received his BS degree in physics from McGill University in 2015. He joined PhotoSound Technologies, Inc., in 2018 and has since been involved in developing and characterizing photoacoustic imaging systems and data acquisition devices. Sergey A. Ermilov holds his master of biomedical engineering degree from Bauman Moscow State Technical University, Russia, and his PhD in bioengineering from Rice University, Houston, Texas, United States. His scientific expertise and interests are in the field of biomedical photoacoustics and include development of theoretical aspects, system modeling, development of multimodality instruments, and data processing. He is an author of 24 peer-reviewed publications and 19 patents related to the field of photoacoustic imaging. Mark A. Anastasio is a Donald Biggar Willett professor in engineering and the head of the Department of Bioengineering, University of Illinois at Urbana–Champaign, Urbana, Illinois, United States. He received his PhD from the University of Chicago, Chicago, Illinois, United States, in 2001. His research broadly addresses computational image science, inverse problems in imaging, and machine learning for imaging applications. He has contributed broadly to emerging biomedical imaging technologies that include photoacoustic computed tomography, ultrasound computed tomography, and x-ray phase-contrast imaging. His work has been supported by numerous NIH grants, and he served for 2 years as the chair of the NIH EITA Study Section. He is a fellow of SPIE, the American Institute for Medical and Biological Engineering, the International Academy of Medical and Biological Engineering, and the Institute of Electrical and Electronics Engineers. Umberto Villa is a research scientist at Oden Institute for Computational Engineering and Science, University of Texas at Austin, Texas, United States. He received his BS and MS degrees in mathematical engineering from Politecnico di Milano, Milan, Italy, in 2005 and 2007, respectively, and his PhD in mathematics from Emory University, Atlanta, Georgia, United States, in 2012. His research interests lie in the computational and mathematical aspects of large-scale inverse problems, imaging science, and uncertainty quantification. |
Image restoration
Photoacoustic tomography
Tomography
Imaging systems
Transducers
Data acquisition
Biomedical optics