Dense speed-of-sound shift imaging for ultrasonic thermometry

Objective. Develop a dense algorithm for calculating the speed-of-sound shift between consecutive acoustic acquisitions as a noninvasive means to evaluating temperature change during thermal ablation. Methods. An algorithm for dense speed-of-sound shift imaging (DSI) was developed to simultaneously incorporate information from the entire field of view using a combination of dense optical flow and inverse problem regularization, thus speeding up the calculation and introducing spatial agreement between pixels natively. Thermal ablation monitoring consisted of two main steps: pixel shift tracking using Farneback optical flow, and mathematical modeling of the relationship between the pixel displacement and temperature change as an inverse problem to find the speed-of-sound shift. A calibration constant translates from speed-of-sound shift to temperature change. The method performance was tested in ex vivo samples and compared to standard thermal strain imaging (TSI) methods. Main results. Thermal ablation at a frequency of 2 MHz was applied to an agarose phantom that created a speed-of-sound shift measured by an L12-5 imaging transducer. A focal spot was reconstructed by solving the inverse problem. Next, a thermocouple measured the temperature rise during thermal ablation of ex vivo chicken breast to calibrate the setup. Temperature changes between 3 °C and 15 °C was measured with high thermometry precision of less than 2 °C error for temperature changes as low as 8 °C. The DSI method outperformed standard TSI in both spatial coherence and runtime in high-intensity focused ultrasound-induced hyperthermia. Significance. Dense ultrasonic speed-of-sound shift imaging can successfully monitor the speed-of-sound shift introduced by thermal ablation. This technique is faster and more robust than current methods, and therefore can be used as a noninvasive, real time and cost-effective thermometry method, with high clinical applicability.


