Uncertainty estimation for temperature measurement with diagnostic ultrasound
 Tina A. Fuhrmann^{1}Email author,
 Olga Georg^{2},
 Julian Haller^{2},
 KlausV. Jenderka^{1} and
 Volker Wilkens^{2}
DOI: 10.1186/s403490160071x
© The Author(s) 2016
Received: 1 July 2016
Accepted: 1 October 2016
Published: 1 December 2016
Abstract
Background
Ultrasound therapies are promising, noninvasive applications with potential to significantly improve, e.g. cancer therapies like viro or immunotherapy or surgical applications. However, a crucial step towards their breakthrough is still missing: affordable and easytohandle quality assurance tools for therapy devices and ways to verify treatment planning algorithms. This deficiency limits the safety and comparability of treatments.
Methods
To overcome this deficiency accurate spatial and temporal temperature maps could be used. In this paper, the suitability of temperature calculation based on timeshifts of diagnostic ultrasound backscattered signals (echotimeshift) is investigated and associated uncertainties are estimated. Different analysis variations were used to calculate the timeshifts: discrete and continuous methods as well as different frames as a reference for temperature calculation (4 s before, 16 s before the frame of interest, base frame). A sigmoid function was fitted and used to calculate temperatures. Twodimensional temperature maps recorded during and after therapeutic ultrasound sonication were examined. All experiments were performed in agargraphite phantoms mimicking nonfatty tissue, with highintensity focused ultrasound being the source of heating.
Results
Continuous methods are more accurate than discrete ones, and uncertainties of calculated temperatures are in general lower, the earlier the reference frame was recorded. Depending on the purpose of the measurement, a compromise has to be made between the following: calculation accuracy (early reference frame), tolerance towards small movements (late reference frame), reproducing large temperature changes or cooling processes (reference frame at a certain point in time), speed of the algorithm (discrete (fast) vs. continuous (slower) shift calculation), and spatial accuracy (interval size for indexshift calculation). Within the range from 20 °C to 44 °C, uncertainties as low as 12.4 % are possible, being mainly due to medium properties.
Conclusions
Temperature measurements using the echotimeshift method might be useful for validation of treatment plan algorithms. This might also be a comparatively accurate, fast, and affordable method for laboratory and clinical quality assessment. Further research is necessary to improve filter algorithms and to extend this method to multiple foci and the usage of temperaturedependent tissue quantities. We used an analytical approach to investigate the uncertainties of temperature measurement. Different analysis variations are compared to determine temperature distribution and development over time.
Keywords
Ultrasound therapy HIFU Dosimetry Ultrasound thermometry Accuracy Uncertainty Echotimeshift method Quality assessment SafetyBackground
Ultrasound therapies
Ultrasound therapies have increasingly gained importance over the past few decades [1–4]. They aim to complement or replace standard therapies [5, 6], for instance, chemo or radiotherapy for cancer treatment [7, 8], thrombolysis of thrombosis [9, 10], or the removal of kidney and gall stones [11, 12]. Additionally, they have evolved for new therapies like viro, immuno, or gene therapy [13, 14], drug delivery through the bloodbrain barrier [15], or sonophoresis [16–18]. The safety of ultrasound therapies and the accuracy of their treatment plans and of their devices are crucial [19] for the wider acceptance and usage of the new therapies, for patient safety, and for administrative approval.
Safety assurance before the treatment
Considerable effort has been made towards objective safety criteria during treatment in the past few years, for example, based on the mechanical and thermal index [20], the maximum allowed temperature on the transducer surface [21], or realtime temperature control.
There is a lack of standards for implementation and execution of quality assurance processes before the therapy or a single treatment session. This is partially because there are not many possibilities to easily, but at the same time reliably, e.g. validate the algorithms of a treatment planning programme, calculated doses, and estimated side effects or to test therapy devices before treatment in medical facilities [22]. It is a major shortcoming for the reproducibility of treatment sessions and the comparability of clinical studies [23]. Tools needed for the tasks mentioned above must be easytohandle, reproducible, accurate, reliable, and inexpensive.
One possibility to test the functionality and the constancy of therapeutic ultrasound devices could be to sonicate a tissuelike phantom with standardized parameters and to monitor temperatures in a region of interest during a certain time (e.g. only during or during and after sonication). In this study, only nonfatty tissue phantoms are used to investigate the general feasibility of our methods.
Temperature measurement
Temperature measurement with diagnostic ultrasound is appropriate since it is relatively cheap, easily accessible in medical facilities, capable of use in real time, and allows two or threedimensional imaging. It thus has major advantages when being used for quality assessment compared to magnetic resonance or thermocouple measurements.
The idea of using ultrasound measurement for noninvasive thermometry has been promoted already 20 years ago [24, 25] and has been successively improved since then. Progress was made i.a. in the fields of twodimensional [26], threedimensional [27], and compound imaging [28], improved algorithms [29–33], and the estimation of limitations [34]. Furthermore, the applicability of the method in different media was described: homogeneous [35, 36] and multilayer [30] phantoms, and in tissue in vitro [37–40] and in vivo [41]. A comprehensive review can be found in [42].
What to expect in this paper
In our study, we investigated the suitability and uncertainties of a common method for temperature measurement: the echotimeshift method. Different analysis variations were used to calculate the timeshifts: discrete and continuous algorithms with the reference frame for each frame being either the frame 4 or 16 s before the frame of interest or the base frame. Moreover, we introduce an uncertainty calculation based on mathematical analysis and show specific possibilities to decrease uncertainties.
Methods
A. Temperature calculation with echotimeshift and echoindexshift
Here, c is the speed of sound in the medium at the initial frame at constant initial temperature 𝜗 _{0}, \(\frac {\partial }{\partial z}\) is the spatial derivation along axial depth z (along a scanline), and δt is the timeshift of the backscattered signal. c is assumed to change approximately linearly with temperature with the thermal coefficient of the speed of sound \(\beta = \frac {1}{c(\vartheta _{0})} \frac {\partial c(\vartheta)}{\partial \vartheta }\). α is the linear coefficient of thermal expansion. For our purposes, α and β are both considered to be independent of temperature and of axial depth, because a homogeneous phantom is used. Therefore, they can be combined to a value k (Eq. 2).
k is the proportionality factor between temperature change on the one hand, and relative change of both speed of sound and length on the other hand. In other words, it specifies the proportionality between temperature change and the shift of the signal, as shown later. k is larger than zero for nonfatty tissue and smaller than zero for fatty tissue. Thus, major problems will occur if the diagnostic ultrasound scanline passes fatty and nonfatty tissue and if a single value is presumed for k in these cases. In our study, we focus on a nonfatty tissue phantom. Prospects for the application in inhomogeneous phantoms, for instance, by using a depthdependent k and an iterative calculation procedure, must be investigated in further studies. k only depends on the particular material but neither on temperature nor on depth.
Here, δi accounts for the incremental indexshift, where incremental means additive along axial depth, and \(\frac {\partial }{\partial i}(\delta i)\) is the indexshift. Note that δt and δi are referred to as the incremental timeshift and the incremental indexshift, respectively, since they add up along a scanline. Their spatial derivation is then solely called time or indexshift.
For a given experiment, the absolute values of δi depend on the sampling frequency and interpolation of the original backscattered RF signal. The incremental indexshift is calculated with crosscorrelation as explained in more detail in the “D. Incremental indexshift calculation methods” section.
B. Temperature calculation with a sigmoid function fit
where i is the index of the data point of the preprocessed (i.e. interpolated or frequency filtered) RF signal, e is Euler’s number, and a to d are fit parameters. d accounts for the incremental indexshift that arose on that scanline axially before the beginning of the phantom, a is the maximum incremental indexshift difference to d, b is the slope at the turning point and therefore determines the value of maximum temperature, and c is the location of the turning point and therefore the location of maximum temperature.
C. Uncertainties of calculated temperatures
where x _{ n } are the independent variables (k, fit parameters a, b, and c), N is their number, and u _{ n } is the uncertainty of the associated variable x _{ n }.
Uncertainty of k
Uncertainty of temperature
The fit parameters are not independent on each other but their dependency is not known in detail. This correlation could be examined in the future for instance with Monte Carlo simulations, also taking interval sizes into account. Therefore, the uncertainties due to fit parameters add to a general fit uncertainty budget without squaring.
where the subscript “ref” means the sum over the associated reference frames. If the reference frame is another than the initial frame, uncertainties due to fitting add up and can only stay the same or grow larger. Otherwise, they can also decrease over time since they are independent of any calculation before. A decrease during cooling is possible, because at lower temperatures, less decorrelation due to spatial extension and fewer losses of interferences occur. This leads to less noise and a more suitable fit. The method of calculating uncertainty presented here assumes that incremental indexshifts can be measured continuously and filtered reasonably.
D. Incremental indexshift calculation methods
Discrete crosscorrelation
Continuous crosscorrelation
Here, ∢ denotes the fourquadrant inverse tangent (atan2 in many programming languages). Continuous phaseshifts are detectable within ±π and therefore the indexshifts δi _{ c } within ±δi _{cf} π. The intervals were truncated to [ −1.5, 1.5] because discrete incremental indexshifts were calculated before, so that δi _{ c } should be smaller than or equal to 1.0.
Extension to more dimensions
Applied calculation variations
Some properties can be modified when performing the calculations. These are the reference frame for and the level of continuity of the incremental indexshift calculation as well as the number of scanlines used for the crosscorrelation calculation. All three influence the actual value of the incremental indexshift and its uncertainty but should finally result in the same calculated temperatures.
For the comparison of incremental indexshift calculation variations, the same RF data was used. The RF signal was linearly interpolated from a 40 MHz sampling frequency to a 160 MHz equivalent to obtain smaller time steps between consecutive values in the backscattered ultrasound signal and therefore larger incremental indexshifts. This procedure is especially necessary for discrete methods.
E. Research design
Measurement setup
Procedure of temperature measurement
At the beginning, multiple frames at baseline temperature were recorded. Therapeutic ultrasound was run in discontinuous mode (3.85 s on, 150 ms off, duty cycle 96.25 %) for 80 s and actively heated the phantom. The ultrasound power was 20.8 W, as determined from radiation force balance measurements with the same sonication conditions. In every break (every 4 s), two diagnostic ultrasound frames were recorded. Afterwards, the phantom was allowed to cool down but it was not actively cooled. Diagnostic ultrasound frames were recorded for approximately another 2 min every 4 s. It turned out that the signal quality remains the same in both frames recorded during one HIFU break and therefore acquiring one frame in each break would be sufficient. A larger duty cycle could be implemented. Only one frame of each measuring break was used for the calculations. For temperature analysis, MATLAB was used (lsqcurvefit for sigmoid fitting with lower and upper boundaries and 4000 iterations per fit, nlparci to obtain the 95 % confidence intervals of the fit parameters with use of the residuals and Jacobian matrix).
Procedure of measurement of k
The value k was obtained by actively heating and cooling the water to welldefined temperatures (20 °C to 44 °C in steps of 8 °C) and waiting for the phantom to adjust its temperature. A smaller water bath, a single transducer (5 MHz, V309, 5.0/0.5 170680, Panametrics), and a pulser/receiver (USKey, Lecoeur Electronique) were used instead of the linear array that was used for HIFU measurements. This was due to the long exposure time to high water temperatures and to prevent damage to the linear array. To improve the measurement of k, a smaller phantom and a linear transducer could be used. Smaller temperature steps and multiple scanlines could then be analysed. Because of the long measurement time, the phantom was put into a standard freezer plastic bag to improve the temporal stability of the chemical composition, mainly the isopropanol content. The temperature measuring devices (thermocouple and resistance thermometer) were calibrated at PTB. Before the experiments were performed, manufacturers’ data on the size of the focus zone of the HIFU transducer were rechecked with hydrophone measurements in water.
F. Algorithm
General proceeding
 1.
