Heterogeneous Angular Spectrum Approach for Trans-skull Imaging and Focusing

TL;DR

This paper introduces a heterogeneous angular spectrum approach (ASA) for trans-skull ultrasound imaging, providing accurate, computationally efficient aberration correction for brain disease diagnosis and treatment.

Contribution

The authors develop a general numerical ASA method for heterogeneous media, enabling fast, accurate trans-skull focusing and localization with reduced errors and increased intensity.

Findings

Achieves 75% error reduction in source localization

Increases peak intensity by 40-60%

Operates with millisecond computation times

Abstract

Ultrasound, alone or in concert with circulating microbubble contrast agents, has emerged as a promising modality for therapy and imaging of brain diseases. While this has become possible due to advancements in aberration correction methods, a range of applications, including adaptive focusing and tracking of the microbubble dynamics through the human skull, may benefit from even more computationally efficient methods to account for skull aberrations. Here, we derive a general method for the angular spectrum approach (ASA) in a heterogeneous medium, based on a numerical marching scheme to approximate the full implicit solution. We then demonstrate its functionality with simulations for (human) skull-related aberration correction and trans-skull passive acoustic mapping. Our simulations show that the general solution provides accurate trans-skull focusing as compared to the uncorrected case for clinically relevant frequencies, apertures, and targets, with the effects of skull attenuation and amplitude shading included. In the case of source localization, our method leads to an average of 75 % error reduction and 40 to 60 % increase in peak intensity, evaluated over the range of frequencies as compared to the homogeneous medium ASA. Overall, total computation times for both focusing and point source localization of the order milliseconds can be attained with this approach. Collectively our findings indicate that the proposed phase correction method based on the ASA could provide a computationally efficient and accurate method for trans-skull transmit focusing and imaging of point scatterers, potentially opening new possibilities for treatment and diagnosis of brain diseases.

Tables

Table 1. Table 1: Mean and standard deviation focal location error and spot size with uncorrected [ Eq.   8 ] and the corrected [ Eq.   9 ] focusing delays, averaged over all frequencies and positions shown in Figs.   3 and  4 .

	Uncorrected		Corrected
Aperture	Error [mm]	Size [mm2]	Error [mm]	Size [mm2]
50 mm	2.1 $\pm$ 1.2	6.8 $\pm$ 5.7	1.3 $\pm$ 1.0	12.6 $\pm$ 16.9
100 mm	2.2 $\pm$ 0.7	7.5 $\pm$ 6.4	1.0 $\pm$ 0.4	9.6 $\pm$ 10.9

Table 2. Table 2: Mean and standard deviation PAM localizaion error and peak intensities with uncorrected [ Eq.   3 ] and the corrected [ Eq.   6 ] beamforming, computed over the positions shown in Fig.   7 .

	Uncorrected		Corrected
Aperture	Error [mm]	Intensity	Error [mm]	Intensity
50 mm	3.7 $\pm$ 2.2	0.73 $\pm$ 0.46	1.2 $\pm$ 0.7	1.0 $\pm$ 0.60
75 mm	2.5 $\pm$ 1.7	0.66 $\pm$ 0.36	0.9 $\pm$ 0.5	1.0 $\pm$ 0.47
100 mm	3.5 $\pm$ 1.9	0.62 $\pm$ 0.28	0.8 $\pm$ 0.4	1.0 $\pm$ 0.39

Equations

$$P(k_x,k_y,z)$$

$$\equiv {\iint }_{-\infty }^{\infty }\tilde{p}(x,y,z)e^{-i\left(k_{x}x+k_{y}y\right)}\mathrm{d}x\mathrm{d}y.$$

$$\frac{{\mathrm{d}}^{2}P}{\mathrm{d}{z}^2}+k_z^{2}P=0,$$

$$P$$

$$\frac{{\mathrm{d}}^{2}P}{\mathrm{d}{z}^2}+k_z^{2}P$$

$$P$$

$$P^{n+1}\approx P^n{e}^{i{k}_z\Delta z}+\frac{e^{i{k}_z\Delta z}}{2i{k}_z}\left(P^n\ast \Lambda \right)\times \Delta z,$$

$$M$$

$$\tau (x)=\frac{\arg \tilde{p}(x,d)}{\omega }.$$

$$\tau (x)=\frac{\sqrt{(x-x_0{)}^2+d^2}}{c_0}.$$