Introduction
High-intensity focused ultrasound (HIFU) utilizes concentrated ultrasound waves to generate therapeutic effects in a non-invasive manner.Thermal ablation and hyperthermia are applications where ultrasound is used to precisely target and deliver high focused energy in order to achieve localized heating of the targeted areas.This focused energy can effectively destroy tumors, making HIFU a valuable option for cancer treatment, especially for prostate, liver, and breast cancer (Wu et al 2003, Li et al 2004, Poissonnier et al 2007).Additionally, HIFU shows potential in non-oncological applications, such as reducing chronic pain, ablation of uterine fibroids and treatment of certain neurological disorders (Fry et al 1954, Chang et al 2015, Maloney and Hwang 2015, Moosa et al 2019).Hyperthermia has been used in various applications of brain tumor and carcinoma treatment as well as drug delivery (Guthkelch et al 1991, Staruch et al 2011, Gray et al 2019, Karmacharya et al 2021).HIFU for thermal ablation and hyperthermia offers a precise and targeted approach, minimizing damage to surrounding healthy tissues, and often requiring minimal recovery time.
As energy is released into the tissue and the local thermal dose rises, it is crucial to monitor and guide these treatments, to make sure that the target temperature is reached at the desired locations, with minimal off-target effects to surrounding healthy tissue.Magnetic resonance imaging (MRI) thermometry is the most widely used technique for HIFU-guidance (Kim 2015), (Ilovitsh et al 2019).The use of MRI for guiding HIFU procedures presents challenges in terms of accessibility, equipment compatibility, patient positioning, and treatment time, although such equipment is carefully designed to minimize the interference of HIFU on MRI.MRI machines are often large and expensive, limiting their availability in certain healthcare facilities.The strong magnetic fields generated by MRI can cause distortions in ultrasound beams, requiring specialized MRI-compatible HIFU systems for accurate targeting.Patient positioning inside the MRI scanner can be challenging due to the need for real-time adjustments and monitoring during HIFU treatments.Additionally, the combination of MRI and HIFU procedures can lead to longer overall treatment times, which may impact patient comfort and resource utilization.
Ultrasound-guided HIFU is an emerging adaptable alternative to MRI guidance due to the advantages of ultrasound in accessibility, real-time imaging, cost-effectiveness, and versatility.Ultrasound can provide continuous visualization for precise targeting and adjustments during the procedure.This real-time feedback enhances accuracy and ensures optimal delivery of focused ultrasound energy.One of the main challenges is the need to probe temperature change using ultrasound imaging for accurate HIFU-based thermometry.
Ultrasound pulse-echo imaging, also known as B-mode imaging, works by emitting short pulses of ultrasound waves into the body, which then reflect off tissue interfaces and are detected by a transducer to create an image, assuming a constant speed-of-sound.Heating changes the inherent speed-of-sound in the tissue, causing a measurable latency in the echo over successive ultrasound frames captured by an imaging transducer (figure 1).Originally, latency was measured by comparing the time it took for sound waves to cross the field of view, but this technique requires two transducers (a source and a receiver) (Sehgal et al 1986).Several works applied a simpler method relating the measured echo shift between B-Mode images directly to the induced temperature change (Seip and Ebbini 1995, Simon et al 1998, Liu and Ebbini 2010).This approach, later named thermal strain imaging (TSI), served as the basis for many of the more recent developments in the field of ultrasonic thermometry.
Notable methods to calculate the echo shift are cross-correlation template matching as well as Loupas' Estimator, but these are computationally burdensome calculations (Loupas et al 1995).Additionally, phase aliasing artifacts common to typical ultrasound imaging methods must be filtered after the echo shift is calculated, often incorporating heavy filters such as the median filter or Savitzky-Golay.Such artifacts greatly impede the calculation of thermal strain from the echo shift measurement.Thermal strain calculation is impartial to the underlying ultrasound acquisition method or heating procedure, and different imaging schemes have been implemented based on TSI, but these methods are still computationally expensive due to the calculation of echo shift by the TSI algorithm (Foiret andFerrara 2015, Ding et al 2016).If the calculation of the echo shift can be improved, the frame rate of the ultrasound acquisition can be raised to above 1 kHz, usually settling on a signal-to-noise ratio to frame rate tradeoff that yields 50-500 frames s −1 , by taking advantage of coherent plane wave compounding (Montaldo et al 2009).Methods taking advantage of plane wave compounding (Lee et al 2019) or GPU calculations (Ebbini et al 2018) have already been used to speed up TSI, but they suffer from the same limitations.Other methods also based on the change in tissue properties have also been used in ultrasound-guided HIFU, mostly based on the same algorithmic concepts (Shaswary et al 2021).Here, we propose a dense optical flow algorithm for echo shift calculation.The method offers reduced computation time and increases spatial coherence leading to a precise temperature estimation.Calculating the echo shift is conceptually identical to pixel-tracking, and echo shift has already been calculated (Mehrabani et al 2008) using the Horn-Shunck Dense optical flow algorithm (Horn and Schunck 1981).This approach has several advantages.First, this algorithm incorporates a spatial smoothness prior, forcing the resulting echo shifts to be spatially coherent as expected.In addition, dense optical flow can incorporate more intricate motion models than just translation and in many cases outperforms block matching algorithms (De Boer and Kalksma 2015).Finally, optical flow algorithms utilize pyramid calculation schemes to incorporate both low and high resolution information quickly.Today, Gunnar Farneback's (Farnebäck 2003) algorithm based on polynomial expansion is widely used to attain fast, spatially coherent pixel-motion results.In our work, we will combine fast plane wave imaging with optical flow to calculate the echo shift along several steering angles in real time, providing additional information about the field of view.
Once the echo shift has been measured, it can be used to calculate the temperature change.In TSI, this is done by differentiating the echo shift in the axial direction to produce thermal strain, which can be directly calibrated to the temperature shift.The relationship between strain and temperature shift is assumed to be linear (Miller et al 2002), though much work has been done to assess the validity of this assumption-concluding that the relationship between thermal strain and temperature change is actually quasi-linear (Civale et al 2013) or almost linear (Techavipoo et al 2004, Varghese andDaniels 2004).
Analysis of the axial displacement gradient as described above once again ignores the spatial coherence of the HIFU therapy.Instead, we propose unifying echo shift information from several plane waves to predict the change in speed-of-sound using an inverse problem method regularized with a spatial smoothness prior, which can similarly be calibrated to yield the underlying temperature shift.Regularized inverse problem approaches have shown promising results in ultrasound applications related to sound speed (Stähli et al 2020, Robins et al 2021), in addition to additional methods based on the beam geometry (Brevett et al 2022), passive reflectors (Sanabria et al 2019) or diverging waves (Rau et al 2021).
In this work, we developed a dense speed-of-sound shift imaging (DSI) method that combines the optical flow algorithm and inverse problem regularization to create a temperature estimating algorithm that can incorporate a wide array of information while forcing both the input and output to obey spatial regularization constraints.Utilization of prior knowledge about the spatial coherence of the HIFU surgery provides a faster and more robust result with no need for additional image processing.

