Skip to main content

An open-source phase correction toolkit for transcranial focused ultrasound



The phase correction on transcranial focused ultrasound is essential to regulate unwanted focal point shift caused by skull bone aberration. The aim of the current study was to design and investigate the feasibility of a ray-based phase correction toolkit for transcranial focused ultrasound.


The peak pressure at focal area was improved by 140.5 ± 7.0% on target I and 134.8 ± 19.1% on target II using proposed phase correction toolkit, respectively. A total computation time of 402.1 ± 24.5 milliseconds was achieved for each sonication.


The designed ray-based phase correction software can be used as a lightweight toolkit to compensate aberrated phase within clinical environment.


Transcranial focused ultrasound (tcFUS) using a large-phased array transducer has become an attractive modality to treat many brain diseases. Initial clinical trials have reported about the treatment of brain tumor [1], neuropathic pain [2, 3], essential tremor [4, 5], Parkinson’s disease [6, 7], and blood–brain barrier opening in patients with Alzheimer’s disease [8]. Despite the substantial benefits of tcFUS, transcranial focusing of therapeutic ultrasound waves remains a challenge because of considerable differences between acoustic properties of the skull [9], the induction of acoustic focus distortion, and focal shift, combined with a significant decrease in focal intensity in the brain tissue [10].

The solution to these limitations is the precise modulation of the amplitude and phase of transmitted acoustic waves from each element of array transducers [11]. The time-reversal approach was introduced initially for phase aberration correction; this approach relies on an implantable acoustic reflector [12] or a mono-element transducer [13, 14], and provides the optimal aberration correction as the gold-standard for validating the performance of phase correction approaches during the development stage. However, in clinical applications, an invasive insertion of the reflector or hydrophone into the brain tumor would be required, and the possibility of infection and unwanted tissue damage would be unavoidable. Another proposed approach was using a minimally invasive technique called acoustic stars [15]. A cavitation signal from a natural nucleation site [16] or an injected liquid droplet [15] was studied and the corresponding time delay was computed for aberration correction. However, the invasive placement of an acoustic source and the collapse of the droplet should be carefully controlled to avoid hemorrhage.

Among various compensation techniques, the patient’s computed tomography (CT) image-based numerical simulation provides a clinically-feasible, noninvasive aberration compensation. A full-wave acoustic model implemented on a finite-difference time-domain method based on CT image derived acoustic property provided compensation with the considerable computation cost of a few hours [17]. A fast acoustic model using a hybrid angular spectrum [18] was reported, and the computation time of 15 min was achieved with compromised accuracy. Nevertheless, a series of sonication simulations was typically planned during the tcFUS treatment to enlarge the ablation volume [5] or to treat multiple targets [19]. However, the long computation time of numerical simulation approaches is currently prohibitive from a treatment standpoint. Computation time is crucial while considering the implementation of a compensation strategy for the clinical environment. A trade-off between compensation accuracy and computation time was reported on simulation [20] and experimental [21] studies.

In an attempt to develop a lightweight phase correction toolkit, the current study focuses on the design of an opensource software implemented ray-based phase correction on a graphics process unit (GPU) of a laptop computer. An integration of developed software to the current clinical tcFUS system was carried out and three human cadaver skulls were utilized to result the phase aberration. The performance of phase correction was characterized in terms of focal point location and focal pressure based on hydrophone scanning. The refocusing performance was validated by comparing with the clinical software-based compensation in local clinical tcFUS system.


The ray-based phase correction software introduced in this study was used to compute the phase correction in a sub-second speed, a total speed of 402.1 ± 24.5 milliseconds, given by the GPU’s computation capability. Two lines on the cross-sectional plane, as shown in Fig. 1a, which across the target point were scanned using hydrophone. Figure 1b exhibits the profile of the pressure field. The center of each figure was fixed to the peak point of the pre-scanned pressure map from the free field sonication. As expected, the ex vivo skull created distortion in the pressure pattern without phase correction. Kranion shows consistent compensation performance to keep the focal point on the targeted point based on these 1D pressure field scans.

Fig. 1

The illustration of the coordinate system and the 1D scanning results. a The illustration for coordinate of each axis. b Comparison of the one-dimensional pressure field scans of different aberration correction techniques. The measured peak negative pressure (PNP) values corresponding to the central line crossing the targeted point were collected and plotted. The results from transcranial (Skull C) focal pressure scanning without correction with ExAblate- and Kranion-based correction were illustrated. A sonication target was identified at the geometrical center of the therapeutic ultrasound transducer (Target I), and a 10 mm off-centered target was (Target II) used

