Source-linked AI summary
Accelerated High-Resolution Photoacoustic Tomography via Compressed Sensing
Simon Arridge, Paul Beard, Marta Betcke, Ben Cox, Nam Huynh, Felix Lucka, Olumide Ogunlade, Edward Zhang
TL;DR
High-resolution 3D PAT requires slow, densely sampled acquisition, limiting simultaneous spatial and temporal resolution. The paper combines spatially sub-sampled sensing with sparsity-constrained variational reconstruction and reports good image quality from highly sub-sampled data, while identifying model accuracy and target complexity as important boundaries.
Problem
3D PAT systems generally provide either high image quality or high frame rates, limiting simultaneous high spatial and temporal resolution for dynamic imaging.
Method
The paper combines random or patterned spatial sub-sampling with variational reconstruction using spatial sparsity constraints, including total variation with Bregman iterations.
Results
Variational sparsity-constrained methods produced good-quality PAT images from sub-sampled data, whereas linear and L2+ methods did not produce acceptable images in any setting.
Takeaways & Limitations
The approach offers opportunities to accelerate sequential PAT scanners and reduce channel counts in parallelized detector-array schemes.
Takeaways & Limitations
Performance depends on accurate forward-model alignment and target difficulty; the easy Tumor1 simulation reached Msub = 128, whereas Tumor2 degraded beyond Msub = 16–32.
Abstract
from arXiv · showhide
Current 3D photoacoustic tomography (PAT) systems offer either high image quality or high frame rates but are not able to deliver high spatial and temporal resolution simultaneously, which limits their ability to image dynamic processes in living tissue. A particular example is the planar Fabry-Perot (FP) scanner, which yields high-resolution images but takes several minutes to sequentially map the photoacoustic field on the sensor plane, point-by-point. However, as the spatio-temporal complexity of many absorbing tissue structures is rather low, the data recorded in such a conventional, regularly sampled fashion is often highly redundant. We demonstrate that combining variational image reconstruction methods using spatial sparsity constraints with the development of novel PAT acquisition systems capable of sub-sampling the acoustic wave field can dramatically increase the acquisition speed while maintaining a good spatial resolution: First, we describe and model two general spatial sub-sampling schemes. Then, we discuss how to implement them using the FP scanner and demonstrate the potential of these novel compressed sensing PAT devices through simulated data from a realistic numerical phantom and through measured data from a dynamic experimental phantom as well as from in-vivo experiments. Our results show that images with good spatial resolution and contrast can be obtained from highly sub-sampled PAT data if variational image reconstruction methods that describe the tissues structures with suitable sparsity-constraints are used. In particular, we examine the use of total variation regularization enhanced by Bregman iterations. These novel reconstruction strategies offer new opportunities to dramatically increase the acquisition speed of PAT scanners that employ point-by-point sequential scanning as well as reducing the channel count of parallelized schemes that use detector arrays.
1. Introduction
PAT combines optical absorption contrast with ultrasound spatial resolution, but high-resolution 3D imaging requires dense sampling over large apertures. Sequential scanning therefore requires thousands of detection points, motivating compressed acquisition strategies.
- PAT images optical absorption with ultrasound resolution and can potentially provide spectroscopic information about absorbing molecules.
- High-resolution 3D PAT requires tens-of-micrometres sampling intervals across centimetre-scale apertures, creating thousands of detection points.
- The paper combines spatially sparse variational reconstruction with sub-sampled acquisition systems to accelerate PAT while preserving image quality.
2. Background
PAT generates ultrasound from laser-induced heating and reconstructs the initial pressure distribution from boundary measurements. Complete sampling must satisfy linked temporal and spatial Nyquist criteria.
- Laser illumination produces localized pressure increases through rapid optical absorption and thermalization in tissue.
- The acoustic wave equation propagates the initial pressure distribution p0, which is reconstructed from pressure measurements on the domain boundary.
- Temporal sampling requires δt < 1/(2ω*_t), while spatial sampling must resolve wave vectors up to ∥k*∥2 = ω*_t/c.
- Figure 1 contrasts localized single-point interrogation with distributed binary spatial windowing of the Fabry-Pérot sensor.
3. Compressed Photoacoustic Sensing
The paper accelerates sequential PAT by replacing complete point-by-point spatial sampling with either random single-point measurements or distributed orthogonal patterns. Both schemes are modeled as windowed integration of the acoustic field and implemented with a Fabry-Pérot scanner.
- Spatial sub-sampling accelerates sequential PAT because the measured field is acquired with fewer locations or patterns than conventional scanning.
- Each measurement multiplies the detection-plane field by a spatial window and integrates it over the detection surface before temporal sampling.
- Random single-point sampling selects Mc locations from a conventional grid, yielding an acceleration factor Msub = M/Mc.
- Distributed patterned sampling uses fewer orthogonal binary windows supported across the full interrogation area, also achieving Msub = M/Mc.
- The Fabry-Pérot scanner implements these strategies by interrogating its planar sensor with focused or patterned laser beams.
4. Image Reconstruction from Sub-Sampled Data
The paper reconstructs images directly from sub-sampled PAT measurements using forward-model-based variational methods, with total variation and Bregman iterations addressing sparsity and contrast bias.
- Forward and sensing models: The forward model maps initial pressure p0 to detection-plane pressure, while sensing operators generate noisy sub-sampled measurements from that field.The sensing operator represents conventional scanning, random single-point sub-sampling, or patterned interrogation.
- Reconstruction strategies: Reconstructions can either first recover complete detection-plane data and then apply standard PAT inversion, or directly estimate p0 through a model-based one-step procedure.The study emphasizes differences between linear reconstruction and variational model-based reconstruction rather than a detailed comparison of one-step and two-step procedures.
- Variational reconstruction: Variational reconstruction minimizes data misfit plus a regularization functional, whose parameter balances fidelity to measured data against prior constraints.Under additive i.i.d. Gaussian noise, the data-fidelity term is combined with a regularizer to obtain a stable solution.
- Total variation regularization: Total variation penalizes the ℓ1 norm of gradient amplitudes, providing an edge-preserving spatial sparsity constraint suited to recovering feature edges from highly sub-sampled data.The examined convex optimization problem is solved with an accelerated proximal-gradient-type method.
- Bregman iterations: TV regularization introduces systematic contrast loss, while Bregman iterations compensate for this bias through a sequence of regularized problems with decreasing residuals.This approach targets the contrast distortion that limits TV-based images for quantitative PAT analysis.
5. Simulation Studies
Simulation studies use realistic mouse-brain-derived phantoms to compare reconstruction methods under spatial sub-sampling. At a factor of 128, TV-based methods preserve main phantom structures, and Bregman iterations improve small-vessel contrast, though the easiest simulation requires caution.
- Realistic numerical phantom: The realistic numerical phantom is derived from a segmented mouse-brain micro-CT scan containing gray matter, vasculature, and dura mater.Tumor1 and Tumor2 assign different initial-pressure values to tumor tissue before down-sampling to the target resolution.
- Simulation setup: Tumor1 simulations use a 128^3 reconstruction grid, exact forward-model knowledge, and spatial sampling at the computational-grid spacing.The setup is explicitly described as inverse-crime data, with conventional sampling fine enough to satisfy the Nyquist criterion.
- Reconstruction results: Despite 128-fold sub-sampling, TV+ and TV+Br reconstruct the phantom’s main structures without excessive image noise, unlike the unsatisfactory BP and TR reconstructions.The comparison uses rSP-128 and sHd-128 data against conventional scanning; sub-sampling artifacts differ between the two acquisition strategies.
- Bregman enhancement: Bregman iterations especially improve the contrast of small-scale vessel structures, with a stronger benefit for sub-sampled than conventional data.The comparison uses identical color scaling and difference projections between TV+Br and TV+ solutions.
5.3. Simulation Studies with Tumor2
Tumor2 simulations model sensor and noise uncertainties to approximate experimental conditions, showing that TV+Br retains acceptable reconstructions at substantially higher sub-sampling factors than random single-point sampling.
- Simulation setup: The very high-quality Tumor1 results at sub-sampling factor 128 require caution because the phantom was close to the sensors, high-contrast, and generated with the reconstruction model.These conditions constitute an inverse crime and limit extrapolation to experimental data.
- Simulation setup: Tumor2 simulations incorporate spatially varying sensor sensitivity, uncertain noise variance, and a coarser conventional sampling grid to better reproduce experimental data.The model uses random multiplicative factors for sensor sensitivity and noise variance, with σs = 0.2 and σn = 0.1.
- Reconstruction Results: At Msub = 32, sHd-based TV+Br reconstructions remain acceptable, whereas rSP-based reconstructions clearly degrade beyond Msub = 16.For rSP, image quality only slightly deteriorates up to Msub = 16 before clear degradation from 16 to 32.
- Reconstruction Results: TV+Br again gives the best visual and PSNR results among the variational methods for both rSP-16 and sHd-16 data.The regularization parameter was selected using the discrepancy principle with κ = 1.25.
- Sampling patterns: The random single-point pattern has only a minor influence compared with the inverse method, while optimal dynamic sampling-pattern design remains future work.The comparison was made against a regular coarse-grid pattern.
6. Experimental Data
Experimental evaluation uses dynamic and static phantoms plus in-vivo vasculature, with preprocessing and acoustic-model choices tailored to each dataset. These experiments test whether sub-sampling performance transfers beyond idealized simulations.
- Experimental datasets: Experimental evaluation covers a pseudo-dynamic knot phantom, a patterned-interrogation hair phantom, and static in-vivo mouse vessels.The datasets span conventional single-point scanning, patterned interrogation, and real biological anatomy.
- Knot phantom: The knot experiment acquired a conventional scan in approximately 15 minutes before rotating the motor shaft to create motion between repeated scans.Two ink-filled polythene tubes formed the moving knot, with one end attached to a motor shaft.
- Preprocessing: Preprocessing included baseline subtraction, zero-phase band-pass filtering, variance-based removal of 1% of locations, and temporal clipping.For the knot data, filtering used approximately 0.5–20 MHz and the retained time points were 10–400.
- Acoustic model and sampling: The acoustic reconstruction grid was twice finer than the scanning-location spacing because the conventional data were temporally oversampled but spatially undersampled.For 20 MHz signals and sound speed 1540 m s−1, the Nyquist criterion would require δy/z = 38.5 µm and δt = 25 ns.
- Hair phantom: The hair phantom used a 150 µm artificial hair target approximately 2 mm above the detection plane and 3 mm below the Intralipid surface, interrogated with scrambled Hadamard patterns.A 640 × 640 micromirror area was grouped into 128 × 128 pixels, and rows of a 16384 × 16384 scrambled Hadamard matrix generated the binary patterns.
- In-vivo vessels: For in-vivo vessels, measurements covered a 14 mm × 14 mm region with 142 × 141 locations and 630 temporal samples at 10 ns resolution.The excitation wavelength was 590 nm; the processed data were clipped to 141 × 141 locations and time points 10–630.
7. Discussion, Outlook and Conclusion
The study shows that suitable sparsity-constrained variational reconstruction can recover high-quality PAT images from highly sub-sampled data, but achievable rates depend strongly on target complexity and model fit. Experimental results are promising, while scanner limitations and forward-model mismatch remain important barriers to realizing the full acceleration potential.
- Variational reconstructions with spatial sparsity constraints were essential for obtaining acceptable-quality images from sub-sampled PAT data, unlike linear and L2+ methods.The finding was reported across the simulation studies.
- Tumor1 tolerated sub-sampling up to Msub = 128 without significant image-quality loss, whereas the mismatched Tumor2 case degraded beyond Msub = 16−32.The achievable rate therefore varied strongly with target difficulty and model agreement.
- Bregman iterations mitigated systematic contrast loss in small structures such as blood vessels, supporting their use in quantitative PAT studies.
- Patterned interrogation sHd was slightly more efficient than single-point sub-sampling rSP.
- Experimental results showed promising sub-sampling rates but weaker method differences than simulations, consistent with unmodeled sensor sensitivity and noise mismatches.The experimental data included known preprocessing, but not all model mismatches were accounted for.
- Future improvements require more accurate forward models and spatio-temporal variational methods that exploit temporal redundancy for 4D PAT.Suggested measures include calibration, uncertainty modeling, and incorporating acoustic absorption models.
Appendix A. Discrete Total Variation Energy
The appendix defines the discrete 3D total variation energy using finite forward differences and specifies how boundary conditions modify the discretization.
- The 3D pressure field is discretized over indexed voxels, and total variation is formed using finite forward differences with Neumann boundary conditions.
- Dirichlet boundary conditions require additional jump terms at domain boundaries in the simulations.
- Experimental reconstructions impose Dirichlet conditions on detection-plane voxels and Neumann conditions on the remaining image-cube faces.
Appendix B. Optimization
The optimization appendix formulates reconstruction as a smooth data term plus a convex regularizer and uses proximal first-order methods tailored to expensive forward and adjoint operators.
- The reconstruction minimizes E(p) = D(p) + λJ(p), combining a smooth strictly convex data term with a convex, potentially nondifferentiable regularizer.
- The gradient of D(p) is computed explicitly, while J(p) is handled through its proximal operator.
- For positivity-constrained total variation, the proximal step solves a positivity-constrained TV denoising problem.
- The implementation uses FISTA with step size η = 1.8/L, where L approximates the Lipschitz constant of A^T C^T C A, and restarts acceleration when the objective increases.
- The positivity-constrained TV denoising proximal operator is implemented with a primal-dual hybrid gradient algorithm.
Appendix C. Implementation
The implementation uses a PAT reconstruction toolbox built around k-Wave, with optimized C++ and CUDA routines for 3D wave propagation on CPU and GPU architectures.
- The reconstruction routines rely on k-Wave to compute the forward and adjoint operators using optimized C++ and CUDA code.These computations can run on parallel CPU or GPU architectures.
Tables and table captions
The paper includes tables identifying common abbreviations and listing parameters for the inversion models used.
- Table 1 lists commonly occurring abbreviations used in the paper.
- Table 2 lists parameters for the different inversion models used.