Theory
In the following section, we describe the mathematical relationship between the predicted echo shift or change in round-trip time δτ and the temperature change δT along the wave propagation path.Previously, it was noted that the relationship of this echo shift to temperature and local sound speed is given by the following equation (Simon et al 1998): where the difference in round-trip time δτ is the cumulative sum of the inverse of changes in local sound speed c (z,T) due to a temperature change δT from the initial temperature T 0 in tissue with a thermal expansion coefficient of β(z) along the wave propagation path.
In thermal strain imaging, this integral is differentiated along the propagation path and several assumptions can be made to simplify the equation for calculation of temperature change from echo latency (Simon et al 1998): where α 1 is a medium-specific parameter describing the thermal expansion coefficient and linear relationship of sound speed with temperature, and can also be thought of as a unit conversion parameter from thermal strain to temperature.A method for ultrasonic estimation of temperature change from thermal strain based on equation (2) will start by calibration of the constant α 1 using a known ground truth temperature sensor.In this work we will assume a linear relation between the speed-of-sound and the change in temperature because we are interested in improving the estimation of sound speed shift, which is the underlying cause of thermal strain, rather than the calibration to temperature measurements.For water-based soft tissue, the constant α 1 has been reported to be −0.1% °C−1 change in thermal strain for a 1 °C increase in temperature (Seo et al 2011).
In our method, DSI, we suggest a new approach to solving equation (1).Rather than differentiate the line integral as is common in thermal strain imaging, we assume the problem to be linearizable.Therefore, equation (1) can be vectorized.The integral is represented as a matrix, denoted M below and the following equation is obtained: where the pixel shift Δd is observed in a plane wave B-mode image due to an echo shift of δτ.This Δd is the result of integrating the integrand Δν, which represents the change in sound slowness according to the following equation: Once the integral operand Δν is isolated, it is equivalent to the derivative of (1) such that equation (2) can be rewritten without the derivative to obtain: with α 2 assuming the tissue-dependent role of units conversion previously denoted in equation (2) as α 1 .In temperatures up to 40 °C, the relationship between sound speed and temperature change is mostly linear and α 2 should correspond to a change of up to 1 m s −1 for a 1 °C rise in temperature, depending on the target tissue (Seo et al 2011).
In practice, we want to measure the pixel shift, calculate Δν and estimate the underlying temperature change, so the following inverse problem is formulated: This problem is over-defined since many Δd images can be acquired from a particular field of view representing a single Δν.After transmitting multiple angled plane waves into the target tissue, a Tikhonov pseudo-inverse is utilized to calculate the optimal Δν satisfying the data.While this increases the potential for a reliable result, we further incorporate a smoothness prior to ensure spatial coherence.We chose to regularize the spatial gradient in the axial, lateral, and diagonal directions using L2 regularization.Thus, the detected change in sound slowness Δν is estimated by solving the following least squares problem: The pipeline of our DSI method utilizes the acquisition of ultrasound plane wave B-Mode images.The pixel shift at consecutive time intervals for a particular steering angle is computed via optical flow to produce the measured Δd vector.This vector is inverted as an inverse problem to receive Δν. Finally, the units of Δν are converted to temperature by a pre-calibrated constant (figure 1).

