Accelerating image reconstruction in three-dimensional optoacoustic tomography on graphics processing units

Ultrafast filtered back-projection for photoacoustic computed tomography

The filtered back-projection (FBP) algorithm is widely used in photoacoustic computed tomography (PACT) for image reconstruction due to its simplicity and efficiency. Yet, the real-time processing of high-speed PACT data remains challenging for regular FBP implementations. To enhance the reconstruction efficiency of the FBP algorithm, researchers have developed FBP implementations based on the graphics processing units (GPUs). However, existing GPU-accelerated FBP algorithms either sacrifice accuracy for efficiency or are still inefficient for high-speed, real-time PACT imaging. Herein, we report an ultrafast GPU acceleration-based FBP implementation for PACT image reconstruction without sacrificing accuracy. Firstly, the computation complexity of the filtering part of the FBP algorithm is significantly simplified with a pre-computed filtering matrix to enhance filtering efficiency. Secondly, the computation efficiency of the back-projection part of the FBP algorithm is dramatically increased through parallel programming and GPU acceleration. As a result, the proposed FBP implementation takes only 0.38 ms to reconstruct a two-dimensional image of 512 × 512 pixels, which is 439 times faster than regular FBP implementations. Numerical and experimental results show that the proposed FBP implementation outperforms existing GPU-based FBP implementations in reconstruction accuracy and computation efficiency. To the best of our knowledge, this is the fastest implementation of the FBP algorithm ever reported in PACT. This work is expected to provide an ultrafast and accurate image reconstruction solution for high-speed, real-time PACT imaging.

Spatiotemporal image reconstruction to enable high-frame-rate dynamic photoacoustic tomography with rotating-gantry volumetric imagers

Significance

Dynamic photoacoustic computed tomography (PACT) is a valuable imaging technique for monitoring physiological processes. However, current dynamic PACT imaging techniques are often limited to two-dimensional spatial imaging. Although volumetric PACT imagers are commercially available, these systems typically employ a rotating measurement gantry in which the tomographic data are sequentially acquired as opposed to being acquired simultaneously at all views. Because the dynamic object varies during the data-acquisition process, the sequential data-acquisition process poses substantial challenges to image reconstruction associated with data incompleteness. The proposed image reconstruction method is highly significant in that it will address these challenges and enable volumetric dynamic PACT imaging with existing preclinical imagers.

Aim

The aim of this study is to develop a spatiotemporal image reconstruction (STIR) method for dynamic PACT that can be applied to commercially available volumetric PACT imagers that employ a sequential scanning strategy. The proposed reconstruction method aims to overcome the challenges caused by the limited number of tomographic measurements acquired per frame.

Approach

A low-rank matrix estimation-based STIR (LRME-STIR) method is proposed to enable dynamic volumetric PACT. The LRME-STIR method leverages the spatiotemporal redundancies in the dynamic object to accurately reconstruct a four-dimensional (4D) spatiotemporal image.

Results

The conducted numerical studies substantiate the LRME-STIR method's efficacy in reconstructing 4D dynamic images from tomographic measurements acquired with a rotating measurement gantry. The experimental study demonstrates the method's ability to faithfully recover the flow of a contrast agent with a frame rate of 10 frames per second, even when only a single tomographic measurement per frame is available.

Conclusions

The proposed LRME-STIR method offers a promising solution to the challenges faced by enabling 4D dynamic imaging using commercially available volumetric PACT imagers. By enabling accurate STIRs, this method has the potential to significantly advance preclinical research and facilitate the monitoring of critical physiological biomarkers.

Model-Based 3-D X-Ray Induced Acoustic Computerized Tomography

X-ray-induced acoustic (XA) computerized tomography (XACT) is an evolving imaging technique that aims to reconstruct the X-ray energy deposition from XA measurements. Main challenges in XACT are the poor signal-to-noise ratio and limited field-of-view, which cause artifacts in the images. We demonstrate the efficacy of model-based (MB) algorithms for three-dimensional XACT and compare with the traditional algorithms. The MB algorithm is based on iterative, matrix-free approach for regularized-least-squares minimization corresponding to XACT. The matrix-free-LSQR (MF-LSQR) and the non-iterative model-backprojection (MBP) reconstructions were evaluated and compared with universal backprojection (UBP), time-reversal (TR) and fast-Fourier transform (FFT)-based reconstructions for numerical and experimental XACT datasets. The results demonstrate the capability of MF-LSQR algorithm to reduce noisy artifacts thus yielding better reconstructions. MBP and MF-LSQR algorithms perform particularly well with the experimental XACT dataset, where noise in signals significantly affects the reconstruction of the target in UBP and FFT-based reconstructions. The TR reconstruction for experimental XACT are comparable to MF-LSQR, but takes thrice as much time and filters the frequency components greater than maximum frequency supported by the grid, resulting loss of resolution. The MB algorithms are able to overcome the challenges in XACT and hence are vital for the clinical translation of XACT.

An Investigation of Signal Preprocessing for Photoacoustic Tomography

