PaperlayerView exact revision on arXiv

Source-linked AI summary

Scalable photoacoustic tomography implementations accounting for the spatial impulse response of transducers

Trung-Thai Do, Paul Escande, Caroline Chaux, Jérôme Gateau, Hwee Kuan Lee

arXiv:2609.29253v1eess.IV

TL;DR

Three-dimensional iterative photoacoustic reconstruction requires repeated forward and adjoint evaluations that are costly at realistic scale, while finite transducer surfaces affect resolution. The paper develops GPU-oriented matrix-free operators with an analytical treatment of detector surfaces. These operators outperform point quadrature in accuracy and speed, and the surface method outperforms back-projection in a large-scale reference reconstruction.

  • Problem

    Large three-dimensional model-based reconstruction requires repeated forward and adjoint evaluations, but realistic systems make these computations and their memory demands costly, while ignoring transducer SIR reduces resolution.

    Source: §1 · pp. 2, 4

  • Method

    The authors represent initial pressure with radial functions and develop matrix-free GPU operators that account for transducer surfaces, including a closed-form surface-area treatment.

    Source: §1 · p. 6, §3.1 · p. 12, §6 · p. 38

  • Results

    The surface method outperforms point-based quadrature in accuracy and speed, and outperforms back-projection on the large-scale reference system.

    Source: §6 · p. 38

  • Takeaways & Limitations

    The analytical surface method avoids detector-surface discretization and brings model-based 3D reconstruction at real-system scale within a single GPU’s capabilities.

    Source: §6 · p. 38

    The reconstruction formulation is deliberately simple and is not intended to provide state-of-the-art reconstruction quality.

    Source: §5.1 · p. 28

Abstract

from arXiv · show

Iterative model-based reconstruction in photoacoustic tomography repeatedly applies the forward operator mapping the initial pressure to the transducer signals, and its adjoint. At the scale of current three-dimensional systems, this operator cannot be stored and must be evaluated matrix-free, while accounting for the finite, focused surface of the transducers, whose spatial impulse response degrades the resolution when ignored. Representing the initial pressure by compactly supported radial functions, we show that the measured signal is exactly a temporal convolution between a system kernel gathering the radial function and the electrical impulse response, and a purely geometric quantity accounting for the portion of the transducer surface reached by the wave emitted from a voxel during one time step. Two implementations are proposed, differing only in how this quantity is evaluated: a quadrature over points of the surface, as in existing works, or a closed-form area, which never discretizes the surface. We derive closed forms for planar and cylindrically focused transducers and provide, in the latter case, two accelerations of the resulting elliptic integrals, a lookup table and a trapezoidal approximation, together with the piecewise planar approximation customary in the literature. These implementations reduce the per-voxel geometric computations and are released as an open-source Python package for graphics processing units. The performance of these operators is first demonstrated on a synthetic phantom, where the lookup-table-based operator reaches the accuracy of the exact evaluation ten times faster and outperforms the point discretization on both accuracy and runtime. A second experiment shows that they enable the processing of a realistic vascular phantom at full scale, with a higher peak signalto-noise ratio and a better resolution than the back-projection counterpart. The released implementations are an important step towards the adoption of three-dimensional modelbased photoacoustic reconstructions with finite and focused transducers.

1 Introduction

Iterative three-dimensional PAT requires repeated matrix-free applications of a forward operator and its adjoint, while finite transducer surfaces and their spatial impulse responses must be modeled accurately. The paper develops GPU-oriented operators using radial-basis representations, including an analytical surface treatment, and validates them on synthetic and realistic systems.