Materials and methods
3.1.Ultrasound setup HIFU was performed with a 2 MHz spherically focused single element transducer (H107, Sonic Concepts, Bothell, WA, USA) that was controlled by a transducer power output unit (TPO-200, Sonic Concepts, Bothell, WA, USA).This transducer has a diameter of 64 mm and is focused at a depth of 45 mm.Each transmitted pulse consisted of a sinusoid with a peak negative pressure of 3.25 MPa and a duty cycle of 50%.These parameters matched standard hyperthermia experiments (Lee et al 2019, Park et al 2019).To create a steady increase in temperature, the therapeutic transducer transmitted with a pulse repetition time of 20 μs with bursts that consisted of 20 cycles.HIFU treatment was performed in intervals of 45 s lasting a total duration of 270 s.After every interval, imaging was performed and the results stored before the next HIFU cycle began.
The therapeutic transducer was placed face-up at the bottom of a custom water tank.A holder was located at the focal spot, where the sample was placed.Imaging was performed using an L12-5 50 mm linear transducer (ATL Philips, WA, USA) placed perpendicular to the sample and coupled with ultrasound gel (figure 3(a)).A programmable ultrasound system (Vantage 256, Verasonics Inc., WA, USA) was used for imaging of a 59 mm × 40 mm field of view.A plane steering acquisition protocol was used to acquire nine plane waves per imaging sequence.Three main steering angles were selected at 5°intervals [−5°, 0°, 5°], each composed of three steering angles selected at 1.5°intervals (for example [−6.5°, −5°, −3.5°]).Raw RF data was stored for offline processing.
Two types of samples were used.The first was a tissue mimicking phantom.The phantom was made by heating a combined mixture of 1% agarose powder (A10752, Alfa Aesar, MA, USA) with deionized water until boiling.1% Silicon Carbide (57391, Sigma Aldrich, MO, USA) was added to the cooling mixture as acoustic scatterers, and the result was poured into a mold to congeal.Next, fresh ex vivo chicken breast samples were used.Each chicken breast sample was cut from the thickest part of a chicken breast in a single piece that could fill the entire field-of-view (pieces with a size of 59 mm × 40 mm × 15 mm 3 ).The samples were placed in a 3D printed cartridge at the focal spot of the transducer.To minimize reflections 6 mm rubber was placed around the sample and coupled with ultrasound gel.As a ground truth reference, a thermocouple was placed near the focal spot.The probe was operated from MATLAB during the imaging sequence through an Arduino Uno R3.
The pressure amplitude in this setup was calibrated using a needle hydrophone (NH0500, Precision Acoustics, Dorchester, UK) with an active aperture of 0.5 mm connected to an oscilloscope (MDO3024, Tektronix, OR, USA).The hydrophone was placed at the focal spot of the therapeutic transducer, such that the measured amplitude was maximal in the x, y, and z directions.
3.2.Post-processing and temperature change calculation RF data was beamformed and interpolated onto a 599 × 390 pixel 2 grid, such that the pixel size was 0.01 mm in both the axial z direction as well as the lateral x direction.Each image was then Hilbert transformed, and the 1.5°a ngled plane waves sub-sets were compounded such that each image acquisition provided three main steering angles.For example, the set [−6.5°, −5°, −3.5°] was compounded to yield a single image of −5°.This unique compounding has been used previously to enhance the signal-to-noise ratio in preparation for the pixel-tracking step (Stähli et al 2020).The dense optic flow was calculated between consecutive images at each of the main steering angles using Farneback's method (Farnebäck 2003).The three calculated pixel-shift images acquired at each time interval were down-sampled x5 in the lateral direction and x10 in the axial direction, then stacked to create the data vector.The matrix M described in equations ( 6) and (7) was pre-calculated and stored as a sparse matrix.Given a new data vector Δd, a simple matrix multiplication was performed to receive the slowness deviation.The regularization parameters, which are responsible for the smoothness of the result, were selected empirically to give the best possible resulting image.They can be fine-tuned based on the size and smoothness desired for an application of the algorithm.We chose to use regularization parameters Γx = 10, Γz = Γxz = 1.This means a certain level of smoothness was required for each of a pixel's 8 neighbors independently: Γx controls the smoothness with regard to the left-right neighbors, Γz controls the smoothness with respect to the pixel neighbors above and below the pixel, and Γxz controls the diagonal smoothness.The resulting image of pixel change in sound slowness was up-sampled back to the original image dimensions, then summed and compared to the expected temperature change to fit a linear curve, producing the calibration constant α 2 .
As a reference to DSI, TSI was also computed using the same data sets.Implementation included a 1D crosscorrelation search algorithm along each plane wave's propagation axis based on Lee et al (2019).The pixel shift was smoothed in the lateral direction with a Savitzky-Golay filter of length 13.4 mm and then the axial gradient was estimated with a second, gradient estimating Savitzky-Golay filter of length 15.6 mm.Each of the main angles described previously was analyzed separately, then the results were compounded for each time step.
Both algorithms were implemented in python.DSI's optical flow component made use of the OpenCV library, while the inverse problem solution was implemented as a sparse matrix multiplication using scipy.TSI, which requires a complex cross-correlation search, was implemented in numba as just-in-time compiled python code utilizing parallel processing to improve performance.