Photoacoustic tomography (PAT) is increasingly being used for high-resolution biological imaging at depth. Signal-to-noise ratios and resolution are the main factors that determine image quality. Various reconstruction algorithms have been proposed and applied to reduce noise and enhance resolution, but the efficacy of signal preprocessing methods which also affect image quality, are seldom discussed. We, therefore, compared common preprocessing techniques, namely bandpass filters, wavelet denoising, empirical mode decomposition, and singular value decomposition. Each was compared with and without accounting for sensor directivity. The denoising performance was evaluated with the contrast-to-noise ratio (CNR), and the resolution was calculated as the full width at half maximum (FWHM) in both the lateral and axial directions. In the phantom experiment, counting in directivity was found to significantly reduce noise, outperforming other methods. Irrespective of directivity, the best performing methods for denoising were bandpass, unfiltered, SVD, wavelet, and EMD, in that order. Only bandpass filtering consistently yielded improvements. Significant improvements in the lateral resolution were observed using directivity in two out of three acquisitions. This study investigated the advantages and disadvantages of different preprocessing methods and may help to determine better practices in PAT reconstruction.

A Review of Emerging Electromagnetic-Acoustic Sensing Techniques for Healthcare Monitoring

Conventional electromagnetic (EM) sensing techniques such as radar and LiDAR are widely used for remote sensing, vehicle applications, weather monitoring, and clinical monitoring. Acoustic techniques such as sonar and ultrasound sensors are also used for consumer applications, such as ranging and in vivo medical/healthcare applications. It has been of long-term interest to doctors and clinical practitioners to realize continuous healthcare monitoring in hospitals and/or homes. Physiological and biopotential signals in real-time serve as important health indicators to predict and prevent serious illness. Emerging electromagnetic-acoustic (EMA) sensing techniques synergistically combine the merits of EM sensing with acoustic imaging to achieve comprehensive detection of physiological and biopotential signals. Further, EMA enables complementary fusion sensing for challenging healthcare settings, such as real-world long-term monitoring of treatment effects at home or in remote environments. This article reviews various examples of EMA sensing instruments, including implementation, performance, and application from the perspectives of circuits to systems. The novel and significant applications to healthcare are discussed. Three types of EMA sensors are presented: (1) Chip-based radar sensors for health status monitoring, (2) Thermo-acoustic sensing instruments for biomedical applications, and (3) Photoacoustic (PA) sensing and imaging systems, including dedicated reconstruction algorithms were reviewed from time-domain, frequency-domain, time-reversal, and model-based solutions. The future of EMA techniques for continuous healthcare with enhanced accuracy supported by artificial intelligence (AI) is also presented.

Enhancing Finite Element-Based Photoacoustic Tomography by Localized Reconstruction Method

Iterative reconstruction algorithm based on finite element (FE) modeling is a powerful approach in photoacoustic tomography (PAT). However, an iterative inverse algorithm using conventional FE meshing of the entire imaging zone is computationally demanding, which hinders this powerful tool in applications where quick image acquisition and/or a large image matrix is needed. To address this challenge, parallel computing techniques are proposed and implemented in the field. Here, we present an alternative approach for 2D PAT, which locoregionally reconstructs the region of interest (ROI) instead of the full imaging zone. Our simulated and phantom experimental results demonstrate that this ROI reconstruction algorithm can produce almost the same image quality as the conventional full zone-based reconstruction algorithm; however, the computation time can be significantly reduced without any additional hardware cost by more than two orders of magnitude (100-fold). This algorithm is further applied and validated in an in vivo study. The major vessel structures in a rat’s brain can be imaged clearly using our ROI-based approach, coupled with a mesh of 11,801 nodes. This novel algorithm can also be parallelized using MPI or GPU acceleration techniques to further enhance the reconstruction performance of FE-based PAT.

Normalization of optical fluence distribution for three-dimensional functional optoacoustic tomography of the breast

SIGNIFICANCE

In three-dimensional (3D) functional optoacoustic tomography (OAT), wavelength-dependent optical attenuation and nonuniform incident optical fluence limit imaging depth and field of view and can hinder accurate estimation of functional quantities, such as the vascular blood oxygenation. These limitations hinder OAT of large objects, such as a human female breast.

AIM

We aim to develop a measurement-data-driven method for normalization of the optical fluence distribution and to investigate blood vasculature detectability and accuracy for estimating vascular blood oxygenation.

APPROACH

The proposed method is based on reasonable assumptions regarding breast anatomy and optical properties. The nonuniform incident optical fluence is estimated based on the illumination geometry in the OAT system, and the depth-dependent optical attenuation is approximated using Beer-Lambert law.

RESULTS

Numerical studies demonstrated that the proposed method significantly enhanced blood vessel detectability and improved estimation accuracy of the vascular blood oxygenation from multiwavelength OAT measurements, compared with direct application of spectral linear unmixing without optical fluence compensation. Experimental results showed that the proposed method revealed previously invisible structures in regions deeper than 15 mm and/or near the chest wall.

CONCLUSIONS

The proposed method provides a straightforward and computationally inexpensive approximation of wavelength-dependent effective optical attenuation and, thus, enables mitigation of the spectral coloring effect in functional 3D OAT imaging.

Denoising of pre-beamformed photoacoustic data using generative adversarial networks

