Temperature imaging with speed of ultrasonic transmission tomography for medical treatment control: A physical model-based method
Chu Zhe-Qia), Yuan Jie†a), Pinter Stephen Z.b), Kripfgans Oliver D.b), Wang Xue-Dingb), Carson Paul L.‡b), Liu Xiao-Junc)
School of Electronic Science and Engineering, Nanjing University, Nanjing 210093, China
Department of Radiology, University of Michigan, Ann Arbor, Michigan 48109, USA
School of Physics, Nanjing University, Nanjing 210093, China

Corresponding author. E-mail: yuanjie@nju.edu.cn

Corresponding author. E-mail: pcarson@umich.edu

*Project supported in part by DoD/BCRP Idea Award, BC095397P1, the National Natural Science Foundation of China (Grant No. 61201425), the Natural Science Foundation of Jiangsu Province, China (Grant No. BK20131280), the Priority Academic Program Development of Jiangsu Provincial Higher Education Institutions, China, and the National Institutes of Health (NIH) of United States (Grant Nos. R01AR060350, R01CA91713, and R01AR055179).

Abstract

Hyperthermia is a promising method to enhance chemo and radiation therapy of breast cancer. In the process of hyperthermia, temperature monitoring is of great importance to assure the effectiveness of treatment. The transmission speed of ultrasound in biomedical tissue changes with temperature. However, when mapping the speed of sound directly to temperature in each pixel as desired for using all speeds of ultrasound data, temperature bipolar edge enhancement artifacts occur near the boundary of two tissues with different speeds of ultrasound. After the analysis of the reasons for causing these artifacts, an optimized method is introduced to rebuild the temperature field image by using the continuity constraint as the judgment criterion. The significant smoothness of the rebuilding image in the transitional area shows that our proposed method can build a more precise temperature image for controlling the medical thermal treatment.

PACS: 43.35.Wa; 43.80.Sh; 81.70.Cv
Keyword: temperature imaging; speed of sound; continuous constrain; thermal therapy
1. Introduction

Surgery, radiation therapy and chemotherapy are widely known as the most common types of cancer treatment and have been studied for a long time. Cancer surgery, [1, 2] an operation to repair or remove part of the body to diagnose or treat cancer, remains the foundation of cancer treatment. Radiation therapy and chemotherapy are methods of killing cancer cells respectively by high-energy radiation and drugs. Many studies[36] have proved the significant effects of these two methods.

With the development of electromagnetic[7] and ultrasonic deep heating, hyperthermia becomes an emerging method for breast tumor treatment. Plenty of researches have proved that the combination of hyperthermia and other methods, especially radiotherapy and chemotherapy, often contributes to better treatment results.[811]

In the process of hyperthermia, temperature monitoring is very important to assure safe treatment.[12] Target temperature ranges are typically quite tight, [13] e.g., 42  ° C or 43  ° C ± 1  ° C, to achieve the highest therapeutic enhancement. Insufficient temperature control not only reduces the effectiveness of hyperthermia but also can cause pain and damage.[14] Temperature measurement methods can be classified as invasive ones and noninvasive ones. Invasive temperature measurement methods are those that heat-sensitive sensors are embedded into tissues. They offer many advantages, such as accurate results, easy operations, and low cost. However, sparse spatial sampling provides temperature information of limited spots instead of the whole tissue[15] and there is some pain and the risk of infection. Non-invasive methods can avoid these problems by not disrupting the skin and other organs. Commonly employed or studied non-invasive measurement methods employ magnetic resonance (MR), [16] microwave, [17, 18] and ultrasound[19, 20] techniques. Three-dimensional (3D) mapping of temperature with the MR method based on the relaxation time, the diffusion coefficient, [21] or proton resonance frequency (PRF)[22] of tissue water has been well described previously. However, these methods only consider the relative temperature changes by comparing a reference image with images acquired under identical experimental conditions. An alternative approach to the calculation of absolute temperature is MR spectroscopic imaging by measuring the water resonance shift.[23, 24] Recently, a real-time MR temperature monitoring method based on the relaxation time of tissue is proposed[25] by demonstrating the capability of real-time simultaneous temperature mapping in aqueous tissue and fat during HIFU ablation. However, the entire process of MR is quite expensive for a cycle of multiple treatments. Techniques based on microwaves to measure temperature can be easily affected by the environment and can only penetrate a limited depth.[26] Ultrasound is a non-ionizing, convenient and inexpensive modality with relatively simple signal processing requirements, which makes it an attractive method to be used for temperature monitoring. Methods of using ultrasound as a temperature measurement include those that measure echo-shift changes in tissue thermal expansion and the speed of sound, [27, 28] those that measure the acoustic attenuation coefficient, [29, 30] and those that exploit the change in backscattered energy (CBE)[31] from tissue in homogeneities. However, the prior knowledge of the speed of sound (SOS) and thermal expansion coefficient are necessary if the echo-shift method is employed.[32]