Two orthogonal scanning planes, a cross-sectional plane (Fig. 2a1) perpendicular to the direction of the beam propagation and a longitudinal plane (Fig. 2e1) on the transducer axis, were implemented to capture the focal point in 2D. A white dot was marked on the free field sonication scan to indicate the maximum intensity point of peak negative pressure (PNP) and was utilized as a reference in the subsequent scanning and analysis. The skullcap was placed on the transducer after the focal point was localized in the scanning maps. Figure 2b1 and 6f1 present the particular aberration induced by Skull C.

Fig. 2

Hydrophone 2D scanning maps the human skullcap (Skull C). A 0.25 mm step resolution and 10 mm × 10 mm coverage area were maintained on all of the hydrophone scanned maps. Lateral (XY) and axial (XZ) hydrophone scanning maps based on the focal point (white dot) of the free field sonication were applied. The PNP map was plotted, and the white dashed line that crosses the focal point was illustrated. The peak intensity for each image is normalized for each local peak pixel value. The result for skull A and B were shown in additional files 1 and 2, respectively

The focus was spread over a large area and displaced from the intended target (Fig. 2b1, f1 and b2, f2). Two phase correction methods refocused the distorted focus. Figures 2 and 3 exhibit that Kranion demonstrated favorable performance on shifting the focal point back to the targeted point. The sonications of Target I resulted in average focal shifts of 0.44 ± 0.5, 0.34 ± 0.37, and 0.17 ± 0.15 mm for no correction, ExAblate-, and Kranion-based corrections, respectively. The sonications of Target II resulted in average focal shifts of 0.59 ± 0.44, 0.38 ± 0.34, and 0.16 ± 0.18 mm for no correction, ExAblate-, and Kranion-based corrections, correspondingly.

Fig. 3

The shifted focal point on Target I and Target II. The quantified focal point shift by skull aberration and the corrected focal point by using ExAblate and Kranion method on target I and target II

Table 1 lists the data depicting the peak pressure levels on the PNP map. Both correction methods improved the level of peak pressure by 152.8 ± 29.7% and 140.5 ± 17.0% on Target I and 162.6 ± 35.8% and 134.8 ± 19.1% on Target II for ExAblated and Kranion, respectively.

Table 1 Pressure on the focal point from 2D scanning for different correction methods


This study aimed to evaluate a developed toolkit (Kranion) for the compensation of skull acoustic phase aberration at a clinical applicable rate by hydrophone measurements using three skull cadavers in an ExAblate 4000 system. The phase correction method introduced in this study successfully corrected the aberrated focus to the intended target position. Although the amplitude correction was not implemented in this study, a similar peak focal pressure compared to the ExAblate software-based correction was achieved. Additional improvement of focal quality could be achieved by adding an amplitude correction module [11] on the basis of transmission loss and attenuation term [22].

The local skull density, thickness, and marrow thickness were varied in different locations and directly influenced the transmission amplitude of the therapeutic wave. In particular, the inefficiency of brain tissue heating had been reported on patients with thick marrow layers [23] in their skulls. The proposed phase correction approach in this study could have enhanced compensation performance in thick marrow cases and should be addressed in future studies. Moreover, in the current MRgFUS clinical protocol, the mean skull density ratio (SDR), which is the ratio between the mean values in Hounsfield units for the marrow and cortical bone [23], is used to screen patients for treatment. Recently, a skull bone injury was reported [24], but no correlation was found between the mean SDR and the presence or absence of skull lesion. This indicated that a potential local skull heating event could be present under sufficient sonication power even with an acceptable mean SDR value. This result suggests that the per-element-based local SDR should be considered in the amplitude correction in future studies.

Table 2 lists the different numbers of activated channels from two methods. As we understand, the clinical system is designed to decide if the channel needs to be turned on or off by comparing the incident angle and critical angle. Kranion has also implemented the similar strategy to make the decision. As Table 2 shows there are different activated channel numbers between Kranion and ExAblate. We believe, the main reasons for these differences are: by 1) the difference of incident angle calculated by the image processing algorithm from two methods and 2) differences in acoustic properties, such as sound speed and density, used in two methods to calculate the critical angle. In this study, we did not attempt to mimic the exact image processing approach as that used in the clinical system software to get the same activated channel number because the technical details in the clinical system are proprietary.