1 Introduction

  • Finite-aperture transducers integrate pressure over their surfaces, and ignoring their spatial impulse response causes a systematic loss of resolution.Cylindrically focused surfaces are also used to benefit from acoustic focusing gain.Source: §1 · p. 2
  • Iterative regularized reconstruction can outperform analytical inversion in 3D, but depends on efficient forward and adjoint evaluations.Analytical inversion is fast but relies on idealized detectors and canonical full-view geometries; full-wave solvers can be costly for repeated evaluations.Source: §1 · p. 3
  • Dense quadrature over curved detector surfaces is a computational bottleneck, while existing radial-basis approaches either approximate surfaces or lack quantified error against the continuous model.The paper identifies the absence of a method combining closed-form treatment of extended curved detectors, quantified approximation error, and GPU-friendly arithmetic at modern 3D scale.Source: §1 · pp. 4–5
  • The framework derives matrix-free forward and adjoint operators from the wave equation, incorporating transducer-surface integration, electrical impulse response, and compactly supported radial bases.It recovers point discretization of the detector surface as one implementation.Source: §1 · p. 5
  • The surface implementation expresses measurements as a temporal convolution with coefficients representing detector area between successive wavefront spheres, with closed forms for planar and cylindrical detectors.The paper also establishes error bounds for both implementations; the point method includes a quadrature term that the analytical surface treatment removes.Source: §1 · p. 5
  • Both implementations reduce to per-voxel computations and batched convolutions, making them suitable for GPU architectures.They are released in a Python package, and the large-scale experiment demonstrates reconstruction at the scale of a real system.Source: §1 · pp. 5–6
  • The surface method outperforms point-based quadrature in accuracy and speed, and enables large-scale model-based reconstruction that outperforms back-projection on the reference system.The reduced-size experiment mimics a realistic cylindrical-transducer acquisition; the large-scale demonstration accounts for transducer SIR.Source: §1 · p. 6

2 Modeling of PAT systems

This section specifies the acoustic and measurement model, discretizes initial pressure with radial functions, and reformulates propagation to separate the radial profile from geometry. It also states the modeling assumptions and reference acquisition parameters used in the operator analysis.

2.1 Notations

  • The notation fixes boldface for vectors and matrices, flattens grid indices, and uses standard Lebesgue and Sobolev spaces for analysis.The index set [n] denotes the first n integers.Source: §2.1 · p. 6

2.2 Optical generation of ultrasound waves

  • The sample occupies a bounded box, and its pressure propagates through a homogeneous, lossless medium at sound speed c.The wave field follows the wave equation, whose solution is given by Kirchhoff’s and Poisson’s formulas.Source: §2.2 · pp. 7–8

2.3 Pressure field measurement

  • Each transducer measures pressure integrated over its surface, followed by the transducer’s electrical impulse response and digital sampling.The EIR is modeled as a band-pass filter; the digitized measurements form a vector across transducers and time samples.Source: §2.3 · p. 8

2.4 Discretization of p0

  • The initial pressure is represented by a Kansa expansion of radial functions centered on a uniform Cartesian grid covering the sample.The unknown is a coefficient vector, and the relationship between coefficients and nodal values depends on the radial function.Source: §2.4 · p. 8

2.5 Setting

  • The analysis assumes a compactly supported bounded radial function, a compactly supported integrable EIR, and separation between the sample support and transducers.The EIR’s practical truncation is taken long enough that its truncation error can be neglected.Source: §2.5 · p. 9

2.6 Decoupling interpolation and propagation

  • Radial symmetry lets propagation from each grid point be corrected by a direction-independent convolution kernel, decoupling the basis profile from wave propagation.The kernel structure supports parallel arrival-time calculations and efficient convolution; angular dependence would undermine memory coalescing.Source: §2.6 · pp. 9–10

2.7 Reference experimental setup

  • The reference acquisition uses 8·10^6 grid points and 11,520 signals, making a dense forward matrix impractical to store.The sample grid has 50 µm spacing, and each signal has about 1,000 samples; dense storage would require about 740 TB in double precision.Source: §2.7 · p. 10

2.8 Choice of the radial function

  • For non-interpolatory radial functions, coefficients require solving a Gram-matrix system; compact support makes that matrix sparse under the stated spacing condition.The band-pass EIR guides support-radius choice, while stationary spacing reduces grid size but lacks convergence guarantees for the initial pressure itself; a complete error theory for that regime is outside the paper’s scope.Source: §2.8 · pp. 10–11