The structure of heat dissipating a blood vessel[33, 34] in a breast tumor is rather complex, making high resolution temperature imaging critical for uniform thermal therapy or treatment enhancement. Local changes in the speed of sound, [35, 36] attenuation of ultrasound, [37] and thermal conductivity[38] are caused by changes in temperature. Since the thermal coefficient of the speed of ultrasound varies with tissue, the speed of ultrasound in tissue can be used to calculate the temperature field. Temperature distribution in biomedical tissue should meet the following equation:

where q is the thermal density, T is the temperature, and λ is the thermal conductivity.

According to Eq.  (1), temperature should change monotonically and smoothly. However, when we directly map the value of the speed of sound to obtain the temperature, a bipolar edge enhancement artifact will occur around the boundary of two different tissues. Unlike the ripple artifact, which is often observed around the heated region in the temperature estimations caused by the thermal– acoustic lens effect, [39, 40] the bipolar edge enhancement artifact results from the abrupt change in thermal coefficient of the speed of sound and the finite resolution of the imaging system.

In this paper, we analyze the causes of these artifacts and propose an optimization method to rebuild a correct temperature field image. The rest of this paper is organized as follows. In Section 2, both the direct method and the proposed method of optimally rebuilding the temperature image are presented. In Section 3, the experimental analysis and the quantitative analysis are given for subjective evaluation. Finally, some conclusions and a discussion are given in Section 4.

2. Principle of rebuilding the temperature image
2.1. Principle of heat transfer

Usually, the temperature, T, of an object can be described as a function of spatial coordinates x, y, z, and time t as follows:

Since the variation of temperature is a gradual process, we can neglect the time parameter to simplify the partial differential equation of heat conduction for temperature fields into

where ϕ and λ respectively represent the heat generated by the heat source per unit volume per unit time and the thermal conductivity.

2.2. Calculation of temperature field

2.2.1. Traditional method

By threshold segmentation of the speed of ultrasound image, tissue and its transitional area are separated by the following rule:

where v(x, y) is the speed of ultrasound at (x, y), v1 and v2 are average ultrasound transmission speeds in tissue1 and tissue2 respectively, Sm (m = 1, 2) are the areas of tissue1 and tissue2 and we suppose v1 < v2.

With the profiles between temperature and speed of ultrasound in tissue1 and tissue2, f1(v) and f2(v), temperatures in tissue1 and tissue2 can be acquired through

The SOS image mixes the speeds of sound of the two tissues at their boundaries due to the finite resolution of the imaging system, given by the system’ s line spread function (LSF), that is,

as the speed of sound in mixtures of homogeneous tissues follows a simple linear mixture model. Equation (7) gives a linear combination of the two different speeds at a sharp vertical boundary with a rectangular LSF. Thus, the measured speed v′ (x) in the transitional area would be

where Δ x = x2x1, Δ y = y2y1, and (xm, ym) (m = 1, 2) are the points in the boundary of the transitional area and tissuem.

The temperature at a sharp vertical boundary ideally would be given by a similar equation

We use the data of the speed of sound in the transitional area with a pixel and conduct the pixel calculation from the speed of sound. The temperature in the mixing area is expressed as

This provides distracting edge enhancement when the transitional area is not accurately segmented by Eq.  (4). To illustrate this, suppose that a homogeneous tissue with low temperature is surrounded by a high temperature tissue. When tissues and their transitional areas are well distinguished and equations  (6) and (9) are correctly applied, the temperature curve for the model is represented by the blue line in Fig.  1. However, if the boundary between two tissues is an unknown mixture of the two, equation  (10) must be employed and even a sharp boundary will have the edge artifacts of the red line when the transitional area is segmented into tissues.