Acquire experimental data including a baseline image at homogeneous temperature distribution
 2.Signal preprocessing

Subtract any existing offset

Interpolate RF signal to a 160 MHz equivalent sampling frequency (interp1)

Calculate Hilbert transform (hilbert)

Calculation of k: sort measurements from low (baseline) to high temperature

 3.
Calculate discrete incremental indexshifts (xcorr)
 4.
Calculate continuous incremental indexshift (xcorr, angle)
 5.
Filter incremental indexshifts along the axial direction for every scanline (medfilt1)
 6.
Linear fit (polyfit) for calculation of k or sigmoid fit (lsqcurvefit, nlparci for uncertainties) for temperature calculation
 7.
 8.
Add temperatures and uncertainties (without uncertainty due to k) to those from reference frames
 9.
In case of temperature calculation, add uncertainty of k
The index i was counted in relation to the front wall of the phantom to meet the requirements of the method, where the phantom has to expand away from the transducer. Halfoverlapping intervals were used for crosscorrelation calculations.
Filtering of indexshifts
Filtering of the calculated incremental indexshifts is essential for this method because outliers strongly influence the quality of both, linear fits for the calculation of k and sigmoid fits for temperature evaluation, including the uncertainties of the fit parameters. For temperature calculation, a median filter was applied twice on incremental indexshifts. First, five neighbours (eleven values) were taken into account, then one neighbour (three values).
A better but very challenging option would be a filter algorithm only based on the incremental indexshift difference between axially consecutive incremental indexshifts. If the difference is larger than an allowed value, the axially consecutive incremental indexshift is defined to be an outlier. The allowed difference is dependent on the therapeutic ultrasound power and absorption coefficient of the medium, the size of the focus zone, the degree of interpolation of the RF sampled data, the time within the heating or cooling process, and the reference frame used. We verified the general possibility of such a filter algorithm and obtained good results. However, we did not use it here because of its difficult implementation and up to now unmanageable usage in a clinic. It might nevertheless be useful for constancy testing purposes, where the abovementioned quantities are usually known.
In this study, neither spatial nor temporal filtering was applied to the final temperature maps since the filtering influences uncertainty of the methods in an incalculable way. For validation of treatment planning algorithms or laboratory dosimetry, temperature filtering is possible and recommended. As shown later, our results suggest that filtering of temperature maps could improve overall accuracy and reliability of single calculated temperature values.
Results and discussion
A. Final values and uncertainties of k
Properties of the phantom material are introduced into temperature calculation via the material dependent value k (Eq. 2).
Values and uncertainties of k
Final k  Standard deviation  Gaussian uncertainty  Uncertainty  

