We extend the forward-adjoint operator framework derived in our previous study to photoacoustic tomography (PAT). In that earlier work, the acoustic forward operator included a reception operator that maps, at each time step, the pressure wavefield in free space onto the boundary (receiver surface). It was shown that this reception operator serves as a left-inverse of an emission operator that maps the pressure restricted to the boundary (emitter surface) onto free space, perfectly complying with the reciprocity law of physics. In this study, we define the full PAT forward operator as a composite mapping composed of an acoustic forward operator equipped with a scaled variant of the previously proposed reception operator, and an operator describing the photoacoustic source. Singularities arising both in the reception step (due to the boundary restriction) and in the photoacoustic source (due to its instantaneous nature) are regularized using regularized Dirac delta distributions. The resulting PAT forward-adjoint operator pair satisfies an inner-product relation, which we verify through numerical experiments on a discretized domain. The effectiveness of the proposed operator pair is further demonstrated using an iterative minimization framework that yields both qualitatively and quantitatively accurate reconstructions of an initial pressure distribution from the corresponding Dirichlet-type boundary data.

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 weakly 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 operator, 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 transforming a full-frequency Hessian operator into an identity operator through a specific preconditioning scheme, thereby enabling 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, while yielding robust reconstructions with reduced sensitivity to noise by eliminating the ill-conditioning of the Hessian operator. Furthermore, compared with the author’s previous study, the present work improves the approximation of the geometrical component of the wave amplitude by introducing a paraxial ray-tracing framework, thereby further enhancing both computational efficiency and reconstruction accuracy. The inversion approach proposed in the present study has enabled the first successful translation of acoustic inverse scattering methods relying on first-scattered waves to a clinical setting.

The acoustic wave equation governs wave propagation induced by either volumetric radiation sources, or by surface sources of monopole or dipole type. For surface sources, boundary value problems yield wavefield representations via the Kirchhoff–Helmholtz or Rayleigh–Sommerfeld integrals. This study begins by examining the equivalence between the analytic expressions of the associated monopole and dipole integral formulations and their regularized approximations. Leveraging these regularized formulations, we introduce \emph{reception operators} that map free space pressure wavefields—obtained by solving the wave equation—onto measured fields restricted to the boundary. Building on this trace mapping, we derive the adjoint of the forward operator. We show that, under the common practical assumption of Dirichlet-type boundary data, the adjoint operator coincides—up to a constant factor—with the time-reversed form of the dipole integral formula, evaluated on the receiver surfaces. This study aims to advance the numerical approximation of forward problems and the solution of inverse problems in acoustics, with a particular focus on applications that require accurate amplitude modeling, including attenuation reconstruction and photoacoustic tomography.

This study presents the first experimental validation of a Hessian-free ray-Born inversion technique for quantitative reconstruction of sound speed from transmission ultrasound data. The method combines single-scattering theory with high-frequency approximations, yielding an inversion framework well suited to the frequency ranges used in clinical ultrasound applications. Unlike previous singly-scattered inversion approaches that account for medium heterogeneities only in the scattering potential, the proposed ray-Born method employs Green's functions approximated along ray trajectories determined by high-frequency assumptions. The associated objective function is linearized and minimized sequentially across increasing frequency bands. At each frequency set, the linearized subproblem is solved using a weighting scheme applied to both the solution and data spaces, which diagonalizes the Hessian and enables its inversion in a single step. The method, previously reported and released as an open-source package, was applied to in-vitro and in-vivo datasets provided by the University of Rochester Medical Center. The reconstructed images were evaluated by comparison with those obtained using a full-wave inversion approach based on a frequency-domain Helmholtz solver. The results demonstrate the strong potential of the Hessian-free ray-Born inversion as a computationally efficient and accurate method suitable for clinical translation.