We have trained generative adversarial networks (GANs) to mimic both the effect of temporal averaging and of singular value decomposition (SVD) denoising. This effectively removes noise and acquisition artifacts and improves signal-to-noise ratio (SNR) in both the radio-frequency (RF) data and in the corresponding photoacoustic reconstructions. The method allows a single frame acquisition instead of averaging multiple frames, reducing scan time and total laser dose significantly. We have tested this method on experimental data, and quantified the improvement over using either SVD denoising or frame averaging individually for both the RF data and the reconstructed images. We achieve a mean squared error (MSE) of 0.05%, structural similarity index measure (SSIM) of 0.78, and a feature similarity index measure (FSIM) of 0.85 compared to our ground-truth RF results. In the subsequent reconstructions using the denoised data we achieve a MSE of 0.05%, SSIM of 0.80, and a FSIM of 0.80 compared to our ground-truth reconstructions.

Fast Correction of “Finite Aperture Effect” in Photoacoustic Tomography Based on Spatial Impulse Response

Photoacoustic computed tomography (PACT) is a fast-developing imaging technique, which can provide structural and functional information in biological tissues with high-resolution beyond the depth of the optical diffusion limit. However, the most current PACT reconstruction method generally employs a point detector assumption, whereas in most PAT systems with circular or spherical scanning modes, the transducer is mostly flat and with a finite size. This model mismatch leads to a notable deterioration in the lateral direction in regions far from the rotation center, which is known as the “finite aperture effect”. In this work, we propose to compensate a novel Back-projection (BP) method based on the transducer’s spatial impulse response (SIR) for fast correction of the “finite aperture effect”. The SIR accounts for the waveform change of the transducer for an arbitrary point source due to the geometry of the detection surface. Simulation results showed that the proposed SIR-BP method can effectively improve the lateral resolution and signal to noise ratio (SNR) in the off-center regions. For a target 4.5 mm far from the rotation center, this new method improved the lateral resolution about five times along with a 7 dB increase in the SNR. Experimental results also showed that this SIR-BP method can well restore the image angular blur to recover small structures, as demonstrated by the imaging of leaf veins. This new method offers a valuable alternative to the conventional BP method, and can guide the design of PAT systems based on circular/spherical scan.

LightSpeed: A Compact, High-Speed Optical-Link-Based 3D Optoacoustic Imager

Wide-scale adoption of optoacoustic imaging in biology and medicine critically depends on availability of affordable scanners combining ease of operation with optimal imaging performance. Here we introduce LightSpeed: a low-cost real-time volumetric handheld optoacoustic imager based on a new compact software-defined ultrasound digital acquisition platform and a pulsed laser diode. It supports the simultaneous signal acquisition from up to 192 ultrasound channels and provides a hig-bandwidth direct optical link (2x 100G Ethernet) to the host-PC for ultra-high frame rate image acquisitions. We demonstrate use of the system for ultrafast (500Hz) 3D human angiography with a rapidly moving handheld probe. LightSpeed attained image quality comparable with a conventional optoacoustic imaging systems employing bulky acquisition electronics and a Q-switched pulsed laser. Our results thus pave the way towards a new generation of compact, affordable and high-performance optoacoustic scanners.

Single laser-shot super-resolution photoacoustic tomography with fast sparsity-based reconstruction

Recently, ℓ 1 -norm based reconstruction approaches have been used with linear array systems to improve photoacoustic resolution and demonstrate undersampled imaging when there is sufficient sparsity in some domain. However, such approaches have yet to beat the half-wavelength resolution limit. In this paper, the ability to beat the half-wavelength diffraction limit is demonstrated using a 5 MHz ring array photoacoustic tomography system and ℓ 1 -norm based reconstruction approaches. We used the array system to image wire targets at ≈ 2 - 3 cm depth in both intralipid scattering solution and water. The minimum observable separation was estimated as 70 ± 10 μ m , improving on the half-wavelength resolution limit of 145 μ m . This improvement was demonstrated even when using a random projection transform to reduce data by 99 % , enabling substantially faster reconstruction times. This is the first photoacoustic tomography approach capable of beating the half-wavelength resolution limit with a single laser shot.

Computationally Efficient Forward Operator for Photoacoustic Tomography Based on Coordinate Transformations

Photoacoustic tomography (PAT) is an imaging modality that utilizes the photoacoustic effect. In PAT, a photoacoustic image is computed from measured data by modeling ultrasound propagation in the imaged domain and solving an inverse problem utilizing a discrete forward operator. However, in realistic measurement geometries with several ultrasound transducers and relatively large imaging volume, an explicit formation and use of the forward operator can be computationally prohibitively expensive. In this work, we propose a transformation-based approach for efficient modeling of photoacoustic signals and reconstruction of photoacoustic images. In the approach, the forward operator is constructed for a reference ultrasound transducer and expanded into a general measurement geometry using transformations that map the formulated forward operator in local coordinates to the global coordinates of the measurement geometry. The inverse problem is solved using a Bayesian framework. The approach is evaluated with numerical simulations and experimental data. The results show that the proposed approach produces accurate 3-D photoacoustic images with a significantly reduced computational cost both in memory requirements and time. In the studied cases, depending on the computational factors, such as discretization, over the 30-fold reduction in memory consumption was achieved without a reduction in image quality compared to a conventional approach.

Approximate back-projection method for improving lateral resolution in circular-scanning-based photoacoustic tomography

PURPOSE

In circular-scanning-based photoacoustic tomography (PAT), the effect of finite transducer aperture has not been effectively resolved. The goal of this paper is to propose a practical reconstruction method that accounts for the finite transducer aperture to improve the lateral resolution.