3 Implementation of the forward operator

This section describes a shared discrete-convolution implementation of the forward operator that separates geometric propagation from radial-function and electrical-impulse-response effects. It compares surface-quadrature and analytical approaches to computing geometric coefficients, including the accuracy–cost tradeoffs associated with upsampling and coarser area evaluation.

3.1 Overview

  • The operator encodes propagation in geometric coefficients and applies the radial-function and electrical impulse response through a system-kernel convolution, with sampling handled by convolution and decimation.The coefficients are represented on an upsampled time grid; the kernel is reused across transducers, and the geometric-coefficient computation is the dominant cost in numerical experiments.Source: §3.1 · pp. 12–14
  • The shared convolution and sampling stage uses a reusable kernel and can exploit FFTs and batched transducer convolutions, while the adjoint reuses the geometric routines.The convolution complexity is O(UL log(UL)); the adjoint is described as having the same cost as the forward operator.Source: §3.1 · pp. 13–14

3.2 Point detector approximation

  • The point-detector method approximates the surface integral with a quadrature over Q points, accumulating voxel contributions into time bins at O(NQ) cost.Its error has a Q^-1 quadrature term and a (UFs)^-1 arrival-time quantization term; increasing the user-set upsampling factor U can improve accuracy.Source: §3.2 · pp. 14–16

3.3 Exploiting the surface of the transducer

  • The surface-based method analytically integrates the surface contribution, avoiding surface discretization by computing the area reached by each voxel’s wave during each time bin.Its error arises from arrival-time quantization and approximating inverse distance within each bin; the stated bound applies when c(UFs)^-1 ≤ distmin.Source: §3.3 · pp. 16–18
  • The user-defined upsampling factor U can improve accuracy, but raises complexity; choosing Δl between 1 and U preserves fine-grid arrival-time resolution while introducing approximation error.Δl = 1 recovers Algorithm 3 exactly, while Δl = U evaluates areas at the native sampling frequency; the latter gave the best trade-off in the reported experiments, though this may be system dependent.Source: §3.3 · pp. 18–19

3.4 Comparison of the two implementations

  • Point-detector quadrature applies to any transducer geometry, whereas the surface method requires geometry-specific analytical area formulas and a more involved implementation.The paper derives these formulas for planar and cylindrical transducers.Source: §3.4 · p. 19
  • Point-method coefficient computation costs O(NQ), independent of U, while surface-method cost O(NB CM) grows linearly with U and can be mitigated by coarser area evaluation.The average number of time intervals B also grows linearly with U; evaluating areas every Δl intervals reduces area evaluations by a factor of Δl.Source: §3.3 · pp. 18–19, §3.4 · p. 19
  • Under the same smoothness assumption, the point method’s convergence rate is significantly slower, dropping to Q^-1/2; the Q^-1 rate in Proposition 2 requires stronger smoothness than that in Proposition 3.Source: §3.4 · p. 19
  • Both implementations have identical memory use, dominated by FFT vectors and workspaces at about 4KUL floating-point values, in addition to O(N) values for p0.The point method’s Q^-1 convergence bound requires stronger smoothness than the surface method’s bound.Source: §3.4 · p. 19

4 Implementations for usual transducer geometries

This section derives closed-form geometric area evaluations for planar and cylindrical transducers. For cylindrically focused transducers, it develops and compares lookup-table and trapezoidal approximations of elliptic-integral evaluation with a piecewise-planar surface approximation.

4 Implementations for usual transducer geometries

  • A shared separable-surface framework derives the geometric areas by splitting the parameter domain into quadrants and evaluating its level-set regions.The framework applies to surfaces with the specified separable form and uses geometry-dependent inputs to compute each area.Source: §4 · pp. 20, 21, 23