Discrete frametoframe  748.0 °C  76.6 °C (10.2 %)  35.4 °C (4.8 %)  112.0 °C (15.0 %) 
Continuous frametobase  747.9 °C  80.1 °C (10.7 %)  12.3 °C (1.7 %)  92.4 °C (12.4 %) 
k was much more dependent on temperature than previously reported in the literature (e.g. [26]) and also than the widespread application of the echotimeshift method would suggest. It is not constant within the considered temperature range as is assumed for the method. This is due to the nonlinear dependence of the speed of sound on the temperature which strongly influences the standard deviation. Therefore, the standard deviation of all calculated values of k plus the uncertainty due to the measurement (devices) and the calculation (algorithm including filtering and fitting) were combined to the final uncertainty of k. Ideally, both would be of the same size and due to the same physical and algorithmic effects.
The independence of k on temperature is inherent to nonfatty tissue [43] and therefore a necessary property of phantoms mimicking nonfatty tissue. Up to 44 °C, k can be considered temperature independent with the uncertainties shown in Table 1. Then, small temperature changes are overestimated and large ones underestimated. Depending on the usage of the method, this must be taken into account and possibly be corrected. The temperature dependence, especially at temperatures higher than 44 °C, is not a fundamental limitation for laboratory dosimetry because it is generally possible to implement smaller baseline temperatures. Continuing research should investigate as to what extent results are transferable to higher temperatures. The temperature dependence of k might be a strong limitation for clinical applications [34] because of the higher initial temperatures (36 °C).
To further decrease uncertainty due to the calculation of k, calculations could be optimized, for instance, by using better filter algorithms for outlying incremental indexshifts or by using transducer arrays that allow performing more comprehensive statistics on the calculation of k.
B. Temperature changes: comparison of different calculation methods
To compare different calculation variations and their uncertainties, in this part, we analyse the (1) incremental indexshifts, (2) temperaturetime curves, and (3) temperature maps.
1a. Discrete vs. continuous incremental indexshifts
During the first few seconds of heating, similar results are obtained for all methods for scanlines acquiring data from within the focus zone (Fig. 6). Discrete incremental indexshifts are sufficiently large for a sigmoid fit, which in turn is a good approximation for the data.
Significant differences between discrete and continuous calculations are found for the frametoframe methods (reference frames 4 s and 16 s before the examined frame) at all other times and in all other places: outside the focus zone even during the first seconds of heating (not shown), everywhere after the first seconds of heating (Fig. 7) and especially during cooling (Fig. 8). Discrete incremental indexshifts in these cases were very small and sometimes not even detectable (e.g. the detected indexshift is 0, the “true” indexshift is 0.4). For the frametobaseframe method, discrete calculation is possible for the whole heating and cooling process. The continuous indexshift serves as an additive correction and lowers uncertainty to some extent.
1b. Influence of the reference frame on incremental indexshifts
The reference frame determines the absolute value of incremental indexshifts. They are larger, the earlier the reference was recorded: smallest for the frametoframe method with the reference 4 s before (Figs. 6 a, b, 7 a, b, and 8 a, b), larger for the reference 16 s before (Figs. 6 c, d, 7 c, d, and 8 c, d), and largest for the reference being the base frame (Figs. 6 e, f, 7 e, f, and 8 e, f). The incremental indexshift for a single interval is positive if the examined correlation interval was heated up or negative if it cooled down. Indexshifts which are too big (larger than six wavelengths, 96 indices in our calculation), which could arise for frametobaseframe methods, do not occur.
During cooling and when using a frametoframe method (Fig. 8 a–d), decreasing and negative incremental indexshifts are detectable within the former heated and then cooling zone. Rising shifts occur at the former heated zone’s edges (Fig. 8 b, d) because of heat transfer out of the heated zone to the borders of the phantom. The sigmoid function does not approximate this curve shape appropriately. Neither decreasing nor negative shifts occur in frametobase methods since the phantom at every scan heated up compared to the base frame.
In general, more noise exists in the incremental indexshift curves for the frametobase methods compared to the continuous frametoframe methods. This is due to a loss of coherence because of spatial extension and therefore decorrelation at higher temperatures. It mostly disappears again if temperatures fall below the threshold where the specific decorrelation occurred if no permanent structural changes arose. Structural changes from one frame to a frame recorded shortly before are few and small enough to be overcome by the algorithm. In addition, noise could, in general, be caused by the mechanical movement of the phantom but this is irrelevant under laboratory conditions. A frametoframe algorithm will be more stable against small mechanical movements as they occur, e.g. in living objects. Smaller intervals for incremental indexshift calculations lead to some extent to a better spatial resolution, larger ones to less noise in the calculated data.
2. Temperature curves over time
Calculated temperaturetime curves of neighbouring scanlines (lateral) and intervals (axial direction, not shown) are in good agreement with each other for the continuous methods and the discrete frametobase method. This agreement of calculated curves with each other is an indication of accuracy of the method, next to absolute calculated uncertainties. This is especially true of the lateral neighbouring scanlines since temperatures for every scanline are calculated with an independent sigmoid fit and no filtering in lateral direction was done. The maximum calculated temperatures are alike within the uncertainty.
For the frametobase methods (Fig. 9 e, f), calculated and measured temperatures are consistent within the uncertainty. An offset between calculated and thermocouple temperatures during cooling might be due to the evaporation of alcohol from the phantom during storage and measurement and therefore an effectively smaller value k than used for temperature calculation based on k measurement. This is very likely and inversely means that significantly higher temperatures occur in the thermocouple measurements during the second half of heating. On the other hand, higher temperatures are rather underestimated due to the method used (see “Temperature calculation with echotimeshift and echoindexshift” section). Both effects could eliminate each other. The remaining underestimation of calculated temperatures compared to thermocouple temperatures is due to the viscous heating artefact [47, 48].
For continuous frametoframe methods (Fig. 9 b, d), measured and calculated temperatures are consistent during heating. The cooling process is neither reproduced satisfyingly enough to be used for laboratory dosimetry nor for the validation of algorithms for treatment planning.
Discrete frametoframe methods reveal major shortcomings. These are due to very small and even falsely nondetected incremental indexshifts, large uncertainties of the sigmoid fit, and error propagation due to adding up of the calculated temperature differences over all frames. The discrete frametoframe methods are hence only suitable if a reference frame of some temporal distance is chosen.
3. Twodimensional temperature maps
The heating process can, in general, be seen for all methods, the cooling process only for the frametobase and continuous methods. The latter do not reproduce the circular shape nor a uniform temperature distribution during cooling.
C. Uncertainty of temperature calculation
General thoughts on the uncertainty of temperature calculation
For the following uncertainty evaluation, the influence of k and the fit parameters (a to c) are considered. As mentioned before, spatial uncertainty is not taken into account separately since it strongly depends on measurement settings and can easily be added to the uncertainty budget. Depending on whether qualitative or quantitative temperature measurement is to be performed, uncertainties of different parameters have to be considered.
Qualitative temperature measurement is sufficient for determining the location of the focus, its shape and extension, or temperature distribution as percentiles. The uncertainty of k is irrelevant, whereas uncertainties of fit parameters and spatial resolution are crucial. The maximum acceptable spatial uncertainty depends on the size of the region to be treated and its closeness to risk organs. Percentiles of temperature distribution can be used to adjust the power of the therapeutic ultrasound device to meet treatment requirements (e.g. keep a minimum temperature in the focus zone for a certain time or stay below a maximum temperature at risk organs).
Quantitative temperature measurement resulting in the exact temperaturetime profiles is, among other things, necessary for thermal dose calculation, to validate treatment plans and the algorithms to calculate them, to estimate bioeffects in the target region, and to prevent side effects in organs at risk.
The maximum acceptable uncertainty for calculated temperatures depends on the desired or suppressed effects and their onset temperature intervals. Since these temperatures are different for various types of tissue and dependent on other variables, too, no assessment about the absolute uncertainties will be given.

