Open Access

Rapid full-wave phase aberration correction method for transcranial high-intensity focused ultrasound therapies

Journal of Therapeutic Ultrasound20164:30

DOI: 10.1186/s40349-016-0074-7

Received: 27 November 2015

Accepted: 13 October 2016

Published: 8 December 2016



Non-invasive high-intensity focused ultrasound (HIFU) can be used to treat a variety of disorders, including those in the brain. However, the differences in acoustic properties between the skull and the surrounding soft tissue cause aberrations in the path of the ultrasonic beam, hindering or preventing treatment.


We present a method for correcting these aberrations that is fast, full-wave, and allows for corrections at multiple treatment locations. The method is simulation-based: an acoustic model is built based on high-resolution CT scans, and simulations are performed using the hybrid angular spectrum (HAS) method to determine the phases needed for correction.


Computation of corrections for clinically applicable resolutions can be achieved in approximately 15 min. Experimental results with a plastic model designed to mimic the aberrations caused by the skull show that the method can recover 95 % of the peak pressure obtained using hydrophone-based time-reversal methods. Testing using an ex vivo human skull flap resulted in recovering up to 70 % of the peak pressure at the focus and 61 % when steering (representing, respectively, a 1.52- and 1.19-fold increase in the peak pressure over the uncorrected case). Additionally, combining the phase correction method with rapid HAS simulations allows evaluation of such treatment metrics as the effect of misregistration on resulting pressure levels.


The method presented here is able to rapidly compute phases required to improve ultrasound focusing through the skull at multiple treatment locations. Combining phase correction with rapid simulation techniques allows for evaluation of various treatment metrics such as the effect of steering on pressure levels. Since the method computes 3D pressure patterns, it may also be suitable for predicting off-focus hot spots during treatments—a primary concern for transcranial HIFU. Additionally, the plastic-skull method presented here may be a useful tool in evaluating the effectiveness of phase correction methods.


MRgHIFU Transcranial HIFU Phase aberrations Ultrasound simulation methods


High-intensity focused ultrasound (HIFU) is becoming prominent as a treatment method for a large number of diseases. The lack of ionizing radiation and non-invasive nature of HIFU makes it well suited as a primary or supplementary treatment for uterine fibroids [1], prostate cancer [2, 3], bone metastasis [4], and other afflictions [5]. In the brain, transcranial HIFU is being investigated for neurostimulation [6], to treat movement disorders [7] and gliomas [8], and may be useful in other treatments [9].

Transcranial HIFU, however, faces challenges not present for treatments in other areas in the body. The large attenuation (10 dB on average at 0.5 MHz and 20 dB at 1.5 MHz [10]) caused by the skull requires the use of large-aperture transducers to evenly spread the power over the skull to prevent skin burns [11]. Additionally, the large difference in acoustic properties of the skull and surrounding soft tissues—notably the speed of sound, which is approximately 2740 m/s in the skull [10] and near 1500 m/s in most soft tissues [12]—introduces phase aberrations that cause the focus to be distorted, displaced, and of lower intensity. Nevertheless, aberrations can be corrected by phasing the elements of a phased-array transducer such that constructive interference is achieved at the intended treatment location.

There exist several methods for determining the phases required to correct for aberrations [13]. The methods can be broadly classified into time-reversal, energy-based, and simulation-based methods. One version of the time-reversal method [14] relies on an implantable hydrophone to measure the phases from each element in order to determine the inverse phase for correction. Although this offers excellent correction, the invasive nature makes it inappropriate for HIFU treatments. Energy-based methods use imaging to determine the effect of the skull on the ultrasonic beam. For example, cavitation can create shock waves that can be imaged noninvasively (“ultrasonic stars” [15]), but the energy required to create cavitation and the bubbles themselves may pose a risk to the patient. Methods based on MR acoustic radiation force imaging [16] avoid cavitation but may require more time to converge to the correct phases. Simulation-based methods predict the phases through numerical simulation. Often, these simulations are based on acoustic properties derived from a CT scan [17], although using MR images of the skull obtained using UTE sequences may be possible [18]. Finite-difference time-domain (FDTD) simulations took about 2 h to compute phase corrections [19]. Faster simulations can be achieved by simplifying the acoustic models used [20], but the resulting phases will have reduced accuracy. It has been shown that full-wave simulation methods have an improved accuracy over ray tracing and can be computed within minutes [21]; however, the technique used there was limited to a single point of correction and did not include experimental results. Another work [22] has also shown full-wave methods to produce superior results compared to ray tracing at the expense of computational time (hours vs. seconds). A review of phase correction methods can be found in [13].