Simulated data
To validate our method, a simulated pixel-shift was induced in a set of steered plane wave images acquired from agarose tissue-mimicking phantom by manually introducing a pixel shift into the RF data captured from real images, implemented by translating pixels along each plane wave's propagation axis.A circular lesion mimicking a heated area was simulated within the images by introducing a virtual temperature shift of approximately 10 °C, at a circular area centered in the image with a radius of 5 mm (figure 2).This virtual heating was entirely controlled so that the algorithms could be compared to a known reference that was not affected by environmental factors like density, imaging parameters, or transducer alignment.In addition, agarose is less absorbent of ultrasonic energy than soft tissue, and we wanted to introduce a virtual heat source representative of soft tissue.The ground truth model is presented alongside TSI thermal strain calculation and DSI slowness deviation reconstruction (figures 2(b), (d), (f), respectively).The pixel shift of the TSI method provides a lateral full-width at half-max (FWHM) of 8.19 mm, while the DSI method gives a FWHM of 8.79 mm.The ground truth lateral FWHM is 8.49 mm (figures 2(a), (c), (e)).In terms of lesion dimensions, the lateral FWHM in TSI reconstruction is 10.56 mm while DSI achieves 9.57 mm (figure 2(g)).The ground truth was 9.77 mm in each direction.In the axial direction, TSI gives a FWHM of 7.89 mm while DSI gives 15.5 mm (figure 2(h)).Overall there is a smoother result with the DSI method.Next, the effect of the HIFU duration on the DSI algorithm performance was assessed (figure 4).Larger speed of sound shift was observed with an increase in treatment duration.The detected slowness deviation and focal area for each time step in figure 4 are displayed in table 1.
Treated region area increased by 48.6% for 90 s compared to 45 s, and by 85.3% for 135 s treatment (table 1).

Temperature estimation
The experiment presented in figure 4 was repeated 6 times.The mean temperature for each treatment interval was calculated.The average thermocouple measurements were used as a ground-truth model to represent the expected change in temperature.Sound slowness and thermal strain were extracted from the DSI and TSI algorithms respectively and calibrated to the ground truth model.Calibration of the temperature profiles yielded values of −0.22% °C−1 for α 1 and 88 ηs/m °C−1 for α 2 , corresponding to a change of 0.3 m s −1 for a 1 °C temperature increase.These values were used for calculating the temperature prediction (figure 5).

Performance benchmark
Finally, to emphasize the improvement in temporal resolution of the DSI algorithm, we performed a bench test to determine the average runtime of the algorithm compared to TSI.The cross-correlation parameters used in the search algorithm were chosen to optimize imaging quality as in figure 3. Four random experiments were selected, and each was processed 100 times with each algorithm.The processing time for each iteration was divided by the number of frames to calculate the mean time per frame.The average runtime of the TSI algorithm was 7.21 ± 0.5 s/frame, while DSI ran in 0.33 ± 0.004 s/frame.Overall a speedup of ×20 was achieved with DSI.