Fig.  1. Illustration of abnormal bipolar edge enhancement artifacts.

2.2.2. Proposed method with optimization

To allow for thin layers of mixed tissue and using Eq.  (10), or one with a different LSF, continuity of temperature distribution is used as a constraint to rebuild the temperature image. After rough segmentation using Eq.  (4), the boundary of the transitional area is determined by the following equation:

where dvm (m = 1, 2) is the average of the gradient of the speed of ultrasound in the identified tissue1 or tissue2:

The transitional area is denoted as D, and its boundary as L. Generally, we assume that the whole transitional area is segmented into two parts as Dm (m = 1, 2) by the boundary. Then thermal conductivity λ in each part is a constant. In addition, there is no inner heat source in each transitional area. Thus, equation  (3) in each Dm is simplified into

Since the two transitional parts share the boundary, we can obtain the heat conduction equation of the whole transitional area, which is expressed as

Temperature at L as the boundary condition is calculated by Eqs.  (11) and (12) and Eqs.  (5) and (6) as follows:

The finite element method[41] is a numerical technique for finding approximate solutions to boundary value problems for partial differential equations. We divide the whole area of interest into similar parts called finite elements and use this method from the calculus of variations to solve the problem. In order to solve Eq.  (14) with the boundary condition described in Eq.  (15), the transitional area D is divided into m non-overlapped triangles e1, e2, … , em. For all m triangles, there are n non-overlapped vertices d1, d2, … , dm.

For triangle ei (i = 1, 2, … , m), coordinates of its three vertices are denoted as (xi1, yi1), (xi2, yi2), and (xi3, yi3). Equation (15) is used to calculate the coefficient matrix Kei:

where r, s, and t are the index numbers of the three vertices of triangle ei, and

The n × n matrix K is built by elements kpg from Kei:

Temperature of all vertices is acquired by solving:

where T = (T1, T2, … , Tn)′ is the vector of the temperature of all vertices. The temperature of vertices on boundary L is given by Eq.  (15) as a known boundary condition.

2.2.3. Reconstruction quality evaluation

While rebuilding the temperature field image based on Eq.  (14), we take the absolute difference between the temperature of one pixel and the mean temperature of its adjacent four pixels as a criterion to evaluate the smoothness of the rebuilding image, which is expressed as

Bipolar edge enhancement artifacts are the abnormal results around the boundary. Mean variation of temperature of adjacent pixels in the transitional area, which is defined in Eq.  (21), is calculated as a criterion to evaluate the image quality of the reconstructed temperature field, and is given as

where M is the total number of pixels in the transitional area.

3. Experiment and analysis

In our experiment, we use a gelatin cylinder with a central flow type to mimic the presence of benign and cancerous masses embedded in glandular tissue as a representative x-ray CT scan of the breast phantom shown in the literature.[42] 15% gelatin and cool water in the tube are respectively used to mimic tumor and fat. To record the transmitted and reflected ultrasound energy, a clinical ultrasound ring array scanner for breast cancer diagnosis, termed computed ultrasound risk evaluation (CURE), has been designed and built at the Karmanos Cancer Institute (KCI), Detroit, MI, USA.[43, 44] During scanning, all transducer elements sequentially emit a fan beam of ultrasound signals toward the opposite side of the ring. The scattered (transmission) and backscattered (reflection) ultrasound signals are subsequently recorded by all elements. Data of all wave fields (including both transmission and reflection) are then used to reconstruct the images of acoustic properties.[42]

Reflection mode and speed of sound images (e.g., Figs.  2(c) and 2(b) respectively) of a gelatin cylinder with a central flow tube (Fig.  2(a)), immersed in water are obtained with a prototype SoftVue® breast imaging system (Delphinus Medical Technologies, Inc., Plymouth, MI, USA). The different sections of the image are segmented from reflection images.

Fig.  2. (a) Cylindrical flow phantom of cooled gelatin in a water bath heated to 38  ° C; (b) speed of sound image, with 15% gelatin (white cylinder with centered dark hole from cooled water in a thin walled latex tube) and with a very attenuating thick walled silicone flow return tube outside the gelatin at 8 o’ clock; (c) reflection image; (d) prior relationship of speed of ultrasound with temperature for gelatin and water; (e) transitional area rough segmentation; (f) optimized transitional area segmentation.