METHODS

We for the first time propose to calculate the spatial-temporal response (STR) of the employed finite-sized transducer in a forward model, and then compensate the time delay and the directional sensitivity of the transducer in the framework of the back-projection method. Both simulation and phantom experiments were carried out to evaluate the lateral resolution improvement with the proposed method. The performance of this new method for imaging complicated targets was also assessed by calculating the mean image gradient.

RESULTS

Simulation results showed that with this new method the lateral resolution for off-center targets can be as good as that for the center targets. Phantom experimental results showed that this new method can improve the lateral resolution more than two times for a point target about 5 mm far from the rotation center. Phantom experimental results also showed that many blurred fine structures of a piece of leaf veins at the off-center regions were well restored with the new method, and the mean image gradient improved about 1.3 times.

CONCLUSION

The proposed new method can effectively account for the effect of finite transducer aperture for circular-scanning-based PAT in homogenous acoustic media. This new method also features its robustness and computational efficiency, so that it is a worthy replacement to the conventional back-projection algorithm in circular-scanning-based PAT. This new method can be of great importance to the design of circular-scanning or spherical-scanning-based PAT systems.

Fast 3-D Opto-Acoustic Simulation for Linear Array With Rectangular Elements

Simulation involves predicting the responses of a physical system. In this article, we simulate opto-acoustic signals generated in a 3-D volume due to the absorption of an optical pulse. A separable computational model is developed, which splits processing into two steps, permitting an order-of-magnitude improvement in computational efficiency over a nonseparable model. The simulated signals represent acoustic waves, measured by a probe with a linear transducer array, in a rotated and translated coordinate frame. Light is delivered by an optical source that moves with the probe's frame. Spatiotemporal impulse response for rectangular-element transducer geometry is derived using Green's function solution to the acoustic wave equation. The approach permits fast and accurate simulation for a probe with an arbitrary trajectory, which is useful for modeling freehand acquisitions. For a 3-D volume of n³ voxels, computation is accelerated by a factor of n . This may potentially have application in opto-acoustic imaging, where clinicians visualize structural and functional features of biological tissue for assessment of cancer and other diseases.

Cross-sectional photoacoustic tomography image reconstruction with a multi-curve integration model

BACKGROUND AND OBJECTIVE

In acoustic inversion of photoacoustic tomography (PAT), an imaging model that precisely describes both the ultrasonic wave propagation and the detector properties is of crucial importance. Inspired by the multi-stripe integration model in clinical X-ray computed tomography systems, in this work, we introduce the Multi-Curve-Integration-based acoustic inversion for cross-sectional Photoacoustic Tomography (MCI-PAT).

METHODS

We assumed that in cross-sectional PAT system, the three-dimensional (3-D) wave propagation problem could be reduced to a two-dimensional (2-D) problem in a limited, yet sufficient field of view. Under such condition, the MCI-PAT imaging model is generated by integrating several circular acoustic curves, the centers of which are points evenly distributed on the finite-size ultrasonic transducer surface. In this way, the spatial detector response is taken into account, while its computational burden does not largely increase because the integration process is performed only on a 2-D plane.

RESULTS

As proven by simulation, phantom and in vivo small animal experiments, the MCI-PAT method leads to promising improvement in PAT image quality. Comparing to traditional imaging models that considered only a single acoustic curve, our proposed method successfully improved the visibility of small structures and achieved evident enhancement on signal-to-noise ratio.

CONCLUSIONS

The performance of the MCI-PAT method demonstrates that for cross-sectional PAT, a 2-D simplification of the propagation of multiple photoacoustic waves is feasible. Due to its simplicity, our method can be used as an add-on to current system models considering only a single acoustic curve.

Multi-Operator Minimum Variance Adaptive Beamforming Algorithms Accelerated With GPU

The goal of this work is to design high-resolution, high-contrast and robust MV adaptive beamforming algorithms, which are also implemented in real-time frame rate. Multi-operator optimization is introduced into MV adaptive beamforming in this work to propose a multi-operator MV adaptive beamforming algorithmic optimization framework. Based on the proposed algorithmic optimization framework, the algorithm optimization can be either conducted by activating a single optimization operator, or conducted by activating multiple optimization operators. The multi-operator MV (MOMV) adaptive beamforming algorithms are then derived from this framework. Moreover, in order to promote the real-time imaging capability of MOMV beamforming, a GPU-based parallel acceleration framework is proposed along with the algorithmic optimization framework by exploring the image-level coarse-grained parallelization and pixel-level fine-grained parallelization. GPU computing resource allocation strategy and memory access strategy are both explored to better design the acceleration framework. Comprehensive quantitative simulation evaluations and qualitative in vivo experiments of imaging performance are studied, and the results demonstrate that the proposed MOMV adaptive beamforming algorithms significantly improve the imaging performance as compared with other MV beamforming algorithms, which have high resolution, high contrast, good robustness, and real-time imaging capability with thousands of acceleration speedup. Furthermore, the MOMV beamforming algorithm without eigen-decomposition and projection optimization operator achieves much higher beamforming frame rate with little downgrade of image quality as compared with the MOMV beamforming algorithm with all optimization operators.

