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.


Background
Transcranial focused ultrasound (tcFUS) using a largephased 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.

Results
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.
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.
The focus was spread over a large area and displaced from the intended target (Fig. 2b1, f1 and b2, f2). Two 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 phase correction methods refocused the distorted focus.

Discussion
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 softwarebased 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.
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 reposition 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).

Conclusions
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 preclinical 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. The central coordinate of 1024 elements (n) in a clinical focused ultrasound system, s n, 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-tobrain boundary, respectively. A geometric-centeroriented vector 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 Â d n;w Â h n;w : * ð1Þ In order to search the collision point, the corresponding trilinear interpolated voxel value Q n, 1 on the point s n, 1 was collected and compared with the predefined Hounsfield Unit threshold τ skull = 700. The integer increments d n, w iteratively increase its value by 1 until the following condition was met. The first collision point, s n, 1 , was localized by searching the end point which satisfy Q n, 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, s n, 2 , was defined by following the same detection method of the first collision point. However, the end point which satisfy Q n, 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 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 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.
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.
Where g x , g y and g z 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 n * n;1 ¼ g x ; g y ; g z ð4Þ

Refraction-based aberration phase
The incident angle θ I, ws can be calculated on the basis of the detected normal vector n * n;1 and the vector 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: , where θ I, ws and θ T, ws denote the incident angle and refraction angle on the water-to-skull boundary, respectively. And c w and c s 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 h n;s * that travels inside of the skull bone could be defined. The second collision point was searched based on the abovementioned method and the normal vector n * n;2 and incident angle θ I, sb were calculated accordingly. Then, the refraction angle θ T, sb and refracted vector h n;b * could be defined based on Snell's Law.
, where c b 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 d n, w , traverses the skull d n, s and the length from second collision point to the point that has the closest distance to targeted point d n, b were collected for each transducer element. The aberrated phase was then calculated using the following equation: , where f 0 is the driving frequency of the FUS system (650 kHz), and c n, w , c n, s and c n, 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 c n, b was defined as equal as the wave speed in water c n, 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, https://netbeans.apache. org) and available for free download from GitHub (https://github.com/jws2f/Kranion) 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].

Registration
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 FP Transducer 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 FP CT 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. The transformation function between the transducer coordinate system and CT volume was then calculated using singular-value decomposition (SVD) method [32]. The FP Transducer can be represented as 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 and the translation matrix T is defined as The final transformation between the voxel center in the CT image (P CT ) and the transducer coordinate system (P Transducer ) can be derived using 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.