$${\nabla }^{2}p-\frac{1}{c^2(\mathbf{r})}\frac{{\partial }^{2}p}{\partial t^2}$$

$$c(\mathbf{r})$$

$$\left({\nabla }^2+k_0^2\right)\tilde{p}$$

$${\mathcal{F}}_k\left[{\nabla }^2\tilde{p}\right]=\left(-k_x^2-k_y^2+\frac{{\partial }^2}{\partial z^2}\right)P,$$

$$\left({\nabla }^2+k_0^2\right)\tilde{p}$$

$$={\mathcal{F}}_k^{-1}\left[\left(-k_x^2-k_y^2+\frac{{\partial }^2}{\partial z^2}\right)P+k_0^{2}P\right]$$

$$\displaystyle=\mathcal{F}_{k}^{-1}\,\Bigg{[}\Bigg{(}\underbrace{k_{0}^{2}-k_{x}^{2}-k_{y}^{2}}_{=\,k_{z}^{2}}+\frac{\partial^{2}}{\partial z^{2}}\Bigg{)}P\Bigg{]}\,.$$

$${\mathcal{F}}_k^{-1}\left[\frac{{\partial }^{2}P}{\partial z^2}+k_z^{2}P\right]$$

$$⟹\frac{{\partial }^{2}P}{\partial z^2}+k_z^{2}P=\Lambda \ast P.$$

$$P_h=A{e}^{i{k}_{z}z}+B{e}^{-i{k}_{z}z},$$

$$g(z|z^{\prime })$$

$$P$$

$$P^n=P(k_x,k_y,z_n),$$

$$P^{n+1}=P(k_x,k_y,z_n+\Delta z).$$

$$P^{n+1}\approx P^n{e}^{i{k}_z\Delta z}+\frac{e^{i{k}_z\Delta z}}{2i{k}_z}\left(P^n\ast \Lambda \right)\times \Delta z.$$

Full text

**© 2020 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. Article DOI: 1010.1109/TMI.2019.2953872. **

1 Introduction

Ultrasound has emerged as a novel modality for the treatment and imaging of brain diseases and disorders [1, 2]. When enhanced by circulating microbubble contrast agents—lipid, albumin, or polymer-shelled gas pockets ($110\text{\,}\mathrm{\SIUnitSymbolMicro m}$) that scatter sound and vibrate in response to incident ultrasound—it can enable a range of new therapeutic interventions[3, 4, 5, 6, 7, 8] and open new possibilities for imaging [9, 10, 11]. This has been facilitated by advancements in methods to account for skull aberrations.

In their simplest form, aberrations introduced by the skull may be corrected with knowledge of phase and amplitude at the transmitter due to a point source at the target, that can then be conjugated to achieve trans-skull focusing on transmit. Such phase and amplitude modulations may be determined experimentally via measurement of the field due to an induced cavitation source at the target location inside the skull [12, 13, 14], however cavitation events in the brain are not always desirable, and the relatively high cavitation threshold requires additional corrections or high power arrays [15]. Use of corrections based on analysis of time domain signals measured from contrast agent microbubble emissions has been demonstrated[16], though the location of these cavitation events cannot be controlled precisely. Alternatively, a point source can be placed within the skull to enable effective trans-skull therapy [17, 18, 19, 20]; however this technique requires invasive placement at the target, which is not possible in realistic scenarios. For these reasons, several simulation techniques have been proposed for aberration correction, including finite difference time domain (FDTD) [21, 22] and k-space propagation models[23, 24, 25]. In these methods, the simulation of the acoustic propagation due to a point source inside the brain is used to estimate the required phases for aberration correction. While these methods can be very precise, as they may account for a broad range of physical effects including, e.g., mode conversion, viscosity and nonlinearity, there exists a fundamental trade-off between their accuracy and required computational complexity [26].

Currently an effective trade-off is obtained by the use of a stepwise, locally homogeneous angular spectrum approach (ASA), which employs either coordinate rotations to account for refraction[23] or accounts for spatial heterogeneity in the space domain before spectral propagation[27]. While these ASA implementations currently provide a reasonable compromise between computational complexity and aberration correction, they do not solve intrinsically the full wave propagation problem [27, 28, 29, 26].