Graphics processing unit accelerating compressed sensing photoacoustic computed tomography with total variation

Photoacoustic computed tomography with compressed sensing (CS-PACT) is a commonly used imaging strategy for sparse-sampling PACT. However, it is very time-consuming because of the iterative process involved in the image reconstruction. In this paper, we present a graphics processing unit (GPU)-based parallel computation framework for total-variation-based CS-PACT and adapted into a custom-made PACT system. Specifically, five compute-intensive operators are extracted from the iteration algorithm and are redesigned for parallel performance on a GPU. We achieved an image reconstruction speed 24-31 times faster than the CPU performance. We performed in vivo experiments on human hands to verify the feasibility of our developed method.

Evaluation Algorithm of Root Canal Shape Based on Steklov Spectrum Analysis

In recent years, we have seen more and more interest in the field of medical images and shape comparison motivated by the latest advances in microcomputed tomography (μCT) acquisition, modelling, and visualization technologies. Usually, biologists need to evaluate the effect of different root canal preparation systems. Current root canal preparation evaluation methods are based on the volume difference, area difference, and transportation of two root canals before and after treatment. The purpose of root canal preparation is to minimize the volume difference and ensure the complete removal of the smear layer. Previous methods can reflect some general geometric differences, but they are not enough to evaluate the quality of root canal shape. To solve this problem, we proposed a novel root canal evaluation method based on spectrum and eigenfunctions of Steklov operators, which can be served as a better alternative to current methods in root canal preparation evaluation. Firstly, the ideal root canal model was simulated according to the root canal model before and after preparation. Secondly, the Steklov spectrum of the two models was calculated. Thirdly, based on the spectrum and the histogram of the Gaussian curvature on the surface, the weight of each eigenvalue was computed. Therefore, the Steklov spectrum distance (SSD), which measures shape difference between the root canals, was defined. Finally, the calculation method that quantifies the root canal preparation effect of root canals was obtained. Through experiments, our method manifested high robustness and accuracy compared with existing state-of-the-art approaches. It also demonstrates the significance of our algorithm's advantages on a variety of challenging root canals through result comparison with counterpart methods.

GPU-accelerated Double-stage Delay-multiply-and-sum Algorithm for Fast Photoacoustic Tomography Using LED Excitation and Linear Arrays

Double-stage delay-multiply-and-sum (DS-DMAS) is an algorithm proposed for photoacoustic image reconstruction. The DS-DMAS algorithm offers a higher contrast than conventional delay-and-sum and delay-multiply and-sum but at the expense of higher computational complexity. Here, we utilized a compute unified device architecture (CUDA) graphics processing unit (GPU) parallel computation approach to address the high complexity of the DS-DMAS for photoacoustic image reconstruction generated from a commercial light-emitting diode (LED)-based photoacoustic scanner. In comparison with a single-threaded central processing unit (CPU), the GPU approach increased speeds by nearly 140-fold for 1024 × 1024 pixel image; there was no decrease in accuracy. The proposed implementation makes it possible to reconstruct photoacoustic images with frame rates of 250, 125, and 83.3 when the images are 64 × 64, 128 × 128, and 256 × 256, respectively. Thus, DS-DMAS can be efficiently used in clinical devices when coupled with CUDA GPU parallel computation.

Optimizing photoacoustic image reconstruction using cross-platform parallel computation

Three-dimensional (3D) image reconstruction involves the computations of an extensive amount of data that leads to tremendous processing time. Therefore, optimization is crucially needed to improve the performance and efficiency. With the widespread use of graphics processing units (GPU), parallel computing is transforming this arduous reconstruction process for numerous imaging modalities, and photoacoustic computed tomography (PACT) is not an exception. Existing works have investigated GPU-based optimization on photoacoustic microscopy (PAM) and PACT reconstruction using compute unified device architecture (CUDA) on either C++ or MATLAB only. However, our study is the first that uses cross-platform GPU computation. It maintains the simplicity of MATLAB, while improves the speed through CUDA/C++ − based MATLAB converted functions called MEXCUDA. Compared to a purely MATLAB with GPU approach, our cross-platform method improves the speed five times. Because MATLAB is widely used in PAM and PACT, this study will open up new avenues for photoacoustic image reconstruction and relevant real-time imaging applications.

A review of GPU-based medical image reconstruction

Tomographic image reconstruction is a computationally demanding task, even more so when advanced models are used to describe a more complete and accurate picture of the image formation process. Such advanced modeling and reconstruction algorithms can lead to better images, often with less dose, but at the price of long calculation times that are hardly compatible with clinical workflows. Fortunately, reconstruction tasks can often be executed advantageously on Graphics Processing Units (GPUs), which are exploited as massively parallel computational engines. This review paper focuses on recent developments made in GPU-based medical image reconstruction, from a CT, PET, SPECT, MRI and US perspective. Strategies and approaches to get the most out of GPUs in image reconstruction are presented as well as innovative applications arising from an increased computing capacity. The future of GPU-based image reconstruction is also envisioned, based on current trends in high-performance computing.

Efficient 3-D Model-Based Reconstruction Scheme for Arbitrary Optoacoustic Acquisition Geometries

