RL-963 Optimal Phased Array Pattern Synthesis for Non-Invasive Cancer Ablation of Liver Tumors Using High Intensity Focused Ultrasound Youssry Y. Botros May 1998 RL-963 = RL-963

Optimal Phased Array Pattern Synthesis for Non-Invasive Cancer Ablation of Liver Tumors Using High Intensity Focused Ultrasound by Youssry Y. Botros A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1998 Doctoral Committee: Assistant Professor Emad S. Ebbini, Co-Chair Professor John L. Volakis, Co-Chair Professor Anthony W. England Professor Thomas B. A. Senior Professor Fawwaz T. Ulaby

~ Youssry Y. Botros 1998 All Rights Reserved

This work is dedicated to my parents, Youssef and Nagat, my sister Faten, my brother-in-law Hani, and my nephew Mina.

ABSTRACT Optimal Phased Array Pattern Synthesis for Non-Invasive Cancer Ablation of Liver Tumors Using High Intensity Focused Ultrasound by Youssry Y. Botros Co-Chairs: Emad S. Ebbini and John L. Volakis In recent years, phased arrays have been employed with reasonable success to apply high intensity focused ultrasound (HIFU) beams to tumors surrounded by homogeneous media to achieve a desired thermal treatment. In the case of tumors surrounded by inhomogeneities, such as liver tumors shadowed by the rib cage, the focusing scheme is highly complicated and has not been addressed in the literature. The problem stems from the need to find an efficient focusing scheme capable of generating well behaved interference patterns in the interior region of the rib cage concurrently with minimizing power deposition over the rib surfaces. In this work, high frequency techniques are employed to predict the ultrasonic fields within the computational domain. Two synthesis techniques are developed and implemented for single and multiple focus pattern synthesis. After experimental validation of the

proposed computational and synthesis approaches, realistic modeling of the problem geometrical features and parameters is considered, analyzed, and tested for single and multiple focusing schemes. This realistic modeling employs exact scattered field representations to predict the ultrasonic fields within the computational domain. Utilizing the proposed computational and synthesis approaches, design curves for selecting the optimal treatment frequency have been constructed, allowing practioners to design treatment plans according to the tumor location and the tissue attenuation. This thesis demonstrated the feasibility of achieving deeply localized HIFIJ treatment in the presence of the rib cage while minimizing direct power incidence over the rib surfaces, and achieving optimal phased array pattern features. A major conclusion of the thesis is that high frequency ray approaches can be efficiently employed for predicting the ultrasonic fields within the domain of interest.

ACKNOWLEDGEMENTS One of the most difficult situations I face now is writing acknowledgments, because expressing my gratitude to everyone who helped me is impossible. My story begins at the computer laboratory, Faculty of Engineering, Alexandria, Egypt. I was trying to contact some professors in the United States in order to pursue my Ph.D. degree. One of my emails went to Professor John L. Volakis at the University of Michigan, who replied on the same day informing me that he could accept me. Three months later, on January 1, 1995, I arrived at Detroit metro airport and Dr. Sunil Bindiganavale, a former Volakis student, picked me up at the airport and gave me a ride to Ann Arbor. On January 2, 1995, John Volakis visited me at the hotel to explain my research project and asked me to have some results quickly. He said, "We have to give them results very quickly'. Until now, I do not know who are "them". Around mid-January, 1995, John Volakis introduced me to my other co-advisor, Professor Emad Ebbini. In May, 1998, I finished my thesis, defended and my conclusion is that "It was fun". I was lucky to have Professors John L. Volakis and Emad S. Ebbini as thesis co-advisors. They offered me both guidance and technical help. Their sincere and continuous encouragement, technical and personal assistance, and great friendship iii

had a great impact on my performance. Many thanks go to my committee members Professor Anthony W. England, Professor Thomas B. A. Senior and Professor Fawwaz T. Ulaby whom all have provided the ultimate help and assistance they could offer. I have to recognize Mr. Tariq Abdelhamid who was an excellent friend and fellow engineer. I met him when I came to Ann Arbor and during my stay, we were contacting on a daily basis and this created an excellent environment and perfect atmosphere for me. I discussed every point of my scientific and personal life with Tariq and I always found good, sincere and honest advice. When I wrote the thesis, he read it twice and from his remarks, I was able to significantly improve my writing style. Special thanks go to his family including his wife Arwa Elnahas and his young daughter Noreen. During my three year stay in the Radiation Laboratory, I made excellent friendships with the faculty, staff and students. I learned much from my colleagues in the Radiation Laboratory and they did not hesitate to help me out when I needed. Specifically, I would like to thank the former members of Professor John Volakis group: Dr. Hristos Anastassiu, Dr. Sunil Bindiganavale, Dr. Jian Gong, Dr. Tayfun Ozdemir and Dr. Dan Ross. Also, I would like to express my deepest appreciation to former students of Professor Emad S. Ebbini group, Dr. Osama Haddadin and Dr. Jian Shen for their sincere help, friendship and encouragement. I would also like to thank all the current members of Professor Volakis group: Lars Andersen, Arik Brown, Michael Carr, Yunus Erdemli, Dejan Filipovic, Zhifang Li, Mike Nurnberger, Kubilay Sertel and Misha Smelyanskiy. Special thanks go to Dr. Thomas Eibert, our post doctoral fellow for reading my thesis and giving very 1V

useful remarks. Also, I would like to thank Dr. Erdem Topsakal for our useful discussions. I would also like to thank Dr. Adib Nashashibi for his cooperation, advice, and encouragement. With regard to Professor Ebbini's group, I would like to thank Phillip VanBaren who significantly contributed to this work by performing an experiment that was critical in validating my model, and Claudio Simon for his helpful and constructive discussions. I also made excellent relations with students in other groups through dinners, parties, gatherings, and soccer games, which all helped me enjoy the stay in Ann Arbor. Technical discussions with all Radiation Laboratory students strongly affected my research and the great ideas, help and support they offered me were effective and helpful. I would like to thank Mr. Mark Casciato, "Frank, the Enemy"' who offered me great assistance during the last stage of my research in addition to accompanying me to very late dinner or breakfast trips. Also, I have great appreciation for my discussions with my friend John D. Shumpert who provided great help in the later stages of my research. I would like also to express my gratitude to Mr. Stephane R. Legault whose cooperation benefited me greatly and my useful discussions with him lead to the achievement of better results. Community wise, I would like to thank all the people I knew outside the university during my stay in Michigan. All the Egyptian students in Ann Arbor were more than friends and I had the chance to make excellent relations with all of them. Special thanks go to the Coptic Community at Michigan who considered me as one of them and I will never forget their kindness and friendship. I would like to thank my family in Alexandria, Egypt. Special thanks and ap v

preciation to my mother Nagat Henri who was a wonderful mother and an excellent friend at the same time. My deep gratitude go also to my father Youssef Botros for his support and encouragement. Appreciation goes to my sister Faten and her husband Hani. Finally, I would like to dedicate this work to the new family member, my nephew Mina who was born while I was here. vi

TABLE OF CONTENTS DEDICATION................................... ACKNOWLEDGEMENTS........................... LIST OF TABLES................................ LIST OF FIGURES............................... LIST OF APPENDICES............................ CHAPTERS 1 Introduction............................... 1.1 Motivation............................. 1.2 Overview of Cancer Treatment Modalities........... 1.3 An Introduction to HIFU.................... 1.4 Significance of HIFU...................... 1.5 Problem Statement and Objectives................ 1.6 Suggested Computational and Synthesis Approaches 1.7 Dissertation Overview...................... ii iii x xi X1i 1 1 2 6 8 9 11 12 2 Problem Geometry and Definitions............................. 16 2.1 Basic Relations......................... 17 2.2 Problem Description and modeling............... 19 2.3 Focusing Definitions.............................. 25 2.4 Thermal Dose Computations.........................28 2.5 Discussion and Conclusions..........................32 3 Focusing in Homogeneous Media......................... 3.1 Formulation............................ 3.2 Computational Example......................... 3.3 Numerical Simulations...................... 3.3.1 Single Focus Pattern Synthesis............. 3.3.2 Multiple Focus Pattern Synthesis.............. 34 35 40 42 44 45 vii

3.4 Discussion and Conclusions................... 4 Hybrid Ray Physical Optics (HRPO) Modeling............ 4.1 Computational Description and Formulation......... 4.2 Computational Example.................... 4.3 Rib Fields............................ 4.4 Experimental Validation.................... 4.5 The HRPO Algorithm...................... 4.6 Numerical Examples....................... 4.6.1 Single Focus Scheme.................. 4.6.2 Multiple Focusing Scheme............... 4.7 Comparison with the Focusing in Homogeneous Media.. 4.7.1 Single Focus Scheme................... 4.7.2 Two Focus Scheme................... 4.8 Discussion and Conclusions................... 5 Hybrid 5.1 5.2 Two-Step Virtual Array Ray (VAR) Modeling........ Computational Description and Formulation......... Field Pattern Synthesis Using VAR.............. 5.2.1 STEP 1: Virtual Array Design............. 5.2.2 STEP 2: Array Synthesis and minimization of intensity levels over the solid Ribs.. 5.2.3 Constrained Optimization... 5.3 The Two Step VAR Algorithm... 5.4 Numerical Examples........... 5.4.1 Single Focus Example...... 5.4.2 Multiple Focus Example.. 5.5 CPU Computational Cost Comparison and VAR Techniques.......... 5.5.1 HRPO CPU Cost........ 5.5.2 VAR CPU Cost......... 5.5.3 CPU Cost Comparison..... 5.6 Discussion and Conclusions....... 52 67 68 73 79 82 86 87 88 90 101 101 107 112 117 118 122 122 123 L 24 1:25 1L26 127 L 34 143 145 146 148 149 152 153 155 162 165 171 172 173 180 between the HRPO 6 Improved Geometrical and Computational Modeling of the Problem 6.1 The Simplified Geometrical Modeling and the Approximate High Frequency Computational Techniques.......... 6.2 Modeling and Computational Formulation............ 6.3 Simple Case: One Line Source Radiates in the Presence of a Lossy Elastic Cylinder...................... 6.4 General Case: K Line Sources and N Lossy Elastic Cylinders 6.5 Pattern Synthesis........................ 6.6 Numerical Examples....................... 6.6.1 500 kHz Operating Frequency with no Attenuation 6.6.2 1 MHz Operating Frequency with no Attenuation. viii.. Vlll

6.6.3 Effects of Tissue Attenuation on the Focusing Scheme 6.7 Discussion and Conclusions................... 7 Design 7.1 7.2 7.3 7.4 7.5 and Computational Aspects................... Attenuation Model and Energy Requirements......... Optimal Power Deposition................... Design Curves.......................... Effect of Rib Power Minimization................ Discussion and Conclusions................... 184 186 206 206 209 213 222 231 8 Conclusions and Recommendations for Future Work........... 238 8.1 Chapters Contributions..................... 238 8.2 Recommendations for Future Work................. 244 APPENDICES.............................. BIBLIOGRAPHY............................... 246 259 ix

LIST OF TABLES Acoustic properties of different tissues..................20 Comparison between different focusing schemes.......... 66 Weights and node values required for numerical evaluation of integrations using Gaussian quadrature scheme................ 2258 Table 2.1 3.1 B.1 x

LIST OF FIGURES Figure 2.1 Schematic description for the geometrical features of the problem at hand................................................. 22 2.2 Simplified two dimensional model of the problem at hand. 23 2.3 Focusing mechanism........................... 24 2.4 Simulation results of the threshold for lesion formation based on (2.20) clearly demonstrating its agreement with the experimental intensityduration formula given in 2.22..................... 31 3.1 Schematic modeling of the focusing mechanism in homogeneous medium. 36 3.2 Two dimensional modeling of the focusing mechanism in homogeneous medium........................................... 37 3.3 Complex excitation vector of the feed array (single focus)....... 46 3.4 Two dimensional intensity pattern (single focus)........... 47 3.5 Focal plane cut in the intensity pattern (single focus)......... 48 3.6 Contour plot based on 25% contours in the intensity pattern (single focus).................................... 49 3.7 Complex excitation vector of the feed array (two foci)........ 53 3.8 Two dimensional intensity pattern (two foci).............. 54 3.9 Focal plane cuts in the intensity pattern (two foci).......... 55 3.10 Contour plot based on 25% contours in the intensity pattern (two foci). 56 3.11 Complex excitation vector of the feed array (four foci)........... 57 3.12 Two dimensional intensity pattern (four foci). 58 3.13 Focal plane cuts in the intensity pattern (four foci). 59 3.14 Contour plot based on 25% contours in the intensity pattern (four foci). 60 3.15 Complex excitation vector of the feed array (six foci).......... 61 3.16 Two dimensional intensity pattern (six foci)............... 62 3.17 Focal plane cuts in the intensity pattern (six foci).............. 63 3.18 Contour plot based on 25% contours in the intensity pattern (six foci). 64 4.1 Computational description of the propagation mechanism. 69 4.2 Geometrical description of the propagation model............ 70 xi

4.3 Two dimensional geometrical modeling of the problem......... 74 4.4 Experimental setup............................ 83 4.5 Two dimensional geometry of the array utilized in the experiment. The wooden bars simulating the ribs are positioned in the y=0 plane along the x-axis with their lengths (depths) running parallel to the z-axis.. 84 4.6 Comparison between measured and computed data........... 91 4.7 Complex excitation vector of the feed array shown in Figure 4.3 (single focus)..................................... 92 4.8 Two dimensional Intensity pattern for the Single Focus Case. The simulated ribs are centered at y=0 and x=30 mm, 69 mm, 108 mm, 147 mm, 186 mm and 225 mm....................... 93 4.9 Focal plane intensity cut at y=150 mm............... 94 4.10 Intensity pattern contour plot (25 % contours) (single focus)... e 95 4.11 Intensity pattern over the rib tissue plane (single focus)......... 96 4.12 Complex excitation vector of the feed array shown in Figure 4.3 (two foci)...................................... 97 4.13 Two dimensional intensity pattern (two foci). The simulated ribs are centered at y=0 and x=30 mm, 69 mm, 108 mm, 147 mm, 186 mm and 225 mm.... v...c 98 4.14 Focal plane intensity cuts at y=135 mm and at y=165 mm...... 99 4.15 Intensity pattern contour plot (25 % contours) (two foci)....... 100 4.16 Intensity pattern over the rib-tissue plane (two foci).......... 103 4.17 Complex excitation of the feed array for the homogeneous scheme (single focus).................................. 104 4.18 Complex excitation of the feed array for the inhomogeneous scheme using the HRPO technique (single focus)................. 104 4.19 Contour plot based on (25 % contours) for the homogeneous schenme (single focus)............................... 105 4.20 Contour plot based on (25 % contours) for the inhomogeneous scheme using the HRPO technique (single focus)................ 105 4.21 Focal plane cut for the homogeneous scheme (single focus)...... 106 4.22 Focal plane cut for the inhomogeneous scheme using the HRPO technique (single focus)............................ 106 4.23 Complex excitation of the feed array for the homogeneous scheme (two foci)....................................... 109 4.24 Complex excitation of the feed array for inhomogeneous scheme using the HRPO technique (two foci)...................... 109 4.25 Contour plot based on (25 % contours) for the homogeneous scheme (tw o foci)................................. 110 4.26 Contour plot based on (25 % contours) for the inhomogeneous scheme using the HRPO technique (two foci).................. 110 4.27 Focal plane cuts for the homogeneous scheme (two focus)....... III xii

4.28 Focal plane cuts for the inhomogeneous scheme using the HRPO (two foci)............................ 111 4.29 Illustration of the reciprocity concept utilized in interpreting the HRPO results.................................. 116 5.1 Two dimensional model illustrating the virtual array concept..... 121 5.2 Complex excitation of the feed array using the VAR method (single focus). Array efficiency VtA= 8.28 %................... 130 5.3 Complex excitation of the feed array using the HRPO method ( single focus). Array Efficiency qTA= 21.93 %.................. 130 5.4 Two dimensional intensity pattern using the VAR method (single focus). 131 5.5 Two dimensional intensity pattern using the HRPO method (single focus)................................. 131 5.6 Focal plane intensity cut using the VAR method (single focus).... 132 5.7 Focal plane intensity cut using the HRPO method (single focus).. 132 5.8 Intensity pattern contour plot (25 % contours) using the VAR method (single focus)................................. 133 5.9 Intensity pattern contour plot (25 % contours) using the HRPO method (single focus)...................................133 5.10 Intensity pattern over the rib-tissue plane using the VAR method (single focus). Relative rib power cr= 2.18 %................ 135 5.11 Intensity pattern over the rib-tissue plane using the HRPO method (single focus). Relative rib power c6= 2.88 %o................. 135 5.12 Complex excitation of the feed array using the VAR method (two foci) (a) Relative amplitudes (b) Relative phases. The Array excitation efficiency 77A= 10.43 %........................... 138 5.13 Complex excitation of the feed array using the HRPO method (two foci) (a) Relative amplitudes (b) Relative phases. The Array excitation efficiency qA= 19.96 %o.......................... 138 5.14 Two dimensional intensity pattern using the VAR method (two foci). 139 5.15 Two dimensional intensity pattern using the HRPO method (two foci).139 5.16 Focal plane intensity cuts using the VAR method (two foci)..... 140 5.17 Focal plane intensity cuts using the HRPO method (two foci).... 140 5.18 Intensity pattern contour plot (25 % contours) using the VAR method (two foci).............................................. 141 5.19 Intensity pattern contour plot (25 % contours) using the HRPO method (two foci)............................................ 141 5.20 Intensity pattern over the rib-tissue plane using the VAR method. Relative rib power 6r= 2.47 %............................ 142 5.21 Intensity pattern over the rib-tissue plane using the HRPO method. Relative rib power ~r= 6.79 %............................ 142 5.22 CPU ratio between the VAR and HRPO techniques............. 151 6.1 Representation of the geometrical features of the propagation model. 157 6.2 Two dimensional geometrical representation of the problem at hand.. 158 xiii

6.3 A line source radiates in the presence of a lossy elastic cylinder.... 159 6.4 Schematic representation of the problem described in Figure 6.3... 160 6.5 Two dimensional model illustrating the field computational concept applied to the problem described in Figure 6.3............... 161 6.6 Geometrical details of the problem.................... 166 6.7 Virtual Array Concept............................... 169 6.8 Complex excitation of the feed array for the single focus example at 500 kHz................................. 175 6.9 Two dimensional intensity pattern for the single focus example at 500 kHz.................................. 176 6.10 Focal plane intensity cut for the single focus example at 500 kHz... 177 6.11 Intensity pattern contour plot (25 % contours) for the single focus example at 500 kHz........................... 178 6.12 Intensity pattern over the rib-tissue plane for the single focus example at 500 kHz.......................................179 6.13 Complex excitation of the feed array for the two focus example at 500 kH z..................................... 181 6.14 Two dimensional intensity pattern for the two focus example at 500 kHz................................................. 182 6.15 Focal plane cut for the two focus example at 500 kHz.......... 188 6.16 Intensity pattern contour plot (25 % contours) for the two focus example at 500 kHz...................................... 189 6.17 Intensity pattern over the rib-tissue plane for the two focus example at 500 kHz........................................ 1.90 6.18 Complex excitation of the feed array for the single focus example at 1 M H z.................................... 91 6.19 Two dimensional intensity pattern for the single focus example at 1 M Hz..................................... 192 6.20 Focal plane intensity cut for the single focus example at 1 MHz... 193 6.21 Intensity pattern contour plot (25 % contours) for the single focus example at 1 MHz.............................. 194 6.22 Intensity pattern over the rib-tissue plane for the single focus example at 1 MHz.................................. 195 6.23 Complex excitation of the feed array for the single focus example at 500 kHz with lossy tissue......................... 196 6.24 Two dimensional intensity pattern for the single focus example at 500 kHz with lossy tissue............................ 197 6.25 Focal plane intensity cut for the single focus example at 500 kHz with lossy tissue.......................................... 1.98 6.26 Intensity pattern contour plot (25 % contours) for the single focus example at 500 kHz with lossy tissue.................. 199 6.27 Intensity pattern over the rib-tissue plane for the single focus example at 500 kHz with lossy tissue........................ 200 xiv

6.28 Complex excitation of the feed array for the single focus example at 1 MHz with lossy tissue........................................ 201 6.29 Two dimensional intensity pattern for the single focus example at 1 M Hz with lossy tissue......................... 202 6.30 Focal plane intensity cut for the single focus example at 1 MHz with lossy tissue................................ 203 6.31 Intensity pattern contour plot (25 % contours) for the single focus example at 1 MHz with lossy tissue................... 204 6.32 Intensity pattern over the rib-tissue plane for the single focus example at 1 MHz with lossy tissue......................... 205 7.1 Variation of the deposition index p with both the focal distance d and the operating frequency in MHz........................ 218 7.2 Deposition index as a function of the treatment frequency for a set of focal distances ranging from 6 cm to 10 cm. 219 7.3 Deposition index as a function of the treatment frequency for a set of focal distances ranging from 11 cm to 15 cm.............. 220 7.4 Percentage increase in the treatment frequency as a function of the allowed deterioration or degradation in the deposition index P..... 221 7.5 Complex excitation vector of the feed array with rib power minimization...................................... 225 7.6 Complex excitation vector of the feed array with no rib power minimization........................................... 226 7.7 Two dimensional intensity pattern when rib power minimization is considered......................................... 227 7.8 Two dimensional intensity pattern with no rib power minimization.. 228 7.9 Focal plane intensity cut with rib power minimization........... 230 7.10 Focal plane intensity cut with no rib power minimization....... 233 7.11 Intensity pattern contour plot (25 % contours) with rib power minim ization................................. 234 7.12 Intensity pattern contour plot (25 % contours) with no rib power minimization.......................................... 235 7.13 Intensity pattern over the rib-tissue plane with rib power minimization. 236 7.14 Intensity pattern over the rib-tissue plane with no rib power minimization..................................... 237 xv

LIST OF APPENDICES APPENDIX A Minimum-Norm Least-Squares Solution................ 247 B Gaussian Integration Technique................... 255 xvi

CHAPTER 1 Introduction The high incidence of liver cancers prompts us to develop more efficient cancer treatment modalities. In this chapter, cancer treatment modalities are summarized with their advantages and drawbacks. Hyperthermia as a non-invasive cancer treatment modality is then discussed. An efficient way of performing hyperthermia is to employ High Intensity Focused Ultrasound (HIFU) to ablate cancer tumors. HIFU is briefly explained with its significant advantages over other cancer treatment modalities. Difficulties associated with its clinical application are summarized. The problem statement, thesis objective, and computational modeling are summarized in this chapter. 1.1 Motivation Liver malignancy is the main motivation for this work due to its significant incidence worldwide. In the United States, an estimated 6000 to 10000 patients per year with colorectal liver metastases are candidates for liver resection or other ablative 1

methods. These methods have a 5 year survival rate of 30 %. In the selected group of patients undergoing surgical resection for colorectal metastases, surgical mortality is about 5 % in major cancer treatment centers, and there is 15-30 % morbidity associated with the operation. Ten to fourteen days of hospitalization are common, and recurrence in the liver occurs in up to 25 % of the patients. In the United States, hepatocellular cancer does not have a high incidence, with approximately 5000 annual cases (one million worldwide), of which only 15 % are candidates for tumor resection. Resected hepatocellular cancer has a five year survival rate of 30 %. In patients with hepatocellular cancer, operative mortality is greater, up to 15 %, despite selection of better risk patients, because many patients with hepatocellular cancer have underlying liver cirrhosis. Liver re-resection for recurrent hepatocellular cancer is frequently, but not usually successful. Due to the importance and significance of liver cancer problems and because of the adverse effects associated with surgical treatments, there is a great need to seek and develop alternative treatment methodologies. In general, cancer treatment methods involve both invasive and non-invasive treatments. Several cancer treatment techniques with different features and mechanisms are briefly presented and discussed in the next section. 1.2 Overview of Cancer Treatment Modalities The major modality for tissue ablation in current clinical practice is surgery. Other clinical applications such as radiotherapy, cryotherapy, laser and phototherapy and 2

chemical tissue infiltration, and hyperthermia are more limited. Surgical ablation is always invasive and requires inpatient hospital stays. As mentioned in the previous section, surgical treatment suffers from several drawbacks, First, its associated risk inside the operating room is around 5 %. Second, postsurgery complications and adverse effects have high occurrence among patients undergoing surgery. Third, due to the way of resecting liver tumors, the percentage of cancer reoccurrence is quite high. Finally, surgical operations are not localized to the tumor volume or region and thus lead to significant damage within surrounding normal tissues and other organs. Another cancer treatment modality is radiotherapy treatment which offers two main advantages. Being non-invasive, it dramatically reduces risks associated with the treatment and post-surgery complications. Another advantage of radiotherapy treatment is that it can be accurately focused to the tumor volume, thus avoiding the damage to normal tissues or surrounding organs. However, due to the utilization of ionizing materials and compounds for this treatment, prolonged usage and/or relatively high doses increase its risks significantly. Cryosurgery, where a very cold needle is utilized to treat a certain area, is invasive and relatively cumbersome resulting in tedious treatments. The major advantage of cryosurgery is the ability to selectively ablate the tumor while preserving adjacent tissue. Laser and photodynamic therapy are not applicable for deep lesions for the most part. When laser therapy is applied to deep lesions, the therapy is accomplished through heating of tissues and thus laser treatment is invasive. 3

Chemical injection for tissue ablation is also invasive, and is not readily controlled. It involves the injection of chemical materials which offer cytotoxic effects for cancer tumors. It may be an efficient treatment methodology if used for small doses or short periods. For deep lesions, ethanol has been the most frequently utilized chemical but its usage is limited with few successful treatment cases. Hyperthermia is a cancer treatment modality which aims at the preferential destruction of cancer or malignant cells. This goal is achieved by increasing the tumor(s) temperatures to therapeutic levels for specified periods of time [54]. Hyperthermia, as an efficient treatment modality for cancer tumors, has been supported by recent biological and clinical research showing that mildly elevated temperatures have cytotoxic effects for solid tumors [53]. These findings have enhanced the validity of using hyperthermia as a cancer treatment modality. Hyperthermia can be combined with other cancer treatment modalities such as radiation therapy and chemotherapy. It is suggested that combined heat and radiation therapy has a synergistic effect on the response of several tumors [51]. This is due to the fact that cancer tumors respond differently to the applied treatment modality. For example, if a tumor has a good response to heat treatment, it will respond poorly to the applied radiation treatment. Therefore, combined treatment is recommended for different types of cancer tumors. Hyperthermia can be applied to the whole body (systemic or whole body hyperthermia) or to a specified region of the whole body (regional hyperthermrnia). Regional hyperthermia can include a region which contains both malignant and normal tissue or can be confined to the tumor region (localized hyperthermia). Also, depending on 4

the depth of the tumor, hyperthermia can be classified either as superficial or deep. Hyperthermia can be applied by various techniques including fluid immersion (systemic hyperthermia), electric currents (interstitial, localized hyperthermia), and electromagnetic and acoustic waves (non-invasive, localized hyperthermia). Among these methodologies, electromagnetic and acoustic external beam heating treatments appear to be more attractive. Absorption of electromagnetic waves in the tissue is due to conduction and displacement currents induced in the tissue region. The best localization when electromagnetic waves are utilized is achieved via microwaves. However, they suffer from severe attenuation (absorption) in biological tissues. Therefore, their application is limited to superficial tumors such as the subcataneous ones. The penetration depth of microwave applicators is around 4 cm to 6 cm [57]. On the other hand, acoustic or ultrasonic external heating is based on the absorption of high energy waves by the tissue [48]. For deep focusing and to increase the power deposition at the tumor, there is a need for multiple transducer systems to carry out an efficient treatment. Each transducer is designed and fed in a such a way that the acoustic or ultrasonic beams converge at the treatment region which is usually the tumor location(s). Multiple transducer applicators are typically found in the form of phased arrays [23, 46, 47, 84] to allow electronic beam scanning and thus remove a need for mechanical translocation of the applicator head. High intensity ultrasonic beams can be focused to the prespecified tumor regions in order to deposit energy (in the form of heat) sufficient to reach certain therapeutic levels required for cancer ablation. This cancer treatment modality is referred 5

to as High Intensity Focused Ultrasound (HIFU) treatment and is briefly defined, introduced and discussed in the subsequent section. 1.3 An Introduction to HIFU High intensity focused ultrasound in the frequency range of 500 kHz to 10 MHz has been used in a wide range of therapeutic applications [20, 21, 35, 51, 52, 53, 54, 66]. For example, HIFU has considerable potential for deep localized hyperthermia, primarily in conjunction with other cancer treatment modalities, e.g. radiation therapy. This is because ultrasound can be easily focused (with mm-size heating spots at the focus) at considerable depths (1 cm - 15 cm with focusing) depending on the frequency of operation [26]. The (typically) small heating spot at the focus can be scanned mechanically or electronically to produce the desired heating pattern tailored to the tumor geometry. Clinical hyperthermia systems utilizing HIFU have been used by a number of groups and their investigations continue to demonstrate the promise of this technology. HIFU has become even more attractive with recent advances in phased array applicators [23, 24, 25, 36, 45, 46, 47]. Phased arrays offer the promise of versatile applicator systems providing several advanced field control features [79, 83] such as: * High-speed electronic scanning. * Tissue aberration correction. * Applicator/patient motion compensation. 6