Table 2 Activated amount of phased array channels

Using ray tracing algorithm to improve focal quality was notable in this study. However, this study has the following limitations. The assumptions of the transmitting plane wave from the linearized wave equation in fluid media could not account for any shear wave, which is an important mode of wave propagation in solid skull bone. This limitation could be solved by adding a shear model to the ray method [25, 26]. The nonlinear propagation phenomenon under high-intensity FUS was disregarded in this study because only low-power sonications were tested. Thus, the applicable implementation of the phase correction proposed in this study might be limited to treatments using low-intensity FUS. The number of skulls used in this study is limited to 3 full-size hemisphere skulls which is the maximum amount we could collected for this study. It is desirable to have a greater number of skulls for better validation.

The single ray representation of the acoustic beam could result in inaccurate refraction on uneven and rough skull surfaces. We observed a very few rays which hit on the non-normal voxels on the skull surface and refracted to odd directions. The implementation of a bundle of rays could filter out the influence from the roughness of the skull surface and result in an averaged ray vector.

The hydrophone-based gold-standard phase correction was desirable in validating the performance of the proposed phase correction to determine the amount of energy recovered through the method. This validation should be investigated in future studies. Two static orthogonal 2D planes were scanned to measure focal point shift which may not sufficiently represent the exact focal shift after aberration in 3D. The actual focal point could be localized more accurately by performing small 3D scans near the anticipated focus for each case and re-position the 2D measurement planes. However, it may increase the overall scanning time considerably.

The fixed value of 2900 m/s for skull bone sound speed in this study was utilized to calculate the refraction angle. The first place where sound speed value of the skull bone was utilized in proposed algorithm is refraction angle calculation (5) which is right after the collision point detection. In this stage, the algorithm doesn’t have any knowledge about the skull bone sound speed to calculate the refraction angle. In order to apply the Snell’s law, we utilized the representative sound speed as an initial value for skull bone. Nevertheless, the heterogeneous property of the skull bone was accounted in later stage of the algorithm to calculate the aberrated phase (6).


A ray-based lightweight toolkit to compensate aberrated phase was developed. The performance of the toolkit was validated by hydrophone measurements through three skull cadavers in a clinical configuration. The accuracy and computation time of the proposed phase correction method are reasonable for its usage in pre-clinical transcranial focused ultrasound research.


Ray-skull collision detection

While an acoustic wave emitted from an ultrasound transducer element propagates along a beam in a medium, it reaches the boundary of adjacent medium. In transcranial focused ultrasound, the intersection points between the wave beam and the boundaries of different medium, such as water, skull bone, and brain tissue, should be defined to implement ray tracing. A threshold-based collision point detection was implemented on the pre-scanned CT volume images of the cadaver skull to identify the coordinates of the intersection points. The volume rendering of three cadaver skulls with fixing frame was reconstructed from scanned CT image as shown in Fig. 4.

Fig. 4

Reconstructed volume rendering of three skulls with their fixing frame. The volume of skull and fixing frame were reconstructed from CT volume images. The bottom, front and side view of three skull volume were shown in the figure. The visible fiducial points of skull A (seven points) were marked using red-doted circles. The variation of shape and surface curvature of each skull could be observed

The central coordinate of 1024 elements (n) in a clinical focused ultrasound system, sn, 0, was obtained from the system manufacturer. The notation n represents the element number of the array transducer and 0 is the indicator to represent the transducer coordinate. Similarly, the number 1 and 2 were indicated the coordinate of collision points on water-to-skull boundary and skull-to-brain boundary, respectively. A geometric-center-oriented vector \( \overset{\rightharpoonup }{h_{n,w}} \) represents the vector of the emitted wave illustrated as a black arrow in Fig. 5a. The ray was traversed with a step size of 0.1 mm. The end point of the vector is, which is the first collision point, defined as below:

$$ {s}_{n,1}={s}_{n,0}+0.1\times {d}_{n,w}\times \overset{\rightharpoonup }{h_{n,w}.} $$
Fig. 5