Discussion and conclusion
In this paper we developed a method for imaging a local speed-of-sound shift with a dense algorithm for the purpose of temperature estimation, aiming to improve the calculation time and spatial coherence of the result.The proposed DSI algorithm is robust and could be used to monitor additional applications that induce a speedof-sound shift.The common method for ultrasonic thermometry is TSI which is based on a cross-correlation search and local axial gradient.Our experiments showed several promising advantages for the DSI algorithm in HIFU guidance compared to TSI.First, across all experiments, it is clear that the smoothness priors used in the optical flow calculation and inverse problem regularization greatly improve the spatial coherence of the results.This is captured in figure 3, where the image processing steps that work well in the simulated TSI results do not generalize well to ex vivo, and the TSI pixel shift estimation gives a large focal area.In figure 3, heavy smoothing was necessary to reach an adequately smooth TSI image.One of the main advantages of DSI is not needing the additional filtering.In the future, knowledge of the shape of a HIFU beam in the transverse imaging plane can intuitively incorporate further apriori information into the DSI inverse problem solution to produce a more complex regularization scheme.This improvement can be further capitalized on by increasing the amount of plane waves to more heavily oversample the sound slowness field.Both of these concepts can be used to reduce the blurring caused by trivial L2 spatial regularization (figure 2(h)).
Another notable benefit of the DSI algorithm is its runtime.TSI and other correlation-based tracking methods are insufficient for real-time use when the region of interest is the entire frame.In this case it is preferable to make use of dense algorithms incorporating all of the available information simultaneously, while eliminating the need for image processing.Removal of the search step of TSI and replacement with optical flow provided a large part of the performance boost, but DSI's slowness estimation is also very fast due to the use of largely sparse matrices in solving the inverse problem.In ex vivo samples we detected a slowness deviation in order of 3 μs m −1 , which corresponds to a sound speed change of 7 m s −1 from a base sound speed of 1540 m s −1 .This translates to a speed-of-sound change of 0.45%.The temperature estimation results were similar between the ground truth thermocouple measurements, the TSI and DSI methods (figure 5).
Several considerations should be made when implementing the DSI algorithm.First, there are inherent tradeoffs in selecting the steering angles and regularization constants described in this algorithm.Selecting the sub-angles used for coherent compounding of the main angles is a direct tradeoff between SNR and resolution with frame rate, as in any coherent compounding scheme.The main angles play a much more prominent role in both processing time and resolution.Pixel tracking is performed on each angle individually, and the inverse problem calculation incorporates data from each of the main angles simultaneously, so the complexity of digital processing done in this study is linear with respect to the number of main steering angles.The scanning step of the plane waves must be selected so that high SNR is achieved by the sub-angles while adequately sampling the viewing field, with the goal to maximize the number of pixels in the region of interest being sampled by more than one plane wave.It is the cross-information acquired by separate plane waves on each pixel that improves the accuracy of the inverse problem solution.In addition, the RF sampling parameters must be fine-tuned for imaging during therapy.The most important consideration is the ultrasonic frequency: using a high frequency improves the resolution which is helpful for pixel tracking, but high frequency ultrasound comes with several  fallbacks.Increased tissue absorption associated with high frequencies causes the SNR to drop significantly at depth, and can also increase the thermal dose when used in between consecutive HIFU pulses.More expensive hardware is required to support the RF sampling rate required for such frequencies as well.Still, the ultrasound frequency should be maximized to improve imaging with these considerations in mind.
There are a number of limitations that should be mentioned.First, as in any ultrasound-based technique, the speed-of-sound shift can be measured only in soft tissue.In cases where the ablated region is located within bone, such as in transcranial ablation, MRI thermometry remains superior.The orthogonal alignment of the imaging and therapeutic transducers used in this study will not perform well with significant interior or exterior obstruction of the imaging transducer, for example during imaging through the ribcage or in areas of the body where the alignment of transducers is not possible.Instead, a dedicated USgFUS transducer like the one in Glickstein et al (2022) should be used.In addition, there is variability in temperature measurements that can be attributed to several factors (figure 5).Here the relationship between tissue temperature and sound speed change was modelled as a scalar unit conversion, but in practice it is a tissue-dependent curve that is linear only for certain temperature ranges (Seo et al 2011, Civale et al 2013) however this is not a new issue in our algorithm but rather a known drawback of temperature estimation methods based on equation (1).In the experiments described above we have performed only HIFU-induced hyperthermia treatments and taken care to stay within the linear region of this curve, but in thermal ablation therapy this temperature range is often crossed and the linearity is compromised.So long as the linear temperature region is maintained, we believe our algorithm can also be extended to more aggressive treatments such as thermal ablation therapy, at least as well as TSI.Another issue with temperature calibration are the absorptive and viscous heating artifacts common to thermocouples (Morris et al 2008, Civale et al 2013).As the tissue vibrates under the effect of HIFU, friction is generated between the tissue and the thermocouple, introducing viscous heat that is not present in true surgical HIFU applications.Additionally, the thermocouple itself absorbs heat from the surrounding tissue.Thus, temperature calibration of slowness deviation is forced to overcompensate for this heat in the thermocouple.Moreover, temperature calibration depends on the thermocouple positioning.A shift in positioning will add a built-in error to the process.Last, the heavy noise present in ultrasound imaging often confounds algorithmic applications, and both TSI and DSI are no exception.The estimation of echo delay is greatly impacted by the noise levels, and each algorithm has a different method of compensating.It is possible to further fine-tune these steps (image processing, or dense estimation) but we see the current implementation as a reasonable tradeoff between image quality and runtime.Finally, there is an important balance between spatial regularization and resolution.Here we have opted to introduce strong regularization, and figure 3 shows that this is preferable to the image processing done in TSI, but at the cost of spatial resolution.
Despite these difficulties, error in temperature measurement is as good as that of TSI, largely due to the high accuracy of the slowness deviation estimation.Our algorithm can provide reasonable, robust slowness deviation estimation in real-time applications, and more work can be done on the calibration to temperature for HIFU guidance.The algorithm can be implemented in real-time as shown by the performance benchmark, by storing the pseudo-inverse matrix, beamforming each plane wave angle independently, tracking the pixel motion, and solving the inverse problem with a simple matrix multiplication.We believe that the speed up in the algorithm run time could facilitate its implementation as part of a closed-loop HIFU controller or be used with matrix array transducers to enable real-time 3D thermometry.