Inherently related to the problem of errors in trans-skull FUS targeting due to skull aberrations is that of source localization. This problem has emerged as a priority in microbubble enhanced US therapy and imaging: when exposed to ultrasound, microbubbles scatter diverging pressure waves that, due to their size (several orders of magnitude smaller than the wavelength), act as point sources. The information carried by these waves can enable the spatiotemporal characterization of the microbubble dynamics. During therapeutic interventions, this information is used to ensure that the desired type of oscillation (i.e., stable or inertial) is taking place at the intended location (i.e. treatment monitoring) [30].

Several approaches based on passive acoustic mapping (PAM) techniques have been proposed for monitoring the microbubble dynamics. While direct implementations do not account for aberrations [31, 32, 33, 30], recent work has successfully incorporated human skull aberration corrections into both frequency-domain[34] and time-domain [11, 35, 26, 36, 37, 38] PAM implementations. These methods require either measurement (invasive) [17, 14], simulations to calculate delays for each point in the image, or a representative point to be used for a local region [39, 26], which incurs a computational expense that scales with the field of view. These methods are thus either infeasible or time consuming (up to several minutes) unless GPUs are used to speed up the computations [26]. This problem can potentially be overcome with the ASA method[30], however current implementations do not account for aberrations.

The spectral selectivity inherent to frequency domain methods such as the ASA [33, 30] is important for characterizing the type of oscillation (e.g. harmonics vs broadband). However current implementations assume a homogeneous medium. In addition to spectral selectivity, reconstructions on the order of milliseconds are important for closed-loop control of the microbubble dynamics[40] and improved temporal resolution during microbubble imaging.[35, 10] While time domain methods [31, 32, 38] can readily incorporate corrections and extract for certain frequency content via filtering, they generally incur high computational loads, unless GPU units are used to speed up the computations [38]. Hence, fast and frequency selective reconstructions to visualize the cerebrovascular microbubble dynamics through the skull may have important implications for image guided therapy and imaging [8].

Herein we derive the general solution for the heterogeneous ASA. First, we present the derivation of a fast phase correction method for arbitrarily distributed, weakly heterogeneous medium. Then, through simulated acoustic propagation, we show the numerical implementation of the algorithm for focal aberration correction and frequency selective passive mapping of point sources through a human skull. Next, we evaluate the focal aberration error, focal pressure, and spot size, as well as the point source localization error and intensity, compared with the corresponding uncorrected cases. Finally, we evaluate the computational cost of the different algorithm permutations.

2 Theory

2.1 Angular Spectrum Approach

The angular spectrum $P$ of a monochromatic111A time convention of $p=\tilde{p}\,e^{-i\omega t}$ is used, so that the forward temporal transform uses a kernel of $e^{+i\omega t}$, and that of the forward spatial transform is then $e^{-i(k_{x}x+k_{y}y)}$. The factors of $1/2\pi$ are applied on the inverse transforms. pressure field $\tilde{p}$ with angular frequency $\omega$ is given by its 2D spatial Fourier transform

[TABLE]

Applying the spatial transform to the homogeneous Helmholtz equation $(\nabla^{2}+k^{2})\tilde{p}=0$ yields an ordinary differential equation for the angular spectrum

[TABLE]

where $k_{z}^{2}=(\omega/c_{0})^{2}-k_{x}^{2}-k_{y}^{2}$, and $c_{0}$ is the (constant) small signal sound speed. If $P_{0}$ is known at some reference plane $z=0$, and if there are no backward-travelling waves, then Eq. 2 has the solution

[TABLE]

The acoustic field in any plane may then be reconstructed with Eq. 3 and evaluation of the inverse transform.

2.2 Extension to Heterogeneous Media

In the case of relatively weak heterogeneity—i.e. the sound speed $c(\bm{r})$ changes slowly compared with the wavelength—the governing equation is

[TABLE]

Here, $\Lambda=\mathcal{F}_{k}\left[k_{0}^{2}(1-\mu)\right]$, $k_{0}=\omega/c_{0}$, $\mu=c_{0}^{2}/c^{2}$, $c_{0}$ is a reference (average) sound speed, and $*$ denotes a 2D convolution over the wavenumber components $k_{x}$ and $k_{y}$. From Eq. 2, the heterogeneity appears as a source term in the governing equation. In the general case, the implicit solution of Eq. 4 may be obtained with a Green’s function technique

[TABLE]

A numerical approximation of Eq. 5 may be made to compute $P$ at arbitrary $z$ via

[TABLE]