Optimal optoacoustic tomographic sampling is often hindered by the frequency-dependent directivity of ultrasound sensors, which can only be accounted for with an accurate 3-D model. Herein, we introduce a 3-D model-based reconstruction method applicable to optoacoustic imaging systems employing detection elements with arbitrary size and shape. The computational complexity and memory requirements are mitigated by introducing an efficient graphic processing unit (GPU)-based implementation of the iterative inversion. On-the-fly calculation of the entries of the model-matrix via a small look-up table avoids otherwise unfeasible storage of matrices typically occupying more than 300GB of memory. Superior imaging performance of the suggested method with respect to standard optoacoustic image reconstruction methods is first validated quantitatively using tissue-mimicking phantoms. Significant improvements in the spatial resolution, contrast to noise ratio and overall 3-D image quality are also reported in real tissues by imaging the finger of a healthy volunteer with a hand-held volumetric optoacoustic imaging system.

A Multi-Grid Iterative Method for Photoacoustic Tomography

Inspired by the recent advances on minimizing nonsmooth or bound-constrained convex functions on models using varying degrees of fidelity, we propose a line search multi-grid (MG) method for full-wave iterative image reconstruction in photoacoustic tomography (PAT) in heterogeneous media. To compute the search direction at each iteration, we decide between the gradient at the target level, or alternatively an approximate error correction at a coarser level, relying on some predefined criteria. To incorporate absorption and dispersion, we derive the analytical adjoint directly from the first-order acoustic wave system. The effectiveness of the proposed method is tested on a total-variation penalized Iterative Shrinkage Thresholding algorithm (ISTA) and its accelerated variant (FISTA), which have been used in many studies of image reconstruction in PAT. The results show the great potential of the proposed method in improving speed of iterative image reconstruction.

Impact of nonstationary optical illumination on image reconstruction in optoacoustic tomography

Optoacoustic tomography (OAT), also known as photoacoustic tomography, is a rapidly emerging hybrid imaging technique that possesses great potential for a wide range of biomedical imaging applications. In OAT, a laser is employed to illuminate the tissue of interest and acoustic signals are produced via the photoacoustic effect. From these data, an estimate of the distribution of the absorbed optical energy density within the tissue is reconstructed, referred to as the object function. This quantity is defined, in part, by the distribution of light fluence within the tissue that is established by the laser source. When performing three-dimensional imaging of large objects, such as a female human breast, it can be difficult to achieve a relatively uniform coverage of light fluence within the volume of interest when the position of the laser source is fixed. To circumvent this, researchers have proposed illumination schemes in which the relative position of the laser source and ultrasound probe is fixed, and both are rotated together to acquire a tomographic dataset. A problem with this rotating-illumination scheme is that the tomographic data are inconsistent; namely, the acoustic data recorded at each tomographic view angle (i.e., probe position) are produced by a distinct object function. In this work, the impact of this data inconsistency on image reconstruction accuracy is investigated systematically. This is accomplished by use of computer-simulation studies and application of mathematical results from the theory of microlocal analysis. These studies specify the set of image discontinuities that can be stably reconstructed with a nonstationary optical illumination setup. The study also includes a comparison of the ability of iterative and analytic image reconstruction methods to mitigate artifacts attributable to the data inconsistency.

Regularized Iterative Weighted Filtered Back-Projection for Few-View Data Photoacoustic Imaging

Photoacoustic imaging is an emerging noninvasive imaging technique with great potential for a wide range of biomedical imaging applications. However, with few-view data the filtered back-projection method will create streak artifacts. In this study, the regularized iterative weighted filtered back-projection method was applied to our photoacoustic imaging of the optical absorption in phantom from few-view data. This method is based on iterative application of a nonexact 2DFBP. By adding a regularization operation in the iterative loop, the streak artifacts have been reduced to a great extent and the convergence properties of the iterative scheme have been improved. Results of numerical simulations demonstrated that the proposed method was superior to the iterative FBP method in terms of both accuracy and robustness to noise. The quantitative image evaluation studies have shown that the proposed method outperforms conventional iterative methods.

Real-Time Model-Based Inversion in Cross-Sectional Optoacoustic Tomography

Analytical (closed-form) inversion schemes have been the standard approach for image reconstruction in optoacoustic tomography due to their fast reconstruction abilities and low memory requirements. Yet, the need for quantitative imaging and artifact reduction has led to the development of more accurate inversion approaches, which rely on accurate forward modeling of the optoacoustic wave generation and propagation. In this way, multiple experimental factors can be incorporated, such as the exact detection geometry, spatio-temporal response of the transducers, and acoustic heterogeneities. The model-based inversion commonly results in very large sparse matrix formulations that require computationally extensive and memory demanding regularization schemes for image reconstruction, hindering their effective implementation in real-time imaging applications. Herein, we introduce a new discretization procedure for efficient model-based reconstructions in two-dimensional optoacoustic tomography that allows for parallel implementation on a graphics processing unit (GPU) with a relatively low numerical complexity. By on-the-fly calculation of the model matrix in each iteration of the inversion procedure, the new approach results in imaging frame rates exceeding 10 Hz, thus enabling real-time image rendering using the model-based approach.

Three-Dimensional Optoacoustic and Laser-Induced Ultrasound Tomography System for Preclinical Research in Mice: Design and Phantom Validation