* Optimal pattern synthesis and control features can be developed to allow adequate therapeutic heating of tumors while avoiding overheating of normal intervening tissue [27, 28, 30]. Despite the obvious advantages of HIFU for non-invasive therapeutic applications and the availability of advanced phased array applicator systems capable of precision deep lesion formation, HIFU is not yet a widely accepted cancer treatment modality in the clinic either for normal hyperthermia or surgery. The reasons for this lack of progress are: * Incomplete understanding of the effectiveness of HIFU to adequately heat specified volumes (e.g. tumors) at deep regions in the presence of inhomogeneities and/or shadowing bone structures. An example of this problem appears when treating liver tumors by focusing high intensity ultrasonic beams inside the rib cage. * Lack of clinical real-time non-invasive high resolution temperature estimation feedback. * Lack of quantitative non-invasive measurement techniques for tissue response to HIFU fields (especially important for tissue ablation). In other words, a non-invasive measurement to evaluate the therapeutic levels of exposure is not yet available. Also, non-invasive means to measure thresholds for irreversible thermal tissue damage in-vivo is not available. * Lack of realistic treatment planning software based on patient-specific geomet 7

rical data sets to help with the optimization of the applicator and the treatment strategy. Furthermore, it would be extremely valuable to develop guidance and visualization mechanisms for therapeutic beams based on a well established imaging modality. 1.4 Significance of HIFU As noted earlier, few applications have utilized HIFU in clinical tissue ablation. Nevertheless, HIFU offers significant promise as an alternative treatment planning modality for tissue ablation. There are potential advantages for HIFU over all current clinical tissue ablation methods. Compared to surgical ablation which is always invasive and often requires inpatient hospital stays, HIFU can be developed to be a non-invasive modality and on an ambulatory basis. In addition, HIFU treatment can be integrated with a tissuepenetrating imaging technique. The advantages of HIFU over radiotherapy are its potential to be actively and interactively linked to an imaging feedback system; thus achieving high degree of focusing accuracy. As presently conceived, however, HIFU is not able to treat large tissue volumes with the intension of principally ablating malignant tissue with lesser damage to the surrounding normal tissues. HIFU shares the major selectivity advantage with Cryosurgery where the treatment is only confined to the tumor region or volume without affecting the adjacent normal tissues. Focused ultrasound is unique in its ability to non-invasively target small and deep volumes for thermal coagulative necrosis [5, 7, 15, 20, 36]. It is possible to target 8

tumors 10-15 cm into the body and safely ablate them while preserving the normal intervening tissue. The registration of prior Magnetic Resonance Imaging (MRI) with ultrasound imaging will help in both the treatment planning and image guidance [45]. As a conclusion, HIFU offers the following advantages for cancer treatment: * It is non-invasive. * It utilizes non-ionizing radiation for treating cancer malignancy. * Can be focused to reach deeply seated tumors. * Can be localized to the tumor volume while preserving the surrounding normal tissues. * Can be applied via phased array applicator systems and thus optimal pattern features can be achieved. * Can be integrated with efficient imaging techniques such as the MRI to provide an integrated cancer treatment [45]. 1.5 Problem Statement and Objectives In previous studies [26, 27], HIFU treatment was performed by utilizing phased array applicators for cancer ablation. Several studies, for example [23, 26, 36, 45, 46, 47, 57, 60, 84], were performed to test phased array usage for deep localized hyperthermia in homogeneous media (free of any strongly scattering obstacles). In these studies, phased arrays with different geometrical features, feed mechanisms and 9

focusing characteristics were tested and implemented to demonstrate the feasibility to deeply focus ultrasonic beams inside a homogeneous medium. Because of the aforementioned significance of liver cancers and due to the advantages of applying HIFU to treat deeply seated organs, HIFU treatment of liver cancer using phased arrays is the emphasis of this thesis. In the case of liver malignancy, cancer tumors are deeply residing inside the human rib cage. Therefore, for a noninvasive hyperthermia treatment using HIFU, it is required to deeply focus ultrasonic beams through the human ribs to reach the liver. This can be achieved by guiding the beams through the intercostal spacings between the human ribs. This task is quite complicated and the presence of the human ribs presents a challenging problem which is addressed in this thesis. The objectives of this thesis are: * Test the feasibility of performing deep localized treatment for liver tumors using HIFU. In this thesis, phased array applicators are utilized to apply HIFU to the human liver inside the rib cage. * Develop, formulate, and validate field computational approaches pertinent to the propagation mechanisms that allow sufficient accuracy for the subject treatment requirements. * Develop and implement optimal phased array pattern synthesis techniques capable of producing the desired therapeutic levels at the tumor location(s) without excessive heating at undesired locations. * Develop and generate design curves that aid in choosing the operating frequency 10

required for an efficient and feasible treatment. 1.6 Suggested Computational and Synthesis Approaches The appropriate choice of the employed computational technique depends on the size of the computational domain and the desired accuracy levels. Two major categories of computational techniques may be employed to predict ultrasonic fields within the domain of interest. The first uses exact numerical methods [37] while the second uses approximate high frequency ray methods [22, 64]. According to the frequency ranges utilized for HIFU treatments and based on tissue layer characteristics [53, 54], the operating acoustic wavelength A within the computational domain is within the mm order. Therefore, with the geometrical features and dimensions of the human body, the computational domain is expected to span tens of wavelengths. Despite significant progress made in numerical methods such as Finite Element Method (FEM) [73, 82], Method of Moments (MoM) [4, 43], Finite Difference Time Domain (FDTD) [59] computational methods, the size of the problem at hand is significantly larger than problems routinely handled by these methods, even on highperformance parallel computers. A 2-D 200A x 200A region of interest (typical for ultrasound) results in millions of grid points with any of these methods resulting in large Central Processing Unit (CPU) and memory requirements. Moreover, these numerical techniques make the synthesis process an intractable task. As a conclusion, 11

numerical methods based on fine tessellation of the computational domain are not appropriate for modeling the propagation mechanism. Alternatively, high frequency ray approaches look more attractive in modeling such a geometry. Therefore, the computational formulation in this thesis is carried out using high frequency ray approaches [2, 17, 22, 58, 64, 71, 81]. These techniques have been shown to be sufficiently accurate when the diffraction points (ribs) are sufficiently apart and the observation points (foci) are also at large distances (with respect to the wavelength) from the scatterer or source. Computationally, high frequency ray techniques offer affordable memory and CPU requirements which are dramatically less than those required by numerical techniques. For the synthesis approaches employed to synthesize the desired field patterns with specific heating foci, pseudo-inverse (PI) techniques introduced in [26, 27, 29] are utilized to compute the feed structure that is capable of generating the desired heat map inside the rib cage. This technique will also be modified to achieve optimal pattern characteristics when deep HIFU treatment inside the rib cage is considered. 1.7 Dissertation Overview The thesis consists of 8 chapters and is organized in the following way: Chapter 1 gave introductory information about the significance of liver cancer and its treatment methodologies. Hyperthermia techniques were summarized and the importance of HIFU as an efficient non-invasive modality was briefly explained. Also, in the first chapter, the problem statement and thesis objectives and goals were 12

clearly defined. The major methods employed for problem formulation were reviewed and high frequency ray methods were recommended for problem modeling. Chapter 2 reviews the basic acoustic relations and quantities and introduces a geometrical description of the problem at hand with its associated parameters. Also, focusing parameters as well as theoretical and empirical relations for the necessary thermal doses are summarized. Actual dimensions and real biological tissue parameters are utilized to propose a model for the propagation mechanism. As an initial and crucial step before addressing the problem in this thesis, Chapter 3 analyzes the problem of focusing ultrasonic beams within homogeneous media with no obstacles. This case represents an ideal focusing reference, which will be compared to the focusing obtained in the presence of the ribs. Also, in Chapter 3, the basic synthesis formulation is introduced and discussed in the context of previous investigations [24, 25, 26, 27]. After presenting the basic computational approach and synthesis algorithms for focusing in homogeneous media, several examples for single and multiple focus pattern synthesis are implemented and discussed. In Chapter 4, a high frequency Hybrid Ray Physical Optics (HRPO) computational approach is developed to predict the ultrasonic fields within the computational domain. The ribs are considered as an array of planar non-penetrable strips. The field values are specified over a set of control (field) points and the pseudo-inverse synthesis approach is applied to compute the feed array. To check that the synthesis process yields the desired field specifications, field maps are constructed everywhere inside the rib cage as well as over the rib surfaces. The setup and results of an experiment performed to validate the HRPO approach are discussed. Single and multiple 13

focus pattern synthesis examples are implemented and tested. Because the HRPO technique lacks explicit control of the rib fields and due to its intensive computational demands, another high frequency technique is developed in Chapter 5 for the analysis and synthesis of the field patterns, which is named the Virtual Array Ray (VAR) technique. This technique which involves two main synthesis steps is implemented and validated using the HRPO method before applying it to single and multiple focus pattern synthesis. It offers several advantages such as ease of implementation and convenience for field pattern synthesis and scanning in addition to explicit control of the ultrasonic fields over the rib surfaces. Computationally, its CPU cost is significantly less than that of the HRPO. It should be noted that in both chapters 4 and 5, the model proposed in chapter 2 is utilized to demonstrate the feasibility of achieving well behaved interference patterns in the presence of strongly scattering obstacles. The main reason for analyzing such a simplified model with approximate high frequency techniques [17, 22, 81] is to determine if focusing deeply inside the rib cage is a feasible task or not. Failing to achieve such a goal makes the pattern synthesis process meaningless. After demonstrating the feasibility of generating well behaved patterns inside the rib cage in Chapters 4 and 5, a more realistic geometrical model is proposed in Chapter 6. For this model, the computational formulation is carried out using the exact eigen solution (Neumann expansion for the fields within the computational domain). The eigen expansion of the domain fields [22, 39, 64, 71, 81] coupled with the boundary conditions are employed to predict ultrasonic fields in the interior of the rib cage and over the rib surfaces. This chapter also provides support for the feasibility of 14

achieving deep localized focusing with well behaved patterns within the domain of interest even with realistic modeling of the rib obstacles. Fields maps constructed using the model proposed in Chapter 6 are similar to those predicted by high frequency approaches introduced in Chapters 4 and 5. This shows the robustness of the employed synthesis and field prediction methods. Chapter 6 also demonstrates that tissue attenuation plays a key role in resulting interference patterns. Simply, higher attenuation values deteriorate the focusing scheme in addition to causing excessive heating at undesired locations. This restricts treatment frequency choices and/or the maximum penetration depth. Therefore, in Chapter 7, an efficient treatment plan is proposed taking into account the attenuation factor and the depth of the tumor as well as the required power deposition level at the tumor location(s). Design curves relating the operating frequency and penetration depth to the heat deposition are generated for use in treatment planning. In addition, Chapter 7 discusses the effect of rib power minimization with the realistic modeling of the problem to conclude that rib power minimization which leads to a constrained optimization problem significantly reduces the power deposition over the rib surfaces. Finally, thesis contributions, conclusions and recommended future work are summarized in Chapter 8. 15

CHAPTER 2 Problem Geometry and Definitions This chapter introduces the basic quantities and relations utilized to predict the behavior of acoustic fields when propagating through homogeneous media. It summarizes the boundary conditions associated with field propagation through inhomogeneous regions. After reviewing these principles, a geometrical model for field prediction is proposed. This model is simplified to a two dimensional one employed to test the feasibility of achieving highly focused beams in the presence of the rib cage. The simplified model is extensively utilized in Chapters 4 and 5 to compute the ultrasonic fields within the computational domain. Focusing definitions and parameters are introduced in this chapter as well. Finally, based on previous investigations [21, 34, 54], thermal dose computational analysis is presented to characterize the intensity, time and temperature effects when performing thermal coagulation of tissues. 16

2.1 Basic Relations The acoustic or ultrasonic field behavior is governed by the following relations [41, 49, 61, 75], p- -( t) Vp(rt)=p t (2.1) 1 Op(. t) V U(, t)=- p(Fc2(f) at ' (2.2) ' P PW) C2 at where p(r, t) is the sound pressure (scalar) at time t and location r' in Newtons/rnm2, u'(r,t) is the particle velocity (vector) at the same time and space coordinates in m/sec, p(r) is the medium density in kg/m3 and c(r) is the speed of acoustic waves (sound) at the location r in m/sec. It is more convenient to express both p and u in terms of a scalar potential 4 which is the velocity potential in m2/sec. This potential is related to the basic acoustic quantities p and u by the following equations [49, 75, 86], (r, t) = Vr(r, ) (2.3) 04(', t) p(ft) =-p(r-) at (2.4) For steady state harmonic time variation with eJwt time dependence assumed and suppressed, the wave equation for 4 can be written in the following standard form [41, 86], V (V$(r)) + Y2(Frj(r) = xS,(r, (2.5) where q(r) is the velocity potential at point r' and y is the complex propagation constant. The source term x8(r) denotes the source amplitude distribution in sec-1. 17

In conjunction with (2.1) to (2.5), the appropriate boundary conditions required for a unique solution of 0 are also introduced. These boundary relations describe the behavior of acoustic waves at interfaces between different media [41, 61, 75] and at the boundaries of the domain. The boundary conditions imply the following: * Continuity of pressure across the boundary separating the two media. Therefore, for an interface separating two media, for example 1 and 2, the boundary condition for the acoustic pressure is P1 - P2 = 0, (2.6) where pi and P2 are the interface pressures in medium 1 and 2, respectively. * Continuity of the normal component of the particle velocity across the boundary between two media. Therefore, the boundary conditions for the acoustic particle velocity is given by n' ( -2). 0, (2.7) where n is the unit normal vector outwardly directed from medium 1 and ui and ui2 are the interface particle velocities in medium 1 and 2, respectively. The general solution of (2.5) is given by the following integral form [17, 86], ()= j,x(r )gf(r'lr")dV' - f [gf(1r')(A'- V'O(F"))- A(i")('. V'gf(i?'))]dS', (2.8) where gf(rlr') is the appropriate space Green's function. By definition, gf(r r') represents the field at a distance r'- r"' from a unit source located at r"'. Under the 18

surface integral, O(rF') is the total ultrasound field on the surface S' which encloses the volume V' and n' is the unit normal vector outwardly directed from S'. For sources in homogeneous media, S' can be allowed to go to infinity and be eliminated from (2.8). Alternatively, the surface integral is used when the sources are not known explicitly and instead the field is given over the surface S' as will be seen later. The time average intensity vector (in Watts/m2) of an acoustic wave is defined by [75, 86], () = Real{p(r) u*(r)}, (2.9) where the asterisk denotes complex conjugation and Real implies the real part. When the time average intensity vector is integrated over the surface 5, the time average acoustic power in Watts is determined from Pav= Iav(ru dS, (2.10) where dS is an infinitesimal element of the surface S and Pay is the time average acoustic power in Watts. 2.2 Problem Description and modeling Figure 2.1 shows the suggested geometrical model for the problem at hand and it also displays a general structure of the feed array utilized in focusing. As shown, a large portion of the liver is located inside the rib cage below the ribs which are considered to be strongly scattering inhomogeneities for the ultrasound waves. A water layer (bolus) is typically placed between the human body and the array for 19

Tissue Attenuation Speed of Sound Specific Impedence Np/(cm MHz) m/sec gm/(cm2 sec) x 10-5 Blood 0.017 1570 1.62 Bone 1.64 3183 4.8 Fat (subcataneous) 0.07 1478 NA Fat (peritoneal) 0.24 1490 NA Heart 0.23 1571 1.64 Liver 0.135 1604 1.63 Skin 0.197 1720 1.59 Table 2.1: Acoustic properties of different tissues. coupling purposes. Directly below the ribs, a fat layer exists above the viscera layer and the liver. For non-invasive treatments, the goal is to compute the excitation of the feed array so that highly localized foci are generated at the target (tumor) location without excessive heating at undesired locations. A second goal as part of the design process is to keep acceptable (minimum) intensity levels over the ribs to minimize the possibility of excessive heating in bones. The thermal conductivity of the bones is significantly higher than that of the soft tissues (four to seven times higher) [52, 53]. Therefore, for the same intensity level, temperature increase in the bones is considerably higher than that in soft tissues (such as skin, fat, viscera and liver). For the HIFU frequency ranges, the properties of the biological layers of the human 20

body are presented in Table 2.1. These properties depend on several parameters such as the person age and the laboratory measurements. A careful examination of Table 2.1 indicates that the specific acoustic impedance for the skin, fat, viscera and liver are so close that the propagation medium can be approximated by a homogeneous one (except for the ribs). Layer attenuation varies from water which is almost lossless to fat and liver tissues which are characterized by relatively high losses. Based on previous investigations and tests [26, 53, 54], the effective or average attenuation for soft biological tissues could be assumed to be 1 dB/(cm MHz) which corresponds 0.115 Np/(cm MHz). For the frequency ranges utilized for HIFU treatments, the operating wavelength in water and soft tissue layers is at most 3 mm (since the lower range for the desired treatment is around 500 kHz). The ribs tend to be cylindrical in shape with elliptical cross section. For an average person, the major and minor diameter of each rib are around 1.8 cm and 1.5 cm, respectively while the intercostal spacing between two adjacent ribs is about 2 cm. The ribs average depth in the axial direction is around 15 cm which is much larger than operating wavelength. This simplifies the three dimensional model shown in Figure 2.1 to a two dimensional one as shown in Figure 2.2. As displayed in Figure 2.2, the feed array is assumed linear with elements radiating in a water layer (bolus) which is placed over the rib cage. The ribs are assumed to be thin, solid and non-vibrating plates placed in the y = 0 plane. 21

Feed Array The Liver & the Tumors Figure 2.1: Schematic description for the geometrical features of the problem at hand. 22

Feed Array 1 /2 3 -s K The Liver & the Tumors Y ----- Figure 2.2: Simplified two dimensional model of the problem at hand. 23

I - ----- - - T Al Aw F F f number = D Figure 2.3 Focusing mechanism. Figure 2.3: Focusing mechanism. 24

2.3 Focusing Definitions Phased arrays are multiple transducer mechanisms that consist of similar elements with certain spatial adjustment [3, 32, 76]. Each element is assumed to be fed or driven by a separate source or signal. This feed generates a certain amplitude and a specific phase at each element. Phased arrays are utilized to improve the power handling capability of a certain system or to achieve a specific type of beam forming or directivity [3, 32]. By changing the feed excitation mechanism of each element, the radiation characteristics change. In general, for each array element, the amplitude and phase of the driving signal can be changed in order to achieve desired field pattern specifications. For a phased array, focusing is achieved by driving the array elements with signals shifted with respect to each other in such a way that the individual beams from the array elements add constructively at the desired focal point. Figure 2.3 shows a generic focusing structure with aperture size D and focal spot distant F from the surface. The dimensions of the focal point in the longitudinal and transverse directions are denoted by Al and Aw, respectively and given by [26], (F) D (2.11) and A W = K2 A D' (2. I12) where A is the wavelength and K1 and K2 are constants determined by the geometry of the focusing structure. The ratio D-is called the f number of the array and it is 25

less than unity for typical hyperthermia transducers. The basic terms and definitions that are frequently used when dealing with focusing ultrasonic waves using phased arrays are as follows:. * The spatial-peak temporal-average intensity at the focus of a radiating source (time harmonic variation) is given by ISPTA = 2pcT( )- (2.13) where rf defines the focal point location relative to the source, r and T are the duty cycle and the period of the gating function, respectively and p is the complex acoustic pressure. * The spatial-peak pulse-average intensity at the focus when the source is driven by a periodically gated sinusoid is ISPPA = p*(c )p ) (2.14) where the average is taken over the pulse duration. Therefore, the time average intensity is related to the pulse average intensity by the following relation TISPPA / c ISPTA = TIP (2.15) T * The spatial-average time-average intensity ISATA, of a beam is the time average acoustical signal averaged over the effective cross section of the beam (defined by the 6-dB lateral dimensions). * The spatial-peak temporal-peak intensity, ISPTP is the highest peak value of the 26

instantaneous intensity at a spatial peak. It is given by 2 ISPTP -, (2.16) pc where Pm is the maximum instantaneous pressure. * The intensity gain at the focal point of a transducer is given by G ISPTA (2.17) ISATA,T where ISATA,T is the spatial-average time-average intensity over the transducer surface. A general formula for the intensity gain of a focused transducer is given by D2 G = K3 (2.18) where K3 is a constant determined by the geometrical features of the focusing system. * The array excitation efficiency for a phased array consisting of K elements driven independently by time harmonic signals is given by A u- > 100, (2.19) where u is the array excitation vector consisting of K complex numbers u= [U1,U2,.......UK]T, IUmaxI is the maximum amplitude particle velocity at the surface of the array and <.,. > defines the inner product of two complex column vectors. Therefore, the array excitation efficiency is a measure of the amplitude (velocity) distribution uniformity over the feed array. The array efficiency is 27

higher when the amplitude distribution is more uniform and becomes unity for a uniformly fed array. The array excitation efficiency is of increasing importance when dealing with limited power handling capability. 2.4 Thermal Dose Computations The concept of thermal dose, while not completely developed, has been accepted by most hyperthermia practitioners. It has been verified experimentally that the lesion size and shape is determined by the thermal dose parameter of equivalent minutes at 43~C, defined as T43(F, t)= RT(rT)-43dr, (2.20) where T is the temperature in ~C, t is the time in minutes, r' is a vector representing the coordinate system of the region of interest and R is an empirically determined constant. It should be noted that (2.20) is derived from cell survival rates at different temperatures [67]. An important implication of (2.20) is that the same desired therapeutic end-point (e.g., coagulation necrosis) can be reached by maintaining tissue temperature at 43~C for 120 minutes or, alternatively, maintaining tissue temprature at 55~C for 2 seconds. These two alternatives are equivalent for non-localized heating of non-perfused tissue, however, they differ greatly when localized heating of perfused tissue is considered. In case of long duration heating, heterogeneous blood perfusion and the thermal conductivit the thermal conductivity of the tissue complicate the requirements of temprature control algorithms, even with the most advanced applicators. On the other hand, short duration heating causes thermal coagulation before blood perfusion and con 28

duction begin to affect the heating pattern. In the latter case, the temprature control problem is reduced to a pattern synthesis problem. This siginificantly simplifies the therapeutic procedure to achieve the desired end-point. This conclusion is based on the theoretical analysis of the transient bioheat transfer equation (BHTE) which has the following form 9T pC- t = V (KtVT) - WbCb(T - Ta) + Q, (2.21) where T is the tissue temperature and it is mainly a function of the position r and the time t. The parameter Ta represents the blood temprature, p is the tissue density, C is the tissue heat capacity. In this equation, KIt, Wb and Cb are tissue conductivity, blood perfusion rate and, blood heat capacity, respectively. The excitation term Q in this equation accounts for the external heating by the incoming ultrasonic waves. In spite of its limitations, especially in accounting for blood perfusion, (2.21) has been generally accepted as a model for tissue temperature response to heating fields [26, 74, 79]. In [26], extensive simulations with the BHTE were performed and compared to both in-vitro and in-vivo experimental data. The results encouraged the use of the BHTE for both prediction and control of the temperature response to localized heating patterns in tissue media. Significant amount of work has been done to determine the thermal and cavitational thresholds required for thermal coagulation of biological tissues. The tests were performed in exposed organs of small animals and the results indicate that thermal coagulation of tissue occurs according to a well-defined intensity exposure relationship 29

[34, 21], I Te = Co (2.22) where I is the intensity in Watts/cm2, Te is the exposure time in seconds and Co is a constant that depends on the tissue type [21, 34, 54]. In [34], Frizzell published invivo data supporting this model for values of Te is in the range from 0.1 to 10 seconds. In his experiment, intensity thresholds for thermal coagulation were compared for cat liver, rabbit liver and cat brain in-vivo for exposure times or pulse durations ranging from.001 second to 35 seconds. The following observations can be made from his work: 1. Thresholds for liver damage were the same for both cats and rabbits. 2. For Te < 0.1 sec, intensity thresholds were higher than those predicted by (2.22). This may be due to the change in the coagulation mechanism from transient to mechanical (transient cavitation). 3. For T, > 10 sec, intensity thresholds were lower than those predicted by (2.22) and this may be due to the blood perfusion effects. 4. For 0.3 < T < 10 sec, the constant C was 460 Watts sec for the a e liver and 200 Watts sec'5/cm2 for the brain. Therefore, different types of tissue will have different threshold values. Figure 2.4 illustrates the agreement between the simulation results of the lesion formation threshold based on the solution of (2.21) and the experimental results using 30

W-I )( J X. S.I a) CO CD 0. 1000 900 800 700 600 Computed threshold -- lntensity*sqrt(Time)=Co............................................................................. -. - ~ ~ ~~~~~~~~~~~.............................................. -......................................................................... ~~....:.~~........: -..... 500 400 300 200 2 3 4 Pulse width (seconds) 5 6 7 8 9 10 Figure 2.4: Simulation results of the threshold for lesion formation based on (2.20) clearly demonstrating its agreement with the experimental intensity-duration formula given in 2.22. 31

(2.22). This agreement exists for exposure durations between 1 and 10 seconds. In addition, this agreement is almost independent of the specific value of the blood perfusion. Therefore, if ultrasonic power deposition is controlled for a specified time duration (1 - 10 seconds) at a given (deep) location inside the tissue, then thermal necrosis will occur at all points within the volume reaching the intensity threshold given by (2.22). 2.5 Discussion and Conclusions This chapter is an introduction to the basic acoustic relations, model parameters and focusing definitions. The geometrical model of the problem at hand was presented in this chapter along with its ultrasonic tissue properties. The chapter also introduced a less complicated model that will be extensively utilized to evaluate if focusing through strongly scattering obstacles is feasible. The ultrasonic properties of the soft tissues enhanced the approximation of homogeneous propagation medium (except for the ribs). Also, in this chapter, the good agreement between empirical and theoretical thermal dose computations indicate that focusing high intensity beams to certain tissues will increase the localized intensity and thus achieve thermal necrosis required for lesion formation. Therefore, keeping the temperature within a specific volume at therapeutic levels for short amounts of times (few seconds) is the basic idea for performing HIFU treatment for this type of tissue. To achieve this task, the beam has to 32

be focused precisely, making the cancer treatment problem equivalent to a focusing or beam forming problem. Thus, the problem of synthesizing temprature heating patterns is converted to a field pattern synthesis problem which is considered and addressed in all subsequent chapters. Before addressing the field pattern synthesis problem in the presence of the rib cage, it is more convenient to carefully study a simpler problem involving focusing ultrasonic beams in homogeneous media (free of any strongly scattering inhomogeneities). This will be addressed and implemented in the following chapter. 33

CHAPTER 3 Focusing in Homogeneous Media In this chapter, ultrasonic focusing in uniform homogeneous media using phased arrays is discussed. Focusing in such media is considered since it represents an ideal application for HIFU. Thus, the resulting focusing scheme may be utilized as a basis for comparing results in the subsequent chapters when the ribs are considered. As concluded in Chapter 1, high frequency ray techniques [15, 17, 22, 39, 64, 71, 81] are employed to analyze the computational model. The basic synthesis techniques and algorithms are introduced, implemented and discussed in this chapter. Specifically, the pseudo-inverse (PI) pattern synthesis technique presented in [26, 27, 29, 30, 38] is proposed as an efficient synthesis modality. This synthesis approach is briefly discussed since it will be employed in subsequent chapters. Several computational examples representing single and multiple focusing schemes are simulated and discussed at the end of the chapter. 34