where $P^{n}=P(k_{x},k_{y},n\Delta z)$. See Appendix A for derivations of Eqs. 4, 5 and 6. Previously, Gu and Jing persented a general forward propagation scheme for field simulations that includes nonlinearity and attenuation [41]. In the absence of these effects, Eq. (12) of that reference becomes (with the notation of this paper)

[TABLE]

Substitution of this value of $M$ into their Eq. (11) recovers Eq. 5 obtained here, indicating the consistency of the results.

3 Methods

3.1 Simulations

Acoustic propagation (either to evaluate the transmit focusing achieved with the computed time delays or to acquire the RF data to be used with the corrected beamforming) was simulated in k-Wave [42]. As the incident angles were largely below the critical angle for skull bone ($\sim$26\text{,}\mathrm{\SIUnitSymbolDegree}$$[43]), elastic effects were not included in the simulation. The skull and tissue densities were computed from CT data of a human skull, and the spatially-dependent sound speed and absorption coefficient $\alpha$ were computed from the Hounsfield scale conversion as described previously [21, 37]. Absorption was included via a power law model $\alpha(\bm{r})=\alpha_{0}(\bm{r})\cdot f^{\beta}$ [44]. A power law exponent of $\beta=1.2$ was defined for the entire medium [45], as k-Wave’s implementation does not allow spatial variation of this parameter. Spatial grid spacing $200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and time step $40\text{\,}\mathrm{ns}$ were used for all simulations (CFL number $c\Delta t/\Delta x=0.44$), and 2000 time steps ($80\text{\,}\mathrm{\SIUnitSymbolMicro s}$) were simulated. The reference sound speed was taken to be the average value over the entire grid. All reconstruction and aberration correction scripts were written in MATLAB and run on a standard desktop computer (Intel Core i7, four cores at 2.8 GHz and 16 GB memory) without parallel or graphical processing techniques. The general computational flows for focal aberration correction and source localization are shown in Fig. 1 and Fig. 2, respectively. While the correction can function in 3 dimensions [see, e.g., Fig. 2(c)], due to the large number of simulations required, simulations were performed in 2D ($60\text{\,}\mathrm{m}\mathrm{m}$ by $80\text{\,}\mathrm{m}\mathrm{m}$, with 10 point perfectly matched layer) for computational efficiency. Thus in the reconstructions $k_{y}=0$.

3.2 Numerical Implementation

3.2.1 Transmit Focusing

To determine the improvement in focal aberration correction through human skull with the proposed general solution (heterogeneous ASA), the focusing delays were computed with and without phase corrections for clinically relevant frequencies ($0.251.5\text{\,}\mathrm{MHz}$) [2, 5, 46, 47]. To achieve focusing at a desired depth $d$, $P(k_{x},d)$ was computed for a delta function at the desired focus, i.e., for $P_{0}=\mathcal{F}_{k}\left[\delta(x-x_{0},0)\right]$. To avoid very steep changes on the spatial computation grid for the corrected case, this delta function was approximated as a cosine-windowed Gaussian distribution with full width at half maximum of $1.5\text{\,}\mathrm{m}\mathrm{m}$. Since the field is monochromatic, the excess phase represents the total time delay $\tau$ for an array element as a function of its position along the transducer face

[TABLE]

Phase unwrapping was used to preserve the full phase of $P$ at the transducer location. For reference, the results of focusing using the delays calculated with Eq. 8 were compared with the results obtained without accounting for medium heterogeneity (i.e., geometrical focusing). In this case, the delays were computed with $c_{0}=$1500\text{,}\mathrm{m}\mathrm{/}\mathrm{s}$$ directly from the target point via

[TABLE]

The excitation time series for each element was a 40-cycle sine pulse with $100\text{\,}\mathrm{kPa}$ amplitude at the desired frequency, modulated with a Tukey window with $R=0.1$, and shifted by the amount given by Eq. 8 [or Eq. 9 in the uncorrected case]. The relative amplitude for each element was scaled according to the normalized amplitude of $\tilde{p}(x,d)$ (or unity in the geometric case); the maximum amplitude of the corrected case was set to be $20\text{\,}\mathrm{\char 37\relax}$ larger than the uncorrected case due to the tapering of the outer elements. The field was then simulated with each element of the array transmitting with the calculated phase delay and amplitude.

3.2.2 Passive Acoustic Mapping