In this paper, we introduce a simulation-based full-wave phase aberration correction method. The method leverages the speed of the hybrid angular spectrum (HAS) method of acoustic simulation [23] in order to rapidly determine the phases required. HAS is implemented on a graphics processing unit (GPU) in order to quickly determine the acoustic pressure field in a 3D volume for each element, allowing for phase correction at multiple treatment locations with negligible increase in computational time.

We demonstrate the method’s effectiveness experimentally using two models. First, we correct the phase aberrations for a 3D-printed plastic aberrator model that is designed to simulate the phase aberrations caused by a human skull, but for which the acoustic properties can be accurately determined. Second, phase corrections are performed using a section of an ex vivo human skull. We further expand on the capabilities of this method by combining it with rapid HAS simulations to demonstrate and evaluate the effects of transducer/model misregistration and the effects of steering on pressure levels with reasonable calculation times.


Simulation technique

The phase aberration correction method in this paper simulates acoustic pressure patterns using the HAS method [23]. HAS alternates between the space and spatial-frequency domains as it propagates a steady state wave through transverse sections of a model to simulate 3D ultrasonic beam patterns; it simulates refraction, reflection, and absorption. HAS demonstrates a computational speed increase over other simulation methods. For example, a 301 × 301 × 300-voxel model with 0.15-mm isotropic resolution with three tissue types computed on a 2.67-GHz i7-Core Windows desktop with 12 GB of RAM using MATLAB 7.10 took 46 s in HAS versus 467 min using a FDTD technique using the same input source [23], though this did not exploit parallelism for either HAS or the FDTD technique. The HAS method models radiating boundary conditions, which was also modeled in the FDTD simulation.

The mathematical basis for the HAS approach is a discrete solution for the steady state Helmholtz equation
$$ {\nabla}^2A+{k}^2A=0, $$
where A is a pressure field and k is the wave number. Using the pressure p n-1(x, y), which is the pressure pattern at the entrance to a given slice of the model at layer n, the space-domain step for HAS is calculated as follows:
$$ {p}_n^{\mathit{\hbox{'}}}\left(x,y\right)={p}_{n-1}\left(x,y\right){e}^{j\varDelta b\left(x,y\right)r\mathit{\hbox{'}}}{e}^{-{a}_n\left(x,y\right)r}. $$
Here, r is the oblique distance between the slices taken at angles of the plane wave components and r′ is the perpendicular distance between the angular component and the z-axis. Δb(x, y) is a propagation term that represents the difference between a particular voxel’s propagation phase and the average propagation phase used in the spatial-frequency domain step (next). a n (x, y) is the attenuation for each voxel. After the space-domain step, the Fourier transform of p n (x, y) is computed yielding A′(α/λ, β/λ). In the spatial-frequency domain, the pattern is propagated using a transfer function
$$ {A}_n\left(\frac{\alpha }{\lambda },\frac{\beta }{\lambda}\right)={p}_n^{\mathit{\hbox{'}}}\left(\frac{\alpha }{\lambda },\frac{\beta }{\lambda}\right){e}^{j{b}_n{\textstyle \hbox{'}}\sqrt{1-{\alpha}^2-{\beta}^2}\Delta z}, $$

where α and β are the directional cosines (α = λf x and β = λf y ), b n is an average propagation constant across the plane, and Δz is the propagation distance from one plane to the next. The inverse Fourier transform of A n will be the pressure pattern at the face of the next plane, p n . The process can be repeated for each slice along the model.

HAS is able to account for inhomogeneities at a voxel-by-voxel level and, if needed, is able to account for reflections by calculating the impedance mismatch between voxels during the space-domain step. Although HAS is capable of modeling scattering [24], it was not employed for these simulations. HAS does not currently model nonlinearities or shear waves.

For phase correction, HAS is used to calculate the 3D pressure patterns throughout the entire model volume from each individual element of the phased-array transducer, initially assuming zero phase and uniform amplitude, as shown in Fig. 1a. Phase offsets for aberration correction are created by taking the negative of the phase from each element as found at any given treatment location. When this negative phase is impressed on each individual element, maximum constructive interference occurs at the desired point, as indicated in Fig. 1b. Arbitrary treatment locations can be handled by saving the phases over the entire treatment volume. There is negligible additional computational cost for multiple treatment sites because the 3D pressure pattern for the treatment volume has already been computed. Currently, the algorithm only accounts for the phase differences between elements and does not alter the amplitude of the elements when steering or at the geometric focus. However, since amplitudes can be saved, it is possible to implement more advanced methods such as those presented in [25] or [21].
Fig. 1

