Window size analysis in ultrasonic temperature estimation with timeshifts

In this work, ultrasonic timeshift method for temperature change estimation was investigated for 2D simulated in-silico synthetic ultrasonic signals. Digital phantom tissue was created in MATLAB environment and acoustic simulation was running on k-Wave toolbox for two different temperature conditions. First temperature distribution was assigned to tissue as uniform 37 °C. Second temperature distribution is Gaussian form with peak at tissue center as 45 °C and tails of Gaussian curve is 37 °C. Signal was analyzed with ultrasonic timeshift method for temperature change estimation. This method is based on four steps, calibration with tissue constant, finding timeshift with cross correlation algorithm, find slope of timeshift vector respect to timestep, and multiply tissue constant and slope of local timeshift vector. This multiplication gives temperature change of local point. In this work, window size of smoothing filter of timeshift vector and linear fitting to timeshift—timestep data was analyzed as parametrically with range 3λ to 10λ with a 1λ increment for both windows equally. As a result, window parameters as 5λ to 7λ give best results, maximum absolute error is 0.82 °C, 0.97 °C and 0.92 °C respectively and mean absolute error is ∼0.35 °C. As a verify, different analysis was performed on different temperature distribution with discrete two peak curves.


Introduction
Ultrasonic imaging technique-based temperature change estimation can be used in tumor ablation and hyperthermia treatments, to determine the success of the treatment during application.In heat therapies, tissue must be kept in a temperature range for a certain period which designed treatment protocol assessed by thermal dose models [1].With proper planning and imaging, the region of interest is heated sufficiently decrease size or destroy tumor region, and surrounding tissues are taken minimum damage from heating application.
Ultrasonic imaging technique is easy applicable modality in clinics and already used for monitoring detect ablation induced change in tissue and gas bubble formation in thermal therapies, [2,3].It is possible to estimate temperature change with a software-based method without or less need for additional instrumentation.It is advantageous in terms of cost and mobility compared to non-invasive MR thermometry.However, ultrasound thermometry is measuring temperature variation, not absolute temperature as in MR thermometry.Method is based on measuring the change in the time of flight of points inside the tissue in response to a change in temperature.When there is a temperature change in the tissue, the speed of sound changes and tissue expands with thermal expansion.Considering these two effects, temperature change estimation is performed based change of time-of-flight information in tissue, [4].It has shown that timeshift method is capable of measure with maximum absolute error as 0.53 °C and mean absolute error as 0.2 °C in Abdoulhassani's work, [5].
Beside timeshift method, there are different methods are used for ultrasonic temperature measurement which uses change of different properties of signal, i.e., backscatterer energy, photoacoustic, shear wave and Nakagami imaging, [6,7].Backscatterer energy method evaluates temperature estimation from change of signal energy which scattered signal from subwavelength scatterer in tissue for different temperature, [8].In photoacoustic imaging, laser burst is applied to tissue and tissue absorbs energy of photons.Tissue converts energy to heat and shows rapidly expansion and shrinkage which cause stimulation of wave in tissue which pressure depend on temperature, [9].Shear wave method also relates temperature and pressure which caused by shear wave from mechanical pulse to tissue, [10].In Nagakami imaging thermometry, as a statistical approach, local Nagakami distribution parameter is calculated, and tissue is mapped with this parameter.After temperature change it seen that local Nagakami distribution is changed, [11,12].
In this study, simulation was run for different temperature distributions on the 2D geometric model containing cylindrical scatterers.The temperature change estimations in the signals were obtained by using the timeshift method which run on simulation signals.Furthermore, the effect of the change in window sizes of analysis algorithm was investigated.In this work, it is aimed that testing window size parameters and find optimal value and show that the method is given adequate temperature change measurement.