Next, the effect of focal aberration correction for the passive acoustic mapping of microbubbles was assessed. Bubbles were approximated as Gaussian sources with center frequencies of 0.4, 0.8, and 1.2 MHz (to represent harmonic components of bubbles excited with trans-skull FUS at 200 kHz or 600 kHz) [30, 2, 5]. Sources ($N=150$) were placed randomly within a region (approximately 2 mm wide by 3 cm long) representing a vessel shape, within a sound speed environment defined from the CT data (see Fig. 5). To determine the error as a function of position and array aperture, $1\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$ sources were then individually simulated in a rectangular grid at axial positions of 30 to 80 mm, and transverse locations from $-$20 to 20 mm.

The resulting pressure from these sources was measured by a virtual linear array with $200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ pitch (e.g., 250 elements for 50 mm array). The RF data $\tilde{p}(x,0)$ were then transformed to give $P_{0}$, from which $P(k_{x},z)$ was computed with Eq. 3 or Eq. 6 as indicated, with a step size $\Delta z=$50\text{,}\mathrm{\SIUnitSymbolMicro m}$$. The pressure distribution was then found from $\tilde{p}=\mathcal{F}_{k}^{-1}[P(k_{x},z)]$, and the pixel intensity field (i.e., the PAM) was computed as $I(x,z)=\|\tilde{p}(x,z)\|^{2}$. The convolutions in Eq. 6 were computed in the frequency domain to improve efficiency, since $P*\Lambda=\mathcal{F}_{k}^{-1}[\tilde{p}\cdot\lambda]$. Finally, the recovered source positions were taken to be the maxima of the generated PAMs and compared with the known source positions from the simulation. An algorithm for the reconstruction process is given in the Supplementary Material.

Because $P$ is proportional to $\operatorname{exp}{ik_{z}z}$, back propagation incurs multiplication by $\operatorname{exp}{-ik_{z}z}$. Evanescent components of the angular spectrum, for which $k_{z}$ is pure imaginary, will then grow exponentially. Therefore, all measured angular spectra $P_{0}$ were windowed with a Tukey window with cosine fraction $R=0.25$ to taper these components [48]. Additionally, all initial spectra were zero padded such that their computational extent was four times larger than their physical extent. Sound speed fields $c(\bm{r})$ were padded with their edge values replicated to match the grid size of the padded $P_{0}$. The axial step $\Delta z$ was $50\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and the field was computed from $z=0$ to $90\text{\,}\mathrm{mm}$.

4 Results

4.1 Focal Aberration Correction

First we estimated the effect of the phase correction on transmit focal accuracy for various focal targets as a function of frequency (Fig. 3). For each target focus, the focal delays were computed geometrically with Eq. 9 and then with corrections given by Eqs. 6 and 8. Across all focal positions and frequencies, the focal error was reduced from $2.1\pm 1.2\text{\,}\mathrm{mm}$ without phase corrections to $1.3\pm 1.0\text{\,}\mathrm{mm}$ with the correction. Aberration errors in focal targeting were generally larger at off-axis focal positions. See, for example the highlighted case in Fig. 3(c–e), where the uncorrected focal error was $5.1\text{\,}\mathrm{mm}$, while the error for the corrected case was $0.6\text{\,}\mathrm{mm}$.

We also evaluated the impact of the standoff distance of the array form the skull on focusing. Results indicate that the focal accuracy and improvement were comparable for different standoff distances ($0.6\pm 0.3\text{\,}\mathrm{mm}$ with the correction and $1.6\pm 1.3\text{\,}\mathrm{mm}$ without, see Supplementary Fig. S-1).

Next we assessed the normalized (to free field) focal pressure and spot size (i.e., the total area within $3\text{\,}\mathrm{d}\mathrm{B}$ of the peak pressure) achieved with the proposed aberration correction algorithm as a function frequency and position. Figure 4 illustrates how the spot size and focal pressure vary in the corrected case compared with the water path case as a function of frequency and position. The decrease in maximum pressure with frequency that is observed [Fig. 4(b–c)] reflects absorption losses that increase as $f^{\beta}$ ($\beta=1.2$ was used). Relative to the respective water case, the area of the focal region in the corrected trans-skull case did increase higher frequencies; however the absolute area became smaller, from $22.9\pm 13.8\text{\,}\mathrm{m}\mathrm{m}\mathrm{{}^{2}}$ for $250\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ to $6.5\pm 3.7\text{\,}\mathrm{m}\mathrm{m}\mathrm{{}^{2}}$ for the $1.5\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$ case. Both the spot size and focal pressure were seen to vary more with the focus’ axial position than with its transverse position [compare Fig. 4(b&d) and Fig. 4(c&e)]. The focal error and spot size in the corrected and uncorrected cases averaged over all frequencies and positions are summarized in Table 1. Together, these data demonstrate that the proposed method is able to correct aberrations and improve focal targeting through the human skull.