Phase correction method. a Ultrasonic waves are simulated from each individual element assuming zero phase, producing a 3D array of phases corresponding to any desired treatment location. Three arbitrary elements are shown here. b The inverse phase determined from simulations is applied to each element, resulting in zero phase at the treatment location for maximum constructive interference with greater intensity than in the uncorrected case

For fast computation that would be required in a clinical setting, each element’s pressure pattern is computed in parallel on a NVIDIA Tesla GPU (NVIDIA, Santa Clara, CA) using Jacket (ArrayFire, Atlanta, GA) and MATLAB 2012 (MathWorks, Natick, MA). For small acoustic models (229 × 159 × 182 voxels with six tissue types and 256 transducer elements), pressure patterns and phase corrections can be computed in approximately 45 s. Models such as the more clinically relevant skull flap model (421 × 648 × 170 voxels with 3000 distinct types of bone and 256 elements) take approximately 15 min for pressure patterns and corrections.

Aberrator model

Evaluation of the phase correction method can be difficult in an actual skull where there are uncertainties associated with determining the acoustic properties. To evaluate the effectiveness of the phase correction method independently of these errors, an aberrator model was developed that was made of 3D-printed photopolymer (VeroWhitePlus, Stratasys, Eden Prairie, MN). It was flat on one side and had randomly varying heights on the other side, as shown in Fig. 2. The heights on the aberrating side varied a maximum of 4 mm on top of a 7.5-mm base (giving a range of 7.5–11.5 mm total thickness). The variance in height was chosen to create phase shifts up to 2π, given the speed of sound of the plastic (Table 1). A separate rectangularly shaped block was used to determine the acoustic properties of the plastic by a through-transmission measurement [26].
Fig. 2

Phase aberrator height map variations in height of 3D-printed aberrator model. Heights were randomly generated to create phase shifts varying from 0 to 2π based on the speed of sound of the photopolymer used in printing (see Table 1)

Table 1

Acoustic properties of phase aberrator model

Speed of sound


2492 m/s

4.72 dB/(cm MHz)

Acoustic properties of the 3D-printed aberrator model as measured by a through-transmission test

Ex vivo skull model

To assess what level of corrections could be expected in an actual skull, an ex vivo section of human skull approximately 15 cm × 10 cm in size was evaluated. The skull flap was discarded, de-identified, and obtained from a deceased patient, and therefore IRB exempt. It was cleaned and frozen several months before being imaged in a clinical CT scanner at a resolution of 0.46 × 0.46 × 0.3 mm. An acoustic model of the skull was built using the CT images and a published method of mapping CT Hounsfield units (HU) to acoustic properties (speed of sound, density, and attenuation). The method, described in Pichardo et al. [27], maps a linear relationship between density and HU. The speed of sound and attenuation are determined from HU using a series of curves that were created by optimizing simulations to match experimental data. Before phase correction experiments, the skull flap was rehydrated overnight inside a degasser to remove air.

Hydrophone scans

Experimental verification of the phase correction method was performed using a 1-MHz 256-element phased-array transducer (IMASONIC SAS, Besançon, France) driven by high-power generators (IGT, Bordeaux, France). The transducer had a circular aperture of 14.5 cm in diameter and a focal length of 13 cm. A custom holder was built to hold the aberrator model and the skull flap at a fixed, known distance and rotation angle relative to the transducer. The distance and rotation of the holder were determined by placing the transducer and holder in an MR imager and performing a 3D GRE scan with 1-mm isotropic resolution. The resulting image was zero fill interpolated down to 0.25 mm isotropic voxel spacing, resulting in a misregistration error of 0.25 mm or less in each direction. Pressures were measured by scanning a hydrophone (HNR-0500, Onda Corporation, Sunnyvale, CA) in a raster pattern in the plane perpendicular to the direction of the beam propagation using two stepper motors (NRT150, Thorlabs Inc., Newton, NJ), as shown in Fig. 3. The signal from the hydrophone was passed through a 10-dB preamplifier before being recorded by a digital oscilloscope (PicoScope 4224, Pico Technology, Tyler, TX). These hydrophone scans had a resolution of 0.25 mm and covered an area of 10 × 10 mm. For phase correction, the hydrophone was left at the intended focus while each element was fired individually, allowing the phase to be recorded at that location. Later, the inverse of the recorded phases was applied to achieve the hydrophone time-reversal scans. The scans were performed both with and without phase corrections at the geometric focus.
Fig. 3