Overview of proposed ray-tracing based phase correction. a Some of the 1024 elements of the clinical focused ultrasound (FUS) transducer and reconstructed 3D patient head based on computed tomography (CT) images. A threshold-based collision point detection was applied to define the intersection point between incident vector and skull. b The refraction diagrams on water-to-skull and skull-to-brain boundary. c The edge operator on three axis and the equations to estimate the normal vector of the skull surface on the basis of the defined collision point. d Overview of developed graphic user interface (GUI) of the Kranion software. Illustrating the FUS transducer (d1) and magnetic resonance imaging (MRI)–CT registration (d2 and d3) of the same patient. e Zoomed view of the blue rectangular depicted in (d). The refracted rays (green lines) on sn, 1 and sn, 2 of 1024 channels were rendered in a sub-second speed. Note that the MRI image was not utilized in phase correction experiment. The MRI image shown in this figure is to show the MRI-CT image registration capability of developed software

In order to search the collision point, the corresponding trilinear interpolated voxel value Qn, 1 on the point sn, 1 was collected and compared with the predefined Hounsfield Unit threshold τskull = 700. The integer increments dn, w iteratively increase its value by 1 until the following condition was met. The first collision point, sn, 1, was localized by searching the end point which satisfy Qn, 1 ≥ τskull. Then the normal vector of the water-skull boundary was calculated based on localized collision point. The incident and refraction angles were calculated after the normal vector was defined. The details about incident angle, refraction angle and normal vector calculation were explained in following sections.

The second collision point, sn, 2, was defined by following the same detection method of the first collision point. However, the end point which satisfy Qn, 2 ≤ τskull was searched and defined as the second collision point. The cancellous bone, which has lower Hounsfield unit (HU) value than the threshold, was simply neglected by defining two computation regions based on two HU peaks, inner and outer cortical bones on collected HU profile along the corresponding ray segment.

Surface normal estimation

A normal vector of the skull surface that occurs on the ray–skull collision point is required to compute the incident and refraction angles of the wave beam (Fig. 5b). The edge of the surface is a small area in the image volume, where the local level changes rapidly. Thus, an edge operator was utilized to detect the best-oriented plane at estimated collision points.

The mathematical foundation for obtaining the three basic functions that estimate the local image gradient was introduced by Zucker and Hummel [27, 28]. A set of three 3x3x3 operators (ψi, ψj and ψk) are used to estimate the image gradient at any voxel of interest, each providing the directional gradient value in x, y and z directions respectively.

$$ {\displaystyle \begin{array}{c}{\psi}_i\left(x,y,z\right)=x/\sqrt{x^2+{y}^2+{z}^2}\\ {}{\psi}_j\left(x,y,z\right)=y/\sqrt{x^2+{y}^2+{z}^2}\\ {}{\psi}_k\left(x,y,z\right)=z/\sqrt{x^2+{y}^2+{z}^2\ }\end{array}} $$

The discrete approximation of these three operators were shown in Fig. 5c. Then the surface normal defining the best edge at collision point, i.e. (α, β, γ), was obtained by convolving ψ with the input image I.

$$ {\displaystyle \begin{array}{c}{g}_x=\left\langle {\psi}_i,I\right\rangle ={\iiint}_{\delta }{\psi}_i\left(x,y,z\right)I\left(\mathrm{x}-\alpha, y-\beta, z-\gamma \right) dx\ dy\ dz\\ {}{g}_y=\left\langle {\psi}_j,I\right\rangle \\ {}{g}_z=\left\langle {\psi}_k,I\right\rangle \end{array}} $$

Where gx, gy and gz are gradient component along x, y and z axis, respectively. Finally, the normal vector on the skull surface based on the detected collision point is defined as

$$ {\overset{\rightharpoonup }{\boldsymbol{n}}}_{n,1}=\left({g}_x,{g}_y,{g}_z\right) $$

Refraction-based aberration phase

The incident angle θI, ws can be calculated on the basis of the detected normal vector \( {\overset{\rightharpoonup }{\boldsymbol{n}}}_{n,1} \) and the vector \( \overset{\rightharpoonup }{h_{n,w}} \) of incident wave (see Fig. 5b). Then, the refraction angle θT, ws could be derived based on Snell’s Law as shown in bellow:

$$ \frac{\mathit{\sin}{\theta}_{I, ws}}{C_w}=\frac{\mathit{\sin}{\theta}_{T, ws}}{C_s} $$

