Hessian-free Ray-Born Inversion for Quantitative Ultrasound Tomography of Weakly Heterogeneous Media

TL;DR

The paper introduces a Hessian-free, ray-Born inversion framework for quantitative ultrasound tomography in weakly heterogeneous media, formulated in the frequency domain. By diagonalizing the Hessian with a tailored weighting and employing a paraxial ray-tracing system, it achieves single-step subproblem solves and substantial computational savings while embedding regularization into the forward model. The method leverages a ray-based Green’s function in a lossy Helmholtz setting, incorporating absorption and dispersion to produce robust, high-resolution sound-speed reconstructions that outperform prototype Born and Hessian-based approaches, particularly under noise. Numerical experiments on 2D breast-like phantoms confirm accurate image reconstruction and improved stability, supporting translation toward 3D clinical QUT applications with reduced computational cost.

Abstract

This study presents a frequency-domain, Hessian-free ray-Born inversion method for quantitative ultrasound tomography, extending the author's previous Hessian-based approach. Both approaches model acoustic wave propagation using a ray-based approximation of the Green's function in smoothly varying heterogeneous media, and perform the inversion iteratively in the frequency domain, progressing from low to high frequencies. In the earlier method, each frequency subproblem was solved through iterative inversion of the Hessian matrix, a process that not only increased computational cost but also made the update steps more sensitive to noise. The present work addresses these limitations by diagonalizing the Hessian matrix through a specific weighting scheme, which enables a single-step inversion for each frequency subproblem. This reformulation reduces the computational expense by approximately an order of magnitude relative to the Hessian-based approach. The weighting scheme also functions as a smoothing regularizer that is intrinsically embedded within the forward operator, thereby balancing computational efficiency and spatial resolution while producing robust reconstructions less sensitive to noise. Furthermore, for approximating the geometrical portion of the amplitude, this study introduces a paraxial ray-tracing system, further enhancing computational efficiency and accuracy. The inversion approach proposed in the present study has enabled the first successful translation of acoustic inverse scattering methods to a clinical setting.

Figures (9)

• Figure 1: Ray parameterization at the scattering point $\boldsymbol{x}'$. The wavevectors $\boldsymbol{k}(\boldsymbol{x}', e)$ (blue) and $\boldsymbol{k}(\boldsymbol{x}', r)$ (red) correspond to rays emanating from the emitter $e$ and receiver $r$, respectively, at $\boldsymbol{x}'$. These wavevectors are expressed in polar coordinates as $(k(\boldsymbol{x}'), \gamma(\boldsymbol{x}', e))$ and $(k(\boldsymbol{x}'), \gamma(\boldsymbol{x}', r))$. The vector $\mathbf{\bar{k}}(r,e, \boldsymbol{x}')$ passes through $\boldsymbol{x}'$ and represents the gradient of a two-way isochron, computed as $\mathbf{\bar{k}}(r,e, \boldsymbol{x}') = \boldsymbol{k}(\boldsymbol{x}', e) + \boldsymbol{k}(r, \boldsymbol{x}')$. Its polar representation is given by $(|\mathbf{\bar{k}}(r,e, \boldsymbol{x}')|, \zeta(r,e, \boldsymbol{x}') )$. The direction-reversed wavevector $\boldsymbol{k}(r, \boldsymbol{x}')$ is represented in polar coordinates as $(k(\boldsymbol{x}'), \gamma(\boldsymbol{x}', r) + \pi)$. Additionally, $\theta(r,e, \boldsymbol{x}')$ denotes the scattering angle. The dependence on $\omega$ is omitted for brevity.
• Figure 2: Phantom used for simulation of synthetic UST data using the full-wave approach: (a) sound speed $[ \mathrm{m} \mathrm{s}^{-1} ]$, (b) absorption coefficient $[ \mathrm{dB} \ \mathrm{MHz}^{-y} \ \mathrm{cm}^{-1} ]$, (c) smoothed sound speed $[ \mathrm{m} \mathrm{s}^{-1} ]$, (d) smoothed absorption coefficient $[ \mathrm{dB} \ \mathrm{MHz}^{-y} \ \mathrm{cm}^{-1} ]$. The maps are shown on a grid consisting of $502 \times 502$ points (used for the wave simulation). The maps in (c) and (d) are smoothed by applying an averaging window of size 17 points to the original acoustic maps in (a) and (b), respectively. The original (nonsmoothed) acoustic maps were used for simulating the data for image reconstruction. The smoothed maps were used for simulating the data for benchmarking against the ray approximation to the heterogeneous Green's function, as described in Section \ref{['sec:results-greens']} . The power law exponent was assumed to be $y = 1.4$ and homogeneous.
• Figure 3: Acoustic source used for all excitations (emitters): (a) time domain, (b) frequency domain: normalized amplitude and phase. $f_{\text{max}}$ indicates the maximum frequency supported by the grid used for the wave simulations.
• Figure 4: Pressure time series at all 256 receivers for a single frequency of 1 MHz following the excitation of emitters 1 and 20: (a) phase for emitter 1, (b) amplitude for emitter 1, (c) phase for emitter 20, and (d) amplitude for emitter 20. For the full-wave simulation, point sources in terms of $s$ were used.
• Figure 5: (a) Ground truth. Reconstructed sound speed images from synthetic data with 40 dB SNR: (b) Time-of-flight-based approach (initial guess). Assuming the true $\alpha_0$ (Figure \ref{['fig:1b']}): (c) Born, (d) Hessian-based ray-Born, (e) Hessian-free ray-Born. Assuming $\alpha_0 = 0$: (f) Born, (g) Hessian-based ray-Born, (h) Hessian-free ray-Born. Assuming $\alpha_0 = 0.5 \ \mathrm{dB \ MHz}^{-y} \mathrm{cm}^{-1}$ (homogeneous within the breast): (i) Born, (j) Hessian-based ray-Born, (k) Hessian-free ray-Born.

Theorems & Definitions (2)

• Definition 1
• Definition 2