In this work, we introduce a novel three-dimensional imaging system for in vivo high-resolution anatomical and functional whole-body visualization of small animal models developed for preclinical and other type of biomedical research. The system (LOUIS-3DM) combines a multiwavelength optoacoustic tomography (OAT) and laser-induced ultrasound tomography (LUT) to obtain coregistered maps of tissue optical absorption and speed of sound, displayed within the skin outline of the studied animal. The most promising applications of the LOUIS-3DM include 3D angiography, cancer research, and longitudinal studies of biological distributions of optoacoustic contrast agents.

A Constrained Variable Projection Reconstruction Method for Photoacoustic Computed Tomography Without Accurate Knowledge of Transducer Responses

Photoacoustic computed tomography (PACT) is an emerging computed imaging modality that exploits optical contrast and ultrasonic detection principles to form images of the absorbed optical energy density within tissue. When the imaging system employs conventional piezoelectric ultrasonic transducers, the ideal photoacoustic (PA) signals are degraded by the transducers' acousto-electric impulse responses (EIRs) during the measurement process. If unaccounted for, this can degrade the accuracy of the reconstructed image. In principle, the effect of the EIRs on the measured PA signals can be ameliorated via deconvolution; images can be reconstructed subsequently by application of a reconstruction method that assumes an idealized EIR. Alternatively, the effect of the EIR can be incorporated into an imaging model and implicitly compensated for during reconstruction. In either case, the efficacy of the correction can be limited by errors in the assumed EIRs. In this work, a joint optimization approach to PACT image reconstruction is proposed for mitigating errors in reconstructed images that are caused by use of an inaccurate EIR. The method exploits the bi-linear nature of the imaging model and seeks to refine the measured EIR during the process of reconstructing the sought-after absorbed optical energy density. Computer-simulation and experimental studies are conducted to investigate the numerical properties of the method and demonstrate its value for mitigating image distortions and enhancing the visibility of fine structures.

Expediting model-based optoacoustic reconstructions with tomographic symmetries

PURPOSE

Image quantification in optoacoustic tomography implies the use of accurate forward models of excitation, propagation, and detection of optoacoustic signals while inversions with high spatial resolution usually involve very large matrices, leading to unreasonably long computation times. The development of fast and memory efficient model-based approaches represents then an important challenge to advance on the quantitative and dynamic imaging capabilities of tomographic optoacoustic imaging.

METHODS

Herein, a method for simplification and acceleration of model-based inversions, relying on inherent symmetries present in common tomographic acquisition geometries, has been introduced. The method is showcased for the case of cylindrical symmetries by using polar image discretization of the time-domain optoacoustic forward model combined with efficient storage and inversion strategies.

RESULTS

The suggested methodology is shown to render fast and accurate model-based inversions in both numerical simulations and post mortem small animal experiments. In case of a full-view detection scheme, the memory requirements are reduced by one order of magnitude while high-resolution reconstructions are achieved at video rate.

CONCLUSIONS

By considering the rotational symmetry present in many tomographic optoacoustic imaging systems, the proposed methodology allows exploiting the advantages of model-based algorithms with feasible computational requirements and fast reconstruction times, so that its convenience and general applicability in optoacoustic imaging systems with tomographic symmetries is anticipated.

Discrete imaging models for three-dimensional optoacoustic tomography using radially symmetric expansion functions

Optoacoustic tomography (OAT), also known as photoacoustic tomography, is an emerging computed biomedical imaging modality that exploits optical contrast and ultrasonic detection principles. Iterative image reconstruction algorithms that are based on discrete imaging models are actively being developed for OAT due to their ability to improve image quality by incorporating accurate models of the imaging physics, instrument response, and measurement noise. In this work, we investigate the use of discrete imaging models based on Kaiser-Bessel window functions for iterative image reconstruction in OAT. A closed-form expression for the pressure produced by a Kaiser-Bessel function is calculated, which facilitates accurate computation of the system matrix. Computer-simulation and experimental studies are employed to demonstrate the potential advantages of Kaiser-Bessel function-based iterative image reconstruction in OAT.

Deblurring algorithms accounting for the finite detector size in photoacoustic tomography

Most reconstruction algorithms for photoacoustic tomography, like back projection or time reversal, work ideally for point-like detectors. For real detectors, which integrate the pressure over their finite size, images reconstructed by these algorithms show some blurring. Iterative reconstruction algorithms using an imaging matrix can take the finite size of real detectors directly into account, but the numerical effort is significantly higher compared to the use of direct algorithms. For spherical or cylindrical detection surfaces, the blurring caused by a finite detector size is proportional to the distance from the rotation center (spin blur) and is equal to the detector size at the detection surface. In this work, we apply deconvolution algorithms to reduce this type of blurring on simulated and on experimental data. Two particular deconvolution methods are compared, which both utilize the fact that a representation of the blurred image in polar coordinates decouples pixels at different radii from the rotation center. Experimental data have been obtained with a flat, rectangular piezoelectric detector measuring signals around a plastisol cylinder containing various small photoacoustic sources with variable distance from the center. Both simulated and experimental results demonstrate a nearly complete elimination of spin blur.

Investigation of the far-field approximation for modeling a transducer's spatial impulse response in photoacoustic computed tomography