, where θI, ws and θT, ws denote the incident angle and refraction angle on the water-to-skull boundary, respectively. And cw and cs are the averaged sound speed of water and skull, respectively. The sound speed for water and skull was defined as 1480 m/s and 2900 m/s, respectively [29]. Then, the refracted vector \( \overset{\rightharpoonup }{h_{n,s}} \) that travels inside of the skull bone could be defined. The second collision point was searched based on the above-mentioned method and the normal vector \( {\overset{\rightharpoonup }{\boldsymbol{n}}}_{n,2} \) and incident angle θI, sb were calculated accordingly. Then, the refraction angle θT, sb and refracted vector \( \overset{\rightharpoonup }{h_{n,b}} \) could be defined based on Snell’s Law.

$$ \frac{\mathit{\sin}{\theta}_{I, sb}}{C_s}=\frac{\mathit{\sin}{\theta}_{T, sb}}{C_b} $$

, where cb denotes the sound speed of brain region and it was set to be the same as the water in this study, because ex-vivo experiment was performed by merging skull into the water without any brain tissue or phantom placed inside the skull.

The length of the ray segments from transducer to first collision point dn, w, traverses the skull dn, s and the length from second collision point to the point that has the closest distance to targeted point dn, b were collected for each transducer element. The aberrated phase was then calculated using the following equation:

$$ {\varnothing}_n=2\pi {f}_0\left(\frac{d_{n,w}}{c_{n,w}}+\frac{d_{n,s}}{c_{n,s}}+\frac{d_{n,b}}{c_{n,b}}\right) $$

, where f0 is the driving frequency of the FUS system (650 kHz), and cn, w, cn, s and cn, b are the averaged sound wave speeds which were estimated by averaging the sound speed derived from the Hounsfield unit in the CT image [30] along each segment of the ray. Note that the wave speed in brain medium cn, b was defined as equal as the wave speed in water cn, w because no brain tissue was utilized in the experiment. The aberrated phase n of each element n was obtained, and additional phase unwrapping to define actual excitation phase range (as – π to π) for FUS transducer was performed. Phases were then printed out as a phase correction configuration file to be imported into the control computer (CPC) of the clinical system. Notably, incident angles that are greater than the critical angle (total reflection) were neglected in the phase computation, and the corresponding elements were turned off during the experiment.

Kranion software