4.2 Cylindrical surface transducer

  • For cylindrical transducers, the exact area expression uses incomplete elliptic integrals, with lookup-table, trapezoidal, and piecewise planar alternatives.The lookup-table and trapezoidal approaches reduce the cost of evaluating the elliptic integral; the piecewise planar approximation uses the planar closed-form solution.Source: §4 · p. 23, §4.2.1 · p. 26

4.1 Plane surface transducer

  • For planar transducers, the framework yields closed-form area calculations using an elementary primitive and instantiates Algorithm 4.The planar surface’s area integral is evaluated through an elementary primitive.Source: §4.1 · p. 24

4.2.1 Elliptic integral

  • Each cylindrical time-step integral requires 2 to 8 evaluations of the elliptic integral, making it the method’s primary computational bottleneck.The algorithm remains valid when a voxel lies on the cylinder’s focal line, where the arccos is never evaluated.Source: §4.2.1 · p. 26
  • The lookup table makes the Exact evaluation ten times faster with no measurable loss of accuracy.It uses 1000 × 1000 non-uniform samples concentrated near the singularity; locating adjacent grid points requires two square-root evaluations.Source: §4.2.1 · p. 27
  • The trapezoidal approximation is the fastest of the three proposed methods, avoids elliptic integrals, and introduces approximation error.It approximates the curve by a segment and reduces each subdomain computation to endpoint evaluations requiring two square-root operations.Source: §4.2.1 · p. 27

4.2.2 Piecewise-Planes

  • The piecewise-planar approximation replaces the cylindrical surface’s swept area with a sum over tangent planes, trading computational cost against curvature-approximation error.Increasing the number of planar sectors raises computational cost while improving how closely the planes represent the cylinder’s curvature.Source: §4.2.2 · p. 27

5 Numerical experiments

The numerical experiments compare reconstruction quality and computational cost across forward-operator approximations and test the LUT operator at full scale against a system-tailored back-projection method. The full-scale vascular-phantom evaluation examines vessel-focused image quality and resolution using masked metrics, Fourier shell correlation, and qualitative vessel recovery.

5.1 Common methodology

  • The shared reconstruction protocol uses a deliberately simple formulation to limit regularization effects, so differences can be attributed mainly to each forward operator.Measurements are synthesized with a very fine Exact operator, withheld from reconstruction to avoid the inverse crime; accuracy is assessed with PSNR and SSIM.Source: §5.1 · pp. 28–29

5.2.1 Setup

  • The experiment considers rotational acquisition around the focal point using a 7 × 7 absorbing-rod phantom, reconstructed on fine 81×81×15 and coarse 27×27×5 grids under noiseless and 1% Gaussian-noise conditions.Sources near the rotation axis are weakly affected by the SIR, with stronger effects as radial distance increases.Source: §5.2 · p. 29, §5.2.1 · pp. 30–31

5.2.2 Results

  • On the noiseless fine grid, Trapezoidal, Point, and Piecewise-Planes PSNR decreased by 0.3 dB, 1.2 dB, and 0.7 dB, respectively, versus Exact and LUT.Point and Piecewise-Planes overestimated the central rod’s amplitude despite locating it accurately.Source: §5.2.2 · p. 31, Table 2 · p. 32, Figure 8 · p. 33
  • Across experiments, LUT matches Exact accuracy while ranking second-fastest; Trapezoidal is fastest but can reduce quality, and Piecewise-Planes is costly without leading in quality.The accuracy–cost comparison is presented for the controlled study; the prose reports similar quality across operators except on the noiseless fine grid.Source: §5.2.2 · p. 32, Table 2 · p. 32, Figure 9 · p. 34

5.3.1 Setup

  • The full-scale experiment models a limited-view kidney-vessel acquisition on 8 million voxels and compares LUT reconstruction with a tailored back-projection baseline.Exact reference synthesis takes about 2 h for one application, making it prohibitive inside an iterative solver; the simulated measurements include 1% Gaussian noise.Source: §5.3.1 · pp. 34–35