2D tissue model and simulation
Ultrasonic signals used in this study were obtained in silico with running acoustic simulation by k-Wave [13,14] package program in MATLAB environment, [15].The computational area was modeled as 2D with distributed reflective and circular particles to form a heterogeneous structure.Spatial dimensions of study were designed as wavelength units for 1 MHz sinusoidal wave and 1540 m s −1 of speed of sound.Consequently, one wavelength was calculated as 1.54 mm.Computational area was designed as 128 for depth and 8λ for lateral width.Circular reflective particles were designed with 0.33λ of characteristic mean diameter and were randomly distributed in tissue ensuring acceptable minimum distance of 1.5 λ between center of particles.
The phantom has a uniform temperature at 37 °C initially.Subsequently, it is heated and has a temperature distribution with Gaussian profile which depending on the depth with the peak at 45 °C in the tissue center and tails at 37 °C as in figure 1. Gaussian distribution was created using the MATLAB 'gaussmf' command with a mean value at 64λ and a standard deviation of 0.0659.
Particle distribution on computational area and temperature distribution of second state are given in figure 2.
Thermoacoustic profiles in the computational area were assigned using the values given in table 1.The speed of sound profile for both the tissue and reflective particles depends on the temperature distribution within the computational area.Speed of sound values for the tissue region were selected based on literature data for human liver tissue.For reflective particles, speed of sound profile was created arbitrarily to show slight deviation from the liver speed of sound, while still providing sufficient reflection for the analysis of backscattered data.Accordingly, the speed of sound in the tissue was taken as 1580 m s −1 and 1420 m/s for reflective particles at 37 °C.It is assumed that there is a change in sound velocity of 1 m/s for every 1 °C increase for both regions.In the study, it was not aimed to model the tissue acoustic properties and behaviors exactly, thus it is assumed that approach with applying rough values for thermoacoustic properties is sufficient.Also, it is assumed that the density and attenuation profiles in the tissue and reflective regions are not temperature dependent, but model takes into account of attenuation of ultrasonic signals in simulation.

Acoustic simulation
The simulation was run separately for both temperature distributions using the 2D solver function using the kWave package under MATLAB environment, [14].In the study, no intervention was made in the governing equations of the kWave.In the simulation, the distance between the grids were taken as λ/10 and timestep were taken as 1/100 of the sinusoidal pulse frequency.These values were found in previous simulations to be sufficient to produce the sinusoidal signal waveform.Grids nodes which located at zero depth and 4λ to 6λ of lateral axis were defined as ultrasonic transducer and capable of transmit and receive acoustic signals.These grids transmit sine wave at 1 MHz frequency and 1 GPa amplitude of pressure for duration of 3 Period, (figure 3).

Boundary conditions and initial values
• At the initial moment, the acoustic pressure value in the tissue is zero.
• At depth X = 0, the input signal is applied for 4 to 6 λ at lateral axis for duration of pulse signal.
• A perfectly matched layer is defined on the four sides of the computational area.This layer absorbs waves passed from the edge of the calculation area.
As a result of the simulation, a time-pressure signal series -consisting of 40 lines represent the transducer grids-was obtained for each temperature condition.Gaussian average was applied to this signal series for each column.Thus, single line signal was obtained and given in the figure 4. Time axis was converted to wavelength using the c mean of tissue sound velocity value.
In these equations, T, k, c, α, z are respectively temperature, tissue constant, speed of sound, linear thermal expansion rate and depth.i is timestep represents both variable and indices, δ i is timeshift as timestep unit.In  this study, it is assumed that density is not depended on temperature, thus, thermal expansion rate is equal to zero.In figure 5  In the ultrasonic time-of-flight method, temperature changes between two different temperature distributions could be determined, not the absolute temperature in the tissue.For this reason, at least two ultrasonic imaging data are needed from the tissue, and local timeshifts could be detected and analyzed.Timeshift is the difference between the time of flight of any point in the tissue and is given in equation (3).In this study, timeshift values from two simulated signal which represents two temperature distribution calculated with cross-correlation method, xcorr in MATLAB.In cross-correlation operation, first signal is assigned as reference signal and both signals are equally windowed around interest of point by 3λ window size.Timeshifts are calculated from finding maximum correlation coefficient of correlation operation for two windows around interest of point.Indices of maximus minus window length gives timeshift value of point.In this study for each 300.timestep, timeshift data is calculated, and timeshift vector is evaluated.Since timeshift vector has noises and has outliner values, consecutive moving median and moving mean filter with equal window size -with 6λ width-is applied on vector and evaluated smooth data.Timeshift data graphic is given in figure 7, with timestep axis converted to spatial unit with c medium value for readilability.
Algorithm needs to evaluate d