We present a MATLAB package for reconstructing sound-speed images from transmission ultrasound data. The package is based on two-point ray tracing and implements two complementary inversion strategies for image reconstruction. The first is a time-of-flight (ToF) method that produces low-resolution, low-contrast images with minimal artefacts. The second is a ray-Born inversion method, which integrates high-frequency ray theory with the Born approximation to generate high-resolution sound-speed reconstructions. Early iterations of the ToF reconstruction are used to provide an initial estimate for the more advanced ray-Born approach. The core of this software package consists of four ray-tracing algorithms, whose accuracy is assessed in this study with respect to known analytical trajectories and accumulated acoustic path lengths. Furthermore, both image-reconstruction strategies have been validated numerically with simulated synthetic datasets and experimentally with open-source in-vitro and in-vivo datasets in related parallel studies.

The acoustic wave equation governs wave propagation induced by either volumetric radiation sources, or by surface sources of monopole or dipole type. For surface sources, boundary value problems yield wavefield representations via the Kirchhoff-Helmholtz or Rayleigh-Sommerfeld integrals. This study begins by establishing an equivalence between the analytic expressions of the associated monopole and dipole integral formulations and their full-waveform approximations. Leveraging this equivalence, we introduce reception operators that map free space pressure wavefields-obtained by solving the wave equation-onto measured fields restricted to the boundary. Building on this trace mapping, we derive the adjoint of the forward operator. We show that, under the common practical assumption of Dirichlet-type boundary data, the adjoint operator coincides-up to a constant factor-with the interior-field time-reversed form of the dipole integral formula, evaluated on the receiver surfaces. This study aims to advance the approximation of forward problems and the solution of inverse problems in acoustics, with a particular focus on applications that require accurate amplitude modeling, including therapeutic ultrasound optimization, attenuation reconstruction, and photoacoustic tomography.

Photoacoustic tomography is a contrast agent-free imaging technique capable of visualizing blood vessels and tumor-associated vascularization in breast tissue. While sophisticated breast imaging systems have been recently developed, there is yet much to be gained in imaging depth, image quality and tissue characterization capability before clinical translation is possible. In response, we have developed a hybrid photoacoustic and ultrasound-transmission tomographic system PAM3. The photoacoustic component has for the first time three-dimensional multi-wavelength imaging capability, and implements substantial technical advancements in critical hardware and software sub-systems. The ultrasound component enables for the first time, a three-dimensional sound speed map of the breast to be incorporated in photoacoustic reconstruction to correct for inhomogeneities, enabling accurate target recovery. The results demonstrate the deepest photoacoustic breast imaging to date namely 48 mm, with a more uniform field of view than hitherto, and an isotropic spatial resolution that rivals that of Magnetic Resonance Imaging. The in vivo performance achieved, and the diagnostic value of interrogating angiogenesis-driven optical contrast as well as tumor mass sound speed contrast, gives confidence in the system's clinical potential.

This study proposes a Hessian-inversion-free ray-born inversion approach for biomedical ultrasound tomography. The proposed approach is a more efficient version of the ray-born inversion approach proposed in [1]. Using these approaches, the propagation of acoustic waves are modelled using a ray approximation to heterogeneous Green’s function. The inverse problem is solved in the frequency domain by iteratively linearisation and minimisation of the objective function from low to high frequencies. In [1], the linear subproblem associated with each frequency interval is solved by an implicit and iterative inversion of the Hessian matrix (inner iterations). Instead, this study applies a preconditioning approach to each linear subproblem so that the Hessian matrix becomes diagonalised, and can thus be inverted in a single step. Using the proposed preconditioning approach, the computational cost of solving each linear subproblem of the proposed ray-Born inversion approach becomes almost the same as solving one linear subproblem associated with a radon-type time-of-flight-based approach using bent rays. More importantly, the smoothness assumptions made for diagonalising the Hessian matrix make the image reconstruction more stable than the inversion approach in [1] to noise.

Simulating propagation of acoustic waves via solving a system of three-coupled first-order linear differential equations using a k-space pseudo-spectral method is popular for biomedical applica-tions, firstly because of availability of an open-source toolbox for implementation of this numerical approach, and secondly because of its efficiency. The k-space pseudo-spectral method is efficient, because it allows coarser computational grids and larger time steps than finite difference and finite element methods for the same accuracy. The goal of this study is to compare this numerical wave solver with an analytical solution to the wave equation using the Green’s function for computing propagation of acoustic waves in homogeneous media. This comparison is done in the frequency domain. Using the k-Wave solver, a match to the Green’s function is obtained after modifying the approach taken for including mass source in the linearised equation of continuity (conservation of mass) in the associated system of wave equations.

