PaperlayerView exact revision on arXiv

Source-linked AI summary

On the Adjoint Operator in Photoacoustic Tomography

Simon R. Arridge, Marta M. Betcke, Ben T. Cox, Felix Lucka, Brad E. Treeby

arXiv:1602.02027v2math.NAcs.NA

TL;DR

Incomplete or sub-sampled PAT data motivate variational reconstruction with prior knowledge, but that approach needs an efficient adjoint of the forward operator. The paper derives a continuous analytical adjoint and implements it with a k-space wave-propagation model; it reports sufficient precision against the discretize-then-adjoint approach and illustrates variational reconstruction.

  • Problem

    For incomplete or sub-sampled PAT data, variational reconstruction can incorporate prior knowledge, but first-order optimization requires an efficient numerical adjoint of the forward operator.

    Source: §1.2 · p. 3

  • Method

    The paper derives the continuous analytical PAT adjoint and implements it with a k-space wave-propagation model.

    Source: §5 · pp. 13–14

  • Results

    The adjoint-then-discretize approach was sufficiently precise compared with the discretize-then-adjoint approach, and the adjoint was illustrated in variational reconstruction with a-priori constraints.

    Source: §5 · p. 15

  • Takeaways & Limitations

    The explicit continuous adjoint can be implemented with existing wave-propagation codes that support a time-varying mass source.

    Source: §5 · p. 14

    The optimization schemes used were computationally inefficient and not tuned for optimal performance, so the example iteration counts do not indicate expected computational cost.

    Source: §5 · p. 15

Abstract

from arXiv · show

Photoacoustic Tomography (PAT) is an emerging biomedical "imaging from coupled physics" technique, in which the image contrast is due to optical absorption, but the information is carried to the surface of the tissue as ultrasound pulses. Many algorithms and formulae for PAT image reconstruction have been proposed for the case when a complete data set is available. In many practical imaging scenarios, however, it is not possible to obtain the full data, or the data may be sub-sampled for faster data acquisition. In such cases, image reconstruction algorithms that can incorporate prior knowledge to ameliorate the loss of data are required. Hence, recently there has been an increased interest in using variational image reconstruction. A crucial ingredient for the application of these techniques is the adjoint of the PAT forward operator, which is described in this article from physical, theoretical and numerical perspectives. First, a simple mathematical derivation of the adjoint of the PAT forward operator in the continuous framework is presented. Then, an efficient numerical implementation of the adjoint using a k-space time domain wave propagation model is described and illustrated in the context of variational PAT image reconstruction, on both 2D and 3D examples including inhomogeneous sound speed. The principal advantage of this analytical adjoint over an algebraic adjoint (obtained by taking the direct adjoint of the particular numerical forward scheme used) is that it can be implemented using currently available fast wave propagation solvers.

1. Introduction

PAT uses optical absorption to generate ultrasound pulses that carry tissue-contrast information to surface sensors. For incomplete or sub-sampled data, variational reconstruction can use prior knowledge, motivating an efficient, discretization-independent adjoint of the forward operator.

1. Introduction

  • PAT converts optical absorption into broadband ultrasound pulses that carry contrast information to the tissue surface, supporting high-resolution soft-tissue imaging.The measured acoustic data can be used to infer the initial pressure distribution; recovering optical properties requires a separate light-transport inversion not considered here.Source: §1 · p. 2, §1.1 · p. 3

1.1. The Physics of Photoacoustic Tomography

  • The initial pressure is inferred from measured acoustic time series, while recovery of optical properties requires a subsequent light-transport inversion outside this paper’s scope.The acoustic inversion supplies data for that later optical inversion.Source: §1.1 · p. 3

1.2. Challenges of Photoacoustic Tomography

  • Variational reconstruction can incorporate prior knowledge when data are incomplete or sub-sampled, but first-order optimization requires an efficient numerical adjoint of the PAT forward operator.The regularization functional encodes prior knowledge about the initial pressure.Source: §1.2 · p. 3
  • The paper derives a continuous analytical adjoint independent of discretization, clarifies its relation to time reversal, and implements it with existing wave-propagation codes.The implementation uses a k-space pseudospectral wave-propagation model.Source: §1.2 · p. 4

2. Theory

The theory models PAT as acoustic wave propagation from an initial pressure distribution followed by a limited sensor measurement. It derives the continuous adjoint as a wave-equation solution evaluated at the recording time and distinguishes its source-term formulation from time reversal’s boundary constraint.