Experimental setup. Schematic of the experimental setup of hydrophone scanning system. The transverse plane was orthogonal to the direction of acoustic propagation (i.e., the transverse plane was horizontal)


Aberrator model

Simulations were run using the phase aberration correction method described in the previous section. Figure 4a shows the phases generated by the simulation-based phase correction method as a function of element location on the transducer. For comparison, the phases obtained from the hydrophone scans are shown in Fig. 4b. The time required to compute the simulation-based corrections was approximately 45 min due to the high resolution and size of the model (0.23 × 0.23 × 0.25 mm with 667 × 667 × 104 voxels), which was needed to emulate a smoothly varying aberrator.
Fig. 4

Phase maps. Computed phases obtained a through simulations and b with hydrophone measurements, for correcting aberration through the aberrator model to the geometric focus as a function of transducer element location

Pressure patterns at the geometric focus were obtained both with numerical simulations and with experimental hydrophone scans. As expected, the aberrator model created significant distortions in the pressure patterns in the absence of phase correction. Both the simulated and experimental pressure patterns demonstrated that without phase correction, the beam pressure was spread over a large area and displaced from the intended target, as shown in Fig. 5, upper row.
Fig. 5

Model pressure patterns. The left column displays the pressure patterns at the geometric focus through the aberrator model obtained from the hydrophone scans. The right column displays the results of the simulated pressure patterns using the acoustic parameters listed in Table 1. The rows from top to bottom represent the cases of uncorrected, corrected using simulated phases calculated with known acoustic parameters, and corrected using phases obtained from hydrophone measurements, respectively

Figure 5, middle and bottom rows, shows the corrected experimental and simulated pressure patterns though the aberrator model using phases found from both the hydrophone measurements and the simulation-based method.

Table 2 gives the relative increase in maximum pressure from the hydrophone scans due to the phase correction and the distance from the point of maximum pressure to the intended treatment location. Quantitatively, the observed peak beam pressure was increased by a factor of 1.41 when simulation-based corrections were applied. Phase corrections obtained with the hydrophone measurements resulted in a 1.50-fold peak pressure increase. Moreover, the distance of the maximum pressure from the intended treatment location was reduced from 1.9 to 0.56 mm or less in both cases.
Table 2

Phase correction improvements in the aberrator model


Ratio of maximum pressure corrected to uncorrected

Distance from focus (mm)




Corrected, simulation-based



Corrected hydrophone-based



Ratio of the highest pressure with phase corrections compared to the uncorrected case and target offset distance for the aberrator model as measured by hydrophone scans

Skull flap

Without phase correction, the ex vivo skull flap produced significant aberrations when the intended focus was at the geometric focus, as shown in Fig. 6, upper row. Displacement of the beam from the intended treatment location and spreading of the beam pattern were observed. Additionally, significant secondary side lobes were observed. An arrow in Fig. 6 indicates a side lobe with 84 % of the peak pressure of the main lobe for the hydrophone scan with no corrections.
Fig. 6

Skull flap pressure patterns. The left column displays the pressure patterns obtained with the hydrophone at the geometric focus through the skull flap. The right column displays the simulated pressure patterns. The rows, from top to bottom, represent the patterns with no corrections, corrections using phases obtained from simulations, and corrections using phases obtained from hydrophone measurements, respectively. The arrow in the uncorrected hydrophone scan points to a side lobe with 84 % of the peak main lobe pressure

After applying simulation-based phase corrections, significant improvements in the beam focus were observed. The computation time to compute the phases required for correction using simulations was approximately 15 min.

Quantitative data depicting the pressure levels and distance from the intended target from the hydrophone scans are listed in Table 3. The maximum pressure was increased by a factor of 1.51 for targeting at the geometric focus using simulation-based corrections. Additionally, the distance of the maximum pressure from the intended treatment location was reduced, the beam profile was smaller, and the secondary lobes were diminished. Hydrophone corrections resulted in a 2.17-fold increase in peak pressure compared to the uncorrected case.
Table 3

Phase correction improvements in ex vivo skull


Ratio of maximum pressure corrected to uncorrected

Distance from focus (mm)




Corrected, simulation-based



Corrected, hydrophone-based



Ratio of the highest pressure with phase corrections compared to the highest uncorrected pressure and target offset distance for the skull model as measured by hydrophone scans