An efficient and accurate image reconstruction algorithm for ultrasound tomography in soft tissue is described and demonstrated, which can recover accurate sound speed distribution from acoustic time series measurements. The approach is based on a second-order iterative minimisation of the difference between the measurements and a model based on a ray-approximation to the heterogeneous Green's function. It overcomes the computational burden of full-wave solvers while avoiding the drawbacks of time-of-flight methods. Through the use of a second-order iterative minimisation scheme, applied stepwise from low to high frequencies, the effects of scattering are incorporated into the inversion.

Over the past decade, the range of applications in biomedical ultrasound exploiting 3D printing has rapidly expanded. For wavefront shaping specifically, 3D printing has enabled a diverse range of new, low-cost approaches for controlling acoustic fields. These methods rely on accurate knowledge of the bulk acoustic properties of the materials; however, to date, robust knowledge of these parameters is lacking for many materials that are commonly used. In this work, the acoustic properties of eight 3D-printed photopolymer materials were characterised over a frequency range from 1 to 3.5 MHz. The properties measured were the frequency-dependent phase velocity and attenuation, group velocity, signal velocity, and mass density. The materials were fabricated using two separate techniques [PolyJet and stereolithograph (SLA)], and included Agilus30, FLXA9960, FLXA9995, Formlabs Clear, RGDA8625, RGDA8630, VeroClear, and VeroWhite. The range of measured density values across all eight materials was 1120-1180 kg  · m-3, while the sound speed values were between 2020 to 2630 m  · s-1, and attenuation values typically in the range 3-9 dB  ·  MHz-1· cm-1.

Ultrasound tomography (UST) has seen a revival of interest in the past decade, especially for breast imaging, due to improvements in both ultrasound and computing hardware. In particular, three-dimensional UST, a fully tomographic method in which the medium to be imaged is surrounded by ultrasound transducers, has become feasible. This has led to renewed attention on UST image reconstruction algorithms. In this paper, a comprehensive derivation and study of a robust framework for large-scale bent-ray UST in 3D for a hemispherical detector array is presented. Two ray-tracing approaches are derived and compared. More significantly, the problem of linking the rays between emitters and receivers, which is challenging in 3D due to the high number of degrees of freedom for the trajectory of rays, is analysed both as a minimisation and as a root-finding problem. The ray-linking problem is parameterised for a convex detection surface and two robust, accurate, and efficient derivative-free ray-linking algorithms are formulated and demonstrated and compared with a Jacobian-based benchmark approach. To stabilise these methods, novel adaptive-smoothing approaches are proposed that control the conditioning of the update matrices to ensure accurate linking. The nonlinear UST problem of estimating the sound speed was recast as a series of linearised subproblems, each solved using the above algorithms and within a steepest descent scheme. The whole imaging algorithm was demonstrated to be robust and accurate on realistic data simulated using a full-wave acoustic model and an anatomical breast phantom, and incorporating the errors due to time-of-flight (TOF) picking that would be present with measured data. This method can used to provide a low-artefact, quantitatively accurate, 3D sound speed maps. In addition to being useful in their own right, such 3D sound speed maps can be used to initialise full-wave inversion methods, or as an input to photoacoustic tomography reconstructions.

Quantitative photo-acoustic tomography (QPAT) seeks to reconstruct a distribution of optical attenuation coefficients inside a sample from a set of time series of pressure data that is measured outside the sample. The associated inverse problems involve two steps, namely acoustic and optical, which can be solved separately or as a direct composite problem. We adopt the latter approach for realistic acoustic media that possess heterogeneous and often not accurately known distributions for sound speed and ambient density, as well as an attenuation following a frequency power law that is evident in tissue media. We use a diffusion approximation (DA) model for the optical portion of the problem. We solve the corresponding composite inverse problem using three total variation (TV) regularised optimisation approaches. Accordingly, we develop two Krylov-subspace inexact-Newton algorithms that utilise the Jacobian matrix in a matrix-free manner in order to handle the computational cost. Additionally, we use a gradient-based algorithm that computes a search direction using the L-BFGS method, and applies a TV regularisation based on the alternating direction method of multipliers (ADMM) as a benchmark, because this method is popular for QPAT and direct QPAT. The results indicate the superiority of the developed inexact Newton algorithms over gradient-based quasi-Newton approaches for a comparable computational complexity.