3.1 Formulation Figure 3.1 demonstrates a schematic description of the focusing mechanism in homogeneous medium. As shown in this figure, a focusing array/aperture is employed to focus the ultrasonic beams inside a homogeneous medium. If the homogeneous medium is assumed to have ultrasonic characteristics similar to that of the water and according to the operating wavelength (less than 3 mm) for HIFU treatments, the model presented in Figure 3.1 can be simply reduced to the two dimensional one shown in Figure 3.2. This simplification is valid if the feed structure has cylindrical shape. That is, the array depth is much larger than the operating wavelength. With the elements in the form of cylindrical line sources, the problem will be equivalent to an array of line sources radiating in a homogeneous medium as shown in Figure 3.2. For the geometry shown in Figure 3.1, the radiated velocity potential from the feed aperture in the homogeneous region (under the aperture) is obtained by considering the first term of (2.8), that is U(F) = j u(r')gf(Fr')dS', (3.1) where r' and r' are the location of the source and field points, respectively and u,(r'') is the normal velocity source distribution at r'. As previously mentioned, gf(rlr') defines the Green 's function which gives the field at r' when a unity source is placed at r'. Discretizing u, in (3.1) into K segments converts the integration in (3.1) into a summation over the aperture segments. Thus, the velocity potential can be written 35

Feed ArraY Homogeneous bledium 0 hornOof the f0clIsIng inechanis'" "I I Schematic Tnodel'ng I Figure 3-1.. ous ynediumgene 36

d 0.< - - Feed Array 1 y Figure 3.2: Two dimensional modeling of the focusing mechanism in homogeneous medium. 37

as K (m) = uk gf(r'mrk)dSk, (3.2) k=-1 k where uk is the complex excitation or feed value of the kth segment located at r' while rm is the mth field point location. Also, gf (rm rk) is the appropriate Green's function which gives the field at rm when a unity source is placed at rk. Alternatively, (3.2) can be written in the following form = Hu, (3.3) where u is the the complex excitation vector of the feed array of length K. The kth element of this vector is equal to Uke3ek, where Uk and Ok are the relative feed magnitude and phase of this element, respectively. In this formulation, the matrix H represents the transfer or propagation operator from the aperture surface to the field point and its entries are given by H(m, k)= = gf(rm Irk) dSk. (3.4) k Therefore, if the excitation vector as well as the field (control) points are specified, the corresponding complex velocity potential can be directly obtained from (3.3) and (3.4). However, in the synthesis process, the complex pressure vector or the velocity potential vector is specified over a set of field or control points and the array/aperture that achieves the desired levels at the prespecified locations is computed. Also, this array should minimize the power deposition elsewhere in the system. For example, for a given set of M control points, the complex control field (pressure) vector &c, is specified and the array with an excitation vector u that minimizes i4c - H ull2 is 38

computed by finding the complex excitation vector u. If M is less than K, (3.3) will be representing an underdetermined system of equations. To deal with such systems, a singular value decomposition technique [38] is utilized as in [15, 26] to solve this system based on a minimum-norm least-square solution for (3.3). This solution is given by [29, 31], u = HP' c, (3.5) where HP' is the pseudo-inverse (PI) of H and it is given by [26, 56], Hpi H*t (H H*t)-1 (3.6) where the superscripts * and t denote matrix conjugation and transpose, respectively. More details about the Pseudo-inverse (PI) minimum-norm least-squares solution are given in Appendix A. The minimum-norm least-squares synthesis is utilized here and in all subsequent chapters of the thesis because of its following advantages [26, 27]: 1. The minimum-norm least-square solution guarantees the exact power deposition at the prespecified control points. This allows for a precise and explicit control of the heating patterns especially at the control locations which are the tumor location(s) and any other points where the power deposition is controlled [5, 7, 13, 15]. 2. The field intensity or pressure distribution produced by the PI solution tends to concentrate (deposit) the energy around the control points and thus minimum energy deposition is observed at other locations. The applicability of the 39

pseudo-inverse synthesis technique to generate efficient focusing schemes has been demonstrated in [6, 14, 26]. 3. The PI solution allows for optimization of several quantities such as the intensity gain at the control points and the array excitation efficiency [26, 27, 29]. 3.2 Computational Example Having presented the theoretical formulation for the synthesis technique, the two dimensional model shown in Figure 3.2 is considered. As shown, an array of K line sources or infinite cylinders is placed on the top of a homogeneous medium occupying the region (y > 0). Assuming that the kth array line source is located at (Xk, yk), the velocity potential due to this element at a distance r with coordinates (x, y) will be given by the following equation [86], uo H^(-jyr2 1S = Ho(2)(-jr) F(p), r > a, (3.7) JY H?2)(-jya) where a is the radius of the vibrating line source or array element, UO is the maximum amplitude velocity at the surface of the transducers. In this equation, H(2)(-j-ya) and H 2)(-jya) are the zero and first order Hankel functions of the second kind, respectively. It should be noted that contour integration was used to determine the reference velocity potential at the surface of the vibrating cylinder. The function F(O) represents a general form of the element directivity function and it can be assumed as the conventional sinc function for the case of line source excitation. Finally, the 40

propagation factor y is given by [42], = a +j/3, (3.8) where a represents the attenuation factor of the ultrasonic waves in the medium in Np/m and 3 is the corresponding phase factor in rad/r. For our application, the primary line sources are placed several wavelengths away from the field (control) points and thus the large argument approximation of the Hankel function may be used to rewrite (3.7) as 1= -0 - F(,), (3.9) with -Uov/1 -= yH(2)(_ jHa)/- (3.10) where r is the distance between the source and field points. The total field at the mth arbitrary point (xmrn Ym) in the far zone is computed by adding the contributions of all the array elements at the specified location. This gives K +(r) = Z ls(rkm), (3.11) k=l where rkm is the distance between the mth field point and the kth array element, given by Xkm = \/(Xk - Xm)2 + (Yk - Ym)2 (3.12) Substituting (3.7), (3.9) and (3.10) into (3.11), the field at the mth test point when the kth element of the linear array is driven by the complex excitation UkCk1 is 41

obtained. This field is expressed as K /km (r)=ke e F (tan1 Xm - Xk (3.13) k=1 V/ rkm I ym-YkI In matrix form, (3.13) can be expressed as = H u, (3.14) where u is the array excitation vector with length IK and its kth complex element equals ukej3k where Uk and 0k are the relative velocity amplitude and phase of the kth array element, respectively. The general entry of the matrix H is given by e —y/km X - Xk H(m,k) = o --- F tan-1 Xm k). (3.15) V7 [m |ym - Yk]A' Next, the synthesis problem is carried out by defining a set of M control points and the corresponding complex potential or intensity value at each one. In this process, the matrix operator at each point is computed according to (3.15) and thus the complex excitation vector u of the feed array is obtained using the PI approach. After computing the feed array, the field map in the homogeneous medium is reconstructed using (3.10) to (3.13) to check the performance of the synthesis algorithm. 3.3 Numerical Simulations In this section, the developed approach is validated by considering several computational examples. In the simulation examples presented in this chapter and the following ones, the problem geometry, dimensions, array shape, elements, operating conditions, and layer parameters are kept the same for comparison purposes. In Chapters 3, 4 and 5, the following assumptions are employed: 42

* The simplest linear array is utilized for testing and validating the suggested computational and synthesis approaches. This array can be considered as a reference array or aperture. Several studies [26, 68, 70, 80, 85] discussed the effects of the array shape, elements and dimensions on the focusing scheme. Also, most of these studies optimally selected these parameters for different focusing schemes and applications. * Optimization of driving signal distribution (excitation of the feed array) or maximization for the intensity gain at the focal spots when focusing in homogeneous media were previously studied [26, 27] and are not addressed here. In the computational examples, a 240 element linear array with A/2 uniform interelement separation was utilized. The operating frequency is 500 kHz which corresponds to a 3 mm operating wavelength in the homogeneous medium (based on a speed of sound of 1500 m/sec). In the following subsections, two types of focusing schemes are simulated, implemented and discussed. The first deals with single focus scheme while the second considers multiple focus schemes. For the single focus case, it is required to create a single spot with a prespecified complex pressure value or intensity level at a desired location inside the tissue or the medium. In multiple focus schemes, more than one location are selected and the feed array has to deposit prespecified levels at their locations. In all numerical simulations considered in this chapter, the following steps are carried out to validate the suggested synthesis approach [26, 27]: * The control locations and the corresponding velocity potential vector '4c are 43

specified. It should be noted that the complex pressure is related to the corresponding velocity potential by (2.4). * Based on linear array excitation with array elements in the form of line sources, the matrix operator H is calculated according to (3.15). * Using (3.5) and (3.6), the complex excitation vector of the feed array is obtained based on the PI synthesis technique [27, 26]. Thus, amplitude and phase distributions over the feed array elements may be computed. 3.3.1 Single Focus Pattern Synthesis In this section, single focus pattern synthesis is implemented and tested for the problem described in Figure 3.2. The focus intensity and its location are chosen arbitrarily. Let us consider focusing at a desired location with the coordinates (120 mm, 150 mm) which is equivalent to (40 A, 50 A) at 500 kHz, where A is the operating acoustic wavelength. The focal intensity is selected to be unity (normalized). The simulation results are demonstrated in Figures 3.3 to 3.6. From these figures, the following can be deduced: * The complex excitation vector of the feed array is obtained and plotted in Figure 3.3. This figure shows that PI synthesis yields a nearly uniform fed array (phase synthesis only). This was expected because all array elements have direct access to the focal point. Because of the nearly uniform amplitude distribution over the feed array elements, the array excitation efficiency for this case is around 89.95 %. 44

* Figures 3.4 and 3.5 show a two dimensional intensity pattern and a focal plane intensity cut, respectively. Both figures indicate that, at the control point location, the specified intensity level is exactly fulfilled. This is an important feature of the PI synthesis and was discussed briefly in [7, 13, 15, 26, 70]. They also indicate that the synthesis pattern algorithm tends to concentrate (deposit) the power at the control location keeping minimal power deposition elsewhere. * Figure 3.6; shows a 25 % contour plot in the intensity pattern. It clearly shows that the desired focus is precisely fulfilled in both magnitude and location. It also shows that the synthesis process yielded significantly low interference patterns outside the focus location. This also justifies the fact that the PI pattern synthesis deposits the energy at the prespecified control locations 3.3.2 Multiple Focus Pattern Synthesis Having demonstrated the idea of single focus pattern synthesis, multiple focus pattern synthesis is now considered to demonstrate the capability of the PI technique to generate efficient multiple focus schemes. In this case, a set of control locations is specified along with the corresponding complex velocity potentials. The following multiple focus examples are considered and simulated: 1. Two foci synthesis: Both foci have unity intensity and are located at (60 mm, 165 mm) and (180 mm, 135 mm) or (20 A, 55 A) and (60 A, 45 A), respectively. The simulation results are displayed in Figures 3.7 to 3.10. From these figures, the following 45

0) E 0.01 <: 0.0t ar 0 50 100 150 200 Array Elements 250 4 | u)2 co (U 03 T-.0._ rr2.i ~ l lI -4 -- L I 0 50 100 150 200 250 Array Elements Figure 3.3: Complex excitation vector of the feed array (single focus). 46

1 0.8, 0) a 0.6, 4 -c la N X 0.4-. z 0.2 250 ~. ~....... ---— — --- --- ---- ~~ ~~~ ~~ 4 -~ ':............ s..r..- `5 '' 200 300 150 200 100 150 50 100 50 x-axis (mm) y-axis(mm) Figure 3y-.4: Two dimensional intensity pattern (single focus). Figure 3.4' Two dimensional intensity pattern (single focus). 47

1 - - 0.9 0.8 0.7 i e 0.5 N o 0.4 0.3 0.2 0.1 -............. '.. -.............. I.. -.............. -.................. -.................... -................. -.. '..............1............................... I...............................................................,............................................................................................................................................................................................................................................................................................................................................................................................................................................................................. 0t 50 0 50 100 150 200 x-axis (mm) Figure 3.5: Focal plane cut in the intensity pattern (single focus). 250 48

9An.1 220 200 180 160 140 E.x 120 x 100 80 60 40 20.................................................................................................................................................................................................................................... 50 100 150 200 250 y-axis (mm) Figure 3.6: Contour plot based on 25% contours in the intensity pattern (single focus). 49

can be deduced: * The complex excitation vector of the feed array shown in Figure 3.7 indicates that the uniform nature of the feed amplitude over array elements no longer exists since the feed array has to deposit specific power levels simultaneously at two different locations. This non uniformity in the feed amplitude distribution can be deduced by observing the degradation in the array efficiency from 89.95 % for the single focus case to 45.72 % for the two foci one. * Figures 3.8 and 3.9 show a two dimensional intensity pattern and focal plane cut in this pattern. Both figures clearly demonstrate that the two foci are recovered with a great accuracy in both level and location. This is also apparent from the 25 % contour plot shown in Figures 3.10. * Figures 3.8, 3.9 and 3.10 still prove that well behaved interference patterns are generated within the computational domain. However, compared with the single focus example, a degradation in the interference patterns is observed in the two foci scheme. 2. Four foci synthesis: The four foci have also equal unity intensity and are located at (60 mm, 120 mm), (60 mm, 150 mm), (180 mm, 120 mm) and (180 mm, 150 mm) or (20 A, 40 A), (20 A, 50 A), (60 A, 40 A) and (60 A, 50 A), respectively. The results of the four foci numerical simulation are presented in Figures 3.11 to 50

3.14. The following can be concluded after analyzing the results presented in these figures: 9 The complex excitation vector of the feed array shown in Figure 3.11 indicate that the dynamic range of the excitation vector increases compared to the previous example concurrently with significant degradation in the array efficiency which is approximately 27.26 % for this example. The increase of the dynamic range is related to the increased number of foci and the need to deposit more energy at more locations. Regarding the degradation in the array excitation efficiency, it was concluded before that increasing the number of foci leads to more deterioration in the array amplitude distribution. * Figure 3.12 and 3.13 show a two dimensional intensity pattern and focal plane intensity cut. Both Figures indicate that still four clear foci can be observed in the homogeneous medium. However, there is a significant deterioration in the interference patterns inside the homogeneous region or at the focal plane. * Figure 3.14 which shows a 25 % contour plot also illustrates the deterioration in the interference patterns compared with the single and double focus examples. 3. Six foci synthesis: The six foci have unity intensity and are located at (60 mm, 120 mm), (60 mm, 150 mm), (120 mm, 120 mm), (120 mm, 150 mm), (180 mm, 120 mm) and 51

(180 mmn 150 rnm) or (20 A, 40 A), (20 A, 50 A), (40 A, 40 A), (40 Aq 50 A), (60 A, 40 A) and (60 A, 50 A), respectively. The simulation results are displayed in Figures 3415 to 3.18 and the following can be deduced: * The complex excitation vector of the feed array shown in Figure 3.15, indicating that higher amplitude and phase variations are observed over the aperture plane. The array efficiency becomes 20.98 % for the six foci case. * The two dimensional intensity pattern shown in Figure 3.16 and the focal plane intensity cut shown in Figure 3.17 indicate that higher interference patterns are observed compared with the previous simulations. The degradation in the interference patterns is expected due to the increase in number of foci. However, from these two figures, still the foci are recovered with excellent precision in both level and location. * The 25 % contour plot shown in Figure 3.18 indicates a significant degradation in the interference patterns compared with the previous cases. 3.4 Discussion and Conclusions This chapter presented the computational and synthesis formulation for the problem of focusing in homogeneous media. The simplest feed array was considered and the pattern was synthesized using the PI approach presented in [5, 15, 26, 27]. After constructing the field maps in the lower half space of the feed array and according to 52

V. 0 o. 0.15 E.0.1 0.05 0 0 50 100 150 200 250 Array Elements 4 2 "CD -2 -4 0 50 100 150 200 250 Array Elements Figure 3.7: Complex excitation vector of the feed array (two foci). 53

1.4 1.2. 0.8. =... 0.6-. o 0.4 0.2 -o0 250 '. -. ".-....... n:..:. 2: 5. 200 300 150 200 100 150 50 100 50 x-axis(mm) y-(mm) Figure 3.8y- Two dimensional intensity pattern (two foci). Figure 3.8' Two dimensional intensity pattern (two foci). 54