By implementing the aforementioned technique, Kranion software was developed by the authors in [31], and Fig. 5d exhibits a screen capture of the GUI. Figure 5e illustrates a close-up view of the refracted rays. The software was developed using the Java programing language on the NetBeans integrated development environment (Apache Software Foundation, and available for free download from GitHub ( under MIT license. The main ray tracing computation was constructed as a compute shader executing on a GPU (GeForce GTX 1080 with 12 GB RAM, NVIDIA, CA, US). The software allows a user to manually register the MR and CT images and rotate the scene to any desired observation view. The FUS transducer geometry could be moved to an arbitrary position, and the corresponding ray tracing of 1024 elements was rendered in a sub-second speed and illustrated on the screen of a laptop computer. Additional information on the software is available in [31].


To simulate the skull phase aberration, the registration between transducer elements and the cadaver skull was obtained by using fiducial points through computer-aided design (CAD) fabrication and CT scanning. A cover and skull fixing frame were fabricated on the basis of the CAD sketch and assembled on top of the transducer. Figure 6a depicts a simplified sketch. Eight tantalum fiducial points (2 mm diameter, Bal-Tec, CA, US) were embedded in the known coordinates of the fixing frame during fabrication. The registration between the transducer and fiducial points was defined using an absolute geometry of the transducer and skull fixing frame. The coordinates of the fiducial points FPTransducer based on the transducer element coordinate system were then derived. Three cadaver skullcaps were obtained (University Hospital of Virginia, VA, US) and fastened on the fixing frame. A 3D volume image of the assembled skull frames was obtained using the LightSpeed VCT CT system (GE Healthcare, IL, US). The coordinates of the fiducial points FPCT were defined by binarizing the CT image of each skull using an appropriate threshold. Granulometry was used to identify the fiducial regions, and the center of mass was calculated for each region to obtain the fiducial position.

Fig. 6

The registration of clinical focused ultrasound system and cadaver skull. a Computer-aided design (CAD) sketch of the clinical transducer, skull fixing frame, and reservoir cover. b Fiducial point-based registration of the transducer elements and skull volume in scanned CT images. Notably, the detailed geometry of the actual FUS transducer was excluded in this study. The sketch in (a) and (b) does not represent the actual scale of the FUS transducer. c Top view illustrating the coordinate configuration, including patient perspective coordinate (patient’s anterior-PA & patient’s posterior-PP, patient’s left-PL & patient’s right-PR), CT image coordinate, and coordinate of hydrophone scanning system (front [F] & back [B], left [L] & right [R])

The transformation function between the transducer coordinate system and CT volume was then calculated using singular-value decomposition (SVD) method [32]. The FPTransducer can be represented as

$$ SVD\left({FP}_{Transducer}\right)= U\varSigma {V}^{\ast }, $$

where U (m × m) and V (n × n) are the unitary matrices, and Σ is the diagonal matrix (m × n). V is the conjugate transpose of V. Then, the rotation matrix R is expressed as

$$ R=U{V}^{\ast } $$

and the translation matrix T is defined as

$$ T=-R\bullet {FP}_{Transducer}+{FP}_{CT} $$

The final transformation between the voxel center in the CT image (PCT) and the transducer coordinate system (PTransducer) can be derived using

$$ {P}_{CT}=R\bullet {P}_{Transducer}+T $$

Figure 6b plots the composite figure of the CT skull frame image and the transducer elements. To obtain an improved see-through view, only half of the transducer is shown. Figure 6c depicts the coordinate configuration lookup table between the software and experiment hardware systems. This configuration combines the coordinate systems of the transducer, patient, and hydrophone scanning domain.

Experimental setup and hydrophone scanning

The experiment was implemented on the magnetic imaging guided focused ultrasound (MRgFUS) clinical system at the Focused Ultrasound Center of the University of Virginia (University of Virginia Health System, VA, US). An ExAblate 650 kHz FUS transducer was utilized in this study. The transducer was carefully detached from the treatment bed and placed in an AIMS III scanning tank (AST3-L, ONDA, CA, US). The facing up transducer was mounted on a constructed frame that was fixed inside the scanning tank (Fig. 6c). The transducer filled with degassed water was positioned in an empty scanning tank during the whole experiment. The degassed water (with the remaining oxygen level maintained at less than 1.1 ppm) was prepared using collected tap water and a vacuum membrane degassing system before configuring the experimental setup. Skull degassing was performed for 1 h before attaching the skull frame to the FUS transducer.

The clinical brain transducer consists of 1024 elements, and each element is individually driven on the basis of predefined phase correction files. The format of the phase correction file that is compatible with ExAblate console software was studied and a new configuration file was reconstructed in accordance with the proposed phase correction calculated by the Kranion software. The Kranion- and ExAblate-based phase correction files were fed into the CPC, and the corresponding sonications were scheduled to compare the resulting pressure fields.

Figure 7 demonstrates an overview of the experimental setup. A needle-type hydrophone (HNA-0400, ONDA, CA, US) and a preamplifier (AH2020, ONDA, CA, US) were utilized to measure the pressure variation in the region of interest. The amplified pressure signal was obtained using an oscilloscope that was synchronized to the transducer driving system through a Bayonet Neill-Concelman connector. The raw data of each 2D pressure scan were saved and labeled in the PC that controlled the hydrophone positioning system.

Fig. 7

The diagram of experiment workflow. a The schematic setup for hydrophone scanning experiment. The phase correction file was constructed in the Kranion software on the basis of the pre-scanned CT image of each collected cadaver skull. The software in the workstation of the clinical system was scheduled for phase-corrected sonication on the basis of the imported correction file. A series of hydrophone scanning was performed on the focal plane to map the transcranial pressure field. The 2D pressure field data of each scanning were collected for postprocessing

Availability of data and materials

The opensource toolkit in this paper is available from and the datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.



Transcranial focused ultrasound


Computed tomography


Graphics process unit


Focused ultrasound


Graphic user interface


Magnetic resonance imaging


Hounsfield unit


Control computer


Computer aided design


Singular-value decomposition


Magnetic imaging guided focused ultrasound


Peak negative pressure


Skull density ratio


  1. 1.

    Coluccia D, et al. First noninvasive thermal ablation of a brain tumor with MR-guided focused ultrasound. J Ther Ultrasound. 2014;2:17.

    Article  Google Scholar 

  2. 2.

    Martin E, et al. High-intensity focused ultrasound for noninvasive functional neurosurgery. Ann Neurol. 2009;66(6):858–61.

    Article  Google Scholar 

  3. 3.

    Jeanmonod D, et al. Transcranial magnetic resonance imaging-guided focused ultrasound: noninvasive central lateral thalamotomy for chronic neuropathic pain. Neurosurg Focus. 2012;32(1):E1.

    Article  Google Scholar 

  4. 4.

    Elias WJ, et al. A pilot study of focused ultrasound thalamotomy for essential tremor. N Engl J Med. 2013;369(7):640–8.

    Article  Google Scholar 

  5. 5.

    Lipsman N, et al. MR-guided focused ultrasound thalamotomy for essential tremor: a proof-of-concept study. Lancet Neurol. 2013;12(5):462–8.

    Article  Google Scholar 

  6. 6.

    Magara A, et al. First experience with MR-guided focused ultrasound in the treatment of Parkinson’s disease. J Ther Ultrasound. 2014;2:11.

    Article  Google Scholar 

  7. 7.

    Fasano A, et al. Magnetic resonance imaging-guided focused ultrasound thalamotomy in Parkinson tremor: reoperation after benefit decay. Mov Disord. 2018;33(5):848–9.

    Article  Google Scholar 

  8. 8.

    Lipsman N, et al. Blood–brain barrier opening in Alzheimer’s disease using MR-guided focused ultrasound. Nat Commun. 2018;9(1):2336.

    Article  Google Scholar 

  9. 9.

    Fry FJ, Barger JE. Acoustical properties of the human skull. J Acoust Soc Am. 1978;63(5):1576–90.

    Article  Google Scholar 

  10. 10.

    Sun J, Hynynen K. Focusing of therapeutic ultrasound through a human skull: a numerical study. J Acoust Soc Am. 1998;104(3 Pt 1):1705–15.

    Article  Google Scholar 

  11. 11.

    White J, Clement GT, Hynynen K. Transcranial ultrasound focus reconstruction with phase and amplitude correction. IEEE Trans Ultrason Ferroelectr Freq Control. 2005;52(9):1518–22.

    Article  Google Scholar 

  12. 12.

    Flax SW, O'Donnell M. Phase-aberration correction using signals from point reflectors and diffuse scatterers: basic principles. IEEE Trans Ultrason Ferroelectr Freq Control. 1988;35(6):758–67.

    Article  Google Scholar 

  13. 13.

    Clement GT, Hynynen K. Micro-receiver guided transcranial beam steering. IEEE Trans Ultrason Ferroelectr Freq Control. 2002;49(4):447–53.

    Article  Google Scholar 

  14. 14.

    Pernot M, et al. In vivo transcranial brain surgery with an ultrasonic time reversal mirror. J Neurosurg. 2007;106(6):1061–6.

    Article  Google Scholar 

  15. 15.

    Haworth KJ, et al. Towards aberration correction of transcranial ultrasound using acoustic droplet vaporization. Ultrasound Med Biol. 2008;34(3):435–45.

    Article  Google Scholar 

  16. 16.

    Gateau J, et al. Transcranial ultrasonic therapy based on time reversal of acoustically induced cavitation bubble signature. IEEE Trans Biomed Eng. 2010;57(1):134–44.

    Article  Google Scholar 

  17. 17.

    Pulkkinen A, et al. Numerical simulations of clinical focused ultrasound functional neurosurgery. Phys Med Biol. 2014;59(7):1679–700.

    Article  Google Scholar 

  18. 18.

    Almquist S, Parker DL, Christensen DA. Rapid full-wave phase aberration correction method for transcranial high-intensity focused ultrasound therapies. J Ther Ultrasound. 2016;4:30.

    Article  Google Scholar 

  19. 19.

    Ghanouni P, et al. Transcranial MR-guided focused ultrasound: a review of the technology and neuro applications. AJR Am J Roentgenol. 2015;205(1):150–9.

    Article  Google Scholar 

  20. 20.

    Kyriakou A, et al. 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.

    Article  Google Scholar 

  21. 21.

    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(1):23–36.

    Article  Google Scholar 

  22. 22.

    Pinton G, et al. Attenuation, scattering, and absorption of ultrasound in the skull bone. Med Phys. 2012;39(1):299–307.

    Article  Google Scholar 

  23. 23.

    Chang WS, et al. Factors associated with successful magnetic resonance-guided focused ultrasound treatment: efficiency of acoustic energy delivery through the skull. J Neurosurg. 2016;124(2):411–6.

    Article  Google Scholar 

  24. 24.

    Schwartz ML, Yeung R, Huang Y, et al. Skull bone marrow injury caused by MR-guided focused ultrasound for cerebral functional procedures. J Neurosurg. 2018;130(3):758–62.

  25. 25.

    Pichardo S, Hynynen K. Treatment of near-skull brain tissue with a focused device using shear-mode conversion: a numerical study. Phys Med Biol. 2007;52(24):7313–32.

    Article  Google Scholar 

  26. 26.

    Jones RM, O'Reilly MA, Hynynen K. Transcranial passive acoustic mapping with hemispherical sparse arrays using CT-based skull-specific aberration corrections: a simulation study. Phys Med Biol. 2013;58(14):4981–5005.

    Article  Google Scholar 

  27. 27.

    Zucker S, Hummel R. A 3-dimensional edge operator. In: Proc. IEEE Conf. on pattern recognition and image proocessing; 1979.

    Google Scholar 

  28. 28.

    Ballard DH, Brown CM. Computer vision, 1st edition, Prentice Hall, ISBN-10:0131653164, 1982; p. 81–3.

  29. 29.

    Aubry JF, et al. Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans. J Acoust Soc Am. 2003;113(1):84–93.

    Article  Google Scholar 

  30. 30.

    Marsac L, et al. Ex vivo optimisation of a heterogeneous speed of sound model of the human skull for non-invasive transcranial focused ultrasound at 1 MHz. Int J Hyperth. 2017;33(6):635–45.

    Article  Google Scholar 

  31. 31.

    Snell J, Anders. Kranion software: Focused Ultrasound Foundation; 2017.

  32. 32.

    Strang G. Linear Algebra and its applications. 4th ed. ISBN-10:0030105676; 2006.

    Google Scholar 

Download references


The authors thank Anders Quigg for his inputs on the development of the Kranion software during its early stage.


This publication was partially supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT, No. 2018R1A2B2007997). All the experimental works and data collection and analysis were supported by Focused Ultrasound Foundation (FUSF) when the first author was an intern and the corresponding author was a Merkin fellow at FUSF.

Author information




CJ developed the software, performed validation experiment, analyzed the data, and was a major contributor in writing the manuscript. JS developed the infrastructure of the software and operated the clinical system for validation experiment. DM conducted cadaver skull preparation and CT scanning. DP designed and coordinated the study, and helped to draft the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Dong-Guk Paeng.

Ethics declarations

Ethics approval and consent to participate

All cadavers were sourced from the Virginia State Anatomical Program (VSAP) via the University of Virginia Medical School. Cadavers are not considered as a human subject in VSAP and no IRB approval was required for the current study.

Consent for publication

Not applicable.

Competing interests

The authors declare no conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1.

Hydrophone 2D scanning maps of the free field sonication, skull aberration without correction, with ExAblated- and Kranion-based corrections on the human skullcap (Skull A). A 0.25 mm step resolution and 10 mm × 10 mm coverage area were maintained on all of the hydrophone scanned maps. Lateral (XY) and axial (XZ) hydrophone scanning maps based on the focal point (white dot) of the free field sonication were applied. The PNP map was plotted, and the white dashed line that crosses the focal point was illustrated. The peak intensity for each image is normalized for each local peak pixel value.

Additional file 2.

Hydrophone 2D scanning maps of the free field sonication, skull aberration without correction, with ExAblated- and Kranion-based corrections on the human skullcap (Skull B). A 0.25 mm step resolution and 10 mm × 10 mm coverage area were maintained on all of the hydrophone scanned maps. Lateral (XY) and axial (XZ) hydrophone scanning maps based on the focal point (white dot) of the free field sonication were applied. The PNP map was plotted, and the white dashed line that crosses the focal point was illustrated. The peak intensity for each image is normalized for each local peak pixel value.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Jin, C., Moore, D., Snell, J. et al. An open-source phase correction toolkit for transcranial focused ultrasound. BMC biomed eng 2, 9 (2020).

Download citation


  • Transcranial focused ultrasound
  • Ray-based phase correction
  • Open-source toolkit