5.3.2 Results

  • For the full-scale experiment, Figure 11 presents side-by-side maximum-intensity projections from the ground truth, back-projection, and LUT reconstructions.The supplied passage describes the comparison display but does not report its quantitative or qualitative outcome.Source: §5.3.2 · p. 35
  • Model-based reconstruction produces sharper vessels, preserves vessel continuity, and partially recovers small-z vessels that BP misses, while larger-z vessels remain barely recovered.Its maximum intensity is ≈1.1 for a ground truth of 1; BP amplitudes are not directly comparable because BP restores relative rather than absolute amplitudes.Source: §5.3.2 · p. 35, Figure 11 · p. 36
  • On unmasked metrics, model-based reconstruction has higher PSNR than BP (27.18 vs. 25.46 dB), but BP has higher SSIM (0.849 vs. 0.657).The phantom is 99.7% background, so these metrics are dominated by background voxels; BP metrics use a least-squares scale adjustment, while model-based metrics are unscaled.Source: §5.3.2 · pp. 35–36
  • On a vessel-focused mask, the model-based reconstruction outperforms BP in both PSNR (14.45 vs. 10.65 dB; +3.8 dB) and SSIM (0.321 vs. 0.029).The mask unions ground-truth support with reconstructed voxels above 0.1, then dilates by a ball of radius 3 voxels.Source: §5.3.2 · p. 37
  • Model-based reconstruction takes 15.8 h for 30 iterations on one NVIDIA A100 GPU, about 32 min per iteration, compared with 3 s for BP.Source: §5.3.2 · p. 37
  • At the 1-bit FSC threshold, model-based reconstruction reaches ζ ≈0.34 versus ζ ≈0.085 for BP, a resolution gain of about 4× (≈10 MHz vs. ≈2.6 MHz).The model-based result remains correlated with ground truth down to a spatial period of approximately 150 µm, or three voxels, the finest scale in the phantom.Source: §5.3.2 · p. 37
  • Accounting explicitly for detector SIR improves reconstruction quality but raises computational requirements; this is the first demonstrated model-based reconstruction on such a large-scale system.The demonstration was enabled by the efficient GPU implementations developed in Sections 3 and 4.Source: §5.3.2 · p. 37

6 Conclusion

The paper presents matrix-free forward and adjoint operators that account for transducer SIR with memory costs scaling to full-resolution 3D systems. Its analytical surface treatment outperforms point-based quadrature and supports large-scale model-based reconstruction with SIR.

6 Conclusion

  • The proposed forward operator and its adjoint are fast, matrix-free, account for transducer SIR, and have memory costs that scale to full-resolution 3D systems.Source: §6 · pp. 37–38
  • The two implementations use a Kansa-type radial discretization and differ in their surface integral treatment: point quadrature versus analytical treatment.The analytical method is instantiated for planar and cylindrical elements, with a methodology for extending it to other surface geometries.Source: §6 · p. 38
  • The shared mathematical framework establishes error bounds relative to the continuous wave-equation model and analyzes each implementation’s computational complexity.This makes the accuracy–computation trade-off explicit.Source: §6 · p. 38
  • The analytical surface method outperforms point-based quadrature in both accuracy and computational speed, and enables large-scale model-based reconstruction accounting for transducer SIR.The authors describe this as the first reconstruction demonstrated at this scale that accounts for transducer SIR.Source: §6 · p. 38
  • The released implementation brings model-based 3D photoacoustic reconstruction at the scale of a real acquisition system within the capabilities of a single GPU.The analytical surface method avoids surface discretization while providing improved accuracy and computational efficiency.Source: §6 · p. 38

Funding

The authors acknowledge institutional, national, and international funding, as well as access to high-performance computing resources.

Funding

  • The work received support from CNRS MITI, GENCI through IDRIS HPC access, and the PHC Merlion partner ministries and A*STAR.The IDRIS allocation was 2025-AD011014692R2.Source: Funding · p. 39