2.1. Mathematical Modeling

  • The forward model starts from initial pressure and zero initial pressure derivative, then models propagation in an unbounded domain without explicit boundary conditions.Time-gating is assumed to remove reflections from experimental equipment.Source: §2.1 · pp. 4–5
  • The sensor model first windows the pressure field over the array’s accessible space and finite recording time, then applies a measurement operator M.The work assumes access to the adjoint measurement operator M*.Source: §2.1 · p. 5

2.2. Derivation of the Adjoint PAT Operator

  • The adjoint of the full forward operator factors as A*=P*M*, with P* mapping sensor-space data to a spatial function evaluated from the adjoint wave solution at time T.The adjoint operator is built from the adjoint measurement operator and the adjoint of the windowed pressure-propagation map.Source: §2.2 · pp. 5–7
  • The derivation’s boundary-term argument requires the initial pressure support to lie inside the measurement surface and recording time T to be long enough for waves to leave that domain.The text also notes that finite detection time makes the window vanish at the endpoints.Source: §2.2 · p. 6
  • The adjoint propagation is defined through a wave equation whose time-varying mass source uses the measured data, unlike time reversal, which imposes time-reversed pressure as a boundary constraint.Both constructions evaluate their resulting field at t=T.Source: §2.2 · p. 7

3. Methods and Implementation

The methods implement the PAT operators with a k-space time-domain model and compare iterative time-reversal, least-squares, positivity-constrained, and total-variation reconstruction schemes.

3.1. Numerical Wave Propagation Model

  • The k-space time-domain scheme discretizes the forward, adjoint, and time-reversal operators on a bounded domain enclosed by a PML.The adjoint can also be implemented with other time-domain wave-propagation techniques that support a time-varying mass source term.Source: §3.1 · p. 7

3.2. Image Reconstruction Methods

  • A Neumann-series enhancement of time reversal is reformulated as an iterative image reconstruction algorithm.The study also omits harmonic extension because the chosen final time leaves residual-wave power below sensor noise.Source: §3.2 · p. 8
  • The least-squares reconstruction uses first-order gradient descent, replacing the time-reversal operator with the adjoint and adding a step size.The authors choose this approach for comparison with iterative time reversal, although CGLS is another stated option.Source: §3.2 · p. 8
  • Positivity can be imposed on least-squares reconstruction by re-projecting each gradient-descent update, and a positivity-constrained iterative time-reversal variant is also defined.The latter follows from the similarity between the iterative time-reversal and adjoint-based updates.Source: §3.2 · pp. 8–9
  • The TV method adds a gradient-amplitude penalty to the data-fit objective and solves the positivity-constrained problem with proximal gradient descent.TV(p) is the ℓ1 norm of the gradient-field amplitude; its proximal operator is a positivity-constrained TV-denoising problem.Source: §3.2 · p. 9
  • For the described iterative algorithms, a step size below 2/θ ensures convergence to a minimizer, with θ computed once per experimental setting by power iteration.Here, θ is the largest singular value of A*A.Source: §3.2 · p. 9

4. Computational Studies

The computational studies assess the analytical adjoint in homogeneous 3D and inhomogeneous 2D settings, then examine its precision and its use in inverse reconstruction. Reconstructions show subtle differences among basic methods, while positivity or total-variation constraints bring a more substantial improvement.

4.1. Computational Scenarios

  • The studies use a homogeneous 3D planar-sensor setup and an inhomogeneous 2D setup with irregular sensors, applying simple point sampling in both.The 3D setup models planar Fabry-Pérot scanners; the 2D medium properties vary over a range encountered in the human breast.Source: §4.1 · p. 9

4.2. Precision Assessment

  • The analytical adjoint is sufficiently close to the algebraic adjoint for variational reconstruction; iterate ℓ2 error stays nearly constant while ℓ∞ error grows slowly without indicating instability.The adjoint discrepancy decreases slightly as N decreases, a trend that can be offset by increasing PML size; single and double precision show no observable difference in that measure.Source: §4.2 · p. 10, Figure 2 · p. 12

4.3. Inverse Reconstructions

  • In the limited-view 3D scenario, positivity or total-variation constraints improve reconstruction quality more substantially than differences among time reversal, back-projection, and iterative approaches.The limited-view geometry also changes intensity decay with depth and produces circular limited-view artifacts.Source: §4.3 · pp. 12–13, Figure 3 · p. 13

