UNIVERSITY OF CALGARY Reducing the Dominant Skin Reflection in Radar Based Microwave Breast Imaging by Batoul Maklad A THESIS SUBMITTED TO THE FACULTY OF GRADUATE STUDIES IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING CALGARY, ALBERTA JANUARY, 2010 c Batoul Maklad 2010 Library and Archives Canada Bibliothèque et Archives Canada Published Heritage Branch Direction du Patrimoine de l’édition 395 Wellington Street Ottawa ON K1A 0N4 Canada 395, rue Wellington Ottawa ON K1A 0N4 Canada Your file Votre référence ISBN: 978-0-494-62091-5 Our file Notre référence ISBN: 978-0-494-62091-5 NOTICE: AVIS: The author has granted a nonexclusive license allowing Library and Archives Canada to reproduce, publish, archive, preserve, conserve, communicate to the public by telecommunication or on the Internet, loan, distribute and sell theses worldwide, for commercial or noncommercial purposes, in microform, paper, electronic and/or any other formats. . The author retains copyright ownership and moral rights in this thesis. Neither the thesis nor substantial extracts from it may be printed or otherwise reproduced without the author’s permission. L’auteur a accordé une licence non exclusive permettant à la Bibliothèque et Archives Canada de reproduire, publier, archiver, sauvegarder, conserver, transmettre au public par télécommunication ou par l’Internet, prêter, distribuer et vendre des thèses partout dans le monde, à des fins commerciales ou autres, sur support microforme, papier, électronique et/ou autres formats. L’auteur conserve la propriété du droit d’auteur et des droits moraux qui protège cette thèse. Ni la thèse ni des extraits substantiels de celle-ci ne doivent être imprimés ou autrement reproduits sans son autorisation. In compliance with the Canadian Privacy Act some supporting forms may have been removed from this thesis. Conformément à la loi canadienne sur la protection de la vie privée, quelques formulaires secondaires ont été enlevés de cette thèse. While these forms may be included in the document page count, their removal does not represent any loss of content from the thesis. Bien que ces formulaires aient inclus dans la pagination, il n’y aura aucun contenu manquant. Abstract Radar based microwave breast imaging has been proposed as a complementary breast imaging technique. One challenge with this approach is to subdue the overwhelming reflections that originate from the skin. These reflections are up to 90 dB larger than reflections from any tumours that may be present, and must be reduced in order to examine the interior of the breast. This thesis describes the development of a technique to reduce skin reflections that is applicable to data recorded from realistic, three-dimensional, numerical breast models. Metrics are outlined that assess the algorithm’s performance and the skin reduction technique is subsequently tested on a variety of scenarios. The results demonstrate that the method is consistent and effective, as the skin reflections are reduced on average by 95 dB in data sets simulated using different antennas and breast models with various shapes and glandular tissue distributions. ii Acknowledgements First, I would like to thank my family, for their continual support. I am blessed to have an amazing mother that is the source of my strength and guidance, as well as a brilliant father that made his expertise available at any time. I am also genuinely grateful to Dr Fear, for giving me the opportunity to work on such a great project and for consistently providing me with her understanding, mentorship and support. This thesis would not have materialized had it not been for a handful of key friends. I am extremely thankful for the unfailing hospitality and encouragement provided by the McNeil family. You guys got me through some tough times and I will never forget it. I would also like to express my profound appreciation towards Sarah Manske and Nadia Barnieh for being there when I needed them most. I am deeply grateful for the efforts of Billy Wu, Philip Chu and Kelly Goss towards keeping me from going INSANE, along with Jeff Sill and Jeremie Bourqui for always being friendly faces to confide in. Finally, I thank god, for all of my blessings and for sending me all these wonderful people to help me. iii Table of Contents Abstract ii Acknowledgements iii Table of Contents iv List of Tables vi List of Figures vii List of Abbreviations ix 1 Introduction 1.1 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Background 2.1 Radar-Based Microwave Breast Cancer Detection . . . . . 2.1.1 Realistic Numerical Breast Models . . . . . . . . . 2.1.2 Data Collection Schemes . . . . . . . . . . . . . . . 2.1.3 Comparison of Microwave Breast Imaging Systems 2.2 Skin Reflection Removal . . . . . . . . . . . . . . . . . . . 2.2.1 RLS/Woody . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Modified RLS . . . . . . . . . . . . . . . . . . . . . 2.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Neighbourhood Modified RLS 3.1 Application of Modified RLS to 3-D Breast Model 3.1.1 3-D Breast Model . . . . . . . . . . . . . . 3.1.2 Processing Signals From 3-D Breast Model 3.1.3 Results & Discussion . . . . . . . . . . . . 3.2 Performance Measures . . . . . . . . . . . . . . . 3.2.1 Single Signal Analysis . . . . . . . . . . . 3.2.2 Image Evaluation . . . . . . . . . . . . . . 3.3 Neighbourhood Selection . . . . . . . . . . . . . . 3.3.1 Row Based Neighbourhoods . . . . . . . . 3.3.2 Row & Column Based Neighbourhoods . . 3.3.3 Beamwidth Based Neighbourhoods . . . . iv . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 3 4 . . . . . . . . 6 6 8 14 18 20 23 29 36 . . . . . . . . . . . 37 37 38 39 42 45 45 51 54 55 65 71 3.4 3.5 3.6 3.7 3.3.4 Antenna Distributions & Cross Correlation Neighbourhoods Filter Weight Estimation . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Adaptive Skin-Dominant Window Selection . . . . . . . . . 3.4.2 Adaptive Rank Selection . . . . . . . . . . . . . . . . . . . . Feedback Filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Mathematical Derivation of Feedback . . . . . . . . . . . . . 3.5.2 Feedback Results . . . . . . . . . . . . . . . . . . . . . . . . Neighbourhood Modified RLS Overview . . . . . . . . . . . . . . . Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Performance of Neighbourhood Modified RLS 4.1 Response to Varying Glandular Distributions . . 4.1.1 Models . . . . . . . . . . . . . . . . . . . 4.1.2 Results . . . . . . . . . . . . . . . . . . . 4.2 Application to Different Breast Shapes . . . . . 4.2.1 Model . . . . . . . . . . . . . . . . . . . 4.2.2 Results . . . . . . . . . . . . . . . . . . . 4.3 Application to Different Antennas . . . . . . . . 4.3.1 Model . . . . . . . . . . . . . . . . . . . 4.3.2 Results . . . . . . . . . . . . . . . . . . . 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . 4.4.1 Neighbourhood Comparison . . . . . . . 4.4.2 Clutter Quantification . . . . . . . . . . 4.4.3 Skin Response Suppression . . . . . . . . 4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 . 81 . 82 . 87 . 89 . 90 . 92 . 94 . 100 . . . . . . . . . . . . . . 102 102 102 105 109 110 111 115 116 118 123 123 125 129 131 5 Conclusions 132 Bibliography 135 A MIST:Mapping MR Intensities to Dielectric Properties 144 B Derivation of the Kalman Gain Vector For the RLS Algorithm 151 v List of Tables 2.1 Data Acquisition Summary . . . . . . . . . . . . . . . . . . . . . . . . 20 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 Single Signal Analysis Measures . . . . . . . . . . . . . . . . . . . . . Image evaluation tools . . . . . . . . . . . . . . . . . . . . . . . . . . Neighbourhood Size for Each Row Using Row Selection Scheme . . . Image Quality Using Row Based Neighbourhoods. . . . . . . . . . . . Neighbourhood Size for Each Row Using Row-Column Selection . . . Image Quality Using Row-Column Based Neighbourhoods. . . . . . . Image Quality Using Modified Row-Column Based Schemes. . . . . . Image Quality Using Beamwidth Based Neighbourhoods. . . . . . . . Average Neighbourhood Size Using Normalized Cross Correlation . . Image Quality Using Beamwidth & Cross Correlation Based Schemes. Statistics of Estimated Ranks . . . . . . . . . . . . . . . . . . . . . . 51 54 57 59 67 67 69 74 78 78 88 4.1 4.2 4.3 4.4 4.5 Image Quality of Breast Models with Varying Interior Properties. . . Resulting Image Quality of Two Different Simple Breast Models. . . . Dielectric Properties of Materials Used in Simulations . . . . . . . . . Differences Between the Wu-King Dipole and the BAVA . . . . . . . Resulting Image Quality of Breast Models Scanned Using the BAVA. 109 114 116 117 122 A.1 Corresponding Pixel Intensity and Dielectric Property Groups . . . . 147 vi List of Figures 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11 2.12 2.13 2.14 2.15 Pipeline of Signal Processing Procedures . . . . . . . . . . . . . Simulation Procedure . . . . . . . . . . . . . . . . . . . . . . . . Planes of the Human Body . . . . . . . . . . . . . . . . . . . . . Example of Pixel Intensity Histogram . . . . . . . . . . . . . . . Relative Permittivity Map of Breast Model From TSAR Group MIST 2-D Data Acquisition . . . . . . . . . . . . . . . . . . . . MIST in the Prone Orientation . . . . . . . . . . . . . . . . . . Data Collection Scheme Employed by TSAR System . . . . . . RASOR Data Acquisition . . . . . . . . . . . . . . . . . . . . . Skin Reflections From 3-D Realistic Breast Model . . . . . . . . Woody Averaging Skin Subtraction Flow Chart . . . . . . . . . Comparison of Woody Averaging and RLS/Woody . . . . . . . Modified RLS Input Vector Initialization Flow Chart . . . . . . Flow Chart of Original Modified RLS . . . . . . . . . . . . . . . Application of the Modified RLS to Simulated Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 9 11 12 13 16 16 17 18 21 27 28 32 34 35 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14 3.15 3.16 3.17 3.18 3.19 3.20 3.21 Dielectric Properties of Simple Breast Model . . . . . . . . . . . . Resulting Images When Modified RLS is Applied to a 3-D Model Flow Chart of Modified RLS . . . . . . . . . . . . . . . . . . . . . Calculating the Calibrated Data . . . . . . . . . . . . . . . . . . . Calculation of Energy Ratio . . . . . . . . . . . . . . . . . . . . . Comparison of Tumour Estimate and Known Tumour Response . Reference Image . . . . . . . . . . . . . . . . . . . . . . . . . . . . Focused Images of Row Neighbourhoods . . . . . . . . . . . . . . Tumour-to-Skin Ratios Using Row Neighbourhoods . . . . . . . . Mean Energy Ratio Using Row Neighbourhoods . . . . . . . . . . Mean Offset of Each Row Using Row Neighbourhoods . . . . . . . Examples of Late-Time Clutter . . . . . . . . . . . . . . . . . . . Mean MSE of Each Row Using Row Neighbourhoods . . . . . . . Preserved Tumour Responses Using Row Neighbourhoods . . . . . 5r2c Antenna Distribution . . . . . . . . . . . . . . . . . . . . . . Tumour-to-Skin Ratios Using Row-Column Neighbourhoods . . . Preserved Tumour Responses Using Row-Column Schemes . . . . Selection of Patch Neighbourhood . . . . . . . . . . . . . . . . . . Preserved Tumour Responses Using Patch Neighbourhoods . . . . Comparison of Recorded Signals . . . . . . . . . . . . . . . . . . . Comparison of Potential Neighbours . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 43 44 46 48 50 52 58 60 61 62 63 63 65 69 70 71 72 75 76 79 vii 3.22 3.23 3.24 3.25 3.26 3.27 3.28 3.29 Resulting Images Using the 2prl 90 Neigbourhood . . . . . . . . Flow Chart of Updated Modified RLS . . . . . . . . . . . . . . . Example of Two Signals That Require a Different SDW . . . . . Tumour-to-Signal Ratio Using Blind & Ideal SDW Selection . . Flow Chart of the Enhanced Modified RLS . . . . . . . . . . . . Flow Chart of Feedback Filter . . . . . . . . . . . . . . . . . . . Comparison of Tumour-to-Skin Ratio With & Without Feedback Flow Chart of the Neighbourhood Modified RLS . . . . . . . . . 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15 4.16 4.17 4.18 Permittivity Maps of Breast Models with Varying Tissue Distributions 104 Tumour-to-Skin Ratios of Models with Varying Tissue Distributions . 106 Tumour-to-Signal Ratio of Models with Varying Tissue Distributions 107 Scan Pattern Applied to Simple Breast Models of Different Shapes . . 110 Relative Permittivity Map Simple Breast Models of Different Shapes . 111 Tumour-to-Skin Ratios of Models with Different Shapes . . . . . . . . 112 Tumour-to-Signal Ratio Using Blind & Ideal SDW on MRI3R . . . . 113 Images Generated From Model With Different Geometry . . . . . . . 114 Scan Pattern Applied to Simple Breast Models Excited With BAVA . 117 Relative Permittivity Map Simple Breast Models Excited With BAVA 117 Comparison of Tumour Responses Recorded Using BAVA . . . . . . . 119 Resulting Tumour-to-Skin Ratios with BAVA . . . . . . . . . . . . . 120 Resulting Tumour-to-Signal Ratios with BAVA . . . . . . . . . . . . . 121 Histograms of Resulting Neighbourhood Sizes . . . . . . . . . . . . . 125 Comparison of Signals Acquired Using Different Antennas . . . . . . 126 Selection of Window for Clutter Ratio Measure . . . . . . . . . . . . 127 Clutter Ratios Using Different Antennas . . . . . . . . . . . . . . . . 128 Skin Subtraction Ratios Using Different Antennas . . . . . . . . . . . 130 A.1 A.2 A.3 A.4 A.5 MR Image Used to Create a Numerical Breast Model . . . . . . . . . Dielectric Property Histograms Using Different Mapping Approaches Preprocessed MR Template Used for Dielectric Mapping . . . . . . . Gaussian Mixture Modeling of MR Intensity Histogram . . . . . . . . Linear Mapping of Pixel Intensities to Dielectric Properties . . . . . . viii . . . . . . . . . . . . . . . . . 80 . 81 . 83 . 86 . 89 . 92 . 93 . 101 144 148 149 149 150 List of Abbreviations RLS UWB MIST RASOR TSAR MR FDTD 2-D 3-D MSE R ~p ~q or w ~ CData NoTumData Tactual NoSkData ResSkData p2p SCR prl 2prl 2prl 90 SDW GCV BAVA Recursive Least Squares Ultra-Wideband Microwave Imaging via Space-Time Beamforming Real Aperture, Synthetically Organised Radar Tissue Sensing Adaptive Radar Magnetic Resonance Finite Difference Time Domain Two Dimensional Three Dimensional Mean Squared Error Autocorrelation Matrix Cross Correlation Vector Filter Weight Coefficient Vector Calibrated Data No Tumour Data Known Tumour Response Skin-Subtracted Signal Residual Skin Response Peak-to-Peak Signal-to-Clutter Ratio Patch Row Limit Neighbourhood Double Width Limited Row Patch Double Width Limited Row Patch With 90% Normalized Cross Correlation Skin Dominant Window Generalized Cross-Validation Balanced Antipodal Vivaldi Antenna ix 1 Chapter 1 Introduction Breast cancer is the most prevalent form of cancer in women worldwide, and accounts for approximately 22% of all new female cancer diagnoses. In Canada over 22,000 women are diagnosed and more than 5,000 breast cancer related deaths are reported annually[1]. Currently, mammography is considered to be the gold standard imaging technique; however, it exhibits a few limitations in terms of its sensitivity and specificity [2]. This modality also exposes the patients to potentially harmful x-rays, uncomfortably compresses the breast, has difficulty with imaging dense breasts and only generates two dimensional (2-D) images [3]. These factors have motivated the development of new breast cancer imaging techniques that are relatively inexpensive, avoid ionizing radiation and are capable of generating three dimensional (3-D) images [4]. Microwave breast cancer detection relies on the dielectric property differences between the various tissues in the breast. Past studies suggest that there is a significant difference between the dielectric properties of healthy and diseased breast tissues [5, 6, 7, 8]. A recent large scale study of the dielectric properties of healthy, malignant and benign tissue samples, that spanned a frequency range of 0.5-20 GHz, was conducted collaboratively by the University of Calgary and the University of Wisconsin [9]. This investigation revealed that the properties are strongly dependent on the ratio between the adipose tissue and the glandular or fibroconnective tissue present in the sample [9]. The contrast between malignant and adipose domi- 2 nated tissues can be as large as a ratio of 10:1, while the difference between cancerous and glandular or fibroconnective tissue is 10% or less[10]. In spite of the low property contrast between malignant and fibro-glandular tissues, microwave breast imaging is being actively pursued by an increasing number of research groups. The main advantages of microwave imaging are that it eliminates patient exposure to x-rays and the need to compress the breast. There are three general approaches to microwave imaging; passive, hybrid and active. An example of a passive method is microwave radiometry wherein the temperature of the breast is measured, and tumours are detected on the basis that they are significantly warmer than the rest of the breast [3]. Hybrid imaging selectively heats potential tumours with microwaves and records the resulting compression waves using an ultrasound transducer [11]. The tumour is expected to generate the largest compression wave, as its relatively large conductivity would cause it to expand the most when heated. Active microwave imaging is broken up into two general categories, tomography and radar based approaches. Microwave tomography irradiates the breast using an antenna array and measures the transmitted fields which are then used to generate a dielectric property map of the breast. Radar based methods, such as confocal microwave imaging, record reflected fields from the breast, and focus these signals to localize significant scatterers [12]. Microwave methods are not in wide-spread clinical use, as several challenges need to be addressed. One of the challenges microwave radiometry faces is that it requires high precision low noise receivers because the signals generated are on the order of 10−13 to 10−16 watts [13]. Microwave induced acoustic imaging has the predicament of ensuring the breast is uniformly heated, while microwave tomography is hindered by 3 its computationally intensive nature. Radar based microwave imaging does not face the above mentioned difficulties and also has the potential to generate 3-D images of the breast. However, reliable detection of tumours in a realistic environment requires the detection of small differences in the reflected signals. This obstacle is made more challenging by the fact that the largest reflection is caused by the skin. The skin has been reported to cause reflections that are up to 40-60 dB larger than the tumour reflection in 2-D and simple 3-D numerical breast models [14, 15]. Consequently, this thesis will focus on reducing the reflections from the skin in order to permit the examination of the breast interior with radar based microwave imaging. 1.1 Objectives The goal of this research is to suppress the reflections from the skin found in signals collected from simulations of realistic, 3-D breast models. Reduction of the skin reflection is expected to make tumour detection feasible. To date, methods in literature focus on 2-D or simple 3-D models; consequently, there is a need for an algorithm that can be applied to realistic 3-D models. More specifically the objectives are to: 1. Design a skin reflection suppression algorithm that preserves the tumour reflection and is effective when applied to realistic 3-D breast models. 2. Develop metrics to evaluate the success of skin reflection reduction and tumour response preservation. 3. Evaluate the algorithm’s performance on breast models with different shapes and varying glandular tissue distributions. 4 4. Compare the algorithm’s effectiveness when different sensors are used to collect the data. 1.2 Thesis Overview Chapter 2 contains background information and a literature review of the pertinent findings previously published in regards to microwave breast cancer detection. The development of realistic 3-D numerical breast models is detailed, followed by a description of the currently reported data acquisition systems. Next, the major obstacles involved in removing the skin reflections from the collected signals are presented along with the algorithms currently employed to reduce the skin reflection. Chapter 3 describes the development of a skin reflection reduction algorithm designed for data acquired from 3-D breast models. First, the capability of a previously published technique is appraised. Next, a set of performance measures are outlined to assess whether the skin reflections have been suppressed and if the tumour reflections have been preserved. The previously published approach is then combined with several different data grouping methods and tested using the aforementioned metrics. Finally, the previously published method is further refined in conjunction with the proposed data selection algorithm. Chapter 4 presents the performance of the developed algorithm when applied to a collection of different realistic 3-D numerical breast models. Models of the same shape with a variety of glandular distributions are examined first and the algorithm’s ability to suppress the skin reflection is investigated. Next, the impact of varying the shape of the breast model, while maintaining a very simple interior breast property 5 distribution is explored. The robustness of the algorithm is also tested by applying it to data acquired using a different sensor. Finally, a comparison of the obtained results is conducted to establish the proposed method’s proficiency. Chapter 5 summarizes the contributions of this thesis and outlines alternative potential applications of the algorithm. To the best of our knowledge this is the only algorithm currently reported in literature that is capable of reducing skin reflections recorded from realistic 3-D breast models. The contributions of the work described in this thesis include: 1. providing an understanding of the data recorded from simulation of 3-D realistic numerical breast models, 2. developing and testing a similarity criteria for data that is helpful in several signal processing challenges, is applicable to signals recorded from monostatic and multistatic systems, and is also adaptive to the sensor, and 3. enhancing a pre-existing skin subtraction algorithm so that it is almost entirely adaptive to the data, and is effective on signals recorded from 3-D MR-based numerical breast models. 6 Chapter 2 Background Microwave breast cancer detection has been proposed as a complementary technique to those currently employed for breast imaging. Radar-based imaging is one example of a microwave based approach and is the system explored in this thesis. In this chapter, the basic ideas behind radar-based methods are discussed. The numerical breast models used to validate this technique are summarized, and an overview of the data collection schemes under investigation is provided. A key challenge in radarbased breast imaging is the reduction of the reflection caused by the skin in order to permit the examination of the breast interior. Consequently, the chapter concludes with a comprehensive description of the two prominent methods to remove the skin reflection: the Modified RLS[16] and the Woody/RLS[17]. 2.1 Radar-Based Microwave Breast Cancer Detection Radar-based microwave breast cancer detection involves illuminating a breast at microwave frequencies using an ultrawide band (UWB) pulse and recording the ensuing reflections [12]. These reflections occur as a result of the pulse impinging upon an interface of two tissues with distinct dielectric properties. The goal of radar based methods is to detect scatterers within the imaging volume. This is accomplished by collecting reflection profiles at various sites around the breast, and sending the data through a pipeline of signal processing procedures. 7 The aim of the signal processing steps is to reduce reflections from objects other than tumours. First, the data is calibrated by removing the incident signal and any characteristic reflections of the antenna from the recorded signal. Next, the overwhelming reflections from the skin are reduced in order to reveal the previously masked tumour response [17, 16]. This process is referred to as ”skin subtraction”. Once the skin subtraction is complete, a clutter reduction algorithm may be applied to isolate tumour reflections from responses originating from glandular tissue and other internal structures of the breast[18, 19]. Finally, images are created by focusing the processed signals, as the remaining coherent responses are presumed to originate from tumours[12, 16, 20, 21, 22]. Fig. 2.1 summarizes the steps in the signal processing pipeline. Figure 2.1: Pipeline of Signal Processing Procedures: Four basic steps are required to create a radar based image of a breast tumour. Currently, three radar-based systems have reported experimental and simulated results in literature. These include: 1. Microwave Imaging via Space Time beamforming (MIST) from the University of Wisconsin; 2. Real Aperature Synthetically Organized Radar (RASOR) from the University of Bristol; and 3. Tissue Sensing Adaptive Radar (TSAR) from the University of Calgary, which is the system examined in this thesis. 8 In order to validate these proposed systems, a variety of data collection schemes and breast models have been developed for both simulations and experiments. This thesis focuses on reducing the skin reflections in signals acquired from 3-D breast models in a simulated environment. Consequently, the following sections describe the breast modeling approaches and data acquisition techniques applied, prior to reviewing the previously reported skin subtraction algorithms. 2.1.1 Realistic Numerical Breast Models The propagation of microwave energy through a medium is primarily dependent on its dielectric properties. As a result, numerical models are created by assigning appropriate dielectric properties to spatial distributions representing the various tissue structures found in the breast. Spatial distributions are typically established with magnetic resonance (MR) images, which provide a map of the glandular and fatty tissues in the breast. The resulting breast models are placed in an electromagnetics simulator, and the interactions of the model with incoming microwaves are monitored. The finite-difference time-domain (FDTD) method has been used extensively in the study of radar-based microwave imaging [23, 24, 25, 26]. Fig. 2.2 shows the steps required to simulate data using an MR based numerical breast model. The development of realistic numerical breast models from MR images has been investigated primarily by two institutes: the University of Wisconsin, and the University of Calgary. As the results of the skin subtraction algorithm employed by the MIST system have solely been reported for 2-D cases, only their 2-D models are described. Conversely, the development of 3-D models is outlined for the TSAR system, as these are the models that are used for testing throughout this thesis. 9 Figure 2.2: Simulation Procedure: Data is simulated by illuminating a numerical breast model at several locations with an antenna, and recording the ensuing reflections using the FDTD method. The relative permittivity maps seen on the right hand side of the figure, ranges from 3 (blue) to 50 (red). University of Wisconsin The numerical breast models developed at the University of Wisconsin are based on high resolution MR images acquired with the patient in the standard prone position (i.e. laying face down). For 2-D models, a slice of the high resolution MR image is used as the model template [27]. Subtle inherent MR image artifacts are filtered out and the resolution of these preprocessed images is enhanced from 0.625 x 0.625 mm pixels to 0.5 x 0.5 mm pixels using linear interpolation [28, 29, 30]. The pixels associated with the skin are segmented and replaced with a 2 mm thick artificial skin layer [28, 27]. Finally, the frequency dependent dielectric properties of the pixels representing the interior of the breast are assigned based on the MR pixel intensities[28]. 10 Several methods have been proposed to map MR pixel intensities to dielectric properties. The mapping techniques assume a linear relationship between the MR pixel intensities and the dielectric properties. The main differences between the approaches are the number of ranges used in the linear mapping, and the manner in which the ranges are selected. The proposed methods include: • Uniform: MR pixel intensity histogram linearly mapped to a single range of dielectric properties [29]. • Bimodal : MR pixel intensity histogram is visually segmented and linearly mapped to 2 ranges of dielectric properties, representing fatty and glandular tissues [29]. • Piecewise Linear : MR pixel intensity histogram is visually segmented and linearly mapped to 3 ranges of dielectric properties, representing fatty and glandular tissues, as well as a transitional region [29]. • Piecewise Linear Using Gaussian Mixture Models: MR pixel intensity histogram is segmented using a Gaussian mixture model and linearly mapped to 3 ranges of dielectric properties, representing fatty and glandular tissues, as well as a transitional region [30]. Further details on these techniques can be found in Appendix A. University of Calgary The TSAR group has focused primarily on creating 3-D MR-based breast models. A 1.5 Tesla Siemens Sonota MR scanner is used to acquire sagittal images using a breast coil and T1 weighted sequence [31, 32]. Fig. 2.3 demonstrates the three of the major 11 anatomical planes of the human body: the sagittal, coronal, and transverse planes. The sagittal plane images have a resolution of 0.43 x 0.43 mm and are collected at 1.2 mm intervals resulting in 512 x 512 pixels per plane[31]. Once the MR data is collected, the tissues inside the breast are segmented and assigned a dielectric value. Two approaches have been proposed for tissue segmentation and property mapping. Figure 2.3: Planes of the Human Body: Three of the major planes of the human body are the sagittal, coronal and transerves planes [33]. The first method involves segmenting the MR image to identify the interior of the breast from its contour using an edge detection algorithm [31]. This process is followed by classifying different tissue types within the breast by placing thresholds on the pixel intensity histogram of the MR image. Fig. 2.4 is an example of a histogram generated from an MR image. Using this histogram, four non-dispersive breast models are created [32]. 1. Simple breast where all the pixels inside the breast are mapped to a relative permittivity of 9 and a conductivity of 0.4 S/m. 2. Simple glandular where pixel intensities greater than a particular threshold 12 Figure 2.4: Example of Pixel Intensity Histogram: The MR pixel intensity histogram used to generate realistic numerical breast models. are set to a relative permittivity of 16 and a conductivity of 1 S/m. 3. Piecewise-linear where intensities are linearly mapped between a relative permittivity of 1 to 36 and a conductivity of 1 to 4 S/m. 4. Bimodal where intensities within a range representing glandular tissues are linearly mapped to a relative permittivity and conductivity that span 15.2 to 27.5 and 1.7 to 3.0 S/m respectively. Pixel values below this range are set to a relative permittivity of 9 and conductivity of 0.4 S/m. For all the models described above, an artificial skin layer is introduced that is 2 mm thick and has a relative permittivity of 36 and a conductivity of 4 S/m. A variation on the piecewise-linear model (piecewise linear II) exhibits varying skin properties and thickness achieved by mapping pixel intensities associated with the skin. The second method takes the frequency dependent properties of breast tissues into account and is applied to MR data collected from patients diagnosed with breast cancer [34]. First, MR images acquired before contrast agent injection are exam- 13 ined. A skin layer with uniform thickness and pixel intensity is overlaid on the breast surface. A k-means approach is then applied to the breast interior and 19 different tissues are classified using ISeg (Schmid & Partner Engineering AG,Zurich, Switzerland) [34]. The pixel intensities are mapped to frequency dependent dielectric properties with the use of the three tissue groups identified in a large scale study of breast tissue properties at microwave frequencies: 0-30% adipose, 31-84% adipose, and 85-100% adipose [9, 34]. A threshold of 25 to 75% about the nominal property value of each group defined in [9] sets the linear mapping ranges. Finally, a realistic tumour is added to the dielectric model by subtracting MR images with and without contrast agent to identify the tumour. The tumour properties are also mapped to dielectric values. Fig. 2.5 demonstrates the resulting numerical model at 3 GHz in the coronal plane. Figure 2.5: Relative Permittivity Map of Breast Model From TSAR Group: Cornal slice of an MR based dispersive numerical breast model at 3 GHz [34] Preliminary trials have been conducted to validate the accuracy of the numerical breast models. The trials compare the simulated signals to the reflections recorded using an experimental TSAR prototype from the same volunteer. The numerical breast models are created using an MR image of the volunteer as a template and 14 the simulated antenna placements mimic the relative distance between the antenna and the breast in the experimental prototype. Good agreement is achieved when comparing skin reflections recorded in simulated and experimental settings [35]. The realistic shape and tissue distributions exhibited by 3-D, MR-based models present a new challenge as the resulting signals are significantly different from those recorded from 2-D models. The first major difference is that the shape of the skin reflections recorded at different locations around the 3-D models fluctuates considerably as the curvature of the breast varies. The second challenge is that the duration of the skin response cannot be as easily identified. Therefore, defining an isolated window over which the algorithm operates is not effective as the duration of the skin response is now longer which causes it to overlap with the tumour response. Consequently, more sophisticated signal processing and data collection methods are required to compensate for these differences. 2.1.2 Data Collection Schemes A variety of different data collection schemes are utilized in radar based imaging. In general, an UWB pulse is transmitted towards the breast and the resulting reflections are recorded. The major differences between the proposed data acquisition systems lie in the number of antennas used for the detection process, the placement of these sensors and the frequency ranges of the emitted pulses. There are two basic data acquisition schemes in radar based imaging: bistatic and monostatic. A bistatic system requires a minimum of two antennas, where one element acts as a transmitter and the other acts as a receiver. Multistatic systems are an example of a bistatic data collection approach where more than two antennas are used. A 15 multistatic breast scan is conducted by transmitting a pulse from one antenna and recording the electromagnetic echoes at all or a subset of the other receivers. The pulse is then emitted by another antenna until all the elements in the array have been used as a transmitter [21, 22]. A monostatic system uses one antenna as the transmitter and the receiver. An example of a monostatic data collection scheme involves creating a synthetic array by scanning one antenna around the breast and recording the ensuing reflections [12]. Alternatively, an antenna array can also be used to monostatically scan the breast by exciting one element in the array and acquiring the resulting signal at that same element. This process is then repeated for all antennas in the array[28]. Next, a brief overview of the data collection approaches of the three main radar based systems that are under investigation is given. MIST Two general data acquisition schemes have been explored for the MIST system using FDTD. The first system has the patient in the supine position which naturally flattens the breast. The aim behind this orientation is to gain access to the upper outer quadrant of the breast where 50% of tumours are generally found [28]. The supine orientation is simulated using a 2-D MR-based breast model, as seen in Fig. 2.6. A monostatic scan is conducted using a 17 element antenna array that is placed directly on the surface of the skin. The array spans 8 cm and each antenna is modeled as a current source. The signals are acquired by recording the electric-field response at the antenna that is illuminating the breast [16]. The excitation used is a 110 ps full-width half-maximum differentiated Gaussian pulse that has significant energy 16 between 1-11GHz and a peak at 6 GHz [28, 16, 36, 37]. Figure 2.6: MIST 2-D Data Acquisition: Sagittal plane breast model in the supine position, where the current sources are represented by black dots, and the tumour is the white dot near the center [16] The second system places the patient in the prone position [27]. This scan orientation is simulated using a 2-D coronal breast model in conjunction with a circular 17 element antenna array that is 10 cm in diameter, as shown in Fig. 2.7 [38]. In this scenario, the elements are no longer in contact with the skin, but once again they are modeled using current sources. The data is still collected in a monostatic manner and a differentiated Gaussian with a full-width half-maximum frequency range of 1.67-11.1 GHz is used to excite the sources [38]. Figure 2.7: MIST in the Prone Orientation: Realistic 2-D coronal plane model that contains a spherical tumour illustrated by the black dot seen on the left hand side of the breast interior. The surrounding white circles represent the current sources that are used to scan the 2-D breast model [38]. 17 TSAR The TSAR data acquisition system scans the patient monostatically in the prone position. Using FDTD, simulations are conducted where a resistively loaded (Wu King) dipole antenna illuminates a breast model with a differentiated Gaussian pulse with a full-width half-maximum of 0.17 ns and a center frequency of 4 GHz [39, 40]. The recorded signals are proportional to the voltage associated with the back scattered electric field observed at the feed of the antenna used to illuminate the breast. The breast model is scanned with a synthetic array that consists of 9 or 12 rows with 10 or 32 receivers on every level. The number of rows and antenna elements on a given row is selected such that an adequate estimate of the breast surface is achieved. The accuracy of the surface estimate is directly related to the size of the breast and the beamwidth of the antenna [41]. The rows are rotationally offset so the closest antenna three rows away is aligned to the element of interest [41]. In [41] a grid of antennas is placed beneath the breast to ensure the model is illuminated sufficiently. This grid consists of a 5x5cm plane of 25 elements separated by 1 cm in either direction. Fig. 2.8 shows an example of a TSAR scan pattern. Figure 2.8: Data Collection Scheme Employed by TSAR System: The scan consists of 9 rows of up to 32 antennas which are represented as dots. The grid antennas are shown in red. 18 RASOR The RASOR data acquisition system is multistatic and places the patient in the prone position [42]. The sensor array is comprised of 16 wideband stacked patch antennas that are positioned on a hemispherical surface as presented in Fig. 2.9. The aim of the curved array is to simulate the shape of a pendulant breast and to design a system that can be used in preliminary clinical trials [42]. Using this sensor array the breast is scanned twice, leading to 120 recorded signals each trial [43]. The second scan is conducted at a slightly offset angle from the first, and a pulse with frequency content between 4 and 9 GHz is used to illuminate the breast in both cases [44]. This system is currently being used for preliminary clinical trials. Figure 2.9: RASOR Data Acquisition: Data is collected in the prone position using 16 wideback stacked patch antennas [45]. 2.1.3 Comparison of Microwave Breast Imaging Systems Overall each radar based microwave imaging system has advantages and disadvantages over the others; however, all the systems face the challenge of reducing reflections from the skin. The image generation technique employed by the MIST system 19 appears to be the most advanced as it compensates for frequency dependent propagation effects and is able to discriminate against artifacts and noise [16]. Some of the major draw backs of the findings reported for this system is that the realistic models have not been validated and that the data have been collected primarily in a simulated setting using 2-D models that are excited with current sources [27]. The TSAR system has reported results for 3-D models that are excited using two different antennas and a preliminary clinical prototype has been developed [41, 34, 31]. However, the image generation technique does not account for dispersion, and further processing is required after skin subtraction to isolate the tumour reflection [18]. One of the major advantages of the RASOR system is ther multistatic data acquisition scheme, as it allows collection of more information about the breast. This system has demonstrated quite a bit of success, as initial clinical results on breast cancer patients have been reported[45]. Nonetheless, these clinical results have not been validated using another modality or any other medical means. Multistatic systems also face the challenge of developing and managing more complicated electronics. The RASOR data acquisition system is also plagued with the fact that it is not robustly adaptive to the dimensions of the patient’s breast. The manner in which the data is collected dictates the characteristics of the received signals, which consequently defines the skin subtraction challenges. The first property that must be considered when reducing the skin response is whether the scan is monostatic or multistatic, as this controls the relationship between the recorded signals. The type of antenna and the distance between the antenna and the breast also has a significant effect on the processing requirements. Table 2.1 20 summarizes the basic characteristics of all the data collection schemes discussed. As the TSAR and MIST systems have monostatic data acquisition systems that place the sensors a distance away from the breast, their resulting skin subtraction problems are similar. As a result, more focus is placed on their proposed techniques for reducing the skin reflection. System MIST TSAR RASOR 2.2 Table 2.1: Data Acquisition Summary Antenna Antenna # of Excitation Scan Type Model Position Signals Frequency 1 Monostatic Current On 17 1-11 GHz Source Breast 1 Distance Monostatic Current From 17 1.67-11.1 GHz Source Breast 1 2 cm Monostatic Wu King From 64-288 1-10 GHz Dipole Breast 16 2 cm Multistatic Stacked From 120 4-9 GHz Patches Breast Skin Reflection Removal One of the main challenges of radar based microwave imaging is the fact that the largest reflection recorded is caused by the skin. This reflection occurs as a result of the incident pulse propagating from the immersion medium to the skin surface and traversing into the breast interior, as demonstrated in Fig. 2.10(a). Previous studies have reported that the skin reflection can be on the order of up to 40-60 dB larger than the tumour response in 2-D and simple 3-D numerical breast models [14, 15]. 21 Several factors make the isolation of the skin response difficult. First, the duration of excitation pulse is longer than the time it takes to travel through the skin. For example, the time it takes to travel through a 2 mm thick layer of skin with a relative permittivity of 36 is 0.04 ns while the full-width half-maximum of the TSAR pulse is 0.17ns. Consequently, the reflections caused at the skin interfaces overlap and extend in time. Second, the curved geometry of both the antenna beam and the breast model gives rise small responses that are recorded after the first dominant response. As a result the skin response is comprised of a dominant reflection that occurs early in the signal and a secondary reflection that arises throughout the rest of the signal. The tumour response is consistently masked by the secondary response, and is only embedded in the dominant reflection when the tumour is close to the skin. Dominant Response Secondary Response 4000 70 Recorded Signal Tumour Response 2000 60 Voltage 40 30 −2000 20 −4000 Voltage 50 0 10 −6000 −8000 0 2 4 1 6 −10 x 10 (a) Reflection of Incident Pulse 1.5 2 Time (s) −10 2.5 −9 x 10 (b) Dominant & Secondary Skin Responses Figure 2.10: Skin Reflections From 3-D Realistic Breast Model: The illumination of a simple breast model is shown in (a) while the resulting signal from exciting a simple breast model is seen in (b). The model used to record the signal in (b) is comprised of a 2 mm skin-layer that encloses a 6 mm spherical tumour embedded in homogeneous fatty tissue. These factors demonstrate that the skin response cannot be sufficiently removed 22 using time gating. Fig. 2.10(b) clearly illustrates this phenomenon as a recorded signal acquired by scanning a simple 3-D breast model is plotted, along with its corresponding isolated tumour response. A skin subtraction algorithm that effectively removes the dominant skin reflection and subdues the secondary response is needed. The algorithm should also limit the amount of distortion caused to the tumour response, as another signal processing technique is applied to the skin-subtracted signals to isolate this reflection. Consequently, it may not be possible to successfully image the tumour immediately following skin subtraction. In general, previously proposed algorithms have created estimates of the skin responses and subtracted these estimates from the recorded signals. Currently, there are three major skin subtraction methods that have been documented for radar-based microwave breast imaging: Calibration[46, 43], RLS/Woody [17], and Modified RLS [16]. The calibration technique is employed by the RASOR system. This method simply subtracts the data acquired from two different scans which are obtained by slightly offsetting the antenna positions. While effective for RASOR, this approach may only be applied to monostatic systems with closely spaced antennas. The RLS/ Woody and Modified RLS, on the other hand, are designed for monostatic systems and are also applicable to multistatic systems [47]. These methods use adaptive filter theory to estimate and remove the skin responses from the recorded signals. The algorithms have been successfully applied to a cylindrical 3-D model and to 2-D MR based breast models, respectively. To date, there have been no reports of successfully subtracting the skin reflections from a realistic 3-D breast model. Consequently, these two algorithms will be explored in detail in order to ascertain the feasibility of expanding their utility to 3-D models. For simplicity, the 23 mathematical descriptions of these techniques denote vectors using lower case bold letters and matrices with upper case bold characters. 2.2.1 RLS/Woody The basic concept of this skin subtraction method is to apply the recursive least squares (RLS) algorithm over the dominant segment of the skin response and the Woody averaging technique over the rest of the signal [17]. To fully understand this approach, the RLS and Woody averaging algorithms are detailed, and a description of the combination of these two techniques is given. The motivation behind integrating two signal processing approaches is also outlined. RLS The goal of RLS is to minimize the mean square error (MSE) of the estimate of a target signal. To gain insight into the application of this algorithm to radar-based microwave breast imaging, the ideal recorded signal is modeled by x = s + i + t, (2.1) where s, i, and t correspond to the responses from the skin, the interior of the breast and the tumour. The aim of the RLS is to estimate the portion of the target signal (x) that is dominated by the skin response (s) using the signals from all non-target antennas. Once the estimate is computed, it is subtracted from the target signal so that s is reduced. Eq.(2.2) presents this process mathematically as: x[n] − ~uT [n]w ≈ t[n] + i′ [n] (2.2) ~ is a vector of filter weight coefficients, and ~u[n] is the input vector. Alwhere w though the filter is not designed to remove i[n], it is expected to partially distort the 24 reflections from the interior of the breast leaving i′ [n]. It should also be noted that the filter will dramatically distort the tumour response if it overlaps with the dominant skin reflection seen early in the signal. Consequently, this algorithm cannot be employed when the tumour is in close proximity to the skin. The input vector at a given time step n is defined with the use of the data collected from all other antennas by concatenating the values of the signals at the same time step. Based on adaptive signal processing theory, the filter weights are calculated using the autocorrelation (R) of the input and the cross-correlation (~p) of the input and the target signal. Eq.(2.3), which is known as the Weiner-Hopf or the Normal equation, demonstrates the relationship between the filter weights, the autocorrelation matrix, and the cross correlation vector [48]. R~ w = ~p (2.3) The optimal filter tap weights are found by solving the Normal equation; however it should be noted that terms (R) and (~p) of Eq.(2.3) are statistical measures. These measures are not practical to estimate because of the significant number of scans the same patient would have to undergo to gain enough information to calculate them. Consequently, the RLS algorithm models these statistics with deterministic timeaveraged values, and the filter weights are found at every time step using a factor called the Kalman gain vector [48]. The derivation of the Kalman gain vector is found in Appendix B. Eq.(2.4) is the filter weight update equation which relates the filter coefficients (~ w[n]) to the Kalman gain vector (~k[n]) as follows ~ [n] = w ~ [n − 1] + ~k[n]α[n] w (2.4) ~ [n − 1] is the weight vector of the past iteration while α[n] is the a priori where w 25 error before the current weight. α[n] is characterized by Eq.(2.5), α[n] = x[n] − ~uT [n]~ w[n − 1] (2.5) where x[n] is the value of the target signal at the current time step. The final overall a posteriori error (e[n]) is demonstrated in the following equation: e[n] = x[n] − ~uT [n]~ w[n] (2.6) which is the equivalent to the residual of the skin estimation process from a skin subtraction standpoint [48, 17]. Woody Averaging The Woody Averaging algorithm estimates the skin response as an average of all the signals recorded at non-target antennas. The average is calculated over time shifted and scaled the data so that misalignments and deviation in amplitudes seen between non-target signals are compensated. The time shifts seen between the signals are caused by variations in the distance between the antenna and the breast surface. The variations in the relative distance between the antenna and the skin are also partially responsible for the amplitude fluctuations observed at different signals. Once the signals have been sufficiently aligned and scaled, they are averaged over every time step, and the average is subtracted from the target signal. The resulting signal is an estimate reflection profile of the interior of the breast. The normalized cross correlation is used to establish the appropriate time shift and scaling factor for a given non-target antenna. As the results of the cross correlation are primarily driven by the skin reflection no windowing is required to isolate the skin response. The process starts by calculating the correlation between the 26 target signal and all other signals originating from antennas on the same row. The scaling factor is then calculated by: • normalizing the minimum and maximum correlation values by the minimum and maximum autocorrelation values of the target antenna, and • averaging the normalized minimum and maximum correlation value vectors Following this calculation the signals are scaled and time-shifted to the time step where the maximum correlation is achieved. The correlation between the target antenna and the time-shifted signals is then calculated again and the process is repeated until convergence is achieved. The average of the resulting signals at each time step is taken and subsequently subtracted out of the target signal. Fig. 2.11 shows an example of the application of the Woody Averaging technique in a flow chart. As demonstrated by Fig. 2.11, the Woody Averaging approach generates a good estimate of the first major peak seen in the skin response; however, the estimate of the second peak is not very accurate. Consequently, the skin-subtracted signal has a significant residual skin response that dwarfs the tumour reflection; however, the Woody Averaging does aid in reducing some of the secondary reflections. 27 Figure 2.11: Woody Averaging Skin Subtraction Flow Chart Combining RLS With Woody Averaging The RLS/Woody algorithm first applies the Woody Averaging method outlined above to the target signal. The original target signal is then segmented such that the portion of the signal dominated by the skin reflection is isolated from the reflection profile of the interior of the breast [17]. This isolation is conducted using an estimate of the skin location and thickness relative to the target antenna. Once the segmentation is complete, the RLS algorithm is applied to the skin-dominated portion of the signals. The results of the Woody Averaging algorithm are used as the estimate of the reflection profile of the interior of the breast. Consequently, the 28 resulting Woody Averaged signal for time steps beyond the skin-dominant window is concatenated to the results of the RLS algorithm. The motivation behind this approach is to combine the strengths of both algorithms while overcoming their respective weaknesses. The Woody averaging technique maintains the tumour response, while leaving large residual traces of the skin reflections. Conversely, the RLS eliminates the skin response while completely distorting tumour reflections. Solely applying the RLS to the segment of the data that is dominated by the skin response would result in a large discontinuity in the filtered signal and the tumour response would remain masked. Consequently, both Woody averaging and RLS algorithms are required [17]. Fig. 2.12 further reinforces the need for combining both algorithms by showing the results of applying them to a recorded signal. Woody Averaging 500 RLS/Woody 30 Tumour Reflection Tumour Reflection 25 0 20 Voltage 0.5 1 1.5 2 2.5 3 −9 x 10 RLS 10 Voltage 15 −500 0 10 5 0 5 0 −5 −5 −10 −10 0 0.5 1 1.5 Time (s) (a) 2 2.5 3 −9 x 10 −15 0 Discontinuity 0.5 1 1.5 Time (s) 2 2.5 3 −9 x 10 (b) Figure 2.12: Comparison of Woody Averaging and RLS/Woody: The results of applying the Woody Averaging and the RLS algorithm, to a signal recorded from a simple 3-D realistic breast model are shown in (a) while the results of applying the RLS/Woody algorithm to the same signal are shown in (b). The RLS/Woody has also been tested on experimental data wherein the signals 29 are collected by illuminating a cylindrical breast phantom with a monopole antenna at serveral locations around the model. The diameter of the breast phantom is 10 cm and the skin layer is 2 mm thick. The dielectric properties of the materials used to create the model, mimic those of skin and fat [17]. By examining the peakto-peak ratio between the skin-subtracted data and the original skin response, the effectiveness of the skin reflection removal is measured [17]. The peak-to-peak ratio for the Woody averaging technique is found to be -26.33 dB while the results for the RLS/Woody are far better at a value of -107.87 dB [17]. A brief comparison of the RLS/Woody and the Modified RLS conducted in [49] shows that both algorithms generate similar peak-to-peak ratios. The main disadvantage of the RLS/Woody is that the concatenation of the two algorithms results in a discontinuity in the signal. As the reflections recorded from 3-D models are more complex, the RLS/Woody becomes a less attractive option because it adds a significant signal processing artifact. As a result the application of the RLS/Woody to signals recorded from a realistic 3-D breast models is not explored further. 2.2.2 Modified RLS The Modified RLS skin subtraction algorithm assumes that the skin responses measured at all antenna sites are similar but not identical due to varying skin thicknesses and breast property heterogeneity. Similarly to the RLS, the response calculated for a particular target antenna is modeled using a filtered combination of the skin responses calculated at all the other receivers. The major difference between this approach and the RLS is that the filter is not updated at every time step and it is 30 applied over the entire duration of the signal. In terms of the ideal signal model in Eq.(2.1), the goal of this process is to calculate a set of filter weights, q, and apply them to an input vector, b[n], so that the skin response (s) is adequately estimated. Once this estimate is acquired, it is subtracted from the target signal (x) to reduce the skin reflections as demonstrated by Eq.(2.7). x[n] − qb[n] ≈ t[n] + i′ [n] (2.7) The filter weights are established by minimizing the mean square error of the estimate of the skin-dominated portion of the target signal. A single set of filter weights is established by time averaging the statistics of the input vector over the skin dominant window (SDW). For a given time step n a small window is used to select the data within the SDW that is used to compute the average. Consequently, the target signal at time step n is approximated as an averaged weighted combination of the data observed at all other antennas over the SDW. Before any filtering takes place, several steps are required to define the input vector (b[n]). The first step establishes a time window (m) where the skin response dominates the target signal, as seen in Fig. 2.13(a). Following the definition of window m, a smaller window (2J + 1) is defined and used to form the input vector along with the data collected from all other receivers. At time step n, the contribution to the input from a non-target antenna a is given by 2J + 1 successive samples centered around n, as shown in Fig. 2.13(b). This step is mathematically represented by b~a [n] = [ba [n − J], . . . , ba [n], . . . , ba [n + J]]T . As demonstrated in Fig. 2.13(c), the overall input at the time step of interest (n) is defined by concatenating the input contributions of all other non-target receivers to create one large column vector. 31 Mathematically, the input vector at time step n when filtering the skin response at T ~T antenna 1 is represented by ~b2N [n] = [~bT 2 [n], . . . , bN [n]] . To set the input vector for the following time step, the small window is shifted by one and the concatenation process is repeated (Fig. 2.13(d)). With the input vector established, the minimization problem may be defined using the following cost function ~q = argmin q noX +m−1 n=no 2 T~ ~ x [n] − q b [n] 1 . 2N (2.8) In Eq.(2.8) x1 [n] is the signal at time step n of target antenna 1, ~q is the filter weight vector, ~b2N [n] is the input vector for target antenna 1 and time interval n = n0 to n0 + m − 1 represents the SDW. As ~b2N [n] is generated using the signals from all non-target antennas, the Modified RLS approximates the skin reflection as weighted combination of all other signals, averaged over the SDW. Once again the Normal equation (R~q = ~p) is used to solve the minimization problem; however, the cost function used is solely based on data originating from the skin-dominant window. Consequently, the time-averaged autocorrelation matrix (R) and the crosscorrelation vector (~p) are calculated as follows, n +m−1 1 oX ~ b2N [n]~bT2N [n] R= m n=n (2.9) o n +m−1 1 oX ~ ~p = b2N [n]x1 [n]. m n=n (2.10) o The above expressions are very similar to the autocorrelation matrix and the cross-correlation vector used in the recursive least squares problem, where statistical means are replaced with time averages [48]. The optimal filter weights are only 32 (a) (b) (c) (d) Figure 2.13: Modified RLS Input Vector Initialization Flow Chart 33 achieved if the data input into the above approximation is stationary. Although, the data does not strictly adhere to this condition, it does appear to be locally stationary over the small window (2J + 1). Consequently, the estimated filter weights may not be optimal. Because of, the high degree of similarity between the received signals, matrix R tends to be ill-conditioned. To remedy this issue, the autocorrelation matrix is replaced with a low rank approximation using eigenvalue decomposition as follows, Rg = g X λj ~uj ~uTj . (2.11) j=1 where λj and ~uj represent the eigenvalues and eigenvectors respectively and g denotes the number of significant eigenvalues. This low rank approximation simply discounts eigenvectors associated with small eigenvalues, as these vectors are considered to be associated with noise. Using the low rank approximation of the autocorrelation matrix (Rg ) the filter weights are established. The filter is then applied to the entire target signal, using the same method described above to determine the input. The filtered signal for antenna 1 (r1 [n]) can therefore be represented by the following equation r1 [n] = x1 [n] − ~qT ~b2N [n]. (2.12) The small sliding window incorporates data that may correspond to tumour responses from the non-target antennas into the input vector. Consequently, this windowing may inadvertently corrupt the tumour response of the filtered target signal. In an attempt to correct this distortion, a feedback filter is applied to the resulting signals. The feedback filter assumes that the filtered signals are a good approximation of the tumour response. Consequently, the distortion seen at antenna 1 may 34 be alleviated by applying the same filter weights to the resulting target signals as follows, r˜1 [n] = r1 [n] + ~qT ~r2N [n] (2.13) where ~r2N [n] is an input vector created with the output of the initial filtration process. The Modified RLS algorithm is summarized by the the flow chart shown in Fig. 2.14. Figure 2.14: Flow Chart of Original Modified RLS: Summary of the steps required to skin subtract the data using the Modified RLS. Fig. 2.15 presents the results obtained when applying this skin subtraction technique to the 2-D breast model shown in Fig. 2.6. Examining Fig. 2.15 reveals that skin response has practically been eliminated while the late time response appears to be preserved. In [27] the Modified RLS is found to perform equally as well on supine and prone breast models. This algorithm is also compared to an idealized skin subtraction technique wherein the signals acquired from the same breast model with and 35 Figure 2.15: Application of the Modified RLS to Simulated Data: The algorithm is applied to 17 signals acquired from a 2-D MR-based realistic breast model. The parameters J, m, and g are set to 3, 20, and 15, respectively [16]. without a varying interior breast property distribution are subtracted [27]. The performance of the Modified RLS is found to be superior to this idealized method as less energy is detected by the Modified RLS near the breast surface [27]. This outcome is a direct result of the breast property heterogeneities arising at the interior interface between the skin and the breast tissue. These breast property heterogeneities cause large reflections in the early-time signal that cannot be sufficiently suppressed using the idealized method. Although the Modified RLS has not been tested on signals recorded from realistic, 3-D breast models it appears to have to most potential, for three major reasons. First, the Modified RLS does not cause a discontinuity in the resulting signals, second, it has been demonstrated to be effective on a variety of 2-D breast models and third, it is applicable to multistatic systems [47]. 36 2.3 Summary Radar based microwave imaging has been proposed as a complementary breast cancer detection technique. This modality illuminates the breast with an ultra-wideband pulse, records the ensuing reflections at several locations around the breast and processes the reflections to find the largest scatters within the imaging volume. The significant scatterers are detected by applying a pipeline of signal processing techniques, that isolate the tumour reflection, and focusing the resulting signals to create an image. One of the first steps in the signal processing pipeline is to remove the overwhelming reflection caused by the skin that masks the tumour response. Although previously proposed skin subtraction algorithms have been successfully tested on 2-D scenarios, a method applicable to monostatic data collected from a realistic 3D breast models has yet to be reported. The next chapter explores the feasibility of expanding the Modified RLS to skin subtract signals from complex 3-D breast models. 37 Chapter 3 Neighbourhood Modified RLS This chapter explores the expansion of the Modified RLS algorithm to reduce skin reflections recorded from a simple, realistic, 3-D breast model. First, the performance of the original method in 3-D is examined, followed by an outline of the performance measures used to quantify the effectiveness of potential skin subtraction algorithms. Next, the development of an input data selection technique termed neighbourhood selection is described. The procedure to estimate the filter weights is then explored by investigating the effect of varying the width of the SDW and incorporating a previously reported adaptive rank selection approach into the Modified RLS. Finally a cost/benefit analysis of including the feedback filter in the neighbourhood Modified RLS algorithm is conducted. 3.1 Application of Modified RLS to 3-D Breast Model The performance of the Modified RLS in 3-D is evaluated by examining the results of applying the algorithm to simulated data recorded from a simple realistic numerical breast model. No major modifications are made to the Modified RLS as the goal is to establish the feasibility of applying this algorithm to data recorded from a 3-D model. To gain insight into the limitations of the algorithm, images are produced using the resulting skin-subtracted data and potential areas of improvement are identified. 38 3.1.1 3-D Breast Model The breast model is derived using a Magnetic Resonance (MR) scan by converting the MR pixel intensities to dielectric property values (simple breast model described in section 2.1.1). The model examined in this chapter consists of a 2 mm thick layer of skin that bounds non-dispersive homogeneous fatty tissue. A 6 mm diameter tumour is embedded in the fat. The dielectric properties of the various tissues included in the model can be seen in Fig. 3.1 as a cross sectional slice of the model’s relative permittivity map is shown. Relative Permittivity Map 50 120 45 110 40 100 35 z (mm) 90 30 80 25 70 20 60 15 50 10 40 20 5 40 60 80 x (mm) 100 120 Figure 3.1: Dielectric Properties of Simple Breast Model: The relative permittivity, and conductivity of the skin (36, 4 S/m), fat (9, 0.4 S/m), tumour (50, 4 S/m), and immersion medium (3, 0 S/m), are established using the first method outlined by the TSAR group in section 2.1.1. The data is simulated using the FDTD method. A single simulation consists of a breast model that is being illuminated with a 1 to 10 GHz differentiated Gaussian pulse using a resistively loaded dipole antenna. A complete data set is acquired by translating the antenna to a number of different locations around the breast and recording the resulting simulated fields. An elliptical scan pattern consisting of 245 different illumination sites is used in conjunction with a 5x5 grid of antennas at the 39 base of the scan. The 245 antennas are placed in 9 rows with up to 32 antennas per row, with the minimum number of antennas on a given row being 20. 3.1.2 Processing Signals From 3-D Breast Model To effectively apply the Modified RLS to data recorded from a 3-D model, the differences between the 2-D and 3-D problems must be identified. There are three major distinctions between the data recorded from the 2-D models examined by the MIST system and the data recorded from the 3-D model described in section 3.1.1. First, the number of sensors used to acquire the data has increased significantly, as a result of illuminating a 3-D volume. Second, the distribution of the sensors within the scan pattern has changed, as there is no longer a single row of receivers. Third, the excitation used by the MIST & TSAR systems are marginally different; consequently, the duration of the skin responses are not the same. In effort to apply the Modified RLS to the data recorded from the 3-D model presented in section 3.1.1, while preserving the skin subtraction process described in [16], only one modification is applied. This modification deals with the selection of the SDW as there was no methodical description on how this was previously accomplished. Consequently, the skin subtraction approach, starts by selecting a target antenna τ , that will be filtered using combination of the signals recorded at all other antenna. Next, the skin dominant portion of the target signal (xτ ) is identified using an adaptive skin-dominant window selection technique, which will be described further in section 3.4.1. The input vector (b~τ M [n]) used to filter the target signal xτ at time step n is defined with the use of segments of all other signals. The contribution of any given non-target antenna i to the input vector b~τ M at time step 40 ~i [n]), is n is 2J + 1 successive samples centered about n. This contribution vector (b represented mathematically by b~i [n] = [bi [n − J], . . . , bi [n], . . . , bi [n − J]]T . (3.1) The overall input vector b~τ N at time n is a column vector that is a concatenation of the contribution vectors bi~[n], from all other non-target antennas. This is described mathematically by, T ~T ~ ~ ~T T T b~τ M [n] = [b~T 1 [n], b2 [n], . . . , bτ −1 [n], bτ +1 [n], . . . , bM [n]] , (3.2) where there are M antennas in the entire scan pattern. The input vector for the following time step is simply established by shifting the center of the 2J +1 successive time over by one. To remain consistent with [16] the value of J is set to 3. Once the input vector is established, the filter weights are computed such that the following cost function is satisfied. ~q = argmin q noX +m−1 2 T~ xτ [n] − ~q bτ M [n] n=no (3.3) Consequently, the filter coefficients (~q) are designed to minimize the mean square error between the target signal xτ , and the weighted input vector q~T b~τ M over the skin dominant portion of the target signal. To solve this minimization problem, the autocorrelation matrix (R) of the input vector, and the cross correlation vector (~p) between the target signal and the input vector are computed and averaged over the SDW as follows n +m−1 1 oX ~ R= bτ M [n]~bTτM [n] m n=n o (3.4) 41 n +m−1 1 oX ~ ~p = bτ M [n]xτ [n]. m n=n (3.5) o The variables no and no + m − 1 in Eq.( 3.4- 3.5) are the time steps associated with the first and second boundary of the SDW, respectively. As there is a high degree of similarity between a subset of the recorded signals, the resulting autocorrelation (R) tends to be ill-conditioned. As a result, eigenvalue decomposition is used to estimate a low rank approximation of R. To remain consistent with the procedure outlined in [16] the rank of the averaged autocorrelation matrices calculated for every target antenna is assumed to be 15. Using simple Weiner filter theory, the filter weights are computed as ~q = R−1 p, 15 ~ (3.6) where R−1 15 is the inverse of the low rank approximation of the averaged autocorrelation matrix. The filter weights are then applied to the entire target signal, using the same technique to establish the input vector. Eq.(3.7) demonstrates this filtration process rτ [n] = xτ [n] − ~qT ~bτ M [n] (3.7) where rτ [n] is the skin-subtracted signal at time step n. To alleviate any distortion to the tumour response, a feedback filter is applied; wherein the same filter weights are applied to the resulting filtered signals and added back into the skin-subtracted signal. Consequently, the final skin-subtracted signal (r˜τ ) with feedback is computed as r˜τ [n] = rτ [n] + ~qT ~rτ M [n]. (3.8) 42 After all the signals have been skin-subtracted, a simple delay-and-sum focusing algorithm is applied to the data in order to create an image. The first step of the focusing algorithm is to divide the imaging volume into small cubes called voxels. Once this step is accomplished, a single voxel is selected as the focal point. The time required for a signal to travel from a given antenna to the focal point is then calculated using estimates of the skin location and thickness, as well as assuming dielectric properties for the immersion medium, skin, and interior of the breast. The travel times from the focal point to all antennas are calculated and the corresponding components of the skin-subtracted signals are summed. The focal point is then scanned through the entire imaging volume, and the same delay-and-sum procedure is repeated. Finally, the summed value calculated for every voxel is squared, and the resulting intensities are normalized by the maximum value. The voxels with the largest intensities are assumed to correspond to the largest scatterer in the imaging volume, which in turn should coincide with the tumour location. 3.1.3 Results & Discussion Fig. 3.2 demonstrates some of the images obtained when the above focusing algorithm is applied to skin-subtracted data acquired using the Modified RLS. The image in Fig. 3.2(a), shows the plane where the largest scatterer is found while Fig. 3.2(b) represents the plane where the tumour should be found. The tumour localization error is 26 mm indicating that the tumour was not successfully detected. 43 Plane @ z = 85 mm 1 20 1 40 0.8 40 0.8 0.6 60 0.4 80 0.2 100 0 120 50 100 y (mm) (a) 150 200 x (mm) x (mm) Plane @ z = 98 mm 20 0.6 60 0.4 80 0.2 100 0 120 50 100 y (mm) 150 200 (b) Figure 3.2: Resulting Images When Modified RLS is Applied to a 3-D Model: Images generated when the original Modified RLS is applied to a signals collected from a simple realistic 3-D breast model. The plane where largest scatterer is found is shown in (a) while the plane where tumour should be found is displayed in (b). Failure to detect the tumour is expected, as the skin response is not sufficiently reduced. The magnitudes of the residual dominant skin reflections are slightly larger than the magnitudes of the tumour responses, while the residual secondary responses dwarf the tumour reflections. The inadequate removal of the skin response is further demonstrated by the peak-to-peak ratios between the skin-subtracted data, and the original skin response. Previously, ratios as low as -100 dB have been reported for the Modified RLS; however, a significantly larger ratio of -51.49 dB is found in this scenario [49]. Consequently, the algorithm must be improved in terms of three major facets. First, the amplitude of the residual dominant reflection must be smaller than the amplitudes of the tumour response. Second, secondary skin response should be subdued more effectively in order to increase the potential of successfully detecting the tumour. Third, the approach must pose less of a computational burden. The Modified RLS is comprised of ten steps, as shown in Fig. 3.3. In effort to successfully skin subtract this model, the four highlighted pink steps in Fig. 3.3 are examined and revised. Selecting a subset of antennas as the input to the algorithm is 44 the first modification examined as this may dramatically decrease signal processing time. Currently, it takes over 5 hours to estimate the skin reflections for 245 antennas on a 2.81 GHz machine with 1.96 GB of RAM. Next, the effect of varying the width of the skin-dominant window is investigated to ascertain the algorithm’s dependence on this parameter. The low rank approximation is also explored, for it is hypothesized that by making this process adaptive, a better skin response estimate is produced. Finally, the value of including the feedback filter to remove distortion caused by strong neighbouring tumour responses is evaluated. Figure 3.3: Flow Chart of Modified RLS: The steps highlighted in pink are examined further in this chapter. To determine the effect of altering these four steps, standard performance metrics are developed. Section 3.2 outlines the implications of these measures and how they are calculated. 45 3.2 Performance Measures There are two general approaches to assessing the results generated by applying variations of the algorithm to the same data set. The first is to study the filtered signals, while the second is to analyze the images that are created by focusing the processed signals. 3.2.1 Single Signal Analysis Several signals are used to compute the single signal performance measures. The following is a list of those signals and how they are calculated. 1. The calibrated data (CData) is determined by subtracting the signals recorded from two different simulations. The first simulation is of an antenna that is in an immersion medium and is radiating the differentiated Gaussian pulse used by the TSAR group as demonstrated by Fig. 3.4(a). The second simulation seen in Fig. 3.4(b) is identical to the first except that a breast model is placed in front of the radiating antenna. The breast model is made up of a skin layer that encloses an interior tissue distribution, including a tumour. The second simulation is repeated several times while placing the antenna at different locations relative to the breast. The overall goal is to remove the original signal that excited the antenna, as well as any other characteristic reflections of the antenna from the recorded signals. Consequently, the calibrated data is computed by taking the difference between the signal recorded from the first and the second set of simulations. 2. The no tumour data (NoTumData) is established in the same manner as 46 (a) (b) Figure 3.4: Calculating the Calibrated Data: The results from simulating the scenario seen in (a) are subtracted out from the signal recorded from simulating the layout in (b). This processing is repeated while placing the antenna at several different locations around the breast model. the calibrated data, except that the breast model used in the second set of simulations does not contain a tumour. It is also important to note that the antenna placements used are identical to those used for the calibrated data simulations. 3. The known tumour responses (Tactual ) are computed by subtracting the no tumour data from the calibrated data of the same breast model. The aim is to isolate the reflections caused by the tumour. 4. The skin-subtracted signals (NoSkData) are calculated by applying the skin subtraction algorithm to the calibrated data while the residual skin responses (ResSkData) are established by filtering the no tumour data. 5. The tumour estimates (Test ) are obtained by taking the difference between skin-subtracted signals and the residual skin response. 47 There are three major concerns when establishing the quality of a skin-subtracted signal: skin suppression effectiveness, tumour response preservation, and clutter reduction capabilities. The tumour response preservation indirectly incorporates the clutter reduction and skin suppression abilities of the algorithm. As a result, four measures are used to quantify the tumour preservation while only one measure is used to describe each of the skin suppression and clutter reduction. 1. The tumour-to-skin ratio (T um2SkinRatio) measures the suppression of the dominant skin reflection by comparing the amplitudes of the skin-subtracted signals and their respective known tumour responses. The amplitude is measured by calculating the peak-to-peak (p2p) value of the signal, which is found by taking the voltage difference between the largest peak and trough. The p2p value of the skin-subtracted signal (p2pSkSubSDW ) is computed using the data found within the skin-dominant window, while the p2p of its corresponding tumour response (p2pTactual ) is calculated over the entire duration of the signal. Once the p2p values have been established, a skin-subtracted-totumour ratio is computed as follows. T um2SkinRatio = 20 ∗ log p2pTactual p2pSkSubSW D (3.9) To evaluate the results of this metric, a comparison is made between the original tumour-to-skin ratio and the tumour-to-skin ratio after skin subtraction. The p2p of the original signal is calculated using the calibrated data. The original tumour-to-skin ratios are typically negative, while a successful skin removal trial would generally result in a positive ratio. Previously, tumour-to-skin ratios ranging from 13.38 to 33.85 dB have been reported 48 after filtering data that originally had ratios of -58.57 to -52.37 dB [49]. 2. The energy ratio quantifies both tumour preservation and clutter reduction by comparing the residual skin responses to the no tumour data. This metric measures the suppression of the secondary skin response; consequently, this criterion is only applied to the data recorded after the second edge of the skin-dominant window, as demonstrated by the red rectangle in Fig. 3.5. Eq.(3.10) is used to calculate the measure. An energy ratio greater than 100 indicates that unwanted signal processing artifacts are corrupting the signal. Therefore, a good energy ratio is close to zero. PN ResSkData2 [n] EnergyRatio = PNn=mo × 100 2 [n] N oT umData n=mo (3.10) Figure 3.5: Calculation of Energy Ratio: The energy ratio between the residual skin response and the no tumour data, is calculated over the portion of the signals contained within the red rectangle, shown above. 49 3. The offset is a tumour preservation metric that measures the time shift between the known tumour response and the tumour estimate. The time shift is established by computing the cross correlation between the two signals and determining the ensuing lag. This measure reveals whether there is another response that strongly resembles the known tumour response, within the tumour estimate. Ideally, the time shift value would be zero. 4. The normalized mean square error (MSE) is another tumour preservation criterion that calculates the difference between the known tumour response and the tumour estimate. This is calculated by normalizing both signals to their maxima and averaging the point-by-point error between both signals over time. Once again a MSE of zero is ideal as it indicates that the shapes of the two signals are identical and the alignment is perfect. Eq.(3.11) demonstrates how this measure is calculated 2 N Test [n] Tactual [n] 1 X M SE = − N n=1 max (Test ) max (Tactual ) (3.11) where N is the length of the entire signal, and n is the time step of interest. Fig. 3.6 compares two tumour estimates with their corresponding known tumour responses. The MSE calculated for the signals seen in Fig. 3.6(a) is 0.026, while the MSE for the signals in Fig. 3.6(b) is 0.046. 50 Antenna 21 8 Antenna 214 1.5 Tactual 6 Tactual Test Test 1 0.5 2 Voltage (V) Voltage (V) 4 0 −2 0 −0.5 −4 −1 −6 −8 0 0.5 1 1.5 Time (s) (a) 2 2.5 3 −9 x 10 −1.5 0 0.5 1 1.5 Time (s) 2 2.5 3 −9 x 10 (b) Figure 3.6: Comparison of Tumour Estimate and Known Tumour Response: The tumour estimate in (a) agrees well with known tumour response leading to an MSE of 0.026. The estimate in (b) is not accurate and has a much higher MSE of 0.046. 5. The percentage of tumours preserved is the final tumour preservation metric and it establishes if the tumour estimates are accurate. The tumour estimate is valid if the energy ratio, the offset, and the MSE adhere to certain thresholds. First, the energy ratio must be below 100, indicating that the estimate is not riddled with large signal processing artifacts. Next, a threshold of 28.8 ps is placed on the offset between the estimate and the known response, as this corresponds to a shift of 3 mm in the tumour position. The 3 mm shift is calculated assuming the wave is propagating in a non-dispersive fatty medium with a relative permittivity of 9. A 3 mm shift for this breast model is permissible, as the diameter of the tumour is 6 mm. The MSE is then considered by placing a threshold of 0.027 on its value. This constraint is set by visually inspecting the tumour estimates and their corresponding MSE . Table 3.1 summarizes the single signal analysis measures. 51 Table 3.1: Single Signal Analysis Measures Metric Definition Desired Range Tumour-to-Skin P2P ratio between Tactual & ≥0 Ratio NoSkData Energy Ratio Energy ratio between 0 ResSkData & NoTumData Offset Cross correlation lag between 0 to ≤28.8ps Tactual & Test (3 mm shift) Normalized MSE Mean squared error between 0 to ≤0.027 Tactual & Test Percentage of Tumour estimates with good 100% Tumours Detected energy ratio, offset, and MSE 3.2.2 Image Evaluation To gain further insight into the performance of the skin subtraction algorithm, images are formed using the resulting signals. Generally, the TSAR signal processing pipeline includes a clutter reduction technique; however, this step is omitted when creating the images. The clutter reduction technique computes an estimate of the tumour response using the skin-subtracted data. Consequently, the images generated without this technique are expected to be dominated by clutter. Two major factors are considered when assessing the quality of the images. The first factor is the location of the largest detected scatterer. The second is the pixel intensity difference between the two largest scattering responses in the imaging volume. In order to isolate the baseline performance of the focusing algorithm from the effectiveness of the skin subtraction process, the images produced when focusing the known tumour responses are used as references. Fig. 3.7(a) demonstrates the reference image generated for the x-y plane containing the greatest pixel intensity. To assess the quality of the tumour localization, a comparison is made between the 52 1 40 0.8 0.6 60 0.4 80 0.2 100 Zoom in of Reference Image @ Plane z = 85 mm 0.6 Extent of Tumour 60 0.4 0.2 65 0 120 50 100 y (mm) (a) 150 200 1 0.8 55 x (mm) x (mm) Reference Image @ Plane z = 85 mm 20 0 100 105 110 115 y (mm) 120 125 130 (b) Figure 3.7: Reference Image: Image generated of the x-y plane where the largest scatterer is detected using known tumour responses is shown in (a) while a zoom in of the same image is shown in (b) to demonstrate the extent of the imaged tumour. location of the greatest pixel intensity found in the reference images and the known position of the tumour. The Euclidean distance between these two sites is deemed the localization error. Any reductions in this error are due to unwanted shifts of the tumour estimates caused by the skin subtraction process. In order to calculate the second image evaluation measure, the size of the tumour response in the image domain must be estimated. This parameter is visually established based on the reference images of three planes that contain the maximum pixel value. Fig. 3.7(b) demonstrates the segmentation applied to the x-y plane of interest. As the extent of the tumour in 2-D is represented by a rectangle, its dimensions are dictated by the outermost pixel that is considered to be associated with the tumour. Pixels that have an intensity value greater than 0.1 are deemed part of the tumour, provided that the neighbouring pixel that is closest to the center of the tumour is also greater than 0.1. The center of the tumour is assumed to be located at the pixel with intensity value 1. Using the dimensions of the rectangles from all three planes, the overall extent of the tumour response is denoted by a rectangular prism in 3-D. Once the 3-D extent of the tumour is ascertained, the signal-to-clutter ratio 53 (SCR) is calculated. This metric is computed by zeroing the region associated with the tumour, and finding the next largest response in the image. These intensities are then compared using the following relationship, SCR = 10 ∗ log M axInts ClInts (3.12) where M axInts and ClInts are the largest and second largest pixel intensity observed in the imaging volume, respectively. The SCR of the images generated using skin-subtracted data are calculated by zeroing the tumour extent defined by the reference image and taking the ratio between the two largest pixel intensities. The SCR value calculated using the skin-subtracted data should typically not be higher than the ratio found using the reference image. It should also be noted that the SCR may not directly correlate with the energy ratio measure calculated during single signal analysis as the focusing process may result in coherent or incoherent addition of signal artifacts. To ascertain the success of the skin subtraction process, the first results that are examined are the focused images. Both tumour estimate and skin-subtracted data are focused, leading to two distinct sets of images. The images generated by the tumour estimate data are an indication of whether or not the tumour response is preserved, while the images from the skin-subtracted data reveal if the tumour can be detected without further processing. The image evaluation tools described above are then applied to both images. Following the image evaluations, the data sets that result in the successful detection of a tumour are identified and the individual signals are examined. Finally, the best data set is typically determined based on the SCR and the number of preserved tumour responses. Table 3.2 summarizes the image 54 evaluation parameters. Metric Localization SCR 3.3 Table 3.2: Image evaluation tools Definition Euclidean distance between maximum pixel intensities of reference image & image of interest Ratio between the pixel intensities of the two largest scatterers in image Desired Range 0 to 21 tumour diameter 0 to SCR of reference image Neighbourhood Selection The goal of the neighbourhood selection algorithm is to facilitate the successful skin subtraction of data acquired from a realistic 3-D breast model. The underlying concept is to create a more meaningful input vector for the skin reflection reduction technique. The neighbourhood selection algorithm limits the data incorporated into the input vector of a given antenna by only allowing antennas that are in its vicinity to contribute. This modification is implemented to ensure that the data used to calculate the filter weights originates from a local region of the breast with similar curvature. The degree of curvature did not present as much of a challenge in previous work, as only 2-D breast models were considered [16]. Consequently, the neighbourhood selection process attempts to mitigate the differences seen in the recorded signals caused by the realistic 3-D breast shapes. The algorithm works on the premise that creating the input vector by grouping antennas that are illuminating a similar region of the breast will generate a better skin response estimate. Both cross correlation between signals and the distance between antennas are considered in this process. Several neighbourhood selection schemes are tested, encompassing a variety of antenna distributions and potential 55 cross correlation thresholds. The prospective antenna groupings that were considered are: 1. setting a hard limit on the number of rows included, 2. restricting the number of rows and columns accepted, 3. admitting antennas that illuminate the same region of the breast, 4. a cross correlation criterion in combination with the best of the above 3 options. The combination of the Modified RLS and the neighbourhood selection algorithm is applied to the simple 3-D realistic breast model described in section 3.1. To test the proposed neighbours, the basic parameters of the Modified RLS are set to consistent values throughout the trials. The skin-dominant window and the low rank approximation of the autocorrelation matrix are established using an adaptive approach that is elaborated upon in section 3.4. The small sliding window (2J+1) on the other hand is set to 7 as reported in [16]. For all the trials presented below, the feedback filter is not applied. The performance measures described in section 3.2 are applied to the resulting data sets to establish the best potential neighbourhood. 3.3.1 Row Based Neighbourhoods The first antenna distribution schemes that are explored are based on rows. For example, a 1 row neighbourhood for a given target antenna consists of antennas within the same row. For cases where multiple rows are included, the input data comes from antennas that are in a limited number of rows centered about row of the target antenna. A 5-row neighbourhood consists of the antennas that are in the same 56 row as the target antenna as well as the antennas in the two rows above and beneath the target. When a row above or beneath the target antenna level does not exist, no exceptions are made to increase the amount of data included in the input vector. The row based neighbourhood selection techniques are described mathematically as follows. Let xk,l be the signal recorded at the target antenna that is on the k th row and the lth column of the antenna scan. The input vector for any given target antenna is defined using concatenated segments of the signals recorded at its neighbouring antennas. Let b~r,c [nJ] be a segment of a signal on the rth row and the cth column that is used to create a portion of the input vector at time step n. This vector is defined by, b~r,c [nJ] = [xr,c [n − J], . . . , xr,c [n], . . . , xr,c [n + J]]T (3.13) where 2J + 1 is the size of the small sliding windows used to create the input vector, and vector xr,c , is the signal recorded at a neighbouring antenna. The row based neighbourhoods includes all antennas that are within a certain number of rows (RT ) from the target antenna. For each row, an intermediate column vector κr [n] is defined by T ~T T~ κ~r [n] = [b~T r,0 [nJ], br,1 [nJ], . . . , br,NCr [nJ]] (3.14) where N Cr is the number of antennas on row r, and κ~r [n] is the contribution of the rth row to the input vector. An exception is made in the creation of κ~r [n] when r = k. For this case κ~k [n] is defined by T ~T T~ T~ T~ κ~k [n] = [b~T k,0 [nJ], bk,1 [nJ], . . . , bk,l−1 [nJ], bk,l+1 [nJ], . . . , bk,NCr [nJ]] (3.15) This exception ensures that the signal recorded at the target antenna is not included 57 ~ [n]) at time step n is generated in the input vector. The row based input vector (bRk,l by concatenating the column vectors (~κ[n]) as follows ~ [n] = [κT~ [n], . . . , κ~T [n], . . . , κT~ [n]]T bRk,l k−RT k k+RT (3.16) where RT guides the number of rows accepted into the neighbourhood. The conducted trials included 1-, 3-, 5-, 7-, and all row neighbourhoods, wherein RT is set to 0, 1, 2, 3, and 4, respectively. The number of antennas on every row for this particular scan pattern is also not the same. As a result, a target antenna on a given row has a different number of signals admitted into the vector. Table 3.3 summarizes the number of antennas in each row, as well as the number of antennas included in the input vector for a given level and neighbourhood. To illustrate the impact of the neighbourhood selection algorithm, images created with 1-, 3-, and 5-row neighbourhoods are shown in Fig. 3.8, while Table 3.4 summarizes the results when the image evaluation tools discussed in section 3.2.2 are applied to Fig. 3.8. Table 3.3: Neighbourhood Size for Each Row Using Row Selection Scheme Row 1 2 3 4 5 6 7 8 9 # of Antennas in Row 25 20 20 20 32 32 32 32 32 1 row 24 19 19 19 31 31 31 31 31 Neighbourhood 3-row 44 64 59 71 83 95 95 95 63 5-row 64 84 116 123 135 147 159 127 95 Upon first glance, it is clear that fewer artifacts are present in the tumour estimate images in comparison to the skin-subtracted images. This phenomenon is expected as any residual responses in the skin-subtracted signals are removed when calculating the tumour estimates. When comparing neighbourhoods, no significant differences 58 Skin-Subtracted Data 1 20 1 40 0.8 40 0.8 60 0.6 60 0.6 0.4 80 0.2 100 50 100 y (mm) 150 0.2 0 50 100 y (mm) 150 200 20 1 20 1 40 0.8 40 0.8 60 0.6 60 0.6 0.4 80 120 x (mm) 120 200 0.2 100 5 0.4 80 100 0 x (mm) 3 x (mm) 120 x (mm) 20 50 100 y (mm) 150 0.4 80 0.2 100 0 120 200 0 50 100 y (mm) 150 200 20 1 20 1 40 0.8 40 0.8 60 0.6 60 0.6 0.4 80 0.2 100 120 0 50 100 y (mm) 150 200 x (mm) 1 Tumour Estimate x (mm) # of rows 0.4 80 0.2 100 120 0 50 100 y (mm) 150 200 Figure 3.8: Focused Images of Row Neighbourhoods are seen in the tumour estimate images; however, substantial improvements are seen in the skin-subtracted images. These differences are also documented by the SCR reported in Table 3.4. Images for the all-row and the 7-row neighbourhoods are not included in Fig. 3.8, because a tumour is not successfully found in the skin-subtracted images for either of these cases. Examining Table 3.4 reveals that the highest SCR is achieved when focusing the 59 Table 3.4: Image Quality Using Row Based Neighbourhoods. SCR (dB) Localization Error (mm) Neighbourhood Tumour Skin-Subtracted Tumour skin-subtracted Estimate Data Estimate Data 1 row 3.87 0.56 1 2.24 3-row 4.15 1.01 1 2.24 5-row 4.98 1.95 1 2.24 7-row 5.36 N/A 1 N/A all-row 5.69 N/A 1 N/A tumour estimates from the all-row and the 7-row scenarios; however, these SCR are larger than the ratio computed from the reference image (5.05 dB). This indicates that the tumour estimates generated from the all-row and 7-row neighbourhoods, exhibit more deconstructive interference than the known tumour response when focused. Consequently, the resulting SCR from these neighbourhoods are deemed artificially high. The SCR of the skin-subtracted images for these two cases are not recorded because the tumour is not detected. This illustrates the all-row and 7-row neighbourhoods’ inability to sufficiently suppress the secondary skin response in order to reveal the tumour reflection. It is postulated that these neighbourhoods are unable to reduce the secondary response because they include antennas that have recorded significantly different signals into their neighbourhoods. Based on the image evaluation, the 5-row neighbourhood appears to have the optimal performance. To ensure that this data set is better than the other proposed options, the single signal analysis is conducted on the three trials that successfully imaged a tumour. As there is a significant number of antennas in the scan pattern, the results of the single signal analysis is presented in a condensed manner. The resulting tumourto-skin ratios are summarized using the average and standard deviation of this ratio 60 Tumour−to−Skin Ratio 30 original 1 row 3 rows 5 rows Skin Subtracted Ratio (dB) 55 50 20 10 45 0 40 −10 35 −20 30 −30 25 −40 20 −50 15 −60 10 0 2 4 Row 6 8 Original Ratio (dB) 60 −70 10 Figure 3.9: Tumour-to-Skin Ratios Using Row Neighbourhoods: The average and standard deviation of the tumour-to-skin ratios before and after skin subtraction are plotted for every row. The averages are represented by circles or squares while the standard deviation is demonstrated using errorbars. over every row. Fig. 3.9 demonstrates the average and standard deviation of the skinsubtracted and original tumour-to-skin ratio per row for the three neighbourhoods under consideration. The original ratio typically lies between -65 to -60 dB while the skin subtraction ratios exhibit an improvement of approximately 90 to 95 dB for all neighbourhoods. The 3-row and the 5-row neighbourhoods have the most consistent results, as fewer outliers are observed. However, the 5-row neighbourhood generates slightly better results, as it is the best case for 4 out of 9 rows. The average energy ratios per row for the three neighbourhoods are displayed in Fig. 3.10. The difference in performance is not considerable, but the variation in the typical energy ratio is fairly significant. The energy ratio for all rows except row 1 is 61 Average Energy Ratio per Row 60 1 row 3 rows 5 rows 50 Energy % 40 30 20 10 0 0 2 4 Row 6 8 10 Figure 3.10: Mean Energy Ratio Using Row Neighbourhoods less than 25%, while rows 4 through 8 have a ratio of less than 10%. The large value seen in row 1 can be explained by two factors. First, fewer antennas are used to create the input vector, as demonstrated by Table 3.3. Second, the energy ratios for 7 of the 25 antennas in row 1 are greater than 80%; a ratio greater than 60% is achieved only once in rows 2 through 9. This dramatic decrease in performance is likely due to the extreme variation in curvature seen at the bottom of the breast. At rows 4 through 7, the breast is somewhat cylindrical but at rows 1 to 3 it begins to converge towards the nipple. The effect of this shape variation is seen most prominently in the data of row 1 as a result of the antenna orientation at this level. Fig. 3.11 presents the average over every row of the offset recorded between the tumour estimate and the known tumour response. Once again the average time shift 62 Time Shift 200 1 row 3 rows 5 rows Time Steps ∆t = 0.48e−12s 150 100 50 0 −50 −100 −150 −200 −250 0 2 4 Row 6 8 10 Figure 3.11: Mean Offset of Each Row Using Row Neighbourhoods does not fluctuate significantly between the three neighbourhoods, except for at rows 2 and 7. The outlier seen at row 2 using the 3-row neighbourhood is caused by 6 antennas that have a time shift greater than 200 time steps. Interestingly, the large time shift recorded for 5 of these 6 cases is a result of clutter corrupting the latetime tumour estimate. Fig. 3.12 demonstrates two examples of this phenomenon. Although large differences are not observed between the three investigated cases, the 5-row neighbourhood exhibited no outliers, and typically had a smaller standard deviation. The 5-row neighbourhood also had an average time shift value of less than 60 time steps for 5 out of 9 rows, while this only occurred at 3 rows for the other cases. Consequently, the 5-row neighbourhood appears to be the most desirable candidate. 63 Antenna 27 0.6 Antenna 34 1 Tumour Estimate Known Tumour Response 0.8 0.4 0.6 0.2 Voltage Voltage 0.4 0 −0.2 0.2 0 −0.4 −0.2 Tumour Estimate Known Tumour Response −0.6 −0.8 0 0.5 1 1.5 2 Time (s) −0.4 2.5 −0.6 3 0 0.5 1 1.5 Time (s) −9 x 10 (a) 2 2.5 −9 x 10 (b) Figure 3.12: Examples of Late-Time Clutter: The first half of the tumour estimate is accurate while the second half contains artifacts that can be mistaken as the tumour response. As a result of this artifact the offset calculated for the signals in (a) is 451ps while the offset seen in (b) is 417 ps. Normalized Mean Square Error Over Entire Signal 1 row 3 rows 5 rows 0.04 Normalized MSE 0.035 0.03 0.025 0.02 0.015 0.01 0 2 4 Row 6 8 10 Figure 3.13: Mean MSE of Each Row Using Row Neighbourhoods 64 The averaged normalized MSE results for the row-based groupings of interest is seen in Fig. 3.13. Once again, rows with a larger number of antennas in their neighbourhoods achieved better results (with the exception of row 1). It should be noted that the original tumour-to-skin ratio for row 1 seen in Fig. 3.9 is much higher than all other rows, implying that smaller skin responses are recorded for this row. When comparing the results of the three neighbourhoods of interest, it appears that the 5-row scenario is the most desirable, as it has the best MSE for 6 out of 9 rows and has the worst performance only once. It is also observed that the lowest overall MSE was achieved at row 7. This result is expected as row 7 is perfectly aligned with the tumour elevation, which would lead to larger tumour responses. Fig. 3.14 illustrates the percentage of tumours preserved in every row and in total, for all three neighbourhoods. The differences between the performances of each neighbourhood at every row can also be explained using the three previously examined measures. Generally, the 5-row case preserves more tumours than the 1- or 3-row scenarios, with exceptions at rows 2, 3, 9. The poor performance in these cases can be attributed to the large variation in curvature of the breast at the elevations corresponding to these rows. Data collected from regions of the breast with significantly different curvatures would be more likely to be incorporated into the input vector with a 5-row neighbourhood. This phenomenon is also reflected in the energy ratio results. When considering the image evaluation and the single signal analysis of the three neighbourhoods, it is clear that the 5-row neighbourhood is the optimal candidate of those considered. Tumour Response Preservation 100 1 row 3 rows 80 5 rows 90 70 80 60 70 50 60 50 40 40 30 30 20 20 10 10 0 1 2 3 4 5 Row 6 7 8 9 Total 0 Overall Percentage of Preserved Responses Percentage of Preserved Responses per Row 65 Figure 3.14: Percentage of Preserved Tumour Responses Using Row Neighbourhoods 3.3.2 Row & Column Based Neighbourhoods In an effort to reduce processing time and to test whether the 5-row neighbourhood is optimal, neighbourhoods selected based on the number of rows and columns centered about the target antenna are investigated. As the 5-row candidate proved to be the best option explored, the number of rows tested remains static at 5 and the number of columns is varied between 2, 4 and 8. The row-column based neighbourhood selection techniques are represented mathematically as follows. The input vector for a given target signal xk,l recorded by an antenna that is on the k th row and the lth column of the scan pattern, is defined using data from antennas that are a certain number of rows (RT) and columns (CT) away from the target. The contribution to the input vector at time step n, of a neighbouring 66 antenna that is on the rth row and cth column (b~r,c [nJ]) is defined by Eq.(3.14). Using the b~r,c [nJ] segments of all the neighbouring antennas on a given row (r) an intermediate column vector γr [n] is initialized as T ~ ~T T~ γ~r [n] = [bT r,l′ −CT [nJ], . . . , br,l′ [nJ], . . . , br,l′ +CT [nJ]] (3.17) where CT is the number of accepted columns, l′ is the closest antenna to the target on row r and γ~r [n] is the contribution of the rth row to the input vector. This step effectively limits the number of columns accepted into the neighbourhood for a given row. An exception is made in the creation of γ~r [n] when r = k. In this case γ~k [n] is defined by T ~ T~ T~ T~ γ~k [n] = [bT k,l−CT [nJ], . . . , bk,l−1 [nJ], bk,l+1 [nJ], . . . , bk,l+CT [nJ]] . (3.18) This exception ensures that the signal recorded at the target antenna is not included ~ [n]) at time step in the input vector. The row-column based input vector (bRCk,l n for target signal xl,k is generated by concatenating the column vectors (γ~r [n]) as follows ~ [n] = [γ T~ [n], . . . , γ~T [n], . . . , γ T~ [n]]T bRCk,l k−RT k k+RT (3.19) where RT guides the number of rows accepted into the neighbourhood. For all the trials RT is set to 2, resulting in 5 row neighbourhoods, while the columns are limited by setting CT to 2, 4, and 8. The 5-row 4-column (5r4c) neighbourhood accepts all antennas that are illuminating nearly a quarter of the breast, while the 5-row 8-column (5r8c) case admits sensors that are irradiating half of the breast. Fig. 3.15 demonstrates a 5-row, 2 column (5r2c) neighbourhood as two antennas on either side of the target are admitted, along with the five closest antennas in the two rows above and beneath the 67 target antenna level. Once again no exceptions are made when a row above or below does not exist. Table 3.5 outlines the number of antennas used at any given row for a particular neighbourhood, while Table 3.6 summarizes the performance of the row/column neighbourhoods. Table 3.5: Neighbourhood Size Row 5-row 2-column Neighbourhood 5-row 4-column 5-row 8-column 5-row for Each Row Using Row-Column Selection 1 2 3 4 5 6 7 8 14 19 24 24 24 24 24 19 26 35 44 44 44 44 44 35 50 67 84 84 84 84 84 67 64 84 116 123 135 147 159 127 9 14 26 50 95 Table 3.6: Image Quality Using Row-Column Based Neighbourhoods. # Tumour SCR (dB) Localization Error (mm) Neighbourhood Responses Tumour Skin Tumour Skin Preserved Estimate Subtracted Estimate Subtracted Data Data 5-row 2-column 140 4.17 1.84 1 3.16 5-row 4-column 168 4.90 2.59 1 2.24 5-row 8-column 175 4.87 2.35 1 2.24 5-row 181 4.98 1.95 1 2.24 Examining Table 3.6 reveals that the 5-row neighbourhood preserves the most tumours, while the 5r2c scenario preserves the least number of tumours. This result is expected as the amount of data used to create the input vector has decreased dramatically in the 5r2c case, which would hinder the performance of the filter for its derivation is based on statistical means. The SCR of the tumour estimate images also supports the 5-row neighbourhood as the best candidate; however, it is also interesting to note that the 5r4c neighbourhood has a SCR comparable to the 5-row case. When looking at the SCR of the skin-subtracted images, the best value occurs 68 in the 5r4c scenario, followed by the 5r8c scenario. It is also interesting to note that the 5r8c scenario performs reasonably well over all three measures. In attempt to improve on the performance of the neighbourhoods, a variation on the row-column scheme is examined. Specifically data collected from the corners of the neighbourhood are removed from the input vector. As this modification only has an impact on the furthermost accepted rows, their corresponding γ~r′ [n] vectors are redefined as T ~ ~ ′~ T~ T γk−RT [n] = [bT k−RT,l′ −(CT−1) [nJ], . . . , bk−RT,l′ [nJ], . . . , bk−RT,l′ +(CT−1) [nJ]] (3.20) and T ~ ~ T~ ′~ T γk+RT [n] = [bT k+RT,l′ −(CT−1) [nJ], . . . , bk+RT,l′ [nJ], . . . , bk+RT,l′ +(CT−1) [nJ]] (3.21) The γ~r [n] vectors associated with all other rows (including the target antenna row) are defined using the same method described above. Consequently, the input vector for target signal xk,l maybe defined as follows, T ~T ′T~ ′T ~ ~ [n] = [γ ′T~ [n], γ ′T ~ b′RCk,l k−RT k−(RT−1) [n], . . . , γk [n], . . . , γk+(RT−1) [n], γk+RT [n]] . (3.22) This modification is motivated by the antenna stagger, as demonstrated by Fig. 3.15. In this figure, a 5r2c neighbourhood is highlighted in light purple about a red target antenna at the center. Looking at the antennas in the corners, it is seen that the two antennas highlighted in dark purple are much further away from the target antenna than any other neighbour. Table 3.7 presents the image evaluation results and the number of tumours detected when neighbourhoods without corners 69 are tested. For simplicity any row/column scenarios mentioned beyond this point refers to this selection schemes that does not include corner antennas. Figure 3.15: 5r2c Antenna Distribution: The target antenna is a red ’x’, while the neighbouring antennas are the black ’x’s highlighted in purple. The two corner antennas marked in dark purple, are significantly further away from the target. The above antenna distribution is a planer representation of the actual scan pattern where the antennas are placed to reflect the sensor stager between adjacent rows. Table 3.7: Image Quality Using Modified Row-Column Based Schemes. # Tumour SCR (dB) Localization Error (mm) Neighbourhood Responses Tumour Skin Tumour Skin Preserved Estimate Subtracted Estimate Subtracted Data Data 5-row 2-column 131 4.02 0.48 1 3.32 5-row 4-column 161 5.10 2.07 1 2.24 5-row 8-column 180 4.96 2.64 1 2.24 5-row 181 4.98 1.95 1 2.24 The results for the 5r8c case are very promising: the number of preserved tumour responses and the SCR for the tumour estimate data are almost identical to the 5row case, and the best SCR for the skin-subtracted images is achieved. The margin between the skin-subtracted SCR of the 5r8c and the 5-row neighbourhoods suggest that 5r8c case is the most desirable candidate. This conclusion is also supported by the tumour-to-skin ratios summarized in Fig. 3.16. The ratios of the 2- and 4-column 70 cases are typically smaller than those recorded for the row neighbourhoods while the values for the 5r8c scenario are comparable to the row cases. This figure also exhibits the same pattern as the row neighbourhoods where row 1 has a relatively high value and row 2 drops significantly. As previously outlined, this shape is a direct reflection of the variation in the curvature of the breast, and the number of antennas included in the input vector. Tumour−to−Skin Ratio 30 original 2 col 4 col 8 col Skin Subtracted Ratio (dB) 55 50 20 10 45 0 40 −10 35 −20 30 −30 25 −40 20 −50 15 −60 10 0 2 4 Row 6 8 Original Ratio (dB) 60 −70 10 Figure 3.16: Tumour-to-Skin Ratios Using Row-Column Neighbourhoods: The average and standard deviation of the tumour-to-skin ratios before and after skin subtraction are plotted for every row. The averages are represented by circles or squares while the standard deviation is demonstrated using errorbars. Fig. 3.17 presents the percentage of preserved tumour responses at every row and overall for row/column based neighbourhoods and the 5-row case. The 5r2c scenario is not included for its tumour response preservation is poor as a result of having limited data to create the input vector. The 8 column neighbourhood has the best 71 performance at rows 2, 3, and 9 where the breast exhibits the largest variations in curvature. This case also preserves more tumour responses than the 5-row scenario in rows 6 and 7, which is highly desirable as the strongest tumour responses are found Tumour Response Preservation 4 col. 8 col. 5 rows 100 80 90 70 80 60 70 50 60 50 40 40 30 30 20 20 10 10 0 1 2 3 4 5 Row 6 7 8 9 Total 0 Overall Percentage of Preserved Responses Percentage of Preserved Responses per Row in these rows. Overall, the 5r8c neighbourhood appears to be the best candidate. Figure 3.17: Preserved Tumour Response Percentage Using Row-Column Schemes 3.3.3 Beamwidth Based Neighbourhoods Up to this point, the neighbourhood selection techniques have been strongly dependent on the antenna distribution within the scan pattern. However, the distance between two rows or columns may vary for different scans. It should also be noted that the data employed for these trials are not acquired from an antenna distribution wherein the receivers are evenly spaced. As a result a neighbourhood selection technique that does not depend on the antenna distribution is required. The method 72 Figure 3.18: Selection of Patch Neighbourhood: Antennas that have overlapping half-energy beamwidths are admitted into the neighbourhood. The resulting neighbours for the target antenna represented by the red x, are the antennas illustrated in black while the antennas shown in grey are rejected. proposed in this section establishes the neighbourhoods based on the beamwidth of the radiating antenna. This beamwidth is approximated using a rectangular patch defined by the 3 dB cutoff of the energy radiated 2 cm away from the aperture in the planes parallel and perpendicular to the antenna. As demonstrated by Fig. 3.18, antennas that have overlapping patches are included in the neighbourhood. This selection approach produces what is called patch neighbourhoods. The patch neighbourhood is tested along with two modifications to the idea. The first involves only considering antennas within the 5-row neighbourhood. This selection scheme is called patch row limit (prl) and is developed to account for the variations of the curvature of the breast without having to explicitly calculate it. The second modification is called 2 patch row limit (2prl) and it multiplies the width of the patch by two in order to account for the limited illumination that occurs outside of the antenna patch. This enhancement also increases the number of antennas allowed into the neighbourhood. A mathematical description of the selection process of the 73 three proposed patch based neighbourhoods follows. For a selected target antenna τ , located at xτ , yτ , zτ , all other antennas in the scan are considered as potential candidates for inclusion into the neighbourhood. Two conditions must hold in order for any given antenna υ, to be included in the neighbourhood Npatch,τ . This is expressed mathematically as υ ∈ βk ≥ Npatch,τ if f p (xτ − xυ )2 + (yτ − yυ )2 (3.23) and β⊥ ≥ p (zτ − zυ )2 where βk and β⊥ are the parallel and perpendicular dimensions of the half energy beamwidth of the antenna calculated 2 cm away from the aperature. For the patch and prl neighbourhood βk is set to 4.4 cm, while in the 2prl this threshold is doubled to a value of 8.8 cm. The β⊥ is set to 2.5 cm for the prl and the 2prl neighbourhoods, whereas a value of 3.5 cm is used for the patch neighbourhood. Once the members of the neighbourhood have been established, the contribution of the signal xυ recorded by neighbouring member υ to the input vector at time step n is defined by b~υ [n] = [xυ [n − J], . . . , xυ [n], . . . , xυ [n + J]]T (3.24) where 2J +1 is the width of the small sliding window used to create the input vector. Consequently, the patch based input vector (b~Pτ [n]) at time step n for target antenna τ may be defined as T ~ T ~T b~Pτ [n] = [b~T 1 [n], . . . , bυ [n], . . . , bυNS [n]] (3.25) 74 where N S is the number of antennas included in the patch neighbourhood (Npatch,τ ). Table 3.8 summarizes the performance measure results when applied to the patch based neighbourhoods. Table 3.8: Image Quality Using Beamwidth Based Neighbourhoods. # Tumour SCR (dB) Localization Error (mm) Neighbourhood Responses Tumour Skin Tumour Skin Preserved Estimate Subtracted Estimate Subtracted Data Data patch 161 5.78 1.93 1.41 2.24 patch row lim. 152 4.90 2.17 1 2.24 2 patch row lim. 184 4.70 2.08 1 2.24 5-row 8-column 180 4.96 2.64 1 2.24 5-row 181 4.98 1.95 1 2.24 The poor tumour response preservation rate of the patch and the prl neighbourhoods further reinforces the importance of having sufficient data in the input vector, as the sizes of these neighbourhoods are significantly smaller than the other three presented schemes. The 2prl case preserves the most tumour responses while having a reasonable SCR for the skin-subtracted images and a slightly lower SCR for the tumour estimate images than the 8 column and 5-row cases. Investigation of the resulting tumour-to-skin ratios reveal no striking differences between the patch based neighbourhoods. Furthermore, the ratios lie within the ranges previously observed for the 5-row and 5r8c neighbourhoods. Fig. 3.19 demonstrates the percentages of preserved tumour responses for the patch, 2prl, and 5r8c neighbourhoods. The superior performance exhibited by the 2prl scenario in row 1, is caused by an increased number of antennas forming the input vector. The hemispherical nature of the scan is directly responsible for the larger neighbourhood, as the antennas at the bottom of the array tend to be closer 75 which leads to more instances of beamwidth overlap. Consequently, the input vectors created for row 1 are effectively the same as the ones produced by the 5-row neighbourhood. The 2prl neighbourhood is also very effective in rows 5 to 8, as it preserved all the tumour responses in row 7 and only one less tumour than the other two cases at row 6. A strong preservation rate in this region is highly desirable as these rows are in close vicinity to the tumour. Although the 8 column neighbourhood creates marginally better images, the 2prl scenario is deemed more robust, as it preserves more tumour responses and is an adaptable approach that makes no prior Tumour Response Preservation 100 80 patch 2prl 8 col. 90 80 70 60 70 50 60 50 40 40 30 30 20 20 10 10 0 1 2 3 4 5 Row 6 7 8 9 Total 0 Overall Percentage of Preserved Responses Percentage of Preserved Responses per Row assumption on the antenna spacing in the scan pattern. Figure 3.19: Percentage of Preserved Tumour Responses With Patch Schemes 76 3.3.4 Antenna Distributions & Cross Correlation Neighbourhoods As a final check to ensure that the data included in the input vector is appropriate, a normalized cross correlation is applied to the neighbours selected by the 2prl scheme. The motivation behind this modification is to discount or remove erroneous signals that are recorded as a result of the antenna coupling with the breast, or any other unexpected artifact. Fig. 3.20(a) demonstrates an example of a signal recorded when the antenna has coupled with the breast, while Fig. 3.20(b) shows a typical response. 8000 6000 6000 4000 4000 2000 Voltage Voltage 2000 0 −2000 0 −2000 −4000 −4000 −6000 −6000 −8000 −10000 0 0.5 1 1.5 Time (s) (a) 2 2.5 3 −9 x 10 −8000 0 0.5 1 1.5 Time (s) 2 2.5 3 −9 x 10 (b) Figure 3.20: Comparison of Recorded Signals: The signals recorded at different locations in the scan may have different shapes. For example the signals seen in (a) is an erroneous signal that is caused by the antenna coupling with the breast, while (b) demonstrates a more common signal. The proposed selection process adds the stipulation that signals must have a correlation greater than a predetermined threshold to be included in the input vector. From a mathematical perspective, this neighbourhood is selected by simply updating 77 the conditions outline in Eq.(3.24) to υ ∈ Npatchx,τ if f p βk ≥ (xτ − xυ )2 + (yτ − yυ )2 p β⊥ ≥ (zτ − zυ )2 (3.26) ψ ≤ xcorr(xτ , xυ ) where ψ is the predetermined cross correlation threshold, while xτ and xυ are the signals recorded at the target antenna and the neighbouring antenna, respectively. Once the members of the neighbourhood have been established, the input vector is initialized using Eq.(3.24-3.25). The threshold (ψ) is established by examining the cross correlations between all antennas in the array, and the correlations between the target antenna and its neighbours as selected by the 2prl method. The best cross correlation achieved between two antennas within the scan is 99.57%, while the worst is 59.16%. When the correlations in the 2prl neighbourhoods are considered, the average maximum cross correlation is approximately the same as the highest correlation between any two signals, while the worst correlation improves to 82.31%. This indicates that neighbourhoods that are selected based on antenna proximity increase the correlation between the target and the input vector. The tested cross correlation thresholds are 97%, 95%, 90%, 85% and 80% and are specifically chosen to lie within the range of the average cross correlation values calculated for the 2prl scenario. To further refine the number of tests, the average neighbourhood sizes generated by the 5 proposed thresholds are tabulated in Table 3.9 with reference to the number of antennas in the 2prl. 78 Table 3.9: Average Neighbourhood Size Using Normalized Cross Correlation Normalized Cross Mean # of Antennas Mean % of Accepted Correlation in Neighbourhood Antennas Relative to 2prl 97% 40 58.37 95% 51 75.78 90% 61 90.46 85% 63 94.36 80% 64 95.85 Based on Table 3.9 it is concluded that 97% and the 85% should not be further evaluated. The 85% threshold is removed because its statistics do not vary greatly from those reported for the 80% case. The 97% threshold is also removed, as previous trials of neighbourhoods with 40 antennas or less have lead to tumour response preservation rates between 50-60% (Table 3.3 and Tables 3.5 - 3.7). Table 3.10 summarizes the results of the single signal analysis and the image evaluation tools when applied to data processed using the 2prl neighbourhood in conjunction with the three normalized cross-correlation thresholds of interest. Table 3.10: Image Quality Using Beamwidth & Cross Correlation Based Schemes. # Tumour SCR (dB) Localization Error (mm) Neighbourhood Responses Tumour Skin Tumour Skin Preserved Estimate Subtracted Estimate Subtracted Data Data 2prl 95 177 4.70 1.69 1.41 2.24 2prl 90 181 4.80 1.80 1 2.24 2prl 80 181 4.74 2.19 1 2.24 2prl 184 4.70 2.08 1 2.24 Of the 3 threshold based neighbourhoods presented in Table 3.10, it is clear that the 95% threshold does not perform as well as other cases. The 90% and 80% scenarios on the other hand only differ in terms of their SCR. Examination of both 79 the tumour-to-skin ratios and the percentage of preserved tumour responses per row give no further insight as to which neighbourhood is better. Consequently, the performance of both neighbourhoods is evaluated based on the results achieved on an atypical signal. Target Signal 1500 XCorr 81.9693 2000 1000 XCorr 91.225 1500 1500 1000 500 1000 −500 −1000 −1500 500 500 Voltage Voltage Voltage 0 0 −500 0 −500 −2000 −1000 −2500 −3500 −1000 −1500 −3000 0 0.5 1 1.5 Time (s) (a) 2 2.5 3 −9 x 10 −2000 0 0.5 1 1.5 Time (s) (b) 2 2.5 3 −9 x 10 −1500 0 0.5 1 1.5 Time (s) 2 2.5 3 −9 x 10 (c) Figure 3.21: Comparison of Potential Neighbours: The neighbourhood of the atypical target signal seen in (a) would include the antennas associated with the signals seen in (b) and (c) if the 2prl 80 scheme is employed; however, the 2prl 90 scheme would reject the antenna that recorded the signal seen in (b) from the neighbourhood. Fig. 3.21(a) is an example of a signal recorded at the top row of the scan where the breast meets the chest wall which causes an extreme curvature to the surface. The unexpected shape of this signal is caused by the incident wave interacting with two distinct surfaces simultaneously. Fig. 3.21(b) and 3.21(c) show an example of a signal that would be included in the neighbourhood if the 80% and 90% thresholds are applied, respectively. When the 80% threshold is in place, 36 antennas are used to create the input vector, while only 3 antennas are included when the 90% threshold is used. Although 36 antennas are used to estimate the skin reflection, this signal’s tumour response is still not successfully recovered. The 90% case does not attempt to estimate the skin response as only 3 neighbours are available and a preliminary threshold of 10 antennas has been placed on the amount of data required 80 to estimate the skin reflection. This preliminary threshold is set based on the smallest neighbourhood that is capable of preserving the tumour response. Consequently, the 90% threshold is deemed more appropriate, as the signal seen in Fig. 3.21(a) should not be processed for it does not have reliable data to estimate the skin response. Fig. 3.22 illustrates the tumour estimate and skin-subtracted images generated when 20 1 20 1 40 0.8 40 0.8 0.6 60 0.4 80 0.2 100 0 120 50 100 y (mm) (a) 150 200 x (mm) x (mm) the 2prl 90 neighbourhood is employed. 0.6 60 0.4 80 0.2 100 0 120 50 100 y (mm) 150 200 (b) Figure 3.22: Resulting Images Using the 2prl 90 Neigbourhood: The images generated when the 2prl 90 neighbourhood and the Modified RLS are applied to tumour estimate data (a) and skin-subtracted data (b). The x-y planes shown above contain the largest scatterer within the imaging volume. Overall, the 2prl neighbourhood preserves the most tumour responses and the 5r8c scenario generates the best SCR; however, the 2prl 90 neighbourhood appears to be the most promising candidate. This option is the most desirable as it is independent of the antenna spacing within the scan and intelligently removes anomalous signals from the input vector. The 2prl 90 neighbourhood is also capable of successfully detecting the tumour in both images at a reasonable SCR level, with a typical localization error. As a result the 2prl 90 neighbourhood is the data selection technique that is applied in the Neighbourhood Modified RLS algorithm. Fig. 3.23 demonstrates the first enhancement made to the Modified RLS implemented to suc- 81 cessfully skin-subtract data acquired from a simple three dimensional realistic breast model. The element highlighted in dark blue represents the update that ensued as a result of the findings of this section while the elements in light purple are investigated in sections 3.4-3.5. Figure 3.23: Flow Chart of Updated Modified RLS: The 2prl neighbourhood selection technique is incorporated into the skin subtraction technique. The steps highlighted in light purple are examined further in this chapter while the step with the dark blue background has been updated based on the results of this section. 3.4 Filter Weight Estimation As outlined in Fig. 3.23 the filter weight estimation process encompasses several steps. In this section, techniques to adaptively identify the SDW are outlined, and the effect of varying this parameter is studied. This investigation is followed by the selection of a rank detection technique to further refine the calculated filter weights. 82 All of these trials are conducted on the simple breast model introduced in section 3.1. The 2prl 90 neighbourhood presented in section 3.3.4 is also employed. Just as in section 3.3, the width of the small sliding window (2J+1) is set to 7 time steps and the rank is established using the approach outlined in section 3.4.2. 3.4.1 Adaptive Skin-Dominant Window Selection In [16], a skin-dominant window selection technique is not described as this problem is trivial with 2-D breast models. The 3-D models present the challenge of exhibiting significantly more complex reflections as the secondary skin response is not seen in 2-D. As a result, the SDW can no longer be manually selected. There are two main criteria that the SDW identification technique must follow. First, the window should not include any part of the tumour response, as this will distort the tumour reflection. Second, the selection process must be adaptive. The need for an adaptive approach is illustrated by Fig. 3.24 for the durations of the skin-dominant portion of these two signals are different. The requirement of having an adaptive SDW selection technique is further reinforced by the fact that the tumour response of the signal seen in Fig. 3.24(a) starts at time step 1628 (0.78 ns) while the SDW of the signal illustrated in Fig. 3.24(b) does not end until time step 1805 (0.87 ns). Two skin-dominant window selection schemes that adhere to the above two criteria are developed and tested. 1. The blind selection method assumes no prior knowledge of the location of the skin or tumour, and relies entirely on the calibrated data to establish the window. This process starts by normalizing and squaring the data. The location of all the peaks and troughs throughout the squared signal are then determined 83 4 1 x 10 2000 1500 0.5 1000 Voltage Voltage 500 0 −0.5 0 −500 −1000 −1500 −1 −2000 −1.5 0 0.5 1 1.5 Time (s) 2 2.5 3 −2500 0 0.5 −9 x 10 (a) 1 1.5 Time (s) 2 2.5 3 −9 x 10 (b) Figure 3.24: Example of Two Signals That Require a Different SDW: The signal is illustrated with a solid black line while the SDW is marked with dashed lines. The width of the window in (a) is 1010 time steps (4.85 e-10 s) while the duration of the SDW in (b) is 1297 time steps (6.23e-10 s). with the use of derivatives. As the signals are discretized the derivatives are approximated using differences. So for target signal xτ , the approximation of the first yτ′ [n] and second yτ′′ [n] order derivatives of the squared signal at time step n is yτ′ [n] = x2τ [n + 1] − x2τ [n] ≈ dx2τ [n] dt (3.27) d2 x2τ [n] . dt2 (3.28) and yτ′′ [n] = yτ′ [n + 1] − yτ′ [n] ≈ Consequently, time step n is deemed to be at a peak or a trough of the signal, if the estimate of the first order derivative is approximately equal to zero. The peaks are delineated from the troughs using the approximation of the second 84 order derivative as follows xτ [n] ∈ peak if f (3.29) yτ′ [n] ≈ 0 yτ′′ [n] < 0 or xτ [n] ∈ trough if f (3.30) yτ′ [n] ≈ 0 yτ′′ [n] > 0 The algorithm then searches for at least 3 peaks before the maximum of the signal, and 9 peaks after the maximum, along with their corresponding troughs. If 3 or 9 peaks do not exist before and after the maximum, respectively, the search algorithm simply exits. To ensure the tumour response is not included in the window, significant peaks are selected using the following criteria for the dipole antenna xτ [np ] > 0.003 xτ [nM p ] (3.31) where xτ [np ] is the value of the signal at the time step associated with the peak of interest, and xτ [nM p ], is the value of the signal at the largest peak. As another precaution to ensure that the tumour response is not incorporated into the SDW, peaks are accepted only if they are adjacent to another significant peak with the exception of the first and last detected significant peaks. Once the desired peaks are identified, the front edge of the window is set to be the 85 trough before the first significant peak, while the second boundary is the trough that follows the last significant peak. 2. The ideal selection requires prior knowledge of the tumour response. In this case, the front edge of the window is selected in the same way as the blind selection, but the second edge is selected with reference to the tumour response. Using a similar technique to blind selection, the significant peaks of the tumour response are established and the troughs are identified. The second boundary is then set to the trough that comes before the first significant peak of the tumour response. This approach is designed on the premise that TSAR may be used as a complementary modality to another imaging technique; therefore, an estimate of the location may be available prior to processing the data. A preliminary version of this method is employed for the neighbourhood selection trials discussed in section 3.3. To assess the skin removal capabilities of the proposed SDW selection techniques, a performance measure called the tumour-to-signal ratio is adopted. Much like the tumour-to-skin ratio described in section 3.2, the p2p values of the skin-subtracted data and the known tumour response are compared. The major difference between both measures lies in the calculation of the p2p value for the skin-subtracted data. The p2p value for the tumour-to-skin ratio is established by only examining the data found within the SDW, while the tumour-to-signal ratio includes the entire duration of the signal. Ideally the tumour-to-signal ratio should be 0 dB. The results of this measure are also analyzed by contrasting its value with the original tumour-to-skin ratio. Fig. 3.25 summarizes the tumour-to-signal ratios calculated from the skin- 86 subtracted data generated using both the blind and ideal selection techniques. Tumour−to−Signal Ratio Original Tumour−to−Signal Ratio (dB) 0 Blind Window 40 Ideal Window 30 −5 20 −10 10 −15 0 −20 −10 −25 −20 −30 −30 −35 −40 −40 −50 −45 −60 −50 0 2 4 Row 6 8 Tumour−to−Skin Ratio (dB) 5 −70 10 Figure 3.25: Tumour-to-Signal Ratio of Each Row When Blind & Ideal SDW Selection is Applied: The average and standard deviation of the tumour-to-signal ratios before and after skin subtraction are plotted for every row. The means are displayed using circles or squares while the standard deviations are shown using errorbars. As observed in previous plots, the averaged results exhibit a parabolic shape related to the variation in breast curvature (row 1 is an outlier as before). Both windows exhibit a pronounced improvement from the original average tumour-toskin ratio; however, the average tumour-to-signal ratios for the ideal windows are consistently better than the ones computed for the blind windows. This result is expected as an ideal window is designed to reduce the secondary skin reflections, while maintaining the tumour response. It should also be noted that successful tumour detection without clutter reduction is only achieved with the ideal window. Although the ideal SDW performs considerably better than the blind SDW, there 87 is still a noteworthy difference between the skin-subtracted signal generated using the ideal case and the known tumour responses. This indicates that there are still residual secondary skin reflections. The average difference between the size of the blind and ideal windows is 418 time steps, while the largest and smallest differences are 870 and 210 time steps, respectively. The mean width of the blind window is 904 time steps. 3.4.2 Adaptive Rank Selection Establishing the low rank approximation of the autocorrelation matrix is a critical step in the filter weight estimation process. Underestimating the rank leads to large residual skin responses, while overestimating introduces noise into the system. Several rank estimation algorithms in the context of linear least-squares problems have been proposed. A commonly used and effective technique for solving ill-conditioned linear systems is the generalized cross validation (GCV) method presented in [50]. The GCV operates on the premise that allowing a bias in the estimation process can greatly reduce the variance in the filter coefficient estimates. First, the eigenvalue decomposition of the autocorrelation matrix R is taken: R = UΛU−1 (3.32) where Λ is a diagonal matrix of eigenvalues and U is a matrix of corresponding eigenvectors. The aim is to establish the number of significant eigenvalues (k). The GCV approximates the rank to be the value that minimizes the following (I − RR−1 (k))~p2 g(k) = [trace(I − RR−1 )] (3.33) 88 where I is the identity matrix, R−1 (k) is the inverse of the low rank approximation of R for any given k, and ~p is the cross correlation between the input and the desired response. The low rank approximation is thereby determined by iteratively varying k until the minimum of g is found. Previously, the rank of the autocorrelation matrix was set to a static value of 15 [16]. Integrating the GCV method into the skin subtraction algorithm leads to fluctuations in the estimated rank. These deviations are quantified by calculating the minimum, maximum, mean and standard deviation of the estimates over all signals. Table 3.11 compares the above statistical measures, when the ideal and blind windows are used to select the SDW. Window Blind Ideal Table 3.11: Statistics of Minimum Maximum 10 22 19 35 Estimated Ranks Mean Standard Deviation 15.64 2.15 27.24 3.21 The mean rank generated using the blind window is significantly smaller than the average rank of the ideal window. The resulting standard deviation of the ideal window is also larger. This trend indicates that the duration of the SDW has a considerable bearing on the computed rank. The larger mean value observed in the ideal case is likely due to the increased amount of information that is incorporated into the autocorrelation matrix, as a result of using a longer SDW. It is also interesting to note that the initial rank estimate of 15 used by [16] agrees with the average rank computed using the blind window. As the increase in computational time caused by the adaptive selection of the rank is marginal, the GCV is incorporated into the enhanced skin subtraction algorithm. Fig. 3.26 illustrates the updated skin subtrac- 89 tion technique. The feedback step highlighted in light purple is examined further in section 3.5. Figure 3.26: Flow Chart of the Enhanced Modified RLS: The adaptive SDW selection and rank determination is incorporated into the skin subtraction approach under development. The steps highlighted in light purple are examined further in this chapter while the steps with the dark blue background have been updated based on the results of the previous two sections. 3.5 Feedback Filtering The feedback filter aims to mitigate distortion to the tumour response as a result of skin subtracting the data. To further understand this distortion the filtering process is modeled mathematically and the value of implementing the feedback filter is evaluated. The differences in the results of the signals processed with and without the feedback are also compared using the performance measures outlined in section 3.2. 90 3.5.1 Mathematical Derivation of Feedback To mathematically describe the effect of the feedback filter, the signal of antenna 1 is modeled as: x1 = s1 + i1 + t1 (3.34) where x1 is a calibrated target signal,s1 is the skin response, i1 is the response from the interior of the breast and t1 is the tumour response. In order to establish the contribution of the feedback, the filtered signal at antenna 1 is analyzed. This signal is represented by r1 [n] = x1 [n] − ~qT ~b2N [n] (3.35) where r1 [n] is the filtered signal at time step n prior to feedback, ~qT is the vector of filter weights and ~b2N [n] is the input vector. As described in section 2.2.2, the input vector is comprised of concatenated segments of data that originate from the target signal’s neighbours. Consequently, the input vector can be modeled as follows b2N [n] = s2N [n] + i2N [n] + t2N [n] (3.36) where s2N [n],i2N [n],t2N [n] are the skin, breast interior, and tumour response contributions to the input from the target’s neighbourhood. Combining Eq.(3.36) with the equation for the filtered response at antenna 1, we obtain r1 [n] = s1 [n] − ~qT s2N [n] + i1 [n] − ~qT~i2N [n] + t1 [n] − ~qT ~t2N [n] = s′1 [n] + i′1 [n] + t1 [n] − ~qT ~t2N [n]. (3.37) where s1 [n],i1 [n],t1 [n] are the skin, breast interior, and tumour responses of the target antenna signal while s′1 [n] and i′1 [n] are residual skin and interior breast responses. 91 Because the filter weights are designed to eliminate the skin response, Bond et al. assumed the s′1 [n] term to be approximately equal to zero while i′1 [n] was not accounted for. As Eq.(3.37) is the output of the skin subtraction without feedback, the goal of the feedback filter is to compensate for the last term on the right hand side. Consequently, the feedback attempts to alleviate this distortion by using the output of Eq.(3.37) as an estimate for the tumour responses. The ~qT ~t2N [n] term can thereby be eliminated by applying the same filter weights to the resulting skin-subtracted signals, and adding it back into the signal as follows, r˜1 [n] = r1 [n] + ~qT ~r2N [n] where • r˜1 [n] is the filtered signal with feedback, • r1 [n] is the tumour response estimated from Eq.(3.37), • and ~r2N [n] is an input vector created out of the tumour response estimated from Eq.(3.37). To evaluate the benefit of adding the feedback, the input vector ~r2N [n] is modeled using the ideal signal model as follows: r2N [n] = s′2N [n] + i′2N [n] + t2N [n] − t2Nfilt [n] (3.38) where s′2N [n], and i′2N [n] are the residual skin and breast interior response, while t2Nfilt [n] is the vector of ~qT ~t2N [n] contributions. Substituting Eq.(3.38) and Eq.(3.37) into the feedback equation it is seen that feedback filter can be represented as r˜1 [n] = s′1 [n] + i′1 [n] + t1 [n] + ~qT s′2N [n] + ~qT i′2N [n] − ~qT ~t2Nfilt [n]. (3.39) As expected the ~qT ~t2N [n] is eliminated; however, comparing Eq.(3.39) to Eq.(3.37) it is seen that Eq.(3.39) has one additional term ~qT i′2N [n]. Consequently, the feedback 92 Figure 3.27: Flow Chart of Feedback Filter is only effective if ~qT i′2N [n] − ~qT ~t2Nfilt [n] ≤ ~qT ~t2N [n]. Fig. 3.27 demonstrates a flow chart of this process. To evaluate the benefit of the feedback filter, signals recorded from the simple breast model introduced in section 3.1 are processed with and without the feedback. The 2prl neighbourhood selection technique is used in conjunction with a 7 time step small sliding window (2J+1) and the rank detection algorithm described in section 3.4.2. As the above inequality is highly dependent on the residual of the interior response, the SDW is set using the ideal selection technique to establish the baseline performance. 3.5.2 Feedback Results The single signal performance of the feedback filter is first evaluated with reference to the filter without feedback. The data processed with feedback detected 118 tumour responses, while 181 are preserved without feedback. This large discrepancy is mostly due to the significant difference in the energy ratios calculated for the two trials. The feedback case produced 54 signals with an energy ratio greater than 100%, while this occurred only twice without feedback. The difference between the two ratios is likely due to the feedback adding extra unwanted energy to the signal, indicating that ~qT i′2N [n] − ~qT ~t2Nfilt [n] is in fact larger than ~qT ~t2N [n]. This hypothesis is also 93 supported by the resulting tumour-to-skin ratios displayed in Fig. 3.28. Tumour−to−Skin Ratio 30 50 20 40 10 30 0 20 −10 10 −20 original feedback no feedback 0 −10 −30 −40 −20 −50 −30 −60 −40 0 2 4 Row 6 8 Original Ratio (dB) Skin Subtracted Ratio (dB) 60 −70 10 Figure 3.28: Tumour-to-Skin Ratio With & Without Feedback: The average and standard deviation of the tumour-to-skin ratios before and after skin subtraction are plotted for every row. The averages are represented by circles or squares while the standard deviation is demonstrated using errorbars. The average tumour-to-skin ratios after skin subtraction for the data generated with feedback are all smaller than the ratios calculated without feedback. This suggests that the average p2p value of the skin-subtracted signal within the SDW has increased. This effect is not desirable as the goal is to eliminate the skin response. Intuitively, the feedback effectively adds a scaled version of the residual skin response at all of the neighbouring antennas back into the signal. Overall, it appears that this feedback is not beneficial when applied to 3-D models. Alternate methods may be pursued in order to estimate appropriate filter weights that alleviate the distortion to the tumour response; however, this is a tumour estimate problem. Consequently, 94 feedback is not incorporated in the Neighbourhood Modified RLS algorithm. 3.6 Neighbourhood Modified RLS Overview The overall goal of the Neighbourhood Modified RLS is to estimate the skin response of the target antenna as a filtered combination of the skin reflections measured at neighbouring antennas. Using basic adaptive filtering principles the filter weights are computed such that the mean square error between the target signal and the estimate is minimized over the skin dominant portion of the signal. As only a single set of filter weights are computed for a given target antenna, the statistics of the input vector are time averaged over the SDW to calculate the weights. The data incorporated into the input vector from the SDW for a given time step is selected using a small window. This small window slides across the entire SDW and a segment of every neighbouring signal that is centered about the time step of interest is incorporated into the input vector. Consequently, the skin response at a given time step is approximated as an averaged weighted combination of the data observed within the SDW at the neighbouring antennas. Once the skin response is estimated it is subtracted from the target signal, and this process is repeated for every sensor in the antenna array. More specifically for a given target antenna τ , the input vector at time step n is defined using concatenated segments of the signals recorded at its neighbouring antennas. The neighbourhoods are selected based on antenna proximity and the cross correlation between the signal recorded at the target (xτ ) and the signal acquired at a candidate neighbour (xυ ). The proximity of the sensors is gauged using the half-energy beamwidth of the antennas. The beamwidth is approximated using a 95 rectangular patch and the dimensions of the patch are dictated by the half-energy beamwidth of the antenna computed 2 cm away from the aperature. Given that target antenna τ is located at xτ , yτ , zτ , any given antenna υ is a member of the target antenna’s neighbourhood (Npatchx,τ ) if υ ∈ Npatchx,τ if f p βk ≥ (xτ − xυ )2 + (yτ − yυ )2 p β⊥ ≥ (zτ − zυ )2 (3.40) ψ ≤ xcorr(xτ , xυ ) where βk and β⊥ are the dimensions of the patch, xυ , yυ , zυ , are the coordinates of a potential neighbouring antenna, and ψ is the predetermined cross correlation threshold. The contribution of a neighbouring antenna υ to the input vector at time step n is b~υ [n] = [xυ [n − J], . . . , xυ [n], . . . , xυ [n + J]]T (3.41) where 2J +1 is the width of the small sliding window used to create the input vector. ~ [n]) at time step n for target antenna τ is the Consequently, the input vector (bPxτ following column vector ~ [n] = [b~T [n], . . . , b~T [n], . . . , b~T [n]]T bPxτ 1 υ NS (3.42) where N S is the number of antennas accepted into the neighbourhood. The input vector at the next time step is established by simply shifting the small sliding window over by one and following the same procedure. Once the input vector is established, the filter weights are computed such that 96 the following cost function is satisfied ~q = argmin q noX +m−1 n=no 2 T~ xτ [n] − ~q bP xτ [n] (3.43) where, ~qT is a column vector of filter weight, while no and no +m−1 are the first and second boundaries of the SDW. Two methods for SDW selection are proposed. Both methods establish the boundaries of the window such that the significant peaks and troughs of the skin response are included within the SDW. The significant peaks and troughs are established by looking at the first and second order derivatives of the signal of interest squared. The main difference between both methods is the signal of interest used. The first method called blind selection requires no prior knowledge of the tumour location and response; consequently, the signal of interest is set to be the calibrated data of the target antenna. The second method named ideal selection uses both the calibrated data and the tumour response of the target antenna to establish the boundaries. In order to describe both techniques using general terms the signal of interest is denoted by vτ . As the signals are discretized the derivatives are approximated using differences. Consequently, the first (yτ′ [n]) and second (yτ′′ [n]) order derivatives of the squared target signal of interest (vτ ), at time step n is yτ′ [n] = vτ2 [n + 1] − vτ2 [n] ≈ dvτ2 [n] dt (3.44) d2 vτ2 [n] . dt2 (3.45) and yτ′′ [n] = yτ′ [n + 1] − yτ′ [n] ≈ The portion of the signal at time step n is deemed to be at a peak or a trough, if the estimate of the first order derivative is approximately equal to zero. The peaks are 97 delineated from the troughs using the approximation of the second order derivatives as follows vτ [n] ∈ peak if f (3.46) yτ′ [n] ≈ 0 yτ′′ [n] < 0 or vτ [n] ∈ trough if f (3.47) yτ′ [n] ≈ 0 yτ′′ [n] > 0 Once the majority of the peaks and troughs that are centered about the maximum of the signal are identified, the significant peaks are established using two criterions. First, the amplitude of the peaks are tested as follows peak ∈ SigP eak if vτ [np ] > χ vτ [nM p ] (3.48) (3.49) where vτ [np ] is the value of the signal at the time step associated with the peak of interest, vτ [nM p ], is the value of the signal at the largest peak and χ is a threshold set based on the characteristics of the antenna used. Next, significant peaks are only accepted if they are adjacent to another significant peak with the exception of the first and last detected significant peak. Both stipulations attempt to ensure that the detected significant peaks are in fact associated with the response we are trying to isolate. In the case where the signal is set to be the calibrated data the targeted 98 response is the skin reflection, while the tumour reflection is isolated when vτ is the tumour response. As the blind SDW selection relies solely on the calibrated data, the edges of this window are selected based on the location of the peaks and the troughs of the squared signal. The first edge of the blind SDW (no ) is set to be the time step associated with the trough that precedes the first significant peak. The following boundary (no + m − 1) is the trough that follows that last significant peak. The ideal selection uses both the calibrated data and the tumour response to identify the SDW. Consequently, the significant peaks and their corresponding troughs are identified for the squared calibrated data, and the squared tumour response. The first boundary of this window (no ) is set the same way as the blind selection, while the other edge (no + m − 1) is located at the trough that precedes the first significant peak of the tumour response. The normal equation R~q = ~p (3.50) is used to solve the minimization problem, seen in Eq.(3.43), where R is the autocorrelation matrix of the input vector, and ~p is the cross correlation vector between the target signal and the input vector. These values are computed using a time average over the SDW as follows n +m−1 1 oX ~ R= bτ N [n]~bTτN [n] m n=n (3.51) o ~p = n +m−1 1 oX ~ bτ N [n]xτ [n]. m n=n (3.52) o In order to calculate the filter weights, the matrix R must be inverted. This is a 99 challenge, as this matrix is ill-conditioned as a result of the high degree of similarity between the neighbouring signals over the SDW. Consequently, the autocorrelation matrix is replaced with a low rank approximation. The GCV method is used to estimate the rank of R. This technique starts by applying eigenvalue decomposition to the autocorrelation matrix as shown, R = UΛU−1 (3.53) where Λ is a diagonal matrix of eigenvalues and U is a matrix of corresponding eigenvectors. The GCV approximates the rank to be the value that minimizes the following function (I − RR−1 (k))~p2 g(k) = [trace(I − RR−1 )] (3.54) where I is the identity matrix and R−1 (k) is the inverse of the low rank approximation of R for any given k. Consequently, the low rank approximation is determined by iteratively varying k until the minimum of g is found. The filter weights are thereby computed as ~qT = R−1 p GCV ~ (3.55) where, R−1 GCV is the inverse of the low rank approximation of R as calculated by the GCV method. The filter weights are then applied to the entire duration of the target signal, using the same procedure as described above to establish the input vector. The final skin-subtracted signal using the Neighbourhood Modified RLS method is represented by ~ [n]. rNModRLS,τ [n] = xτ [n] − ~qT bPxτ (3.56) 100 3.7 Summary This chapter outlines the necessary modifications to the Modified RLS in order to successfully skin subtract data recorded from a simple 3-D breast model. First, a neighbourhood identification technique must be applied to create an appropriate input vector. Several neighbourhood identification approaches are tested, and the most robust candidate selects the neighbours based on the antenna beamwidth and spatial distribution, along with the normalized cross correlation. Next, two adaptive skin-dominant window selection techniques are discussed, wherein one technique requires prior knowledge of the tumour location, and the other is completely independent. Both methods achieve the goal of significantly reducing the skin response. The technique that requires prior knowledge also has the positive side effect of producing signals that can generate images of the tumour. This is followed by the successful incorporation of an adaptive rank determination algorithm into the Modified RLS. Finally the effect of including feedback into the filter is examined and ultimately found to be detrimental. All of these enhancements combined form what is now called the Neighbourhood Modified RLS which is represented by the flow chart Fig. 3.29. 101 Figure 3.29: Flow Chart of the Neighbourhood Modified RLS: The steps highlighted in dark blue are the elements of the Modified RLS algorithm that have been examined and updated in order to create the Neighbourhood Modified RLS, while the grayed out step has been removed entirely from the algorithm. 102 Chapter 4 Performance of Neighbourhood Modified RLS To further assess the capabilities of the Neighbourhood Modified RLS developed in Chapter 3, the algorithm is applied to signals collected from simulations of a variety of numerical breast models. The first set of tests examines the response of this skin subtraction technique to the inclusion of glandular tissue in the numerical breast model. The next trial evaluates the robustness of this approach to models of different shapes. The last experiment investigates the application of the algorithm to data acquired using a different antenna. The chapter is concluded with a comparison of the results achieved throughout all the trials. 4.1 Response to Varying Glandular Distributions The Neighbourhood Modified RLS is tested on breast models that include glandular tissue to establish the feasibility of applying this algorithm in increasingly realistic settings. The reduction of the skin response and the preservation of the tumour response are evaluated as the presence of glandular tissue beneath the skin changes the reflected signals. 4.1.1 Models Four MR-based numerical breast models with varying degrees of complexity are tested. The first baseline model is the simple breast model introduced in Chapter 103 3, which is comprised of a 2-mm thick skin layer that bounds fatty tissue and a spherical tumour with a diameter of 6-mm. The next three models have the same shape but include different dielectric property distributions representing glandular tissues. The varying property distributions are established by applying the dielectric mapping approaches outlined in section 2.1.1 as follows: 1. Simple glandular : MR pixel intensities greater than a particular threshold are set to relative permittivity of 16 and a conductivity of 1 S/m while pixels below this value are mapped to fatty tissue. 2. Bimodal : Pixel intensities above a given threshold are linearly mapped to relative permittivity and conductivity ranges of 15.2 to 27.5 and 1.7 to 3.0 S/m, respectively. Pixel values below the threshold are assumed to be fat. 3. Piecewise-linear II : MR pixel intensities are linearly mapped to a relative permittivity and conductivity ranges of 1 to 36 and 1 to 4 S/m, respectively. An artificial 2-mm thick skin layer is applied to the simple glandular and bimodal models. The relative permittivity of the artificial skin layer is set to 36 while the conductivity is 4 S/m. The piecewise-linear II model relies solely on linear mapping to establish the shape and properties of both the skin and glandular tissue. Consequently, the skin layer of this model has varying thickness, and dielectric properties. A 6 mm spherical tumour is included in all models at the same location and is assigned a relative permittivity of 50 and a conductivity of 4 S/m. The relative permittivity of the fatty tissue is set to 9 while the conductivity is 0.4 S/m. Dispersion is not accounted for in these models. Fig. 4.1 illustrates the relative permittivity of all four models along the z-plane that cuts through the center of the tumour. 104 Simple Breast Simple Glandular 50 45 20 50 45 20 40 60 30 25 80 20 100 5 100 y (mm) 150 30 25 80 20 15 120 10 50 35 60 100 15 120 140 40 40 35 x (mm) x (mm) 40 140 200 10 5 50 (a) 100 y (mm) 150 200 (b) Bimodal Piecewise Linear II 50 45 20 50 45 20 40 40 60 30 25 80 20 100 15 120 140 40 35 10 5 50 100 y (mm) (c) 150 200 x (mm) x (mm) 40 35 30 60 25 80 20 100 15 120 10 140 5 50 100 y (mm) 150 200 (d) Figure 4.1: Relative Permittivity Maps of the Same Breast Model with Varying Tissue Distributions: The simplest model is called Simple Breast (a), while the most complicated model is referred to as Piecewise Linear II (d). The two intermediate models are named Simple Glandular (b) and Bimodal (c), where the Bimodal model is more complex of the two. 105 As the dielectric contrast between the fat and the glandular tissue increases, so does the complexity of the signals. An increased dielectric contrast results in larger reflections originating from the glandular structures which interfere with the tumour response. Consequently, the piecewise linear II model generates the most challenging signals while the bimodal model has the next most complex reflections. The data is acquired by conducting FDTD simulations of the illumination of each breast model using the same antenna, scan pattern and excitation described in section 3.1. 4.1.2 Results Two major aspects are investigated when examining the application of the Neighbourhood Modified RLS to the four models described above. First, the algorithm’s ability to significantly reduce the skin reflections is demonstrated. Next, the preservation of the tumour responses is discussed. Both attributes are quantified by applying the performance measures introduced in section 3.2 to the skin-subtracted data. Fig. 4.2(a) demonstrates the original tumour-to-skin ratio of the four tested models, while Fig. 4.2(b) shows the resulting tumour-to-skin ratios using the blind window. As expected, the simple breast model has the highest average original tumourto-skin ratio while the piecewise linear II model has the lowest ratio. The skin reflections for all 4 models are successfully reduced as their mean tumour-to-skin ratios are all greater than zero; however, the extent of the skin response reduction depends on the complexity of the model. The algorithm is the least effective on the 7th row of the piecewise linear II model. This poor performance is caused by the significant number of signals in row 7 that originally have low tumour-to-skin ratios, as a result of their tumour responses exhibiting very small p2p values. 106 Original Tumour−to−Skin Ratio (dB) −50 Original Ratio (dB) −60 −70 −80 −90 Simple Breast Simple Glandular Bimodal Piecewise Linear II −100 −110 0 2 4 Row 6 8 10 (a) Filtered Tumour−to−Skin Ratio (dB) 60 Skin Subtracted Ratio (dB) 50 40 30 20 10 Simple Breast Simple Glandular Bimodal Piecewise Linear II 0 −10 −20 0 2 4 Row 6 8 10 (b) Figure 4.2: Tumour-to-Skin Ratios of Models with Varying Tissue Distributions: The original and skin-subtracted ratios are shown in (a) and (b), respectively. The means are marked by various shapes while the errorbars represent the standard deviations. 107 Simple Glandular 0 −10 −10 −20 −20 −30 −30 −40 −40 −50 −50 −60 −60 −70 −80 −90 −70 Original Ideal Window Blind Window 0 2 −80 4 Row 6 8 Original Tumour−to−Skin Ratio (dB) Tumour−to−Signal Ratio (dB) 0 −90 10 (a) Piecewise Linear II −20 −30 −30 −40 −40 −50 −50 −60 −60 −70 −70 −80 −80 −90 −100 −110 −90 Original Ideal Window Blind Window 0 2 −100 4 Row 6 8 Original Tumour−to−Skin Ratio (dB) Tumour−to−Signal Ratio (dB) −20 −110 10 (b) Figure 4.3: Tumour-to-Signal Ratio of Models with Varying Tissue Distributions: Blind & ideal SDW selection is applied to Simple Glandular (a) & Piecewise Linear II (b). 108 Fig. 4.3 summarizes the tumour-to-signal ratios generated using the ideal and blind windows for the simple glandular and piecewise linear models, respectively. In both cases, the ratios generated using the ideal window are higher, which is anticipated as making the SDW larger calibrates the filter to remove a larger portion of the unwanted signal. The best average tumour-to-signal ratio for the simple glandular model is approximately -10 dB, while the highest mean ratio achieved by the piecewise linear model lies at about -30 dB. Next, the single signal analysis described in chapter 3 is conducted and images are formed with the known tumour signals and the resulting skin-subtracted data. Table 4.1 summarizes the number of preserved tumour responses, the SCR and the tumour localization for the skin-subtracted data generated using the blind SDW. Intuitively, as the glandular distribution of the breast models becomes more complicated, fewer tumour responses are successfully preserved. The SCR and tumour localization are compared to the images generated using the known tumour response. The SCR of the tumour estimate images are within 15% of the SCR values for the references images with the exception of the piecewise linear II case. This discrepancy is likely due to inaccuracies in the skin thickness estimate, which are subsequently used in the focusing algorithm. The tumour localization of both types of images is generally in agreement with the largest difference observed in the bimodal case with an error of 2.24 mm. It is also noted that more tumour responses are preserved using the blind SDW than the ideal window presented in Chapter 3. This observation is attributed to the fact that the ideal window is larger which increases the chances of including parts of the tumour response into the SDW. Consequently, the potential of distorting the 109 Table 4.1: Resulting Image Quality of the Same Breast Model with Varying Interior Properties: Each model has 245 signals and blind SDW selection is applied to the data. # Tumour SCR (dB) Localization (mm) Breast Model Responses Tumour Known Tumour Known Preserved Estimate Tumour Estimate Tumour Simple Breast 212 5.16 5.05 (60,110,84) (60,110,85) Simple Glandular 178 3.25 3.88 (60,110,85) (61,110,85) Bimodal 151 3.22 3.11 (59,110,84) (60,110,86) Piecewise Linear II 124 5.50 7.13 (61,110,82) (61,110,83) tumour estimate is increased. Finally, the filtered signals are focused and images are created. The tumour is not detected in the images generated for the models that contain glandular tissue, as the residual glandular reflections are to large. However, the primary skin reflection is successfully reduced in all four cases, indicating that the algorithm is robust to varying skin thickness and properties. 4.2 Application to Different Breast Shapes The Neighbourhood Modified RLS is tested on two breast models with different geometries to assess the robustness of the neighbourhood selection algorithm to variations in breast shape. The algorithm’s success is appraised based on its primary and secondary skin response reduction capabilities when the blind and ideal windows are applied. 110 (a) (b) Figure 4.4: Scan Pattern Applied to Simple Breast Models of Different Shapes: The models are derived off of MR images of the same patient; consequently, the model in (a) is called MRI3L, while the model in (b) is named MRI3R. 4.2.1 Model Two relatively simple breast models are examined for these tests. Both models are derived from MR images and are comprised of a 2 mm thick skin layer that encases a spherical tumour embedded in a homogenous fatty medium. Dispersion is not accounted for in both models, and the dielectric properties of the tissues remain consistent with previous trials. The first model is the simple breast model described in section 4.1.1 which is generated using an MR scan of the left breast of the patient; consequently, this is referred to as MRI3L for the remainder of this section. The second model is based on MR images of the right breast of the same patient and is thereby labeled MRI3R. There are two major differences between these models. First, the tumour diameter of MRI3L is 6 mm while the tumour in MRI3R spans 8 mm. Second, the overall shape of both models is substantially different, as illustrated by Fig. 4.4 and Fig. 4.5. Fig. 4.4 displays the skin models and scan patterns for both models. In attempt to 111 MRI3L MRI3R 50 20 20 40 35 60 30 25 80 20 100 15 120 x (mm) x (mm) 40 140 50 45 40 40 30 60 80 20 100 10 10 5 50 100 y (mm) (a) 150 200 120 50 100 y (mm) 150 200 (b) Figure 4.5: Relative Permittivity Map Simple Breast Models of Different Shapes: The z plane shown above intersects with the center of the tumour in both cases, and is located at the same elevation (z = 80 mm). MRI3L is shown in (a) while MRI3R is seen in (b). make an unbiased comparison of the geometry of both models, the MRI3R is rotated by 180 degrees to gain a similar view of both breasts and the limits of all three axes are the same. To maintain consistency between the trials, the data is acquired using the same antenna and excitation outlined in section 4.1.1. The scan pattern for both models is comprised of the same number of rows and antennas so that an objective analysis of the algorithm’s performance is conducted. These antennas are placed approximately 2 cm away from the breast model. Fig. 4.5 clearly demonstrates the difference in the shape between both models as the z-planes shown in both images are at z = 80 mm. 4.2.2 Results The performance measures outlined in section 3.2 are used to characterize the skin reduction capabilities of the algorithm. First, the single signal analysis of the results 112 generated using the ideal and blind window is conducted. Next, the focusing algorithm is applied to the skin-subtracted, tumour estimate, and known tumour signals and the resulting images are evaluated. The average and standard deviation of the tumour-to-skin ratios over every row, for both models, using the blind window are presented in Fig. 4.6. The overall trends of the original and processed tumour-to-skin ratios are generally the same within a respective model; however, a difference is observed when the trends of the two models are compared. This deviation is another demonstration of the shape variation between the MRI3R and MRI3L models. The performance of the Neighbourhood Modified RLS remains consistent, as the average improvements in p2p ratios between the skin and the tumour responses are 99 dB and 97 dB for the right and left breasts, respectively. 60 0 50 −10 40 −20 30 −30 20 −40 10 −50 0 −60 Original MRI3R Original MRI3L Blind Window MRI3R Blind Window MRI3L −10 −20 −30 0 1 2 3 4 Row Original Ratio (dB) Skin Subtracted Ratio (dB) Tumour−to−Skin Using Blind Window −70 −80 −90 5 6 7 8 9 Figure 4.6: Tumour-to-Skin Ratios of Models with Different Shapes 113 MRI3R 0 −10 −10 −20 −20 −30 −30 −40 −40 −50 −50 −60 −70 −80 −60 Original Blind Window Ideal Window 0 1 2 3 −70 4 Row 5 6 7 8 9 Original Tumour−to−Skin Ratio (dB) Tumour−to−Signal Ratio (dB) 0 −80 Figure 4.7: Tumour-to-Signal Ratio Using Blind & Ideal SDW on MRI3R Fig. 4.7 summarizes the tumour-to-signal ratios calculated for MRI3R model using the ideal and blind windows. As expected the ideal window has superior performance to the blind window in terms of reducing unwanted late-time responses. The maximum ratio achieved with the ideal window lies at about -10 dB while the rest of the ratios are generally about -15 dB. These values are comparable with the results seen for the simple glandular case presented in section 4.1.2; however, the tumour is successfully detected when images are formed using the skin-subtracted signals. Fig. 4.8 illustrates the tumour estimate and skin-subtracted images generated when the Neighbourhood Modified RLS is applied to MRI3R using the ideal window. 114 Tumour Estimate 20 0.8 40 0.6 60 0.4 80 0.2 100 0 50 100 y (mm) (a) 150 200 1 20 x (mm) x (mm) Skin Subtracted 1 0.8 40 0.6 60 0.4 80 0.2 100 0 50 100 y (mm) 150 200 (b) Figure 4.8: Images Generated From Model With Different Geometry: The site of the tumour is readily apparent amidst several sources of clutter in (b), while limited clutter sites are seen in (a). Table 4.2 summarizes the SCR and the tumour localization of the skin-subtracted, tumour estimate, and known tumour images generated using the ideal window. Reasonable agreement is observed between the resulting measures of the tumour estimate and reference images, as the localization error between the two is less than 1.5 mm and the SCR are within 15% of each other. The skin-subtracted images for both models have nearly identical attributes, as their localization error with reference to the known tumour image is the same and their SCR are almost equivalent. Table 4.2: Resulting Image Quality of Two Different Simple Breast Models. SCR (dB) Localization (mm) Breast Skin Tumour Known Skin Tumour Known Model Subtracted Estimate Tumour Subtracted Estimate Tumour Data Data MRI3L 1.80 4.80 5.05 (60,111,83) (60,110,84) (60,110,85) MRI3R 1.82 2.56 3.03 (50,113,81) (50,112,81) (50,111,82) The percentage of tumour responses preserved using the blind window for the 115 right and left models are 82.7 and 85.5, respectively. The ideal window preserves 70.5% and 73.2% of the tumour responses for the right and left breast models, respectively. These results are deemed acceptable for it is not realistic to measure a reliable tumour response at every antenna location. Once again the same trend is observed where the ideal window generates fewer accurate tumour estimates than the blind window. Although fewer tumour responses are preserved by the ideal window, the tumour is detected in images, because the late-time clutter is significantly reduced. The blind window does not suppress these unwanted responses as effectively; consequently, its signals require further processing in order to successfully focus an image. Overall, the robust nature of the Neighbourhood Modified RLS is demonstrated as it is able to suppress the skin reflections acquired from breast models with distinct geometry using both blind and ideal windows. The application of the ideal window also produces signals that are capable of generating an image of the tumour, when processed by the focusing algorithm. 4.3 Application to Different Antennas The algorithm’s performance is further explored by examining if it is effective when the data is acquired using a different antenna. Two breast models are investigated and the Neighbourhood Modified RLS skin response removal capabilities are assessed. As the radiating element has changed, the neighbourhood selection algorithm is updated to account for the beamwidth difference. 116 4.3.1 Model Both models are derived from MR images of the same patient and are comprised of skin, fat, and a 6 mm tumour. Once again dispersion is not included; however, the dielectric properties of the three tissues are slightly different then those previously employed. Table 4.3 summarizes the properties of the materials used in the models. Table 4.3: Dielectric Properties of Materials Used in Simulations Material Permittivity Conductivity Immersion 2.4 0.04 Medium Skin 36 4 Fat 8.18 0.4 Tumour 54.1 0.79 The scan pattern used for both trials is made up of 11 rows of 16 antennas as demonstrated in Fig. 4.9, and the antennas are placed approximately 2 cm away from the skin. The main difference between the two models is the location of the tumour. In one case the tumour is positioned towards the center of the model (MRI316R Central), while in second scenario the tumour is placed in close proximity to the skin (MRI316R Edge). Fig 4.10 illustrates the cross sections of the breast models at the tumour locations in the z-plane. A Balanced Antipodal Vivaldi Antenna (BAVA) is used to illuminate the breast for the data sets used in this section. This antenna is designed specifically for microwave breast cancer detection and is currently integrated into the TSAR system [51]. To instill consistency between the trials, the same differentiated Gaussian pulse is used to excite the antenna. The main differences between the BAVA and the (Wu-King) dipole antenna previously employed are highlighted in Table 4.4. 117 Figure 4.9: Scan Pattern Applied to Simple Breast Models Excited With BAVA Table 4.4: Differences Between the Wu-King Dipole and the BAVA Antenna Beamwidth (mm) Notes Wu King Dipole 44x35 - inefficient due to resistive loading - bandwidth undetermined (poor S11) BAVA 34x45 - implemented experimentally - bandwith of 2.4-18GHz MRI316R Central MRI316R Edge 55 55 50 50 50 50 45 45 40 35 30 150 25 40 100 x (mm) x (mm) 100 35 30 150 25 20 200 15 20 200 15 10 250 5 50 100 150 y (mm) (a) 200 250 10 250 5 50 100 150 y (mm) 200 250 (b) Figure 4.10: Relative Permittivity Map Simple Breast Models Excited With BAVA: The z planes shown above intersect with the center of the tumour in both cases. The elevation of the z plane in both cases are different, in (a) z = 93 mm while z = 93 mm in (b). 118 4.3.2 Results Using the single signal analysis and the image evaluation tools outlined in section 3.2, the performance of the Neighbourhood Modified RLS is appraised. First, the difference between the magnitudes of the tumour responses in both models is explored. Next, the resulting tumour-to-skin and tumour-to-signal ratios are examined. Last, the percentage of accurate tumour estimates is presented along with the SCR and tumour localizations of the generated images. Fig. 4.11 illustrates the largest and median tumour responses of both models in terms of their p2p values. The biggest tumour response is two orders of magnitude larger than median tumour response in the edge model, while the difference observed for the central case is only a factor of 5. The significant amplitude deviation observed for the edge model is directly related to propagation losses. Consequently, substantial variations are observed in the tumour responses of the edge case, because the distance between the antenna and tumour varies considerably. The consistent nature of the magnitude of the central model’s tumour responses is also reflected in the original tumour-to-skin ratios presented in Fig. 4.12, as the standard deviations of the mean values are significantly smaller than those observed in the edge case. The mean tumour-to-skin ratios after subtraction are also shown in Fig. 4.12 for both models. The ratios generated using the blind windows are all generally above 0 dB for the central case, indicating that the skin response is successfully reduced. The negative ratios recorded at rows 8, 9, and 11 are caused by the limited size of the neighbourhoods in that region, in conjunction with the small p2p values exhibited by many of the tumour responses. The reduction in tumour- 119 MRI316R Central Largest Tumour Response −6 5 0 −5 0 1 2 3 MRI316R Edge Largest Tumour Response −5 Voltage Voltage 5 x 10 4 x 10 0 −5 0 5 1 2 3 4 −9 0 1 2 Time (s) 3 Median Tumour Response −7 5 Voltage Voltage x 10 −1 0 x 10 Median Tumour Response −6 1 4 5 −9 x 10 (a) 5 −9 x 10 x 10 0 −5 0 1 2 Time (s) 3 4 5 −9 x 10 (b) Figure 4.11: Comparison of Tumour Responses Recorded Using BAVA: The amplitude difference between the largest and the median tumour response recorded from the MRI316R Central model (b) is 5 times, while an amplitude variation of 2 orders magnitude is observed with the the MRI316R Edge model (b). to-skin ratio observed in the edge case, results from the larger number of tumour responses with very small p2p values. The ideal window normally generates lower tumour-to-skin ratios, because the SDW is typically larger which means it contains more complex signals. As a result, the estimate of the filter weight coefficients is less accurate. Fig. 4.13 demonstrates the mean tumour-to-signal ratios computed for the edge and central cases using the blind and ideal windows. Once again the ideal window has better average ratios than the blind window and the effect of having smaller neighbourhoods in rows 8, 9, and 11 is seen. The blind and ideal ratios lie between -70 to -40 dB for the central case and -80 to -50 dB for the edge data. The difference between the two models is once again due to the number of tumour responses with very small p2p value in the edge dataset. 120 20 −40 10 −50 0 −60 −10 −70 −20 −80 −30 −90 Original Ideal Window Blind Window −40 −50 −60 0 2 4 −100 Skin Subtracted Ratio (dB) Original Ratio (dB) Tumour−to−Skin Ratio MRI316 Central Tumour −110 6 Row 8 10 −120 12 (a) Original Ratio (dB) −30 20 −40 0 −50 −20 −60 −40 −70 −60 −80 −80 −90 −100 −100 −120 −110 Original Ideal Window Blind Window −140 −160 0 2 4 −120 Skin Subtracted Ratio (dB) Tumour−to−Skin Ratio MRI316 Edge Tumour 40 −130 6 Row 8 10 12 (b) Figure 4.12: Resulting Tumour-to-Skin Ratios with BAVA: The same model is test with two different tumour placements. The tumour is positioned towards the center of the model in (a) & at the edge of the model in (b). 121 −40 −50 −50 −60 −60 −70 −70 −80 −80 −90 −90 −100 Original Ideal Window Blind Window −100 0 2 4 Original Tumour−to−Skin Ratio (dB) Tumour−to−Signal Ratio (dB) Tumour−to−Signal MRI316 Central Tumour −110 6 Row 8 10 12 (a) −40 −30 −50 −40 −60 −50 −70 −60 −80 −70 −90 −80 −100 −90 −110 −100 −120 −110 Original Ideal Window Blind Window −130 −140 −150 0 2 4 −120 −130 6 Row 8 10 Original Tumour−to−Skin Ratio (dB) Tumour−to−Signal Ratio (dB) Tumour−to−Signal MRI316 Edge Tumour −140 12 (b) Figure 4.13: Resulting Tumour-to-Signal Ratios with BAVA: Blind & ideal SDW selection is applied to MRI316 Central in (a) & to MRI316 Edge in (b) 122 Table 4.5 summarizes the results of applying the single signal and image evaluation metrics on skin-subtracted data using the blind window. The localization of Table 4.5: Resulting Image Quality of Breast Models Scanned Using the BAVA. % Tumour SCR (dB) Localization (mm) Breast Model Responses Tumour Known Tumour Known Dectected Estimate Tumour Estimate Tumour MRI316R Central 69 3.22 0.52 (133,134,92) (133,134,91) MRI316R Edge 48 2.27 2.07 (161,117,89) (163,116,96) the tumour estimate and the known tumour images are in agreement for the central case, while an error of 7.35 mm is observed in the edge scenario. The significant error found in the edge case is due to the limited number of tumours that are successfully detected. The decreased detection rate exhibited by this model is a result of two factors. First, as the tumour is in close proximity to the skin, the tumour response and the skin response overlap in the recorded signals acquired from antennas close to the tumour. This leads to the inclusion of part of the tumour response into the SDW, thereby causing distortion to the tumour estimates. Second, the antennas that are positioned on the opposing side of the breast relative to the tumour record such small tumour responses that the reflections almost resemble noise. Consequently, the algorithm is unable to recover these responses. It is also observed that the SCR of known tumour image for the central case is very low. This phenomenon is caused by the level of noise present in the known tumour responses relative to the p2p value of the tumour reflection. As the tumour is central in this model, the p2p values of its tumour responses are fairly small. Consequently, the focusing algorithm is unable to effectively resolve the tumour response from noise present in the signal. This problem does not arise with the edge model because the 123 p2p values of the tumour responses recorded from antennas in close proximity to the tumour are significantly larger than the noise. Ultimately, the primary skin response is effectively removed for both models indicating that the Neighbourhood Modified RLS has the potential to be applied to data collected from a variety of antennas. The tumour cannot be imaged when focusing the skin-subtracted data, as significant artifacts are seen in the signals recorded by the BAVA. However, it is proven that the tumour responses are preserved, as the tumour is successfully imaged using the tumour estimates. The signal artifacts in the BAVA signals also cause the extremely small tumour-to-signal values computed and are discussed further in section 4.4.2. 4.4 Discussion The overall performance of the algorithm is further assessed by comparing the results of the breast models previously investigated. The adaptive nature of the neighbourhood selection algorithm is examined first, followed by an evaluation of the signal artifacts present prior to skin subtraction. These artifacts are characterized and related back to the challenge of skin subtraction. The skin reduction effectiveness is then quantified and analyzed. 4.4.1 Neighbourhood Comparison The neighbourhoods created by processing the models described in this Chapter are evaluated by comparing the sizes of the resulting groupings. This comparison is conducted by creating a histogram of the generated neighbourhood sizes (quantified in 124 terms of the number of antennas) presented in Fig. 4.14. As the two models scanned by the BAVA antenna have identical neighbourhoods, only one of the cases is illustrated. This result is expected as the same scan pattern is used to illuminate a breast of the same shape with a homogeneous interior tissue distribution. Furthermore the grid scan applied to the models introduced in section 4.1.1 is not included in the analysis to ensure that an unbiased comparison is made between the other models. First, the neighbourhoods of models that have the same shape but different dielectric properties and skin thicknesses (MRI3L and MRI3L Piecewise Linear II) are compared. The resulting neighbourhood sizes deviate slightly due to the variations in the skin thickness and dielectric properties. Conversely, when comparing models of different shapes (MRI3R and MRI3L) a significant difference is seen in the neighbourhood distribution even though the same number of antennas is used in both scans. This indicates that the shape of the breast has a large bearing on the antenna groupings chosen by the neighbourhood selection algorithm. It is also interesting to note that the neighbourhoods selected for the models illuminated using the BAVA are significantly smaller than those generated for all the other cases. This phenomenon is caused by the narrower beamwidth exhibited by the BAVA in combination with the less dense scan pattern used to acquire the data. 125 Neighbourhood Sizes Generated Using MRI3R 70 60 50 40 30 20 Neighbourhood Sizes Generated Using MRI3L 80 Quantity of Receivers with Given Neighbourhood Size Quantity of Receivers with Given Neighbourhood Size 80 10 70 60 50 40 30 20 10 0 0 10 20 30 40 50 60 70 80 90 Size of Neighbourhood (# of Antennas) 0 0 100 10 20 30 (a) 70 60 50 40 30 20 10 0 50 60 70 80 90 100 (b) Neighbourhood Sizes Generated Using MRI316R 80 Quantity of Receivers with Given Neighbourhood Size Quantity of Receivers with Given Neighbourhood Size 80 40 Size of Neighbourhood (# of Antennas) Neighbourhood Sizes MRI3L Piecewise Linear II 70 60 50 40 30 20 10 5 10 15 20 25 30 35 Size of Neighbourhood (# of Antennas) (c) 40 0 0 10 20 30 40 50 60 70 80 90 Size of Neighbourhood (# of Antennas) 100 (d) Figure 4.14: Histograms of Resulting Neighbourhood Sizes: The models used for (b) & (d) are introduced in Section 4.1.1 (MRI3L &MRI3L Piecewise Linear II), while the models in (a) and & (c) are decribed in Section 4.2.1 (MRI3R) & 4.3.1 (MRI316R), respectively. 4.4.2 Clutter Quantification Previously, the primary skin reflection was defined to be the dominant skin response that is generally captured within the SDW, while the secondary reflection includes all other responses beyond the window. However, the reflections recorded by the BAVA antenna contain signal artifacts, referred to as clutter, that were not previously seen using the dipole antenna. Fig. 4.15 illustrates a typical signal recorded by the BAVA 126 and dipole antennas and their corresponding tumour only signals. Calibrated Data 0 −5000 −10000 0 0.5 1 1.5 2 Calibrated Data 0.1 Voltage Voltage 5000 2.5 3 0 −0.1 0 1 2 3 4 −9 Tumour Only 1 0 −20 0 0.5 1 1.5 Time (s) (a) 2 2.5 3 −9 x 10 x 10 Tumour Only −6 Voltage Voltage 20 5 −9 x 10 x 10 0 −1 0 1 2 Time (s) 3 4 5 −9 x 10 (b) Figure 4.15: Comparison of Signals Acquired Using Different Antennas: Typical calibrated data, and tumour only signals recorded by the dipole (a) antenna and the BAVA (b). The additional clutter seen in the BAVA signal is caused by internal reflections within the antenna, and the fact that the simulations are run for a longer period of time. Consequently, a measure is developed to quantify the unwanted reflections known as clutter. This metric is the ratio in dB between the p2p value of the calibrated data within a window to the p2p value of the known tumour response. The purpose of placing the window on the calibrated data is to make an unbiased comparison of the clutter in the datasets recorded using the BAVA and dipole antennas. As a result, the window applied to the reflections acquired with the dipole is simply the rest of the signal beyond the SDW, while the responses collected with the BAVA are segmented based on a large late-time response seen in all the signals, as shown in Fig. 4.16. The large late-time response is postulated to be multiple reflections between the internal structures of the BAVA. This phenomenon is not seen in the dipole signals because the Wu-King Dipole is a lossy antenna that less complicated 127 Clutter Window Selection for BAVA Signals 0.04 Calibrated Data Clutter Window 0.03 Voltage 0.02 0.01 0 −0.01 −0.02 −0.03 0 0.5 1 1.5 2 2.5 Time (s) 3 3.5 4 4.5 −9 x 10 Figure 4.16: Selection of Window for Clutter Ratio Measure internal structures. Fig. 4.17 illustrates the difference in the tumour-to-clutter ratios calculated for the simple breast models illuminated with the BAVA and dipole antennas. As the known tumour response is used in this metric, the standard deviation of the MRI316R edge ratios are significantly larger than in any other model because the tumour is in close proximity to the skin. Overall, the clutter ratios for the signals recorded using the BAVA antenna are significantly lower, indicating that the tumour response is more heavily masked. This result explains the poor tumour-to-signal values calculated for the MRI316R models in section 4.3. 128 Simple Breast MRI3 (Dipole) −10 MRI3R MRI3L Tum. to Clutter Ratios (dB) −15 −20 −25 −30 −35 −40 −45 −50 −55 0 1 2 3 4 Row 5 6 7 8 9 (a) Simple Breast MRI316R (BAVA) −50 Tum. to Clutter Ratios (dB) −60 −70 −80 −90 −100 −110 Central Tum Edge Tum 0 2 4 6 Row 8 10 12 (b) Figure 4.17: Clutter Ratios Using Different Antennas: The clutter ratio from skin– subtracted data recorded from the simple breast models excited using the dipole (a) antenna & the BAVA (b). The blind SDW is used in the skin subtraction process. 129 4.4.3 Skin Response Suppression To assess the reduction of the skin response a new metric called the skin subtraction ratio (SkinSubRatio) is employed. This measure calculates the p2p value of the segment of the signal within the SDW before and after skin subtraction and takes the ratio of the two values as shown by equation 4.1: SkinSubRatio = 20 ∗ log p2pSkSubSW D . p2pCDataSDW (4.1) As the primary goal of the Neighbourhood Modified RLS is to reduce the skin reflection, focus is placed on the results generated using the blind window. The ideal window is not examined because it was developed for tumour detection purposes and is typically marginally less successful at reducing the primary skin response. Fig. 4.18 presents the mean and the standard deviation of the skin subtraction ratio over every row, for the simple breast datasets acquired using the BAVA and dipole antennas. The dominant skin response is reduced by about 96 dB on average in all four cases. The algorithm’s consistent nature is also seen with the models that include glandular tissue as their average skin subtraction ratio is 94 dB. Virtually no variations are seen in the average skin subtraction ratio values calculated for datasets associated with the BAVA while marginally larger deviations are seen in the other examined cases. This result is expected as the shape of the breast remains the same for the BAVA trials, while the geometry of the other models varies. Furthermore, the extremely small differences between the MRI316R models occur on the rows that are closest to the tumour, rows 7 to 11. The consistent reduction of the primary skin response in all the tested cases demonstrates the efficient and robust nature of the Neighbourhood Modified RLS. 130 Simple Breast MRI3 (Dipole) 106 MRI3R MRI3L Skin Subtraction Ratio (dB) 104 102 100 98 96 94 92 90 0 1 2 3 4 Row 5 6 7 8 9 (a) Simple Breast MRI316R (BAVA) 115 Skin Subtraction Ratio (dB) 110 105 100 95 90 Central Tum Edge Tum 85 80 0 2 4 6 Row 8 10 12 (b) Figure 4.18: Skin Subtraction Ratios Using Different Antennas: The skin subtraction ratio metric is applied to the simple breast datasets recorded with the dipole (a) antenna & the BAVA (b). The blind SDW is used in the skin subtraction process. 131 4.5 Summary The performance of the Neighbourhood Modified RLS technique is evaluated by testing its ability to reduce the skin response in breast models with varying glandular distributions, different shapes, and on data acquired using a different antenna. The adaptive nature of the neighbourhood selection technique is presented, as significant differences are seen when the breast geometry varies or when the antenna spacing and beamwidth change. The resulting groupings exhibit only a marginal deviation when the glandular tissue distributions of the breast models differ. The advantages and disadvantages of applying the blind and ideal windows during the filter weight estimation process are also examined. Applying the ideal window to datasets that are acquired by the Wu King Dipole antenna and do not include glandular tissue leads to successful tumour detection; however, the same results are not achieved with the BAVA, as too many artifacts are present in the recorded signal. The blind window, on the other hand, attains the primary goal of reducing the primary skin response and subduing the secondary skin response on a more consistent basis, as an improvement of 94 to 96 dB is observed from the variety of tested breast models. 132 Chapter 5 Conclusions TSAR is a radar based microwave imaging approach that has been proposed as a complementary breast cancer detection technique. This modality relies on the dielectric property contrast between tissues to detect malignant lesions. Data is collected by illuminating the breast with an ultra-wideband pulse and recording the ensuing reflections at several locations. The recorded signals are then processed to isolate the reflections from the tumour. One of the major challenges faced by this technique is removing the overwhelming skin reflections that mask the tumour responses. Several skin subtraction algorithms have been successfully tested; however, a method applicable to signals acquired from a realistic 3-D breast model that is compatible with the TSAR data collection system has not been reported. This thesis focuses on the development of a skin subtraction algorithm for monostatic data collected from realistic, 3-D MR-based, numerical breast models. The new skin reduction technique is a fusion of a neighbourhood selection algorithm with an updated version of the Modified RLS called the Neighbourhood Modified RLS. The neighbourhood identification technique selects an appropriate group of antennas to be used to create the input vector for the Modified RLS. The selection is conducted based on the beamwidth of the antenna and the proximity of the sensors, in conjunction with the normalized cross correlation between the recorded signals. The updates to the Modified RLS include the integration of an adaptive SDW selection approach and a rank determination technique. The previously proposed feedback 133 filter is also removed from the Modified RLS. Two SDW selection techniques are explored in depth: the blind and the ideal methods. The ideal technique requires prior knowledge of the tumour location, while the blind window does not require any prior information. The tumour responses are more consistently preserved and the dominant skin response is more effectively reduced using the blind window. However, when the neighbourhood is appropriately selected, the ideal window has the positive side effect of producing signals that allow the tumour to be imaged without any further signal processing. The performance of the Neighbourhood Modified RLS is tested on a variety of realistic, 3-D MR-based breast models. The first set of models exhibit varying glandular distributions, while the second set incorporates different breast shapes. Both of these trials scan the breast models using a (Wu-King) dipole antenna. The data for the third set was acquired using the BAVA on two models with the same shape and properties, but with different tumour locations. Based on these trials, the neighbourhood selection technique is found to be robust for three reasons. First, the antenna groupings do not change significantly when the glandular distribution is changed. Second, the primary skin response is practically eliminated in all tested cases. Third, the tumour is successfully imaged in a simple breast model with a different shape using the ideal window. The adaptive nature of the neighbourhood identification approach is also demonstrated, as the size of the neighbourhoods change significantly when the breast geometry is varied or when the antenna scan pattern and beamwidth changes. Overall, the Neighbourhood Modified RLS proved to be an effective algorithm that significantly reduces the primary skin response by approximately 95dB on a consistent basis. 134 Previously, skin-subtraction algorithms have been reported, that reduce the skin response by up to 100dB; however, these techniques were not tested on realistic 3-D models [49]. To the writer’s knowledge, the Neighbourhood Modified RLS is the first algorithm to have successfully skin-subtracted signals acquired monostatically from realistic, 3-D MR-based, numerical breast models. The main contributions of this thesis include: 1. providing an understanding of the data recorded from simulation of 3-D realistic numerical breast models, 2. developing and testing a similarity criteria for data that is helpful in several signal processing challenges, is applicable to signals recorded from monostatic and multistatic systems, and is also adaptive to the sensor, and 3. enhancing a pre-existing skin subtraction algorithm so that it is almost entirely adaptive to the data, and is effective on signals recorded from 3-D MR-based numerical breast models. Several aspects of the algorithm can be explored further. First, the Neighbourhood Modified RLS should be tested on data acquired in experimental settings, as realistically shaped 3-D breast models are now available and pre-clinical trials are becoming more prevalent[45]. Next, the neighbourhood identification technique could potentially be used to enhance image generation as information about the antenna beam overlap may help decipher the origin of responses. Last, the neighbourhood selection algorithm may also be applied to multistatic systems, as the relationship between the signals recorded at two different sites may be reconciled using the antenna beamwidth [47]. Bibliography [1] “Canadian cancer statistics 2007,” Canadian Cancer Society/National Cancer Institute of Canada, Tech. Rep., 2007. [2] A. Smith, P. Hall, and D. Marcello, “Emerging technologies in breast cancer detection.” Radiology Management, vol. 17, no. 6, pp. 16–24, Jul. 2004. [3] K. Foster, “Thermographic detection of breast cancer,” IEEE Engineering in Medicine and Biology Magazine, vol. 17, no. 6, pp. 10–14, 1998. [4] S. Nass, I. Henderson, and J. Lashof, Mammography and Beyond: Developing Techniques for the Early Detection of Breast Cancer. Institute of Medicine, National Academy Press, Washington, DC, USA, 2001. [5] A. Campbell, “Dielectric properties of female human breast tissue measured in vitro at 3.2 GHz,” Physics in Medicine and Biology, vol. 37, no. 1, pp. 193–210, 1992. [6] S. Chaudhary, R. Mishra, A. Swarup, and J. Thomas, “Dielectric properties of normal & malignant human breast tissues at radiowave & microwave frequencies.” Indian Journal of Biochemistry & Biophysics, vol. 21, no. 1, pp. 76–79, 1984. [7] J. Choi, J. Cho, Y. Lee, J. Yim, B. Kang, K. Keun Oh, W. Hee Jung, H. Jung Kim, C. Cheon, H. Lee, et al., “Microwave detection of metastasized breast cancer cells in the lymph node; potential application for sentinel lym- 135 136 phadenectomy,” Breast Cancer Research and Treatment, vol. 86, no. 2, pp. 107– 115, 2004. [8] W. Joines, Y. Zhang, C. Li, and R. Jirtle, “The measured electrical properties of normal and malignant human tissues from 50 to 900 MHz,” Medical Physics, vol. 21, pp. 547–550, 1994. [9] M. Lazebnik, L. McCartney, D. Popovic, C. Watkins, M. Lindstrom, J. Harter, S. Sewall, A. Magliocco, J. Booske, M. Okoniewski, et al., “A large-scale study of the ultrawideband microwave dielectric properties of normal breast tissue obtained from reduction surgeries,” Physics in Medicine and Biology, vol. 52, no. 10, pp. 2637–2656, 2007. [10] M. Lazebnik, D. Popovic, L. McCartney, C. Watkins, M. Lindstrom, J. Harter, S. Sewall, T. Ogilvie, A. Magliocco, T. Breslin, et al., “A large-scale study of the ultrawideband microwave dielectric properties of normal, benign and malignant breast tissues obtained from cancer surgeries,” Physics in Medicine and Biology, vol. 52, no. 20, pp. 6093–6116, 2007. [11] L. Wang, X. Zhao, H. Sun, and G. Ku, “Microwave-induced acoustic imaging of biological tissues,” Review of Scientific Instruments, vol. 70, pp. 3744–3748, 1999. [12] E. Fear, S. Hagness, P. Meaney, M. Okoniewski, and M. Stuchly, “Enhancing breast tumor detection with near-field imaging,” IEEE Microwave Magazine, vol. 3, no. 1, pp. 48–56, 2002. 137 [13] R. Tipa and O. Baltag, “Microwave thermography for cancer detection,” Romanian Journal of Physics, vol. 51, no. 3/4, pp. 371–377, 2006. [14] P. Kosmas and C. Rappaport, “FDTD-based time reversal for microwave breast cancer detection-localization in three dimensions,” IEEE Transactions on Microwave Theory and Techniques, vol. 54, no. 4 Part 2, pp. 1921–1927, 2006. [15] T. Williams, E. Fear, and D. Westwick, “Tissue sensing adaptive radar for breast cancer detection-investigations of an improved skin-sensing method,” IEEE Transactions on Microwave Theory and Techniques, vol. 54, no. 4 Part 1, pp. 1308–1314, 2006. [16] E. Bond, X. Li, S. Hagness, and B. Van Veen, “Microwave imaging via spacetime beamforming for early detection of breast cancer,” IEEE Transactions Antennas and Propagation, vol. 51, no. 8, pp. 1690–1705, Aug. 2003. [17] J. Sill and E. Fear, “Tissue sensing adaptive radar for breast cancer detectionexperimental investigation of simple tumor models,” IEEE Transactions on Microwave Theory and Techniques, vol. 53, no. 11, pp. 3312–3319, 2005. [18] D. Kurrant, E. Fear, and D. Westwick, “Tumor Response Estimation in Radarbased Microwave Breast Cancer Detection,” IEEE Transactions on Biomedical Engineering, vol. 55, no. 12, pp. 2801–2811, 2008. [19] P. Kosmas and C. Rappaport, “A matched-filter FDTD-based time reversal approach for microwave breast cancer detection,” IEEE Transactions on Antennas and Propagation, vol. 54, no. 4, pp. 1257–1264, 2006. 138 [20] ——, “Time reversal with the FDTD method for microwave breast cancer detection,” IEEE Transactions on Microwave Theory and Techniques, vol. 53, no. 7, pp. 2317–2323, 2005. [21] R. Nilavalan, A. Gbedemah, I. Craddock, X. Li, and S. Hagness, “Numerical investigation of breast tumour detection using multi-static radar,” Electronics Letters, vol. 39, no. 25, pp. 1787–1789, 2003. [22] Y. Xie, B. Guo, L. Xu, J. Li, and P. Stoica, “Multi-static adaptive microwave imaging for early breast cancer detection,” IEEE Transactions on Biomedical Engineering, vol. 53, pp. 1647–1657, 2006. [23] S. Hagness, A. Taflove, and J. Bridges, “Two-dimensional FDTD analysis of a pulsed microwave confocal system for breast cancer detection: Fixed-focus and antenna-array sensors,” IEEE Transactions on Biomedical Engineering, vol. 45, no. 12, pp. 1470–1479, 1998. [24] P. Kosmas, C. Rappaport, and E. Bishop, “Modeling with the FDTD method for microwave breast cancer detection,” IEEE Transactions on Microwave Theory and Techniques, vol. 52, no. 8 Part 2, pp. 1890–1897, 2004. [25] I. Rekanos and A. Raisanen, “Microwave imaging in the time domain of buried multiple scatterers by using an FDTD-based optimization technique,” IEEE Transactions on Magnetics, vol. 39, no. 3 Part 1, pp. 1381–1384, 2003. [26] X. Chen, D. Liang, and K. Huang, “Microwave imaging 3-D buried objects using parallel genetic algorithm combined with FDTD technique,” Journal of Electromagnetic Waves and Applications, vol. 20, no. 13, pp. 1761–1774, 2006. 139 [27] X. Li, E. Bond, B. Van Veen, and S. Hagness, “An overview of ultra-wideband microwave imaging via space-time beamforming for early-stage breast-cancer detection,” IEEE Antennas and Propagation Magazine, vol. 47, no. 1, pp. 19– 34, 2005. [28] X. Li and S. Hagness, “A Confocal Microwave Imaging Algorithm for Breast Cancer Detection,” IEEE Microwave and Wireless Component Letters, vol. 11, no. 3, pp. 130–132, 2001. [29] M. Converse, E. Bond, B. Veen, and C. Hagness, “A computational study of ultra-wideband versus narrowband microwave hyperthermia for breast cancer treatment,” IEEE Transactions on Microwave Theory and Techniques, vol. 54, no. 5, pp. 2169–2180, 2006. [30] E. Zastrow, S. Davis, M. Lazebnik, F. Kelcz, B. Van Veen, and S. Hagness, “Development of Anatomically Realistic Numerical Breast Phantoms with Accurate Dielectric Properties for Modeling Microwave Interactions with the Human Breast,” IEEE Transactions on Biomedical Engineering, vol. 55, no. 12, pp. 2792–2800, 2008. [31] T. C. Williams, “Radar based breast cancer detection: Skin estimation,” Ph.D. dissertation, University of Calgary, Dept. Elect & Comp. Eng., Calgary, AB, 2008. [32] J. Sill, T. Williams, E. Fear, R. Frayne, and M. Okoniewski, “Realistic Breast Models for Second Generation Tissue Sensing Adaptive Radar System,” in EuCAP 2007-The Second European Conference on Antennas and Propagation, 140 2007, 2007, pp. 1–4. [33] “Anatomical terminology.” [Online]. Available: http://training.seer.cancer. gov/anatomy/body/terminology.html [34] D. Kurrant and E. Fear, “An Improved Technique to Predict the Time-of-Arrival of a Tumor Response in Radar-Based Breast Imaging.” IEEE Transactions on Biomedical Engineering, vol. 56, no. 4, pp. 1200–1208, 2009. [35] J. M. Sill, personal communication, 2009. [36] S. Davis, E. Bond, X. Li, S. Hagness, and B. Veen, “Microwave Imaging via Space-Time Beamforming for Early Detection of Breast Cancer: Beamformer Design in the Frequency Domain,” Journal of Electromagnetic Waves and Applications, vol. 17, no. 2, pp. 357–381, 2003. [37] S. Davis, H. Tandradinata, S. Hagness, and B. Van Veen, “Ultrawideband microwave breast cancer detection: a detection-theoretic approach using the generalized likelihood ratio test,” IEEE Transactions on Biomedical Engineering, vol. 52, no. 7, pp. 1237–1250, 2005. [38] D. Winters, E. Bond, B. Van Veen, and S. Hagness, “Estimation of the frequency-dependent average dielectric properties of breast tissue using a timedomain inverse scattering technique,” IEEE Transactions on Antennas and Propagation, vol. 54, no. 11 Part 2, pp. 3517–3528, 2006. [39] E. Fear, X. Li, S. Hagness, and M. Stuchly, “Confocal microwave imaging for breast cancer detection: Localization of tumors in three dimensions,” IEEE 141 Transactions on Biomedical Engineering, vol. 49, no. 8, pp. 812–822, 2002. [40] E. Fear and J. Sill, “Preliminary investigations of tissue sensing adaptive radar for breast tumor detection,” in Proceedings of the 25th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 2003., vol. 4, 2003, pp. 3787–3790. [41] T. Williams, J. Sill, and E. Fear, “Breast surface estimation for radarbased breast imaging systems,” IEEE Transactions on Biomedical Engineering, vol. 55, no. 6, pp. 1678–1686, 2008. [42] I. Craddock, A. Preece, J. Leendertz, M. Klemm, R. Nilavalan, and R. Benjamin, “Development of a hemi-spherical wideband antenna array for breast cancer imaging,” in EuCAP 2006-First European Conference on Antennas and Propagation, 2006., 2006, pp. 1–5. [43] M. Klemm, I. Craddock, J. Leendertz, A. Preece, and R. Benjamin, “Breast cancer detection using symmetrical antenna array,” in EuCAP 2007-The Second European Conference onAntennas and Propagation, 2007., 2007, pp. 1–5. [44] I. Craddock, M. Klemm, J. Leendertz, A. Preece, and R. Benjamin, “Development and application of a UWB radar system for breast imaging,” in LAPC 2008-Loughborough Antennas and Propagation Conference, 2008., 2008, pp. 24– 27. [45] M. Klemm, I. Craddock, J. Leendertz, A. Preece, and R. Benjamin, “RadarBased Breast Cancer Detection Using a Hemispherical Antenna Array— 142 Experimental Results,” IEEE Transactions on Antennas and Propagation, vol. 57, no. 6, pp. 1692–1704, 2009. [46] R. Nilavalan, J. Leendertz, I. Craddock, A. Preece, and R. Benjamin, “Numerical analysis of microwave detection of breast tumours using synthetic focussing techniques,” in IEEE Antennas and Propagation Society International Symposium, 2004, vol. 3, 2004, pp. 2440–2443. [47] M. O’Halloran, E. Jones, and M. Glavin, “Quasi-Multistatic MIST Beamforming for the Early Detection of Breast Cancer.” To appear in IEEE Transactions on Biomedical Engineering, 2009. [48] S. Haykin, Adaptive filter theory. Prentice-Hall, Inc. Upper Saddle River, NJ, USA, 1996. [49] J. M. Sill, “Second generation experimental system for tissue sensing adaptive radar,” Master’s thesis, University of Calgary, Dept. Elect & Comp. Eng., Calgary, AB, 2004. [50] G. Golub, M. Heath, and G. Whaba, “Generalized Cross Validation Technique as a Method for Choosing a Good Ridge Parameter,” Technometrics, vol. 21, no. 2, pp. 215–223, 1979. [51] J. Bourqui, “Balance antipodal vivaldi antenna and dielectric director for breast cancer detection,” Master’s thesis, University of Calgary, Dept. Elect & Comp. Eng., Calgary, AB, 2008. 143 [52] E. Zastrow, S. Davis, and S. Hagness, “Safety assessment of breast cancer detection via ultrawideband microwave radar operating in pulsed-radiation mode,” Microwave and Optical Technology Letters, vol. 49, no. 1, pp. 221–225. 144 Appendix A MIST:Mapping MR Intensities to Dielectric Properties Three basic methods have been proposed to map MR pixel intensities of the breast interior to dielectric properties: uniform mapping, piecewise linear mapping, and bimodal mapping [29]. All of these techniques use a first order Debye dispersion model to simulate the frequency dependent behaviour of breast tissue; however, the correlation of pixel intensities to an appropriate set of Debye parameters differs [28, 16, 27, 29, 52, 30]. Before any of these techniques are applied, a pixel intensity histogram of the MR template is created, as shown in Fig. A.1(b). (a) (b) Figure A.1: MR Image Used to Create a Numerical Breast Model: The MR template is acquired in the coronal plane as shown in (a) its corresponding pixel intensity histogram is illustrated in (b) [29]. Uniform mapping was developed first and does not effectively account for the large 145 contrast in dielectric properties between fatty and fibroglandular tissue found by the large scale breast tissue study [28, 29]. This method involves applying a smooth linear map between the MR pixel intensity range to a range of Debye parameters [29]. The spread of the Debye parameters span a certain percentage above and below the published median values for the tissue of interest. In [28, 16] this percentage was 10% about a relative permittivity of 9.8 and a conductivity of 0.4 S/m, while in [29] a threshold of 50% about a relative permittivity of 15.66 and a conductivity of 1.03 S/m was employed. Fig. A.2(a) demonstrates the resulting histogram of the dielectric constant at 6 GHz, when uniform mapping is applied to the MR image in Fig. A.1(a). Piecewise linear mapping accounts for the disparity between the dielectric properties of fatty, and fibroglandular tissues by assigning a distinct spread of Debye parameters for each tissue type [29]. This process starts by visually segmenting the histogram seen in Fig. A.1(b) into three regions: fatty, fibroglandular, and transitional [29]. The median intensity in the fatty region is mapped to the median Debye parameters documented for fat. The maximum and minimum pixel intensities in that region are then mapped to Debye parameters that are 10% about the median values [29]. The same procedure is applied to the pixels associated with fibroglandular tissue, while using the appropriate Debye values. The transition region is mapped to parameters that lie within the range of maximum fatty tissue and minimum fibroglandular [29]. Fig. A.2(b) demonstrates the resulting histogram of the dielectric constant at 6 GHz, when piecewise linear mapping is applied to the data in Fig. A.1(a). Similar to the piecewise linear mapping, the bimodal method assumes that two 146 distinct tissue types are present; however during the segmentation of the MR histogram, a single threshold is used thereby eliminating the transitional region [29]. Once the two regions have been identified, the same method as described for the piecewise linear technique is used to map the pixel intensities to frequency dependent dielectric values [29]. The median Debye parameters and ranges used for fatty and fibroglandular tissues in the piecewise linear technique are also applied here. This approach ultimately leads to a sharper dielectric contrast between fatty and fibroglandular tissue [29]. Fig. A.2(c) demonstrates the resulting histogram of the dielectric constant at 6 GHz, when bimodal mapping is applied to the data from Fig. A.1(a). Most recently in [30], the filtering, segmenting, and dielectric property mapping methods employed to create realistic three dimensional breast models have been updated. In these MR models, the skin is artificially introduced as a 1.5 mm thick layer that overlays the breast and a 0.5 cm chest muscle layer that is covered by 1.5 cm thick subcutaneous fat layer at the top of the model [30]. Fig. A is an example of a preprocessed MR template. The same previously described piecewise linear mapping scheme is used to create the numerical models; however, the intensity histogram is no longer segmented by inspection and the dielectric property values are assigned based on the three tissue groups outlined [9]. By fitting the intensity histogram to a Gaussian mixture model, 7 parameters are identified and used to define distinct regions of the histogram. Each region is associated with a different tissue type [30]. The 7 parameters are characterized using the minimum(mg ) and maximum (Mf ) intensity values, as well as the means(µf , µg ) and standard deviations (mσg , Mg , mf , mσf ) of the two Gaussians 147 Intensity Range mg − mσg mσg − µg µ g − Mg Mg − m f m f − µf µf − mσf mσf − Mf Intensity Group Ig1 Ig2 Ig3 Ig4 Ig5 Ig6 Ig7 Property Group Pg1 Pg2 Pg3 Pg4 Pg5 Pg6 Pg7 Table A.1: Corresponding Pixel Intensity and Dielectric Property Groups found in the mixture model [30]. Fig. A.4(a) shows an example of an intensity histogram, along with the resulting Gaussian mixture model while Fig. A.4(b) presents the Gaussians used to create it. The upper and lower limits of the dielectric properties of a given tissue group are defined in Fig. A.5(a). Once the property range is set, the pixel intensities in the corresponding regions are linearly mapped to dielectric values [30]. Table A.1 outlines the pixel intensity ranges as well the corresponding property maps while Fig. A.5(b) visually demonstrates the linear mapping at 6 GHz. Once all the preprocessed slices of the MR image are mapped, the planes are stacked on top of each other to create a three dimensional breast model [30]. 148 (a) (b) (c) Figure A.2: Dielectric Property Histograms Using Different Mapping Approaches: The results of applying the uniform method is displayed in (a), while the dielectric property histograms for the piecewise linear and bimodal techniques are demonstrated in (b) and (c), respectively [29]. 149 Figure A.3: Preprocessed MR Template Used for Dielectric Mapping [30] (a) (b) Figure A.4: Gaussian Mixture Modeling of MR Intensity Histogram: The original histogram and its Gaussian mixture estimate are shown in (a), while (b) illustrates the Gaussians used to create the estimate [30]. 150 (a) (b) Figure A.5: Linear Mapping of Pixel Intensities to Dielectric Properties: The linear map seen in (b) is used at 6 GHz to establish which dielectric property group demonstrated in (a) is associated with a given pixel intensity [30]. 151 Appendix B Derivation of the Kalman Gain Vector For the RLS Algorithm The RLS algorithm models the statistical measures encountered when trying to minimize the MSE with deterministic time-averaged values. Consequently, the MSE is replaced with a weighted sum of squared errors, while the autocorrelation and cross correlation are replaced with deterministic time averaged values [48]. Eq.(B.1) shows the time averaged equivalent of the normal equation, F[n]~ w[n] = ~z[n] → ~ [n] = F−1 [n]~z[n] w (B.1) where F[n] and ~z[n] are the deterministic time average values of the autocorrelation and cross correlation, respectively. These time-based measures can be calculated using the following relationship, F[n] = λF[n − 1] + ~u[n]~uT [n] → ~ = λ~z[n − 1] + d[n]~uT [n] z[n] (B.2) where F[n-1] and ~z[n−1] are the deterministic time averaged autocorrelation matrix and cross correlation vector of the last iteration and λ is the forgetting factor. The forgetting factor is the weight that is applied to the previously mentioned sum of squared errors and is a value between and including 0 and 1 [48]. The filter weights are established using Eq.(B.1); however, this requires the inversion of a matrix which is computationally intensive. As a result, matrix inversion 152 lemma is applied to Eq.(B.2) to create the following analytical expression for the inverse of F[n] [48]. λ−2 P[n − 1]~u[n]~uT [n]P[n − 1] 1 + λ−1 ~uT [n]P[n − 1]~u[n] = λ−1 P[n − 1] − λ−1~k[n]~uT [n]P[n − 1] F−1 [n] = P[n] = λ−1 P[n − 1] − (B.3) ~k[n] is known as the Kalman gain vector and can be calculated as follows: ~k[n] = P[n−1]~ u[n] . λ+~ uT [n]P[n−1]~ u[n] The gain vector is used to update the filter tap weights in order to avoid calculating the time averaged cross-correlation vector ~z[n] at every time step. Using a combination of the Kalman gain vector and Eq.(B.1-B.3), the following update equation for the weight vector is computed ~ [n] = w ~ [n − 1] + ~k[n]α[n] w (B.4) ~ [n − 1] is the weight vector of the past iteration while α[n] is the a priori where w error before the current weight [48]. α[n] is characterized by Eq.(B.5), α[n] = x[n] − ~uT [n]~ w[n − 1] where x[n] is the value of the target signal at the current time step. (B.5)

1/--страниц