Figure 1 .
Figure 1.A schematic diagram describing the DSI algorithm.Ultrasound B-Mode images are acquired at various plane wave steering angles and the pixel shift at consecutive time intervals is computed.This produces the measured Δd vector which is inverted to receive Δν. Finally, the units of Δν are converted to temperature by a pre-calibrated constant.

4. 2 .
Ex vivo results and comparison to TSI Next, ex vivo chicken breast experiments were conducted.HIFU was applied at the center of the sample, and an ablated region was visible at the center of the images of the sliced sample (figure 3(b)).A thermocouple was located at the focal spot to serve as the ground truth measurement.The imaging array acquired plane waves on the transversal plane of the sample.The recorded pixel shift and the reconstructed image were compared between the TSI and DSI methods (figure 3), and the focal area was defined as the number of pixels above half the maximal sound slowness deviation multiplied by the resolution in each axis.The lesion detected with TSI has a focal area of 226 mm 2 (figure 3(d)) while the DSI lesion recorded a focal area of 177 mm 2 .

Figure 2 .
Figure 2. Simulation results.Average axial pixel shift for (a) ground truth.(c) TSI, and (e) DSI methods.The reconstruction of the heated area is for (b) ground truth, (d) TSI, and (f) DSI methods.The reconstruction cross-sections are analyzed in the lateral (g) and axial (h) directions for the three methods.Axes are common to subfigures (a)-(f).Colorbar is common to subfigures (a), (c), (e).

Figure 3 .
Figure 3. Experimental comparison between the TSI and DSI methods.(a) Experimental setup illustration.(b) Image of an ex vivo chicken breast sample following a 270 s ablation treatment.Pixel shifts result for the (c) TSI and (e) DSI methods after the first treatment.Treated region reconstruction of the pixel shift maps using (d) TSI and (f) DSI methods.Axes are common to subfigures (c)-(f).Colorbar is common to subfigures (c), (e).

Figure 5 .
Figure 5.Comparison of Temperature estimation as a function of HIFU duration (N = 6 samples).Ground truth thermocouple measurements (blue line), TSI (orange line), and DSI (green line).Results are presented as mean ± STD.

Table 1 .
Detected slowness deviation and focal area of the treatment spot based on figure 4.