4.2 Passive Acoustic Mapping

To assess the ability of the proposed method for PAM we estimated the errors in the peak location compared with the source position for three different frequencies (Figure 5). The total (radial) error using the heterogeneous ASA was reduced compared to that for the PAMs based on the homogeneous medium (i.e., uncorrected) ASA. Specifically, using Eq. 6, the localization errors were $1.6\pm 0.9\text{\,}\mathrm{mm}$, $0.6\pm 0.4\text{\,}\mathrm{mm}$, and $0.6\pm 0.5\text{\,}\mathrm{mm}$ for $400\text{\,}\mathrm{kHz}$, $800\text{\,}\mathrm{kHz}$, and $1.2\text{\,}\mathrm{MHz}$ sources, respectively, compared with $2.8\pm 1.5\text{\,}\mathrm{mm}$, $2.0\pm 1.1\text{\,}\mathrm{mm}$, and $2.2\pm 1.8\text{\,}\mathrm{mm}$ when Eq. 3 was used (i.e., when the presence of the skull was not accounted for). Importantly, across all simulated point sources, there were no outliers in the corrected case [e.g., in Figure 5(d), where the PAMs for the indicated sources are shown at right]. Further, the intensities of the peaks in the reconstructions were larger and more uniform in the corrected case overall: normalized to the mean corrected peak intensity for each case, the uncorrected PAMs had peak intensities of 0.94$\,\pm\,$0.09, 0.78$\,\pm\,$0.12, and 0.59$\,\pm\,$0.14 for the $400\text{\,}\mathrm{kHz}$, $800\text{\,}\mathrm{kHz}$, and $1.2\text{\,}\mathrm{MHz}$, cases respectively.

Next we assessed the localization accuracy as a function of the aperture size. For larger apertures with the same pitch (i.e., with additional elements), the localization is further improved for the corrected case. Figure 6 shows the effect of the aperture size on localization accuracy for the trans-skull case with the uncorrected ASA (top row) and corrected ASA (bottom row) for sources at 400 kHz. For the corrected case, the localization error decreases with increasing aperture size ($1.5\pm 0.9\text{\,}\mathrm{mm}$, $0.4\pm 0.3\text{\,}\mathrm{mm}$, and $0.5\pm 0.2\text{\,}\mathrm{mm}$ for $50\text{\,}\mathrm{mm}$, $75\text{\,}\mathrm{mm}$, and $100\text{\,}\mathrm{mm}$ apertures, respectively). For the uncorrected case, a larger aperture does not reduce the error [$2.7\pm 1.6\text{\,}\mathrm{mm}$, $1.2\pm 0.8\text{\,}\mathrm{mm}$, and $1.5\pm 0.6\text{\,}\mathrm{mm}$ for $50\text{\,}\mathrm{mm}$, $75\text{\,}\mathrm{mm}$, and $100\text{\,}\mathrm{mm}$ apertures, respectively; see Fig. 6(b–c)].

Finally, we estimated the localization error as a function of position and aperture using a rectangular grid of $1\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$ sources within the skull (Figure 7) Across all positions, the localization error was reduced by $6080\text{\,}\mathrm{\char 37\relax}$, and the localization error in the corrected positions decreased monotonically with aperture for the phase corrected case. Further, the peak intensity in the uncorrected case was on average $3040\text{\,}\mathrm{\char 37\relax}$ smaller than in the corrected case. Table 2 summarizes the mean and standard deviation of the error over all points in the grid in Fig. 7. Full details, including the axial and transverse error and spot size are included in the Supplementary Tables S-1–S-3.

4.3 Computational Efficiency

