Source-linked AI summary
Full-Wave Iterative Image Reconstruction in Photoacoustic Tomography with Acoustically Inhomogeneous Media
Chao Huang, Kun Wang, Liming Nie, Lihong V. Wang, Mark A. Anastasio
TL;DR
PACT reconstruction in acoustically heterogeneous media is limited by model assumptions, computational burden, and sensitivity to incomplete or noisy data. The paper introduces an exact-wave-equation discrete model with matched operators and evaluates a k-space-based iterative methodology on simulated and experimental data, where the iterative results generally show lower RMSE and fewer artifacts than time reversal.
Problem
Existing heterogeneous-media PACT reconstruction methods are limited by weak-variation assumptions, computational burden, or inadequate mitigation of incomplete and noisy measurements.
Method
The paper develops a discrete imaging model based on the exact PA wave equation, with matched forward and backprojection operators implemented using the k-space pseudospectral method.
Results
The methodology mitigated artifacts from data incompleteness, noise, finite sampling, and modeling errors in computer-simulated and experimental PACT data, with iterative RMSE lower than time reversal in reported cases.
Takeaways & Limitations
The reconstruction method has potential for preclinical and clinical PACT applications.
Abstract
from arXiv · showhide
Existing approaches to image reconstruction in photoacoustic computed tomography (PACT) with acoustically heterogeneous media are limited to weakly varying media, are computationally burdensome, and/or cannot effectively mitigate the effects of measurement data incompleteness and noise. In this work, we develop and investigate a discrete imaging model for PACT that is based on the exact photoacoustic (PA) wave equation and facilitates the circumvention of these limitations. A key contribution of the work is the establishment of a procedure to implement a matched forward and backprojection operator pair associated with the discrete imaging model, which permits application of a wide-range of modern image reconstruction algorithms that can mitigate the effects of data incompleteness and noise. The forward and backprojection operators are based on the k-space pseudospectral method for computing numerical solutions to the PA wave equation in the time domain. The developed reconstruction methodology is investigated by use of both computer-simulated and experimental PACT measurement data.
I. INTRODUCTION
PACT combines optical contrast with ultrasound detection, but acoustic heterogeneity, attenuation, and incomplete measurements challenge existing reconstruction methods. The paper develops a full-wave iterative approach based on an exact discrete PA wave equation and matched operators.
- PACT background: PACT generates acoustic wavefields from pulsed optical illumination and reconstructs absorbed optical energy density from ultrasonic measurements.Optical absorption is related to tissue hemoglobin concentration and molecular constitution, supporting anatomical, functional, and molecular imaging applications.
- Limitations of existing methods: Most existing PACT reconstruction algorithms assume lossless, acoustically homogeneous media, although bone, gas pockets, and skull can strongly perturb PA wavefields.These violations can produce significant image distortions and artifacts.
- Limitations of existing methods: Ray-based methods compensate for weak speed-of-sound variations but require those variations to occur on length scales large relative to the effective acoustic wavelength.Those assumptions can be violated in preclinical and clinical PACT applications.
- Limitations of existing methods: Full-wave alternatives model the exact PA equation but can impose substantial computational burdens or practical sampling requirements.The cited limitations include FEM computational cost, time-reversal dependence on dense enclosing measurements, and unestablished robustness of Neumann-series methods to sparse sampling.
- Proposed approach: The proposed approach uses an exact-equation discrete imaging model with matched forward and backprojection operators implemented through the k-space pseudospectral method.The operators support iterative algorithms intended to mitigate artifacts from incomplete data, noise, finite sampling, and modeling errors.
- Proposed approach: The model describes wave propagation in media with inhomogeneous speed of sound and density plus frequency-power-law acoustic attenuation.The continuous formulation represents absorption and dispersion with lossy derivative operators based on the fractional Laplacian.
B. Photoacoustic wavefield propagation: Discrete formulation
The discrete propagation formulation samples pressure and particle velocity on a 3D Cartesian grid and advances them using Fourier-domain spatial derivatives. This k-space pseudospectral implementation reduces computational and sampling demands relative to real-space methods.
- Propagation method: The k-space pseudospectral method computes numerical PA wavefield solutions to the coupled propagation equations in space and time.It uses FFT-based spatial derivatives and has less restrictive spatial and temporal sampling requirements than real-space finite-element and finite-difference methods.
- Grid and field representation: The discrete formulation represents sampled pressure and particle-velocity fields as time-indexed 3D matrices on a Cartesian grid.The grid contains N=N1N2N3 vertices, and lexicographically ordered vector representations are used for matrix quantities.
- Medium representation: Ambient density and squared speed-of-sound distributions are represented discretely with diagonal matrix forms.The formulation defines diagonal matrices whose entries correspond to sampled medium properties.
- Spectral differentiation: One-dimensional spatial derivatives are computed in the Fourier domain using forward and inverse 3D discrete Fourier transforms and Hadamard products.The derivative construction uses wave-number matrices Ki defined over the spatial grid dimensions.
- Discrete operators: The discrete coupled equations incorporate absorption and dispersion through vectorized terms associated with the sampled pressure and particle-velocity fields.The relevant powers and matrix operators are applied element-wise or through lexicographically ordered vector representations.
C. The image reconstruction problem
The reconstruction problem estimates the initial pressure distribution from time-resolved measurements and known acoustic properties. Its discrete system matrix enables regularized iterative inversion using forward and adjoint operators.
- Data and assumptions: Measured pressure data are collected at L ultrasonic transducer positions outside the object over M temporal samples.The model treats transducer responses as negligible and transducers as point-like for simplicity.
- Reconstruction target: The reconstruction seeks the initial pressure distribution, or equivalently absorbed optical energy density, from measured data and known medium parameters.The known parameters include speed of sound, ambient density, attenuation coefficient, and frequency-power-law exponent.
- Discrete imaging model: The discrete imaging model maps the N×1 initial-pressure vector to measured pressure data through an LM×N system matrix H.The measured-data vector contains all transducer locations and temporal samples.
- Iterative reconstruction: Regularized inversion can be performed with iterative methods that minimize a penalized least-squares cost function using H and its adjoint H†.The action of both operators must be computed to apply these iterative methods.
III. EXPLICIT FORMULATION OF DISCRETE IMAGING
The paper formulates a discrete PACT imaging model by propagating wavefield variables with a matrix operator and relating initial pressure to measured transducer data. It then identifies the system matrix and its transpose for matched forward and backprojection.
- System matrix: The explicit system matrix H is obtained by combining the initial-state mapping, time propagation, and measurement interpolation operations.H maps the N × 1 discrete initial pressure vector to the measured pressure samples.
- Discrete wavefield propagation: The k-space formulation represents all wavefield variables at each time step with a 7N × 1 state vector and propagates them using a 7N × 7N matrix W.The state is advanced from t = 0 through the final discretized time using block propagation matrices defined from W.
- Discrete wavefield propagation: The propagator matrix W encodes the coupled updates of the discretized wavefield variables, including spatial-derivative and material-dependent blocks.Its block structure contains identity, derivative, and material-related matrices such as Di and E.
- Initial conditions: The initial pressure distribution p0 determines the initial wavefield state through the equation of state and the prescribed initial conditions.The resulting vector places the initial pressure information into the wavefield representation before time propagation.
- Measurement model: Computed wavefields are converted to measured pressure data through interpolation because transducer locations generally do not coincide with Cartesian-grid vertices.The interpolation uses the four nearest grid nodes and barycentric coordinates obtained from Delaunay triangulation.
- Matched backprojection: Because H has real-valued elements, the matched backprojection operator H† is its transpose HT, with a corresponding explicit matrix form.The two-dimensional formulation is included as a special case of the three-dimensional construction.
IV. DESCRIPTIONS OF NUMERICAL AND EXPERIMENTAL
The numerical studies evaluate the proposed discrete imaging model and its adjoint in iterative reconstruction from incomplete 2D and 3D PACT data. Reconstructions use a TV-regularized PLS objective and are compared with a TR-based method.
- Numerical studies: Numerical studies assessed the effectiveness and robustness of the proposed discrete imaging model for iterative reconstruction from incomplete 2D and 3D PACT data.
- Iterative reconstruction: The system matrix and its adjoint were incorporated into an iterative algorithm minimizing a PLS cost function with a total variation penalty.
- Comparison method: Performance was compared against an existing time-reversal reconstruction algorithm.
A. Implementation of the forward and backprojection operators
The operators are implemented with k-Wave-based k-space propagation and adjoint backpropagation, then used in constrained TV-regularized reconstruction. The implementation includes absorbing boundaries and GPU-accelerated numerical studies.
- Forward operator: The forward operator H is implemented with the k-space pseudospectral method in the MATLAB k-Wave toolbox.A perfectly matched layer surrounds the computational grid to prevent waves from leaving one side and re-entering on the opposite side.
- Backprojection operator: The backprojection HT is computed by applying transposed measurement interpolation and reverse-time transposed propagation operations.Sparse interpolation matrices and the nontrivial blocks of W are transposed directly; the PML modification is applied to both W and WT.
- Reconstruction algorithm: The iterative reconstruction solves a nonnegative optimization problem combining data fidelity with a TV regularization term.The regularization parameter λ is empirically selected and fixed across the studies.
- Comparison algorithm: The TR algorithm reconstructs images by solving the discretized acoustic equations backward in time with the same PML parameters used for system-matrix construction.
- Computational setup: The studies use 512 × 512 grids for 2D simulations and 256 × 256 × 128 grids for 3D simulations, with computations performed on a CPU-GPU workstation.
C. Computer-simulation studies of 2D PACT
The studies evaluate reconstruction across full-view, few-view, and limited-view geometries, using numerical and experimental PACT data with heterogeneous media and controlled map errors. Iterative reconstruction generally produced fewer artifacts and less noise than time reversal when measurements were incomplete or noisy, and achieved lower RMSE in the reported disc-phantom studies.
- Scanning geometries: Three scanning geometries tested robustness to different types and degrees of data incompleteness: full-view, few-view, and limited-view.Full-view used 180 transducers on a circle; few-view used 60, while limited-view used 90 on a semicircle.
- Numerical phantoms and data: The simulations used blood-vessel and low-contrast-disc numerical phantoms to investigate reconstruction robustness under incomplete data and noise.Pressure data were computed with the k-space pseudospectral method at the transducer locations.
- Errors in SOS and density maps: The simulations modeled SOS and density maps representative of a monkey skull and introduced 1.3% Gaussian noise plus a 7-pixel, 1.4 mm registration shift.These perturbations were used to investigate robustness to inevitable map-estimation errors.
- Different scanning geometries: With few-view and limited-view scanning, iterative reconstructions contained fewer artifacts and less noise than time-reversal results, with image values closer to the phantom.Under full-view scanning, both methods produced accurate reconstructed images.
- Noisy data: RMSE was 0.026 for time reversal and 0.007 for iterative reconstruction for the noisy low-contrast-disc phantom with error-free SOS and density maps.The reported results suggest stronger noise reduction by the iterative algorithm in this setting.
- Errors in SOS and density maps: RMSE was 0.086 for time reversal and 0.034 for iterative reconstruction when SOS and density maps contained errors.The iterative result also had a cleaner background and sharper disc boundaries, attributed to TV regularization.
C. 3D simulation results
The iterative method reconstructed 3D PACT images more accurately and with greater robustness to noise, incomplete measurements, and modeling errors than the TR method, but required substantially more computation.
- 3D simulation: The iterative method achieved an RMSE of 0.003 versus 0.018 for TR in the 3D blood-vessel phantom.Both reconstructions were visualized using maximum intensity projection in the same grayscale window.
- Computational cost: TR required approximately 6 minutes, whereas the iterative method with 10 iterations required 110 minutes.The reported computational comparison highlights a substantial runtime cost for iterative reconstruction.
- Experimental results: The iterative method produced higher-resolution images than TR in the few-view experiment and mitigated artifacts and distortions from incomplete measurement data.The paper attributes the higher resolution to TV regularization, which mitigates model errors such as neglected shear-wave and finite-aperture effects.
- Experimental results: In experimental few-view and limited-view scans, the iterative method achieved RMSE values of 0.002 and 0.003 versus 0.005 and 0.007 for TR.The experiments used an acrylic shell to introduce acoustic heterogeneity and incomplete measurement geometries.
- Overall assessment: The proposed methodology was investigated using both simulated and experimental data to mitigate artifacts from data incompleteness, noise, finite sampling, and modeling errors.The approach uses a full-wave discrete imaging model and matched forward and backprojection operators.
APPENDIX-A: MODELING TRANSDUCER IMPULSE
The appendix incorporates finite transducer impulse responses into the discrete PACT imaging model by modeling spatial averaging and electrical filtering within the system matrix.
- System-matrix formulation: The transducer impulse responses, including SIR and EIR, can be incorporated into the system matrix.This extends the discrete imaging operator beyond point-like transducer measurements.
- Spatial impulse response: The SIR models spatial averaging over each finite transducer surface by dividing it into small patches and approximating the surface integral by a sum.The patch area is chosen much smaller than the acoustic wavelength.
- Spatial impulse response: The averaged pressure at each transducer is formed from pressure samples at its surface patches using a matrix representation.The formulation assumes equal-area patches for simplicity but can be extended to transducer-specific patch counts and areas.
- Measured-data representation: The complete measured-data vector is assembled across all transducers and temporal samples before applying the system-matrix formulation.The construction relates point-like pressure data, spatially averaged pressure, and electrically filtered measurements.
- Electrical impulse response: The EIR models the piezoelectric transducer’s electrical response as convolution of the averaged pressure sequence with the EIR.The transducer is treated as a linear shift-invariant system, with identical EIRs assumed for simplicity.
ALGORITHM FOR PACT
The algorithm solves the TV-regularized optimization problem iteratively using a three-dimensional extension of FISTA and a separate denoising step. Its discrete operators use finite differences with specified zero boundary conditions.
- Optimization solver: The optimization problem in Eqn. (43) is solved iteratively using Algorithm 1.The solver uses a zero initial guess in all studies.
- Optimization solver: Algorithm 1 extends FISTA to three dimensions for the reconstruction optimization problem.The extension is attributed to the FISTA algorithm described in Ref..
- Denoising step: A denoising subproblem defined by Eqn. (60) is solved within the iterative procedure using Algorithm 2.The paper states that this subproblem can be solved efficiently.
- Boundary conditions: The discrete operators impose zero boundary values on the outer faces of the computational grid.The assumed boundary entries for a, b, and c are set to zero at the relevant index boundaries.
- Discrete operators: The discrete operator l maps a matrix y to three directional finite-difference matrices a, b, and c.The differences are defined along the first, second, and third dimensions over their respective interior index ranges.
- Discrete operators: The operator Pp maps the directional matrices (a, b, c) to corresponding output matrices (d, e, f), with component expressions specified for the discrete grid.The operator definitions use grid indices n1, n2, and n3 and complementary boundary assumptions.