We present an optimisation framework for photo-acoustic tomography of the brain based on a system of coupled equations that describe the propagation of sound waves in linear isotropic inhomogeneous and lossy elastic media with absorption and physical dispersion following a frequency power law using fractional Laplacian operators. The adjoint of the associated continuous forward operator is derived, and a numerical framework for computing this adjoint based on a k-space pseudo-spectral method is presented. We analytically show that the derived continuous adjoint matches the adjoint of an associated discretised forward operator. We include this adjoint in a first-order positivity constrained optimisation algorithm that is regularised by total variation minimisation, and show that the iterates monotonically converge to a minimiser of an objective function, even in the presence of some error in estimating the physical parameters of the medium.

A class of sparse optimization techniques that require solely matrix–vector products, rather than an explicit access to the forward matrix and its transpose, has been paid much attention in the recent decade for dealing with large-scale inverse problems. This study tailors application of the so-called Gradient Projection for Sparse Reconstruction (GPSR) to large-scale time-difference three-dimensional electrical impedance tomography (3D EIT). 3D EIT typically suffers from the need for a large number of voxels to cover the whole domain, so its application to real-time imaging, for example monitoring of lung function, remains scarce since the large number of degrees of freedom of the problem extremely increases storage space and reconstruction time. This study shows the great potential of the GPSR for large-size time-difference 3D EIT. Further studies are needed to improve its accuracy for imaging small-size anomalies.

Existing total variation (TV) solvers that have been applied in electrical impedance tomography (EIT) smooth the TV function in order to cope with its non-differentiability around the origin, and thus imposes some numerical errors on the solution. Furthermore, these solvers require storage of Hessian, and are thus very impractical for large-scale computations, especially 3D EIT. These shortcomings were addressed by TV solvers that are based on first-order optimization methods. However, the application of these solvers to EIT remains scarce. In this manuscript, we propose an accelerated version of a gradient-based TV solver based on augmented Lagrangian and alternating direction method of multipliers, referred to as TVAL3, and apply it to EIT. The results demonstrate the superiority of the accelerated algorithm over existing TV solvers in EIT with regard to both accuracy and speed. (C) 2016 Elsevier Inc. All rights reserved.

Inspired by the recent advances on minimizing nonsmooth or bound-constrained convex functions on models using varying degrees of fidelity, we propose a line search multigrid (MG) method for full-wave iterative image reconstruction in photoacoustic tomography (PAT) in heterogeneous media. To compute the search direction at each iteration, we decide between the gradient at the target level, or alternatively an approximate error correction at a coarser level, relying on some predefined criteria. To incorporate absorption and dispersion, we derive the analytical adjoint directly from the first-order acoustic wave system. The effectiveness of the proposed method is tested on a total-variation penalized Iterative Shrinkage Thresholding algorithm (ISTA) and its accelerated variant (FISTA), which have been used in many studies of image reconstruction in PAT. The results show the great potential of the proposed method in improving speed of iterative image reconstruction.

Inspired by the recent advances on minimizing nonsmooth or bound-constrained convex functions on models using varying degrees of fidelity, we propose a line search multi-grid (MG) method for full-wave iterative image reconstruction in photoacoustic tomography (PAT) in heterogeneous media. To compute the search direction at each iteration, we decide between the gradient at the target level, or alternatively an approximate error correction at a coarser level, relying on some predefined criteria. To incorporate absorption and dispersion, we derive the analytical adjoint directly from the first-order acoustic wave system. The effectiveness of the proposed method is tested on a total-variation penalized Iterative Shrinkage Thresholding algorithm (ISTA) and its accelerated variant (FISTA), which have been used in many studies of image reconstruction in PAT. The results show the great potential of the proposed method in improving speed of iterative image reconstruction.