The phase correction method can also be used to target multiple treatment locations with a negligible increase in computational costs. Figure 7 shows the phase correction applied to a location 5 mm transverse to the direction of the ultrasound propagation. Uncorrected, the pressure pattern showed a side lobe with 54 % of the peak pressure of the main lobe, as indicated by an arrow in the figure. Additionally, the distance of the maximum pressure from the intended treatment location was 2.34 mm. After correcting using phases from simulations, a peak pressure 119 % larger than the uncorrected case was found at a distance of 0.9 mm from the intended focus.
Fig. 7

Steering pressure patterns. Hydrophone scans of the phase correction method applied to the ex vivo skull flap with beam steering at a location 5 mm transverse to the geometric focus. The intended treatment location is at 0 mm in y and 5 mm in x. a Hydrophone scans with no phase correction. b Corrections using simulated phases. c Correction using hydrophone-measured phases. The arrow indicates a primary side lobe with 54 % of the main lobe pressure

Further simulation-based analysis

The HAS method combined with phase correction allows for rapid evaluation of a number of different parameters. Of particular clinical interest are the patient-specific responses to treatments—including the ability to electronically steer the beam while maintaining adequate heating capability and the sensitivity of the system to errors in misregistration.

Figure 8 shows the point-wise improvement when the beam is electronically steered for both the aberrator model and the skull flap. For each point, a simulation was performed using uncorrected and corrected phases when steering. This figure compares the ratio of the maximum pressure in the focal plane of the corrected to the uncorrected simulations as a function of the extent of steering.
Fig. 8

Ratio of steering improvement. Point-wise ratio of improvement of maximum pressures in the focal plane for the aberrator model (left) and the skull flap (right). Each point represents the ratio of the maximum corrected to uncorrected pressures in the focal plane when steering to a particular x, y location

Figure 9 shows the sensitivity of both the aberrator model and skull flap to misregistration using simulations with phase correction to the geometric focus. For each point, the model was translated by a specified amount in the plane perpendicular to the direction of the beam propagation while the phases were kept for the focus to mimic the effects of a misregistration. The maximum pressure in the focal plane for the shifted case was then compared to the maximum pressure assuming no shifting (i.e., at shift position 0, 0). In Figs. 8 and 9, the ratios were chosen to demonstrate the improvement independent of the power of transducer.
Fig. 9

Simulated misregistration sensitivity. Simulated sensitivity of the aberrator model (left) and skull flap (right) to translational misregistration. Each pixel represents the maximum pressure in the focal plane when the model has been shifted in a plane perpendicular to the direction of beam propagation by the specified distance

Although the results here cannot be generalized to other skulls, they provide a proof-of-concept of the ability of HAS in combination with phase correction procedures to rapidly simulate possible treatment scenarios. Each graph in Figs. 8 and 9 represents 529 individual simulations and can be created in less than 4 h using a Windows computer with an i5-core processor and 16 GB of RAM using MATLAB 2015a. Note that these simulations were not run in parallel—doing so would provide additional speed. Each individual simulation took approximately 20 s for the skull flap and 25 s for the aberrator model.

As a demonstration of the spatial variation in the propagation phases is related to the amount of misregistration, Fig. 10 shows the autocorrelation of 2D phase-length patterns for both the aberrator model and the skull flap. For this figure, the phase-length patterns were calculated by accumulating the phase shift encountered along the paths parallel to the axis of the ultrasound propagation. (For the skull flap, the phases were smoothed using a Gaussian kernel with a 2-pixel standard deviation to account for the CT measurement noise.) The phase lengths were then converted to phasor notation to best represent the potential for interference. A circular section of each pattern equal in size to the beam profile when it encountered the model was correlated with the untruncated pattern and normalized. This is representative of the overall spatial correlation between the phases required for phase correction.
Fig. 10

Autocorrelation of path lengths. The normalized autocorrelation of the phase lengths for the aberrator model (left) and skull flap (right) along the axis of the ultrasound propagation. Contour lines are displayed for values of 25 % (white), 50 % (red), and 75 % (black) of the maximum correlation