When ultrasonic transducers with large detecting areas and/or compact measurement geometries are employed in photoacoustic computed tomography (PACT), the spatial resolution of reconstructed images can be significantly degraded. Our goal in this work is to clarify the domain of validity of an imaging model that mitigates such effects by use of a far-field approximation. Computer-simulation studies are described that demonstrate the far-field-based imaging model is highly accurate for a practical 3D PACT imaging geometry employed in an existing small animal imaging system. For use in special cases where the far-field approximation is violated, an extension of the far-field-based imaging model is proposed that divides the transducer face into a small number of rectangular patches that are each described accurately by use of the far-field approximation.

Model-based optoacoustic image reconstruction of large three-dimensional tomographic datasets acquired with an array of directional detectors

Image quality in 3-D optoacoustic (photoacoustic) tomography is greatly influenced by both the measurement system, in particular the number and spatial arrangement of ultrasound sensors, and the ability to account for the spatio-temporal response of the sensor element(s) in the reconstruction algorithm. Herein we present a reconstruction procedure based on the inversion of a time-domain forward model incorporating the spatial impulse response due to the shape of the transducer, which is subsequently applied in a tomographic system based on a translation-rotation scan of a linear detector array. The proposed method was also adapted to cope with the data-intensive requirements of high-resolution volumetric optoacoustic imaging. The processing of 2 · 10 (4) individual signals resulted in well-resolved images of both ~ 200 μm absorbers in phantoms and complex vascular structures in biological tissue. The results reported herein demonstrate that the introduced model-based methodology exhibits a better contrast and resolution than standard back-projection and model-based algorithms that assume point detectors. Moreover, the capability of handling large datasets anticipates that model-based methods incorporating the sensor properties can become standard practice in volumetric opto acoustic image formation.

Acoustic Inversion in Optoacoustic Tomography: A Review

Optoacoustic tomography enables volumetric imaging with optical contrast in biological tissue at depths beyond the optical mean free path by the use of optical excitation and acoustic detection. The hybrid nature of optoacoustic tomography gives rise to two distinct inverse problems: The optical inverse problem, related to the propagation of the excitation light in tissue, and the acoustic inverse problem, which deals with the propagation and detection of the generated acoustic waves. Since the two inverse problems have different physical underpinnings and are governed by different types of equations, they are often treated independently as unrelated problems. From an imaging standpoint, the acoustic inverse problem relates to forming an image from the measured acoustic data, whereas the optical inverse problem relates to quantifying the formed image. This review focuses on the acoustic aspects of optoacoustic tomography, specifically acoustic reconstruction algorithms and imaging-system practicalities. As these two aspects are intimately linked, and no silver bullet exists in the path towards high-performance imaging, we adopt a holistic approach in our review and discuss the many links between the two aspects. Four classes of reconstruction algorithms are reviewed: time-domain (so called back-projection) formulae, frequency-domain formulae, time-reversal algorithms, and model-based algorithms. These algorithms are discussed in the context of the various acoustic detectors and detection surfaces which are commonly used in experimental studies. We further discuss the effects of non-ideal imaging scenarios on the quality of reconstruction and review methods that can mitigate these effects. Namely, we consider the cases of finite detector aperture, limited-view tomography, spatial under-sampling of the acoustic signals, and acoustic heterogeneities and losses.

Volumetric real-time tracking of peripheral human vasculature with GPU-accelerated three-dimensional optoacoustic tomography

Optoacoustic tomography provides a unique possibility for ultra-high-speed 3-D imaging by acquiring complete volumetric datasets from interrogation of tissue by a single nanosecond-duration laser pulse. Yet, similarly to ultrasound, optoacoustics is a time-resolved imaging method, thus, fast 3-D imaging implies real-time acquisition and processing of high speed data from hundreds of detectors simultaneously, which presents significant technological challenges. Herein we present a highly efficient graphical processing unit (GPU) framework for real-time reconstruction and visualization of 3-D tomographic optoacoustic data. By utilizing a newly developed 3-D optoacoustic scanner, which simultaneously acquires signals with a handheld 256-element spherical ultrasonic array system, we further demonstrate tracking of deep tissue human vasculature rendered at a rate of 10 volumetric frames per second. The flexibility provided by the handheld hardware design, combined with the real-time operation, makes the developed platform highly usable for both clinical imaging practice and small animal research applications.

Optoacoustic imaging and tomography: reconstruction approaches and outstanding challenges in image performance and quantification

This paper comprehensively reviews the emerging topic of optoacoustic imaging from the image reconstruction and quantification perspective. Optoacoustic imaging combines highly attractive features, including rich contrast and high versatility in sensing diverse biological targets, excellent spatial resolution not compromised by light scattering, and relatively low cost of implementation. Yet, living objects present a complex target for optoacoustic imaging due to the presence of a highly heterogeneous tissue background in the form of strong spatial variations of scattering and absorption. Extracting quantified information on the actual distribution of tissue chromophores and other biomarkers constitutes therefore a challenging problem. Image quantification is further compromised by some frequently-used approximated inversion formulae. In this review, the currently available optoacoustic image reconstruction and quantification approaches are assessed, including back-projection and model-based inversion algorithms, sparse signal representation, wavelet-based approaches, methods for reduction of acoustic artifacts as well as multi-spectral methods for visualization of tissue bio-markers. Applicability of the different methodologies is further analyzed in the context of real-life performance in small animal and clinical in-vivo imaging scenarios.