A Reciprocity-Law-Compliant Photoacoustic Forward-Adjoint Operator

1. Problem Statement

Photoacoustic Tomography (PAT) is a hybrid imaging modality that combines optical contrast with acoustic resolution. The core mathematical challenge lies in solving the acoustic inverse problem: reconstructing the initial pressure distribution ($p_0$) inside a tissue from time-series pressure measurements recorded on a boundary surface.

While various reconstruction algorithms exist (e.g., time-reversal, backprojection), many lack a rigorous mathematical framework that strictly adheres to the reciprocity law of physics when implemented in a discrete, computational setting. Specifically, defining a forward operator (mapping source to boundary data) and its corresponding adjoint operator (mapping boundary data back to the source) that satisfy an exact inner-product relation is difficult, especially when dealing with:

• Singularities: The instantaneous nature of the photoacoustic source and the restriction of fields to a boundary surface.
• Discretization: The mismatch between continuous physical laws and discrete grid-based numerical simulations (e.g., receivers not aligning perfectly with grid points).
• Reciprocity: Ensuring the adjoint operator is the true mathematical transpose of the forward operator, which is crucial for iterative optimization methods (like gradient descent) to converge correctly.

2. Methodology

The author extends a previously proposed acoustic forward-adjoint framework to the specific context of PAT. The methodology involves three main components:

A. Mathematical Formulation

The PAT process is modeled as a composite operator $A = A_{acoustic} \circ Q$:

1. Source Operator ($Q$): Maps the initial pressure distribution $p_0$ to a temporally regularized source term $s_Q$. Since the physical source is instantaneous (Dirac delta in time), it is regularized using a smoothed temporal delta function $\delta_q(t)$ to avoid numerical instability.
2. Acoustic Forward Operator ($A_{acoustic}$): Solves the wave equation in free space with the source $s_Q$. Crucially, it employs a reception operator ($R^D_{G,r}$) that maps the continuous pressure field in free space to the boundary.
   • This reception operator is defined as a left-inverse of an emission operator.
   • It treats the boundary restriction as a "smeared singularity" in the normal derivative of the field, regularized using a spatial regularized Dirac delta distribution ($\delta_b$).
   • A geometric scaling factor ($a_R$) is introduced to account for the Jacobian between volumetric and surface elements, ensuring the operator is scale-invariant.

B. The Adjoint Operator

Using the inner-product definition, the adjoint of the composite operator is derived as $A^* = Q^* \circ (A_{acoustic})^*$.

• Adjoint Source ($Q^*$): Derived via integration by parts, it acts as a time-reversed operation on the test function, evaluated at $t=0$.
• Adjoint Acoustic Operator: Corresponds to a time-reversed double-layer potential integral. In the discrete domain, this is implemented as a source term in the adjoint wave equation, where the boundary data acts as a force source distributed over the receiver surface using the regularized delta function.

C. Discretization and Numerical Implementation

• Grid: A regular Cartesian grid with staggered time steps is used.
• Wave Solver: A k-space pseudospectral method is employed to solve the first-order coupled wave equations (velocity, density, pressure) with high accuracy.
• Boundary Handling: Perfectly Matched Layers (PML) are used to simulate free-space propagation on a finite grid.
• Regularization: Both the temporal source and the spatial boundary restriction are regularized using sinc-based Dirac delta approximations. A thresholding technique is applied to the delta function to reduce computational cost by ignoring points far from the receiver nodes.
• Optimization: An iterative projected gradient descent algorithm is used to minimize a scaled least-squares objective function, reconstructing $p_0$ from boundary data.

3. Key Contributions

1. Reciprocity-Compliant Framework: The paper defines a forward-adjoint pair that strictly satisfies the reciprocity law. The reception operator is mathematically proven to be the left-inverse of the emission operator, ensuring the adjoint is the true transpose.
2. Regularized Singularity Handling: The authors introduce a rigorous method to handle singularities arising from:
   • The instantaneous photoacoustic source (time regularization).
   • The restriction of the field to a boundary surface (spatial regularization via $\delta_b$).
3. Scaled Operator Formulation: By introducing a scaling factor ($1/a_R$), the authors eliminate the dependency of the adjoint operator on the bandwidth parameter of the regularization kernel, resolving issues of non-uniqueness found in previous formulations.
4. Generalization to Off-Grid Receivers: The framework is explicitly designed to handle receivers that do not align with the computational grid points (off-grid), using the regularized delta distribution to interpolate fields between the grid and the physical receiver surface.

4. Results

The study validates the methodology through numerical experiments in 2D:

• Inner-Product Consistency Test:

  • On-Grid Receivers: When receivers align with grid points, the relative difference (RD) between the left-hand side and right-hand side of the inner-product identity was $9.32 \times 10^{-4}\%$.
  • Off-Grid Receivers: When receivers were placed on a circular boundary (not aligned with grid points), the RD was $2.07 \times 10^{-5}\%$.
  • Significance: These near-zero errors confirm the mathematical consistency of the derived forward-adjoint pair.

• Image Reconstruction:

  • A vessel phantom was used as the ground truth.
  • The reconstruction was performed using an iterative minimization framework with the proposed operators.
  • Avoidance of Inverse Crime: The simulation (data generation) and reconstruction used different grid resolutions and receiver discretizations to ensure realistic validation.
  • Outcome: The algorithm converged monotonically (25 iterations). The reconstructed image showed high qualitative and quantitative accuracy, closely matching the ground truth with low relative error.

5. Significance

This work provides a robust theoretical and computational foundation for iterative PAT reconstruction.

• Algorithmic Stability: By ensuring the forward and adjoint operators are true mathematical transposes, gradient-based optimization methods are guaranteed to follow the correct descent direction, leading to stable and accurate convergence.
• Physical Fidelity: The strict adherence to the reciprocity law ensures that the numerical model respects fundamental physical principles, reducing artifacts in the reconstructed images.
• Practical Applicability: The ability to handle off-grid receivers and arbitrary boundary geometries makes this framework suitable for realistic clinical imaging setups where transducer placement is rarely perfectly aligned with a computational grid.
• Future Impact: The authors note that this framework is a stepping stone toward clinical applications, offering a path to higher-resolution and more accurate photoacoustic imaging.