The phase correction method introduced here was able to rapidly calculate the phases required to correct for aberrations. In the case where the acoustic properties were fully known (for the 3D-printed phase aberrator model), the simulation-based correction was able to recover 95 % of the peak pressure achieved with hydrophone time reversal for targeting the geometric focus (as seen in Table 2). The relatively small difference between the simulation-based and hydrophone-based corrections for the aberrator model (95 % of the peak pressure was recovered) is most likely due primarily to registration errors. The simulations require close matching between the simulated and experimental placement of the aberrator with respect to the transducer elements. Misregistration will lead to poorer corrections. Figure 9 demonstrates this effect with translational errors in a single plane. Relatively small misregistrations can cause a significant decrease in the pressures seen. For example, simulations on the misregistration for the aberrator model show that a 1-mm translation of the model results in a 1.1–7.8 % loss in the peak pressure. Misregistration causes less of an effect for the skull model, which has smoother variations over the area through which the beam travels. Figure 10 further demonstrates this fact. When there is longer correlation distance for the phase lengths, one would expect misregistration in that direction to be less disruptive, a factor for phase correction. The correlation pattern for the skull flap shows longer correlation distances along the y-axis, which agrees with the data in Fig. 9 showing more tolerance to misregistration along that axis. Furthermore, the skull flap displays broader areas of correlation than the aberrator model, again agreeing with the results in Fig. 9 that show that the skull flap is more tolerant to misregistration overall. The 50 % contour in the results for the skull flap encompasses an area of 205 mm2 while the area for the aberrator model covers 24 mm2.

The hydrophone time-reversal method is performed experimentally with the aberrator in place and is not subject to the registration requirements. However, there may be slight errors in the location of the transducer element positions or the power output of each element. We estimate a misregistration error of less than 0.25 mm in each direction for these experiments. Data from the aberrator model, where the acoustic properties are homogeneous and known, agree with this: 95 % of the pressure is recovered, which also agrees with the simulated misregistration data presented in Fig. 9.

The 3D-printed aberrator model introduced here is a useful tool for assessing the ability of phase correction methods in the absence of acoustic parameter uncertainties. In the skull flap, where acoustic parameter uncertainties are present, simulation-based corrections were only able to recover 70 % of the pressure achieved using hydrophone time-reversal when targeting to the geometric focus (and 61 % when steering) compared to 95 % in the aberrator model. The aberrator model could be used to evaluate other possible methods of phase correction (ray tracing or FDTD time-reversal, for example) and to compare using different transducers or setups without needing to account for possible errors in acoustic modeling.

The difference between the percentage of pressure recovered using the simulation-based and hydrophone time-reversal methods of phase correction in the aberrator model (95 %) and the human skull flap (70 %) demonstrates the necessity for accurate acoustic parameter estimation of the skull. In transcranial HIFU, a large barrier to treatments is the necessary power required to achieve sufficient heating at the focal location without causing damage to unintended locations. Although the acoustic model of the skull flap used in this paper showed a 1.52-fold increase in peak power over uncorrected results, the hydrophone time-reversal phase correction demonstrated a 2.17-fold improvement. This suggests that there may be benefit in further research into acoustic parameter estimation of the skull. The aberrator model parameters were determined via a through-transmission test on a rectangular block of the same material used to 3D print the model, while the bone was modeled using the method presented in [27], which assigned values based on CT images. However, there is some evidence that the method presented in [27] may not be accurate in all cases as simulations using the method deviated in temperature by 24 % compared to the observed data, and the peak focal point distance between the simulations and experimental measurements was off by an average of 1.6 mm [28].

Other non-invasive methods have shown recovery of 14–58 % of peak intensity of hydrophone-based corrections using human skull flaps [29] compared to this method’s recovery of 49 % of the peak intensity (70 % of the peak pressure) at the geometric focus. However, caution must be used in directly comparing results for different transducer setups and skulls. Larger aperture size will usually cause the ultrasound to travel through a larger, more inhomogeneous region, yielding a greater need for phase correction. Additionally, a change in frequency of the transducer will change the precision required for phase correction. Unquantified variations in acoustic properties between human skulls also make direct comparisons difficult. Ideally, phase correction methods should be evaluated on the same system using an aberrator model similar to the one used here to evaluate effectiveness absent these disparities.

The method introduced here is more computationally efficient when compared to other methods presented in the literature. Simulations on a GPU for the skull flap modeled with clinically relevant resolution were performed in approximately 15 min. This computational time is for the full 3D volume for each individual element, allowing for corrections at multiple treatment locations. This computational time represents an eightfold improvement over a similar FDTD-based time-reversal simulation that was performed for only a single treatment location [19].

The method described has the benefit of simulating a full 3D volume. While we have used the phase information to demonstrate the possibility for multiple treatment locations, it would also be possible to use the amplitude information for predicting off-focus hot spots. For transcranial treatments, this could allow risk assessment of skin burns or heating in important brain structures. Combining the phase correction method with rapid simulations could allow for many clinically valuable insights, including determining which patients may not be suitable for transcranial HIFU treatments, treatment envelope evaluation, establishing limits for the maximum allowable power, or treatment planning.