The first is that the fit model must be correct.
The sigmoid fit model represents the incremental indexshift data reasonably for continuous methods, the discrete frametobase method, and the other discrete methods during the heating process within the heated zone. This is because the focus region is small and the transducer was not moved.
Problems arise for frametoframe methods during cooling: shifts fall in the heated zone due to cooling and rise on its borders due to thermal conduction out of the heated zone (Fig. 8 b, d; heating and cooling occur along a scanline). The sigmoid model is also insufficient for discrete frametoframe methods since neither discretization of the underlying data nor the discretization error is taken into account for the fit. If the maximum incremental indexshift is only a few indices, incorrect steep or flat sigmoid fits result. These uncertainties do not exist for the continuous methods since any shift is continuously calculable.
Depending on the type of measurement, other fit functions must be found, for instance, for validation of treatment plan algorithms with a moving transducer and multiple heated regions or if the focus is nonsymmetric. Calculations can then be applied as shown in this paper.

The second prerequisite is that variables have to be uncorrelated.
Parameters k, a, b, and c have to be uncorrelated to each other (Eq. 7) for uncertainty calculation. Only k is independent. It has not been possible to express the other dependencies in mathematical terms up to now (Monte Carlo methods could be a solution). We assume that the correlation between fit parameters a and b cannot be neglected for uncertainty calculations and that both reinforce each other because of the turning point symmetry of the sigmoid curve. Correlations with fit parameter d can probably be neglected because its uncertainty is very small. Any correlation between a and c is probably very small, whereas b and c could influence each other.