A Table of main notations.

The notation table collects symbols for the model, surface methods, transducer geometries, and numerical experiments.

A Table of main notations.

  • The table defines grid and radial-basis notation, including grid points covering the domain, the supported radial function, and initial-pressure coefficients.It also lists the minimum grid-to-transducer distance.Source: A Table of main notations · p. 40
  • Model and convolution notation covers the system kernel, arrival-time distributions, upsampling, computed measurements, and the forward and adjoint operators.Source: A Table of main notations · p. 40
  • The surface-method notation identifies quadrature points and areas, the surface portion reached during a time interval, and geometric parameters used in surface integration.The table also lists the piecewise-plane count and parameters for planar and cylindrical transducers.Source: A Table of main notations · p. 40
  • The remaining entries define surface parametrization and geometry variables, followed by notation for numerical experiments and reconstruction parameters.These include lookup-table sampling counts, the number of piecewise planes, the ground-truth phantom, and the regularization parameter.Source: A Table of main notations · p. 40

B Proof of Proposition 1

The proof reduces the relevant surface integral to a one-dimensional integral and uses that reduction to establish Proposition 1. The appendices also state the mathematical framework and distribution notation used in the justification.

B Proof of Proposition 1

  • The proof substitutes the discrete initial-pressure model into the propagation expression, evaluates the resulting contribution, and concludes with the announced identity.The stated distance assumption makes one term vanish in the derivation.Source: B Proof of Proposition 1 · pp. 41–42
  • The proof parametrizes the wavefront sphere and changes variables using the distance to the grid point, reducing the surface integral to a one-dimensional integral.The distance is independent of azimuth, and the radial function’s support and oddness properties are used in the ensuing integration.Source: B Proof of Proposition 1 · pp. 41–42

C Mathematical justification of Section 3

C.1 Quick reminders on distributions

  • The reminders define distributions as continuous linear forms on smooth compactly supported test functions and introduce the pairing and Dirac mass notation.They also describe convolution with distributions and weak integration of parameter-dependent distributions.Source: §C.1 · pp. 42–43

C.2 Justification of the derivations

This section justifies the derivations in the distributional sense and establishes that the sampled signal can be computed by discrete convolution followed by downsampling. It covers both quadrature-based and surface-integral representations.

C.2 Justification of the derivations

  • The measure formulation covers both quadrature over surface points and integration with the surface measure.The measure is positive, supported on Σ, and has total mass |Σ|.Source: §C.2 · p. 43
  • The derivation establishes equality in D′, which holds almost everywhere and everywhere when g is continuous.The proof uses integrability of g and establishes the required distributional construction.Source: §C.2 · pp. 43, 45

D Discrete convolutions

  • Convolving a piecewise constant function with a kernel and sampling it is computed exactly by discrete convolution followed by downsampling.The discrete-convolution derivation applies to arbitrary functions.Source: D Discrete convolutions · p. 45
  • On the upsampled grid, the signal coefficients are h ⋆ f; the original-grid coefficients follow by downsampling by U.The kernel h is compactly supported and belongs to L1(R).Source: D Discrete convolutions · pp. 45–47

E.1 Overview of the proofs

The proof overview decomposes each method into successive approximations and bounds their errors separately. The point-detector and surface methods differ in which stages they perform exactly.

E.1 Overview of the proofs

  • The point-detector method first approximates the surface integral by quadrature, then approximates the resulting function as piecewise constant.For this method, the piecewise-constant approximation introduces no further error in the stated decomposition.Source: §E.1 · p. 47
  • The surface method evaluates the surface integral analytically, then approximates the piecewise-constant coefficients using a midpoint rule.The diagram’s error bounds sum contributions along each method’s approximation chain.Source: §E.1 · p. 47
  • Both methods use convolution identities to express the signal through h and their respective approximations of fΣ.The identities follow from associativity of convolution, and the discrete convolution is exact.Source: §E.1 · p. 48
  • The error in each method is bounded independently using a decomposition of the signal difference.The decomposition is evaluated on a finer grid.Source: §E.1 · p. 48