a. n - o z 1.4 1.2 1 0.8 0.6 0.4 0.2............................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................... 0 I- `L:) 50 100 150 200 250 x-axis (mm) 1.4 1.2.1 0.0 -c 0.8 - 0.6 o 0.4 z 0.2 A,,,,I........................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................ n.; u - 0 50 100 150 200 250 x-axis (mm) Figure 3.9: Focal plane cuts in the intensity pattern (two foci). 55

vJll[ I I 220 200 180 160 -140 E E - 120 x 100 80 60 40 20.................................................................. ~........................................................................................................! 50 100 150 200 250 y-axis (mm) Figure 3.10: Contour plot based on 25% contours in the intensity pattern (two foci). 56

tU) V.t Q) 4 -V. 0.3 E >0.2 0.1 Qo, 0 50 100 150 200 250 Array Elements 4 1) 2 CD > 0 (d -2 -3 0 50 100 150 200 250 Array Elements Figure 3.11: Complex excitation vector of the feed array (four foci). 57

1.4........:. 1.2 1............ 0.8 x-axis (mm) y-axis (mm) Figure 3.12: Two dimensional intensity pattern (four foci). - 58 50

1 a I co 0 z v 0 10 20 30 40 50 60 70 x-axis (mm) Figure 3.13: Focal plane cuts in the intensity pattern (four foci). 80 59

'OA - f4 r I I I 2 2 0 -....................................................................... -........................ 180............................................................................ 200 1 60....-..........160 g 1400 - *-...........................................................................................................................40 E E 2 0 L - * - - - * * * * * - - * - - - * - -..- -. - - -............... -.......:..50 100 200 250 80 4 0......................................................................................... 20............................................................ 50 100 150 200 250 y-axis (mm) Figure 3.14: Contour plot based on 25% contours in the intensity pattern (four foci). 60

0.7 l 0............. 0.2 0.1 0 0 50 100 150 200 250 Array Elements 2................................................. 2 - O0 100 150 200 250 Array Elements Figure 3.15: Complex excitation vector of the feed array (six foci). e-2 -4 61

1.4 1.2 C c0.6 o 0.4 0.2 0O 250............ q..................... * fi...........................................: * * *.: * * * *.: * * * *:................................ *... *.. *........................................ fi god 200 150 100 200 150 50 100 50 x-axis (mm) 0 y-ax(mm) Figure 3.16: Two dimensional intensity pattern (six foci). Figure 3.16: Two dimensional intensity pattern (six foci). 62

1.4 1.2................................................................................ - c l * '* 1' '' '' ',0................. 1.. 1 CD NO0-6. - 1 I A zO 4t l.......... I............. -.- ----.. — l 0.2 - f1........................ l\.............................. co4. z 0.2........ 0 0 50 100 150 200 250 x-axis (mm) 1, > 0.8........................... -.. -.. -- - 0) c 0.6 -.... l z0.2.... 0 0 50 100 150 200 250 x-axis (mm) Figure 3.17: Focal plane cuts in the intensity pattern (six foci). 63

240 220 200 180 160 140 E E.- 120 I x 100 80 60 40 20 NC=93 — - 4................................................... i....................................... ^^ ~ 4, 94, 00 c. 0 ^:.._..................^.. -..-.................. o 0................................................ 9............ I.........................'0.......... ~..,...................................!- I I m 50 100 150 y-axis (mm) Figure 3.18: Contour plot based on 25% contours tern (six foci). 200 250 in the intensity pat 64

the data presented in Figures 3.3 to 3.18, the following may be concluded: 1. In all of the presented examples, well behaved interference patterns were observed in the region of interest. This is due to the fact that the PI technique deposits the power at the prespecified foci locations with minimal power deposition elsewhere. 2. For both single and multiple focusing pattern synthesis, the complex pressures or the intensity levels were achieved with an excellent accuracy in both value and location. This is also an essential feature of the PI algorithm. 3. As the number of foci increases, array excitation efficiency defined by (2.19) tends to decrease due to the high amplitude fluctuations over the feed array elements. Also, higher interference patterns are observed since the feed array deposits certain intensity levels at different locations. The corresponding values of these parameters as a function of the number of foci are tabulated in Table 3.1. 4. With the increase in the number of foci, the total power at the surface of the phased array increases. This can be deduced by calculating the power level at the array surface in each case. 65

Number of Foci Array Excitation Efficiency Control Point Gain M r A% G(m)/G(l) 1 89.95 1 2 46.33 0.897 4 27.26 1.131 6 20.98 1.143 Table 3.1: Comparison between different focusing schemes. 66

CHAPTER 4 Hybrid Ray Physical Optics (HRPO) Modeling After presenting the application of phased array pattern synthesis techniques and their utilization in generating highly localized fields in homogeneous media, there is a need to include the ribs which are considered strongly scattering inhomogeneities. This complicates the computational as well as the synthesis processes. To address this problem, the computational formulation is first carried out. Then, phased array pattern synthesis techniques can be invoked to compute the feed array. As stated in Chapter 1, high frequency ray techniques [17, 22, 39, 64] are employed for modeling the propagation mechanism. For computing ultrasonic fields at any location within the computational domain, a ray method is employed to propagate the ultrasonic fields from the feed array to the frontal plane of the rib cage. Then, the induced sources due to these fields are calculated along the intercostal spacings between the ribs. These sources are considered as the secondary sources for the region interior to the rib cage. Subsequently, the Physical Optics (PO) integration method [17, 64, 71] is employed to predict the 67

ultrasonic fields everywhere in the interior region of the rib cage. The feed array is computed using PI technique [15, 26, 27]. To validate the suggested approach, a real experiment was constructed to test a specific focusing scheme. Comparisons between computed and measured data were performed in order to validate the suggested computational approach. Several computational examples for single and multiple focus pattern synthesis are analyzed and discussed at the end of the chapter. 4.1 Computational Description and Formulation Figure 4.1 represents a schematic diagram that describes the idea of the computational formulation employed for predicting ultrasonic fields in the presence of the rib cage. A hybrid methodology is considered for the simulation of this problem. The key steps of the suggested approach can be deduced from Figure 4.1 and summarized in the following points: 1. A ray method is used to calculate the fields over the sub-apertures (intercostal spacings between the ribs) with the feed elements of the array being uniform line sources. The ribs are assumed solid (hard) and are considered non-penetrable by the incoming ultrasonic beams. 2. The sub-aperture fields represent the excitation to the medium below the ribs. They are integrated over the sub-apertures using the PO technique which is sufficiently accurate for large apertures. The contributions of the sub-apertures are added to give the total field at any point in the medium below the ribs. Therefore, the above computational procedure will be referred to as the Hybrid 68

Feed Array Ribs PO Technique) \ Sub-aperture Intercostal Spacing on of the propagation nechanisln Figure 4.1'.COynputational descriptio 69

Feed Array B / z r, ~ * --- —---— __z2L/ —~^ / // / A: Feed Array point;rd B: Upper Region Point. d C: Sub-aperture Point. D D: Field Point. Figure 4.2 Geometrical description of the propagation model. 70

Ray-Physical Optics (HRPO) technique. 3. The relation between the velocity potential and the array excitation vector is formulated in matrix form to employ the PI technique for pattern synthesis as indicated in Chapter 3 and in [15, 26, 27]. In formulating the problem at hand, continuous quantities in the upper region (above the rib plane) will be referred to with the subscript u and the corresponding discrete quantities with the subscript k. Over the rib plane which includes both the solid rib and the intercostal spacings, continuous and discrete quantities are referred to as a and p, respectively. Finally, for the lower region (below the ribs), continuous and discrete quantities will be referred to as d and m, respectively. Figure 4.2 demonstrates the geometrical distances and locations for the different field points that will be used in the computational formulation. Starting from (2.8), the velocity potential for a source radiating in homogeneous media is given by ru= juS(U)gfu(ruIru')dS', (4.1) where (rFu) is the radiated velocity potential in m2/sec at location ru from a unity source with normal velocity distribution us(ru') located at rFu', gfu(rFujru') is the appropriate Green's function and S' is the feed aperture surface area. Typically, the feed array is placed approximately at 5 cm above the rib plane. This distance is equivalent to 17 A at 500 kHz. Since the feed array is far away from the rib cage, -the Green's function can be approximated by the half space one. As shown in Figure 4.2, the radiated fields will illuminate the tissue spacings between the solid ribs and hence induce secondary sources within the intercostal 71

spacings. The normal velocity distribution over the aperture rib plane is given by Ua(r) = n V ) (4.2) where ~(ra ') is the velocity potential at r ' located within the sub-apertures or intercostal spacings between the ribs and n' is the unit normal to the sub-aperture surface. Substituting from (4.1) into (4.2), the normal velocity distribution can be written as Ua(r' = [V' us(r gfu(, )dS') (4.3) The sources, given in (4.3), excite the ultrasonic field in the lower region. The total velocity potential in the region interior to the rib cage due to these sources is calculated by using the aperture velocity distribution as the secondary source in the first term of (2.8) to give 4d(rd)= Ua(ra )gfd(FdIra')dSa, (4.4) where rid is the field point location. It should be noted that (4.4) is in the context of the PO approximation. Replacing Ua(ra') in (4.4) with its value from (4.3) results in bd(rd) = j [n (s( )g fu( gIrfru'))]gfd(rdIra')dS'dSa- (4.5) Discretizing the feed aperture to K elements, each with uniform field (feed) distribution, the integration in (4.5) will be transformed to a summation over all the source segments. Thus, the discretized version of (4.5) can be written as K m = ukej k H(m,), (4.6) k=l where kek is the complex excitation velocity of of t kth element of the discretized source which is the linear array in this case. As described in Chapter 3, H(m, Vk is 72

the effect of the kth element at the mth field point under the solid ribs and its general element H(m, k) is given by H(m, k) = j n * V'[gu(aIu')gfd( dIa')] dSj (4.7) In matrix form, (4.6) can be rewritten in the following compact form ad = Hu, (4.8) where uu is the complex excitation vector of the feed array with length K. As described in Chapter 3, the first step towards the patterns synthesis is performed by finding the excitation vector u, that generates the desired field pattern. The synthesis procedure will be described in the following section. 4.2 Computational Example The problem of ultrasound field propagation through the intercostal spacings between the solid ribs and the surrounding tissue was approximated by the two dimensional model shown in Figure 4.3. As indicated in this figure, the rib cage was modeled as a planar strip array located at the interface of the bolus and tissue (simulated as a stratified medium). The feed array is located on top of the bolus (water) and although shown planar in Figure 4.3, this is not required in the formulation. The parameters of interest are: * K: total number of the elements in the linear array (feed array). * d: inter-element separation between array elements. 73

1 2 d t le et k element K y Figure 4.3: Two dimensional geometrical modeling of the problem. 74

* Uk amplitude (driving amplitude) of the kth element of the array. * Ok phase shift (driving phase) of the kth element of the array. * Yk height of the kth element of the feed array over the rib plane. * N: number of solid (hard) ribs, each of width 11. * 1: separation (intercostal spacing) between ribs. It is often referred to as the sub-aperture width. Additional parameters are required to define the stratified medium below the ribs (fat, viscera and liver tissues). However, as seen from Table 2.1, the specific acoustic impedance of these layers are so close that the medium is nearly homogeneous (except for the solid ribs). Nevertheless, this assumption is made for simplicity and does not represent a restriction on the employed analysis. The HRPO approach for the problem shown in Figure 4.3 will be demonstrated in accordance to the formulation presented in Section 4.1. The velocity potential at the general point x' on the nth sub-aperture (intercostal spacing) due to the feeding of the kth element of the array can be formulated by including the effect of the complex excitation factor Ukejok in (3.7) and rewrite it as [86], jo -,rk n =s(rkn -)=-/Ouke - k F(O), (4.9) where 41s(rkn) is the radiated velocity potential from a certain line source and rkn is the distance between the kth source and the general point (x', 0) on the nth subaperture. To calculate the total field at the general point (x', 0) due to all K elements 75

of the array, array element contributions are superimposed to give K a(X') =- Z ls(rkn), (4.10) k=l where 4a(X') is the overall or total velocity potential at a general point (x', 0O) on the nth sub-aperture. The potential due to the nth sub-aperture at any point rn inside the rib cage can be obtained by using the induced aperture normal velocity distribution as the excitation in (4.4). Namely, (rn) = ua(x')gfd(rnir')dS (4.11) where Ua(x') is the normal velocity at a general point (x', 0) on the nth sub-aperture and gfd(rnIr ) is the appropriate Green's function for the lower half space. As indicated by (4.2), ua(x') can be obtained from the following relation Ua(X ) = n'. V' a(X'). (4.12) For the half space below the rib surfaces, Green's function gfd(rnIrn ) can be expressed in the form [17, 39, 86], gfd(rnr) = ((- ' ). (4.13) To get the expression of the velocity potential at the general point rn, (4.12) and (4.13) are utilized to replace the secondary sources and the Green's function in (4.11), respectively. Thus, 4(rn) can be rewritten as K(rn) = -2 In'. V' +(x')]H )(-jr')dx'. (4.4) It can be seen that (4.14) represents the radiated fields based on the PO approximate method [17, 64, 71]. Substituting from (4.9) in (4.14), the velocity potential at 76

a certain location inside the rib cage can be expressed as 2(rn) = uke ( H )(-r)dx. (4.15) Rearranging (4.15) and noting the independence between the linear integration over the sub-aperture and the summation over the feed array elements, the radiated velocity potential from the nth sub-aperture at the field point can be expressed as — 'Y^rkn <(rn) = -j Eukek Ho2)(-jhyr )dx'. (4.16) 2 Lk= k^1 Ix'' V/7'^kn By imposing the superposition principle, the total velocity potential at the specified field point due to all sub-apertures will be N Obtotal(r) = Z. r(r), (4.17) n=l where Vtotal(r) is the total velocity potential at the field point r. Substituting from (4.16) into (4.17), 4total(r) will be N K J, rk /total(r) = 2 uk H 2)(- j'yr' )dx'. (4.18) 2 n=l k==l Jl rkn Because of the independence between the two summations presented in (4.18), the total velocity potential can be expressed as K A total(r) = 2 N Ho2)(-jr')dx' (4.19) 2 (4.19) tk=1 n=1 J2l k- n To simplify the analysis, (4.19) can be reduced to the following convenient form. K N total(r) = 7' E uke^ E hT(k, n, r), (4.20) k=l n=l where rl e-v/Q+S(') (2) hT(k, n, r)= / H42)(-jr')d', (4.21) Jy/Q + S(x') 77

S(x') = (X')2 + 2x'(1 + I)(n - 1) + 2x'd(k - 1) + x'l, (4.22) Q = y + x4 + (n - 1)2(1 + 1)2 - 2Xk(n - 1)(1 + 11), (4.23) and r' = y2 + ( - -(n -1)(+ ))2 (4.24) in which x and y are the coordinates of the field point. In order to apply the pattern synthesis technique, (4.19) is expressed in matrix form. To accomplish this, (4.20) is tested at a general field point with Xm and ym coordinates. The total field at this point 1total(rm) is achieved by replacing r by rm in (4.20), viz K total(rm) = EukeJOH(m k), (4.25) k=1 where N H(m, k) = Z hT(k, n, m). (4.26) n=1 In this form, H(m, k) can be identified as the contribution of all sub-apertures at the mth field point due to the kth line source. Testing (4.25) M points gives the following matrix system 4ctotal =H u, (4.27) where uu is the complex excitation vector of the feed array with length K with its kth entry equal to UkejOk and H is the total matrix response (of size M x K) where its 78

entries are given by (4.26). Also, ctota, is a vector of length M, each of its elements represents the velocity potential at the control (field) point rm, m = 1, 2,...., M. For HIFU treatment, it is required to control the field or intensity levels at a certain number of control points in the interior region of the rib cage. Typically, the number of array elements K is significantly higher than that of the control points M, that is K >> M. Therefore, (4.27) represents an underdetermined system to be solved for the desired response u,. Following the procedure in [15, 26], (4.27) can be solved to obtain Uu = HPiIJctotal, (4.28) where HPi is the pseudo-inverse of H given by (3.6). The PI pattern synthesis technique is briefly discussed in Appendix A. 4.3 Rib Fields After formulating the computational and synthesis methods, field levels over the rib surfaces are computed. The computational model for the rib field assumes the following: 1. The ribs are of infinite depth (two dimensional) and this is shown by the model presented in Figure 4.3. 2. The incident waves at the rib locations are locally planar. 3. As displayed in Table 2.1, the specific acoustic impedance of the hard ribs is significantly higher than that of water, thus unity reflection coefficient may be 79

assumed as the worst case scenario. Based on these assumptions, the total velocity potential at the general point x' over a certain rib along the rib-tissue plane can be expressed as [17, 44, 86], K 1r(X') = E uke z kgfu(r*sFrk)dS' (4.29) k=l Sr where S' is the surface area of the solid rib and rks is the vector connecting the kth element of the feed array and the sth point over the rib surface. In matrix form, (4.29) can be written in the following way cr= Hr Uu, (4.30) in which the general element Hr(s, k) of the matrix Hr can be expressed according to the following relation Hr(s k) = gfu(krk)dSr. (4.31) After the feed array is computed using the suggested approach, the ratio between the power incident upon the rib surfaces and that reaching the frontal plane of the rib cage is calculated. This ratio gives a measure for the power intercepted and blocked by the solid ribs to that coming from the feed array. The lower this ratio is, the more power reaches the interior locations of the rib cage. This ratio can be defined in the following way r p x 100, (4.32) Pr + PS where Er is referred to as the relative rib power ratio, Pr and PS are the total power reaching the ribs and intercostal spacings, respectively. Pr and Ps can be calculated 80

using the basic relations described in Chapter 2, namely N,r PrIi=(rr )dS'r (4.33) ir=1 -Jihrib and Ns '= E I (rs') dS;, (4.34) i -=l ish1 spacing where Nr and Ns are the number of ribs and intercostal spacings, respectively, rr ' and rs ' represent two general points on the rib and intercostal spacing, respectively. The quantities Ii7 and I, represent the intensities over the ribs and intercostal spacings, respectively. Also, dSr' and dSs' represent the elemental areas for the ribs and intercostal spacings, respectively. The model shown in Figure 4.3 is utilized to examine the performance of the suggested HRPO approach by computing the fields over the rib surfaces. Based on (4.30) and (4.31), the total surface velocity potential at the general point x' along the rib-tissue plane is given by K e-yr Or(I) = -O uke' (4.35) k=l1 where rk' is the distance between the general point (x', O) on the rib-tissue plane and the kth element of the array. It is given by rk = yVy, 2 (4.36) with x1 = (n - 1)(l + 11) + (1/2) - (k - 1)d + x'. (4.37) 81

The corresponding expression in matrix form can be deduced in the same way as (4.30), where Hr is a vector of length K, each of its entries is given by e-a-/rks Hr(s,k)= -o =, (4.38) where Hr(s, k) represents the effect of the kth line source of the feed array at the s th general rib point. Using (4.38), the corresponding intensity pattern over the rib surfaces can be evaluated to examine whether rib levels exceed the threshold required for safe treatments. After achieving the intensity levels using (4.33) and (4.34), the rib power ratio Er can be evaluated from (4.32). 4.4 Experimental Validation After formulating both field computations and pattern synthesis approaches required for HIFU treatment in the presence of the rib cage, it is required next to validate the HRPO suggested approach. An efficient way to achieve this validation is by constructing a real experiment. After measuring the fields experimentally, the HRPO technique is applied to the same set up in order to compare the measured data with the computed one. The experimental setup is shown in Figure 4.4 and this set up was used to focus through a strip array represented by 6 wooden bars each 18 mm wide and separated by 19 mm. The homogeneous medium is simulated by water and the goal is to deeply focus the HIFU beams through the wooden bars to a specific location inside the water tank. As displayed in Figure 4.5, four flat subpanels, each having 16 elements are used to construct the feed array. The element width was 3 mm and the elements were separated by 3.5 mm within each panel. Also, the panel 82

Figure 4.4: Experimental setup. 83

30 -150 20 -100 -10 50 100 -20 Z(mm) X (mm) 150 -30 Figure 4.5: Two dimensional geometry of the array utilized in the experiment. The wooden bars simulating the ribs are positioned in the y=O plane along the x-axis with their lengths (depths) running parallel to the z-axis. 84

width was 61.5 mm and the inter-angles between the adjacent subpanels were 162.5~ and 145~ for the middle and end subpanels, respectively. The panel depth was 60 mm which is large enough for the 2-D approximation considered here. The excitation frequency was 485 kHz and the feed array was adjusted to focus at a location which is 10 cm far from the wooden bars. The experimental results are obtained by measuring the intensity levels at the focal plane. The HRPO approach was applied to the same geometry and focal plane fields are compared to the experimental data. The experimental and computational results are given in Figure 4.6. Clearly, a good agreement exists between the calculated and measured heating patterns especially in the region of interest which includes the focus location. This focal plane cut is used to validate the proposed approach. As observed in this figure, the measured and computed levels of the main lobe are almost identical and the desired focus is accurately fulfilled in both magnitude and location. On the other hand, away from the main lobe, the measured data deviate from the calculated due to three reasons. The first is the uncertainty in the bar locations and dimensions (the associated resolution was around 1 mm which equals one third of a wavelength). Second, multiple reflections from the boundaries of the water tank caused a standing wave which affected the region outside the main lobe of the focal plane. Third, the lack of accurate alignment between the feed array and the wooden bars and this inaccuracy resulted in phasing errors which slightly deteriorate the focusing scheme away from the focus. Despite the experimental limitations given above, the results presented in this section support the use of the hybrid model in predicting the quality of phased array 85

field patterns in the presence of strong obstacles. That is, the hybrid model offers a computationally efficient tool for analysis and synthesis of phased array heating patterns for localized treatment of deep targets in the presence of strongly scattering obstacles. 4.5 The HRPO Algorithm According to the aforementioned field computational approach and based on the PI patterns synthesis techniques developed and implemented in [26, 27], feed array design algorithm using the HRPO method may be summarized as follows: 1. First, define the field control points (M points) and the corresponding complex pressure or intensity values at these points (heating levels). The respective field points are usually the tumor location and points at which minimum power deposition must be kept as low as possible. 2. For each point, evaluate the response vector using (4.19) to (4.26) to construct the system response matrix with size M x K. 3. Using the pseudo-inverse technique and following the procedure discussed in Chapter 3 and this chapter, the complex excitation vector of the feed array uu is computed. After computing the feed array, the region below the ribs is scanned to obtain field intensity map to ensure the fulfillment of the desired intensity levels at the specified control locations. Finally, the evaluation of the heating levels over the solid ribs is 86

carried out to ensure that the employed pattern synthesis approach has maintained relatively low levels over the rib surfaces. Otherwise, step one must be revisited to redefine the M control points by including the overheated rib locations and repeat steps 1, 2 and 3 of this algorithm. 4.6 Numerical Examples To examine the performance of the suggested HRPO technique, both single and multiple focus pattern synthesis are carried out. The simulations employed a typical biological medium associated with a sound speed of 1500 m/sec and a density of 1000 kg/m3. The operating frequency was set to 500 kHz (3 mm wavelength) unless otherwise noted. The geometrical features of the feed array and tissue parameters are given by * 240 planar, uniformly spaced excitation elements with interelement separation equal to half a wavelength. * The array elements heights from the rib plane is 50 mm which is approximately to 17A. * Number of solid ribs is 6. * Rib width is approximately 18 mm or 6 A. * Sub-aperture (tissue between ribs) width is 21 mm or 7A. 87

4.6.1 Single Focus Scheme To test the suggested HRPO approach, a single focus pattern with a unity intensity focus located at (120 mm, 150 mm) or (40 A, 50 A) is synthesized. Typically, for therapeutic phased array focusing, the ratio between the focal spot distance from the array surface and the maximum aperture width (f-number) should be kept around unity. This restriction controls the choice of the phased array distance from the rib plane. More details about the spot size and the 6-db beamwidth are given in [50]. The design goal is to compute the complex excitation vector of the feed array uu by supplying the focal specifications to the HRPO synthesis algorithm. Then, field maps inside the rib cage and over the rib-tissue plane are obtained using (4.27) and (4.38), respectively. As demonstrated in [26], if the field control points are distant enough from each other, then, the matrix H becomes well conditioned and thus PI can be directly applied to compute the complex excitation vector of the feed array. Also, in calculating the integral in (4.21), fourth order Gaussian integration scheme was employed [1, 19], (see Appendix B). Because of the highly oscillating integrand, the main integral was subdivided until convergence was achieved. The numerical results are illustrated in Figures 4.7 to 4.11. Specifically, Figure 4.7 shows the complex excitation vector of the feed array with its relative amplitudes and phases as computed by the PI method [5, 7, 15, 26, 27]. The amplitude distribution demonstrates the performance of the employed PI technique since array elements which have insignificant effect on the control points are minimally fed. These array elements are those facing the rib obstacles. Thus, excessive amount of heating in the 88

undesired tissues is avoided. Different field maps for the synthesized single focus pattern are shown in Figures 4.8 to 4.10. These maps were constructed using one wavelength scanning step in both the x and y directions. As shown in Figure 4.8, ultrasound energy is almost concentrated at the focus location with minimal interference patterns outside the focal plane. Figure 4.9 displays a contour map based on 25 % contours for the field intensity pattern inside the rib cage. Figures 4.8 and 4.9 show that ultrasound energy passes through the intercostal spacings or sub-apertures between the solid ribs and converges towards the focus location. As observed in the focal plane cut shown in Figure 4.10, the intensity level at the focus location is precisely fulfilled. This is because the PI algorithm achieves the exact power deposition at the control point locations(s) and minimizes power deposition elsewhere. Finally, Figure 4.11 shows the normalized fields (with respect to the maximum level at the focus location) over the rib plane. It is indeed pleasing to observe that the intensity levels over the ribs and the intercostal spacings are vanishingly small compared with that at the focus. As shown in the figure, the maximum value for the normalized intensity over the rib-tissue plane is less than 0.06. This directly proves that the pattern synthesis algorithm used in the analysis of this model passes most of the ultrasound energy through the intercostal spacings between the solid ribs. Moreover, the rib power ratio defined in (4.32) is negligible (less than 3 % in this case). As a result, the suggested HRPO algorithm eliminates concerns relating to the pain level on the rib/bone surfaces. 89

4.6.2 Multiple Focusing Scheme The proposed HRPO is also tested for synthesizing multiple focus patterns. In this example, focusing was simultaneously performed at two different locations inside the rib cage. These locations are given by the coordinates (60 mm, 165 mm) which is equivalent to (20 A, 55 A) and (180 mm, 135 mm) which is (60 A, 45 A). The results are displayed in Figures 4.12 to 4.16. As observed in Figure 4.12, the amplitude distribution of the complex excitation vector of the feed array has the oscillatory nature observed when single focus was tested. The array elements with low power contribution at the foci locations are minimally fed while those with direct access to the foci locations are excited with high power or amplitude. The two dimensional intensity pattern shown in Figure 4.13 indicates that ultrasonic beams propagate through the intercostal spacings between the ribs to add constructively at the prespecified locations. The contour plot shown in Figure 4.15 shows high field localization at the foci locations coupled with minimal interference patterns in the interior region of the rib cage. Figure 4.14 shows that the intensity levels at the foci locations are obtained precisely while Figure 4.16 shows the HRPO synthesis technique yields minimal intensity levels over the solid ribs. From Figures 4.13 and 4.15 there is a higher interference pattern observed in the near bone region compared to the single focus case. This is because the PO technique fails to predict the near zone field of a source or radiator. 90

- Measured Data \ 0.9 x. Computed Data....................... 0.8........................................... 0.8 -......;......................................................... - 0.7. 6..............1 1...... g 0.6 0.5........ | 0.4...... 0.3............. 0.2.. --- - - I.......................x.. 0.1. /xx x - x XX X: V x X x x ' X Xx X 0 X II XI X Ix 0 5 10 15 20 25 30 35 40 45 x-axis Figure 4.6: Comparison between measured and computed data. 91

0.7 0.6 (c (0 v 0.5 0.4 E 0.3.' 0.2 aT 0.1 0 C - - -~IIP-.N....... I 25 200 250 I 50 100 150 Array Elements (1 a) If) (0 cn'. A. --- - ' 2 OL - 0 50 100 150 200 250 Array Elements Figure 4.7: Complex excitation vector of the feed array shown in Figure 4.3 (single focus). 92

C -o~ 0.... 0.... *.:........ ------: 0 0 250............................................................................. *............ N....... no -...... K........ v.... z t * ---s:: * —.: 5..s.o | so-vs = -- 250 200 300 150 200 100 150 50 100 x-axis (mm) 0 0 y-axis (mm) Figure 4.8: Two dimensional Intensity pattern for the Single Focus Case. The simulated ribs are centered at y=0 and x=30 mm, 69 mm, 108 mm, 147 mm, 186 mm and 225 mm. 93

0) o O Z z 50 100 150 200 x-axis (mm) Figure 4.9: Focal plane intensity cut at y=150 mm. 250 94

va~l_ I II 220 200 180 160 -140 E E.-~ 120 x 100 80 60 40 20...................................................................................................................................................I.......... I..,...................................................................................................................................... I................................................................................................................................ II I I 50 100 150 y-axis (mm) Figure 4.10: Intensity pattern contour plot (25 focus). 200 250 % contours) (single 95

u) c-.c ' N E z. z 0 50 100 150 200 2 Rib-Tissue axis (mm) Figure 4.11: Intensity pattern over the rib tissue plane (single focus). 50 96

() V.O U) E cD 0- 0.2 0 0 50 100 150 200 250 Array Elements Array Elements 4 ure 4.3 (two foci). -c -2 -4 0 50 100 150 200 250 Array Elements Figure 4.12: Complex excitation vector of the feed array shown in Figure 4.3 (two foci). 97

1.4 1.2 -z 1 ~ 0.8 -a) ~ 0.6 -z 0.4 0.2 O 250 ~ ' ~ '..... 4 II I.................... *........... t.......:-...: *.,:................ *.....: *.::.,:........ *..........................................:.: '.,:.,: l......,......,..........,. *: 1;s:.: ';.:... L.... -............. 300 50 x-axis (mm) 0 y-axis (mm) Figure 4.13: Two dimensional intensity pattern (two foci). The simulated ribs are centered at y=0 and x=30 mm, 69 mm, 108 mm, 147 mm, 186 mm and 225 mm. 98

21 c 0.8..........................................................................- 0.6 ' -................................... * 0 i E ~ o............................................................................................ z o0.4... ~ l 0.2 0 50 100 150 200 251 x-axis (mm) 1.4 1.2.. 2.1 ). I: 0 ~ ~ aC 0.2 0 50 100 150 200 251 x-axis (mm) Figure 4.14: Focal plane intensity cuts at y=135 mm and at y=165 mm. 0 99

240 220 200 180 160 -140 E.- 120 x 100 80 60 40 20.......................................................................................................................................................................................... -..... -..................................................................................................................................................................................................! 50 100 150 200 y-axis (mm) Figure 4.15: Intensity pattern contour plot (25 % contours) 250 (two foci). 100

4.7 Comparison with the Focusing in Homogeneous M\edia Having demonstrated the single and multiple focusing schemes through the rib cage, it is required to compare these schemes with the ideal case when focusing was considered within a medium free from strongly scattering obstacles. For both single and two foci schemes, the complex excitation vectors of the feed array, contour maps and focal plane cuts are compared. 4.7.1 Single Focus Scheme In order to understand the focusing scheme limits, it is essential to compare the homogeneous and inhomogeneo~us schemes when applied to synthesize a certain pattern with a desired focal intensity at a specific location. For the single focus scheme, a unity intensity is intended at (120 mm, 150 mm) or (40 A, 50 A). The results for both cases are compared and displayed in Figures 4.17 to 4.22 and from these figures, the following can be concluded: * The complex excitation vector of the feed array for the homogeneous scheme tends to be more uniform than that of the inhomogeneous case. This is apparent when comparing Figures 4.17 and 4.18. This is due to the fact that the elements of the feed array in the homogeneous case contribute equally to the focal point because of their direct access to the focal location. On the other hand, for the inhomogeneous scheme, the focal spot is located in the shadow region of the 101

elements that are facing the solid ribs and thus according to the reciprocity principle, these elements have weak contribution or effect at the focal location. Therefore, they are minimally fed or excited and this directly results in lower array efficiency for the inhomogeneous scheme. * The interference patterns in the inhomogeneous case tend to deteriorate when compared to the homogeneous case. Figures 4.19 and 4.20 display the contour plots for both cases. Clearly, the introduction of the solid ribs in the computational model reduces the focusing mechanism capability of producing highly localized fields inside the rib cage. This reduction is due to the decrease in the effective aperture used for focusing when the ribs are considered. The effective aperture is represented by the intercostal spacings between the ribs and the smaller this aperture becomes, the lower the aperture gain will be and thus interference pattern deterioration is expected. * The focal plane intensity distribution is better when homogeneous focusing is considered. This fact is clearly demonstrated by comparing Figures 4.21 and 4.22. In both graphs, the focal spot is achieved in both location and level with high degree of accuracy. However, interference patterns over the focal plane are significantly higher in the inhomogeneous focusing case and this is also related to the reduction of the effective aperture when the ribs are considered. 102

0 50 100 150 200 Rib-Tissue axis (mm) Figure 4.16: Intensity pattern over the rib-tissue plane (two foci). 103

250 Array Elements. inn rk '40.1 7M I — 4' 0 so 100 150 200 25 s0 Array Elements Figure 4.17: Complex excitation of the feed array for the homogeneous scheme (single focus). 0.7.. 0.6 a_ I0 8 0.5 'I 0.4 a 0.3 aD 0.2 0.1 0 Array Elements 4 a. DID 0. ID a: - 2 ll l 1 klil ktk li 1.' L,ilj ii J Jli llJJI lIIi I 0 n' -l IL._-._ -- _ -i- - -i nn - LJ iLiJlJuJlUL j 7 VI I I I I-... r r ll I l -d _,- i 0 50 100 150 200 250 Array Elements Figure 4.18: Complex excitation of the feed array for the inhomogeneous scheme using the HRPO technique (single focus). 104

24011 220 200 180 160 - 140 E.- 120 100 80...........;........................................... 60 - 40 20.. I 50 100 150 y-axis (mm) 200 250 Figure 4.19: Contour plot based on (25 % neous scheme (single focus). contours) for the homoge 240. I 220 200 180 - 160 - - 140 E.I 120 100.......... I.................I <Z>................ 80 - 60 - 40 20...I 50 100 150 y-axis (mm) 200 250 Figure 4.20: Contour plot based on (25 % contours) for the inhomogeneous scheme using the HRPO technique (single focus). 105

1 0.9 0.8 0.7 * 0.6.G 0.5 a o 0.4 z 0.3 0.2 0.1 I I I. ~........ 100 15 200 250 C'I 0 Focal 50 plane cut for 1 00 150 200 250 x-axis (mm) the homogeneous scheme (single focus). Figure 4.21: 1,,, T I I 0. 9........................................:..................... 0.89............................................................ 0.6..................:....................i.................................. 0.7 0.5............................................................................ z 0,.3................................................ 0.2.............................................................................. 0.1 0 50 100 150 200 250 x-axis (mm) Figure 4.22: Focal plane cut for the inhomogeneous scheme using the HRPO technique (single focus). 106

4.7.2 Twol Focus Scheme After comparing the homogeneous and inhomogeneous schemes for the single focus pattern synthesis, the two focus mechanism is now tested for both cases. As shown in the single focus comparison, complex excitation vector of the feed array, contour maps and focal plane cuts are compared when a two focus pattern is synthesized and tested for the homogeneous and inhomogeneous cases. The foci locations are fixed at (60 mm, 165 mm) which is equivalent to (20 A, 55 A) and (180 mm, 135 mm) which is equivalent to (60 A, 45 A). Both foci are desired to have unity intensity. The comparison results are demonstrated in Figures 4.23 to 4.28. From these figures, the following conclusions are reached: * The amplitude of the complex array excitation for the homogeneous scheme is still more uniform than that for the inhomogeneous one. Figures 4.23 and 4.24 demonstrate this conclusion. The difference between the two cases is significantly reduced compared to the single focus case. This is due to the fact that in the two foci pattern synthesis for the homogeneous focusing, the feed array deposits certain amounts of power at different locations and directions. Therefore, each array elements is driven with restricted access to each focal point. On the other hand, for the inhomogeneous focusing, the rib existence still dominates the array feed and thus array elements facing the solid ribs are fed with small power. * The interference patterns for the homogeneous scheme is dramatically better than those for the inhomogeneous one. Figures 4.25 and 4.26 show that for 107

the same test case, interference schemes are deteriorated when the ribs are considered. The main reason for this deterioration is due to the reduction of the available aperture or acoustical window when the ribs are considered. According to the inhomogeneous interference pattern shown in Figure 4.26, it is not feasible to achieve more than two foci with acceptable interference patterns when considering the ribs. * Focal plane fields are deteriorated with inhomogeneous focusing. This can be directly deduced from Figures 4.27 and 4.28. The degradation is mainly because of the loss of effective aperture area when the ribs are included in the geometry. 108

Array Elements 'W'11i'4P -4 0 50 100 150 200 250 Array Elements Figure 4.23: Complex excitation of the feed array for the homogeneous scheme (two foci). 1.. IA.0.8-.0.6 0.4 a 0.2 0 Array Elements 4, 1 —~- - 1 1 I I.ua is mI:-A 2-L 2L~.j,,III,!~~I ]jili~ 'Jr]l,!l, _4.. 0 50 100 150 200 250 Array Elements Figure 4.24: Complex excitation of the feed array for inhomogeneous scheme using the HRPO technique (two foci). 109

"An. z-v I - 220 -200 180 160 - 140 - E.- 120 -100 80 60 40 20 Figure 4.25:..........~ ~ ~ ~ ~~ ~~ ~- ~~ ~...........................................-~ ~ ~- ~ ~~-~.~-~:..........................................................~...................................................................................................................................................................................................................................; I I, I 50 100 150 y-axis (mm) Contour plot based on (25 % neous scheme (two foci). 200 250 contours) for the homoge 240 ---- I 220 200 180 160 E- 140 E. 120 1 00 80 60 40 20 -........................................................................ %-. i - * ~ <' * - - ':?5* o4 -.a........................................................................... - ' * 0 o c... >..................................................... - - -........................... K. *...................,......................,.-........................................-................................. -................. I I - - - 50 100 150 200 250 y-axis (mm) Figure 4.26: Contour plot based on (25 % contours) for the inhomogeneous scheme using the HRPO technique (two foci). 110

x-axis (mm) 0 50 100 150 200 250 x-axis (mm) Figure 4.27: Focal plane cuts for the homogeneous scheme (two focus). 1.4. 1. C 0. z O..2 8........ 50 100 150 200 2! S0 x-axis (mm) 250 x-axis (mm) Figure 4.28: Focal plane cuts for the inhomogeneous scheme using the HRPO (two foci). 111

4.8 Discussion and Conclusions In this chapter, a hybrid ray physical optics computational approach for analysis and synthesis of phased array multiple focus patterns in the presence of strong scatterers was developed, implemented and experimentally validated. This hybrid approach ignores certain diffraction terms from the rib edges and is therefore expected to be less accurate than a solution based on a full wave numerical implementation. However, this loss of accuracy is not expected to have significant effects on the validity of the results obtained in the examples shown since the neglected contributions provide correction terms that primarily improve the solution in the shadow region near an obstacle or an edge. At distances several wavelengths away from the obstacle, ray theory gives an acceptable solution and therefore the hybrid method is sufficiently accurate at foci located far away from the ribs. Furthermore, the use of the Rayleigh-Sommerfeld diffraction integral [17, 75] eliminates any difficulties relating to the prediction of fields near caustics. Therefore, theoretical foundation of the hybrid approach presented in this chapter is quite sound. It is interesting to note that Figure 4.7 and 4.12 indicate that the synthesis procedure utilizes the full aperture of the array, including elements that do not have a direct path to the focus. In contrast, previously proposed methods relying strictly on ray tracing (see [50] for an example) utilize only a fraction of the array aperture. This could severely limit the array gain and, consequently, limit the ability to generate useful heating patterns. A fully developed diffraction analysis model with uniform theories [58] or a full 112

wave numerical model will provide a better simulation of the intensity patterns in the presence of highly scattering obstacles. However, even with a more accurate model, uncertainties in the rib locations with respect to the feed array are likely to give phase errors which may offset any gains in accuracy. Both simulation and experimental data shown in the chapter clearly demonstrate that focused phased array field patterns can be generated in the presence of strongly scattering obstacles. Thus, phased array systems are capable of utilizing all acoustical windows (intercostal spacings between the ribs) to focus on targets partially shadowed by these obstacles. This is a unique feature to phased array applicator systems which cannot be achieved by other focusing systems. That is, phased arrays continue to offer the promise of non-invasive treatment of tumors and other target tissues considered inaccessible to other focusing systems, e.g., lens-focused and shell transducers. The results shown here and experimentally obtained in other investigations, for example [50], demonstrate the potential for large-aperture ultrasound arrays for tissue targets in the liver and the heart. Another point worth discussing in the simulation results is the low intensity levels observed over the solid ribs which are much less than those reaching the rib plane (less than 6 % or 7 %) despite that no rib field minimization criteria were imposed in the synthesis process. An interpretation to this phenomena can be concluded from Figure 4.29 where by invoking reciprocity, the following can be deduced: * Invoking reciprocity, the problem may be looked as a line source placed at the tumor location. 113

* With the ribs considered as solid impenetrable objects, there will be no partial transmission through them. * Thus, array elements facing the rib obstacles are minimally excited because of their location in the shadow zone of the radiating line source (see Figure 4.29). * Because of their direct access to the solid ribs, these elements provide minimal power deposition over the ribs. * In the same time, by recalling reciprocity, the feed array will produce the desired level at the tumor location. With regard to the overall performance of the HRPO modeling scheme, three major drawbacks are associated with the utilization of the HRPO for pattern synthesis and field computations. They can be summarized in the following points: 1. The presence of strongly scattering obstacles calls for a modification of the basic synthesis algorithm. This modification involves a rib power minimization criteria in order to ensure minimal power deposition over the solid obstacles. 2. The use of the physical optics PO integration scheme for computing ultrasonic beams in the interior of the rib cage introduced some limitations due to the inaccuracy of the PO in predicting the fields in the near zone of apertures or at grazing angles. 3. For field maps reconstruction after computing the feed array, the HRPO approach is time consuming because of the fine discretization of the PO integration for each field point inside the rib cage. 114

To overcome these drawbacks, a two-step synthesis procedure is presented and tested for this purpose in the next chapter. 115

IL I Feed Array Plane I I I I I I I I I I I I I I I j/ H L I --- —- ~~. —,Shadc,,v! I I I I I I I I I II I I I I I.I, I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I ) I 'W Regions II I I '! i: R / w Solid I Ribs w Source L -- Lower Power H — Higher Power Figure 4.29: Illustration of the reciprocity concept utilized in interpreting the HRPO results. 116

CHAPTER 5 Hybrid Two-Step Virtual Array Ray (VAR) Modeling Although the HRPO technique developed in Chapter 4 deposits acceptable power levels over the solid ribs, it lacks the explicit control of intensity levels over these obstacles. Also, because of the employed PO approximation, HRPO did not accurately predict the fields in the near bone vicinity. Moreover, from a computational point of view, field maps reconstruction after computing the feed array was an extremely time consuming task. To overcome these drawbacks, a two step hybrid computational method is developed, implemented and validated in this chapter. In the first step, a virtual array is introduced over the intercostal spacings between the solid ribs to generate the prespecified intensity levels at a set of control points within the target region. The second step is the design of the actual feed array that induces the required virtual sources along the intercostal spacings. The rib power minimization and design optimization are carried out via the PI technique and by a constrained preconditioned PI method. The VAR method computes the required 117

primary sources while maintaining minimal power deposition over the solid obstacles. Several examples representing single and multiple focusing schemes are presented in this chapter to provide a brief comparison with the HRPO approach developed and implemented in Chapter 4. 5.1 Computational Description and Formulation The first step of the suggested technique involves a pattern synthesis procedure performed in the half-space beneath the ribs. A set of fictitious line sources separated by A spacings is introduced along the intercostal spacings between the ribs. These sources constitute a virtual array (VA) covering the available acoustical windows into the target region. If the first step does not yield a useful pattern, then no array configuration can be designed to achieve the treatment objective. In other words, the first step determines if the available acoustical window puts a fundamental limit on the ability to focus in the presence of a given obstacle. The second step is another synthesis problem with the fictitious sources within the intercostal spaces defined as control points and the real array elements defined as sources. Standard ray frequency techniques are employed to propagate the fields from the feed array to the intercostal spacings between the solid ribs. Therefore, the suggested two step approach may be referred to as the Virtual Array Ray (VAR) technique. The computational formulation of the problem is first introduced and discussed so that field patterns can be synthesized in the subsequent sections. Fields inside the 118

rib cage are computed using the second term of (2.8) where the fields over the rib aperture are employed to compute the radiated velocity potential according to the following equation 0( = ji [(")gf (rl' ) - (')(n'. V'gf (rr'))]dS', (5.1) where S' is the surface over the aperture of the rib cage, b(r) is the radiated velocity potential in m2/sec at a distance r' from a certain source. This source has a normal velocity distribution u(r') = nii' V'4(r') and it is located at r' and gf(rr/') is the appropriate Green's function. Because of the impenetrable ribs, image theory [22, 71] can be invoked to eliminate the second part of the integral. Thus, after discretizing the aperture into P segments, the total field at any location r'd below the ribs can be expressed as P p Otd(rd) = p gfdrd p), (5.2) p=l where up is the complex excitation of the pth virtual element located at r'. It should be noted that rd is the field point location and gfd(rdrFp) is the half space Green's function of the interior region of the rib cage. A similar expression is also valid for the fields generated by the source array/aperture which has K sources/segments. Specifically, K tu(U = - Uk gfu(rurk), (5.3) k=l where Uk is the complex excitation of the kth array element located at rk. In this equation, ru is a general field point in the upper region. The function gfU(ru\rk) is the appropriate Green's function for the region above the solid ribs (upper region). 119

To proceed with the synthesis problem, it is important to express (5.2) and (5.3) in matrix form, where (5.2) can be written as Ad = Hd up, (5.4) where the subscript d indicates that the calculations are performed in the interior region (below the ribs). In this case, up is a P x 1 complex excitation vector of the virtual array and Hd(m, p) represents the field at the location rm due to a unity source at r-. That is, Hd(m,p) = gfd(FmIrp). (5.5) For the upper region, the field matrix representation as deduced from (5.3) will take the form u = Hu u,, (5.6) where the subscript u signifies the upper region (above the ribs), and uu is the discretized source vector of the feed aperture with length K. Also, Hu is the matrix operator whose row dimension is equal to the number of field points at which the velocity potential is computed and its column dimension is K. The entries of this operator are the fields of the kth discrete source evaluated at the jth field point. Thus, Hu(i,k) = gfui(ri rk), (5 7) where r' is the Zth field point location. 120

1, Feed Array k Element K I:-,.:.. _....a..:.. - x A-A y Rib Intercostal Rib bpacing X/2.4 IN 9 Virtual Array Elements A - A Section Figure 5.1: Two dimensional model illustrating the virtual array concept. 121

5.2 Field Pattern Synthesis Using VAR According to Figure 5.1, the feed array (primary sources) is used to illuminate the tissue spacings between the solid ribs and hence induce secondary fields. These secondary fields can be considered as a virtual array responsible for field excitation in the interior region of the rib cage. As mentioned earlier, by introducing the virtual array concept, feed array design can be carried out in two steps. First, the virtual array is designed to generate the desired pattern in the interior of the rib cage. Subsequently, the actual array is computed to generate the desired fields across the intercostal distances between the ribs concurrently with minimizing direct power incidence over the rib surfaces. In both steps of the design, the PI approach (constrained or unconstrained) is employed to compute the virtual array feed. 5.2.1 STEP 1: Virtual Array Design In the first synthesis step, it is required to focus the ultrasound energy at a certain set of locations (spots) within the rib cage while minimizing power deposition elsewhere. Since the foci locations and the corresponding velocity potential vector ad at M field points are known a priori, the synthesis procedure becomes that of specifying the complex excitation vector of the virtual array up such that the corresponding underdetermined system in (5.4) is satisfied. The PI technique (minimum norm least square solution) is employed to determine the virtual sources for a given vector 9Pd [26, 27]. The application of the PI method to (5.4) yields the complex excitation 122

vector of the virtual array = Hp"' d, (5.8) where Hd' represents the PI of Hd and is given by (3.6). It is understood that Hd must be well-conditioned for an accurate solution. This requires that the field points be at some distance from each other [26]. This is the case for the system given in (5.5) which is utilized for computing the virtual array. However, because of the Nyquist sampling requirements, the system in (5.7) will be highly ill-conditioned [5, 7, 26]. Thus, preconditioning is required in order to synthesize a specific field pattern using the PI approach and this point will be addressed in the next section. 5.2.2 STEP 2: Array Synthesis and minimization of intensity levels over the solid Ribs The second step of the synthesis technique involves computing the feed array that minimizes direct power incidence over the rib surfaces while still pursuing a design to recover the specified fields over the intercostal spacings as dictated by STEP 1. According to (5.1) and as shown in Chapter 4 the total localized potential at a general point on a specific rib is given by r = Hr Uu, (5.9) where the general element of Hr is given by Hr(s, k) = gfu(rksrk), (5.10) 123

where rk is the kth segment vector location, rks is the vector connecting the kth array element or aperture segment to the sth field point over the rib. This matrix entry represents the potential of the kth array element at the sth point on the rib plane. As noted earlier, the HRPO synthesis discussed in Chapter 4, did not explicitly minimize the rib fields in the synthesis process. On the other hand, the two-step procedure presented here explicitly accounts for the field levels over the rib surfaces. Therefore, it allows for a constrained optimization analysis to achieve the desired focal intensity levels while minimizing the filed pressure on the solid ribs. This optimization problem is described below. 5.2.3 Constrained Optimization As stated earlier, Hr is the propagation matrix from the surface of the feed array to the frontal plane of the rib surfaces and Hu is a similar operator from the feed array surface to the intercostal spacings. The goal is to solve the following optimization problem: min |IIHuu112 subject to Ip = Huuu,Vu,. (5.11) where 1'p is the velocity potential vector at the locations of the virtual array elements. This optimization problem can be solved using Lagrange multipliers [38]. However, a more elegant solution is obtained via the projection theorem in the form of a weighted minimum norm solution given by u, = WPiH (HuWPiHut)PiPp, (5.12) 124

where W is a positive definite weighting matrix, given by W = HtHr. (5.13) To improve the condition of the matrices Hr and Hu, singular value decomposition (SVD) analysis is employed [56]. The condition number in the preconditioning process specifies the dynamic range of the velocity distribution, but as a rule of thumb, higher condition numbers will lead to higher velocity dynamic ranges and thus lower array efficiencies. 5.3 The Two Step VAR Algorithm The VAR pattern synthesis algorithm can be summarized as follows: 1. Specify the x and y coordinates of a set of M control points. For example, location(s) within the tumor(s) volume or other locations where the heat deposition is kept within certain levels. Also, the corresponding complex velocity potential values are defined at these locations. 2. Calculate the matrix operator Hd and use (5.8) to compute the virtual array over the intercostal spacings. 3. Calculatehe matrix operator Hr which represents the propagation operator from the surface of the feed array to the sampled locations over the solid ribs. As noted earlier, this operator is poorly conditioned due to the closeness of the sample points over the ribs. 125

4. Obtain the weighting matrix W that preconditions the system from (5.13) and perform singular value decomposition [38] to improve the condition of the operator HU. 5. Use equation (5.12) to compute the complex excitation vector of the feed array by finding; uu. After computing the complex excitation vector feed array, field maps interior to the rib cage and over the solid ribs are constructed. It should be noted that in order to compare between the HRPO approach developed in Chapter 4 and the VAR approach developed in this chapter, similar field patterns are implemented and tested with both techniques. 5.4 Numerical Examples Having demonstrated the formulation for field computations and pattern synthesis, a computational example is presented, analyzed and discussed to test the proposed two step algorithm. The model suggested in Figure 5.1 is considered to synthesize single and multiple focus pattern schemes using the VAR approach. The VAR validation is carried out by comparing the results of the VAR approach to those of the HRPO implemented in Chapter 4. To formulate the matrix equations, Green's functions utilized throughout this chapter are replaced with their 2-D expressions [22, 39]. Namely, gf(lrEk) = - H-J ( —jyrkp) (5.14) 2 126

where rkp is the distance between the primary and secondary source points, viz. rkp = y2 + ( -- Xk)2 (5.15) with (Xk, Yk) denoting the kth element coordinates of the array and xp is the location of the pth virtual line source. In (5.14), H(2)(.) represents the zeroth order Hankel function of the second kind. Similarly, for the region below the ribs, the half space Green's function is given by fd(rml rp) = H2) (-jrmp) (5.16) where rmp = /y2 + (p - Xm)2, (5.17) in which (Xm, Yn) are the field point coordinates. Substituting from (5.14) to (5.17) into (5.1) to (5.8), the system matrices for the upper and lower regions are constructed. The synthesis equations (5.8) and (5.12) are then employed to extract the complex excitation vector of the feed array. To test the design, the excitation vector is utilized in reconstructing the field maps under the rib cage and over the rib surfaces. The fields over the solid ribs can be obtained from (5.9) and (5.10). 5.4.1 Single Focus Example Using the aforementioned steps, a single focus pattern synthesis is first carried out. The x and y coordinates for the focus are chosen arbitrarily to be (120 mm, 150 mm) or (40 A, 50 A). The synthesis process is carried out using both HRPO and VAR techniques. In the subsequent maps, intensity values are normalized with respect to 127

the maximum value which usually exists at the focus locations. The synthesis results are shown in Figures 5.2 to 5.9, from which, the following can be concluded: 1. The oscillatory nature of the amplitude distribution of the complex excitation vector of the feed array is observed when the HRPO technique was used. This can be clearly observed from Figure 5.2 where the array elements facing the rib surfaces are minimally fed because of their weak contribution at the focal point. On the other hand, as shown in Figure 5.3, the complex excitation vector when the VAR technique was employed does not have this property, since the array was fed in such a way that minimizes the power deposition over the rib surfaces. 2. At the desired focus location, the specified levels are accurately fulfilled for both approaches as shown in the two dimensional intensity patterns displayed in Figures 5.4 and 5.5. Also, the same fact can be deduced by considering the focal plane cuts in Figures 5.6 and 5.7. 3. Minimal intensity levels are observed everywhere interior to the rib cage, as clearly displayed in the contour plots (25 % contours) in Figures 5.8 and 5.9. This is also apparent from the rib-tissue field patterns presented in Figures 5.10 and 5.11. 4. In the target region below the solid ribs, well behaved interference patterns are observed. This is apparent from Figures 5.4 and 5.5 as well as Figures 5.8 and 5.9. 5. Most of the ultrasound energy propagates through the intercostal spacings 128

avoiding the solid ribs to converge at the desired focal point. This can be observed in Figures 5.10 and 5.11. 6. The VAR technique tends to minimize the power deposition over the rib surfaces compared with the HRPO technique. This is clear from Figures 5.10 and 5.11, where the rib intensity levels in the VAR case are significantly lower than those predicted with the HRPO technique. For this single focus design, there are some advantages associated with the VAR over the HRPO technique which are summarized as follows: 1. The VAR approach appears to lead to a better usage of the intercostal spacings between the ribs. As shown in Figures 5.10 and 5.11, the intensity distribution over the intercostal spacings for the VAR case tends to be more uniform than that of the HRPO. As predicted by these two figures, the maximum intensity level over the intercostal gaps when the VAR approach was employed was approximately 0.04 while in the HRPO, this value was around 0.06 which is about 50 % percent more. Furthermore, the ratio of the power over the ribs to that in the intercostal spacings is about 2.18 % for the VAR technique compared to 2.88 % for the HRPO case. As will be shown, this is even more true for the multiple-focus case. 2. The VAR technique eliminates field singularities near caustics (present in the near bone regions). In the HRPO technique, these singularities were due to the approximate PO evaluation [15]. 129

0 50 100 150 200 Array Elements 2-'" " [. ', JIl lt hi2i - 4.. 0 50 100 150 200 250 Array Elements Figure 5.2: Complex excitation of the feed array using the VAR method (single focus). Array efficiency TA= 8.28 %. 0.7 0.6 0.5 -a f10.4 g 0.3 0a 3.2 0.1 0...I 250 Array Elements 4.l.|. 24 ~~..........." I 0 50 100 150 200 25 50 Array Elements Figure 5.3: Complex excitation of the feed array using the HRPO method ( single focus). Array Efficiency rA= 21.93 %. 130

1. 0.8.. t c 0.6. C 0.4, o z 0.2 - I 250. 2!50 "' 300 250 x-axis (mm) y-axis (mm) Figure 5.4: Two dimensional intensity pattern using the VAR method (single focus). 1 -0.8, a 0.6,. 0.4, z 0.2, [ 0>l[ 250.- *300 250 0 0 x-axis (mm) y-axis (mm) Figure 5.5: Two dimensional intensity pattern using the HRPO method (single focus). 131

1 I 0.9 -................ 0.8 -............................................................... 0.7.................. -...... 0.6............................................................................................ x-as (mm)0.6 e 0.5: Focal plane intensity cut using the VAR method zfocus). 0.2 0.1 0 so50 100 150 200 250 x-axis (mm) Figure 5.6: Focal plane intensity cut using the VAR method (single focus). 1r I. 1 l 0................................................................................8 -............................ 0.7 - 8.......................................................................................................................... 0.2 -............................................................ 0.5 -............................................................................. 0 50 100 150 200 x-axis (mm) 0.3. x-axis (mm) Figure 5.7: Focal plane intensity cut using the HRPO focus). 250 method (single 132

5odn I 2 2 0 -................................... 2 0 0 -................................. 1 8 0 -................................. 1 6 0...................................................... - 1 4 0.............................................. E. 12 0 -................ O1 0 0 O-................................................... 8 0 -... *.:............................... 80 6 0 - * * * *....:........................................ 60: 4 0.............................................. 2 0 - - -.................................................... 20 50 100 150 y-axis (mm) Figure 5.8: Intensity pattern contour plot VAR method (single focus). 2440;220................................................ 2 0 0..............................................!18 0 ~ ~...........................................:..... '180 1160......................... 1:4 0 -................................................................. E 8 0 -....................................... 60 0 -...... * *:............................................................... 80.6 0 -.......................... 40............................................. 20..................................................... & 200 250 (25 % contours) using the................................................................................................................................................................................ 50 100 150 y-axis (mm) Figure 5.9: Intensity pattern contour plot HRPO method (single focus). 200 250 (25 % contours) using the 133

5.4.2 Multiple Focus Example Application of the of the VAR synthesis technique is also extended for synthesizing multiple focus pattern. In this case, focusing is simultaneously performed at two different locations. Both foci are assumed to have unity intensity and are located at (60 mm, 165 mm) and (180 mm, 135 mm) or (20 A, 55 A) and (60 A, 45 A), respectively. The simulation results are illustrated in Figures 5.12 to 5.21. From these figures, the following can be deduced: 1. Figures 5.12 and 5.13 show the complex excitation vector of the feed array for both the VAR and HRPO approaches, respectively. As seen from the amplitude distribution, the array excitation efficiency using the HRPO technique is much higher compared to that using the VAR technique because preconditioning included in the second step of the synthesis process increases the dynamic range of the excitation vector. An algorithm for increasing the excitation efficiency has been proposed in [26]. However, the array efficiency which is a measure of the field uniformity over the array elements does not have to be high (close to 100%). This is due to the fact that some array elements are excited with low power in order to achieve signal cancellation over the strongly scattering obstacles. It is worth noting that in the simulations, an optimization technique employed to increase the array excitation efficiency was tested for several multiple focusing schemes [26]. It was deduced that field levels over the ribs increase dramatically when the array becomes nearly phase driven. Thus, for minimizing the rib fields, the efficiency becomes of less importance. 134

0I.0. l I I 0.04 0.035 -... 0.03.c _D 0.025 0.02 z -.....................;.:....:...:.............1 A::: A A, 1:.......i 0.0151- ~ 0.01 -.. o.0o5 - I... -. —....-. -- I V0 50 100 150 Rib-Tissue axis (mm) 200 250 Figure 5.10: Intensity pattern over the rib-tissue plane using the VAR method (single focus). Relative rib power cr= 2.18 %. 100 150 Rib-Tissue axis (mm) 250 Figure 5.11: Intensity pattern over the rib-tissue plane using the HRPO method (single focus). Relative rib power cr= 2.88 %. 135

2. Figure 5.1]4 shows the intensity pattern inside the rib cage after computing the linear array via the two step VAR technique. This field map demonstrates that the prespecified intensity levels are accurately recovered at the foci locations. Also, low interference pattern levels are observed inside the rib cage. Figure 5.15 shows the same example using the array synthesized using the HRPO technique. As seen, interference patterns from the HRPO are more deteriorated than those achieved when the VAR approach was employed for pattern synthesis. This is particularly true especially in the region directly below the ribs where the physical optics technique has known limitations. 3. It is clear from Figures 5.16 and 5.17 that at the focal plane, both technique delivers the desired power at the specified locations. The VAR technique as seen in Figure 5.16 produces a slight distortion in the focal amplitude due to the preconditioning of the system matrix in the second step of the synthesis procedure. After testing several examples, the observed distortion in the intensity level at the focal point was less than 1 % and thus, this distortion can be neglected. Over the focal plane, similar interference patterns are observed for both techniques. 4. The contour plots given in Figures 5.18 and 5.19 demonstrate well behaved interference patterns inside the rib cage when the VAR technique was employed for pattern synthesis. It should be noted that both contour plots were generated based on 25 % contours in the intensity patterns. The efficient utilization of the acoustical window presented by the intercostal spacings minimizes the 136

interference patterns inside the rib cage. 5. Figures 5.20 and 5.21 show the intensity level over the rib-tissue plane. It is clear that the rib fields generated when the VAR technique was employed are significantly lower than those resulting from the HRPO technique. The rib minimization constraint included in the second synthesis step minimizes direct power deposition over the rib surfaces In summary, one of the most important goals in array design is to avoid illumination of the solid ribs and to propagate most of the ultrasound energy through the intercostal spacings. Both techniques achieve this goal but in the HRPO case, the intensity levels over the ribs are not explicitly controlled. On the other hand, for the two step VAR algorithm, minimization of the rib intensity levels is part of the synthesis process. Also, by controlling rib sampling and condition of the matrices, rib fields can be more precisely controlled. For the examples considered here, the VAR technique results in a rib power ratio 6r ranges from one to three percent whereas the corresponding ratio for the HRPO ranges from five to ten percent. 137

4 u2 -4 0 50 100 150 200 250 Array Elements Figure 5.12: Complex excitation of the feed array using the VAR method (two foci) (a) Relative amplitudes (b) Relative phases. The Array excitation efficiency 7-A= 10.43 %. 0 u.o X0.6 > 0.4 c 0.2 Array Elements 0 a. CD co 4 4 llll lJ i.-2........ -4 0 50 100 150 200 250 Array Elements Figure 5.13: Complex excitation of the feed array using the HRPO method (two foci) (a) Relative amplitudes (b) Relative phases. The Array excitation efficiency 7A= 19.96 %. 138

0.8. a. 0.6. N 0.4. z 0.2,. 0. 250....,.,...,.....,..:,.... -;- '' I s b x-axis (mm) y-axis (mm) Figure 5.14: Two dimensional intensity pattern method (two foci).,~~ 300 250 using the VAR 300 250 11.4..20 o 0.8. 0.6, 0.4,. 0.2.. 0.. 250 I x-axis (mm) y-axis (mm) Figure 5.15: Two dimensional intensity pattern using the HRPO method (two foci). 139

I I I - I 0^0.8 C 0.6 0 0.4 0.2 Z 0.2........'........................... 0 5 100 150 2 I 2 0 50 100 150 200 250 x-axis (mm) Figure 5.16: Focal plane intensity cuts using the VAR method (two foci).,.4, 1.2 c" 1 N 0.6 S E I 0.4 Z z 0.2.....................................................................................................................................................:........... 0, - I, 0 50 100 x-axis (mm) 150 200 2 0 0 50 100 150 200 250 x-axis (mm) Figure 5.17: Focal plane intensity cuts using the HRPO method (two foci). 140

Z4U,,I I r I I 220 200 180 160 -- 140 E E. 120 100 80 60 40 20 Figure 5.18:....................................................... *.......... 0.:i -............................................... c 0o -..................................................................................................................... I I I... 50 100 150 200 y-axis (mm) Intensity pattern contour plot (25 VAR method (two foci). 250 % contours) using the 240. I 2 2 0............................................................................... - 200 7~~ - o.. 20...................................... - 180 - 0.................................. 1 4 0o:.................................................... E I._120 I=?,-: oC> o 6 0........................... 4 0 -.............................................................................................. 2 0 '............................................ 0 100 150 200 250 y-axis (mm) Figure 5.19: Intensity pattern contour plot (25 % contours) using the HRPO method (two foci). 141 141

0 50 100 150 200 250 Rib-Tissue axis (mm) Figure 5.20: Intensity pattern over the rib-tissue plane using the VAR method. Relative rib power 6r= 2.47 %. 0.25 1 ----------------- i ----- i ----- 0.25 oo.1........,,.... M.... 1.......... 0.0 0 50 100 150 200 250 Rib-Tissue axis (mm) Figure 5.21: Intensity pattern over the rib-tissue plane using the HRPO method. Relative rib power cr= 6.79 %. 142

5.5 CPU Computational Cost Comparison between the HRPO and VAR Techniques One of the main advantages of using high frequency methods is their affordable memory and CPU requirements. Compared to numerical methods, high frequency techniques require significantly lower CPU and memory demands. However, for calculating the field within a domain that spans tens or hundreds of wavelengths, CPU and memory issues become of increasing importance. Therefore, in this section, the CPU and memory costs are compared for both HRPO and VAR techniques. Before comparing between the HRPO and VAR costs, the following two points are explicitly addressed: * For both techniques, memory may not be considered as important since the storage needed for both methods is usually proportional to multiplication of the number of array elements (either actual or virtual) by the field control points. Practically, the number of foci that can be generated is limited to one or two as maximum (with deep focusing) and that of the array elements is around tens or hundreds. As a result, there is no memory problem associated with either approaches. * Regarding the CPU cost, the main part of the CPU is consumed in the scanning process. This is due to the fact that after the synthesis is completed and the feed array is computed, fields within the computational domain are evaluated at each point. This process has intensive CPU requirements because of the 143

high number of domain grid points. Therefore, in this section, only the CPU requirements for field scanning in both techniques are discussed. To characterize both approaches, it is essential to introduce the basic parameters that contribute to the CPU cost. These parameters are: * Number of feed array elements is K. * Number of control points (foci) is M. * Number of ribs is N. * Number of intercostal spacings between the ribs is N - 1. * Interelement separation between the virtual array elements is A * Number of virtual array elements is 21(N-1) where 1 is the intercostal spacing A between the ribs. * Number of segments per wavelength for the PO integration utilized in the HRPO approach is Ipo. * Number of field points to be scanned in the x-direction after the feed array is computed is N,. * Number of field points to be scanned in the y-direction after the feed array is computed is Ny. 144

5.5.1 HRPO CPU Cost When the HRPO technique is employed to predict the fields within the entire computational domain, several processes are involved requiring different number of operations. These processes will be considered separately as follows: * The PO Integration Scheme: Assuming that the number of PO segments within each intercostal spacing per wavelength is Ipo and by applying fourth order Gaussian integration scheme to evaluate the resulting PO integration, then (4Ipo) operations will result per intercostal spacing per wavelength. This yields (4Ipo l(N-1) operations for calA culating the total number of operations required to perform the PO integration. * Effect of Feed Array: The calculation performed in the previous item assumes single incidence over the rib-tissue plane. However, in the actual case, K array elements are utilized to excite the rib-tissue fields. That is, the total number of operations is 4Ipo l() K. This number represents the total number of operations required to predict the field at any given location in the interior of the rib cage. * Field Map Reconstruction: For a two dimensional region of interest composed of NzNy points, the total number of operations will be Nhrpo = 4Ipo ( 1)KNNy, (5.1 8) A 145

where Nhrpo is the total number of operations required by the HRPO scannring scheme. As indicated in (5.18), extremely high number of operations is expected when the scanning process is performed using the HRPO computational method. With the increase of the operating frequency, the number of array elements increases (assuming fixed interelement separation of half wavelength) and also the number of field points scanned in the domain (scanning step is one wavelength). This results in a dramatic increase in the executed number of operations and consequently the CPU time. 5.5.2 VAR CPU Cost For the VAR case, the following processes are included in the scanning process: * Virtual Array Construction According to the half wavelength interelement spacing of the virtual sources, the number of virtual sources per intercostal spacing is given by j2. Considering all the intercostal spacings, a number of 21(N-1) virtual elements exist within the rib-tissue plan. * Effect of Feed Array: To construct a certain virtual source when K elements of the feed array are excited, a number of K operations is required. Therefore, a total number of K 2(-1) operations is required to construct the whole virtual array. * Field Map Reconstruction: 146

After the virtual array is computed over the intercostal spacings, it is considered as the feed array for the lower region. The number of operations required to scan a region composed of NxNy field points can be obtained by Nvar = K A ( ) + ( )NN (519) or Nvar = 2(A )(K + NxNy). (5.20) It should be noted that the main difference between (5.18) and (5.19) is that the in the VAR case, the synthesis process is divided into two complete subprocesses. Thus, when scanning the domain field after extracting the excitation vector of the feed array, field points under the rib surfaces have no contribution on the process of computing the virtual array. This process is represented by the first term of (5.19). On the other hand, after the virtual array excitation is deduced from the first scanning step, the problem will be equivalent to a linear array radiating in half space and this reduces the complexity of operations. The scanning process under the rib surfaces is represented by the second term of (5.19). As a conclusion, the HRPO scanning process may be regarded as equivalent to the multiplication of both terms while the VAR scanning is the summation between these terms. 147

5.5.3 CPU Cost Comparison To compare between the CPU cost for both techniques, a CPU ratio factor Q is introduced and defined as Nvar Q = Nr (5.21) Nhrpo Substituting from (5.18) and (5.19) into (5.21), the CPU ratio factor will be K + NNy = K~NnN (5.22) 2 IpoK NN (5.22) It should be noted that, for a fixed problem, all the quantities are dependent on the operation wavelength. Thus, when the treatment frequency changes, K, Nx and Ny varies and consequently the operations ratio Q. Therefore, the CPU cost for both approaches will be a frequency dependent quantity and as a result, it is recommenced to define a frequency factor and express all parameters as functions of this factor. This factor may be defined as F f(kHZ) (5.23) 500 where f is the operating frequency and 500 represents the 500 kHz excitation case. Based on the frequency factor definition and assuming the geometrical dimensions mentioned in the aforementioned examples, the frequency dependent values for all array and rib parameters are given by * K=240 F,. * NA= 80 Fn. 148

. N,= 90 Fn. The numbers 80 and 90 correspond to the domain size in wavelengths at 500 kHz operating frequency assuming one wavelength step size in the scanning process. In the HRPO simulations, convergence is achieved when the number of PO segments per wavelength is around 3 to 5. An intermediate value of 4 is chosen, i.e., Ipo=4. The result of the CPU comparison is plotted in Figure 5.22 showing that the VAR computational approach provides dramatically better performance. For all operating frequencies, the VAR is much faster in scanning the fields in the interior region of the rib cage with a speed gain of at least 2000. Therefore, in all subsequent simulations, the VAR computational technique will be employed for field maps reconstruction. It should be noted that differences between the HRPO and VAR field values for a given excitation vector are negligible and both approaches yield almost the same field maps everywhere. 5.6 Discussion and Conclusions In this chapter, a two-step hybrid approach was presented for performing deep focusing of ultrasonic waves in the presence of strongly scattering objects, e.g., the rib cage. This method was specifically chosen because it provided more control on both upper and lower regions of the rib plane (see Figure 5.1). Most importantly, this two step VAR approach achieved better interference patterns coupled with rib power minimization over the rib surfaces. In addition, the 149

VAR approach required much less CPU time to generate the intensity map in the region interior to the rib cage. This is because Nyquist sampling rate was used over the intercostal spacings and thus transforming the complicated continuous source integration performed in the HRPO approach to a simple reduced summation. The proposed technique demonstrates the feasibility of focusing high intensity ultrasonic beams deeply in the liver vicinity. The efficient utilization of the acoustical window between the solid ribs made it possible to treat organs placed in the shadowed region of strong inhomogeneities such as the ribs cage. As displayed in the simulations, focusing HIFU beams to distances of order 10 cm to 15 cm with a reasonable spot sizes (less than 1 cm in diameter) is possible. The size of the focal spot is primarily determined by the array f number and the operating wavelength A. In addition, a shape factor accounting for the array shape is also an important factor [26, 60]. The developed VAR technique eliminates the utilization of the PO integration method which has poor prediction accuracy in certain parts of the computational domain in addition to its enormous CPU requirements. Also, the VAR approach explicitly accounts for the rib power minimization. 150

750 E Frequency (kHz) Figure 5.22: CPU ratio between the VAR and HRPO techniques. 151

CHAPTER 6 Improved Geometrical and Computational Modeling of the Problem In Chapters 4 and 5, high frequency approximate techniques were employed to predict the ultrasonic fields within the entire computational domain. A simplified geometrical model was proposed in both chapters to address two main issues. The first was to examine the feasibility of focusing in the presence of strongly scattering obstacles such as the rib cage. The second was to study the behavior of field synthesis algorithms in terms of the field patterns quality inside the rib cage and direct power deposition over the rib surfaces. Pattern synthesis techniques capable of producing highly localized fields in the tumor (target) vicinity were developed and validated in the earlier chapters. These chapters demonstrated the feasibility of generating deeply focused patterns inside the rib cage. However, a more realistic computational model with actual geometrical representation and tissue characteristics is needed. Concurrently, more accurate field computational techniques are necessary for ultrasonic field prediction within the computational domain. 152

In this chapter, a more realistic model is proposed based on the actual geometry, layer parameters and operating conditions. This model employs the eigen function representation [22, 39, 81] for modeling the ultrasound scattering from the ribs. Finally, field pattern synthesis is carried out via the two step VAR technique developed in Chapter 5 to compute the feed array. Several examples are studied and implemented followed by a detailed discussion at the end of the chapter. 6.1 The Simplified Geometrical Modeling and the Approximate High Frequency Computational Techniques According to the results presented in the last two chapters, it can be concluded that focusing in the presence of the rib cage is feasible. However, with the geometrical model suggested in Chapter 2 and high frequency techniques employed in Chapters 4 and 5, the proposed model has the following drawbacks: 1. It lacks realistic modeling of the rib shapes, dimensions and characteristics. As displayed in Figure 2.2, the ribs were modeled as perfectly solid (non-penetrable) thin plates. In reality, they tend to be cylindrical in shape with an elliptical or circular cross section. Also, they are not perfectly solid but lossy elastic. 2. The attenuation in all tissue layers including the ribs was ignored and as presented in Chapters 4 and 5, numerical simulations were tested for lossless cases 153

only. Regarding the accuracy of the employed computational techniques, rib scattering and partial transmission through the ribs were ignored. The logic behind ignoring both factors was based on the following: 1. At any stage of the synthesis process, if focusing through such strong inhomogeneities was not feasible, there is no need to improve the computational accuracy. Fortunately, well behaved patterns in the presence of the rib obstacles were achievable [5, 7, 15]. 2. In the region of interest (tumor locations), terms resulting from the rib obstacles introduce weak contribution and can thus be neglected. 3. Exact methods such as the eigen field expansion (Neumann Expansion) approach [17, 39, 64, 71] have intensive analytical as well as computational requirements especially for such a large geometry. Thus, a simple model was initially suggested to examine the feasibility of focusing in the presence of the rib obstacles. In conclusion, Chapters 4 and 5 dealt with developing efficient and reliable pattern synthesis schemes for a simplified model of the problem. Also, these chapters presented introductory steps towards the solution of the general problem. They simply resulted in an answer to the following questions: * Can we focus high intensity ultrasonic beams in the presence of strongly scattering obstacles such as the rib cage? The answer was YES and as described 154

in Chapter 4, a computational model was developed, tested and then validated with a real experiment [15]. * Can optimal, well behaved interference patterns be synthesized inside the rib cage with minimal power deposition over the rib surfaces? Again, the answer was YES and as shown in Chapter 5, the two step VAR synthesis approach was developed and tested. If the answer to either of the previous two questions was negative, there would have been no point in proposing a more complicated model or accurate field prediction techniques. Thus, the goal in this chapter is to achieve the following: * Propose a more accurate model that matches the realistic geometrical features and tissue properties involved in the focusing problem. * Develop better field computational techniques based on exact approaches and employ them to predict, with high precision, the total field within the computational domain including the interior of the rib cage, rib surfaces and near bone area. 6.2 Modeling and Computational Formulation To account for the geometrical features of the rib cage, the model displayed in Figure 6.1 is proposed. As shown, the ribs are modeled as an array of N circular cylinders with infinite depth (two dimensional). Also, Figure 6.1 shows a linear feed array composed of K line sources. This array excites the ultrasonic fields within the 155

domain and as mentioned in Chapter 2, successful pattern synthesis should yield well behaved patterns inside the rib cage concurrently with minimal power deposition over the rib surfaces. The two dimensional geometrical presentation of the model proposed in Figure 6.1 is displayed in Figure 6.2. As shown in this figure, a linear array composed of K identical line sources is placed on the top of a water layer that separates it from the human body. The ribs are simulated as circular cylinders and assumed to be lossy elastic. To analyze this problem efficiently, it is more convenient to start by a simple case where a line source is radiating in the presence of a lossy elastic cylinder of infinite depth in the axial direction. This simple case is shown in Figure 6.3, and its associated parameters and description are presented in Figures 6.4 and 6.5. The computational analysis for this simplified case is considered as the basic step towards the formulation of field prediction mechanisms applied to the general problem modeled by Figures 6.1 and 6.2. Then, the two step synthesis approach presented in Chapter 5 and in [5, 7] is employed to synthesize the field patterns inside the rib cage. The computational formulation and synthesis technique are validated by testing several examples with realistic data. Finally, the attenuation effects are characterized as a function of the operating frequency. 156

Feed Array The Ribs The;iver fof the prPea6.1 Representation of the eoetrical Figure gation model. 157

1 2 k K-1 K 0 - ------ rk 1 n N Field Point y Figure 6.2: Two dimensional geometrical representation of the problem at hand. 158

Line Source - 1' Incident Wave Scattered Wave. Wave Lossy Elastic Cylinder Figure 6.3: A line source radiates in the presence of a lossy elastic cylinder. 159

Line Source Y \.....R..y r T Lossy Elastic Cylinder Field Point Figure 6.4: Schematic representation of the problem described in Figure 6.3 160

Line Source,, A Cylindrical Waves 1) x y Elastic Cylinder inc i I J i q i I / Field Point * sea tot inc + sea Figure 6.5: Two dimensional model illustrating the field computational concept applied to the problem described in Figure 6.3. 161

6.3 Simple Case: One Line Source Radiates in the Presence of a Lossy Elastic Cylinder Figure 6.3 shows a line source illuminating a lossy elastic cylinder, with infinite depth along the axial direction. This is a standard scattering problem and the typical approach to solve this problem was developed and briefly discussed in [22, 39, 81]. However, the solution is only feasible for cylindrically shaped scattering objects. Since the ribs are assumed to be circular cylinders, cylindrical basis functions can be employed for the field expansion. Let Red i represent a regular cylindrical basis function and Ou ci represent an outgoing cylindrical basis function [39]. These functions can be mathematically expressed as [22, 81], Re bi(7y, r') = Ji(7, Ir) ei, (6.1) and Ou =(, ' ) = H)(ye, Ir) e3 (6.2) where Ji is the ith order Bessel function, Hi2) is the ith order cylindrical Hankel functions of the second kind, i is an integer in the range [-oo, oo] and - is the complex propagation constant defined in (3.8). The point r' in cylindrical coordinate system has a radial distance of Ir1 and an angle of 0 with the origin located at the scatterer (cylinder) center. The mathematical representation of (6.1) and (6.2) takes 162

the following form Z(7, |F|) = Z(-j7I|FI), (6~3) where Zi represents either a Bessel or a Hankel function. As indicated in [17, 22, 39, 71], the incident field can be written as = ainc() ai Re qi(^7, F), (6.4) where ai is the incident field expansion coefficient for the ith term of the expansion. Similarly, the scattered field can be written as a function of the outgoing cylindrical eigen functions, viz ksca() = bi Ou /i(^, r), (6.5) where bi is the scattered field expansion coefficient for the jth term of the expansion. Also, the total transmitted field in the elastic cylinder of radius ac can be expressed as 4tra( )- ei Re, (7, r), RI _< ac, (6.6) where ei is the total transmitted field expansion coefficient for the ith term of the expansion and yc is the propagation constant of the cylinder. The incident field coefficients can be evaluated by assuming half space Green's functions of the radiating element which is considered as a line source placed at the location r"', that is in (r)- 2 (t r-rl) (6.7) 163

Applying the addition theorem to the right hand side of (6.7) allows us to write the incident field as >inc(r) = Ou (y, r ) Re i(y, ), (6i.8) where r' is the smaller of r' and r', and r'i is the larger of r' and r'. Assuming that r'< F', then, (6.8) becomes inc (r) = Ou i(y, ') Re i(7 ). (69) Comparing (6.4) and (6.9), the ith expansion coefficient can be deduced and expressed as -j ai - Ou y,). (6.10) Substituting from (6.2), the coefficient ai takes the form a = -J H()(, I (6.11) 2 where the angle 0' is the line source angle considering the origin located at the scatterer center. For a solution of Osca, boundary conditions discussed in Chapter 2 are applied. That is, continuity of the pressure and normal velocity over the surface of the scatterers are enforced, giving p(oinc(r) + 4sca(r)) r=ac = Pcktra(r) tr=ac, (6.12) and a(inc(r ) + 4sca(r ))Ir=ac a= tra(FIr=ac, (6.-3) Or d 164

Substituting (6.1) to (6.5) into the boundary conditions and applying the orthogonality of the cylindrical basis functions, results in two equations for bi and ei and solving yields p(aiJi(yac) + biH2)(yac)) ei =, (6.14) PcJi(7cac) and b PcyJi(%cac)J(Tac) - pyJ(cya,)J(yac) (6 b-) ai. (6.15) p7cJi'(acac)H2) (-ac) - pc-Ji(cac)(H2) '(/ac))' More compactly bi = bai, (6.16) where i pc=PcJi( '(Ycac)J'(ac) - pycJ' (cac)Jiac) (6.17) pcJi'(7cac) H(2) (-ac) - pcyJi(,.cac)(H(2) (ac))' The basic eigen formulation in this section will be extensively employed when the feed array is excited to illuminate the rib cage which is the source of inhomogeneity within the domain. 6.4 General Case: K Line Sources and N Lossy Elastic Cylinders In this section, results developed in the previous section are employed to obtain the total field when K line sources illuminate N circular elastic cylinders (see Figures 6.1 and 6.2). The total field at any point within the computational domain (except 165

th k Array Element A / \ rkn \. X n rko YO,i A^ ko xto X0 N I - if kn Xn I — ' x On to / o Rib Ii th n Rib v y Field Point Figure 6.6: Geometrical details of the problem. 166

for the inhomogeneities) can be expressed as tot (Fr) = qinc(r) + sca (r), (6..18) where 'tot(r) is the total field at the location rf. To calculate each field separately, (6.4) and (6.5) are employed. Due to the fact that the incident field value should be the same regardless of how it is calculated, any rib can be chosen arbitrarily for calculating the incident field. After accounting for contributions of all array elements, the incident field can be expressed as K =inc(r) akJ(y, Iro )ejio, (6.19) k=l i where ro and 00 are the location and angle of the field point, respectively (see Figure 6.6). The scattered field is obtained by superimposing scattered fields from all cylinders, giving K N (r) E= bkH(2)(y, I)e, (6.20) k=1 i n=1 where rn and On are the field location and angle, respectively, referenced to the center of the nth rib as shown in Figure 6.6. Combining (6.19) and (6.20), the total field is given by tot () = E E r( 1 ]io + bE H 2)(y, |r I) j * (6.21) k=l i n=l Considering each element of the feed array is fed separately by a continuous wave signal of relative amplitude Uk and phase Ok, the series coefficients in (6.21) can be 167

replaced by their equivalent from (6.11) and (6.15) to get K r kue) = 2 E I I E H2) )eJ1) y, Iro \)eJ?3 k=l i N + E, Iknl)H2)(?n )ejknein, (6.22) n=1 where Uke~3k is the excitation of the kth array element, rk0 is the distance between the kth element of the feed array to a certain rib or cylindrical scatterer and r'k' is that between the kth element of the feed array and the nth rib, and Oko and 6kn are the corresponding angles, respectively as displayed in Figure 6.6. In the two step formulation described in Chapter 5, a virtual array is introduced along the intercostal spacings between the ribs. As shown in Figure 6.7, the virtual sources are separated by the minimum Nyquist spatial distance which corresponds to A/2 spacings between array elements. To apply this approach, the matrix equation that relates both actual and virtual array feeds should be formulated first. This can be achieved by computing the normal velocity distribution over the virtual aperture using the following equation [17, 86], p(rp)= n.(VtO0t) =rp, (6.23) where ri is a general point on the virtual plane. Substituting for qtot from (6.22) into (6.23) and using Bessel and Hankel functions identities [2], the normal particle velocity up(rp) at the general point r- over the virtual plane is given by Up(-p) = - E ukejk {H 2)(', Irkoejk~ J+(Y o l) + Ji(^Y, Ir Pol) 2 r i r n=1 11 PI 168

x y 4 Rib th virtual p Element i ~n-1 / I Rib Acoustical Window Virtual Elements Figure 6.7: Virtual Array Concept. 169

where rpo is the vector connecting the pth element of the virtual array to the center of an arbitrary rib and 0O is the corresponding angle. In matrix form, (6.24) can be expressed as p = Hui u,, (6.25) where Hu1 is the matrix operator that relates the velocity distribution over the virtual array to that over the actual array (above the ribs). From (6.24), the entries of Hu1 are given by H1(pk) = E {H 2)(l, Irkko ) Ji+(, P'o i,) i Hl J (P, k) D JJo 2. "1rpI + Z k o(2 )(,r I|l)eikn[H() (-y I\ I) - H2)(Y, Ipo I)] eiin n=1 (6.26) As indicated in Chapter 5, the next step is to use the virtual array as the excitation for the approximate half space below its surface. Therefore, the total velocity potential at the general mth point below the virtual array is given by P (rm) = Up(Fp)gfd(rmrp), (6.27) p=l where P is the number of virtual sources, rm is the field point and rp is the virtual source location while gfd(rm rp) is the appropriate Green's function. In matrix form, (6.27) can be expressed as ad = Hd up, (6.28) where itd is the velocity potential vector at the location r and up is the velocity distribution of the virtual array between the intercostal spacings. The general entry 170

of the matrix Hd is given by Hd(m, p) = gfd(rm rp). (6.29) It should be noted that (6.25) and (6.28) are similar to those obtained in Chapter 5 when the two step VAR formulation was introduced for pattern synthesis. The only difference is that (6.25) relates the velocity distribution of the actual array to that of the virtual one via the matrix HU1, while (5.9) relates the velocity potential over the virtual plane to the velocity distribution over the actual one. 6.5 Pattern Synthesis By employing the computational formulation discussed in the aforementioned section, the two step synthesis approach presented in Chapter 5 and [5, 7] can be utilized to compute the complex excitation vector of the feed array. The key steps for this approach can be summarized as follows: 1. Obtain the excitation vector of the virtual array that is capable of producing highly localized fields inside the rib cage. This can be performed by applying the direct PI approach to (6.28) to get up = HdP d, (6.30) where ad is the complex velocity potential vector at a set of M control field points. The PI of Hd is given by (3.6). 2. Formulate the constrained optimization problem by evaluating the incident 171

fields over the solid ribs which are given by (5.9) and (5.10). Hence, the optimization problem can be stated as min ||IIHruu2 subject to up = HuuIu V uu (6 31) The solution to (6.31) will take the following form UU = WPHt(HuiWPiH*t)Pup (6.32) where W is a positive definite weighting matrix given by W = HrtHr. (6.33) After computing the feed array, the total fields are scanned inside the rib cage and over the ribs to verify the pattern synthesis algorithm. The scanning process is performed by first evaluating the matrices Hu1, Hr and Hd. Then, the virtual array excitation vector is computed according to (6.25). Subsequently, the fields inside the rib cage are calculated by substituting in (6.28) using the virtual array excitation vector. 6.6 Numerical Examples After presenting the basic synthesis and analysis approaches for the realistic model presented in the previous sections, several examples for pattern synthesis are considered. The same single and multiple focus examples presented in Chapters 4 and 5 are simulated with 500 kHz to validate the computational approach for the realistic model. In addition, a higher excitation frequency of 1 MHz is used to study focusing 172

with higher frequencies. Both single and multiple focusing schemes will be tested for lossless cases as in Chapters 4 and 5 for comparison purposes. Finally, attenuation and frequency dependence effects are examined. 6.6.1 500 kHz Operating Frequency with no Attenuation In Chapters 4 and 5, the proposed synthesis approaches were tested at 500 kHz. The lossless nature of the soft tissues is maintained here to validate the suggested computational and synthesis approaches by comparing field patters with those generated in Chapter 5. The simulation has the following parameters: * Operating frequency is 500 kHz which corresponds to a 3 mm wavelength in the soft layers. * 240 planar, uniformly spaced excitation elements with interelement separation equal to a half wavelength. * Height of the array elements from the rib plane is 17A. * Number of ribs is 6. * Rib diameter is 6 A which corresponds to 1.8 cm. This value was chosen since its represents the actual rib diameter for an average person. * Aperture (tissue between ribs) width is 7A, which corresponds to 2.1 cm separation on the average. * Rib attenuation is 1.64 Np/(cm MHz) which corresponds to 0.366 Np/A. 173

* Sound speed in ribs is almost double that in soft tissues as indicated in Table 2.1. * Rib density is approximately one and half times that of soft tissues as shown in Table 2.1. * Proposed target location is (120 mm, 150 mm) or (40 A, 50 A). This set of data is applied to the VAR synthesis approach described in the previous section. After applying the feed array excitation vector, field maps inside the rib cage and over the rib surfaces are constructed to test the synthesis algorithm. The simulation results for this example are displayed in Figures 6.8 to 6.12. Figure 6.8 displays the magnitude and phase distribution of the feed array excitation vector and as can be seen, both distributions follow the same patterns observed in Figure 5.2 for the same example. Figures 6.9 and 6.11, indicate that ultrasonic beams propagate through the intercostal gaps between the solid ribs to add constructively at the desired focal location. The focal plane cut shown in Figure 6.10 proves that the prespecified level was achieved accurately at the focus location, a feature that exists in the synthesis approaches presented in Chapters 4 and 5. Figure 6.12 demonstrates that the power deposition over the rib surfaces with the realistic modeling is minimized and not affected either by the geometrical features of the ribs or their ultrasonic characteristics. The relative rib power ratio Er for this test case is about 6.8 %o. Regarding multiple focus pattern synthesis, a two foci pattern is synthesized and discussed. Two foci each has unity intensity are simulated at the locations (60 mm, 174

I) 0 0 1. -' E 0 a: fron 0 50 100 150 200 250 Array Elements 4 i,,i0 i ).... c-2 J.... 0 50 100 150 200 250 Array Elements Figure 6.8: Complex excitation of the feed array for the single focus example at 500 kHz. 175

1.4. 1.2 /) c | 0.86. 0..: z 0.4 0.2. 250 '' ~ '... ~ ~ ~. ~ ~~ ~ ~ ~ ~ ~: ~ ~ ~ - ~ ~ ~ ' ~ ' ~........ ~... ~-,.. ~ *. 250 200 300 150 200 100 150 50 100 50 x-axis (mm) y-axis (mm) y-axis (mm) Figure 6.9: Two dimensional intensity pattern for the single focus example at 500 kHz. 176

.t au c - 0 Z 0 50 100 150 200 250 x-axis (mm) Figure 6.10: Focal plane 500 kHz. intensity cut for the single focus example at 177

240n 220 200 180 160 E E 140. 120 x 100 80 60 40 20.................................................................................................................................................:~:..............................., O~:...................................................................................................... ~......................................... - -...........................................-....................................................... n...................................... -............................................L.........**********................................................................................ - I a. 50 100 150 200 250 y-axis (mm) Figure 6.11: Intensity pattern contour plot (25 % contours) for the single focus example at 500 kHz. 178

0.025 0.02 a) c" C 'O 0.015 N 0 z 0.01...................................................................::............................................................................................::............................................................ h: 0 50 100 Rib-Tissue axis (mm) 150 200 250 Figure 6.12: Intensity pattern over the rib-tissue focus example at 500 kHz. plane for the single 179

165 mm) and (180 mm, 135 mm) or (20 A, 55 A) and (60 A, 45 A), respectively. The complex excitation vector of the feed array is plotted in Figure 6.13, indicating higher amplitude variations over the array elements. This is expected when multiple focus synthesis is considered. The scanned field intensity patterns and focal plane cuts displayed in Figures 6.14 and 6.15 show that deep focusing with acceptable interference patterns is still achievable. This was also deduced from the contour plot given in Figure 6.16. The differences between the field patterns obtained in this chapter and those obtained in Figures 5.14, 5.16 and 5.18 are negligible and in both cases, well behaved interference patterns were generated. Finally, Figure 6.17 which shows the intensity distribution over the rib-tissue plane, indicates that the rib field minimization approach employed in the synthesis process achieved minimal rib power deposition. The relative rib power ratio c, in this case is about 7.4 %. 6.6.2 1 MHz Operating Frequency with no Attenuation Utilizing higher frequencies in HIFU treatment is motivated by the need to obtain higher power deposition at the tumor location coupled with shorter periods to achieve the desired thermal coagulation of malignant tissues. This motivates the need to study focusing as a function of the operating frequency. To address this issue, single focus scheme is carried out using 1 MHz excitation frequency. The input data for this simulation are: * Operating frequency is 1 MHz which corresponds to a 1.5 mm wavelength in the water, fat, viscera and liver. 180

6... 6.................. 35...................................................................................... 50 100 150 200........ 250.!)2....................................... a0 0 50 100 150 200 250 Array Elements 0 50 100 150 200 250 Array Elements Figure 6.13: Complex excitation of the feed array for the two focus example at 500 kHz. 181

0.8 S 0.6. N 0.4 0 0.2' 250 i......... Q..... *.::.::....... ^..::....... * '.:.: *. -::.:.... * s....... -. -....... -.. v. -.. -......... A.,.,,.. A.,*. -......,:., -:.::.: I' 200 300 150 200 1820 a' 100 50 x-axis (mm) y-axis (mm) ample at 500 kHz. 182

* 480 planar, uniformly spaced excitation elements with interelement separation equal to a half wavelength. * Height of the array elements from the rib plane is 34 A. * Number of ribs is 6. * Rib diameter is 12 A. * Aperture (tissue between ribs) width is 14 A. * Rib attenuation is 1.64 Np/(cm MHz) which corresponds to 0.366 Np/A. * Sound speed in the ribs is double that in the soft tissue as indicated in Table 2.1. * Rib density is approximately one and half times that of the soft tissue as shown in Table 2.1. * Proposed target location is (120 mm, 150 mm) or (80 A, 100 A). It should be noted that for both the 500 kHz and the 1 MHz simulations, the interelement separation of the feed array was fixed at half wavelength to enable fair comparisons between the two simulations. The results for the single focus case at 1 MHz are given in Figures 6.18 to 6.22. Figure 6.18 shows the complex excitation vector of the feed array and as shown, feed distribution over the feed array is similar to that obtained at 500 kHz. Figures 6.19 shows a two dimensional intensity pattern inside the rib cage showing an improved focusing scheme compared to that at 500 kHz. This is mainly due to the increased 183

gain at the focal point due to the increase in the excitation frequency. Figures 6.20 and 6.21 display the a focal plot cut and a contour plot (25 % contours), respectively. They show that at 1 MHz, the intensity level at the focal point is precisely fulfilled in both magnitude and location coupled with acceptable interference patterns in the interior of the rib cage. Figure 6.22 shows the intensity pattern over the rib-tissue plane and it can be seen that the intensity levels over the rib surfaces are lower than those over the intercostal spacings. The relative rib power cr is 9.15 % in this case. However, some ribs are exposed to higher intensity levels compared with the 500 kHz excitation case due to the increase of the scattered power. In summary, the 1 MHz test case leads to two important advantages over the 500 kHz case. The first is the improved focusing scheme with better interference patterns and the second is the higher power deposition at the tumor location(s). However, this operating frequency can not be utilized in treating deeply seated tumors because of the higher losses/attenuation within the medium, as will be seen in the following section. 6.6.3 Effects of Tissue Attenuation on the Focusing Scheme To address the effect of soft tissue attenuation on the focusing scheme, the single focus example discussed previously is reconsidered. For the worst case scenario, attenuation of 1 dB/(cm MHz) may be introduced for all soft tissue layers. Two excitation frequencies, 500 kHz and 1 MHz are tested with a unity intensity focus desired at the point (120 mm, 150 mm). 184

Results of the synthesis approach are displayed in Figures 6.23 to 6.27 for excitation at 500 kHz. Figure 6.23 shows the complex excitation vector of the feed array. As shown, the pattern distribution over the feed array elements is identical to that given in Figure 6.8 except that higher amplitudes are observed over the array elements. This scaling is expected because of the increased array power required to compensate for power loss within the soft tissues. The phase distribution over the feed array in the presence of losses is similar to that obtained when the lossless case was considered. Regarding the feasibility of focusing, Figures 6.24 and 6.25 show that acceptable focusing scheme with relatively good interference patterns is achieved. However, intensity levels over the virtual plane as well as below the rib surfaces are higher than those observed for the lossless case and this is clearly displayed in the intensity pattern given in Figure 6.27 and the contour plot shown in Figure 6.26. Although the interference pattern is deteriorated, the synthesis algorithm tends to deposit the desired levels at the foci locations and this is shown explicitly in the focal plane cut in Figure 6.25. As expected, higher field intensity patterns are observed over the rib-tissue plane when compared to those in Figure 6.12. However, the relative rib power Er is still around 6.9 % close to the lossless case. It is worth noting that a small fraction of the total power is deposited at the focal location. In this example, the ratio between the focal spot power and that at the frontal plane of the rib cage is around 9.28 %. This ratio is a function of the focal location and tissue attenuation. The results for the 1 MHz case are presented in Figures 6.28 to 6.32. Figure 6.28 shows the complex excitation vector of the feed array and as stated for the 500 kHz 185

excitation, it has the same scaling feature compared with the corresponding lossless case. The interference patterns shown in Figure 6.29 are too deteriorated and become unacceptable for treatment. This can be also deduced from Figure 6.31 where high interference patterns are also present in the interior of the rib cage. Although the prespecified intensity level was accurately achieved at (120 mm, 150 mm) as shown in Figure 6.30, the higher levels within the soft tissue region cause significant temperature increases in this area. As a result, the advantage of achieving well behaved patterns inside the rib cage disappeared. Figure 6.32 demonstrates the intensity distribution over the rib-tissue plane and as can be seen, intensity levels over the rib surfaces are dramatically higher than those observed in Figure 6.22 when the tissue loss was ignored. However, the relative rib power r, is around 8.5 % which is close to that with lossless tissues. The ratio between the focal spot power and that at the frontal plane of the rib cage is around 2.15 %. 6.7 Discussion and Conclusions The main goal in this chapter was to propose, implement and test a more realistic model for the problem at hand. First, a model with improved rib geometrical features and ultrasonic characteristics was proposed as shown in Figure 6.1. To analyze this model, exact field prediction techniques [22, 39, 81] based on eigen field expansion were developed and the computational formulation was implemented. Important conclusions were drawn from the realistic model simulations. In particular, focusing HIFU beams even when realistic modeling of the rib cage is considered is still feasi 186

ble. This conclusion is valid irrespective of the computational model or the employed field computational techniques. In addition, this chapter demonstrated the robustness of the two step synthesis algorithm developed in Chapters 5, since by using this approach, well behaved interference patterns concurrently with minimal rib intensity levels were achieved. However, regardless of the employed field computational methods, similar pattern features and desired intensity levels were observed when the VAR approach was employed for pattern synthesis in the examples discussed in Chapter 5 and 6. Introducing tissue attenuation resulted in limitations for the actual treatment since the focusing scheme was deteriorated when the frequency was increased due to the presence of higher interference patterns. Therefore, depth of penetration (target location) should control the choice of the operating frequency. This characterization is addressed and discussed in the following chapter in the context of a design approach to allow for an optimal treatment plan. Finally, regarding the computational cost of the employed exact eigen expansion [17, 22, 39, 71, 81], it should be noted that with the huge size of the computational domain, it was necessary to perform both synthesis and scanning problems in an efficient way. For example, the Bessel and Hankel functions were tabulated before starting the synthesis and table look up was carried out when needed. Also, for calculation accuracy, sufficient terms (around 100 term) from the series expansion were considered to ensure that convergence is achieved. 187

1,^ 0.8 -_ () E 0.6 N X 0.4 Z Z 0.2 ^ 0.8... -... c 0.6 0 Z 0.2.... 0 50 100 150 200 250 x-axis (mm) Figure 6.15: Focal plane cut for the two focus example at 500 kHz. 188

240 220 200 180 160 E 140 CS3 ' 120 X x 100 80 60 40 20 -..:.................................... -...................................................... ^................................... -. ~- i...............____ _....^ ^.................................... -.... 9~................................................. 9~ 4i'....Q ~.................................................... ~~~.........,:................. ~.... '................................................................................................~.................-.... o; <> ^, I I I ^...................................................... i i i i 50 100 150 200 y-axis (mm) Figure 6.16: Intensity pattern contour plot (25 % contours) two focus example at 500 kHz. 250 for the 189

01 1 L I I L ___- ___ -50 0 50 100 150 200 250 Rib-Tissue axis (mm) Figure 6.17: Intensity pattern over the rib-tissue plane for the two focus example at 500 kHz. 190

5i 4-' 0 0 50 100 150 200 250 Array Elements 300 350 400 450 500 A I I I I I I I I I I It I I I I dII I I I I1 I I I I U) 4 U1) CO) co ( 0~ U)0 U1) tE,.% -2 '1 1 11 1 1.1III II. I 1 - - - I I I'.. 1. -I I. III I I A A. 0 50 Figure 6.18: 100 150 200 250 300 Array Elements Complex excitation of the feed example at 1 MHz. 350 400 450 array for the single 500 focus 191

1 0.8. ~ 0.6. "o N I 0.4. o z 0.2.. Z O> 250 200... ~....... ~.~,....." ".. 25 300 150 200 100 - - 150 50 100 50 x-axis (mm) 0 y-axis(mm) Figure 6.19: Two dimensional intensity pattern for the single focus example at 1 MHz. 192

1 0.9 0.8 0.7 1 0.6 C 0.5 N E o 0.4 z 0.3 0.2 0.1 I a - I............................................................................................................................................................................................................................................................................! r-.............................. I --.I................................................................... . I........-........... -........... -........... -. I.........-............ -........I........... -............ n - % — --- 0 Figure 6.20: 50 100 150 200 250 x-axis (mm) Focal plane intensity cut for the single focus example at 1 MHz. 193

240 220.. 200............................................................................. 180............................................................................. 1 6 0................e ~ -. -................................................................. E140............................................................................. co 80F --- -- ---- - -.............................................. c 6 12 0................................. I.: 100............................................................................. 6 0 -.................................................................... 0....... 50 100 150 200 250 y-axis (mm) Figure 6.21: Intensity pattern contour plot (25 % contours) for the single focus example at 1 MHz. 194

0' -50 0 50 100 150 200 250 Rib-Tissue axis (mm) Figure 6.22: Intensity pattern over the rib-tissue plane for the single focus example at 1 MHz. 195

0 50 100 150 200 250 Array Elements I I I ca2 -4 0 50 100 150 200 250 Array Elements Figure 6.23: Complex excitation of the feed array for the single focus example at 500 kHz with lossy tissue. 196

1.4 1.2 U):tC I 0.8s (1) N 0.6 E CO z 0.4 0.2. 0 250 ' '' ' ' ~ ~ ~ ~ ~ ~. ~~....= ~.. ~.~,.~.... ~....:.:.,: ' *.................. 200 150 200 100 150 50 100 50 x-axis (mm) 0 y-axis (mm) Figure 6.24: Two dimensional intensity pattern for the single focus example at 500 kHz with lossy tissue. 197

1.4 - 1.2 1-...................................................................I.............1 E E x 0.8 - ~ 0.6 H 0.4..... I I.............................I..................................-.....................I... IXI 0.2 -..............: —.. L A 1 It. - -- - - 0 50 100 150 y-axis (mm) 200 250 Figure 6.25: Focal plane intensity cut for the single 500 kHz with lossy tissue. focus example at 198

220 2 2 0 -............................................ 200 l 1 8 0.............. 1 6 0............. XE140 L * - - -. 1 2 0 -..................................................... e............ ~.......................d 100.... 80 -------— l 60 L><. <.. 40 W 50 100 150 200 250 y-axis (mm) Figure 6.26: Intensity pattern contour plot (25 % contours) for the single focus example at 500 kHz with lossy tissue. 199

0 o '1 ' '' -l - ' - -50 0 50 100 150 200 250 Rib-Tissue axis (mm) Figure 6.27: Intensity pattern over the rib-tissue plane for the single focus example at 500 kHz with lossy tissue. 200

t) a) a 20 0. E <15.> a) no 0 50 100 150 200 250 300 350 400 450 500 Array Elements U) a. a) -Z -2 -- -- 1... — I 0 50 100 150 200 250 300 350 400 450 500 Array Elements Figure 6.28: Complex excitation of the feed array for the single focus example at 1 MHz with lossy tissue. 201

1.8 1.6 1.4, 1.2 c 0 E 1 0 z I 0.8 O0 250... ~. ~. ~~ a a ~~ a ~r a 2 '~ ~D.'.':..~....:. a.,:;r 300 n r x-axis (mm) y-axis(mm) y-axis (mm) Figure 6.29: Two dimensional intensity pattern for the single focus example at 1 MHz with lossy tissue. 202

1 0.9 0.8 0.7 1 'a 0.6 c ~:C 0.5 N E o 0.4 z 0.3 0.2 0.1 I I -................ -........... -..... -................. -...................... -,...................................................................................................................................................................................................................................................................................................................................................................................... -................. I -- - - I.......... -......... I......... I. I........ 0..................... I.... - rA l l -.L -M - - 0 Figure 6.30: 50 100 150 200 y-axis (mm) Focal plane intensity cut for the single focus 1 MHz with lossy tissue. 250 example at 203

f3A A;LLatl III I I ' --- —1- V 2 2 0.......................................................................................... 20O 160............. 1 2 0..-........................... -=................ 80 coo 6 0......................................................... 20 50 100 150 200 250 y-axis (mm) Figure 6.31: Intensity pattern contour plot (25 % contours) for the single focus example at 1 MHz with lossy tissue. 204

0.9 i 0.8 D...-...........2 0.7.......... 0.7..-............................................................ 0.6.............. C 4................. -50 0 50 100 150 200 250 Rib-Tissue axis (m m) Figure 6.32: Intensity pattern over the rib-tissue plane for the single focus example at 1 MHz with 0ossy tissue. 205

CHAPTER 7 Design and Computational Aspects In the previous three chapters, phased array pattern synthesis algorithms for HIFU treatment of liver tumors have been formulated, implemented and validated. Regardless of the employed synthesis algorithm, practioners who use phased arrays need to design an optimal treatment plan for liver tumors. This means that for a specific tumor, treatment frequency should be optimally selected. In this chapter, parameters affecting the focusing scheme such as operating frequency, maximum penetration depth and soft tissue attenuation are discussed and characterized. In addition, effect of the rib power minimization is carefully addressed and analyzed to demonstrate the importance of the constrained solution obtained in the second step of the VAR synthesis procedure. 7.1 Attenuation Model and Energy Requirements In Chapters 3, 4 and 5, the propagation model was assumed to be lossless with soft tissues speed of sound of 1500 m/sec and density of 1000 kg/m3. As indicated 206

in [77], an appropriate model of the attenuation can take the following frequency dependence form a = ao f, (7.1) where ac, is the attenuation coefficient at frequency 1 MHz, and v is a constant, and according to [77], its value is 1.1. The thermal model described in Chapter 2 will be utilized in all simulations discussed in this chapter. The tumor and normal tissues are assumed to have the same thermal conduction and propagation features when excited by ultrasonic waves. In both malignant and normal tissues, the attenuation coefficient can be considered to be equal [63]. The power deposition within the volume of a certain tumor can be evaluated using the following relation Pt = IKtTh, (7.2) where Pt is the power deposition within the tumor volume, I is a constant that depends on the dimensions and geometrical features of the tumor, Kt is the thermal conductivity and Th is the temperature increase in ~C. Assuming that the selected points inside the tumor correspond to the focal points along the scan trajectory, then the total power within the tumor can be approximated by [26], Nf Pt= QiVi, (7.3) m=l where Qi is the power deposition level at the ith focal point, Vi is the effective volume of the focal spot, and Nf is the number of focal spots along the scan trajectory. To 207

simplify (7.3), the following assumptions may be assumed: * The effective size of the focal point does not change as they traverse the scan trajectory. * Equal power deposition is assumed at each focal point. Using these two assumptions and (7.3), the power deposition at each point can be expressed as Qav= Pt (7.4) -NfV' where V is the effective size of any focal point and Qav is the average power deposition level. After evaluating the average power deposition at each focal point, the intensity relations have to be invoked. Using (2.13), the spatial-peak time-average intensity at a focal point is given by Qav ISPTA = (7.5) 2cC where a is the tissue absorption coefficient which is assumed to be equal to the attenuation factor as indicated by [63]. To get the spatial peak pulse average intensity at the focal point, the duty cycle ratio - should be evaluated and since, equal dwell time is assumed at each focal point, then T = NfT. (7.6) 208

Therefore, the spatial-peak pulse-average intensity at the focal point defined by (2.14) can be expressed in the following form ISPPA = NQav (7.7) 2a According to (7.7), the average power deposition within the tumor volume is related to the applied ultrasonic intensity of the incoming beam and the tissue absorption coefficient at this location. In the following section, the treatment frequency is chosen such that it maximizes the power deposition at the focal points which represent the tumor(s) locations. 7.2 Optimal Power Deposition Because of the need to generate a lesion in the tumor vicinity, HIFU beams are focused at the tumor location(s) for a certain period of time to achieve the necessary thermal necrosis. As displayed in Figure 2.4, the time required to reach this thermal stage decreases with the increase of the applied ultrasonic power. Therefore, for more power deposition at the focal spots, higher frequency ranges for HIFU are recommended. This increases the path attenuation [26, 53, 54] and thus higher interference patterns are expected inside the rib cage especially in the near bone area and over the rib surfaces. This fact was demonstrated in the previous chapter when the operating frequency was increased from 500 kHz to 1 MHz. The intensity levels in the near bone area at 500 kHz were significantly lower than those generated at 1 MHz as displayed in Figures 6.24 and 6.29. Thus, it is required to choose the treatment frequency based 209

on the focal spot location and the acceptable intensity levels in the near bone vicinity. The following factors affect the HIFU treatment: 1. The operating frequency f, which plays a major role by directly affecting the power deposition at the focus location. It should be chosen in such a way that ensures higher power deposition at the tumor location without causing excessive overheating at undesired locations. 2. The maximum penetration depth d which is the focal spot location (target location). This distance affects the interference patterns inside the rib cage since, at a given operating frequency, going deeper inside the rib cage increases the overall path attenuation and hence interference patterns in the interior of the rib cage are significantly deteriorated. 3. The tissue attenuation ac is assumed, according to [63, 77], to have the frequency dependence described by (7.1). The power deposition at the focal point(s) is expressed in terms of the operating frequency, tumor location, and tissue attenuation, then standard maximization technique will be employed to extract a simple criterion that controls the selection of the operation conditions according to tumor and tissue characteristics. From (7.7), the average power deposition at a point located at a distance d under the ribs can be expressed as Qav(d) = ISPPA(d). (7.8) Nf 210

Assuming plane wave unfocused system, the power deposition of the ultrasonic beams propagating from the frontal plane of the rib cage to the tumor location can be expressed as 2cQav(d) = -Ioe-2 (7.9) Nf(7.9) where Io is the intensity level at the frontal plane of the rib cage. Substituting from (7.1) into (7.9) yields 2aoIo v 2,of 'd Qav(d) = f-2f. (7.10) Equation (7.10) can be reduced to Qav(d) = Qofe-2,of'd, (7.11) where Qo = 2 (7.12) It is now required to select the operating frequency that maximizes the power deposition at the tumor location, i.e. select f that maximizes Qav(d) given by (7.10). The first derivative of (7.10) with respect to the operating frequency f is given by Q(d= -Q[f e-2aodfd(-2vaodfv-1) + e-2afdvfv-l] (7.13) af Equating the first derivative with zero yields, 2aofv d = 1 (7.14) or f= 2d (7.15) 211

Therefore, the optimal operating frequency can be expressed as -1 fopt(d) = (2aod)-v, (7.16) where fopt(d) is the optimal treatment frequency that maximizes the power deposition level when a certain treatment is desired at distance d under the skin. From (7.16), it is concluded that to reach deeper tumors, the operating frequency has to be reduced. It should be noted that although the formulation was carried out and implemented for the single focus scheme, it can be extended to cover the multiple focus one. This can be directly achieved by replacing Nf by Nj, where M is the number of foci. At the frequency given by (7.16), the value of the power deposition can be deduced by substituting from (7.16) into (7.11). This substitution yields 1 Qavopt (d) = Qo e-, (7.17) 2a0d where Qavop is the power deposition level at the optimal treatment frequency. It is worth noting that this value of the power deposition level represents the maximum that can be achieved at a certain location d. For convenience purposes, a power deposition index may be defined by normalizing the power deposition level in (7.11) with respect to QO. This normalization introduces the parameter tu(d) = fve-2aof"d, (7.18) where 11(d) is the power deposition index when performing HIFU treatment at distance d under the skin. From (7.18), the power deposition index for the optimal case 212

is given by 1 -1 =,opt( ) 2aod (7.19) This equation can be rewritten in terms of the optimal treatment frequency as Popt(d) = fopt(d)-1. (7.20) 7.3 Design Curves In this section, curves relating the power deposition level at the focal plane to the operating frequency and penetration depth are generated and utilized to design the appropriate treatment plan. In practice, liver tumors exist at depths from 6 cm to 15 cm under the skin. A set of design curves is generated by scanning the operating frequency from 100 kHz to 1.6 MHz at each specific location to generate graphs that can be utilized in selecting the optimal treatment planning. Figures 7.1 shows the variation of the power deposition index with both the operating frequency and depth when the factor v is chosen to be 1.1. To study the effect of the operating frequency on the heat deposition index, the distance d is fixed and the frequency is scanned from 100 kHZ to 1.6 MHz and another set of graphs is generated and plotted in Figure 7.2. As shown in this figure, the power deposition index,a has a certain maximum at a frequency value that matches the optimal frequency deduced from (7.16). Figure 7.3 shows the same data at distances from 11 to 15 cm which indicate that significant degradation in the power deposition level occurred when the frequency and/or distance increases. 213

From the two sets of design curves generated and plotted in Figures 7.2 and 7.3, it is observed that the power deposition index,u is fairly flat around the optimal value. This motivates the thinking to increase the operating frequency allowing certain degradation in the power deposition level at the tumor location. For example, if the tumor resides at 10 cm under the skin, then according to (7.16), the optimal operating frequency is 450 kHz and, from (7.19), the corresponding power deposition index is 0.1585. If operating frequency is increased to 550 kHz, the power deposition index will be 0.1561. This means that increasing the operating frequency by approximately 20 % lead to a power deposition reduction of less than 2 %. It is then recommended to increase the operating frequency by a certain percentage concurrently with keeping the power deposition at the tumor location within certain range from the maximum value given by (7.19). To formulate this problem, the operating frequency is related to the optimal one by the following equation rf(d)= f (7.21) fopt (d) where rff (d) is the ratio between the utilized operating frequency and that corresponding to maximum power deposition at the tumor location. Similarly, the power deposition index at a given frequency is related to its value when the maximum heat deposition is intended by using the following equation Af(d)= - (d) (7.22) where Af(d) is the power deposition index ratio and it indicates the operating conditions deviation from the optimal situation. It is clear that at optimal frequency of operation, both f (d) and Af(d) have the unity value. 214

Using (7.16), (7.18), (7.19), (7.20) and (7.21), it can be shown that rf(d) and Af(d) can be related by the following equation Af(d) = F(d)e-()e. (7.23) For a certain degradation in the power deposition, the factor Af(d) is specified and the solution of (7.23) should lead to operating frequency. Thus, (7.23) can be rewritten on the following form a, = xe-X, (7.24) where ar is a constant related to the required ratio of the power deposition index, and is given by ^f(d) ar = (), (7.25) e and x is the parameter we solve for given by x = rF(d). (7.26) Practically, the value of x is close to unity and in order to be able to find an appropriate solution to (7.24), a new parameter y is defined according to the following equation y = x-, (7.27) where y is smaller than unity. Substituting from (7.27) into (7.25) gives the following equation (y + l)e-Y = b,, (7.28) 215

where br is a constant given by br = are. (7.29) The series expansion for the exponent eY for small values of y can be approximated by the following e-Y 1 i - y. (7.30) Substituting with this expansion in (7.28) yields (1 + y)(l - y) = br,. (7.31) Equation (7.31) can be reduced to 1 - y2 = b, (7,32) which results in the following equation, Y = ~V1-b, (7.33) According to (7.27), the corresponding solutions for x will be x= V1-b+1, (7.34) or x = -V+/1 - +1. (7.35) Since the goal is to increase the operating frequency, then, the value of x given in (7.34) is chosen. Substituting from (7.25) and (7.26) in (7.34), the following relation can be deduced r(d) = / - Af(d) + 1, (7.36) 216

from which, the proposed operating frequency that can be employed without significant degradation of the power deposition is given by the following relation f fopt(V/ - A(d) + 1)>. (7.37) which means that if 10 % was allowed as degradation in the heat deposition level, then the treatment power deposition level is 90 % of its maximum and as a result, the value of Af(d) is 90 %. Substituting into (7.38) gives f = 1.32 fpt (7.38) This means that a small degradation in the power deposition level leads to a significant increase of the operating frequency. Figure 7.4 indicates the gain in the operating frequency versus the degradation in the power deposition level at the focal point. 217

~~ ~ ~I ~.~. I 0.05 2 1.5 10;, ^-0. 0.5 14 16 ~ Operating Frequency MHz Focal Distance d Figure 7.1: Variation of the deposition index y with both the focal distance d and the operating frequency in MHz. 218

35, --, I. - 0. 0.2 x 0. _o c o.0.1 ' d=6 cm::: d=7 cm 3........................................................... d= 8 cm. 3 -:... - I d=8 cm.:. * K d=9 cm d=10 cm 5 2 5 0. 0.0 0 0.2 0.4 0.6 0.8 1 Operating Frequency MHz 1.2 1.4 1.6 Figure 7.2: Deposition index as a function of the treatment frequency for a set of focal distances ranging from 6 cm to 10 cm. 219

n 16 rII-I 0.14 I I I..:................. ~- * d=11 cm ---- d=12 cm...... I- d=13 cm x --- d=14 cm d=15 cm.......:..................... 0.121 0.1 - x 'a C 0 ).08 - 3.06 - 0.04 F 0.02 F I I I I I --- -- -- A n t I I ~ 0 0.2 0.4 0.6 0.8 1 Operating Frequency MHz 1.2 1.4 1.6 Figure 7.3: Deposition index as a function of the treatment frequency for a set of focal distances ranging from 11 cm to 15 cm. 220

/(U 60F 0.c8 C a) =0.50 0) c 0 C 0 0 0) CD........................................................................................................................................................................................................................................................................................................ 20 if. I I I I 0 0.05 0 0.05 0.1 0.15 0.2 0.25 Degredation in the power deposition index 0.3 0.35 0.4 Figure 7.4: Percentage increase in the treatment frequency tion of the allowed deterioration or degradation position index p. as a funcin the de 221

7.4 Effect of Rib Power Minimization According to the VAR synthesis approach presented in Chapter 5, the second step of the synthesis process explicitly accounts for the minimization of intensity levels and power deposition over the ribs. An important question arises, What would be the situation if minimization of the rib power was not considered? In other words, what is the effect of the constrained optimization on field patterns inside the rib cage as well as over the rib surfaces? To answer these questions, the following two cases are considered: * Considering minimization of the rib power in the second step of the two step VAR synthesis. * Ignoring the minimization of the rib power levels in the synthesis. Thus, no constraints are enforced while computing the feed array that is responsible for inducing the desired virtual sources over the intercostal spacings. To test these two cases, the following was assumed: 1. A single focus pattern synthesis is performed. 2. Soft tissue attenuation is neglected. 3. Focal intensity of unity value is desired at (120 mm, 150 mm). 4. Operating frequency is set to 500 kHz. 222

In the following, complex excitation vector of the feed array, interference patterns, focal plane intensity pattern and rib-tissue intensity patterns are compared for both cases. * Excitation vector of the feed array: Figures 7.5 and 7.6 show the complex excitation vector of the feed array for both cases. From these figures, the following can be observed: 1. Figure 7.5 shows the excitation vector over the feed array when rib power minimization was applied and as can be seen, the pattern amplitude of this vector is similar to that resulting in Chapter 6 when the minimization of the rib power was invoked. 2. Figure 7.6 shows the same excitation vector when no rib power minimization criterion was enforced and it can be seen that the amplitude distribution over the feed array follows the same pattern observed when the HRPO technique is employed for pattern synthesis. 3. Figure 7.6 implies also that in the absence of rib minimization constraint, array elements which have no direct access to the focal point (tumor location) are excited with low power (amplitude) and this explains the semi periodic behavior of the velocity amplitude distribution. 4. Figure 7.5 displays the case when rib power minimization was considered and it is observed that the periodic nature of the feed array amplitude vanishes since elements with no direct access to the focus are excited in such a way so that signal cancellation is achieved over the rib surfaces. 223

5. The array efficiency in the case when rib power minimization was considered is approximately 10.33 %, while this efficiency was about 22.92 % with no rib power minimization. This was also observed when comparisons between the HRPO and VAR approaches were performed in Chapter 5. * Interference patterns: Figures 7.7 and 7.8 show the two dimensional intensity patterns for both cases. The following can be concluded: 1. Similar interference patterns are observed when single focus scheme is considered. This feature was also observed in Chapter 5 when the VAR and HRPO techniques were compared for the single focus case. 2. This similarity is due to the fact that the difference between the two cases appears in the rib-tissue intensity fields since the PI approach utilized for both cases tends to minimize power deposition everywhere inside the rib cage. 224

/) 1.5..... a: 0 0 50 100 150 200 250 Array Elements u)2...... co IL-4 0 50 100 150 200 250 Array Elements Figure 7.5: Complex excitation vector of the feed array with rib power minimization. 225

0.5 (o0.4......... 0 *60n:3-. 0.3............. E ~0.2..... co |o.3- i: 0:. Cc 0.i1.....w-, 0 50 100 150 200 250 Array Elements 4 -41 0 50 100 150 200 250 Array Elements Figure 7.6: Complex excitation vector of the feed array with no rib power minimization. 226

1.4 -1.2-. 1 c | 0.8s J 0.6 -Co - 0.4 -0.2 O 0> 250 '' ''. A....................... s............................ *.. *. *,.:........................... *... * *.:..................................................... z 9.,'""""" I' 200 300 150 200 100 150 100 50 50 x-axis (mm) y-axis (mm) Figure 7.7: Two dimenona intensity pattern when rib power minimization is considered. 227

1.4 1.2. 1 _' 0.8 - 'a 0.6. Z 0.4 -0.2 O. 250 200 ~ ' ~ ' ~.............. ~ ~ ~. ~. ~. LJ;n......,......................................, b W *: -....:........................................ 5. *.......................... i N.:... -...... 250 300 150 200 100 150 100 50 50 x-axis (mm) 0y-axis (mm) Figure 7.8: Two dimensional intensity pattern with no rib power minimization. 228

* Focal Plane Intensity Distributions and Contour Patterns: The focal plane and contour cuts are displayed in Figures 7.9 to 712 and from these figures, the following can be deduced: 1. There is a similarity between the focal plane intensity distribution for both cases as shown in Figures 7.11 and 7.12. This can be attributed to the fact that the minimization criterion is enforced over the rib-tissue plane. 2. Regarding the contour plots, similarity is also observed as shown in the 25 % contour plots in Figures 7.11 and 7.12. This similarity is due to the same fact explained in 1. 229

0 50 100 150 200 250 x-axis Figure 7.9: Focal plane intensity cut with rib power minimization. 230

* Rib-Tissue Pattern Distribution: Figures 7.13 and 7.14 display the intensity pattern over the rib-tissue plane for the constrained and unconstrained cases. From these figures, the following can be deduced: 1. It is clear that the intensity levels over the rib surfaces are significantly higher when no rib power minimization was enforced in the second step of the VAR approach. This is clearly demonstrated in Figure 7.13 where the rib power percentage cr is around 6.4 %. On the other hand, this percentage dramatically increases to 13.78 % when no rib power minimization is invoked. 2. Due to the realistic modeling of the propagation mechanism that accounts for tissue parameters, rib shapes and attenuation, ranges for or have dramatically increased compared to those obtained in Chapters 4 and 5. In these chapters, a simplified model with thin solid ribs was assumed with no partial transmission through the ribs. 7.5 Discussion and Conclusions This chapter dealt with design issues and computational aspects for deep HIFU treatment inside the rib cage. The design issues introduced the main parameters and factors affecting the focusing scheme and it was demonstrated that a treatment plan can be designed based on a simple set of design curves. The loss in the soft tissue 231

layers (fat, viscera and liver) turned out to be the main parameter when selecting the operating conditions of the HIFU treatment. As illustrated by the results, focusing up to 10 cm or 12 cm depths can be achieved concurrently with well behaved interference patterns. However, more penetration depths lead to higher intensity levels over the frontal plane of the rib cage. These levels, especially in the near bone region, produce excessive amount of heat in the near bone area and over the rib surfaces. The effect of rib power minimization over the field maps within the computational domain was also addressed. The examples presented in this chapters indicated that rib power minimization dramatically affects the fields over the rib surfacers as well as the linear array feed. Using the minimization, the relative rib power was significantly reduced with almost no change in the intensity pattern over the focal plane or inside the rib cage. 232

1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 50 100 150 200 25C x-axis Figure 7.10: Focal plane intensity cut with no rib power minimization. 233

Z4L) I I 220 - 200 - 180.... 160. 140. 120 -100 - 80 - 60 - 40 - 20 - Figure 7.11:................. _. _ _ _..... 50 100 150 200 250 y-axis Intensity pattern contour plot (25 % contours) with rib power minimization. 234

240 220 200 180 160 140.a 9 120 x 100 80 60 40 20 -......................................................................................................................................................................................................................................................................................................................................................................................r................................................. l 50 100 150 200 y-axis Figure 7.12: Intensity pattern contour plot (25 % contours) rib power minimization. 250 with no 235

100 Rib-Tissue axis (mm) Figure 7.13: Intensity pattern over the rib-tissue minimization. plane with rib power 236

0.015.......................................... o -50 0 50 100 150 200 250 Rib-Tissue axis (mm) Figure 7.14: Intensity pattern over the rib-tissue plane with no rib power minimization. 237

CHAPTER 8 Conclusions and Recommendations for Future Work This chapter summarizes the main conclusions and contributions of the thesis work. It also presents some recommendations for future work. 8.1 Chapters Contributions The problem of focusing high intensity ultrasonic beams through the rib cage was considered in this thesis. Phased arrays were utilized for focusing because of their ability to produce efficient heating patterns capable of achieving the desired therapeutic levels required for HIFU treatments. The main components of the research conducted are as follows: 1. Testing the ultimate or ideal case of focusing in homogeneous media using multiple transducer phased array applicators. This study was performed in Chapter 3 to show the best case for focusing when no obstacles exist within the medium. 238

The focusing in homogeneous media gives a brief introduction to the basic synthesis algorithm and parameters. 2. A hybrid ray physical optics (HRPO) computational technique was developed and implemented in Chapter 4 to show that focusing through strong inhomogeneities such as the rib cage is feasible. This model was validated with an experiment and good agreement between measured and simulated data was found. The HRPO algorithm proved the following: * High frequency ray techniques are more attractive for modeling the problem of focusing through the rib cage. Numerical methods such as the finite element method have huge CPU and memory requirements. Also, such techniques can not be easily employed for patterns synthesis. For a problem with a region of interest that spans hundreds of wavelengths, the accuracy of these numerical methods is not expected to have significant difference from that of high frequency schemes [17, 22, 71, 81]. * As a minimum-norm least-squares solution, the PI pattern synthesis technique coupled with the HRPO algorithm deposits the prespecified power at the desired foci locations concurrently with acceptable power deposition everywhere within the computational domain which is required for HIFU treatments. * The HRPO synthesis algorithm [13, 14, 15] guides most of the ultrasonic energy through the intercostal spacings between the ribs and hence acceptable power deposition levels were observed over the rib surfaces. However, 239

the HRPO technique has the following three drawbacks: (a) It does not account explicitly for power minimization over the rib surfaces which may lead to excessive heat accumulation over the rib surfaces. (b) The PO technique employed to predict the ultrasonic fields in the interior regions of the rib cage has poor prediction accuracy in the near zone (i.e. near bones). (c) The computational cost of the HRPO technique in constructing field maps in the interior region of the rib cage is unaffordable. 3. Motivated by the aforementioned drawbacks, another synthesis approach was developed and implemented in Chapter 5. Namely, the two step hybrid virtual array ray (VAR) technique was formulated and validated by comparing it to the HRPO technique. This method offers the following advantages: * It accounts for the geometrical features of the rib cage by making efficient utilization of the acoustical windows represented by the intercostal spacings between the ribs. * It takes into account minimization of the incident power deposited upon the rib surfaces. This minimization is considered as part of the main design by adding a constraint in the synthesis approach. * Computationally, it is significantly less intensive than the HRPO technique. As described in Chapter 5, the CPU cost ratio between the two techniques is at least 2000. Thus, after computing the feed array using 240

either HRPO or VAR, fields scanning was performed using VAR. The field maps constructed by the VAR and HRPO scanning techniques are identical for a given array excitation and this gives an additional validation for both approaches. * The VAR technique removed field uncertainty associated with the PO approach especially in the near bone area. 4. A realistic geometrical model and exact field computational techniques were proposed and analyzed in Chapter 6 to test the feasibility of focusing inside the rib cage with the ribs considered as lossy elastic cylinders and the attenuation of all layers (soft and solid) was taken into account. By analyzing this model using exact eigen expansions [17, 22, 71, 81], it turned out that the synthesis algorithms are still capable of generated well behaved interference patterns. Testing this realistic model with accurate computational techniques lead the following conclusions * demonstrated the robustness of the developed synthesis techniques even with the inclusion of the rib and tissue shapes and parameters. The field patterns generated using simple computational models did not differ much from those generated with more realistic geometrical features and more accurate field prediction techniques. This lead to the conclusion that regardless of the rib shapes and features, well behaved interference patterns can be generated in the presence of these inhomogeneities. * Using the geometric and computational model developed and implemented 241

in Chapter 6, an optimal treatment plan was proposed in Chapter 7. The parameters that affect this plan are the attenuation in soft tissues, operating frequency, penetration depth and power deposition levels at the tumor location(s). Design curves assuming plane wave unfocused system were generated to obtain the optimal choice of the operating frequency that maximizes the power deposition at a given location. The optimal frequency is a function of tissue attenuation (absorption coefficient) and penetration depth. Curves for optimal treatment planning were generated and discussed. In Chapter 7, the treatment frequency was allowed to exceed the optimal one because of the need to increase the treatment frequency to achieve faster tissue coagulation. The percentage of increase is controlled by the desired degradation for the power deposition at the tumor location. To conclude, the basic contributions of this dissertation are: * A hybrid ray physical optics (HRPO) computational technique was developed, implemented and tested to focus high intensity ultrasonic waves in the presence of strongly scattering obstacles such as the rib cage. The HRPO was experimentally validated. * A two step virtual array ray (VAR) technique was formulated, implemented and tested to focus high intensity ultrasonic waves in the presence of strongly scattering obstacles such as the rib cage. This technique takes into account the minimization of power deposition over the rib surfaces. 242

* Exact computational methods based on eigen function representation for rib scattering were implemented and tested. The implementation accounted for the geometrical features of the rib obstacles and tissue parameters such as the attenuation. Pattern synthesis was also carried out using the VAR technique. * Design curves for selecting the treatment frequency subject to tissue attenuation and absorption, tumor location and heat deposition level were generated and tested. The conclusions of this thesis can be summarized in the following points: * Focusing high intensity ultrasonic beams deeply inside the rib cage is feasible. This conclusion means that HIFU can reach deep regions inside the rib cage and this facilitates the treatment of tissue and organs previously considered non reachable by ultrasonic waves. * High frequency ray techniques provide accurate means for predicting ultrasonic fields within the computational domain for a given excitation or feed. They offer ease of implementation in addition to their affordable CPU and memory requirements. * Optimal field pattern synthesis can be developed and implemented for performing efficient HIFU treatment with well behaved interference patterns and minimal power deposition over the ribs. * Optimal treatment plan can be designed according to patient data, tumor specification and tissue properties (typically from MRI scans). The treatment fre 243

quency can be selected in such a way that achieves maximum power deposition within the tumor area and shorter treatment times. 8.2 Recommendations for Future Work This work demonstrated that HIFU treatments of deeply seated tumors, partially shadowed by the rib cage, is possible using high intensity focused ultrasound. The synthesis process assumed that the rib cage geometry, tumor locations, layer dimensions and characteristics are known a priori. In the future, accurate and efficient imaging techniques should be integrated with the focusing mechanism to identify the geometrical features and layer parameters needed for pattern synthesis. This could be done by using the feed array for both imaging and therapy, in such a context that the following tasks can be performed: 1. First, the feed array is used to illuminate the whole domain by low power wave to identify the geometrical features of the rib cage, rib shapes, locations and dimensions. 2. Use MRI imaging to identify the tumor location and dimensions. 3. Supply the collected geometrical data to the patterns synthesis algorithm to design the focusing array 4. Temperature feedback algorithms [68, 69] can be employed to provide real time feedback on the treatment progress (lesion formation). 244

Instead of using phased array which are fed by single continuous sinusoidal signals, code fed arrays may be employed for sensing purposes. In this case, array elements are fed with different wave forms or signals to achieve optimal imaging or sensing features [40, 72]. Also, several geometric modeling and computational techniques may be suggested for future work. In this thesis, the ribs modeling as lossy elastic cylinders with infinite depth along the axial direction does not match the real case precisely. Therefore, finite ribs with elliptical cross section should be considered. Due to the rib minimization scheme employed in computing the feed array, changing the rib shapes is not expected to have significant effect on the synthesis process. Future work may also include the use of different array configurations to generate the desired focusing scheme. In [26, 27], several feed schemes were characterized and implemented for focusing ultrasonic wave in homogeneous media. To conclude, the proposed future work may include one or more of the following schemes: 1. Improve the accuracy of the computational model. This will require three dimensional analysis for the propagation mechanism and field prediction techniques [42]. 2. Integrate efficient imaging schemes with the therapeutic schemes developed in this thesis. 3. Include application of non-invasive temperature estimation techniques to area time focusing problem. One way of achieving this is by making use of the backscattered ultrasound RF-data [74]. 245

APPENDICES 246

APPENDIX A Minimum-Norm Least-Squares Solution This appendix provides the mathematical details of the minimum-norm leastsquares solution employed for pattern synthesis. A generalized solution of the resulting system is provided followed by a description of the different situations that may appear when physical problems are considered. A.1 Minimum-Norm Solution In all chapters of the thesis, the employed computational formulation yielded the following system of equation = H u, (A.1) where 4 represents the complex field values at a set of M field or control points, u is the complex excitation vector of the feed array and H represents the forward propagation operator from the aperture surface to the field points. The pattern synthesis problem can be summarized as follows: 247

For a given field or control vector 4 E CM, where CM is the set of M-tuples of complex numbers which defines the control point space, it is required to find the minimum energy excitation vector of the feed array u E CK, where C0 is the set of K-tuples of complex numbers which defines the feed excitation space, in such a way that minimizes R given by R = IIHu -_ l\2. (A.2) To solve the aforementioned pattern synthesis problem, it is essential to employ the minimum-norm least-squares solution presented in [26, 27, 56]. This solution yields the complex excitation vector of the feed array u denoted u. An efficient solution for this problem is achieved by finding u e CK that satisfies the following optimization problem mmi J(u) = IIH u- _||2 + ~ U, (A.3) where the parameter ( is referred to as the regularization parameter [78]. This is due to the fact that in some problems, the matrix H is ill-conditioned. To solve (A.3) for u, it is required to fulfill the following equation J(u + 6u) - J(u) = 0. (A.4) The solution of (A.4) is given by u = ((I + H*tH)-lH*', (A.5) where I is the identity matrix and H*t is the conjugate transpose (adjoint) of H. As discussed in [27, 26], the matrix H*t represents the backward propagation operator 248

from the control point space to the source space. The regularization parameter g can be deduced using the method of Lagrange multipliers [62]. A more elegant or generalized solution is obtained using the minimum-norm leastsquares solution developed in [56]. This solution yields the Pseudo-Inverse (PI) of H denoted by HP'. This inverse is called also the generalized inverse and can be obtained using the singular value decomposition (SVD) of the matrix H according to the following equation H =U E V, (A.6) and HP= V S-1 U*t, (A-7) where U is an m x m orthogonal matrix, V is an n x n orthogonal matrix and S is a diagonal elements matrix composed of the singular values of the matrix H. The matrix E can be written as E = diag {Jf, a2.......... O-min(m,n)}, (A.8) where ci is the ith singular value of the matrix H. It should be noted that the singular values are usually ordered in descending order, that is 1 > 2 > **......* mirn(mn) > 0. (A.9) From the singular value decomposition of the matrix H, important information about the system can be deduced such as: 249

~ The System Rank The system rank is determined by the number of nonzero singular values. * The System Condition The condition number of a certain system is defined as the ratio between the maximum and minimum singular values. The information about the condition of the matrix is important before analyzing the stability of the solution. The condition status of the matrix defines the way employed to deal with the system and as indicated in Chapters 3 and 4, well conditioned systems lead to a stable solution when the the PI technique is applied. On the other hand, as shown in Chapter 5, badly conditioned systems require preconditioning as discussed in the constrained PI approach. Using the aforementioned singular value decomposition, the minimum-norm leastsquares solution of (A.1) is given by the following equation ui = HP. (A.10) The solution given in this equation is rigorously unique since it is the only solution that minimizes I H u - 4 12 which has the minimum-norm I ll 12 [18, 55]. It should be noted that this feature is always fulfilled regardless of the rank of the matrix H. If the matrix H has full rank, the PI solution takes the following form: 1. The minimum-norm solution If M is less than K, then (A.10) represents an underdetermined system of equations which has an infinite number of solutions. The PI solution in this 250

case represents the minimum-norm solution given by the following equation i= H*t (HH*t)-~l,. (A.11) 2. The exact inverse solution If M equals K, then (A.10) represents a square system of equations and the PI solution in this case will be the matrix inverse of H which is given by u = H-1. (A.12) 3. The least-square solution If M is greater than K, then (A.10) represents an overdetermined system of equations which has a unique solution that minimizes the function R given in (A.2). This solution is given by = (H*t H)- H*t i. (A.13) It should be noted that the least-squares solution is the the one which satisfies (A.5) with 5 equals zero. In many cases of practical interest, it is desired to evaluate the vector u based on a minimum number of control points. In all chapters of the thesis, the control points represented the focal points and there is a possibility to include also the locations where the power deposition is maintained within certain levels. Therefore, the minimumnorm solution becomes of importance in these practical cases and in addition to the practical issues, the minimum-norm solutions offer the following advantages when utilized for synthesizing focusing schemes required for HIFU treatments. 251

* In the practical situation, it is desired to control the field intensity levels at a set of M control points. If an array of K elements is utilized, then K is greater than M. Thus, the system becomes underdetermined with an infinite number of solutions. The minimum-norm solution is the one that minimizes the energy at the feed array surface concurrently with depositing the desired power levels at the control locations which is required for HIFU treatment. * As described throughout the thesis chapters, the minimum-norm solution reproduces the desired field values at the control field locations with an excellent precision. There is no distortion or scaling in the focal point levels and therefore, it is possible to precisely control the field pressure or intensity levels at the tumor location(s). * The minimum-norm solution allows for optimization of specific parameters such as the array excitation efficiency and the control gain. This is achieved by employing a modified version of the minimum-norm solution which is the weighted minimum-norm. This version can be achieved by replacing the H*t with WH*t in (A.11). Thus, the weighted minimum-norm will take the following form w= W H*t (H W H*t)-l D. (A. 14) where fiw is the weighted minimum-norm solution and W is a positive definite matrix. An application of this weighted PI is the case discussed in [26, 27] where the complex excitation vector of the feed array was computed in such a way that achieves maximum array excitation efficiency. Also, in Chapters 5 and 6 of this thesis, a 252

weighted version of the PI was employed to compute the feed array that minimizes the direct incidence over the rib surfaces. A.2 Computational Aspects of the MinimumNorm Solution The minimum-norm solution requires computing of the inverse of the matrix HH*t in (A.11). This is followed by a matrix multiplication with the operator H*t. An efficient way to perform the two operations is by using normal inverse and matrix multiplication routines. However, singular value decomposition subroutines may be also utilized since the singular values and singular vectors obtained from the SVD analysis are important in the analysis of the matrix HH*t. The aforementioned way of computing the minimum-norm solution is referred to as the direct way. Alternatively, iterative schemes can be employed to solve (A.). These iterative schemes ae auseful when the size of the propagation matrix is too large for direct solvers, [8, 10, 16]. The successive approximation, steepest descent and conjugate gradient techniques [33, 65] are guaranteed to converge to the minimum-norm leastsquares solution obtained by the PI technique [56]. Also, minimal residual schemes such as the Generalized Minimal Residual (GMRES) can be efficiently employed to solve the systems which are badly conditioned such as the systems generated in the second step of the VAR techniques [9, 11, 12]. It is worth noting that preconditioning can be coupled to any of the previously mentioned schemes in order to speed 253

up convergence. Preconditioners vary from the simplest Diagonal Preconditioner (DPC) which has trivial memory and CPU requirements to the more complicated Approximate Inverse Preconditioner (AIPC) which requires intensive memory and CPU demands [16, 65]. 254

APPENDIX B Gaussian Integration Technique Numerical integration schemes are utilized to compute integrals when the closed form analytical evaluation for these integrals becomes non feasible. In Chapter 4 of the thesis, the physical optics (PO) integration method yielded a one dimensional integration with no analytical solution. An efficient and accurate way to calculate the PO integral in Chapter 4 is to use Gaussian Quadrature Integration scheme [1, 19]. The basic formulation and parameters of this method are introduced and summarized in this appendix. B.1 Gaussian Integration Formulation It is required to numerically compute the following integration rb Iab= f(x)dx (B.1) where lab is the integral value, [a, b] is the integration interval within which the integration function f(x) is bounded. For more convenience, it is recommended to 255

transform the integration variable x according to the following equation b-a b+ a x= -2 t (B.2) where t is a new integration variable. Applying this transformation to (B.1) yields the following integration b-a 11 b-a b+a 2 2 (B.3) According to (B.3), the problem now becomes the numerical evaluation of the integral I, = f(alt + bl)dt, 1 (B.4) where b-a al= 2 ' (B.5) and b+a bl = 2 (B.6) The original integration I is related to I, by the following equation I = al I. (B.7) Gaussian integration scheme chooses n points for evaluating the integral given in (B.4), thus n I1 = S w-f(att + bl), i=l (B.8) 256

where wi is the ith coefficient of the expansion and ti is the corresponding function point and node located within the interval [-1, 1]. To evaluate the coefficients and nodes that give accurate results for higher degree polynomial, Lagrange polynomials can be utilized [19] to extract the coefficients and the nodes. Table B.1 shows both nodes and coefficients for several values of n starting from 1 to 5. 257

n Node t(n,j) Coefficient w(n,j) 2 0.57735 1.0000 0.57735 1.0000 3 0.77460 0.55556 00000 0.88889 -0.77460 0.55559 4 0.86114 0.34785 0.33998 0.65214 0.33998 0.65214 - 0.86114 0.34785 5 0.90618 0.23692 0.53847 0.47863 00000 0.56889 - 0.53847 0.47863 - 0.90618 0.23692 Table B.1: Weights and node values required for numerical evaluation of integrations using Gaussian quadrature scheme. 258

BIBLIOGRAPHY 259

BIBLIOGRAPHY [1] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, National Bureau of Standards, 1972. [2] C. A. Balanis, Advanced Engineering Electromagnetics, John Wiley & Sons, 1989. [3] C. A. Balanis, Antenna Theory, John Wiley & Sons, 1997. [4] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz, "AIM: Adaptive Integral Method for Solving Large-Scale Electromagnetic Scattering and Radiation Problems," Radio Science, pp. 1225-1251, September 1996. [5] Y. Y. Botros, E. S. Ebbini, and J. L. Volakis, "Optimal Synthesis of Phased Array Field Patterns in the Presence of the Rib cage," in the IEEE International Ultrasonics Symposium, Toronto, Canada, October 1997. [6] Y. Y. Botros, E. S. Ebbini, and J. L. Volakis, "Phased Array Pattern Synthesis for Hyperthermia Treatment of Liver Cancers," in the Joint Annual Meeting between the Radiation Research Society and the North America Hyperthermia Society, Rohde Island, Providence, May 1997. [7] Y. Y. Botros, E. S. Ebbini, and J. L. Volakis, "Two Step Hybrid Virtual ArrayRay (VAR) Technique for Focusing through the Rib Cage," IEEE Transactions on Ultrasonics, Ferroelectronics and Frequency Control, July 1998. [8] Y. Y. Botros, J. Gong, S. R. Legault, J. L. Volakis, and T. B. A. Senior, "Implementation Aspects of the Perfectly Matched Absorbing Layers for Finite Element Method Applications," in the Fourteenth National Radio Science Conference (invited paper), Cairo, Egypt, March 1997. [9] Y. Y. Botros, J. Gong, and J. L. Volakis, "PML Study for Finite Element Modeling of Antenna and Microwave Circuits," in the Thiteenth Applied Computational Electromagnetics (ACES) Conference, Naval Postgraduate School (invited paper), Monterey, March 1997. [10] Y. Y. Botros and J. L. Volakis, "Approximate Inverse Preconditioning of GMRES Solver Applied to PML Applications," in the International Workshop on Finite Elements for Microwave Engineering (invited paper), France, July 1998. 260

[11] Y. Y. Botros and J. L. Volakis, "Convergence Improvements of Iterative Solvers for Poorly Conditioned Antenna Applications," in the AP-S International Symposium and URSI Radio Science Meeting, June 1998. [12] Y. Y. Botros and J. L. Volakis, "Detailed Convergence Study for Perfectly Matched Layer (PML) Applications," in the Fourteenth Applied Computational Electromagnetics (ACES) Conference (invited paper), Naval Postgraduate School, Monterey, March 1998. [13] Y. Y. Botros, J. L. Volakis, and E. S. Ebbini, "Analysis and Synthesis of MultipleFocus Phased Array Heating Patterns Through the Rib Cage: A simulation Model," in the Thirteenth National Radio Science Conference (invited paper), Cairo, Egypt, volume 1, March 1996. [14] Y. Y. Botros, J. L. Volakis, E. S. Ebbini, and P. D. VanBaren, "Deep Localized Hyperthermia Through the Rib cage Using Ultrasound Heating Patterns," in the Third Joint Meeting Between the Acoustical Society of America and the Acoustical Society of Japan, Hawaii, December 1996. [15] Y. Y. Botros, J. L. Volakis, E. S. Ebbini, and P. D. VanBaren, "A Hybrid Computational Model for Ultrasound Phased Array Heating in Presence of Strongly Scattering Obstacles," IEEE Transactions on Biomedical Engineering, November 1997. [16] Y. Y. Botros and J. Volakis, "On the Convergence of the FEM Systems Terminated by the Perfectly Matched Layer (PML) Absorbers," IEEE Microwave and Guided Wave Letters (Submitted). [17] J. J. Bowman, T. B. A. Senior, and P. L. E. Uslenghi, Electromagnetic and Acoustic Scattering by Simple Shapes, Hemisphere Publishing Corp., 1987. [18] W. Brogan, Modern Control Theory, New Jersey: Prentice-Hall, 1985. [19] R. L. Burden and J. D. Faires, Numerical Analysis, PWS Publishers, 1985. [20] J. Chapelon, P. Faure, D. Cathignol, R. Souchon, F. Gorry, and A. Gelet, "The Feasibility of Tissue Ablation Using High Intensity Electronically Focused Ultrasound," in the IEEE International Ultrasonics Symposium, volume 2, pp. 1211-1214, 1993. [21] J. Chapelon, J. Margonari, D. Cathinol, A. Gelet, C. Guers, and Y. Theillere, "Thresholds for Tissue Ablation by Focused Ultrasound," IEEE Transactions on Ultrasonics, Ferroelectronics and Frequency Control, vol. 2, pp. 1653-1656, 1990. [22] W. C. Chew, Waves and Fields in Inhomogeneous Media, Van Nostrand Reinhold, 1992. 261

[23] C. J. Diederich and K. Hynynen, "The Feasibilty of Using Electrically Focused Ultrasound Arrays to Induce Deep Hyperthermia via Body Cavitation," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 38, no. 3, pp. 207-219, 1991. [24] C. Dorme and M. A. Fink, "Ultrasonic Beam Steering Through Inhomogeneous Layers with a Time Reversal Mirror," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 43, no. 1, pp. 167-175, January 1996. [25] J. Driller and F. L. Lizzi, "Therapeutic Applications of Ultrasound: A Review," IEEE Engineering in Medicine and Biology Magazine, vol. 6, no. 4, pp. 33-40, 1987. [26] E. S. Ebbini, Deep Localized Hyperthermia with Ultrasound Phased Arrays Using the Pseudoinverse Pattern Synthesis Method, PhD thesis, University of Illinois at Urbana-Champaign, 1990. [27] E. S. Ebbini and C. Cain, "Multiple-Focus Ultrasound Phased-Array Pattern Synthesis: Optimal Driving Signal Distribution for Hyperthermia," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 36, no. 5, pp. 540-548, 1989. [28] E. S. Ebbini and C. Cain, "A Spherical-Section Ultrasound Phased Array Applicator for Deep Localized Hyperthermia," IEEE Transactions on Biomedical Engineering, vol. 38, pp. 634-643, 1991. [29] E. S. Ebbini, F. Ngo, and C. A. Cain, "The Pseudoinverse Pattern Synthesis Method: Experimental Verification Using a Prototype Cylindrical-Section Ultrasound Hyperthermia Phased-Array Applicator," in the IEEE International Ultrasonics Symposium, November 1989. [30] E. S. Ebbini, S. Umemura, M. S. Ibbini, and C. Cain, "A Cylindrical-Section Ultrasound Phased-Array Applicator for Hyperthermia Cancer Therapy," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 35, no. 5, pp. 561-572, September 1988. [31] E. S. Ebbini, H. Wang, M. O'Donnell, and C. A. Cain, "Acoustic Feedback for Hyperthermia Phased-Array Applicators: Aberration Correction, Motion Compensation and Multiple Focusing in the Presence of Tissue Inhomogeneities," in the IEEE International Ultrasonics Symposium, pp. 1343-1346, November 1991. [32] R. S. Elliott, Antenna Theory and Design, Prentice Hall, Inc., 1981. [33] R. B. et. al., Templetes for the Solution of Linear Systems: Building Blocks for Iterative Solvers, SIAM, 1994. [34] L. A. Frizzell, "Thresholds Dosages for Damage to Mammalian Liver by High Intensity Focused Ultrasound," IEEE Transactions on Ultrasonics, Ferroelectronics and Frequency Control, vol. 35, no. 5,, 1989. 262

[35] W. J. Fry, J. W. Barnard, F. J. Fry, R. F. Krumins, and J. F. Brennan, "Ultrasonic Lesions in the Mammalian Central Nervous System," Science, vol. 122, pp. 517-518, 1955. [36] A. Gelet, J. Y. Chapelon, J. Margonari, Y. Theilliere, F. Corry, R. Souchon, and R. Bouvier, "High Intensity Focused Ultrasound Experimentation on Human Benign Prostatic Hypertrophy," in Eur. Urology, volume 23, pp. 44-47, 1993. [37] D. Givoli, Numerical Methods for Problems in Infinite Domains, Elsevier, New York, 1992. [38] G. H. Golub and C. F. VanLoan, Matrix Computations, The Johns Hopkins Uinversity Press, 1985. [39] 0. S. Haddadin, Ultrasound Inverse Scattering for Tomographic Imaging and Self-Focusing Arrays, PhD thesis, The University of Michigan, Ann Arbor, 1997. [40] 0. S. Haddadin and E. S. Ebbini, "Self-Focusing Arrays for Imaging and Therapy Through Inhomogeneous Media," in the IEEE International Ultrasonics Symposium, November 1996. [41] D. E. Hall, Basic Acoustics, Krieger, Melbourne, Florida, 1993. [42] R. F. Harington, Time Harmonic Electromagnetic Fields, McGraw Hill, 1994. [43] R. F. Harrington, Field Computation by Moment Methods, IEEE Press with Oxford University Press, 1993. [44] F. Hetzel, Hyperthermia in vivo, in Physical Aspects of Hyperthermia, G. Nausbaum Ed., New York: AAPM, G. Nausbaum Ed., 1982. [45] K. Hynenen, A. Chung, T. Fjield, M. Buchanan, D. Daum, V. Colucci, P. Lopath, and F. Jolesz, "Feasibility of Using Ultrasound Phased Arrays for MRI Monitored Noninvasive Surgery," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 43, no. 6, pp. 1043-1053, November 1996. [46] K. Hynynen, "Hot Spots Created at Skin-Air Interfaces During Ultrasound Hyperthermia," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 36, no. 2, pp. 241-2248, 1989. [47] K. Hynynen and K. L. Davi, "Small Cylindrical Ultrasound Sources for Induction of Hyperthermia via Body Cavities," Int. J. Hyperthermia, vol. 9, no. 2, pp. 263 — 274, 1993. [48] R. Jian, Bioheat Transfer: Mathematical Models of Thermal Systems, ch2 in Hyperthermia and Cancer Therapy, F. Storm Ed., Boston, MA: G. K. Hall Medical Publishers, 1983. [49] L. E. Kinsler, A. R. Frey, A. B. Coppens, and J. V. Sanders, Fundamentals of Acoustics, John Wiley & Sons, 1982. 263

[50] J.-U. Kluiwstra, Y. Zhang, P. D. VanBaren, S. Strickbereger, E. Ebbini, and C. Cain, "Ultrasound Phased Array for Noninvasive Myocardial Ablation: Initial Studies," in the IEEE International Ultrasonics Symposium, volume 2, pp. 1605 -1608, 1995. [51] F. Kremkau, "Cancer Therapy with Ultrasound: A Historical Review," J. Clin. Ultrasound, vol. 7, pp. 287-300, 1979. [52] P. Lele, Local Hyperthermia by Ultrasound, in Hperthermia in Cancer Therapy, G. Nausbaum Ed., New York: AAPM, 1982. [53] P. Lele, Physical Aspects and Clinical Studies with Ultrasonic Hyperthermia, in Hperthermia in Cancer Therapy, F. Storm Ed., Boston, MA: G. K. Hall Medical Publishers, 1983. [54] P. Lele, Effect of ultrasound on solid mammalian tissues and tumors in vivo, in Ultrasound, M. Rapacholi, M. Grandolfo, and A. Rindi Ed., New York: Plenum Publishing Corp., 1987. [55] S. Leon, Linear Algebra with Applications, New York: Aacmillan, 1980. [56] D. Luenberger, Optimization by Vector Space Methods, John Wiely & Sons, 1976. [57] R. Magin and A. Peterson, "Noninvasive Microwave Phased Arrays for Local Hyperthermia: A Review," Int. J. of Hyperthermia, vol. 5, no. 4, pp. 429-450, 1989. [58] D. A. McNamara, C. W. I. Pistorius, and J. A. G. Malherbe, The Uniform Geometrical Theory of Diffraction, Artech House, Boston, MA, 1990. [59] M. A. Morgan, Finite Element and Finite Difference Methods in Electromagnetic Scattering, Elsevier, New York, 1990. [60] K. B. Ocheltree and L. Frizzel, "Sound Field Calculation for Rectangular Sources," International Journal of Hyperthermia, vol. 6, no. 6, pp. 267-279, 1990. [61] D. Pierce, Acoustics: An Introduction to its Physical Principles and Applications, Acoustical Society of America, Woodburg, New York, 1989. [62] L. Potter and A. Arun, "Energy Concentration in Bandlimited Extrapolation," IEEE Trans. Acoust. Speech, Signal Processing, vol. 37, no. 7, pp. 1027-1041, 1989. [63] R. Roemer, W. Swindell, S. Clegg, and R. Kress, "Simulation of a Focused Scanned Ultrasonic Heating of Deep Seated Tumors: The Effects of Blood Perfusion," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 31, no. 5, pp. 457-466, 1984. 264

[64] G. T. Ruck, D. E. Barrick, W. D. Stuart, and C. K. Krichbaum, Radar Cross Section Handbook, Plenum Press, 1970. [65] Y. Saad, Iterative Methods for Sparse Linear System, PWS Publishing Co., 1996. [66] N. Sanghvi, F. Fry, R. Foster, and G. C. R. Chua, "System Design Considerations for High Intensity Focused Ultrasound Device for the Treatment of Tissue Invivo," Med. Biol. Eng. Comp., vol. 29, pp. 748, 1991. [67] S. Saparato and W. Dewey, "Thermal Dose Determination in Cancer Therapy," Int. J. Rad. Onc. Biol. Phys., vol. 10, pp. 787-800, 1984. [68] R. Seip and E. S. Ebbini, "Invasive and Non-invasive Feedback for Ultrasound Phased Array Thermometry," in the IEEE International Ultrasonics Symposium, November 1995. [69] R. Seip, E. S. Ebbini, M. O'Donnell, and C. A. Cain, "Non-Invasive Detection of Thermal Effects due to Highly Focused Ultrasonic Fields," the IEEE International Ultrasonics Symposium, vol. 2, pp. 1229-1232, 1993. [70] R. Seip, P. VanBaren,, and E. S. Ebbini, "Dynamic Focusing in Ultrasound Hyperthermia Treatments Using Implantable Hydrophone Arrays," IEEE Transactions on Ultrasonics, Ferroelectronics and Frequency Control, vol. 41, no. 5, pp. 706-713, Sept. 1994. [71] T. B. A. Senior and J. L. Volakis, Approximate Boundary Conditions in Electromagnetics, IEE press, London, 1995. [72] J. Shen and E. S. Ebbini, "A New Coded Excitation Ultrasound Imaging System: Part II: Operator Design," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 43, no. 1, pp. 141-148, January 1996. [73] P. P. Silvester and R. L. Ferrari, Finite Elements for Electrical Engineers, 2nd ed., Cambridge Univ. Press Cambridge, U.K., 1990. [74] C. Simon,, P. D. VanBaren, and E. S. Ebbini, "Qunatitative Analysis and Applications of Non-Invasive Temperature Estimation using Diagnostic Ultrasound," in the IEEE International Ultrasonics Symposium, Toronto, Canada, October 1997. [75] D. H. Staelin, A. W. Morgenthaler, and J. A. Kong, Electromagnetic Waves, Prentice Hall, Inc., 1994. [76] W. L. Stutzman and G. A. Thiele, Antenna Theory and Design, John Wiley & Sons, 1981. [77] W. Swindell, "A Theoretical Study of Nonlinear Effects with Focused Ultrasound in Tissues: An Acoustic Bragg Peak," Ultrasound in Med. and Biology, vol. 11, no. 1, pp. 121-130, 1985. 265

[78] A. Tikhonov and V. Arsenin, Solutions of Ill-Posed Problems, Washington DC: Winston, 1977. [79] P. VanBaren and E. S. Ebbini, "Multi-point Temperature Control During Hyperthermia Treatments: Theory and Simulation," IEEE Transactions on Biomedical Engineering, vol. 42, no. 8, pp. 818-827, August 1995. [80] P. VanBaren, R. Seip, and E. S. Ebbini, "Real-Time Dynamic Focusing Through Tissue Inhomogeneities During Hyperthermia Treatments with Phased Arrays," in the IEEE International Ultrasonics Symposium, pp. 1815-1819, November 1994. [81] V. V. Varadan, A. Lakhtakia, and V. K. Varadan, editors, Field Representations and Introduction to Scattering, Elsevier Science Publishers, 1991. [82] J. Volakis, T. Ozdemir, and J. Gong, "Hybrid Finite Element Methodologies for Antennas and Scattering," IEEE Transactions on Antennas and Propagation, pp. 493-507, March 1997. [83] H. Wang, E. S. Ebbini, and C. A. Cain, "Computationally Efficient Algorithms for Control of Ultrasound Phased-Array Hyperthermia Applicators Based on a Pseudo Inverse Methos," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 37, pp. 274 —277, 1990. [84] H. Wang, E. S. Ebbini, M. O'Donnell, and C. A. Cain, "Phase Aberration Correction and Motion Compensation for Ultrasonic Hyperthermia Phased Arrays: Experimental Results," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 41, no. 1, pp. 34-43, January 1994. [85] H. Wang, E. S. Ebbini, M. O'Donnell, R. Seip, and C. A. Cain, "Adaptive 2-D Cylindrical Section Phased Array System for Ultrasonic Hyperthermia," in the IEEE International Ultrasonics Symposium, volume 2, pp. 1261-1264, 1992. [86] L. Ziomek, Acoustic Field Theory and Space-Time Signal Processing, CRC Press, 1995. 266