For focal aberration correction, calculation of the focusing delays and amplitudes from Eqs. 6 and 8 required $25.2\pm 8.5\text{\,}\mathrm{ms}$. Note that registration of the the 2D sound speed map with the computational grid required $1.63\pm 0.03\text{\,}\mathrm{s}$, but only needs to be done once during registration. For PAM, the correction using the heterogeneous ASA ([Eq. 6]) the reconstructions took approximately $166\pm 37\text{\,}\mathrm{ms}$, for a step size $\Delta z=$50\text{,}\mathrm{\SIUnitSymbolMicro m}$$ over all calculations. This is approximately four times slower than the uncorrected ASA [i.e., with Eq. 3], where the computation of the PAMs required approximately $44\pm 4\text{\,}\mathrm{ms}$ per frequency. For a 3D PAM example case, $540\pm 18\text{\,}\mathrm{ms}$ was required per frequency to evaluate Eq. 6 for data from a 24-by-24 element array and $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ step size (173 376 voxels) and 3054 time samples. All reported times are for a standard desktop computer with no parallel processing or GPU acceleration.

5 Discussion

This work augments the ASA for sound propagation in heterogeneous media with a fast phase correction technique (Figs. 1 and 2). The methods were derived under assumptions of relatively weak heterogeneity, specifically that $\left|\nabla\rho_{0}/\rho_{0}+2\nabla c/c\right|$ is small compared with the wavelength (see Supplementary Material). While these conditions are not strictly met by the impedance contrast represented by the skull, our results from simulations with clinical CT data have demonstrated more than $50\text{\,}\mathrm{\char 37\relax}$ reduction in localization and transmission focusing error at clinically relevant frequencies and geometries. This improvement in accuracy is on the scale of a half wavelength with little additional computational burden.

Via phase extraction and amplitude shading, the error in focal targeting was reduced by approximately $60\text{\,}\mathrm{\char 37\relax}$ as compared to the uncorrected case (from $2.2\pm 0.7\text{\,}\mathrm{mm}$ without phase corrections to $1.0\pm 0.4\text{\,}\mathrm{mm}$ with the full correction—over a range of positions and frequencies, Fig. 3). The improvement in accuracy is more pronounced at higher frequencies, since the derivation assumes $c(\bm{r})$ varies slowly compared with the wavelength; thus the approximation is more valid for smaller wavelengths (higher frequencies). This is not considered as a major impediment, as aberration is more pronounced at higher frequencies, where the performance of the proposed method is favorable (Fig. 4).

Further, the method was seen to reduce the error in trans-skull source localization from computed PAMs by $6080\text{\,}\mathrm{\char 37\relax}$, (from $2.5\pm 2.3\text{\,}\mathrm{mm}$ to $0.6\pm 0.5\text{\,}\mathrm{mm}$) for a $1.2\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$ sources (Fig. 5). Localization accuracy was seen to increase with larger apertures (Fig. 6). A larger aperture imparts higher spatial frequency resolution (i.e., smaller $\Delta k_{x}$ and $\Delta k_{y}$ in the discrete transform), and thus improved coherence of the reconstruction at the correct source location. Larger apertures in the uncorrected case do not in general improve the accuracy of the localization in the presence of aberrations.

Compared with previously proposed methods for adaptation of ASA to heterogeneous media, our method has the advantage of requiring a single propagation step, and does not require transformations between each axial position [27, 23, 49]. Its inherent efficiency should allow its use in conjunction with real-time control methods based on PAM that allow for location-specific cavitation control [40]. Moreover, the proposed algorithm is faster than full simulation methods for phase correction (e.g., 166 ms, compared with 2 minutes for the corresponding 2D k-Wave simulation). Three dimensional heterogeneous PAM examples required 520 ms per frequency for $1.7\times 10^{5}$ voxels with 576 channels and 3054 time steps, which is comparable to GPU-optimized homogeneous time domain algorithms run on high performance computers (e.g., 85 ms for $10^{3}$ voxels and $3\times 10^{5}$ time samples in Ref. 38).

Computationally efficient methods may also lead to improved aberration correction accuracy by allowing iterative focal improvement. Such approaches are further supported by clinical experience that has shown that one still needs to correct for small errors ($12\text{\,}\mathrm{m}\mathrm{m}$) during targeting. [50, 2] Currently, these corrections are calculated from CT data and confirmed via MRI–guided methods, wherein low-level focal heating (a few degrees Celsius) is measured with MR temperature imaging (MRTI) [51]. In experimental settings, aberration correction has also been confirmed by small tissue displacements (a few microns) that are induced by radiation force, using MR acoustic radiation force imaging (MR-ARFI) [52, 53]. Based on these methods, optimization routines can be developed to refine the phase of each element by maximizing focal displacement [54, 28] or focal heating for a given excitation pulse. To ensure that such adaptive focusing methods are performed in timescales similar to current CT-based methods for a single correction (of the order of seconds) the computation of the phase delays needs to take place on time scales of the order of a few milliseconds. Hence the proposed methods can also provide fast PAM for spatio-temporal control of the cerebrovascular microbubble dynamics.