E.2 Proofs

The proofs attribute each method’s error to specific approximation steps and bound those terms under stated regularity and distance assumptions. Supporting lemmata establish the required bounds and regularity properties.

E.2 Proofs

  • For the point-detector method, err3 is zero; its remaining errors arise from the other terms in the decomposition.The proof bounds the quadrature error using smoothness over surface cells and treats arrival-time quantization separately.Source: §E.2 · pp. 48–50
  • The point-detector error analysis uses Taylor expansion over convex surface cells, with a mean-value bound when h is only in W 1,∞(R).The Taylor-based bound uses derivatives of the surface integrand; the weaker regularity case uses the mean value inequality.Source: §E.2 · pp. 49–50
  • For the surface method, err1 is zero; the remaining terms arise from arrival-time quantization and approximation of the inverse distance.The proof combines the bounds on err2 and err3.Source: §E.2 · pp. 50–51

E.3 Supporting lemmata

  • The supporting lemmata establish regularity and bounds for g and the surface integrand under assumptions including positive minimum distance and regularity of h.The proofs also use the padding convention to keep arrival times within the upsampled grid.Source: §E.3 · pp. 51–53

F Elliptic integral evaluation of the arc integral

This section derives a closed form for the cylindrical transducer’s arc integral using elliptic integrals. It also identifies when a reciprocal-modulus transformation is needed to evaluate them.

F Elliptic integral evaluation of the arc integral

  • The closed form for the cylindrical transducer’s arc integral uses a separable function ψ(u) = ρ2 + R2 − 2Rρ cos u.The derivation reduces the integral to elliptic-integral expressions.Source: F Elliptic integral evaluation of the… · p. 54
  • The derivation’s elliptic-integral expressions are well-defined under the stated bound on ν(r) sin2 t.The proof obtains the bound from ψ(u) ≤ r2 on the integration interval.Source: F Elliptic integral evaluation of the… · p. 54
  • When r < R + ρ, ν(r) exceeds 1, so the reciprocal-modulus transformation is required before evaluating the elliptic integrals.This is the regime before the wavefront reaches the far edge of the circle of radius R in the plane of xi.Source: F Elliptic integral evaluation of the… · p. 55

G Adjoint operators

The adjoint operators transpose the forward operator’s linear stages in reverse order, preserving the inner-product identity up to floating-point accuracy. Their geometry stage gathers contributions, and the coarser area grid retains the forward pass’s reduction in area evaluations.

G Adjoint operators

  • The adjoints transpose the forward stages in reverse order and satisfy the inner-product identity up to floating-point accuracy.The forward operator’s stages are geometric assembly, convolution, and downsampling.Source: G Adjoint operators · p. 55
  • Transposing downsampling replicates each time sample across the U upsampled indices it averages, at O(UL) operations.Source: G Adjoint operators · p. 55
  • The adjoint convolution is correlation with the reversed kernel, computed by the same batched FFT at O(UL log(UL)) operations per transducer.Because the kernel is real-valued, reversal corresponds to conjugating its spectrum.Source: G Adjoint operators · p. 55
  • The adjoint changes geometric assembly from scattering voxel contributions to gathering time samples, while reusing the forward geometric quantities and routines.The padding convention keeps scatter and gather indices aligned; the Piecewise-Planes operator sums tangent-plane contributions in both directions.Source: G Adjoint operators · p. 55
  • With ∆l = U, the forward spreads each area over ∆l time intervals, and the adjoint gathers those contributions with the same reduction in area evaluations.Source: G Adjoint operators · p. 55
  • The adjoint matches the forward operator’s complexity and memory footprint; voxel loops avoid write conflicts, while summing over K transducers requires a reduction.Its difference from the forward operator is its access pattern: it accumulates into the output entry z[i].Source: G Adjoint operators · p. 56
Loading 2609.29253v1…