Filtered incremental indexshift data having to be representative is the third prerequisite.
That means that filter algorithms should change or delete outliers in a representative way. Filtering is easiest for continuous frametoframe methods, quite challenging for frametobase methods and largely impossible for discrete frametoframe methods, because of small and discrete shifts. It might be laborious or not possible to correct the following effects: loss of correlation because of small mechanical movements, decorrelation due to thermal expansion and the thermoacoustic lens effect, signal instability of the ultrasound transducer, and change of chemical composition over time (i.e. due to healing of tissue or chemical instability of phantoms).
Noise from the thermoacoustic lens’ effect could be reduced by methods like spatial compound imaging [28].
Uncertainty analysis of temperature calculation
The three prerequisites mentioned above are not fulfilled for the discrete variations. Therefore, only continuous temperature calculation variations are considered in the following.
Uncertainty of the sigmoid fit parameters was calculated with the residuals and Jacobian matrix (leads to 95 % intervals) and divided by two to obtain the standard deviation (67 % intervals). The uncertainty of k was determined in experiments performed by us.
Quantiles of uncertainty values
δ𝜗>0.1 °C  δ𝜗>6.0 °C  

Q _{50 %} /°C  Q _{5 %} /°C  Q _{95 %} /°C  Q _{50 %} /°C  Q _{5 %} /°C  Q _{95 %} /°C  
Heating and cooling  
Frametoframe 4 s before  3.92  0.20  377.30  24.77  1.57  389.81 
Frametoframe 16 s before  0.47  0.05  79.10  8.06  0.97  112.44 
Frametobaseframe  0.11  0.02  0.97  1.26  0.78  4.73 
Only heating  
Frametoframe 4 s before  0.40  0.07  2.31  2.59  1.20  6.20 
Frametoframe 16 s before  0.09  0.026  1.01  1.33  0.85  2.59 
Frametobaseframe  0.08  0.018  2.02  2.04  0.84  5.74 
Summary of uncertainty considerations
We investigated influences of the algorithmic evaluation (uncertainties of sigmoid fit parameters) as well as material properties (uncertainty of k which is due to the nonlinear dependence of the speed of sound on temperature and thermal expansion). Depending on the purpose of the uncertainty consideration, frametobase methods or frametoframe methods, where the reference frame has a sufficiently large temporal distance to the frame of interest, are suitable. If heating and cooling processes are to be monitored, the reference frame generally has to have a larger temporal distance.
The main contribution to uncertainty results from k. It could be significantly lowered if the nonlinearity were understood in more detail, if measurements were performed in other media or tissue where the relation of the speed of sound to temperature is more linear, or to a small amount if the measurement of k were enhanced. For tasks regarding human tissue, it is not possible to use different phantoms because the nonlinear dependence of the speed of sound on temperature is inherent to nonfatty tissue.
Measuring errors are thought to be negligible compared to the uncertainties of the value k. The ‘true’ uncertainty might be larger than what is shown in Table 2, because correlations between fit parameters a and b are not taken into account and are not thought to be negligible.
Conclusions
Discrete and continuous analysis variations of the echotimeshift method with reference frames both 4 s and 16 s before the frame of interest along with the base frame were investigated. Special attention was paid to uncertainty considerations of the method in general and concrete measurements on a tissuemimicking phantom.
The best results with regard to uncertainties, the reproducibility of the shape of the heated zone during heating and cooling, and in accordance with thermocouple measurements were obtained with the continuous frametobaseframe method. Continuous frametoframe methods provide good results during the heating process and additionally, a small robustness against small mechanical movements. Based on our findings, a minimum uncertainty of 12.4 % is intrinsic to the timeshift method within the range 20 °C to 44 °C due to tissue properties (nonlinear dependence of the speed of sound on temperature). Algorithmic evaluation also adds to the uncertainty budget but with further processing (e.g. spatial and temporal filtering of temperature maps) and additional physical and biological assumptions it could significantly be reduced.
No uncertainty statement is possible for a single calculated temperature value. An uncertainty calculation was only reasonable for the continuous methods because of incalculable influences on uncertainties of the discrete methods. As far as possible, calculated temperatures were compared to and similar to thermocouple measurements.
Continuous timeshift calculations should be used not only for the quality assessment of therapeutic ultrasound devices but also for realtime therapy control because they outstandingly improve measurement compared to discrete methods. Discrete frametoframe methods have major shortcomings for laboratory dosimetry. The most significant ones are small incremental timeshifts from one frame to the reference frame and the additivity of uncertainties.
Methods using a reference frame a few seconds before the frame of interest are more robust against small mechanical movements and correlation loss due to thermal expansion and the change of interferences. Therefore, they seem most promising for various tasks.
Whether or not the achieved uncertainties are sufficient for laboratory dosimetry, validation of treatment plan algorithms, or decisions whether or not a therapeutic ultrasound device can be used for therapy, must be determined in further investigations. It is dependent on the biological responses of human tissue on different temperatures. To verify complex treatment plans, this method could be extended by altering or replacing the sigmoid fit function.
Abbreviations
 HIFU:

Highintensity focused ultrasound
 RF signal:

Backscattered highfrequency signal
Declarations
Acknowledgements
Not applicable.
Funding
This work was funded through the EMRP project HLT03 DUTy–Dosimetry for Ultrasound Therapy. The EMRP is jointly funded by the EMRP participating countries within EURAMET and by the European Union. The funding body did not influence the design of the study, collection, analysis, or interpretation of the data.
Availability of data and materials
The datasets during and/or analysed during the current study are available from the corresponding author on reasonable request.
Authors’ contributions
All authors designed and redesigned the experiments and contributed their experience and ideas to the evaluations as well as to the uncertainty calculation methods. TF and OG conducted all experiments and performed the analysis of the frametoframe method. Frametobase frame methods were investigated by TF. TF compiled this paper and all authors reviewed and revised it.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Tyshlek D, Aubry JF, ter Haar G, Hananel A, Foley J, Eames M, Kassell N, Simonin HH. Focused ultrasound development and clinical adoption: 2013 update on the growth of the field. J Ther Ultrasound. 2014; 2(1):1–7.View ArticleGoogle Scholar
 Mason TJ. Therapeutic ultrasound an overview. Ultrason Sonochem. 2011; 18(4):847–52.PubMedView ArticleGoogle Scholar
 Kennedy JE. High intensity focused ultrasound: surgery of the future?Br J Radiol. 2003; 76(909):590–9.PubMedView ArticleGoogle Scholar
 Paliwal S, Mitragotri S. Therapeutic opportunities in biological responses of ultrasound. Ultrasonics. 2008; 48(4):271–8.PubMedView ArticleGoogle Scholar
 ter Haar G. Therapeutic ultrasound. Eur J Ultrasound. 1999; 9(1):3–9.PubMedView ArticleGoogle Scholar
 Jenne JW, Preusser T, Günther M. Highintensity focused ultrasound: principles, therapy guidance, simulations and applications. Zeitschrift für Medizinische Physik. 2012; 22(4):311–22.PubMedView ArticleGoogle Scholar
 AlBataineh O, Jenne J, Huber P. Clinical and future applications of high intensity focused ultrasound in cancer. Cancer Treat Rev. 2012; 38(5):346–53.PubMedView ArticleGoogle Scholar
 Shehata IA. Treatment with high intensity focused ultrasound: secrets revealed. Eur J Radiol. 2012; 81(3):534–41.PubMedView ArticleGoogle Scholar
 Francis CW. Ultrasoundenhanced thrombolysis. Echocardiography. 2001; 18(3):239–46.PubMedView ArticleGoogle Scholar
 Wright C, Hynynen K, Goertz D. In vitro and in vivo highintensity focused ultrasound thrombolysis. Investig Radiol. 2012; 47(4):217–25.View ArticleGoogle Scholar
 Yoshizawa S, Ikeda T, Ito A, Ota R, Takagi S, Matsumoto Y. High intensity focused ultrasound lithotripsy with cavitating microbubbles. Med Biol Eng Comput. 2009; 47(8):851–60.PubMedView ArticleGoogle Scholar
 Duryea AP, Hall TL, Maxwell AD, Xu Z, Cain CA, Roberts WW. Histotripsy erosion of model urinary calculi. J Endourol. 2011; 25(2):341–4.PubMedPubMed CentralView ArticleGoogle Scholar
 Sorace AG, Warram JM, Mahoney M, Zinn KR, Hoyt K. Enhancement of adenovirus delivery after ultrasoundstimulated therapy in a cancer model. Ultrasound Med Biol. 2013; 39(12):2374–81.PubMedPubMed CentralView ArticleGoogle Scholar
 Newman CMH, Bettinger T. Gene therapy progress and prospects: ultrasound for gene transfer. Gene Therapy. 2007; 14(6):465–75.PubMedView ArticleGoogle Scholar
 Aryal M, Arvanitis CD, Alexander PM, McDannold N. Ultrasoundmediated bloodbrain barrier disruption for targeted drug delivery in the central nervous system. Adv Drug Deliv Rev. 2014; 72:94–109.PubMedView ArticleGoogle Scholar
 BorSengShu E, De Carvalho Nogueira R, Figueiredo EG, Evaristo EF, Conforto AB, Teixeira MJ. Sonothrombolysis for acute ischemic stroke: a systematic review of randomized controlled trials. Neurosurg Focus. 2012; 32(1):1–6.View ArticleGoogle Scholar
 Ogura M, Paliwal S, Mitragotri S. Lowfrequency sonophoresis: current status and future prospects. Adv Drug Deliv Rev. 2008; 60(10):1218–23.PubMedView ArticleGoogle Scholar
 Mitragotri S, Kost J. Lowfrequency sonophoresis. Adv Drug Deliv Rev. 2004; 56(5):589–601.PubMedView ArticleGoogle Scholar
 Miller DL, Smith NB, Bailey MR, Czarnota GJ, Hynynen K, Makin IRS. Overview of therapeutic ultrasound applications and safety considerations. J Ultrasound Med. 2012; 31(4):623–34.PubMedPubMed CentralGoogle Scholar
 IEC 62359:2010. Ultrasonicsfield characterization—test methods for the determination of thermal and mechanical indices related to medical diagnostic ultrasonic fields. Edition 2.0, International Electrotechnical Commission. Geneva; 2010.
 IEC 60601237:2007. Particular requirements for the basic safety and essential performance of ultrasonic medical diagnostic and monitoring equipment. International Electrotechnical Commission; 2007.
 Shaw A, Martin E, Haller J, ter Haar G. Equipment, measurement and dose—a survey for therapeutic ultrasound. J Ther Ultrasound. 2016; 4(7):1–9.Google Scholar
 Shaw A, ter Haar G. Telling it like it is. J Ther Ultrasound. 2013; 1(4):1–2.Google Scholar
 Seip R, VanBaren P, Cain CA, Ebbini ES. Noninvasive realtime multipoint temperature control for ultrasound phased array treatments. IEEE Trans Ultrason Ferroelectr Freq Control. 1996; 43(6):1063–73.View ArticleGoogle Scholar
 MaassMoreno R, Damianou CA. Noninvasive temperature estimation in tissue via ultrasound echoshifts. part 1. analytical model. J Acoust Soc Am. 1996; 100(4):2514–21.PubMedView ArticleGoogle Scholar
 Simon C, VanBaren P, Ebbini ES. Twodimensional temperature estimation using diagnostic ultrasound. IEEE Trans Ultrason Ferroelectr Freq Control. 1998; 45(4):1088–99.PubMedView ArticleGoogle Scholar
 Anand A, Kaczkowski P. J. Noninvasive measurement of local thermal diffusivity using backscattered ultrasound and focused ultrasound heating. Ultrasound Med Biol. 2008; 34(9):1449–1464.PubMedPubMed CentralView ArticleGoogle Scholar
 Pernot M, Tanter M, Bercoff J, Waters KR, Fink M. Temperature estimation using ultrasonic spatial compound imaging. IEEE Trans Ultrason Ferroelectr Freq Control. 2004; 51(5):606–15.PubMedView ArticleGoogle Scholar
 Liu D, Ebbini ES. Realtime 2D temperature imaging using ultrasound. IEEE Trans Biomed Eng. 2010; 57(1):12–6.PubMedView ArticleGoogle Scholar
 Teixeira CA, Pereira WCA, Ruano AE, Graça Ruano M. On the possibility of noninvasive multilayer temperature estimation using softcomputing methods. Ultrasonics. 2010; 50(1):32–43.PubMedView ArticleGoogle Scholar
 Huang SM, Liu HL, Li ML. Improved temperature imaging of focused ultrasound thermal therapy using a sigmoid model based crosscorrelation algorithm. In: 2012 IEEE International Ultrasonics Symposium (IUS). IEEE: 2012. p. 2766–8.
 Chen BT, Shieh J, Huang CW, Chen WS, Chen SR, Chen CS. Ultrasound thermal mapping based on a hybrid method combining physical and statistical models. Ultrasound Med Biol. 2014; 40(1):115–29.PubMedView ArticleGoogle Scholar
 Dutta D, Mahmoud A, Leers S, Kim K, Motion artifact reduction in ultrasound based thermal strain imaging of atherosclerotic plaques using timeseries analysis. IEEE Trans Ultrason Ferroelectr Freq Control. 2013; 60(8):1660–1668.PubMedPubMed CentralView ArticleGoogle Scholar
 Miller NR, Bamber JC, Meany PM. Fundamental limitations of noninvasive temperature imaging by means of ultrasound echo strain imaging. Ultrasound Med Biol. 2002; 28(10):1319–33.PubMedView ArticleGoogle Scholar
 Daniels MJ, Varghese T, Madsen EL, Zagzebski JA. Noninvasive ultrasoundbased temperature imaging for monitoring radiofrequency heatingphantom results. Phys Med Biol. 2007; 52(16):4827–43.PubMedView ArticleGoogle Scholar
 Bazán I, Vazquez M, Ramos A, Vera A, Leija L. A performance analysis of echographic ultrasonic techniques for noninvasive temperature estimation in hyperthermia range using phantoms with scatterers. Ultrasonics. 2009; 49(3):358–76.PubMedView ArticleGoogle Scholar
 MaassMoreno R, Damianou CA, Sanghvi NT. Noninvasive temperature estimation in tissue via ultrasound echoshifts. part 2. in vitro study. J Acoust Soc Am. 1996; 100(4):2522–30.PubMedView ArticleGoogle Scholar
 Miller NR, Bamber JC, ter Haar GR. Imaging of temperatureinduced echo strain: preliminary in vitro study to assess feasibility for guiding focused ultrasound surgery. Ultrasound Med Biol. 2004; 30(3):345–56.PubMedView ArticleGoogle Scholar
 Souchon R, Bouchoux G, Maciejko E, Lafon C, Cathignol D, Bertrand M, Chapelon JY. Monitoring the formation of thermal lesions with heatinduced echostrain imaging: a feasibility study. Ultrasound Med Biol. 2005; 31(2):251–9.PubMedView ArticleGoogle Scholar
 Civale J, Rivens I, ter Haar G, Morris H, Coussios C, Friend P, Bamber J. Calibration of ultrasound backscatter temperature imaging for highintensity focused ultrasound treatment planning. Ultrasound Med Biol. 2013; 39(9):1596–612.PubMedView ArticleGoogle Scholar
 Varghese T, Zagzebski JA, Chen Q, Techavipoo U, Frank G, Johnson C, Wright A, Lee FT. Ultrasound monitoring of temperature change during radiofrequency ablation: preliminary invivo results. Ultrasound Med Biol. 2002; 28(3):321–9.PubMedView ArticleGoogle Scholar
 Arthur RM, Straube WL, Trobaugh JW, Moros EG. Noninvasive estimation of hyperthermia temperatures with ultrasound. Int J Hyperthermia. 2005; 21(6):589–600.PubMedView ArticleGoogle Scholar
 Bamber JC, Hill CR. Ultrasonic attenuation and propagation speed in mammalian tissues as a function of temperature. Ultrasound Med Biol. 1979; 5(2):149–57.PubMedView ArticleGoogle Scholar
 JCGM 100:2008. Evaluation of measurement data—guide to the expression of uncertainty in measurement. Joint Committee for Guides in Metrology. 2008.
 Loupas T, Powers JT, Gill RW. An axial velocity estimator for ultrasound blood flow imaging, based on a full evaluation of the Doppler equation by means of a twodimensional autocorrelation approach. IEEE Trans Ultrason Ferroelectr Freq Control. 1995; 42(4):672–88.View ArticleGoogle Scholar
 Madsen EL, Zagzebski JA, Banjavie RA, Jutila RE. Tissue mimicking materials for ultrasound phantoms. Med Phys. 1978; 5(5):391–4.PubMedView ArticleGoogle Scholar
 Waterman FM, Leeper JB. Temperature artifacts produced by thermocouples used in conjunction with 1 and 3 MHz ultrasound. Int J Hyperthermia. 1990; 6(2):383–99.PubMedView ArticleGoogle Scholar
 Hynynen K, Martin CJ, Watmough DJ, Mallard JR. Errors in temperature measurement by thermocouple probes during ultrasound induced hyperthermia. Br J Radiol. 1983; 56:969–70.PubMedView ArticleGoogle Scholar