We have presented a simulation-based method for phase correction that allows for improved focusing in transcranial HIFU applications. The method is rapid; phases can be computed in approximately 15 min for the resolution required for transcranial applications, and it creates corrections for a 3D volume allowing for treatment planning. The method was tested on both a 3D-printed acoustic aberrator and an ex vivo human skull flap where it resulted in increased pressure, decreased beam width, and reduction in side lobes. Corrections were achieved for both targeting to the geometric focus and steering to a target 5 mm away from the focus. Additionally, the method of corrections was combined with fast ultrasound simulations to evaluate the effects of steering and transducer/model misregistration on pressure levels. Since simulations are performed for a 3D volume, future work will include predicting maximum safe power usage and off-focus hot spots produced during treatments.



Finite-difference time-domain


Hybrid angular spectrum


High-intensity focused ultrasound


Hounsfield units



The authors wish to thank Henrik Odéen and Allison Payne for their insight and comments relating to this work as well as Robb Merrill for his help in designing and manufacturing the holder used for registration.


Funding for this study was provided by the Focused Ultrasound Foundation and NIH R01 grants EB013433 and CA172787.

Availability of data and materials

Please contact author for data requests.

Authors’ contributions

SA performed the experiments and drafted the manuscript. SA and DAC developed and implemented the phase correction method. SA, DAC, and DLP conceived of and designed the study. All authors read and approved the manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable

Ethics approval and consent to participate

Tissues used in this study were de-identified and obtained from a deceased patient, and therefore IRB exempt.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