Figure  2(e) shows the original segmented transitional area obtained using Eqs.  (4)– (6) and Eq.  (10). After applying our proposed method (using Eqs.  (11) and (12)), the optimized transitional area is segmented as shown in Fig.  2(f). It is the inaccurate traditional segmentation that is responsible for the nonphysical bipolar edge enhancement, according to our previous analysis. By comparison, it is obvious that the transitional area is optimized with our proposed method.

Curves of sound speeds of water and 15% gelatin were interpolated in Refs.  [45] and [46], respectively, to make the temperature image and profiles. Figure  3(a) shows the un-optimized three-dimensional (3D) temperature image of the whole gelatin, water and their boundary, obtained by using the traditional direct mapping method. The optimized 3D temperature image of the whole area is shown in Fig.  3(b), calculated by Comsol® Multi-Physics 3.5a with 25 repeated simulations through the finite element method. Figure  3(c) displays the comparison of temperature profile between with and without optimization.

From the 3D temperature image calculated directly by the speed of sound image and a profile of it, we can see that in the transitional area between gelatin and water there are nonphysical bipolar edge enhancement artifacts and the temperature changes sharply. Comparing Fig.  3(a) with Fig.  3(b) and also comparing the blue temperature profile with the red one in Fig.  3(c), after optimization, the edge enhancement artifacts are effectively alleviated and the temperature in the transitional area is much smoother.

Fig.  3. Experiment results comparison. (a) 3D temperature field built by the traditional method without optimization; (b) 3D temperature field built by the proposed method with optimization; (c) comparison of rebuilt temperature profiles from different angles with and without optimization.

Absolute differences of temperature ∇ T (x, y) from the transitional area with and without optimization in all rebuilt temperature images rows are shown in Fig.  4. The value of ∇ T (x, y) without optimization is much larger than that with optimization, which proves that our proposed method can rebuild a better temperature image which fits Eq.  (12) very well. Mean variation of temperature Δ T in the transitional area is shown in Fig.  5. With our proposed optimization method, the mean variation of temperature decreases tremendously, which means that the temperature in the transitional area is much smoother.

Fig.  4. Comparison of the absolute difference between with and without optimization.

Fig.  5. Comparison of mean variation between with and without optimization.

4. Conclusions and discussion

The structures of breast tissues and in particular breast cancers are rather complex, producing quite non-uninform heating patterns in microwave and ultrasound hyperthermia adjunctive to repeated chemo and radiation therapy sessions. This makes low cost, high speed, high fidelity temperature imaging essential for routine achievement of the benefits of hyperthermia in a cost efficient health care system. The strong bipolar edge enhancement artifacts in the temperature imaging derived from pixel by pixel mapping from the speed of sound will make this task much more difficult. The smoothing achieved with our proposed optimization method is helpful and adaptive. Experimental results show that after optimizing, bipolar edge enhancement artifacts are effectively reduced and the temperature in the transitional area is much smoother. Temperature imaging with transmission tomography for localized feedback of hyperthermia seems quite feasible where stable ultrasound transmission imaging can be achieved.

Acknowledgment

We thank N. Duric of the Karmanos Cancer Institute, and M. L. Scarpelli of the University of Michigan for contributions to conceptualization and experiments, respectively.