D Di
i term for finding temperature changes in tissue.For finding this term, timeshift-timestep graph is evaluated and the window with 6λ width, as in timeshift calculation, is created around point of interest in graph.With linear fitting for values inside in windows local trend is calculated as d D Di i term, (figure 8).In results seen that maximum absolute error is found 0.97 °C and mean absolute error in heating region is found 0.34 °C.At peak of heated region, difference of maximum defined values and maximum values of estimation is 0.07 °C.

Further ınvestigatıon
Window sizes in filtering operation and linear fitting is important since window size affects error in temperature measurement.A window parameter of 6λ resulted in a maximum absolute error of ∼1 °C and a mean error of ∼ 0.35 °C.To achieve a finer solution, a parametric sweeping was also studied for window size.Window sizes were varied in the range of 3λ to 10λ with a 1λ increment for both the smoothing filter and linear fitting.Temperature errors were calculated for each parameter analysis, (figure 10).
According to the results, while using a window size of 5λ to 7λ, the maximum error falls within the range of 0.84 °C to 0.97 °C.Within the same size interval, the mean absolute error is approximately 0.34 °C, and the mean absolute error at the heating center is below 0.1 °C.Findings indicate that narrower window widths result in higher maximum absolute errors, particularly at 3λ and 4λ, exceeding 1 °C which shown in figure 11.Also,  increasing the window width leads to over smoothing in temperature measurement and causes an increase in peak temperature error at the heating center, although the mean absolute error remains at same levels.By evaluating with these results window size of 5λ to 7λ are given best results with less and acceptable errors in this study.For these window sizes, the ratio of error measurements which equal 0.5 °C or lower, to the total measurements was calculated as %91.4,%94.9 and %88.3 for 5λ, 6λ and 7λ respectively.
In further study, for verify the results, another analysis was held with signals obtained for different temperature distributions but same geometric model.In this study, first temperature state is uniform with 37 °C and second temperature distribution is defined with temperature change.Temperature change was modelled as two Gaussian curve which curves are located far from each other in geometry.First peak is located at near 40λ, and temperature change is 8 °C, and second peak is located at near 90λ, and temperature change is 4 °C.In analysis, window size for filter was selected as 6λ and temperature estimation is given in figure 12.
After analyzing second simulation signal which contain two temperature increase regions, it is found that maximum absolute error is 0.92 °C and mean absolute error is 0.29 °C.At each temperature peak region, difference of maximum values is below 0.1 °C.The ratio of measurement error which equal 0.5 °C or lower, to the total measurements was calculated as %92.7.In this study, window size parameter of smoothing of timeshift and linear fitting process are studied for finding optimal temperature measurement errors.Windows sizes are varied 3λ to 10λ with 1λ increment and parametric analysis are conducted.According to the results 5λ to 7λ give best results, maximum absolute error doesn't exceed 1 °C and mean absolute error is ∼0.35 °C.For these window interval, minimum ratio for error measurement equal or less than 0.5 °C to total measurements is 88.3%.At the heating center, measurements of peak value of estimation and peak temperature change value is found below than 0.1 °C.