School of Computing, University of Utah
Utah Center for Advanced Imaging Research, Department of Radiology, University of Utah
Department of Radiology, University of Utah
Department of Bioengineering, University of Utah
Department of Electrical and Computer Engineering, University of Utah


  1. Okada A, Morita Y, Fukunishi H, Takeichi K, Murakami T. Non-invasive magnetic resonance-guided focused ultrasound treatment of uterine fibroids in a large Japanese population: impact of the learning curve on patient outcome. Ultrasound Obstet Gynecol. 2009;34:579–83.View ArticlePubMedGoogle Scholar
  2. Murat F-J, Poissonnier L, Rabilloud M, Belot A, Bouvier R, Rouviere O, Chapelon J-Y, Gelet A. Mid-term results demonstrate salvage high-intensity focused ultrasound (HIFU) as an effective and acceptably morbid salvage treatment option for locally radiorecurrent prostate cancer. Eur Urol. 2009;55:640–9.View ArticlePubMedGoogle Scholar
  3. Poissonnier L, Chapelon J-Y, Rouvière O, Curiel L, Bouvier R, Martin X, Dubernard JM, Gelet A. Control of prostate cancer by transrectal HIFU in 227 patients. Eur Urol. 2007;51:381–7.View ArticlePubMedGoogle Scholar
  4. Li C, Zhang W, Fan W, Huang J, Zhang F, Wu P. Noninvasive treatment of malignant bone tumors using high-intensity focused ultrasound. Cancer. 2010;116:3934–42.View ArticlePubMedGoogle Scholar
  5. Jolesz FA, McDannold N. Current status and future potential of MRI-guided focused ultrasound surgery. J Magn Reson Imaging. 2008;27:391–9.View ArticlePubMedGoogle Scholar
  6. King RL, Brown JR, Newsome WT, Pauly KB. Effective parameters for ultrasound-induced in vivo neurostimulation. Ultrasound Med Biol. 2013;39:312–31.View ArticlePubMedGoogle Scholar
  7. Elias WJ, Huss D, Voss T, Loomba J, Khaled M, Zadicario E, Frysinger RC, Sperling SA, Wylie S, Monteith SJ, Druzgal J, Shah BB, Harrison M, Wintermark M. A pilot study of focused ultrasound thalamotomy for essential tremor. N Engl J Med. 2013;369:640–8.View ArticlePubMedGoogle Scholar
  8. McDannold N, Clement G, Black P, Jolesz F, Hynynen K. Transcranial MRI-guided focused ultrasound surgery of brain tumors: initial findings in three patients. Neurosurgery. 2010;66:323–32.View ArticlePubMedPubMed CentralGoogle Scholar
  9. Monteith S, Sheehan J, Medel R, Wintermark M, Eames M, Snell J, Kassell NF, Elias WJ. Potential intracranial applications of magnetic resonance–guided focused ultrasound surgery: a review. J Neurosurg. 2013;118:215–21.View ArticlePubMedGoogle Scholar
  10. Fry FJ, Barger JE. Acoustical properties of the human skull. J Acoust Soc Am. 1978;63:1576–90.View ArticlePubMedGoogle Scholar
  11. Clement GT, White J, Hynynen K. Investigation of a large-area phased array for focused ultrasound surgery through the skull. Phys Med Biol. 2000;45:1071.View ArticlePubMedGoogle Scholar
  12. Fairhead AC. ICRU report 61: providing reference data for tissue properties. J Acoust Soc Am. 1999;105:1324.View ArticleGoogle Scholar
  13. Kyriakou A, Neufeld E, Werner B, Paulides MM, Szekely G, Kuster N. A review of numerical and experimental compensation techniques for skull-induced phase aberrations in transcranial focused ultrasound. Int J Hyperthermia. 2014;30:36–46.View ArticlePubMedGoogle Scholar
  14. Fink M, Cassereau D, Derode A, Prada C, Roux P, Tanter M, Thomas J-L, Wu F. Time-reversed acoustics. Rep Prog Phys. 2000;63:1933.View ArticleGoogle Scholar
  15. Pernot M, Montaldo G, Tanter M, Fink M. “Ultrasonic stars” for time-reversal focusing using induced cavitation bubbles. Appl Phys Lett. 2006;88:34102-1-34102–3.View ArticleGoogle Scholar
  16. Vyas U, Kaye E, Pauly KB. Transcranial phase aberration correction using beam simulations and MR-ARFI. Med Phys. 2014;41:32901.View ArticleGoogle Scholar
  17. Aubry J-F, Tanter M, Pernot M, Thomas J-L, Fink M. Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans. J Acoust Soc Am. 2003;113:84–93.View ArticlePubMedGoogle Scholar
  18. Miller GW, Eames M, Snell J, Aubry J-F. Ultrashort echo-time MRI versus CT for skull aberration correction in MR-guided transcranial focused ultrasound: In vitro comparison on human calvaria. Med Phys. 2015;42:2223–33.View ArticlePubMedGoogle Scholar
  19. Marquet F, Pernot M, Aubry J-F, Montaldo G, Marsac L, Tanter M, Fink M. Non-invasive transcranial ultrasound therapy based on a 3D CT scan: protocol validation and in vitro results. Phys Med Biol. 2009;54:2597.View ArticlePubMedGoogle Scholar
  20. Clement GT, Hynynen K. Correlation of ultrasound phase with physical skull properties. Ultrasound Med Biol. 2002;28:617–24.View ArticlePubMedGoogle Scholar
  21. Kyriakou A, Neufeld E, Werner B, Székely G, Kuster N. Full-wave acoustic and thermal modeling of transcranial ultrasound propagation and investigation of skull-induced aberration correction techniques: a feasibility study. J Ther Ultrasound. 2015;3:11.View ArticlePubMedPubMed CentralGoogle Scholar
  22. Jones RM, Hynynen K. Comparison of analytical and numerical approaches for CT-based aberration correction in transcranial passive acoustic imaging. Phys Med Biol. 2016;61:23.View ArticlePubMedGoogle Scholar
  23. Vyas U, Christensen D. Ultrasound beam simulations in inhomogeneous tissue geometries using the hybrid angular spectrum method. IEEE Trans Ultrason Ferroelectr Freq Control. 2012;59:1093–100.View ArticlePubMedGoogle Scholar
  24. Christensen D, Almquist S. Incorporating tissue absorption and scattering in rapid ultrasound beam modeling. Proc SPIE. 2013;8584:85840X–1–7.View ArticleGoogle Scholar
  25. Gélat P, ter Haar G, Saffari N. A comparison of methods for focusing the field of a HIFU array transducer through human ribs. Phys Med Biol. 2014;59:3139–71.View ArticlePubMedGoogle Scholar
  26. Kak AC, Dines KA. Signal processing of broadband pulsed ultrasound: measurement of attenuation of soft biological tissues. IEEE Trans Biomed Eng. 1978;BME-25:321–44.View ArticleGoogle Scholar
  27. Pichardo S, Sin VW, Hynynen K. Multi-frequency characterization of the speed of sound and attenuation coefficient for longitudinal transmission of freshly excised human skulls. Phys Med Biol. 2011;56:219.View ArticlePubMedGoogle Scholar
  28. Pulkkinen A, Werner B, Martin E, Hynynen K. Numerical simulations of clinical focused ultrasound functional neurosurgery. Phys Med Biol. 2014;59:1679.View ArticlePubMedPubMed CentralGoogle Scholar
  29. Clement GT, Hynynen K. A non-invasive method for focusing ultrasound through the human skull. Phys Med Biol. 2002;47:1219–36.View ArticlePubMedGoogle Scholar


© The Author(s). 2016