Reference
1 Nelson H, Petrelli N, Carlin A, Couture J, Fleshman J, Guillem J, Miedema B, Ota D and Sargent D 2001 J. Natl. Cancer Inst. 93 583 DOI:10.1093/jnci/93.8.583 [Cited within:1]
2 Merkow R P, Bentrem D J, Cohen M E, Paruch J L, Weber S M, Ko C Y and Bilimoria K Y 2013 J. Am. Coll. Surgeons 217 685 DOI:10.1016/j.jamcollsurg.2013.05.015 [Cited within:1]
3 Birgisson H, Påhlman L, Gunnarsson U and Glimelius B 2007 Acta Oncol. 46 504 DOI:10.1080/02841860701348670 [Cited within:1]
4 Okuma K, Yamashita H, Niibe Y, Hayakawa K and Nakagawa K 2011 J. Med. Case. Rep. 5 111 DOI:10.1186/1752-1947-5-111 [Cited within:1]
5 Bosset J F, Calais G, Mineur L, Maingon P, Radosevic J L, Daban A, Bardet E, Beny A, Briffaux A and Collette L 2005 J. Clin. Oncol. 23 5620 DOI:10.1200/JCO.2005.02.113 [Cited within:1]
6 Rodin G and Ahles T A 2012 J. Clin. Oncol. 30 3568 DOI:10.1200/JCO.2012.43.5776 [Cited within:1]
7 Krauss A, Nill S, Tacke M and Oelfke U 2011 Int. J. Radiat. Oncol. Biol. Phys. 79 579 DOI:10.1016/j.ijrobp.2010.03.043 [Cited within:1]
8 Zagar T M, Oleson J R, Vujaskovic Z, Dewhirst M W, Craciunescu O I, Blackwell K L, Prosnitz L R and Jones E L 2010 Int. J. Hyperther. 26 612 DOI:10.3109/02656736.2010.487194 [Cited within:1]
9 Huilgol N G, Gupta S and Sridhar C R 2010 J. Can. Res. Ther. 6 492 DOI:10.4103/0973-1482.77101 [Cited within:1]
10 Cherukuri P, Glazer E S and Curley S A 2010 Adv. Drug. Deliv. Rev. 62 339 DOI:10.1016/j.addr.2009.11.006 [Cited within:1]
11 Issels R D, Lindner L H, Verweij J, Wust P, Reichardt P, Schem B C, Abdel R S, Daugaard S, Salat C and Wendtner C M 2010 Lancet. Oncol. 11 561 DOI:10.1016/S1470-2045(10)70071-1 [Cited within:1]
12 Ma G and Jiang G 20103rd International Conference on IEEE Biomedical Engineering and Informatics (BMEI) October 16–18, 2010 Yantai, China 1357 DOI:10.1109/BMEI.2010.5639382 [Cited within:1]
13 Hildebrand t B, Wust P, Ahlers O, Dieing A, Sreenivasa G, Kerner T, Felix R and Riess H 2002 Crit. Rev. Oncol. Hematol. 43 33 DOI:10.1016/S1040-8428(01)00179-2 [Cited within:1]
14 Dewhirst M W, Viglianti B L, Lora-Michiels M, Hanson M and Hoopes P J 2003 Int. J. Hyperther. 19 267 DOI:10.1080/0265673031000119006 [Cited within:1]
15 Yang C L, Zhu H, Wu S C, Bai Y P and Gao H J 2010 J. Ultras. Med. 29 1787 [Cited within:1]
16 Dickinson R J, Hall A S, Hind A J and Young I R 1985 J. Comput. Assist. Tomogr. 10 468 [Cited within:1]
17 Stauffer P R, Rodriques D B, Salahi S, Topsakal E, Oliveira T R, Prakash A, D’Isidoro F, Reudink D, Snow B W and Maccarini P F 2013 SPIE BiOS. International Society for Optics and Photonics85840R DOI:10.1117/12.2003976 [Cited within:1]
18 Maruyma K, Mizushina S, Sugiura T, Van L G M J, Hand J W, Marrocco G, Bardati F, Edwards A D, Azzopardi D and Land D 2000 IEEE Trans. Microw. Theory Tech. 48 2141 DOI:10.1109/22.884206 [Cited within:1]
19 Wadley H N G, Norton S J, Mauer F, Droney B, Ash E A and Sayers C M 1986 Philos. Trans. R. Soc. Lond. A 320 341 DOI:10.1098/rsta.1986.0123 [Cited within:1]
20 Islam N, Hale R, Taylor M and Wilson A 2013 Physiol. Meas. 34 1103 DOI:10.1088/0967-3334/34/9/1103 [Cited within:1]
21 Le Bihan D, Delannoy J and Levin R L 1989 Radiology 171 853 DOI:10.1148/radiology.171.3.2717764 [Cited within:1]
22 Ishihara Y, Calderon A, Watanabe H, Okamoto K, Suzuki Y, Kuroda K and Suzuki Y 1995 Magn. Reson. Med. 34 814 DOI:10.1002/(ISSN)1522-2594 [Cited within:1]
23 Kuroda K, Oshio K, Chung A H, Hynynen K and Jolesz F A 1997 Magn. Reson. Med. 38 845 DOI:10.1002/(ISSN)1522-2594 [Cited within:1]
24 Shmueli K, Dodd S J, Li T Q and Duyn J H 2011 Magn. Reson. Med. 65 35 DOI:10.1002/mrm.v65.1 [Cited within:1]
25 Diakite M, Odéen H, Todd N, Payne A and Parker D L 2014 Magn. Reson. Med. 72 178 DOI:10.1002/mrm.24900 [Cited within:1]
26 Leroy Y, Bocquet B and Mamouni A 1998 Physiol. Meas. 19 127 DOI:10.1088/0967-3334/19/2/001 [Cited within:1]
27 Maass Moreno R and Damianou C A 1996 J. Acoust. Soc. Am. 100 2514 DOI:10.1121/1.417359 [Cited within:1]
28 Maass Moreno R, Damianou C A and Sanghvi N T 1996 J. Acoust. Soc. Am. 100 2522 DOI:10.1121/1.417360 [Cited within:1]
29 Rahimian S, Tavakkoli J and Meairs S 2012 J. Acoust. Soc. Am. 131 3211 DOI:10.1121/1.4707975 [Cited within:1]
30 Bamber J C and Hill C R 1979 Ultrasound. Med. Biol. 5 149 DOI:10.1016/0301-5629(79)90083-8 [Cited within:1]
31 Arthur R M, Straube W L, Starman J D and Moros E G 2003 Med. Phys. 30 1021 DOI:10.1118/1.1570373 [Cited within:1]
32 Miller N R, Bamber J C and Meaney P M 2002 Ultrasound. Med. Biol. 28 1319 DOI:10.1016/S0301-5629(02)00608-7 [Cited within:1]
33 Song C W, Lokshina A, Rhee J G, Patten M and Levitt S H 1984 IEEE T Bio-Med. Eng. 1 9 DOI:10.1109/TBME.1984.325364 [Cited within:1]
34 Gautherie M 1980 Ann. N. Y. Acad. Sci. 335 383 DOI:10.1111/nyas.1980.335.issue-1 [Cited within:1]
35 Liu D and Ebbini E S 2010 IEEE T, Bio-Medical. Eng. 57 12 DOI:10.1109/TBME.2009.2035103 [Cited within:1]
36 Kinsler L E, Frey A R and Coppens A B et al. 1999 Fundamentals of Acoustics4th edn. New York Wiley-VCH [Cited within:1]
37 Bamber J C and Hill C R 1979 Ultrasound. Med. Biol. 5 149 DOI:10.1016/0301-5629(79)90083-8 [Cited within:1]
38 Bhattacharya A and Mahajan R L 2003 Physiol. Meas. 24 769 DOI:10.1088/0967-3334/24/3/312 [Cited within:1]
39 Simon C, VanBaren P and Ebbini E S 1998 IEEE T. Ultrason. Ferr. 45 1088 DOI:10.1109/58.710592 [Cited within:1]
40 Le Floch C and Fink M 1997Proceedings of the 1997 Ultrasonics Symposium October 5–8, 1997Toronto, Canada 1301 DOI:10.1109/ULTSYM.1997.661816 [Cited within:1]
41 Dhatt G, Lefrançois E and Touzot G 2012 Finite Element Mehtod New York John Wiley & Sons [Cited within:1]
42 Li C P, Duric N, Littrup P and Huang L J 2009 Ultrasound. Med. Biol. 35 1615 DOI:10.1016/j.ultrasmedbio.2009.05.011 [Cited within:2]
43 Li C, Duric N and Huang L 2008International Conference on BioMedical Engineering and InformaticsMay 27–30, 2008Hainan, China 708 DOI:10.1109/BMEI.2008.303 [Cited within:1]
44 Duric N, Littrup P, Poulo L, Babkin A, Pevzner R, Holsapple E, Rama O and Glide C 2007 Med. Phys. 34 773 DOI:10.1118/1.2432161 [Cited within:1]
45 Del Grosso V A and Mader C W 1972 J. Acoust. Soc. Am. 52 1442 DOI:10.1121/1.1913258 [Cited within:1]
46 Parker N G and Povey M J W 2012 Food. Hydrocoll. 26 99 DOI:10.1016/j.foodhyd.2011.04.016 [Cited within:1]