Discussion
To verify findings, new analysis is conducted with 6λ window width with another temperature distribution.Similar results are given as 0.92 °C for maximum absolute error, 0.29 °C for mean absolute error and below than 0.1 °C for peak value difference.In these analysis, ratio of error measurements equal or below 0.5 °C to total measurements is found as %92.7.
Results shows that, temperature measurement for any point could be performed with 1 °C maximum error and ∼85% of measurements are performed with 0.5 °C or below temperature error with 5λ to 7λ window size.One of main aims of study is creating and testing algorithm in silico 2D model for reduce experiment costs with using before in vitro and clinical phases.With these results it seen that the algorithm could be tested on in vitro experiment.Source of errors in temperature estimation can be explained as simulation errors, speckle nature of ultrasound and filter type operations of algorithm.These sources are related some process in analysis of algorithm, i.e., simulation error causes different speckle pattern.Hence, instead of analysis each of part of algorithm, as an initial study, error analysis is held for whole algorithm.
Timeshift detection and temperature estimation is held on very high spatial frequency (for each ∼0.15-0.25λ)for assessing measurement capability of algorithm.It costs memory and time load, and time cost could probably limit real time applications.Running algorithm on region of interest, instead whole signal, this cost could be reduced.Also, limitation of instrumentation parameters as sampling frequency and frequency value, could change spatial frequency of measurement.It should be noted that, since estimation spatial frequency is ∼0.25λ, spatial resolution is window size parameter i.e., 6λ in this study because of study couldn't find temperature value of any point but interval defined by window size.Another limitation of method is rising from assessing of tissue constant in clinical studies.For each person, tissue constant could show different thermoacoustic profile due to structure, fat-water profile, or health condition.In thermal treatment, need to assess tissue constant value for each application.But it is easier in in-vitro and in-silico experiments due to calibration process.
The aim of this study is to investigate the effectiveness of temperature distribution on a 2D signal containing a basic speckle model and the variation of window parameters in the algorithm.Experimental studies for validation of algorithm will be conducted in future work and are not included in the scope of this study.For further studies, it is aimed that conducted 2-Dimensional simulation for performing temperature mapping of distribution instead of A-line data and conduct phantom experiments on tissue for testing algorithm.In real life thermal applications, due to bioheat behavior of tissue, two-dimensional temperature distribution is occurred as maximum at heating center and contours around the center with gradient to ∼37 °C.Hence, it is important to capture two-dimensional temperature distribution for assessing thermal damage to ROI and healthy surrounding area.

Figure 1 .
Figure 1.Gaussian temperature distribution in the phantom.

Medium 3 .
Speed of sound (37 °C) dc/dT (m/s•K) Density (37 °C) Acoustic Attenuation (db/cm) Reference Medium 1 (Liver) Estimation of temperature change In the case of a temperature change in the tissue, the flight time of the ultrasonic wave arriving at any point in the tissue varies depending on the change in the sound velocity and the thermal expansion of the tissue.In figure 5(a), timeshifts of regularly placed scatterers are shown in phantom geometry for arbitrary positive speed of sound change respect to temperature and no thermal expansion conditions after positive temperature change.By detecting the ultrasonic flight time and take into account change of sound of velocity and thermal expansion of tissue, the temperature change in the tissues could be detected.The temperature change in the biological tissue is formulated by Ebbini for continuous time reference, [4].Since computational domain needs to use discrete time domain, discrete equations have been derived at equations (1) and (2).
(b), timeshift-timestep graph is demonstrated for same scatterers which used in figure 5(a).

3. 1 .
Algorithm for analysisAnalysis algorithm is based on finding parts of equation(1).The first term on the right side of the equation is the k , medium the tissue constant.In this study, it is assumed that it does not change with temperature and determined using thermoacoustic properties of the tissue and reflective particle distribution.k medium is calculated with equation (2) as −1575.7 K, for geometry which minimum distance between center of particles is 1.5λ.While finding the second term dD Dii on the right side, firstly needs to find the d i term and secondly evaluating the timeshift-timestep graph.Determining the local trend of any point in graph with linear fitting using local window, gives us d D Di i term.

Figure 5 .
Figure 5. (a) Timeshift for regularly placed scatterers respect to different temperature distribution (b) Demonstration of timeshifttimestep graph for regularly placed scatterers.

Figure 6 .
Figure 6.Time of flight any point in signal and timeshift.

Figure 8 .
Figure 8. Derivation of local slope with linear fitting.

Figure 9 .
Figure 9. Temperature change estimation and absolute error.

Figure 10 .
Figure 10.Measurement error by window size.

Figure 11 .
Figure 11.Estimation of temperature change for different window size.
Ultrasonic timeshift method for temperature change estimation is based on two processes: finding timeshift and finding local slope of timeshift-timestep graphic.Finding timeshift is based on cross correlation operation.Since timeshift values show outliners and noises caused by acoustic interference and speckle, smoothing filter is applied.On the other hand, for finding local slope of timeshift-timestep graph, linear fitting method is applied with window around point of interest and values inside of window.Local slopes of graph are multiplied by tissue constant k medium and results are recorded as temperature change estimation of region of interest.