5. Discussion, Conclusions and Outlook

The discussion presents the continuous PAT adjoint, explains its numerical implementation and reconstruction use, and identifies assumptions and directions for further work.

5. Discussion, Conclusions and Outlook

  • The continuous adjoint is connected to a wave equation with a time-varying mass source, and the proof uses a measurement window and the Green’s function representation.The window makes boundary terms from integration by parts in time vanish.Source: §5 · p. 13
  • The derivation covers a more general choice of sensor region Γ and recording time T than the standard boundary-measurement setting.In the standard setting, the boundary terms vanish only if p0 is supported inside the measurement surface and T is long enough for all waves to leave the domain.Source: §5 · pp. 13–14
  • The adjoint-then-discretize approach is sufficiently precise compared with the scheme-specific algebraic adjoint, which must be derived separately for each wave-propagation code.The algebraic approach provides a discrete, purely algebraic characterization of the adjoint.Source: §5 · p. 15
  • The adjoint can be used in variational reconstruction problems that incorporate a-priori constraints on the initial pressure.Source: §5 · p. 15
  • The paper contrasts its approach with an alternative adjoint based on a wave-transmission problem solved by a computationally intensive BEM-FEM method [5].That work focuses on short recording times and full boundary measurements.Source: §5 · p. 15
  • Future study is needed on how time reversal and the adjoint differ across sensor geometries, near-sensor targets, limited coverage, and short recordings.The distinct artifacts can cause different deficiencies in in-vivo image interpretation.Source: §5 · p. 15
  • The reconstruction schemes used here are computationally inefficient and untuned, so example iteration counts do not indicate expected computational cost.More sophisticated optimization methods are proposed for future work, especially for large-scale experimental models.Source: §5 · p. 15
  • Example scripts to reproduce the results are planned for a PAT reconstruction toolbox intended to complement k-Wave.Source: §5 · p. 15

Appendix A. k-Space Time Domain Method

The appendix describes the k-space time-domain collocation scheme and its discretization of wave propagation and measurement operators.

Appendix A. k-Space Time Domain Method

  • The scheme embeds the spatial domain in a rectangular volume discretized on a regular grid, with a PML absorbing boundary to damp outgoing waves.Time is discretized over [0, T] with steps Δt = T/Nt.Source: Appendix A · p. 16
  • The k-space method interpolates between collocation points using a truncated Fourier series, enabling efficient gradient calculations with the FFT.Source: Appendix A · p. 16
  • The wavenumber-domain correction κ compensates for the finite-difference time-derivative approximation; it is exact in homogeneous media and reduces errors such as numerical dispersion in heterogeneous media.Source: Appendix A · p. 16
  • The implementation uses a homogeneous reference sound speed chosen for stability and interleaves time updates with staggered-grid gradient calculations.Grid staggering incorporates spatial translations of ±Δξ.Source: Appendix A · p. 16
  • The discrete window operator W implements the measurement window function ω used in the continuous model.Source: Appendix A · p. 16

Appendix B. Pseudospectral Implementation of the PAT Operators

The appendix outlines how the k-Wave iteration implements discrete forward and adjoint PAT operators, including their scaling and relation to time reversal.

Appendix B. Pseudospectral Implementation of the PAT Operators

  • The forward operator uses the k-Wave iteration, applies smoothing Q to the initial pressure to suppress high spatial frequencies, and maps the measured pressure time series through M.In the simulations, M is the identity matrix.Source: Appendix B · p. 17
  • The adjoint implementation uses M* on the data and reverses the sensor data in time before the iteration, then applies Q to the terminal pressure.The self-adjoint smoothing Q is needed for adjointness with the forward implementation.Source: Appendix B · p. 17
  • The adjoint output is rescaled by c0^2 and ρ̃0 because the implementation solves the adjoint first-order system rather than the second-order wave equation.The time-reversal implementation instead introduces reversed sensor pressure as a Dirichlet source and omits the adjoint step that sets sensor density to zero.Source: Appendix B · p. 17

Appendix C. Discrete Total Variation Functional

The appendix gives the finite-difference discretization of total variation for 2D images with Neumann boundary conditions and notes that the 3D implementation is analogous.

Appendix C. Discrete Total Variation Functional

  • The 2D total-variation seminorm uses forward finite differences with Neumann boundary conditions, implemented by copying boundary pixel values beyond the image edges.The 3D implementation works in the same way.Source: Appendix C · p. 18
Loading 1602.02027v2…