Our study has some limitations. First, only a single skull has been used on our simulations; in future work we aim to evaluate the proposed approach using different skulls and skull segments. Additionally, experimental validation using different skulls, microbubbles, and phased arrays is needed for assessing the robustness of the proposed method for focal aberration correction and imaging (passive or active). Further, the algorithm is derived under the assumption of forward propagation only; extension of the method to account for reflection is nontrivial, but could improve the localization and focal targeting accuracy. The reported computation times have strong dependence on the hardware, size of the domain, and the choice of reconstruction parameters (e.g., step size, number of pad channels, etc.). However, in all cases presented it could be reduced to a few milliseconds with more sophisticated implementations, including more powerful hardware, parallel computing, and graphical processing techniques. Finally, our formulation does not include reflection, attenuation, and nonlinearity; however, the proposed method can be readily augmented to account for the effects of attenuation,[41] and nonlinearity[55, 56] potentially enabling the quantification of the acoustic emissions and the investigation of nonlinear sound propagation though heterogeneous media, albeit at the price of some additional computational burden.

6 Conclusion

We derived the general solution for the heterogeneous ASA. Numerical data showed that the general solution provides accurate trans-skull transmit focusing and point source localization. Sub-millimeter errors were attained with only a modest increase in computational complexity both for focusing and point source localization using clinically relevant frequencies ($0.251.5\text{\,}\mathrm{MHz}$) and array apertures ($50100\text{\,}\mathrm{mm}$). The computation times on a standard computer were comparable to those reported for GPU-optimized uncorrected time domain algorithms on high performance computers. Collectively, our findings indicate that the heterogeneous ASA may create new possibilities for treatment, treatment monitoring, and diagnosis of brain diseases.

Appendix A Derivation of ASA Results

A.1 Governing Equation

Provided the medium is weakly heterogeneous (see Supplementary Material), propagation may be described by

[TABLE]

The sound speed may be written as the sum of a reference sound speed $c_{0}$ and a spatially-varying part $c^{\prime}(\bm{r})$[57]

[TABLE]

Because Eq. A.1 is most valid for small changes in sound speed [58], the mean value of $c(\bm{r})$ is the natural choice for the reference sound speed $c_{0}$. Defining $\mu(\bm{r})=c_{0}^{2}/c^{2}(\bm{r})$ and taking the temporal Fourier transform of Eq. A.1 gives

[TABLE]

where $k_{0}=\omega/c_{0}$. Note that for a uniform medium, then $\mu=1$, and Eq. (A.3) reduces to the homogeneous Helmholtz equation as expected. Now defining an auxiliary function $\lambda(\bm{r})\equiv k_{0}^{2}\left(1-\mu\right)$, Eq. A.3 may be written

[TABLE]

where the convolution theorem has been used, and $\Lambda(k_{x},k_{y},z)=\mathcal{F}_{k}[\lambda(x,y,z)]$. Then, since

[TABLE]

the left hand side of Eq. A.3 becomes

[TABLE]

Then Eqs. A.4 and A.6 yield Equation 4

[TABLE]

A.2 Implicit Solution and Approximation

Equation A.7 may be approached in a general way by treating the right hand side as a source term and using a Green’s function technique [55]. The homogeneous solution is

[TABLE]

and the appropriate Green’s function for the 1D homogeneous Helmholtz equation is[59, 60]

[TABLE]

Then the full solution is then given by

[TABLE]

If it is assumed that there are no backward travelling waves, then $B=0$, and from the boundary condition $P|_{z=0}=P_{0}$ at the source plane, the implicit solution is found [Eq. 5, repeated here]:

[TABLE]

The approximate marching scheme suggested by Jing et al. [56] defines the angular spectrum at some discrete axial position $z_{n}$ to be

[TABLE]

and

[TABLE]

Then, with knowledge of the initial condition $P^{0}=P_{0}$, the field may be approximated at arbitrary $z$ via

[TABLE]