DIGITIZATION AND ANALYSIS OF MAMMO GRAPHIC IMAGES FOR EARLY DETECTION OF BREAST CANCER By Farzin Aghdasi B. Sc. (Eng. Hon.) Imperial College of Science and Technology, University of London, London, U.K., 1977 M. B. A. University of Portland, Portland, Oregon, U.S.A., 1978 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES ELECTRICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 1994 © Farzin Aghdasi, 1994 ___________________________ In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of IL E The University of British Columbia Vancouver, Canada Date DE-6 (2188) -tC# iX I91L- Abstract X-ray mammography is the proven method for early detection of breast cancer. Digital pro cessing and analysis of mammographic images can potentially assist in improved performance of radiologists in earlier detection and recognition of abnormalities. In this work a novel image acquisition system based on an area scanning CCD array has been developed for the digitization of mammograms at high spatial and photometric resolutions. The system characteristic parameters were measured. The quality of the resulting images in terms of sharpness and noise content is comparable with that obtained by the more expensive and slower drum laser-scanning microdensitometer. The clinical application of soft-copy display of digitized images are evaluated. To further improve the quality of the images, restoration algorithms were applied to restore the images from the degrading effects of the system’s blur and noise. Performance of three filtering techniques was compared. A new method for the reduction of boundary truncation artifacts in image restoration was suggested and studied. The process of radiographic image formation was modeled and two locally adaptive smooth ing filters were employed to counter signal-dependent radiographic noise before application of restoration filters. The results of the restored images show a marked improvement in detectabil ity of smallest particles of microcalcifications when judged by a human observer. Image segmentation routines were developed to separate microcalcifications from the back ground parenchymal pattern. Performances of two algorithmic approaches to segmentation and two artificial neural networks were compared. Over 100 numerical features were automatically extracted from the clusters of microcalcifications. These features were evaluated for their ability to separate the benign and malignant formations. Using a database of 68 digitized mammo grams a discriminant function was calculated. The sensitivity and specificity of this approach in recognition of malignant microcalcification clusters is shown to be comparable to that of trained radiologists. 11 Table of Contents Table of Contents iii List of Tables viii List of Figures x Acknowledgment 1 2 3 xiii Introduction 1 1.1 Motivation 1 1.2 Human Female Breast 2 1.3 Breast Cancer 3 1.4 Mammographic Presentation 3 1.5 Objectives 1.6 The Structure of This Thesis 5 1.7 Claims to originality 7 State of the Art in Digital Mammogram Processing and Analysis 8 2.1 Image enhancement 8 2.2 Risk assessment 9 2.3 Automated detection of abnormalities 9 2.4 Differential Diagnosis 10 Radiographic Film Digitization 3.1 11 Image Acquisition Hardware 3.1.1 12 Microdensitometers 12 111 3.1.2 Video Cameras 3.1.3 CCD Sensors 13 3.2 Performance Characteristics 19 3.3 Noise Reduction 23 3.4 Comparison of the Three Film Digitization Systems 25 3.5 Image Acquisition at Multiple Scales 29 3.6 The Software Environment 30 3.6.1 MAMPRO 31 Mammogram Data Base 32 3.7 4 5 13 Clinical Evaluation of the Film Digitization System 36 4.1 Working Hypothesis 37 4.2 Materials and Methods 38 4.3 Results 39 4.4 Discussion 42 Image Restoration 46 5.1 Problem Description 46 5.2 Image Processing 47 5.3 Image Restoration 48 5.4 The System Transfer Function 50 5.5 The Wiener Filter 50 5.5.1 6 . Iterative Restoration 52 5.6 The Constrained Least Squares Filter 54 5.7 Results 55 Restoration in the Presence of Signal-Dependent Noise 62 6.1 Image Formation System 63 6.1.1 63 The X-ray Screen-Film Imaging System iv 6.1.2 7 The CCD Camera 67 6.2 Image Observation Model 69 6.3 Image Restoration 77 Local adaptive Wiener smoothing filter 78 6.3.2 Bayesian smoothing filter 79 6.3.3 The Deconvolution Filter 80 6.4 Results 80 6.5 Summary 83 Reduction of Boundary Artifacts in Image Restoration 85 7.1 Image Restoration Artifacts 86 7.2 Problem description 88 7.3 7.4 8 6.3.1 7.2.1 Analysis of errors in using the conventional linear de-convolution 7.2.2 Restoration by circular de-convolution . 92 . 97 Restoration by image extension and periodic convolution • 100 7.3.1 • 104 Computational Load Experiments • 107 Segmentation 118 8.1 Image Analysis • 118 8.2 Detection and Segmentation of Microcalcification Clusters • 119 8.3 Local Histogram Thresholding Algorithm 8.3.1 . . . 119 . Object Labeling • 124 8.4 Edge Detection and Region Growing Algorithm • 126 8.5 Neural Networks Techniques • 127 8.6 8.5.1 The Modified Hopfield Network 8.5.2 The Feed Forward Neural Network Results and Discussions • . 129 • 132 • V • • . 133 8.7 9 Summary 135 Classification 137 9.1 Data Base 138 9.2 Features 139 9.3 Analysis • . 9.4 Results • . 9.5 Summary • . 10 Conclusions, and recommendations for future work 143 146 147 10.1 Objectives • . 10.2 Summary of the work • . 10.3 Suggestions for future work 143 147 147 150 . Bibliography 152 A Mammography 165 A.1 Mammographic Image Formation System 165 A.2 Sources of Image Blur In Mammography 167 A.2.1 Motion Blur 167 A.2.2 Geometric Distortion 167 A.2.3 The Radiographic Receptor . 168 A.3 The System Characteristic Parameters 171 . A.3.1 The Modulation Transfer Function 171 A.3.2 Contrast 171 A.3.3 Noise 176 A.3.4 The X-ray Scatter 176 A.3.5 Digitization 178 A.4 Mammographic Image Interpretation . vi 181 A.5 Radiographically Occult Cancers . 182 B Questionnaire For Clinical Evaluation of Images 188 C Features 190 C.1 Features from Individual Microcalcifications 190 C.2 Features from Microcalcification Clusters 190 VI’ List of Tables 3.1 Main Menu Functions in the Software package MAMPRO 33 3.2 Image Restoration Functions in the Software package MAMPRO 34 3.3 Other Image Processing Functions in the Software package MAMPRO 35 4.1 Confusion matrix for Benign (B) vs. Malignant (M) classification by conventional (CMM) and digital (DMM) magnification mammography 4.2 40 Confusion matrix for Benign (B) vs. Malignant (M) classification by conven tional (CMM) and digital (DMM) magnification mammography; ground truth is radiologists’ majority opinion 4.3 41 Inter-observer (dis)agreements for Benign vs. Malignant classification in conven tional magnification mammography 4.4 41 Comparison of classification by conventional (CMM) and digital (DMM) magni fication mammography for each observer 4.5 42 Comparison of image features by conventional (CMM) and digital (DMM) mag nification mammography for the first observer 43 5.1 Comparison of restoration filter performances 56 6.1 Filter performance factors 83 8.1 Classification matrix for algorithm 1 133 8.2 Classification matrix for algorithm 2 134 8.3 Classification matrix for Perceptron 135 9.1 Discriminant Power of Features 144 9.2 Classification performance of computer algorithms on training images . vii’ 146 9.3 Jackknife classification performance of computer algorithms A.1 Radiographic characteristics of breast masses ix 146 183 List of Figures 3.1 The image acquisition system 15 3.2 The image acquisition hardware block diagram based on a linear CCD 16 3.3 Spatial distribution of illumination in the light box 18 3.4 The image acquisition hardware block diagram based on a 2-D CCD 20 3.5 System Contrast Transfer Function 22 3.6 System Modulation Transfer Function (a) in two dimensions w, and w; (b) a one dimensional slice along w, 3.7 24 Modulation Transfer Function of the sampling aperture of the three digitizers compared with the MTF of the screen-film combination: (a) 2-D CCD; (b) Linear CCD; (c) Matrix Laser Scanner in high resolution mode; (d) Kodak Min-R Ortho M screen-film 3.8 27 Noise standard deviation for (a) the screen-film scanned by a 70jtm aperture; (b) Matrix Laser scanner; (c) One frame of the 2-D CCD; (d) Average of 40 frames of the 2-D CCD with correction for fixed pattern noise 28 5.1 The digital image formation model 49 5.2 The Modulation Transfer Function of (a) the system; and the restoration filter with (b) noise to signal power ratio a = 0.01; and (c) a = 0.001; (d) and (e) show the improvement in MTF due to the application of the restoration filters Hr corresponding to (b) and (c) 5.3 Restoration of test images: f is 53 the latent image, fb is the blurred image, g and 1 are the observed images; all other images are the results of restoration by g 57 various filters x 5.4 A portion of a standard test phantom imaged at a far distance near the res olution limit of the camera. Results of restoration by the two Wiener filters: Wienerl with NSR=.01 & Wiener2 with NSR=.001, and enhancement by 3x3 convolutional masks sharpl & sharp2 as defined in the text 5.5 59 Restoration of a mammogram: (a) The observed image; (b) Restored image using Wiener filter (with NSR=.001) compensating for the camera MTF; (c) Restored image using Wiener filter (with NSR=.001) compensating for the combined sys tem MTF 60 6.1 The radiographic image formation model 64 6.2 The Modulation Transfer Functions 68 6.3 The digitizing camera image formation model 70 6.4 A model of radiographic image formation 71 6.5 A simplified model of radiographic image formation 73 6.6 Generation of the noise source n 2 75 6.7 Generation of the noise source n 1 76 6.8 Adaptive smoothing and deconvolution of a digitized mammogram (a) observed image, (b) smoothed by LLMMSE filter, (c) smoothed by MAP filter, (d,e,f), result of Wiener deconvolution of images (a,b,c), respectively 84 7.1 The image degradation and restoration model in the absence of noise 87 7.2 Modulation Transfer Function of the digitizing camera 90 7.3 Restoration of a step edge acquired by a digitizing photographic camera: a) The observed image of the step edge; b) The restored image using linear deconvolu tion, notice the boundary truncation artifact; c)The restored image of the step edge, a trapezoidal was applied to the data prior to its deconvolution; d)The 7.4 restored image of the step edge using the proposed approach 91 Extension of the image by mirroring in the x and y directions 101 xi 7.5 The extended signal under the conditions necessary for exact recovery 7.6 Spatial support of the boundary artifacts around the margins of the image for 105 the case of L,. <M 7.7 106 Restoration of the blurred image of Bayan: a) original; b) a 64x64 section cut from the original larger image; c) a 64x64 section of the original image blurred by the out-of-focus model of the digitizing camera; d) result of restoring by linear deconvolution; e) result of restoring by linear deconvolution after windowing the data by a trapezoidal window; f) result of restoring using image extension and circular deconvolution; g) the difference image of b) and c); h) the difference image of b) and d); i) the difference image of b) and f) 109 8.1 Segmentation using the Local Histogram Thresholding algorithm 8.2 The subimage grid displacement 122 8.3 Effects of area and gradient tests on segmentation 123 8.4 Examples of segmentation using the Local Histogram Thresholding algorithm 8.5 Comparison of two segmentation algorithms . . . . . . 121 125 128 A.1 X-ray mammography system 166 A.2 The heel effect 169 A.3 MTF curves for three Kodak screen-film combinations [131] . . . . . 172 A.4 Variation of contrast with photon energy [126j 174 A.5 Screen-film characteristic curves 175 A.6 Noise power spectrum; film noise (solid line); quantum noise (dash line); and total noise [130] 177 A.7 Ratio of scattered to primary radiation as a function of radiation field [128] A.8 Contrast sensitivity of the eye [126] . . 179 185 xli Acknowledgment I wish to express my gratitude to my research supervisor Dr. Rabab K. Ward for her exemplary and continued support, care and guidance. I owe Dr. Branko Palcic, Professor of Pathology and head of Cancer Imaging at British Columbia Cancer Research Center, in whose laboratories .1 spent the last five years, a personal debt of gratitude for his unbounded enthusiasm and foresight throughout the course of this research. His generous financial support is also most gratefully acknowledged. This work would not have been possible without the constant help and encouragement of Dr. Jacqueline Morgan-Parkes, professor of radiology at the British Columbia Cancer Agency. In all the clinical aspects of this work her collaboration has been indispensable. I wish to acknowledge the support of Cancer Imaging, British Columbia Cancer Agency, and Xiffix Technologies Corporation of Vancouver BC. In particular thanks are due to chief engineer Mr. Bruno Jaggi P.Eng. without whøse support the scope of this work would have been much limited. I am indebted to my friends and colleagues Dr. Calumn MacAulay, Mr. Steven Poon P.Eng., Mr. Brian Pontifex P.Eng., Mr. John Fengler P.Eng., and Mr. Daniel Nesbit for helpful discussions and for the use of hardware or software tools that they had designed. Thanks are also due to many other staff of Cancer Imaging and Xiffix for their technical support. The moral support and encouragement of other friends, too many to be mentioned by name here, will be always remembered. Special thanks to my wife Fariba for her patience, love, faith, and helpful ideas, and to my son Bayan, much of whose companionship has been sacrificed for the completion of this study; to them I dedicate this work. This study was supported financially by the Canadian Commonwealth Scholarship, the University Graduate Fellowship, the British Columbia Cancer Agency and Xillix Technologies Corporation. xl” Chapter 1 Introduction 1.1 Motivation For women in the developed countries, breast cancer is the most frequently diagnosed cancer and a leading cause of cancer deaths. Over 10 percent of all women in North America will develop breast cancer during their lifetime [1]. In Canada [2] breast cancer accounts for approximately 30% of all new cancer cases and for 20% of cancer deaths in women. In British Columbia the incidence rate is about 105 women in 100,000 with a mortality rate of 32 in 100,000 [3]. Breast cancer therefore is a major health threat to women of epidemic proportions. The probability of success in curing breast cancer is directly related to the stage at which it is detected. The earlier the detection, the higher are the chances of a successful treatment. Breast preserving surgical therapy is only possible at the early stages of the disease. There are several techniques for detection of asymptomatic breast cancer. These include regular breast self-examination, clinical examination by a physician and imaging of the breast. Clinical examination of the breast can detect tumors, but not infrequently bigger than 1 or 2 centimeters in diameter. They may have spread to the axillary lymph nodes and require systemic as well as local treatment and more importantly the patients may not be cured. Imaging the breast has therefore emerged as an indispensable tool for early detection. Several types of breast imaging has been investigated by researchers. These include ther mography, diaphanography, ultrasonography and mammography. Thermography is based on the hypothesis that cancerous cells are more active than normal cells. A heat sensing device is used to detect “hot spots”. The clinical reliability of this method is very low and therefore it is rarely used. 1 Chapter 1. Introduction 2 In diaphanography the breast is transilluminated with strong light in the visible spectrum. This method is generally thought to be of limited usefulness. Ultrasonic imaging technology is emerging as a useful tool though as yet it lacks the nec essary resolution. It is most effective in differentiating cystic from solid masses. Magnetic resonance, impedance measurements and several other imaging modalities are also currently under investigation. The X-ray mammography holds the greatest promise of effective early detection of breast cancer. Special low dose X-ray equipment has been developed reducing the danger of radiation induced cancer to negligible levels. It has been well documented that by screening postmenopausal women using X-ray mam mography, the mortality rate can be reduced significantly [4]. Well over one quarter of a million women in the USA have participated in a “Breast Cancer Detection Demonstration Project” which illustrated reduced mortality rate for the screened population [4]. Of equal importance is the fact that detecting the cancer at carcinoma in situ stage, when treatment with minimal surgery followed by ionizing radiation is still possible, results in a much higher quality of life for the patient. Mass screening of asymptomatic women however will generate a vast number of mammo grams. This overload is sufficiently large to place a thorough screening program beyond the means of most communities. A helper tool such as a computerized prescreening device which could aid radiologists in the detection and recognition of subtle signs of abnormality would greatly speed up this process and make it more economical. Additionally, digital analysis of mammographic images can assist in more objective and therefore more accurate interpretation and diagnosis. 1.2 Human Female Breast The breast or mammary gland is a modified sweat gland that has the specific function of milk production. An understanding of its basic anatomy, physiology and histology is important in Chapter 1. Introduction 3 the interpretation of mammography. Briefly, the adult breast is composed of three basic structures: the skin, the subcutaneous fat, and the breast tissue, which includes the parenchyma and the stroma. The stroma is composed of fat and other connective tissues. The parenchyma is divided into 15 to 20 segments, each drained by a lactiferous duct. The ducts converge beneath the nipple, with 5 to 10 major ducts draining into the nipple. Each duct drains a lobe composed of 20 to 40 lobules. The composition of the breast changes with life cycle and is a function of both the diet and hormonal variation [5]. 1.3 Breast Cancer Although the origin and causes of breast cancer are not yet known, it has been established that early detection is the most effective strategy in its management. There are over eighty histologically different diseases of the breast covering the spectrum from mild benign cases to highly aggressive invasive malignancies. There are great differences in appearance of these abnormalities on a mammogram, with many similarities between the benign and malignant cases. Additionally mammographic appearance of normal breasts are quite variable. Many small breast cancers are non-palpable, and some are also mammographically occult. Differential diagnosis of breast cancer is therefore both important and difficult. 1.4 Mammographic Presentation X-ray mammography is used for two purposes: as a diagnostic tool, and as a mass screening method. In screening centers mammography is used to identify suspicious cases, while as a diagnostic tool it is used for the detection of abnormalities, recognition of malignant tumors and preoperative localization of such tumors The increased use of mammography has resulted in an increased incidence of detection of smaller lesions. In British Columbia, for example, at screening, 15-20% of all detected cancers are found in the carcinoma in situ stage. The normal presentation before mammography screening was 2-3%. Early detection of pathological tissues depends on the mammographic Chapter 1. Introduction 4 image quality. Radiologists typically scan mammograms for any signs of abnormality which may include 1) clusters of microcalcifications, 2) well-defined and rounded or ill-defined and spiculated masses, and 3) other signs such as skin thickening or architectural distortion of the breast etc. Ra diologists usually compare a mammogram with other views of the same breast, those of the other breast or previous mammograms of the patient to check for asymmetries or developing densities. The radiographic images are quite complex. Despite a highly evolved human visual system a radiologist requires many years of training to detect subtle abnormalities in a complex parenchymal pattern. To aid radiologists in this complex task researchers have reported different image and vision processing techniques. A summary of this work is given in chapter 2. Despite these efforts there is much room for improvement and this remains an active field of research. Computer vision methods have in fact been used successfully in another area of cancer imaging, namely that of cell classification in cervical cancer screening programs [6]. In digital mammography the efforts in applications of digital image processing have so far concentrated on image enhancement for visual interpretation by a radiologist. Image analysis for detection and classification have also been reported. Two important aspects however have not received attention. The first is that in the process of X-ray image formation the latent image suffers degradation due to the system blur and noise. The process of digitization of the image further contributes to both the blur and the noise. The first step should therefore be the restoration of the image from these effects. The second aspect is the need for accurate and rapid digitization of the film images if any image processing algorithm is to be used in a clinical setting. In this work these two issues are addressed before application of pattern recognition techniques for detection and recognition of mammographic abnormalities. Chapter 1. Introduction L5 5 Objectives The objective of this study is to develop a system to assist in the detection and recognition of microcalcifications which are frequently present in a mammogram. The most subtle signs of abnormality are not visible in the digitized mammograms the original films — and may not be visible in because of the noise and blur in the image. These are due to the X-ray machine, the camera and the digitizing equipment. If the characteristics of the noise and the different blur functions are known or obtained then these effects may be removed or reduced by image restoration techniques. Different image enhancement techniques have been reported in the literature but image restoration has not yet been applied on mammographic images. While image enhancement has its own advantages such as removing noise and highlighting some features, it does not compensate for the blur. Thus it would be of great value to restore the image to aid radiologists in their analysis, prior to applying the automatic edge detection, feature extraction and classification routines. Following image acquisition and restoration, segmentation and recognition of microcalcifi cation clusters will be developed and evaluated. 1.6 The Structure of This Thesis Chapter 2 contains a review of the state of the art in processing and analysis of the digitized mammograms. In chapter 3, the hardware and software development of a system for the digitization of mammograms using a Charge Coupled Device (CCD) sensor are described. The acquisition system characteristics are reported in terms of its spatial and photometric response. The Modulation Transfer Function (MTF) of the system is calculated from its Square Wave Response Function (SWRF’) and also from its Edge Spread Function (ESF) in both spatial and frequency domains. It is shown that the combined effects of fixed pattern and random noise can be reduced to within round-off noise, having variance of 0.25 gray levels over the whole image. Chapter 1. Introduction 6 In chapter 4, the clinical evaluation of the film digitization system is reported. It will be shown that in most cases radiologists can extract essentially the same diagnostic information from electronically magnified images as from the magnification mammography procedure. In chapter 5, the effects of image restoration on the digitized mammograms are considered. Two simple Wiener image restoration filters which counter the high frequency attenuation due to the MTF and sharpen the mammographic image details are tested. The results of these procedures are the removal, to a great extent, of the blurring effects of the different components of the image formation system and the generated noise. It has been postulated that this image restoration will aid in the quantitative analysis of mammographic images and assist the radiologists in improving their diagnosis. Chapter 6 extends the restoration of mammograms to include the effects of radiographic signal-dependent noise. To accomplish this a comprehensive image formation model is pre sented. In chapter 7, a novel method of carrying out filtering in the frequency domain using image extension and circular convolution are presented. This approach eliminates much of the bound ary artifacts associated with linear, frequency domain restoration of truncated images. we also present a mathematical analysis of this technique. Chapter 8 presents algorithms for automated segmentation of microcalcifications from the background parenchyma of the breast. Chapter 9 describes the extraction of over 100 quantitative features from the segmented clusters of microcalcifications and their application in classification of mammographic abnor malities. Finally, chapter 10 presents conclusions and suggestions for further research. The appendix gives many details and physical principles of mammographic image formation and interpretation. Chapter 1. Introduction 1.7 7 Claims to originality Materials presented in chapter 2 are gleaned from the relevant published literature. CCD cameras have been used by others to digitize X-ray films. The system developed in chapter 3 however is novel in construction and in many of its characteristics as described in the body of this thesis. The software package MAMPRO is a new research tool that contains all of the algorithms developed in this work. The clinical investigations reported in chapter 4 are novel. The restoration filters employed in chapter 5 are well known, but their application to digi tized mammographic images is new. The signal-dependent restoration filters described in chapter 6 are known. Their derivation for the radiographic noise is new, as are the radiographic image formation models. The image extension technique presented in chapter 7 has been reported in the literature in the context of spatial domain filtering. The analysis of its effects in frequency domain filtering for restoration problems, however, is new. The first segmentation algorithm described in chapter 8 is taken from the literature. The next two algorithms were developed by me for this work. Most of the features used in chapter 9 are reported in literature in other contexts. Evaluation of their utility in classification of mammograms is new. The data base of digitized mammograms and the resulting multivariate discriminant functions are also new. At present, no practical (commercial) system exists which provides means to radiologists for accurate viewing, interactive manipulation, or objective (automated and quantitative) analysis of mammograms. The radiological workstation described here forms the basis for such a system. It provides a helper tool (and a second opinion) to aid radiologists in their task of primary diagnosis of mammograms. Chapter 2 State of the Art in Digital Mammogram Processing and Analysis Most early efforts in digital analysis of mammograms were concentrated on processing of Xero radiographs, poor quality radiographs, or have relied on manual measurement of features. Since screen-film mammography has only recently produced acceptable radiographic image quality, we will give a partial list of published work in this area since 1980. This is a very brief chapter that merely points to the relevant literature in the field; details of the published studies that have a direct bearing on this work are given in the subsequent chapters. Different attempts to digitally process mammograms have had different goals in mind. These can be divided into the following four groups: 2.1 Image enhancement Digital image enhancement techniques either employ global manipulation of grey levels or locally adapt such manipulations to image features. The input image is modified using a set of usually heuristic rules to enhance the visibility of certain desired features [7]. In one approach a linear combination of smoothed images and a non-linear contrast transformation is used to obtain enhancement [8]. Alternatively a global estimate of the background breast structure is employed to bring out pathologic abnormalities [9]. Image neighborhoods that are locally adapted to the spatial extent of image features are selected. In this way enhancement techniques such as contrast manipulation, histogram equalization, etc. respond to the local image detail [10, 11, 12, 13, 15]. Finally a non-linear mapping is used to encode the image gray levels in an attempt to equalize the system noise which is signal-dependent [14]. In other approaches to smoothing of mammographic images adaptive order statistic filters 8 Chapter 2. State of the Art in Digital Mammogram Processing and Analysis 9 and tree-structured wavelet transforms are employed [16, 17]. These filters operate on the mammogram at multiple resolutions and therefore are potentially useful in preserving image features such as edges while smoothing the image [18]. 2.2 Risk assessment It has been suggested that a woman’s risk of developing breast cancer can be determined by the pattern of parenchymal densities on the mammogram [133]. The goal of computer aided risk assessment is to correlate the mammographic parenchymal pattern with the risk of developing breast cancer using objective and repeatable measures commonly based on texture [19, 20, 21] or density [22, 23]. 2.3 Automated detection of abnormalities Two types of abnormalities have been investigated namely microcalcifications and masses. Most detection schemes enhance the conspicuity of abnormalities as an intermediate stage, before selecting candidate pixels belonging to abnormalities. The classical method of unsharp masking for the detection of abnormalities is evaluated in [24]. To detect the presence of microcalcifications, linear filtering techniques have been employed. Researchers at the University of Chicago have used matched filtering to enhance the signal, while a box-rim filter is employed to suppress the signal. The presence of microcalcifications is detected from the difference of these signal enhanced and signal suppressed images [25, 26, 27, 28, 29]. Neural network techniques are then employed to reduce the number of false positive detections [30, 31, 32]. Many other techniques have been used for the extraction of microcalcification images from the background. These include: i) local area thresholding [35, 36], ii) morphological opera tions [37, 38, 39], iii) stochastic image models [40, 41], iv) features derived from contour plots for signal peak detection [42, 58], v) multiresolution approaches [44, 45], and vi) wavelet trans forms [46, 47]. A battery of tests involving local contrast, shape, size, gradient, and proximity Chapter 2. State of the Art in Digital Mammogram Processing and Analysis 10 to other microcalcifications have been commonly used to reduce false positive detections of clusters of microcalcifications [33, 34]. Research into the detection and segmentation of masses has also been reported in the lit erature. Asymmetry between the right and left breasts have been exploited to detect possible masses [49, 50, 51, 52]. Template matching has been applied for the detection and segmentation of circumscribed masses [48]. The concept of multiresolution image processing based on fuzzy pyramid linking is used to detect and subsequently classify masses [53, 54]. Scale space filtering is used to extract closed contours associated with mass boundaries [55]. Texture features have also been used to detect stellate lesions [56, 57], and Markov Random Field image model has been utilized to segment tumors [66]. 2.4 Differential Diagnosis The aim of differential diagnosis is to apply pattern recognition methods to differentiate benign and malignant lesions, or to separate benign or malignant subgroups. In one study radiographs of biopsy specimen were used to extract features and calculate a discriminant function to classify clusters of microcalcifications [59]. Different sets of features were used to approach the same task for screen-film mammograms [60, 61, 62]. In a different approach to the same problem, individual objects were not extracted from the background. Instead structural features were computed from the image. Subsequently artificial neural networks were applied to these feature vectors to classify the whole mammogram without the need for the prior segmentation of microcaicifications [67, 68, 69]. Shape features have been calculated for the classification of masses [63, 64, 65, 66]. Fractal dimension has also been used to classify lesions based on image texture [70]. Other related developments are in the areas of film digitization [71, 72] and development of smart workstations for computer aided diagnosis [73, 74, 75, 76, 77]. The most recent advances in almost all of the above areas are reported in the proceedings of the two international workshops that have so far been held in this field [159, 160, 161]. Chapter 3 Radiographic Film Digitization Accurate digitization of the film is the first essential step to a successful automated analy sis. The quality of digitized images is of fundamental importance as any information lost in the acquisition and digitization process can not be recovered by further processing. Accurate interpretation of digitized X-ray screen-film mammograms has been limited by the lack of a ‘specialized’ image acquisition system with high speed and accuracy. Several researchers, e.g. [19], attempted to extract numerical feature descriptions for tumor classification. They report however that a major llmitation on the success of these efforts is the need for accurately digitized images. The mammographic features of interest range from several centimeters across, as in architectural distortions of the breast, to fine microcalcifications of less than 0.1 mm across. Very high resolutions are therefore required for at least some portions of the image. Digitizing the entire mammogram at the highest possible resolution requires a prohibitive amount of memory. Thus my approach is to develop an acquisition system that facilitates the digitization of mammogram images at the various scales needed for the different stages of analysis. The digitization should be rapid and within the required accuracy. Such a system would be useful for both stages of pre-screening in the mass screening of non-symptomatic women as well as for the more detailed analysis of suspicious lesions. Two novel image acquisition systems based on a linear and an area scanning scientific grade Charge Coupled Device (CCD) arrays were developed. Both systems will be described in this chapter and their performance will be compared. Although both systems are useful, the second system is superior and meets the specified requirements in that it offers a fast method of digitizing mammograms with high spatial and photometric resolutions. The system is capable 11 Chapter 3. Radiographic Film Digitization 12 of acquiring 6 frames per second where each frame consists of over 1.3 miffion pixels digitized to 12 bits per pixel. The fixed pattern and the random noise (of optical and electronic origin) are minimized using background subtraction and signal averaging techniques. The quality of the resulting image in terms of sharpness and noise content is comparable with that obtained by the more expensive and slower drum laser-scanning microdensitometer. 3.1 Image Acquisition Hardware There are three methods of film digitization, namely using: 1) microdensitometers; 2) video cameras; 3) one-dimensional and two-dimensional CCD digital cameras. 3.1.1 Microdensitometers X-ray film digitization has traditionally been performed using a microdensitometer. A nar row beam of light is transmitted through the film, and a light sensitive device measures the transmittance which is then converted to optical density using a calibration curve. Lighting conditions approximate specular reflection. The film is moved past the beam by either a rotat ing drum or a flat bed mechanism. The sampling aperture of these systems is limited to the spot size. The scan-time is proportional to the number of sampling points or sampling lines depending on the scan mechanism, and may be quite long. Furthermore the calibration and operation of the system is an involved process requiring a skilled operator. For these reasons microdensitometers have not been widely used to digitize radiograplis in clinical settings. X-ray film digitization using laser scanning microdensitometer has been widely reported in the literature. A comparative study of digitized film radiography systems and the associated diagnostic and operational advantages are given in references [78, 79]. The performance charac teristics of laser digitizers are reported in [80, 81, 82, 83]. A recent assessment of these systems for clinical applications concludes that they are generally slow and expensive [84]. Chapter 3. Radiographic Film Digitization 3.1.2 13 Video Cameras Video cameras have dynamic ranges and band width limitations that affect their performance at the spatial and photometric resolutions required. As an example the Panasonic WV-BD400 video camera was evaluated. The video signal was captured and digitized to 512 x 480 pixels by the frame grabber. The contrast transfer function of this camera falls off to 5% at 400 TV lines on standard video resolution charts when the chart-camera distance is set for a full field of view. This corresponds to a spatial resolution of 1.25 line pairs per mm. This is clearly inadequate for detection of submilimeter microcalcifications. Also the useful gray scale resolution is only about 7 bits [851 and the camera response as a function of gray scale is non-linear with non-unity gain factor. The analogue output of the camera has to be digitized resulting in further quantization errors and the integration times are restricted to be less than 1/30 seconds for interlaced frames at 60 Hz mains frequency. Therefore, video cameras were not considered any further. 3.1.3 CCD Sensors Two new image acquisition systems which take advantage of the recent advances in Charge Coupled Device (CCD) technology were developed. The first system is based on a linear array CCD, and the second system uses a two-dimensional CCD array. Linear CCIJ The first system uses a digital scanner (Datacopy 612, Mountain View, CA) which has a linear array solid state detector (Fairchild CCD 122). There are 1728 pixels spaced contiguously every 13 m over a 22.5 mm length. Two-dimensional images can be acquired by moving the sensor in 13 im steps across the image. Scanning 2846 lines per image results in image files of 4.9 Megabytes where each pixel is digitized to 8 bits. The amount of optical aberration at the corners however is severe and therefore the image was limited to 2500 lines. The integration time is fixed at 3.5 ms. The camera output is fed to the Datacopy image processing interface board model 110 Chapter 3. Radiographic Film Digitization 14 housed in an IBM AT computer. The camera interface board communicates with the frame grabber (FG100AT, Imaging Technology Inc., Woburn, MA) and a digital signal processing board (SKY32O, SKY computer Inc., Boston, MA) via the host computer bus (IBM-AT). The arrangement is shown in Figures 3.1 and 3.2. The digital signal processing board based on TMS32O signal processing chip is used to speed up the scanning operation. The resulting image file is written onto the disk. The frame grabber board is used to process and display the image on a Sony RGB monitor. Since this board can only handle 512 x 512 image files the input image was subsampled prior to display. Six mammograms were digitized using this system. Visual inspection of these images show them to be of inferior quality compared with images obtained from a microdensitometer. This system provides a large field of view but the scanning mechanism imposes two limit a tions on its performance. The first limitation relates to the illumination source and is discussed below. The second performance limitation of this system is the scan time. This speed limitation is characteristic of all point and line scanning mechanisms such as drum laser scanners and linear array CCD cameras. Nonetheless the system is typically at least five times faster than a microdensitometer and is considerably less expensive. Several other researchers have also reported the use of a linear array CCII camera [13, 14]. As discussed later, the random noise contributed by the illumination source and the sensor electronics may be reduced by averaging several frames. This approach however is impractical for this system due to the scan time which is several minutes per frame. To further improve the speed and accuracy of the system the second acquisition system was built. Illumination The camera was mounted above a light box transilluminating a mammogram. This light box is a common viewing box used by radiologists. It has two a.c. powered fluorescent lamps under a light diffuser. Images obtained by the linear CCD displayed a regular pattern of light and ! Mammo gram Fig. 3.1. The Image Acquisition System Mouse Stabilize Illumination Source ‘crocomputer with Interface and Display Peripherals High Resolution Display Camera Stand High Resolution CCD Sensor Cfl, -mm w Keyboard rwwn Computer Monitor THE IMAGE ACQUISITION SXSTEM QT1 CTI I. C.) T Buffer F I ii Green LJBlueLr1: DAC FI1 I I A JIILLTj FG100 Frame Grabber Fig. 3.2. The image acquisition hardware based on ci linear COD I I I Ti I I F]UflDAJI I[U! r-’ Frame U I: rIMUXLUTI I THE IMAGE ACQUISITION HARDWARE Chapter 3. Radiographic Film Digitization 17 dark bands parallel to the CCD array. This artifact was interpreted as being due to light flicker produced by the a.c. source. The 60 Hz ripple in the power supply means that the illumination flux varies over time. Consequently, different lines of the image will have different intensities. For the case of the linear CCD, the camera’s fixed integration time of 3.5 ms is not long enough to overcome this problem. This artifact is sufficiently disturbing as to render the image essentially unusable for quantitative work. This artifact may be minimized by use of a high frequency ballast. The standard ballast, having output frequency of 60Hz was replaced by a Plaser ballast of 600Hz. The disturbing artifact was reduced but was still perceptible. Therefore a new light box (Gordon Instruments TVS, Orchard Park, NY) was acquired with four d.c. driven, 400 watts, quartz lamps. This unit is air cooled and can transilluminate transparencies of sizes up to 14” x 17”. The spatial non-uniformity of light emitted through the top surface of the light box was measured using an optical power meter. A radiometric filter was coupled to the light power detector to make broad-band measurements in the range of 450-950 nm wavelengths. The incident light power is given in Figure 3.3 as a function of distance from the centre. It can be seen that a hot ring exists about 10 cm from the center. The maximum variation is about 20%. Narrow band radiometric and photometric measurements were also performed which confirmed the existence of hot spots. The light non-uniformity was also measured using the CCD camera. the background illumi nation has a variance of 50 gray levels at mid range on an 8-bit scale. However this ‘hot spot’ non-uniformity can be accounted for as it is of fixed pattern. The fluorescence lighting provides a more uniform illumination (variance of 10 on an 8-bit scale). It is also cooler and therefore does not require the fan associated with the d.c. driven source. The problems of flicker and constant (i.e. non-adjustable) output light can be solved by use of variable integration times and area scanning CCD arrays. Chapter 3. Radiographic Film Digitization 18 7 6 I 5 4 3 2 1 0 0 20 30 Distance (Centimeters) Figure 3.3: Spatial distribution of illumination in the light box 40 Chapter 3. Radiographic Film Digitization 19 2-D CCD The second system is an adaptation of the detector originally developed for a solid state microscope (SSM) employing Micro-Imager 1400 (XilJix Technologies Corp., Vancouver, BC, Canada) [86, 87]. This detector uses a two-dimensional CCD array (Kodak KAF-1400) to cap ture 1320 x 103b pixels with square sensor elements of 6.8 im per side. Up to six frames per second can be acquired with variable integration times. The CCD output is digitized to 12 bits, 8 of which can be displayed at any one time. The contents of frame memory is displayed on a high resolution monitor (1280x1024 pixels) and the system work station (Apollo DN4500) has access to the frame memory for image transfer and subsequent analysis. The system block dia gram is given in Figure 3.4. A highly interactive C program with an easy to operate graphical user interface was developed for mammographic image analysis. The following is the description of the performance characteristics of this system. The performance characteristics of the light source, the lens, and the transfer functions of the camera are described in section 3.2. Noise reduction is discussed in section 3.3 and a comparison of film digitization systems is presented in section 3.4. 3.2 Performance Characteristics The following six parameters characterize the quality of an image acquisition system [88] spatial resolution, photometric resolution, photometric linearity and spectral response of the transducer, spatial distortion due to non uniformity of illumination and size/shape variations of transducer pixels, temporal distortion due to illumination fluctuations and electronic noise. The photometric range of the detector is measured to be between 0 and 3.2 Optical Densities. This dynamic range is sufficient to cover all radiographic densities present in a properly exposed and developed mammogram. The CCD response was measured as a function of the intensity of incident light and was found to be linear. Since the CCD measures the optical transmittance of the mammogram at each pixel location, a logarithmic transformation was programmed in the acquisition system’s look-up table (LUT) to enable the display of the image in the optical Chapter 3. Radiographic Film Digitization 20 rj) j) z C 0 C) 0 C 0 -D U) (6, U -n E 0) a -U - 0 0 -Q a 0 -I C,, a 0 U U) U) U E U) -c I CV) U) LL Chapter 3. Radiographic Film Digitization 21 density domain. Using a Nikon 55 mm f/2.8 Micro Nikkor lens and the 2-11 CCD sensor it is possible to digitize a mammogram with a minimum sampling interval of 12.6 um. This high density sampling is required for the task of detection of minute microcalcifications which may be present in the vicinity of the primary tumor site. These “satellite” microcalcifications are indicative of histologic multifocality and are believed to be of prognostic value [89]. The Contrast Transfer Function (CTF) of the camera using the above lens is measured and given in Fig. 3.5. The limit of resolution was measured to be 770 TV lines when imaging a standard resolution chart. Although the CTF is commonly used to specify the spatial resolution of a camera, it is not in general a linear measure and therefore the overall CTF can not be readily computed from the CTF of cascaded elements. Therefore the Modulation Transfer Function (MTF) which is the Fourier Transform of the Point Spread Function (i.e. blur) was used. Two methods were used in measuring the MTF of the camera 1) the Square Wave Response Function and 2) the Edge Spread Function. The Square Wave Response Function (SWRF) has been used by others [90] and can be directly measured using bar pattern test objects. Each sinusoidal component of the square wave with frequency nf will have its amplitude MTF(nf) so that: S(f) > = MTF(nf) (3.1) n1 3 and we can compute the MTF from MTF(f) = [S(f) + S(3f) - S(5f) + (3.2) where S is the SWRF. From equation (3.2) it can be seen that the values of the MTF at high frequencies are dependent on the accuracy of the SWRF measurements at much higher multiples of these frequencies. At these higher frequencies the SWRF values are extremely difficult to measure Chapter 3. Radiographic Film Digitization 22 1 .0 0.9 0.8 0.6 U H 0 0.5 0.4 0.3 0.2 0.1 0.0 0 10 20 30 40 Frequency (Cyc’es/mm) Figure 3.5: System Contrast Transfer Function 50 60 Chapter 3. Radiographic Film Digitization 23 since they are limited by noise. Thus the results of this method were not further employed in this work. In the second method a large black test object with a straight edge is imaged. The Edge Spread Function (ESF) is obtained by averaging 128 one-dimensional step edges. The Line Spread Function (LSF) is obtained by differentiating the ESF. Finally one slice of the MTF is obtained by the one-dimensional discrete Fourier Transform of the LSF [91]. The differentiation may be done in space domain (x, y) or equivalently in the spatial frequency domain (u, v): MTF(u, v) 6(v) . = 2?rju G(u, v) (3.3) where G(u, v) is the Fourier Transform of ESF, and 6 is the unit impulse function. The twodimensional MTF is a separable function formed by the product of one-dimensional MTFs in the u and v directions. Furthermore the Point Spread Function (PSF) of the camera is found by a two-dimensional inverse transform of the MTF. A plot of the MTF is given in Figure 3.6. The results of this method was selected to be further used in this work. Noise Reduction 3.3 The aim here is to correct the noise introduced in the system. Specifically we are concerned with the noise element added by the film digitizing process. This noise is due to the following four sources: a) the illumination source; b) the lens; c) the CCD and its associated amplifier; and d) the analog to digital converter. We combine these four sources and divide the overall noise present in the data into two types: the fixed pattern noise such as the optical shading of the light source and aberrations of the lens, and random noise of optical and electronic origin. The fixed pattern noise can be corrected by image calibration using [92]: ‘cat. = k ‘raw bright where ‘raw is the raw image, ‘dark — ‘dark T T — (3.4) dark 1 is the image with the shutter closed, Ight is the image of the background without the mammogram and k is a positive constant. For rapid processing, Chapter 3. Radiographic Film Digitization 24 C q h 0 0 U- LI) 0 f r eO..— / n—%. fl (a) 1 .0 0.9 0.8 0.7 0.6 U F— 0.5 0.4 0.3 0.2 0.1 0.0 0 10 20 30 40 I I 50 60 • I 70 80 Frequency (Cycles/mm) (b) Figure 3.6: System Modulation Transfer Function (a) in two dimensions w and w; (b) a one dimensional slice along w Chapter 3. Radiographic Film Digitization 25 integer arithmetic can be used. In order to avoid discontinuities created by division by an integer, the following approximate relationship was implemented. cal. = ‘raw T - — ‘bright + bias (3.5) Applying this method to the full image, the noise standard deviation a, was reduced from 3.4 to 1.7 gray levels on the 8-bit scale at mid-range. The random noise can be reduced by averaging multiple frames. If N frames are averaged, the standard deviation of the noise will be reduced by a factor of After subtraction of the fixed pattern noise 25 frames were averaged. This operation takes approximately 10 seconds and reduces the standard deviation of the noise to 0.5 gray levels on the 8-bit scale at mid range. This low level of noise was found to be independent of the integration time or the size of the image and is primarily due to the arithmetic round-off errors introduced in equation (3.5). After correction for both types of noise, a 1.3 megabyte image with a maximum background variation of +2 gray levels was obtained. 3.4 Comparison of the Three Film Digitization Systems A wide range of microdensitometers are in use and their characteristics have been widely re ported in the literature. We have presented the above two film digitizers based on a linear CCD and a two-dimensional CCD. We will now compare the performance of these systems with a typical microdensitometer. We will use a ‘Matrix Laser Digitizer’ as a typical unit of modern design [81]. Pixel Size: Microdensitometers usually provide the option of a few fixed pixel sizes, in the range of 50200 sum. Each pixel is circular and its minimum size is limited by minimum beam spot size. For CCD cameras the pixel is square and therefore congruent pixels cover the image without overlaps. My first system (linear CCD) has 13 um pixels and the second system (2-D CCD) Chapter 3. Radiographic Film Digitization 26 has 6.8 um pixels. Additionally, since the two-dimensional CCD has 100% fill factor there are no gaps in the light sensitive area of the detector. This can not be achieved by the circular scanning spot of microdensitometers. Modulation Transfer Function: The MTF of each digitizer is the combined effect of the MTF of the lens and the MTF of the sensor. All three detectors can be designed such that the result ant MTF of the sensor is close to the sampling aperture function with little degradation introduced by the sensor electronics. We can therefore compare the aperture functions of the three systems. Fig. 3.7 gives the plot of the MTFs of the three sensor apertures as compared with the MTF of a high resolution screen-film combination. We have used the Matrix laser digitizer in the high resolution mode [81] and the Kodak Min-R screen and the Ortho-M film. The effect of finer sampling by the CCD in obtaining higher spatial resolutions is evident. Signal to Noise Ratio: The temporal noise of the digitizer is not a function of the sampling interval and may be measured as a function of the film density. Fig. 3.8 is a plot of digitizer noise as a function of film density compared with screen-film noise. It can be seen that temporal noise in the microdensitometer and the (JCD sensor, after averaging, are comparable and lower than screen film noise for densities less than 2.5. Since the area scanning CCD grabs an image at least 2 orders of magnitude faster than the other two systems, it is the only digitizer that allows noise reduction by averaging of multiple frames within acceptable time scales. After noise reduction, the system noise is primarily due to quantization noise in the analog to digital ( A/D ) converter. Dynamic Range: The useful density range of a film (defined as the density range with sensitometric slope is between 0.2 and 3.2 optical density units. Therefore a dynamic range of > 0.5) in the signal Chapter 3. Radiographic Film Digitization 27 (a) 2—D CCD (b) Linear CCD (c) Laser Scanner (d)Screen—Filrn 0.6-\ H 0.5-i’. ‘% 0.4 (d)\ 0.30.2 0.1 00 ‘ \ ‘\ - - 0 10 I I 20 30 • 40 50 I I 60 70 80 Frequency (Cycles/mm) Figure 3.7: Modulation Transfer Function of the sampling aperture of the three digitizers compared with the MTF of the screen-ifim combination: (a) 2-D CCD; (b) Linear CCD; (c) Matrix Laser Scanner in high resolution mode; (d) Kodak Min-R Ortho-M screen-film. Chapter 3. Radiographic Film Digitization 28 . 0.18 0.16 (c)/ 014 0.12 -r • “o 0.1 0 --- (1) G) ---a--- 0.08 (a) (b) (c) (d) Screen—Film Microdensitometer CCD single frame CCD calibrated if / /1 i / 0.06 (a) 0.5 0.7 0.9 1 .1 1.3 1 .5 1 .7 1 .9 2.1 2.3 2.5 Optical Density Units Figure 3.8: the screen-ifim scanned by a 7Opm aperture; (b) the 2-D CCD; (d) Average of 40 frames of the 2-D CCD Noise standard deviation for (a) Matrix Laser scanner; (c) One frame of with correction for fixed pattern noise. Chapter 3. Radiographic Film Digitization 29 (transmission) space is required. This can be obtained by all the above digitizers and a 12-bit A/D. Image Acquisition Time The clinical application of digital radiology requires processing of hundreds of images per day in an average sized radiology department. Film digitization time is therefore significant and should preferably be less than the film processing time of one second per inch. Two-dimensional CCD scanning has a major advantage in this regard and is between a hundred and a thousand times faster than the other two systems. Summary In summary, both our CCD based systems can give higher resolutions than that of the densit ometer and the two-dimensional CCD based system is significantly faster than the other two systems. It should also be mentioned that a CCD based system could be manufactured at much lower costs than the densitometer. 3.5 Image Acquisition at Multiple Scales It is often necessary to acquire several images at different spatial resolutions from the same mammogram. To facilitate rapid image acquisition the camera was mounted on a belt driven vertical linear drive. The mechanism was coupled to a three stack stepper motor. The lens was also connected to a rotary gear system that was driven by a smaller stepper motor. In this way the camera could provide real time automated optical zoom capabilities. The drive controller for the lens subassembly was programmed with a spline corresponding to the required amount of rotation of the focus ring. In this way the focus ring would automatically maintain focus while the camera moves vertically. Since the mammogram inherently contains a two dimensional image at a fixed location atop the light box, this form of obtaining autofocus is adequate. An x-y stage, controlled via RS232 port from the computer could be added to the Chapter 3. Radiographic Film Digitization 30 system to provide pan capability. In initial clinical use however, practicing radiologists have found that a manual pan of the film under the camera is just as acceptable when the system is being used interactively. This form of zoom and pan provide magnification with increased spatial resolution and is superior to the software zoom and pan provided by pixel replication in commercially available software for imaging boards. Once a region of interest (ROT) had been identified, either through detection of microcal cification clusters or masses, its center coordinates were marked and new images with much higher resolutions were obtained within a radius of up to 2.5 cm. These newly acquired images were subject to various signal processing algorithms in an effort to identify minute calcification formations associated with “satellite tumor sites”, around the central reference tumor. These sateffites are normally mammographically occult to the naked eye but it has been pathologically demonstrated that the presence or absence of these satellite tumors are highly prognostic [891. The question which was addressed is whether the provision of higher spatial and photometric resolution by itself is sufficient to enable the detection of these cancers if, in fact, they exist. For this latter analysis a data base of mammograms is required on which subsequent pathological analysis has shown the existence of satellite cancers. This issue will be further investigated in chapter 6. 3.6 The Software Environment None of the currently available commercial image processing software libraries are of direct use since they are for general purpose use and do not conform to the particular requirements of mammogram analysis such as image size, etc. These systems have not been designed for object detection or pattern recognition work, are of a “turn key” style and generally unsuitable for development work. A software environment therefore was developed to facilitate the acquisition, storage and retrieval, processing and display of images of various size and format. Development work was done on the following three hardware platforms as they became available to me during the course of this project: Chapter 3. Radiographic Film Digitization 31 1. an IBMPC compatible computer with 80386 CPU and a FG100AT Imaging Technology monochrome imaging board. 2. an IBMPC compatible computer with 80486 CPU and a 1280 Matrox colour imaging board 3. an Apollo 4500 computer with a Univision monochrome frame grabber. The C programming language was chosen for all algorithm development. It is a medium level language that combines ease of programming of high level languages with the speed of machine specific languages. For the PC environment, elementary routines were developed for manipulating image files. The visual display unit available was a 480 x 640 pixel VGA monitor. Therefore routines were developed to select subimages of 320 x 320 pixels from the input images. Also, routines for subsampling the original image were written to create smaller image files at lower resolutions. With a subsampling rate of 1:2 in each direction, it takes only 30% more disk space to store all the low resolution ancestors of a given image in a pyramid linked data structure. Routines were also written to enable writing of image data to the video memory for display. Two images of 320 x 320 can be simultaneously displayed on the monitor. This enables rapid visual comparison of input and processed images. Routines were also written for calculating the global histogram and for applying a user defined global threshold to the image gray level. When the Matrox imaging board became available a software package named WINMAM was developed for the acquisition and display of images. This set of routines use the windows graphical interface and a secondary large screen (1280 x 1024 pixels) is then used for image display. 3.6.1 MAMPRO The C programming language operating under UNIX was used in software development for the Apollo workstation. A user friendly, menu driven graphical interface was employed in Chapter 3. Radiographic Film Digitization 32 this system and the routines are callable from this interface. The resulting package is called MAMPRO and a listing of its menu structure developed for this study is given in tables 3.1 to 3.3. 3.7 Mammogram Data Base An image data base was accumulated from mammograms obtained from the British Columbia Cancer Agency. To minimize the effects of image degradation associated with the X-ray mam mographic unit, all images were collected from similar units. The size of the data base was limited by the availability of films with related radiologic and pathologic reports and the digital storage facilities, however over 500 mammograms were examined and over 200 digitized images were obtained. Each image is 1280 x 1024 pixels in size. Most mammograms were digitized at either 50 m or 100 m at 12 bits using our 2D CCD sensor. Further details are given in the following chapters. Additionally, the image data base consists of five images of 1728 x 2500 pixels obtained from the Datacopy linear scanner and 39 images of 1024 x 1024 pixels from the Royal Marsden Hospital in London, U.K., digitized on a drum laser scanner. Chapter 3. Radiographic Film Digitization [ Main Menu Camera Control Functions Set Input Look Up Table (LUT) Set Integration Time Shutter On/Off Screen Facilities Set Display LUT Set Pseudo Color Graphics / Text Overlay Zoom / Pan Image Acquisition Scan Average Subtract Background Image Storage/Retrieval Main Memory / Frame Memory File Format Image Statistics Intensity Profile Linear / Log. Histogram Image Processing Image Restoration Image Enhancement Edge Detection Morphological Processing Image Arithmetic Images Arithmetic Image Segmentation Features Extraction Label Calcifications Calculate Object Features Calculate Cluster Features Table 3.1: Main Menu Functions in the Software package IVIAMPRO 33 Chapter 3. Radiographic Film Digitization Add Noise Clean Noise Add Blur Form Pyramid Power Spectrum Restore Camera Restore Mammogram 34 Add Gaussian noise Add Poisson noise Add photonic noise Add radiographic noise Wiener smoother MMSE filter MAP filter Add uniform blur Add Gaussian blur Add camera blur Add radiographic blur Subsample by /2 Image extension Spatial / Freq. domain Inverse filter Wiener filter CbS filter Adaptive / Signal-dependent Table 3.2: Image Restoration Functions in the Software package MAMPRO _____________________________ Chapter 3. Radiographic Film Digitization Image Enhancement Edge Detection Morphological Proc. Image Arithmetic with a constant 35 Image convolution Rank filters Average filter Sharp filters Unsharp mask filters Enhance contrast Equalize histogram Optical Density Sobel filter Roberts filter Kirsch filter Prewitt filter Laplacian filter Morphological edges Marr & Hildreth Dilation Erosion Closing Opening Morphological edges Detect salts Detect peppers Skeletonizing Image inversion Image threshold Image + x / and or xor Image clip clamp Images ± x / and or xor images mini/max Images rms error Set threshold Increase SNR Segment calcifications Label objects Create binary masks - Images Arithmetic between two images Image Segmentation - Table 3.3: Other Image Processing Functions in the Software package MAMPRO Chapter 4 Clinical Evaluation of the Film Digitization System A prototype unit of the film digitization system was evaluated subjectively from a clinical point of view by experienced radiologists from the BC Cancer Agency. As a result of their initial feedback it was postulated that such a film digitization unit will be of clinical utility even without any further image processing and analysis software. This chapter reports on a clinical investigation that was designed and implemented to assess this aspect of the device. Radiologists normally scan a mammogram for a number of abnormalities including the presence of any microcalcification clusters. Often a hand-held magnifying glass is used to ensure that the very small and faint microcalcifications are detected. If a definite assessment still cannot be made the patient is recalled and a magnification mammogram is taken. Magnification mammography is a well established conventional procedure that is used as a diagnostic tool in evaluation of microcalcifications. It achieves 1.5 to 2 times magnification of a selected portion of the breast image by introducing an air gap between the breast and the screen-film receptor. The air gap also increases the unsharpness of the image and therefore larger magnification ratios are not practical. A superior image to the conventional mammography is obtained through exposing the patient with several times more ionizing radiation. The noise reduction is primarily due to lower amounts of Poisson-distributed quantum noise at higher X-ray doses in the primary beam. The procedure normally involves patient recall with the associated cost and anxiety. The film digitization system also produces magnified images by means of electro-optics. Digital magnification of mammograms results in images that are superior to the conventional 36 Chapter 4. Clinical Evaluation of the Film Digitization System 37 mammogram. Digital magnification is obtained by spatial magnification and photometric map ping via an optimized camera Look Up Table (LUT). The optical arrangement generally pro duces a demagnification, mapping the field of view to the CCD chip size. For example, for the case of our prototype unit, a square 50 micron sampling interval on the X-ray film maps a 64 mm x 51.2 mm portion of mammogram to 8.7 mm x 7 mm CCD area, i.e. a reduction of over 7.3 times. When viewed on an 18” display monitor the 1280 x 1024 image will be 360 x 288 mm in dimension, i.e. a net magnification of over 5.6. This is at least three times the magni fication produced by the conventional magnification mammography. Although this can almost be achieved using hand held magnifying glass, no improvement in contrast or conspicuity is possible in the latter case. It is the combination of magnification with contrast enhancement that provides the unique capabilities of Digital Magnification Mammography (DMM). For this study we performed no post-processing or enhancement of the acquired images. The acquisition camera was used with a linear LUT and a suitable choice of minimum and maximum grey levels. This effectively provides for an implied adjustment of window and level for the input grey levels. This operation is necessary since the acquisition is in 12 bits but the display is in 8 bits. This operation, coupled to the inherent contrast properties of the display unit leads to a pronounced improvement in the image contrast and conspicuity of microcalcifications. To assess the clinical utility of DMM a study was designed and carried out in collaboration with three radiologists from British Columbia Cancer Agency (BCCA). This was a preliminary clinical investigation of the application of DMM for the evaluation of mammograms and primary diagnosis of early breast cancer by experienced practicing radiologists. 4.1 Working Hypothesis: We tested the hypothesis that practicing radiologists can extract essentially the same diagnostic information from the original mammogram in regards to microcalcifications by using DMM in place of obtaining an extra magnification view, thus sparing the patient the added exposure to ionizing radiation. Chapter 4. Clinical Evaluation of the Film Digitization System 4.2 38 Materials and Methods: Mammograms of patients who had been referred to British Columbia Cancer Agency BCCA for further radiographic follow-up were obtained. Each of the patients had a pair of mammograms (craniocaudal and mediolateral oblique view) and one magnification mammography performed on her. Subsequently a biopsy was carried out on each of these patients, and we obtained the relevant pathology report on the excised tissue. All cases were reviewed by an experienced radiologist from BCCA, to ensure their suitability for this study. Specifically, each suitable case should have a visible suspicious cluster of microcalcifications. The films represented hard cases that conventionally would require magnification mammography. We selected 35 cases (three films per case i.e. 105 images) involving difficult-to-diagnose microcalcifications (as judged by a radiologist) without associated masses or any other signs of abnormality. Due to these requirements, we examined and rejected many more cases which either had a visible lesion present on the film, or the original mammograms were unavailable. Whenever the original mammograms were performed at an outside laboratory (normally within BC), the films were requested through the BCCA department of Diagnostic Imaging. Each film was subsequently digitized twice, at 50 um and at 100 ,um sampling intervals. The magnification views each needed only to be digitized at 100 1 um. In this way, each case is comprised of 5 images only one of which was used in the present study. Three radiologists participated in the study. The study consisted of two parts. In the first part, each radiologist individually reviewed the original mammograms and the digitally mag nified images on a monitor. In the second part, each participant again separately reviewed the original mammograms as well as the extra magnification films as in a conventional reading. In each case a detailed questionnaire, given in Appendix B, was filled by each observer. The questionnaire involved 18 questions in 6 categories pertaining to the features of microcalcifica tions. The questionnaire quantifies five attributes of the microcalcifications, namely: number in a cluster, shape, density, margination, and spatial arrangement. The questionnaire also records the overall clinical assessment. Chapter 4. Clinical Evaluation of the Film Digitization System 39 The two parts of the study were kept separate by both a time interval and random shuffling of the cases so that the observers had no recollection of their previous analysis. The three observers did not discuss the cases among themselves and no other clinical data or patient history was supplied to them. 4.3 Results: Two types of results were obtained namely, a qualitative assessment of the images by the radiologists and a quantitative comparison of the data obtained from the questionnaire. Qualitative assessment of the images by the radiologists: All the three radiologists expressed that they were able to see the details of the microcalci fications better on the monitor than on the original mammograms. Two other experienced radiologists with special interest in mammography evaluated the digitized images using the mammograms that they brought along with themselves. One of these radiologists is associated with BCCA, while the other is affiliated with a community hospital. All the clinicians expressed that they can see more microcalcifications on the monitor and they find the soft-copy display easier to use for primary diagnosis than the films. Although the radiologists were free to refer to the original mammograms during the evaluation of the digitized images, they chose not to do so. They would first look at the films to determine the relative location of the microcal cifications (i.e. in which quadrant of the breast they are) and then use the soft-copy display exclusively for detailed description of the microcalcifications. Quantitative comparison of the data: Comparative data for the three observers are given in table 4.1 in the form of a classification confusion matrix. CMM refers to the conventional magnification mammography and DMM refers to the digital magnification mammography. The last question in the questionnaire quan tifies the overall diagnostic impression of the microcalcifications and can be viewed as a two Chapter 4. Clinical Evaluation of the Film Digitization System Obeerver #1 CMM CMM Fbbserver #3 CMM 0MM Conipicuity 67 77 74 49 67 77 74 49 7 2 6 4 22 88 69 3? 44 73 66 31 2 7 1 8 78 38 49 -3 27 43 -14 Snitivity 20 6 20 6 3 6 3 8 B M B 22 4 19 7 B M B M to 16 7 19 B M 0MM Accuracy Malignath Pathology M Obeerver #2 Specificity Bcnign Pathology 0MM 1 1. 40 % 89 T T T 1 T Table 4.1: Confusion matrix for Benign (B) vs. Malignant (M) classification by conventional (CMM) and digital (DMM) magnification mammography. class classification problem, with the two classes labeled as either benign or malignant. In table 4.1 we have also included calculations of “Sensitivity”, “Specificity”, and “Accu racy” as used in the literature and reviewed in the Appendix. We also introduce a new metric called “Conspicuity” defined as (4.1) C=2A—1 where C is the conspicuity, and A is the accuracy. The rationale for this metric is this: a random classification for a two class problem, under the assumption of equal probabilities, gives an accuracy of 50% . If the image of any abnormality is so inconspicuous as to render the classification essentially random then Cz0. A highly conspicuous abnormality, with C1, should ideally result in 100% accuracy in classification. It is well known that some cancers are pathologically evident but radiologically occult. Therefore if ground truth, so far as the images are concerned, is taken to be the 2/3 majority of radiologists responses we obtain the results of table 4.2. Table 4.3 gives the extent of inter-observer agreements for the conventional reading of mam mograms. Observer # 1 agreed in 77% of cases with observer # 2, and 57% with observer # 3. The second and third observers only agreed in 40% of the cases. It can be seen that a considerable disagreement exists. This may be a typical example, and a natural consequence of the subjective and non-quantitative method of interpretation currently in practice. ___________ Chapter 4. Clinical Evaluation of the Film Digitization System Berign features Observer #1 CMM 0MM Observer #2 0MM 0MM Observer #3 0MM DMM 1 j MaJiguant features 1 j Sensitivity Specificity % % 41 ] j Accuracy % Conspicuity B M B •___ 23 1 23 1 0 11 0 11 100 96 97 94 100 96 97 94 B M 6 5 3 8 46 96 80 60 73 83 80 60 M 23 1 20 4 B M B M 11 13 8 16 1 10 0 11 91 46 60 20 100 33 54 9 W 1 j Table 4.2: Confusion matrix for Benign (B) vs. Malignant (M) classification by conventional (CMM) and digital (DMM) magnification mammography; ground truth is radiologists’ majority opinion. Observer # 2 Benign Malignant Observer # 1 Observer # 3 Benign Malignant Benign Malignant 22 7 10 19 1 5 2 4 Observer # 3 Benign Malignant [ 10 2 13 10 Table 4.3: Inter-observer (dis)agreements for Benign vs. Malignant classification in conventional magnification mammography. Chapter 4. Clinical Evaluation of the Film Digitization System 42 CMM Malignant DMM Benign Observer # 1 Benign Malignant 20 3 3 9 Observer #2 Benign Malignant 23 6 1 5 Observer #3 Benign Malignant 3 9 5 18 Table 4.4: Comparison of classification by conventional (CMM) and digital (DMM) magnifica tion mammography for each observer. It is instructive to compare the results of conventional and digital magnification for each individual observer. These are given in table 4.4. The classification problem is carried out in two steps: detection and visualization of micro calcifications, and determination of probability of malignancy. Both of these factors affect the final outcome. The second step however is not a function of image quality. We should therefore compare the radiologists responses for visibility of image features. Typical results from the first observer are given in table 4.5. 4.4 Discussion Observer # 1 had some preliminary exposure to the prototype unit, while the second and third observers had never used soft-copy images for primary diagnosis of mammograms. It appears from the results that some familiarity with this form of image display improves the outcome. We also observe that in both the conventional and digital arms of the study, the first observer was consistently closer to the “truth”. Considering the conventional analysis we assume that all three radiologists are equally well experienced, and yet their performance in predicting the probability of malignancy is very different. In particular the third observer is extra cautious and calls for far too many cases of malignancy. The accuracy of this observer is in fact worse than a random classification of 50% C j. CD CD 1+ c2 C CD C CD o ‘Cl) CD Cc aq :1 Chapter 4. Clinical Evaluation of the Film Digitization System 44 leading to negative values of image conspicuity. It is clear that this observer is influenced by her past experience and the particular group of patients that she normally examines. Our data base is a particularly difficult one to evaluate, and the probability distribution of malignant cases is different than the normal expectations of this observer. Concentrating on the results of the first observer, we note from tables 4.1 and 4.2 that the hit rate for both the conventional and digital magnification is the same. This indicates that within the statistical validity of the experiment the working hypothesis is confirmed. This result should be interpreted with caution. The conventional magnification procedure results in images that inherently have higher signal to noise ratios. There is however some loss of sharpness in the magnified images. The digital magnification has optical magnification and electronic contrast enhancement, but amplifies both the signal and the recorded noise in the original mammogram. This study does not compare these images from a technical view i.e. by signal to noise ratio, contrast, etc. but rather from a clinical point of view. As an example referring to table 4.5 we note that out of the 35 cases we have 9 true positive cases (TP=9), 20 true negative cases (TN=20), 3 false positive cases (FP=3) and 3 false negative cases (FN=3). This gives a false positive ratio of FPR=3/23=13%, or specificity of 87%, a false negative ratio of FNR=3/12=25%, or sensitivity of 75%. The overall accuracy is 29/35=83%. We have defined a conspicuity scale (or C-scale) as (2 x accuracy -1). Since a random diagnosis is likely to lead to 50% FPR and 50% FNR it will give a C-value of 0, i.e. the conspicuity of the abnormality has not increased by use of DMM. Alternatively if DMM could completely obviate the use of conventional magnification mammography we would have no false positives or false negatives, an accuracy of 100% and a C-value of 1. Any C-value above zero quantifies the contribution of DMM. In this case C = 0.66, i.e. a 66% increase in conspicuity. The above analysis can also be carried out in respect of individual questions in the ques tionnaire to quantify the contribution of DMM in increasing the conspicuity of various features of abnormalities. For the feature in question 1, the number of micro calcifications were classified as few ( less 2 than 10 per cm ) ). 2 or many (equal to or greater than 10 per cm Chapter 4. Clinical Evaluation of the Film Digitization System 45 Hence out of the 35 cases correct estimates were made in 24 cases giving an accuracy of 69%. In 8 cases DMM has underestimated the number of microcalcifications while in 3 cases it has overestimated it. The C-value is 0.37, indicating a 37% increase in conspicuity. Question 5 relates to spatial arrangement of microcalcifications and its confusion matrix is given in table 4.5. From table 4.5 we have an accuracy of 77% and C-value of 0.54, i.e. a 54% increase in conspicuity. Question 4 relates to margination of microcalcifications and classifies them into smooth or irregular as given in table 4.5, with an accuracy of 57% and a modest increase in conspicuity of 14%. Question 3 measures the density of calcifications, and classifies them into uniformly dense and poorly defined or smudgy. The confusion matrix is again given in table 4.5, from which we see an accuracy of 0.63 and C-value of 0.26. Question 2 relates to the shape of microcalcifications and is a seven class classification problem with multiple labels being allowed. The confusion matrix for this feature has 49 entries, i.e. larger than the total number of cases. For this reason this method of analysis was considered inappropriate. In summary the most important of the above results is the number of cases that a radiolo gist can correctly identify as benign or malignant without the use of conventional magnification mammography. For the three radiologists the average ratio was 26 out of 35 cases or 74%. For observer # 1, who had prior exposure to the system, this ratio was 83%. We have therefore shown that within the parameters of the experiment, the number of magnification mammogra-. phies may be reduced significantly by the use of digital magnification using the proposed image acquisition device. Chapter 5 Image Restoration In this chapter we employ filtering techniques in an attempt to reduce the image degradations and improve the detectability of abnormalities in digitized mammograms. 5.1 Problem Description Many more breast cancers are now detected in earlier stages as a result of greater participation of women in screening programs. The majority of early carcinomas of the breast are indicated by the presence of one or more clusters of microcalcifications on a mammogram. Detection of subtle, small, and low contrast microcalcifications, is therefore gaining increased significance. Over 30% of all early carcinomas of the breast are detected solely on the basis of the pres ence of a cluster of subtle microcalcifications. Usually a biopsy is performed and the presence of microcalcifications together with the associated malignancy is confirmed pathologically. A radiograph of the biopsy specimen is often taken employing contact radiography where the parameters are optimized for the best possible image. The specimen radiograph has a higher quality image due to lower noise associated with a much higher X-ray dose, and reduced scat tering in the thin specimen. The specimen radiograph almost always shows a larger number of microcalcifications than was visible in the original mammogram. Clearly then the mammogra phy imaging process misses much needed information and introduces image degradations that may play a critical role in the detection of early breast cancer. A closely related problem is that mammography usually underestiniates the extent of a le sion. From several correlated pathologic-radiologic studies it has been shown that the smallest 46 Chapter 5. Image Restoration 47 microcalcifications, normally in the periphery of the lesions, are not visible on the mammo gram. This problem is particularly significant when multifocality is involved and “satellite” microcalcifications are present in the vicinity of the primary tumor. These satellite microcal cifications can be found as far away as 4 cm from the primary lesion and are believed to be highly prognostic [89]. Two factors contribute to this phenomenon namely the observation system noise and the system blur. Image processing techniques may be used to overcome some of these degradations. 5.2 Image Processing Image processing can be considered as a preprocessing stage in digital analysis of mammograms. This preprocessing is used to calibrate the image, remove or reduce the random and fixed pattern noise and counter the effects of non-linear illumination and camera response. It is also used to correct the degrading effects of the transfer functions of the X-ray source, the imaging path and the screen-film receptor as well as to correct for the blurring effects of the camera Modulation Transfer Function (MTF). There is much published work in the literature in the areas of image enhancement and restoration, for example [94], [95], and [96]. We have treated noise reduction, image calibration and some general image enhancement techniques in chapter 3. While image enhancement techniques have their own advantages they do not consider the process of image degradation in the derivation of algorithms. We have already noted that the most subtle signs of abnormality are not visible in the raw mammogram or in its digitized version because of the noise and blur in the image. Degradations caused by the blurring ef fects and noise introduced by the X-ray machine and the screen-film detector appear in the mammogram, and similar degradations of the camera and the digitizing equipment affect the digitized mammogram. We postulate that if the characteristics of the noise and the different blur functions are known or obtained then these effects may be removed or reduced by im age restoration techniques. The restored image may then be directly viewed on the monitor, Chapter 5. Image Restoration 48 subjected to further image enhancement or processed for automated diagnosis. In this chapter We will consider image restoration from the degradiiig effects of the system transfer function. A prerequisite of restoration is the knowledge of the image formation model and the transfer functions of the system components as well as the noise model. The overall Modulation Transfer Function is considered to be due to both the screen-film receptor and the film digitization system. The most common restoration methods use linear filtering techniques such as Wiener filtering. In the following chapter We will extend this work to include the effects of signal-dependent noise. 5.3 Image Restoration Image restoration is a mathematical process in which operations are performed on an observed image so as to estimate the original object that would be observed if no degradations were present in the image formation system used. Basically the procedure is to model the image degradation effects of the system and then find and perform appropriate operations to ‘undo’ these degrading effects. Thus in order to effectively design a digital image restoration procedure, it is necessary to first quantify or characterize the image degradation effects of the physical imaging system, the image digitizer, and the image display. Due to the statistical nature of the degradations, ideal restoration is not possible. However> some degree of improvement may be feasible. We seek such improvements in the X-ray imaging of the female breast. To restore an image the blur functions and the noise strengths of the system must first be measured. The two components of the overall system are the X-ray imaging system and the two-dimensional CCD digitizing system, as shown in Fig. 5.1. These are modeled as a cascade of linear and shift invariant systems. This siniple model serves as a starting point and has the advantage of mathematical tractability. The Modulation Transfer Function (MTF) of the camera was measured as discussed in chapter 3. The MTF of the screen-film combination is normally available from the manufacturer. The noise components from the camera and the noise recorded on the film are here assumed to be additive and their strengths are estimated. Chapter 5. Image Restoration Matrix Film Object czzzzz> 49 X-ray Mammography Digitizing Camera Figure 5.1: The digital image formation model rj- 50 Chapter 5. Image Restoration The System Transfer Function 5.4 The relative importance of the camera MTF and the screen-film MTF needs careful attention. The values of the MTF of the best available X-ray screen-film combination at frequencies above 16 cycles/mm are known to be less than 5% of the MTF’s maximum value [97]. This implies that image details associated with these frequency components are severely attenuated resulting in blurring of the image. The mammographic image can therefore be considered to be band limited. The camera MTF’s value is over 50% at this limit (16 cycles/mm) indicating that the screen-film (and not the digitizing camera) limits the detection of microcalcifications. This relationship depends of course on the geometrical magnification provided by the lens, and holds only when the object and its image are of equal size. If the mammogram is imaged at lower magnifications the effect of the camera MTF increases. A combined impulse response can be found by convolving the impulse response of the indi vidual elements in the imaging path. In the Fourier domain this translates into multiplication: H(f,f) = (5.1) j(f,f) H(f,f,) 8 H f is the screen-film MTF and 3 where H is the overall system modulation transfer function, H H is the digitizing camera MTF. The system is modeled as (5.2) g(x,y)= h(x,y)* f(x,y)-i- n(x,y) where * is the convolution operator, h(x,y) is the overall system blur function, and n(x,y) is the additive noise, and g(x,y) is the observed image. 5.5 The Wiener Filter The reconstruction ifiter is designed to minimize a certain estimation error e which may be defined in a variety of forms. In the Wiener filter formulation = E is defined as: dy} [f(x,y)f(x,y)j d Ex{JJ 2 (5.3) ________ Chapter 5. Image Restoration 51 where Exp is the expectation operator, f(x, y) is the original, undegraded, and unknown image and f(x,y) is the estimate of f(a,y). The expectation in equation (5.3) is performed over the entire spatial support of the image. The restored image F(f, f!,) in the transform domain is P(f, f!d) = H(f, f) where H(f, fr,) is the restoration filter, and G(f, G(f, fr,) is f) (5.4) the Fourier Transform of the observed image. For the assumed linear shift invariant system the Wiener filter in the Fourier domain is known to be x,jy f— ’ 1 H rJx,J!J) I Jy) — 2 W(f,f,) where H is the Fourier Transform of h, and WN and WF are the power spectrum of the noise and signal respectively. Note that since the acquisition system blur function, h(a, y), is commonly symmetric, the transfer function, H(f, fr,) is a real function, and thus the above restoration filter is also real. To calculate equation (5.5) we are faced with the fact that the power spectrum of the latent image WF, in practice, is unknown. In our first attempt to calculate (5.5) we estimate WF from the degraded images as follows. The Wiener filter (5.5) can be re-written as: H*.WF — H — where the (fr, fr,,) TT 11 TIT 2 I’Vp + TXT VVN (5.6) are implied. From (5.2) the observed image in the frequency domain is given by: G=H-F+N (5.7) Therefore the power spectrum of the observed image WG can be calculated as: 0 W = H 2 WF + WN (5.8) Chapter 5. Image Restoration 52 and hence substituting for the unknown WF we get H—’ —ff T WG-WN W (5.9) From equation (5.9) we note that the Wiener filter can be decomposed into a cascade of a smoothing filter, WGWw and the inverse filter, . The computations were carried using the Discrete Fourier Transform (DFT), and WG was replaced by I G 2 In this implementation, in order to minimize the effects of the zeros of H, we selected only the first L values of the blur function as being significant and set the rest to a fixed, small but non-zero constant H(L). L can be treated as a design parameter and its value may be chosen empirically. The noise power spectrum was estimated from selected smooth regions of the noisy and blurred images. The spectral energy of the relatively low contrast soft tissue images on a mammogram are generally concentrated near the zero frequency. The Wiener smoothing filter WGW_WN was also implemented separately and proved to be quite effective in smoothing mammographic images as judged visually. The alternative method of calculating equation (5.5) assumes the noise to signal power ratio a = to be equal to a constant. The value of this constant is adjusted empirically to obtain the best restored image, as judged visually. The filter takes the form of a sharpening high pass filter and assists in better visualization of details of microcalcifications. The filter transfer function is given in Figure 5.2. Although the first method of deriving WF from the observed image has been suggested in the literature, we found that the second method of assuming a constant noise to signal power ratio produced better results. 5.5.1 Iterative Restoration The Wiener filter of equation (5.6) is the optimum linear restoration filter in the Minimum Mean Square Error sense. As stated earlier the problem is that normally the power spectrum of the latent image Wp’ is not known a priori. In section 5.5 above we tried to estimate WF from the power spectrum of the observed image WG. This resulted in equation (5.9). Alternatively ______________________ Chapter 5. Image Restoration 53 16 14 (a) System MTF 1 2 .‘ (b) Hr with NSR=.O1 ,‘ .‘ (c) Hr with NSR=.OO1 1 0 Li_ . (c) (d) MTF * Hr(.Q1) / (e) MTF * Hr(.0O1) 8 / / / 6 ,. 4. Frequency (Cycles/mm) Figure 5.2: The Modulation Transfer Function of (a) the system; and the restoration filter with (b) noise to signal power ratio o = 0.01; and (c) c = 0.001; (d) and (e) show the improvement in MTF due to the application of the restoration filters Hr corresponding to (b) and (c). Chapter 5. Image Restoration 54 we assumed that WF has the same shape as WN and took to be a constant. In this section we will estimate WF iteratively. The first estimate of Wp’ is obtained from the observed image, by :1:1 c 2 (5.10) From this we formulate a restoration filter = (5.11) F +WN IHl W 2 This filter is used to obtain an estimate of the restored image = H$) G (5.12) Now a new estimate of WF can be obtained =1 E() 12 (5.13) The iteration may be stopped based on some convergence criteria. It was found effective enough to arbitrarily stop the procedure after a few iterations. 5.6 The Constrained Least Squares Filter The classical Constrained Least Square deconvolution procedure seeks to minimize the objective function: II QE 112 +) II where Q (G — HF) 112 (5.14) is a two dimensional high pass filter, and A is the Lagrange multiplier. The first term in the objective function imposes a smoothness criterion by minimizing the high frequency components of the estimate of the restored image F. The second term aims at a least squares fit to the observed image G. The Lagrange multiplier which is related to the extent of the noise in the observation is treated as a design parameter and adjusted empirically for ‘best’ Chapter 5. Image Restoration 55 results. It is well known that minimization of the above objective function yields the following Constrained Least Squares deconvolution filter [96]: H* Hr H 2 . IQ (5.15) 12 where y is a constant inversely proportional to the Lagrange multiplier. Following the conven tional approach the Laplacian operator was chosen, approximated by 0 q(x,y)= 1 0 10 —4 1 10 to enforce the smoothness constraint. Therefore Q(w,w) = —1 + 0.5cos(w) + 0.5cos(w) (5.16) Results 5.7 In order to compare the performances of the various restoration filters quantitatively it is necessary to have access to the ‘true’ latent image. For this purpose first a phantom resolution chart was imaged. The phantom image contains both high and low frequency objects of varying contrasts. This image was considered to be our ‘latent’ image f(x,y). The image was blurred with the system transfer function and independent white Gaussian noise of standard deviation a,-. = 5 or 10 grey levels was added to it, to form the ‘observed’ blurred and noisy (r, y) respectively. Each image was then filtered with the inverse filter, two 1 image g(x, y) or g implementations of the Wiener filter, the iterative filter, and the Constrained Least Squares filter. In each case we used the error metric e defined as ])2 e = (5.17) where images are of size M x M. While it is true that this metric does not correctly capture the visual quality of the images, it is mathematically simple to implement, and is in common Chapter 5. Image Restoration 56 Filter Parameter None Inverse Wiener Wiener2 Iterative Wiener Constrained Least Squares j Restoration Error e SNR=l SNR=80 SNR=100 SNRz1000 o, = 5 u = 5 ‘y = 0.001 ‘y 0.01 = 0.1 1 y= 36.5 98.1 42.4 27.1 26.7 53 73.3 44.2 45.2 24.7 31.8 37.6 Table 5.1: Comparison of restoration filter performances use. Table 5.1 gives a summary of the results. In this table Wiener2 refers to the cascade of a Wiener smoothing filter WGW_WN and an inverse filter as described in section 5.5 above. A sample of images is given in Figure 5.3. In this figure Wiener3 refers to the iterative implementation of the Wiener ifiter. It can be seen that the restoration filters can be optimized to sharpen the image detail, with little degradation due to the amplified noise. Since a large portion of the phantom consists of fiat areas with nearly constant grey levels, the amplification of the noise causes the error e to be larger in some cases than the error in the unrestored image. Judging visually, five of the output images show considerable improvement in visualization of image detail, with the Constrained Least Squares filter with = 0.01 producing the best result. This observation is consistent with the numerical results of table 5.1. Figure 5.4 shows the effect of the Wiener filters on the actual image of a test phantom obtained by the two-dimensional CCD digitizer. The object is placed at a distance far from the camera such that the optical image formed at the CCD plane contains substantiai energy at frequencies near the limit of the resolution capability of the camera. A blurred image is obtained in this way while the camera is maintained at its best focus setting. No additional degradation was introduced this time, and the impact of processing the image with the Wiener Chapter 5. Image Restoration = 00 q inverse of f’ oriqin.il. 57 91 of I I FI( I 9 UI I 0(1 q of uienii 1)1 bI’ii H fijij iiod Iv I sli 10055 lOll. 8 hI s_d. I • nil 0 5 1101 f foist sqi of t (If 001 d ..[‘ 0 I ii 10 0 11111 101 UI qI 10 1 Uk 11(1 5 1 001 t1 1 .q 1 5111’ 1 01 of U lOner 1 1 (j I Ii wiener sor 80 10 I ulIrul az)”,’ lJlt((1. ((f (1 Uk Ill I e fl 1 — 5 S. - hill I ref 1 S I • I -a- 1111111 III H qIof 11mm_I if t fI II I fIPIPII inli IfIJIJIJIJI 0 Iii 1 IJIJIJIJI ( I 1 00001 q p 1101 I i e i-immo . I —, j i 1 fiMflI IlsiSt ‘I 9 sql (PIll I I PIll I i — t inmi of 91 9 0 III i of jl (1111 1110f9 II I I 41 I Figure 5.3: Restoration of test images: f is the latent image, fj, is the blurred image, g and g 1 are the observed images; all other images are the results of restoration by various filters. Chapter 5. Image Restoration 58 filter was tested. A fixed noise to signal ratio (NSR) is assumed in the design of the Wiener filters. The value of this constant was adjusted empirically to obtain the best restored image, as judged visually. In Figure 5.4 the outputs labeled Wienerl and Wiener2 are produced by application of Wiener filters with constant noise to signal ratios of NSR=0.0l and NS]Iii.00l respectively. To compare the performance of the deconvolution filter with enhancement filters we include the effects of two high pass filters shown in Figure 5.4. Both of these filters are 3x3 convolutional filters with the following mask values: o —1 0 —1 5 —i o —1 0 ;sharp2 = —1 —1 —1 —l 9 —1 —1 —1 —1 These edge sharpening filters were implemented in the spatial domain and therefore were computationally faster and more economical in the storage requirements. It can be seen that the Wiener filter with NSR=.00l produces the sharpest image in which the four separate lines at the top half of the image are clearly distinguished. The same filters were used on a data base of 30 mammograms selected from diagnostic films of patients referred to the British Columbia Cancer Agency’s Vancouver clinic. The mammo grams represent typical cases containing clusters of low-contrast microcalcifications present in the normal breast parenchyma. The films were digitized with contiguous square pixels of 10Om x lOOjim. Figure 5.5 shows the effect of the Wiener filter with NSR=.00l on a mammogram. In Figure 5.5a the observed digitized image is shown. Figure 5.5b gives the result of the application of a different Wiener filter which restores the image from the effects of the digitizing system only, i.e. only the MTF of the camera is used in the derivation of the filter. Figure 5.5c is the result of a Wiener filter which restores the image from the combined effects of the digitizing system and the screen-film combination. The sharpening effect of the filter is clearly visible as are the characteristic textured patterns in the background due to the amplification of the system noise. The microcalcifications associated with the primary tumor can be resolved Chapter 5. Image Restoration 59 original Wienerl Wiener2 sharp2 sharpi Figure 5.4: A portion of a standard test phantom imaged at a far distance near the resolution limit of the camera. Results of restoration by the two Wiener filters: Wienerl with NSR=.O1 & Wiener2 with NSR=.OO1, and enhancement by 3x3 convolutional masks sharpi & sharp2 as defined in the text. Chapter 5. Image Restoration (a) 60 (b) (c) Figure 5.5: Restoration of a mammogram: (a) The observed image; (b) Restored image using Wiener filter (with NSR=.001) compensating for the camera MTF; (c) Restored image using Wiener filter (with NSR=.001) compensating for the combined system MTF. Chapter 5. Image Restoration 61 better, and several ‘satellite’ microcalcifications are now visible. It can be seen from Figure 5.5 that the effects of camera MTF at the sampling interval of 100 1 um can not be neglected. In summary in this chapter we have reported on our implementation of a number of restora tion algorithms. The performances of these filters with various parameters were compared. The effects of both the screen-film combination and the digitizing equipment were considered. The restored images are sharper and assist in better visualization of the image detail. Chapter 6 Restoration in the Presence of Signal-Dependent Noise The work described in Chapter 5 assumes a globally stationary image with signal-independent additive noise. In this chapter the restoration of mammographic images is extended to include the signal-dependent nature of radiographic noise. We consider a non-stationary image model and signal-dependent noise of photonic and film-grain origins. Both the camera blur arid the MTF of the screen-film combination are considered. The camera noise is minimized through averaging and background subtraction as described in chapter 3. The signal-dependent nature of the radiographic noise is modeled by a linear shift-invariant system and the relative strengths of various noise sources are compared. We investigate the application of two locally adaptive image smoothing filters to improve the signal to noise ratio of digitized mammogram images. To minimize the effects of the system blur a deconvolution filter is then applied in conjunction with these smoothing filters resulting in better visualization of image details. The deconvolution filter is based on the Minimum Mean Squared Error (MMSE) criteria, while the smoothing filters utilize the Bayesian and the Wiener criteria. Of the two smoothing filters the Bayesian estimator is found to outperform the adaptive Wiener filter. The filters are implemented in a real time processing environment as part of my MAMPRO mammographic image acquisition and analysis system. The objective of this approach is to facilitate the detection of microcalcifications at the earliest stage of their formation. In section 6.1, We characterize the image formation system and derive an image observatiomi model in section 6.2. In section 6.3 different image restoration procedures are discussed and designed so as to minimize or reduce the effects of the different 62 Chapter 6. Restoration in the Presence of Signal-Dependent Noise 63 blurs and noise degradations. 6.1 Image Formation System The system degradations are considered to originate from the cascade of two imaging systems as described in chapter 5. In chapter 5 we assumed a globally stationary image with additive Gaussian noise independent of the signal. In this chapter we study and develop a model to describe the formation of digitized mammographic images with special attention to the signaldependent nature of the noise. In order to find a more realistic image degradation model, the different components of the image formation system are first examined. Firstly there is the X ray image formation system which produces the image recorded on the film. The X-ray system is followed by the digitizing optoelectronic camera that produces the digital image matrix. We will first consider these two systems separately and then in section 6.2, we combine them to formulate an overall image observation model. 6.1.1 The X-ray Screen-Film Imaging System Radiographic images suffer from a number of degradations which are inherent in the image formation system. A block diagram of radiographic image formation is given in Figure 6.1. These degradations may be broadly divided into two categories of blur or unsharpness, and noise or statistical fluctuations in image intensity. The Image Blur: The image blur is due to the following four sources: a) the X-ray source; b) the geometry of imaging; c) the beam scattering by the subject; and d) the image detection and display components. The illumination received at the screen is non-uniform due to the X-ray source and the geometry of the imaging. The largest amount of exposure is received along the central beam, with the quantum fluence decreasing proportional to 3 cos O , where 0 is the angle of inclination from the normal [98]. This effect is modulated by the ‘heel’ effect of the anode. Fig. 6.1 The radiographic image formation model a) Chapter 6. Restoration in the Presence of Signal-Dependent Noise 65 Additionally, the focal spot on the anode of the X-ray tube has a finite size creating penumbral shadows in the image. The extent of this shadow depends on the object-film distance, and it affects the spatial resolution limit of the image. It is also this effect that limits the extent of useful magnification views to about two times only [99j. We can model this effect by a twodimensional convolution of the point spread function corresponding to the effective focal spot aperture with the image. The X-ray scatter in the breast reduces the contrast of the image. This is a major source of degradation. Breast compression devices and vibrating anti-scatter grids are used to reduce this effect. Nevertheless, X-ray scatter remains a major limitation. The overall effect of these limitations can be partially evaluated practically by comparing a pre-operative mammogram of a patient who has undergone biopsy, with the specimen radiograph. In specimen radiography the object-film distance is reduced and much of the breast mass responsible for beam scatter is absent. Additionally higher doses are employed leading to better contrast and reduced input quantum noise. The specimen is also centrally located reducing the effects of geometrical distortions. The screen-film combination as a detector has been studied extensively [97]. The intensifying fluorescent screen absorbs the X-ray photons and radiates many more photons in the visible range of wavelengths which expose the film. The X-ray absorption and re-radiation efficiency of the screen , as well as the photon absorption efficiency of the film , 2 i contribute to the overall contrast for a given dose to the patient. The amplification and scattering mechanisms of the screen are stochastic in nature. Rabbani [100] has shown that the uncorrelated component of the quantum noise passes through the screen unaltered while the correlated component is filtered by the system contrast transfer function. The amplification m is described by its probability distribution function Pr(m), having mean ñi, and variance u. The scattering process can be modeled by a two-dimensional linear convolution operation and forms the major component of the screen film Modulation Transfer Function (MTF). For contrasts below about 6%, the system noise is the limiting factor in the visibility of the Chapter 6. Restoration in the Presence of Signal-Dependent Noise 66 image details while the system MTF has little effect on it [101]. For higher contrasts however, the MTF contributes significantly to these limitations. This fact indicates that in digitized mammograms restoring the image from the effects of the system MTF may result in better visualization of smaller objects. Finally the response of the film to the incident photons is non-linear. This non-linearity is described by the D D = — log E characteristic curve. -y log E + /3 we note that both Although this is commonly written as and /3 are functions of the exposure level E. The contrast transfer function of the screen-film combination is both a function of exposure (due to the characteristic curve) and frequency (due to MTF). Image Noise The process of X-ray image formation is also associated with noise which is generated by four sources. These noise sources are: a) the quantum noise due to the discrete nature of the X-ray photons; b) the screen mottle due to the stochastic nature of amplifications and scattering; c) the screen structure noise due to its inhomogeneous phosphor coating; and d) the film grain noise of the emulsion coating. The screen structure noise is generally considered to contribute less than 2% to the overall noise [97]. In Figure 6.1, X-ray quantum fluence Q, in the (i, j)th pixel contains spatial non-uniformities due to the geometry of imaging and the inherent photon noise which has a Poisson distribution. The effect of focal spot size is modeled by convolution with its aperture function, and that of the X-ray scatter is represented as an effective low pass filter. The effects due to cos 3 0 term are quite small since for typical geometries involved 0 is less than 5°. We use the Kodak Min-R screen together with a Kodak Ortho-M film. The densitometric data, MTF, and other parameters of this screen-film combination have been reported by Bunch [97]. Details of this model are further described in section 6.2. Chapter 6. Restoration in the Presence of Signal-Dependent Noise 6.1.2 67 The CCD Camera In addition to the X-ray, the digitizing camera also blurs the image. The finite size of each pixel gives the aperture function and provides the theoretical limit to the MTF of the camera. The lens contribution to the blur can at best be limited to the diffraction properties of the optical components employed. As described in chapter 3, we use a two-dimensional CCD (Kodak KAF 1400) with 100% fill factor and 1035 x 1320 square pixels of 6.8 um per side. This gives an MTF of the form ‘sinc(w)sinc(w)’. This function is multiplied by the MTF of the lens. We used a Nikon Nikkor 55 mm lens which, for the best focus conditions has a minimum MTF of 0.6 at the Nyquest sampling frequency of our sensor. The MTFs of the camera and that of the screen-film are shown in the same plot in Figure 6.2. The spatial frequency axis refers to the frequency content of images formed on the film and on the CCD planes, respectively. When comparing these two MTFs we note that these two curves should refer to the same imaging plane. If the geometry of imaging and the focal length of the lens are chosen such that a ‘life size’ image is formed at the CCD then the two MTF curves are directly comparable. Under these conditions we note that the screen-film combination (and not the digitizing camera) is the limiting factor in spatial resolution. The field of view is now restricted to the area of the CCD, i.e. a mere 7 mm x 9 mm. As we increase the field of view the effective camera MTF referred to the object (mammogram) plane deteriorates. In addition to the noise present on the developed film the digitizing system adds another noise element. This noise element is due to the following four sources: a) the illumination source; b) the lens; c) the CCD and its associated amplifier; and d) the analogue to digital converter. Figure 6.3 gives a block diagram of the CCD image formation from the mammogram. We combine these four sources and divide the overall noise present in the data into two types, the fixed pattern noise such as the optical shading of the light source and the aberrations of the lens, and the random noise of optical and electronic origin. The fixed pattern noise can generally be corrected by image calibration, while the random component of the camera noise may be reduced using averaging. The final calibrated image therefore will contain a small Chapter 6. Restoration in the Presence of Signal-Dependent Noise 68 Modulation Transfer Functions 1 .0 0.9 0.8 0.7 - 0.6 H 0.5 0.4 0.3 0.2 0.1 0.0 0 8 16 24 32 40 48 56 Frequency (Cycles/mm) Figure 6.2: The Modulation Transfer Functions 64 72 80 Chapter 6. Restoration in the Presence of Signal-Dependent Noise 69 amount of observation noise, which is due to the camera and uncorrelated to the image, and a more significant radiographic noise which is signal-dependent. 6.2 Image Observation Model It is customary in image restoration literature to consider the image observation model to be linear. In a commonly used model the observed image, g, is considered to be the result of linear convolution of the latent image, f, with a blur function, h, and addition of independent, zero- mean, white, Gaussian noise. While this model has the advantage of mathematical tractability, it is an over simplification of the actual situation. We have employed this model and have shown in chapter 5 that improvement in the appear ance of digitized mammographic images is possible. In particular the blur and noise contributed by the digitizing camera can be modeled in this form. The X-ray image formation however, can not be modeled in this way. Specifically, radiographic noise is strongly signal-dependent and the film density vs. log exposure characteristic curve introduces non-linearity in the convolution term. We formulated an image observation model for the X-ray system as shown in Figure 6.4 based on the physical description of the system given in Figure 6.1. In this model the weakening of the off-axis rays may be compensated for by a pointwise normalization of each pixel value by the cos 3 0 term. The input is the number of X-ray quanta received at the screen per pixel, and the output is g , the observed optical density of each pixel. We have included three 3 sources of blur in this model: hfocalspot, hscatter, and hscreen. There are also three scaling factors: i is the mean screen amplification ratio; and and 2 are the absorption efficiencies of the screen and the film respectively. The three noise sources represent n 1 the correlated, and 2 the uncorrelated, components of the input quantum noise, and n n 3 the associated film-grain noise. This model may be simplified by making the following observations. The scalar factors may be taken out of the system block diagram and reflected at the input. The blur due to the focal Fig. 6.3.The digitizing camera image formation model I I ft Fig. 6.4. A model of radiographic image formation 2 n 3 n I. Chapter 6. Restoration in the Presence of Signal-Dependent Noise 72 spot size is a function of the relative distance of the lesion of interest and the screen. This is of course not known a priori and for the cases where the object of interest is in contact with the screen there are no blur contributions from the focal spot. Finally we assume that the X-ray scatter in the subject is largely absorbed by the vibrating grid. An effective grid system ensures that the scattered photons are not absorbed by the screen and therefore the scatter can be modeled by a gain factor independent of spatial frequency. Based on these considerations a simplified version of the above model is given in Figure 6.5. In this model we have ignored the effects of the weakening of the off-axis rays, the focal spot size and the X-ray scatter in the subject. In Figure 6.5 the ‘ideal’ image is the number of light quanta, , 23 absorbed by the film q in each pixel (i,j) f where Q, = q,, (6.1) = is the number of X-ray quanta received at the screen per pixel, amplification ratio, and = ?]12, f is the mean screen i.e. the combined absorption efficiencies of the screen-film combination. Note that in Figure 6.5 the images while ñi f, fi, and f2 are in the exposure domain and g are in the optical density domain. P is the non-linear, D — log E, characteristic function of the film. The noise sources in Figure 6.5 n , n 1 3 are zero-mean additive signal-dependent , and n 2 white noise sources uncorrelated with each other. The effects of the scalar factors of Figure 6.4 are now incorporated in the magnitudes of these sources. Specifically, n 2 is the uncorrelated component of the Poisson noise of the X-ray source and is generated as in Figure 6.6. = ./j.n’ (6.2) where n’ is a zero-mean unit-variance Gaussian random variable. fi is a random variable with Poisson probability distribution whose expected value is f. The noise component n 1 is the amplified noise due to the screen amplification fluctuations: = ..,/jn’ 1 k (6.3) Chapter 6. Restoration in the Presence of Signal-Dependent Noise 1 g )+n 1 g=f[bf*(f+n J 2 3 +n * = 13 q * + fl, 11 iii Q Figure 6.5: A simplified model of radiographic image formation 73 Chapter 6. Restoration in the Presence of Signal-Dependent Noise 74 It is generated according to the block diagram of Figure 6.7. The gain factor is = where E iJth (i + -i-) (6.4) is the excess Poisson noise. The film-grain noise is represented by n : 3 3 n = k f 3 n’ (6.5) and k = 3 a.logioe (6.6) where a is the average film-grain area, and / is the sampling interval, i.e. the pixel size. For optimally exposed film 3 is taken to be 0.5. From Figure 6.5 we can write the following image observation equations g = 31 P[h *(f + f [h 31 * n + ] ni)+ 2 f * ni 3 f + (h n ) ]+n 3 +2 where h f is the point spread function of the screen-film and 3 the constant parameters for our screen fihn combination (6.7) 713 ( * (6.8) signifies linear convolution. Using 0.58, ñz = 284, = 112, and MTF and densitometric data as published in [97]) we note that n 2 is at least two orders of magnitude smaller than either n 1 or n . The latter two quantities are of comparable magnitudes. We will 3 therefore ignore n 2 from now on and write: 3 * g=P[h ) 1 j+n (f+n (6.9) The problem of film non-linearity may be handled in any one of the following three ways: i) the function I’ may be explicitly incorporated in the restoration filter; ii) the processing may be done in the exposure domain where a linear relationship exists; or iii) a small-signal model may be used to derive a linear equation. Since mammographic images are of very low contrast we will use the small-signal analysis assumption. If a linear approximation to the above non-linear Chapter 6. Restoration in the Presence of Signal-Dependent Noise poisson generator ii’ = 2 n Gauss(O,1) Figure 6.6: Generation of the noise source n 2 75 Chapter 6. Restoration in the Presence of Signal-Dependent Noise 1” ni 1c n, = k, = k, ‘JT i.[ n’ (1 + / ffl) c = excess Poisson noise ffl = mean screen amplification Figure 6.7: Generation of the noise source n 1 76 Chapter 6. Restoration in the Presence of Signal-Dependent Noise 77 equation is made then 9 4 n = hf * f 4 +n (6.10) * 1 +n n 3 (6.11) where all variables are now in the density domain and the appropriate conversion constants are incorporated in them. n 4 is now the total radiographic noise. This relation can be combined with the linear convolution model of digitization by the CCD camera. We therefore consider that first the ‘ideal’ image, f(i,j), is blurred by the system impulse response, h(i, j), and a small amount of uncorrelated camera observation noise, flc(,j) is added to it to produce the blurred image fb(i,j). Subsequently the signal-dependent radiographic noise n (i,j) is added to it, which results in the final observed image 9 4 (i,j): f&(i,j) = h(i,j) g(i,j) = (i,j) 4 f&(i,j) + n * f(i,j) + n(i,) (6.12) (6.13) The system impulse response h is due to the combined effect of the camera and the screen-film system, i.e.: h(i,j) = f(i,j) 3 h * h(i,j) (6.14) Note that the additive noise model above is not restrictive since any multiplicative noise can be reformulated as additive signal-dependent noise. We consider n, to be a white noise field with zero-mean Gaussian distribution. 6.3 Image Restoration We will use the image formation model of equations (6.12) and (6.13). According to this model the observed image was formed in two steps. The latent image f, was first blurred (with the addition of n) to form fb, and then contaminated by signal-dependent noise n . We 4 therefore divide the restoration problem into two steps. In the first step we apply smoothing Chapter 6. Restoration in the Presence of Signal-Dependent Noise 78 techniques using local statistics to ‘clean-up’ the image from the signal-dependent noise n . In 4 the second step we apply a classical Wiener filter to deblur the image from the combined effects of the system’s Modulation Transfer Function and the noise n. The following two criteria of optimality are used in the first step in deriving two smoothing filters: the minimum mean squared error (MMSE), and the maximum a posteriori (MAP) joint probability. The required image smoothing may be achieved using either the global image or the local image approach. In the first approach the image is assumed to be a stationary random process. In chapter 5 a Wiener smoothing filter is employed, assuming the system noise is a white Gaussian process independent of the image grey levels. This is a non-adaptive global approach in which the MMSE criterion of optimality is applied over the whole image. In this chapter we study the more realistic of the two approaches (i.e. the one which uses local processing). We postulate that the breast image is non-stationary due to the presence of structure in the parenchymal pattern. This is particularly true of mammograms containing microcalcifications or masses. Additionally, the blur process is a local operation and therefore it is reasonable to expect that the restoration should also be performed locally. The grey level histogram of a complete mammogram is commonly not Gaussian. If we subtract the local mean from each pixel however, the resulting image grey levels have a nearly normal distribution. Therefore in this work we consider the image model to be Gaussian with non-stationary mean and non-stationary variance (NMNV). 6.3.1 Local adaptive Wiener smoothing filter For the NMNV model [102] the image is considered to be Gaussian only locally in small neigh borhoods. Using the MMSE criterion locally will lead to an adaptive local linear minimum mean square error (LLMMSE) filter. The estimated pixel value fb [103] is: fg+ where and c.r U4(g_g) 2 U; (6.15) are the local mean and variance of the observed image, and u 4 is the noise variance. The required local statistics are estimated from the observed image and a priori Chapter 6. Restoration in the Presence of Signal-Dependent Noise 79 knowledge about the image is not required. The noise power is first estimated for each pixel and this knowledge is used to calculate a new estimate of the signal according to equation (6.15). If the noise power is negligibly small (i.e. u >> u ) then the estimated pixel value is very close to its observed value. At the other 4 extreme if all of the observed power is due to noise (i.e. u ) then the best estimate of 4 u the signal is the local average. The noise power is estimated from a knowledge of the noise model. The total radiographic noise n 4 is the sum of the film-grain noise 3 n and the noise due to the quantum mottle. The , noise power U4 can therefore be readily estimated for each pixel. For the quantum mottle n , we observe that the underlying noise process is due to the 1 discrete nature of the X-ray photons and therefore has a Poisson probability distribution. The grey level fb(i,j) for the pixel (i,j) is related to the number of photons incident on it. The grey level has a code value between zero and 255 which is a linear quantization of the film density D(i,j), and the film density is a known function of the exposure. Therefore the number of incident photons can be calculated for each pixel. Since the mean and variance of a Poisson distribution are equal, the exposure va’ue will directly determine the noise power a . 1 Finally, since the two components of the radiographic noise are independent of each other the combined noise power can be obtained as 1 u 6.3.2 MTF 2 (6.16) Bayesian smoothing filter The maximum a posteviori (MAP) filter for the above case of the non-stationary mean and non-stationary variance (NMNV) image model and Poisson noise has the form [104]: - fb = 2 where the average power of the blurred image u 6 is obtained from the observed image 1 (6.7) Chapter 6. Restoration in the Presence of Signal-Dependent Noise = 6.3.3 max[(u — 80 4), 0]. (6.18) The Deconvolution Filter After the application of one of the above smoothing filters, the resultant image of the blurred image obtain f, f& f, is an estimate (equation 6.13). A modified Wiener restoration filter was selected to an estimate of the ideal image: P(w,.,w) H,.(w,LJ) = (6.19) ) 71 H,.(w,w).Fb(w,w I II*(w w) H(w,w,) (6.20) where (w,w), and H(w,w) are Fourier Transforms off, fb, and h respectively, 6 F * is the complex conjugation. H,. is the reconstruction filter in the frequency domain and and is a measure of the noise to signal power ratio. 6.4 Results We implemented both the LLMMSE filter (equation 6.15) and the MAP filter (equation 6.17), followed by the deconvolution filter (equation 6.20). We note that the application of these filters may or may not be required depending on our operating conditions interval ) ( such as the sampling and also depending on the application in mind. For example in any reading of a mammogram where attention to fine spatial and photometric details is not required, image deconvolution will not be necessary. An example of this is when mammograms of the left and right breast of the same woman are being examined for the detection of bilateral asymmetry. In such cases large pixel sizes ( e.g. 200 m per side ) may be used which leads to smaller, and therefore more manageable images. The radiographic noise will also be sniafler in these images obviating the need for image smoothing. We evaluated each of the above three filters individually and also as combinations. The two smoothing filters may be used individually in cases where radiographic noise is judged Chapter 6. Restoration in the Presence of Signal-Dependent Noise 81 to be the limiting factor in interpreting the films. The deconvolution filter may be employed alone in cases where the image blur is the principal consideration and oniy a small amount of uncorrelated noise is present. This would be the case when a relatively large sampling interval (i.e. pixel size ) is utilized. The cascade of adaptive smoothing and deblurring filters should be applied in cases where both radiographic noise and image blur are present. We have chosen visual assessment of processed images to determine image quality. Various quantitative measures, such as MMSE have also been proposed in the literature. Calculation of MMSE however requires knowledge of the ideal image and may be performed using simulated degradations. To evaluate the effectiveness of each one of the two smoothing filters we considered a portion of a digitized mammogram containing a cluster of microcalcifications suspicious of malignancy. Radiographic noise was then added to this image. Although this will exaggerate the total amount of noise present in a mammogram it enables us to assess, more readily, the perfor mance of the smoothing filters. It also represents a inamniogram obtained under less than the ideal imaging conditions. A 5x5 square window was used to calculate local statistics. The performance factor P for each filter is defined as the square root of the ratio of average noise power before smoothing , to the average noise power after smoothing = 4 = [104]: (6.21) - (fb(i,j) — fb(i,))) (6.22) (6.23) Table 6.1 shows the advantages of our noise compensation procedure and gives a summary of the performance factor for each filter. The extent of the noise was controlled by a constant multiplier factor, ), in these experiments. The cases of A = 1, 0.5, & 0 correspond to ‘severe’, ‘moderate’, and no noise respectively. Finally, to examine the effect of the combined operation we processed 30 images from our Chapter 6. Restoration in the Presence of Signal-Dependent Noise 82 data base. Following smoothing, each smoothed image was deblurred by the Wiener filter. In this implementation, we again considered the signal to noise ratio to be a constant independent of the spatial frequency of the image. We used a signal to noise power ratio of 20 in the deconvolution filter. The images were displayed on a high resolution (1024 x 1280 pixels, colour) monitor, and a radiologist and an image processing engineer reviewed the images. In all cases the processed images were sharper and revealed greater amounts of image detail. Generally the noise content of the images was also increased. This, however did not interfere with identification of microcal cification clusters. The opinion of the reviewers was that application of these processing steps normally assisted in better visualization of image detail. Figure 6.8 presents a typical image at various processing stages. Here the mammogram contains ‘real’ ( i.e. not simulated ) radiographic noise. The results of the deconvolution both before and after smoothing of the mammograms are shown. Figure 6.8a is the observed noisy and blurred mammogram, and Figure 6.8d is the restoration of it using the Wiener filter but without any prior smoothing. The amplification of the noise obscures much of the detail of this image. Figure 6.8b is the result of processing the image of Figure 6.8a with the LLMMSE smoothing filter, and Figure 6.8e is its restoration using the Wiener filter. Clearly the image detail has been sharpened without any significant gain in the noise. Figure 6.8c is the result of smoothing of Figure 6.8a with the MAP filter, and Figure 6.8f is the deblurring of Figure 6.8c with the Wiener filter. The beneficial effects of these noise smoothing and detail sharpening filters are clearly visible. The visual appearance of these images are consistent with the measured values of the filter performance factors. Since the MAP filter shows higher performance factors, this filter was implemented as part of the real-tinie mammographic image acquisition and analysis system. Chapter 6. Restoration in the Presence of Signal-Dependent Noise 83 Table 6.1: Filter performance factors LLMMSE MAP 6.5 Moderate Noise 1.63 1.76 Severe Noise 1.47 1.63 Summary In summary radiographic images suffer from both signal-dependent noise and system blur. We have designed and implemented locally adaptive smoothing filters to reduce the effect of the noise. The smoothed images are then subjected to a deblurring algorithm to reduce the effects of system blur. The resulting images improve the visibility of subtle signs of abnormality and thus help the earlier detection of breast cancer. Chapter 6. Restoration in the Presence of Signal-Dependent Noise 84 Figure 6.8: Adaptive smoothing and deconvolution of a digitized mammogram (a) observed image, (b) smoothed by LLMMSE filter, (c) smoothed by MAP filter, (d,e,f), result of Wiener deconvolution of images (a,b ,c), respectively. Chapter 7 Reduction of Boundary Artifacts in Image Restoration In this chapter we investigate the problem of restoration of large images. When a mammogram is digitized to high resolutions, it becomes computationally very expensive to process the whole image at once. When the image is subdivided into smaller images artifacts are introduced due to boundary truncation effects. In this chapter we analyse these effects and suggest a simple novel approach to overcome these effects. The abrupt boundary truncation of an image introduces artifacts in the restored image that may be visually objectionable. These artifacts are particularly severe when the restoration filter contains significant high frequency components which is usually the case. The traditional solution is to smooth the image data using special window functions such as Hamming or trape zoidal windows, before applying the restoration filter. This method improves the results but still distorts the image, especially at the margins. Instead of the customary ‘linear’ convolution of the image with the restoration filter we examine a different procedure. This procedure is simple and exploits the natural property of ‘circular’ or periodic convolution of the Discrete Fourier Transform. Instead of padding the image by zeros, it is padded by a reflected version of it. this is followed by ‘circular’ convolution with the restoration filter. This procedure is shown to lead to much better restoration results than the windowing techniques. The computational effort is also improved since our method requires half the number of computations required by the conventional linear de-convolution method. 85 Chapter 7. Reduction of Boundary Artifacts in Image Restoration 7.1 86 Image Restoration Artifacts To remove certain well characterized degradations from images, image restoration techniques are employed. These however generally introduce various artifacts of their own in the image. Several authors have studied these artifacts. Tekaip and Sezan have identified and detailed four types of image restoration artifacts, namely the filtered noise, the filter deviation, the Point Spread Function error and the boundary truncation artifacts [105]. In particular the artifacts associated with the image boundary truncation can dominate the restored image under certain conditions. Woods was the first to discuss the boundary truncation artifact [106]. White and Brzakovic have considered the problem of image extension to minimize errors in convolution of images with spatial masks [1071. Tan, Lim and Tan have discussed the boundary artifacts in image restoration, but they have not considered extending the image to reduce these artifacts [108, 1091. We analyze the boundary artifact problem in section 7.2, and proceed to suggest a solution to it in section 7.3. In section 7.2.1 we discuss the problems arising when linear deconvolution is used, and in section 7.2.2 We will discuss the related problem of aliasing due to periodic convolutions. Section 7.4 considers the computational load of the new approach and we report the experimental results in section 7.5. The image degradation model is usually assumed to be that of a linear convolution and additive noise. In this chapter we do not consider the noise and concentrate on the deblurring. Even though the omission of the noise may not hold in actual situations, this analysis serves as an appropriate starting point due to its mathematical tractability. This image restoration model is schematically depicted in Figure 7.1. Consider an image of a scene obtained by a photographic camera. The observed image (x, y) is cut from the background due to the finite aperture of the image acquisition device. 1 g Since the reconstruction filter is usually a convolutional filter, and not a point-wise operation, the restored image f(x,y) will contain artifacts due to the abrupt truncated boundaries in (x, y). 1 g This problem is of special relevance in the restoration of large images, where a large image is Chapter 7. Reduction of Boundary Artifacts in Image Restoration f(x,y) g(x,y).. .g(x,y) Figure 7.1: The image degradation and restoration model in the absence of noise. 87 f(x,y) Chapter 7. Reduction of Boundary Artifacts in Image Restoration 88 subdivided into subimages and each subimage is restored separately. It is well known that the distribution of intensities in a large image is generally not Gaussian. Typical images of common scenes have non-stationary mean and variance. It has been suggested that the removal of the local mean from the image will render the image Gaussian [112]. It is for this reason that the observed image is divided into smaller size subimages and then each of these subimages is restored separately. Additionally the computational load of restoration is reduced by image subdivision. The boundaries of each subimage introduced by dividing the large image lead to undesirable artifacts. The traditional solution to the image truncation problem has been to employ special win dowing functions, such as Hamming or Hanning, to smooth the effect of truncation. This is followed by zero padding the observed image to the length of the restoration filter, and DFT operations are then used for deconvolution. The zero padding is necessary to achieve linear deconvolution (convolution) since such DFT processes are by nature circular or periodic ones. Here we examine the effects of some simple measures to minimize the effects of these bound ary truncation artifacts. We show that for this problem, the natural circular deconvolution (convolution) process of a DFT can be used to advantage in lessening the image truncation effects. We find that when the image is globally stationary, then padding with the image itself (instead of zeros) leads to better expected results. We then relax the condition of global sta tionarity to that of local stationarity and propose a simple padding method by which the results are further improved. Whenever convenient, and without loss of generality, we will illustrate our approach using one-dimensional images. 7.2 Problem description Consider a long one-dimensional signal f(x) whose length >> M, and its blurred version g(x). Let 1 g ( x) be the observed truncated section of g(x). gi(n), the discretized version of 1 g ( x), is limited to a width of M pixels. Assume the impulse response of the reconstruction filter hr(n) has L,. non-zero terms. After processing gi(n) by h(n) the resulting image f(m) is of length Chapter 7. Reduction of Boundary Artifacts in Image Restoration L + M — 1. The first and last L,. artifact and only the middle M — — 89 1 terms of f(n) are affected by the boundary truncation L,. + 1 terms are true estimates of f(n). If L,. << M and if h is a smoothing function i.e. a low pass filter, then the boundary artifact has the effect of smoothing a band of L,. — 1 pixels wide around the image. This effect is usually of little visual objection. It is common however for the degrading function to be a blurring one i.e. to be a low pass function and therefore hr() tends to be a high pass filter i.e. the transfer function corresponding to hr() will contain high frequency components. Moreover, the reconstruction filter is commonly designed and applied in the frequency domain such as in inverse filtering, Wiener filtering and other associated techniques. Under these circumstances and particularly due to sharp transitions in the transfer function of the reconstruction filter, the length of the impulse response hr(n), L, is likely to be large. The effects of a high pass filter on an image with sharp boundaries are the introduction of highly objectionable artifacts. If L,. M + 1 these artifacts spread over the whole image. Fig. 7.3 gives an example of this effect. In Fig. 7.3 the degraded image is obtained by imaging a step edge with a CCII camera using a square aperture. The camera sensor has 1317 x 1035 pixels with 100% fill factor such that there are no gaps between adjacent pixels. Each pixel is 6.8 zm x 6.8 im giving a Nyquist bandwidth of 73 cycles/mm. The aperture function sets the theoretical limit of the Modulation Transfer Function (MTF). A further reduction to the camera bandwidth is contributed by the lens. We have determined the MTF of the digitizing camera experimentally and this is shown in Fig. 7.2. This MTF has positive non-zero values at all frequencies below the Nyquist limit (64 pixels in our case). The reconstruction filter is a modified smooth Wiener filter which was designed from the knowledge of the camera Modulation Transfer Function. We used the following filter Hr I H-f-c (7.1) where c = 0.01 gave the visually best restored image (Fig. 7.3b). The ringing artifacts produced by frequency domain filtering are clearly visible in the restored image of Fig. 7.3b. Chapter 7. Reduction of Boundary Artifacts in Image Restoration 90 Modu’ation Transfer Function 1 .0 0.9 0.8 a 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0 8 16 24 32 40 48 56 64 72 Frequency (Cycles/mm) Figure 7.2: Modulation Transfer Function of the digitizing camera. 80 Chapter 7. Reduction of Boundary Artifacts in Image Restoration (a) (b) (c) 91 (d) Figure 7.3: Restoration of a step edge acquired by a digitizing photographic camera: a) The observed image of the step edge; b) The restored image using linear deconvolution, notice the boundary truncation artifact; c)The restored image of the step edge, a trapezoidal was applied to the data prior to its deconvolution; d)The restored image of the step edge using the proposed approach. Chapter 7. Reduction of Boundary Artifacts in Image Restoration 7.2.1 92 Analysis of errors in using the conventional linear de-convolution In this section we investigate the sources and nature of errors arising from the abrupt truncation of the image. To do that we will first consider the restoration of a blurred image in the absence of picture truncation. Assume there is no boundary truncation and the object f(n) is surrounded by a background of grey levels equal to zeros: {0, 0, 0, fo ..., fRi, 0, 0, 0} Assume f(n), the finite impulse response h(n) and thus the resulting g(n) are of finite lengths R,LandMrespectively(M = R+L—1). Assume that g(n)= {0,0,0,go,...,gR_1,gR,...,gM,0,0,0} is completely observed, then g(n) where * = h(n) = 0 * f(n) 0 < n < M — 1 (7.2) otherwise represents the linear convolution operation. To obtain f(n) we shall operate on g(n) by h(n) whose length is Lr. Since the length of g(n) is M, L is usually chosen to be In the rest of this paper we assume L,. M. )(n)= h(n)*g(n) J(n) is of length M + L — 1 2M — M. 0< ii < M+ Lr 1 (7.3) 1. If we perform our computation using a DFT then we pad h,f,g and hr by zeros (for example, as recommended by [110]). We take the DFT of size (M + Lr — 1)-point. Without loss of generality form now on we will take the N-point DFT, where N = 2M. We get G(k) = H(k) F(k) (7.4) P(k) = Hr(k) H(k) F(k) (7.5) Applying Hr(k) we get . . Chapter 7. Reduction of Boundary Artifacts in Image Restoration and choosing H(k) as the inverse filter (assuming H(k) P(k) 93 0 Vk) we obtain F(k) (7.6) i.e. F(k) is recovered exactly. In the above we studied the case when f(n) is of length Rand g(n) is of length M and we took N-point DFT (N = = R+L—i 2M). Suppose we now increase the length of f(n) but we still apply N-point DFTs. We have two cases: Case 1: Length of f(n) is > R but N L + 1: here the length of g(n) is — N and taking < N-point DFT will result in equation (7.4) as above. Case 2: Length of f(n) is truncate g(n) N by gN(n) = = {fo,...,fUN_l}. > {0, 0,0, go, {go, N ..., — L + 1. The resulting image g(n) is of length , gN, 1 gj.r ...} > N. Let us and denote the segment of g(n) whose length is Let us denote the corresponding N pixels of f(m) by fN(n) = for 0< n< N—i we get gN(n) h(n) = * fN(n) 0 n < N — 1 (7.7) Let the N-point DFT of gN(n) and fN(n) be GN(k) and FN(k) respectively. We get H(k) FN(k) (7.8) H(k) . FN(k) + EN(k) (7.9) GN(k) but rather GN(Ic) = where BN(k) is an error sequence. Although the image formation of g(n) in (7.2) and (7.7) are similar, their N-point DFT’s (7.4) and (7.9) are not the same because the length of g(n.) is now greater than N. This is the first problem which we encounter when we truncate the image g(n). The second error arises from circular convolution as follows: The finite aperture of the image acquisition device modulates the observed image. Since the aperture window w(n) is of size Chapter 7. Reduction of Boundary Artifacts in Image Restoration 94 M(< N), then the observed truncated image g (n) is 1 gi(n) = g(n) w(n) = gN(n).w(n) (7.10) where gi(n) is of length M(< N). To restore the observed g (n) we apply the restoration filter, hr(fl), of size L,. = M. To 1 apply hr(n), we pad each of gi(n) and hr(n) by zeros and use the N-point (N = 2M) DFT. In what follows we let G (k) and W(k) be the N-point DFT of the zero-padded gi(n) and 1 w(n) respectively. Since gN(n) is of length N(> M), the N-point DFT of (7.10) results in the circular convolution Gl(k)=GN(k)® W(k) where ® is the periodic (circular) convolution operator. (7.11) Thus, from (7.9) and (7.11) we conclude that in general (k) 1 G = The application of the inverse filter [H(k). FN(k) + EN(k)j ® W(k) (7.12) to G (k) in (7.12) will not result in E(k) = F(k) as in 1 (7.6) due to EN(k) and the circular convolution with W(k). Even if EN(k) = 0 (case 1 above), we will still have a circular convolution with W(k). Thus, it is clear from (7.12) that the exact recovery of F(k) via linear deconvolution is usually impossible. The classical approach to the boundary problem has been to employ different smoothing windows, i.e. multiplying g (n) in (7.10) by w(n)’s of different shapes so as to lessen the 1 effects of the abrupt boundaries and the literature has many examples of such windows [111]. Although this approach restricts the influence of the boundary effects, the results may still be objectionable. In Fig.7.3c we illustrate the use of a trapezoidal window. The same Wiener filter as in Fig.7.3b is employed to restore the observed noisy step edge. The original image was multiplied by a trapezoidal window such that a band of 16 pixels wide around the image was attenuated Chapter 7. Reduction of Boundary Artifacts in Image Restoration 95 linearly while the central region of 32 x 32 pixels were unaltered. Clearly the step edge itself becomes sharper, reversing the blurring effect of the camera; however the margins of the image are severely distorted. So far we have analyzed the sources of errors associated with The first problem is that of f(n) image truncation in restoration. being too long resulting in length of g(n) greater than N. When the N-point DFT is used this results in the model (7.9). The second problem is that of circular convolution (7.11) with W(k). In what follows we model the problem differently so that the effects of the second problem are avoided. Specifically let us select the M points {fo,fl,...,fM—1} from the long sequence f(7i) — and assume that the observed (n) (of size 1 g M) is due to the convolution of h(n) with {fo, plus the modeling errors. Let us pad the M selected points of f(n) ..., fM—1} with zeros to length N to obtain the set fi(n) thus f(n) — fi(n) = {...,f_2,f1,O,O,O,...,fM,fM±1,...} Padding g (n) with zeros to size N, the resulting 1 (m) 1 g = h(n) * fi(n) + ei(n) (7.13) gi(n) is for 0 < n < N — 1 (7.14) Taking the N-point DFT we get (k) 1 G = (k) 1 H(k) F(k) + E . (7.15) Chapter 7. Reduction of Boundary Artifacts in Image Restoration 96 Thus unlike (7.12), by using the present model utilizing fj(n) of (7.13), we now have a simple additive error sequence 1 E ( k). We show in appendix 7A that E ( 1 k) = h(j)W (f(_n) - (_1)kf(M - n)) W (7.16) for 0 < k < N-i, where WN = e_2111 In the space domain the error sequence ei(n) is ei(n) = I Jh(j)[f(n-j)-fi(n-j)] ç (n_j) 1 (_Jh(j)f for 0<n<M-i for M<n<N—i Please note that in the above we assumed h(n) to be in the form {h ,h 0 , 1 h(n) is symmetric and is represented in the form {h_, ..., ho, . .., h. (7.17) ..., } hL_}. If however then the equation (7.17) becomes e(n) (Ll)/ h 2 (J)[f(fl_J)_ f(n-j)] _(Ll)/ h 2 (i)fl(i) Applying the restoration filter Hr(k) E(k) for 0 <n < M- 1 for n<0 or n>M results in = = H(k).G(k) (7.18) = Fi(k)+Hr(k)Ei(k) (7.19) And the error in the restoration is E ( t k) = P(k) = H(k) — F ( 1 k) (7.20) E ( 1 k) (7.21) i.e. the error in restoration is given by (7.21) and (7.16). Chapter 7. Reduction of Boundary Artifacts in Image Restoration 7.2.2 97 Restoration by circular de-convolution Consider again the degraded image g(n). Truncating g(n) to length M to obtain the observed gi(n) and padding the latter by zeros creates severe edges within the signal 1 g ( n) which causes the errors in (7.12) or (7.14). Thus instead of assuming the observed gi(n) is a truncated signal let us assume 1 g ( n) as a period in a periodic signal and let us investigate the effects of the periodic or circular de-convolution of this period with the reconstruction filter hr(n). Let us denote two periods of g (n) as 2 1 (n) is thus the N-sample sequence 2 g ( n). g g ( 2 n) {go,...,gM—1,go,...,gM—i} = i.e. g ( 2 n) = g ( 1 n) + gi(n — M) 0 <n < N 1 — (7.22) This model implies that the object resulting in g2(n) (compare with (7.13)) is formed of the two periods f ( 2 n) = {fo, fi, ..., fM—i, fo, fi, ..., fM—i} (7.23) i.e. f(n) — f ( 2 n) = {..., f—a f-i — — fM-i, 0,0,0,..., f — fo, fM+i — fi, ...} (7.24) i.e. we assume 2 g ( n) is the result of the circular convolution of f (n) with h(n) plus errors 2 g2(n) = h(n.) f2(n) + e2(n) 0 n K N — 1 (7.25) where ®is the circular convolution operator. Since g was assumed periodic with period M, studying two periods instead of one will have no effect on the results of computations, but it facilitates the comparison of the error sequences which are now all of size N as in section 7.2.1. Applying the restoration filter h(n) whose length is Lr M we get f2(fl) = hr(fl) ® g ( 2 m) 0 n K N — 1 (7.26) Chapter 7. Reduction of Boundary Artifacts in Image Restoration Since g(n) was originally formed by linear convolution of f 98 and h, and f2 in (7.26) is produced by circular deconvolution, thus for the considered pixels, f2(n) 0 < n < M —‘ f(n) (7.27) in general due to the ‘wraparound’ effect of the periodic convolution. Let 2 G ( k) be the N-point DFT of 2 g ( n), then from (7.24) we obtain G ( 2 k) H(k). 2 F ( k) + E (k) 2 0 <m < N — 1 (7.28) and applying the inverse filter Hr(k) we obtain F ( 2 k) = Hr(k) G(k) 0 < . fl < N 1 (7.29) F ( 2 k) + Hr(k) E(k) (7.30) . and thus the error in the restoration is Ernv(k) = F ( 2 k) F ( 2 k) (7.31) = ll(k).E ( 2 k) (7.32) — i.e. it is exactly the same as (7.21) except that here the term 1 E ( k) is replaced by E (k). Thus 2 if 2 E ( k) < E (k) then we expect the error in restoration for this circular deconvolution case to 1 be less than that of the previous linear deconvolution case. Appendix 7B shows that E ( 2 k) for 0 < k < N e2(n) = — = h(j)W [f(-n) - f(M n)j [i + (_i)k] W (733) 1, and in the space domain we can show that 2 e ( n) is 2 1Zo’h(j)[f(n( n-j)1 j)-f for 0<n<M-1 ( for M<n<N—1 (n—M) 2 e (7.34) Comparing (7.33) and (7.16) we note that in this case the error 2 E ( k) is a function of the quantity b ( 2 n,k) = [f(-n) - f(M - n)j [1 + (_i)k] (7.35) Chapter 7. Reduction of Boundary Artifacts in Image Restoration 99 while in the case of linear convolution the error 1 E ( k) (7.16) is the same function as in (7.33) but of the quantity bi(n,k) = [f(_n) - (_1)kf(M n)] - (7.36) We restate (7.36) and (7.35) as (k) is a function of 1 E I f(—n)—f(M—n) forkeven I f(—n)+f(M—n) forkodd and (k) is a function of 2 E I 2[f(—n)--f(M—n)} 10 forkeven forkodd For an image of constant grey level, note that E(k) = 0 only when k is even, but E (k) 2 = 0 Vk. In this special case E (k) < E 2 (k) Vk. Also if the image is “globally stationary”, i.e. 1 Exp{f(n)} = K ‘v/n (7.37) where K, the image mean is a constant and Exp represents the expected value, then Exp{Ei(k)} 0 only when k is even, but not when k is odd. However from (36) Exp{E (k)} 2 case we also have E (k) 2 = 0 Vk. In this (k) Vk. The same results also hold when the grey levels around 1 E the boundary of the image are all equal to the same constant. An alternative way to see that is to note that e (n) in (7.34) is a function of f(n)—f 2 (n)which 2 is given in (7.24) and ei(n) is a function of since Exp{f (n)} or Exp{f(n) 1 However the Exp{f(n.) — — fi(m) or f(n) — fi(n) given in (7.13). In the latter fi(n)} are never equal to zero then Exp{e (n)} is never zero. 1 (n)} may be equal to zero when the areas surrounding the picture’s 2 f boundaries are of similar gray level values. For these cases we therefore expect the error in restoration to be smaller for the case of circular convolution than that of the linear convolution. In particular if the image is truly periodic with period M then E (k) 2 = 0 Vk i.e. the error vanishes when circular convolution is used. This fact suggests how structural features within the image can be exploited to reduce boundary truncation artifacts. = Chapter 7. Reduction of Boundary Artifacts in Image Restoration 100 In this section we have modeled the degradation process by a circular convolution process and therefore, as is well known, it leads to wraparound errors which are included in e (n). The 2 question now is how we can use these errors of the circular convolution to our advantage so that they lessen the effects of the boundary truncation artifacts. In the next section we propose a simple procedure which enables the use of periodic convolution and reduces the impact of boundary truncation artifacts. 7.3 Restoration by image extension and periodic convolution Our proposed solution to the boundary problem is to form a new image 3 g ( x, y) as in Fig.7.4. The new image is N x N(N — 2M) pixels. The top left quadrant of g (x, y) is the observed 3 M x M image g (x,y) and the other three quadrants are extensions of g 1 (x,y). 1 The top right quadrant is the mirror image of g (x, y) about the AA aids and the bottom half is the 1 mirror image of the top half about the BB axis. The method of extending an image at the boundaries albeit by reflecting only a few of the rows and columns has been earlier suggested for performing convolutions of the image with spatial masks (such as in smoothing or edge detection applications) [107]. The extension of the image by reflecting all its rows and columns and its application with circular deconvolution in restoration problems are new. If it is desired to completely remove the boundary effect in the most general case we need to apply signal prediction techniques to extend an M x M image by M/2 pixels on each side. This is clearly a non trivial task and our proposed method of image extension by mirror imaging effectively achieves a first order prediction in a straight forward manner. We have seen in section 7.2 that using circular convolution, via N x N DFT, results in smaller expected errors than using linear convolution when the image is globally stationary. We therefore use circular convolution here also, i.e. we shall not pad the newly formed N x N image by zeros before applying the restoration filter hr(,y) in the transform domain. Considering again the untruncated one-dimensional image g(n), g(n) = g(n). w(n) 0< n M —1 __ ___ ___ ___ ___ Chapter 7. Reduction of Boundary Artifacts in Image Restoration A’ AL B 11 A Figure 7.4: Extension of the image by mirroring in the x and y directions. ioi Chapter 7. Reduction of Boundary Artifacts in Image Restoration 102 The newly formed image in the one dimensional case is g ( 3 n) = (7.38) {go,gi,..,g_i,gM—i,...,go} i.e. g ( 3 n) = gi(n)+gi(N—n —1) 0 <n <N—i (7.39) We assume that the point spread function is symmetric and that the image g3(n) is formed by a circular convolution of the symmetric h(n) with 3 f ( n) where now 3 f ( n) is given by f ( 3 n) = {fo, fi, ..., fM—i, fM—i, ..., (7.40) fo} thus f(n)— (f 3 n) {...,f2 — fi,f-i — fo,0,0,0,...,fM — fM—i,fM+i — fM-2,...} (7.41) In the frequency domain, using N-point DFT and circular convolution we have G ( 3 k) H(k) 3 F ( k) + E (k) 3 = (7.42) . We now perform the circular convolution of the extended image 3 g ( n) with the reconstruction filter h(n) in the DFT domain ‘ ( 3 k) = (k) 3 Hr(k) G = F ( 3 k)+E(k) . 0 < k < N 1 (7.43) where E ” T (k) = (k) 3 H(k) E . (7.44) Note that (7.44) is the same as (7.21) and (7.32) except for the 3 E ( k) term. We have shown in the previous section that when the image is globally stationary the expected error using circular deconvolution is less than that using linear deconvolution. For the present case we show below that the expected error (using mirroring and circular convolution) Chapter 7. Reduction of Boundary Artifacts in Image Restoration 103 is less than that using linear deconvolution if the image is locally stationary. Local stationary is a more relaxed condition than global stationary [112]. We show in Appendices 7A and 7C that for symmetric h(n) the error term E (k) is given 1 by (L—1)/2 (k) 1 E > = M—1 —1 > h(j)W jzzl f(m)W n=—j (L—1)/2 j—i h(j)W + nM—j M—1+J f(n)W + — f(n)Wj (7.45) n=M for 0 < k < N — 1, while the error term E (k) is given by 3 (k) 3 E = Ei(k) + WE (—k) 1 (7.46) and hence (L—1)/2 (k) 3 E —1 {f(n) h(j)W = - f(—n — 1)}Wj n=—j j=1 M-1 {—f(n)+f(N—n—1)}Wj — nM-j (L—1)/2 -1 h(j)W + jzzl {—f(n)+f(—n—1)}Wj n=O M—1+j {f(n)—f(N—n—1)}Wj + (7.47) n=M for 0 < k < N e3(n) = f — 1. In the space domain e (n) is 3 Jh(j)[f(n-j)-f ( 3 n-j)] (N_n_1) 3 (e for 0<n<M-1 for M<n<N—1 (7.48) Comparing the expressions for E (k) (7.45) and E 1 (k) (7.47) we note that they have the same 3 form except that in E (k) every pixel value 3 f(n) has been replaced by the difference of two pixel values, where the two pixels lie in the same local neighbourhood. For example, inside the first term in (7.45) and (7.47), f(n) in (7.45) is replaced by the difference of f(n) and f(—n— 1) where the possible values of n lie between -1 and — and where L is the width of the blurring Chapter 7. Reduction of Boundary Artifacts in Image Restoration function. Thus the maximum distance between f(n) and f(—n — 104 1) is less than L. The same applies for the other 3 terms in (7.45) and (7.47). Thus when the image is locally stationary in the truncation boundaries then 3 Exp{E ( k)} = 0 Vk and thus Exp{Ea} < Exp{Ei} (7.49) We therefore conclude that our proposed method of extension of the observed image by reflection followed by circular deconvolution leads to smaller expected errors in the restored image when the truncated image is locally stationary at each of its boundaries. Also, as illustrated in the Fig. 7.5, if the grey levels in the neighbourhood of the boundary f(0), up to L,. pixels are all equal to a constant K , and those from M 1 boundary f(M — — L,. — 1 up to the 1) are also equal to another constant K , then E 2 (k) = 0 Vk. In the case of 3 circular deconvolution (section 7.2.2) E (k) = 0 under the added condition that K 2 1 = K . 2 In the two dimensional case, if the restoration function hr(fl, m) has L,. non-zero terms in both directions then the wrap-around effects due to circular convolution affects a band, (L,. — 1)-pixels wide, along the margins of the image. Referring to Fig. 7.6 we observe that if the objects of interest fall inside region A then they will be free from distortion due to wrap around, but the margins are not. If this margin is a smooth background region in which the energy is concentrated at a d.c. level, then the wrap-around errors introduced by the circular convolution operation is such that they will maintain this d.c. level exactly and the object is exactly recovered. We also note that if the support of the blur function is along one axis only (as for example in motion blur) then to obtain exact recovery, the grey-levels of the margins of different lines of g (x, y) need not be all equal; additionally, the grey levels may be different for the right and 1 left L,. 7.3.1 — 1 pixels of the same line. Computational Load For our proposed method (section 7.3) we use the N(= 2M)-point DFT. In the two-dimensional case we can show that the computational load compared to the circular convolution (i.e. without Chapter 7. Reduction of Boundary Artifacts in Image Restoration 105 b a 0 L r M-L r M M÷L r N-Lf 1 Figure 7.5: The extended signal under the conditions necessary for exact recovery. N-i Chapter 7. Reduction of Boundary Artifacts in Image Restoration 106 L? 1 A L1 Figure 7.6: Spatial support of the boundary artifacts around the margins of the image for the case of Lr <M. Chapter 7. Reduction of Boundary Artifacts in Image Restoration 107 any zero padding) will increase from O(M 2 M) to O(4M 2 log 2 log 2 2M). However since the resulting extended image g (n) is real and symmetric considerable gains can be made. Taking 3 N-point DFT of (7.39), after a spatial shift of one half a pixel, we get (k) 3 G = 2Re{Gi(k)} (7.50) M-1 = 2Re{ gi(m)e_32N} for 0< k < N — 1 (7.51) M-1 = 2 gi(n) cos(’ffkn/M) (7.52) Therefore we only need to compute the real DFT transform of an M-point signal. The trans formation (7.52) is the real part of the Discrete Fourier Transform and FFT algorithms can be adopted in its computation. For conventional linear deconvolution, to calculate G 1 (k) we need to calculate the real part as well as the imaginary part of G (k). The same applies to the inverse DFT computations. Our 1 proposed technique therefore requires half the computational load of the traditional approach. 7.4 Experiments Example: Consider a one-dimensional image M pixels long, g (n). We extend this image to N 1 2M pixels by reflecting it about a point half a pixel away from the last pixel to form g (N 1 We now form a periodic function 2 (n), one period of which is g (n) 2 the first L,. — 1 values of gi(n) be all equal to a, and that its last L gi(n) + g (N 1 — — — = n). n). Let 1 values be all equal to b. The signal g (n) is shown in Figure 7.5. 2 Assume that the reconstruction filter has Lr(<M-i) non-zero values , 1 , 0 h h ...,hrL_1 and that these values are normalized to have a sum of 1. Now perform a circular convolution of (n) and hr(n) to compute f(n). 2 g )(n)= a.hr(p)=aforO<n<Lr_1 b.hr(p)bforN_Lr+1nN_1 Chapter 7. Reduction of Boundary Artifacts in Image Restoration 108 Thus the aliasing error cancels the boundary truncation error and therefore the margins remain undistorted. In the above example we obtained the exact result because we made the aliases compensate for the truncation errors. To understand why this is so, let us consider the case of g(n) C, for 0 < n < 4, i.e. a constant gray level and h = = a, b, c. Using the linear deconvolution, the recovered image is f(n)=C{a,a+b,1,1,1,b+c,c} 0< n< M+L —1 where we have normalized the area a+b+c under h,. to be equal to 1. If we use image extension, we get the length of g = 2M = 10 and because of circular convolution the new f(n) is composed of the sum of the restoration results of the signal plus its aliases. Thus, 1(n) = C{a,a+ b,1,1,1,1,1,1,b+ c,c} + C{b+ c,c,0,0,0,0,0,0,0,0}+ C{0, 0, 0, 0, 0,0,0,0,a,a+ b} where the second and third quantities are due to the aliases to the left and right of the signal. Thus, )(n) C{1, 1, 1,1,1,1,1,1,1, 1} and after deleting the second half of f(n) we obtain )(n) = C{1, 1, 1,1, 1} Thus we have used the aliases to our advantage. Fig.7.3d illustrates the effect of this technique in restoring the image of the blurred step edge. In the general case, i.e. for any image, this method of image extension preserves a first order continuity at the image boundaries and therefore greatly reduces the undesirable effects of boundary discontinuity artifacts. The only requirement is that within each horizontal or vertical lines the margins should be reasonably smooth. Chapter 7. Reduction of Boundary Artifacts in Image Restoration 109 Figure 7.7: Restoration of the blurred image of Bayan: a) original; b) a 64x64 section cut from the original larger image; c) a 64x64 section of the original image blurred by the out-of-focus model of the digitizing camera; d) result of restoring by linear deconvolution; e) result of restoring by linear deconvolution alter windowing the data by a trapezoidal window; f) result of restoring using image extension and circular deconvolution; g) the difference image of b) and c); h) the difference image of b) and d); i) the difference image of b) and f). Chapter 7. Reduction of Boundary Artifacts in Image Restoration 110 To show how the new method performs on any image even if the boundaries are not smooth, we apply it on the 64x64 pixel image of a smiling child, Bayan. Fig.7.7a shows the 128x128 pixel image of Bayan and Fig.7.7b is a 64x64 section of his face. Fig.7.7c is the appropriate 64x64 section of the result of the linear convolution of Fig.7.7a with a 64x64 blurring mask h(n). The blur function of Fig.7.2 was used in this experiment. Fig.7.7d is the result of linear deconvolution of the inverse restoration filter with the blurred and zero-padded image of Fig.7.7c. The effects of the boundary artifacts are clearly visible. Fig.7.7e was obtained by first applying a two-dimensional trapezoidal window on Fig.7.7c in the space domain before zero-padding and linear convolution with the same restoration filter. Fig.7.7f is the result of the application of the restoration filter according to our proposed method. Figs.7.7g, 7.7h, and 7.7i, are the difference images of the original image in Fig.7.7b and Figs.7.7c, 7.7d, and 7.7f respectively. A bias of 128 grey-levels has been added for display purposes. Summary We have discussed the problem of image truncation in image deconvolution. The commonly practiced method of frequency domain image filtering firstly involves the application of a window (such as hamming, trapezoidal, etc.) in the space domain to smooth the effects of truncation. Since by its nature, DFT performs circular deconvolution (convolution), after applying the window, the image is zero-padded to the required DFT size so that linear deconvolution is achieved. This approach partially corrects for the boundary truncation artifacts but at the cost of reducing the useful part of the image. In this chapter we have exploited the natural circular deconvolution process to our ad vantage. We first show that if the blurred image is globally stationary, then using circular deconvolution (i.e. instead of padding with zeros we pad by the image itself) leads to reduced errors arising from boundary truncation effects. We have then proposed a simple image exten sion technique by which the observed image is mirrored in both spatial directions. We have Chapter 7. Reduction of Boundary Artifacts in Image Restoration 111 shown that when the blur is a symmetric function, applying circular deconvolution on the njir.. rored image leads to reduced errors when each area around the boundaries of the image is only locally stationary. Our proposed method thus leads to much improved restored images. This technique is also of advantage in the restoration of large blurred images and requires half the computational effort of the traditional linear deconvolutional approach. Chapter 7. Reduction of Boundary Artifacts in Image Restoration 112 Appendix 7A In this appendix we will derive an expression for the boundary truncation error associated with the use of linear convolution. Consider a long image f and its linear convolution with an L( M)-point unit sample response h to produce the blurred image g. We observe M points of g, denoting the observed sequence g 1 O<i<M-1 = We take the N(= 2M)-point DFT of g 1 G ( 1 k) gi(i)W = Ok<N—1 M—1 L—1 G ( 1 k) h(j)f(i—j)W = =O j=U L—1 M—1 [(i h(j)W = — j=O Now let n = i — j M-1-j L-1 G ( 1 k) h(j)W = j=O L—1 Gi(k) h(j)Wjjc = j=O L—i —i ( M—1 f(n)W + jzzO ( M—1 f(n)Wj — n0 nM—j 3 fM—i h(j)W = / f(n)Wj n=—j f(n)W + 3 f(—n)W \n=O G ( 1 k) f(M — — n)Wj” n=1 = H(k)F(k) + E (k) 1 where B ( 1 k) = h(j)W (f(_n) - (_1)kf(M - n)) W Chapter 7. Reduction of Boundary Artifacts in Image Restoration If we assume a symmetric blur function h(n) in the form {h_, 1 113 ..., ho, ..., - } then we 1 h have (L—1)/2 Ei(k) —1 M—i h(j)W = f(n)Wj— j=1 f(n)W n=—j (L—1)/2 j—1 M—1+, h(j)Wc —f(n)W+ + n=O 3=1 f(n)W n=M Appendix 7B In this appendix we will derive an expression for the boundary truncation error associated with the use of circular convolution. We will use two periods of 1 g ( n) in order to obtain an N-sample sequence. We have g ( 2 n) = {go,...,gM—1,go,...,gM—;} = gi(n)+gi(n—M) We also define the sequence f ( 2 n) {fo, = ..., fM—i, fo, ..., fM—1} fi(n)+fi(n—M) Taking N-point DFT we have G ( 2 k) = (k) 1 (1 + (_1)Ic)G and F ( 2 k) (1+(—l)’)F ( 1 k) In terms of the image formation model we have G ( 1 k) = (k)-t-E 1 H(k)F ( k) G ( 2 k) = H(k)F ( 2 k) + E (k) 2 Chapter 7. Reduction of Boundary Artifacts in Image Restoration 114 and hence (k) 2 E (k) 2 E = (k) 1 (1 + (_1)ic)E h(j)W [f(-n) - f(M - n)j [1 + (_i)k] W Chapter 7. Reduction of Boundary Artifacts in Image Restoration 115 Appendix 7C In this appendix we will derive an expression for the boundary truncation error associated with the use of image extension and circular convolution. The original image is assumed to be (n) 3 f = {fof1,...,fM—1,fM—1,...,fo} and the corresponding observed image after extension is (n) 3 g = {go,gl,...,gM—1,gM—1,...,go} We will assume a symmetric blur function h(n) in the form {h_L, ...h , 0 have (L—1)/2 g3(i) (i—j) 3 h(j)f = OiN—1 j=—(L—1)/2 We take the N(= 2M)-point DFT of g3 (k) 3 G (i)W 3 g 0 <k < N We also have (n) 3 g = gi(n)+gi(N--n—1) giving (k) 3 G = (k) 3 G Gi(k) + W/’G (—k) 1 H(k) 3 F ( k) + E (k) 3 . (k) 1 G = (k)-j-E 1 H(k).F ( k) (k) 3 E = E(k) We conclude that -j- (—k) 1 WE — 1 ...,hL_I} then we Chapter 7. Reduction of Boundary Artifacts in Image Restoration 116 The expression for 1 E ( k) for the case of a symmetric h(n) is given in appendix 7A above. Substituting for Ei(k) we have (L—1)/2 E ( 3 k) —1 h(j)W = f(n)Wj f(n)Wf — j1 n=M—j (L—1)/2 >Z + —1 M-—i 1 f(n)W ç 1 h(j)W f(n)Wk — j=rl n=M—j (L—1)/2 j—1 JkI—1-)-j h(j)W7 —f(n)W+ + f(n)W n=O (L—1)/2 n=M j—1 h(j)Wjc + Wi/c M—1+j f(n)W/ — j=1 + W f(n)W/ nO Substituting the dummy variable n’ for n + 1 in the second and fourth square bracket we get (L—1)/2 E ( 3 k) —1 h(j)W f(n)W — n=—j jzl (L—1)/2 0 M h(j)W/ + f(n’ jz1 — 1)W’Ic j—1 h(j)W + — 1)W7” nzzM-j-F1 n’=—j+l (L—1)/2 f(n’ — M—1--j f(n)Wj + — f(n)Wj j=1 j (L—1)/2 h(j)W + M+j f(m’ — — f(n’ 1)W” + j1 — 1)W’Ic nM+1 With new dummy substitutions n = N — n’ or n —n’ as appropriate in the same terms we obtain (L—1)/2 E ( 3 k) M—i —1 h(j)W f(n)Wj — f(n)Wj n=—j (L—1)/2 M—1+j —1 h(j)Wk + f(—n—1)Wj— jzrl (L—1)/2 + j—1 h(j)W? —>f(n)Wf+ (L—1)/2 —1 h(j)W] j=1 f(n)W n=O 3=1 + f(N—n—1)W nM f(—n—1)Wj+ — n=—j f(N—n_1)Wj n=M—j Chapter 7. Reduction of Boundary Artifacts in Image Restoration 117 We can now combine the terms in the first and fourth square brackets together, and similarly the terms in the second and third square brackets to obtain the final expression (L—1)/2 E ( 3 k) h(j)W = j=1 {f(n)—f(—n--1)}Wj n=—j M-1 — > {—f(n)+f(N—n—1)}W n=M—j (L—1)/2 k {—f(n)+f(—n—1)}Wi 1 h(j)W + n=O M-1+j + {f(n)_f(N_n_1)}W] Chapter 8 Segmentation This chapter describes several approaches to the segmentation of microcalcifications from the background parenchyma. Segmentation is a prerequisite step in the analysis of mammographic images. 8.1 Image Analysis Image analysis is concerned with the extraction of semantic descriptions from images. It differs from image processing in that its output is not another image. Image analysis techniques can be employed in detection and localization of abnormalities and in extracting quantitative features. Pattern recognition methods are then used for recognition and classification of abnormalities, as described in chapter 9. Common approaches in image analysis include matching, segmentation, shape analysis and description [113]. Matching normally involves template matching using cross correlation techniques or matched ifitering. These techniques may be applied in a transform domain to render the match invari ant to size, rotation or translation [94]. Template matching is a common technique in machine vision for segmentation for example of manufactured parts. Due to great variability in shapes and sizes of microcalcifications however, this approach is not likely to be useful for the detection of abnormalities. Segmentation methods involve thresholding, edge detection and texture based segmenta tion [1141. Thresholding may be applied to gray level or some other property such as first or second 118 Chapter 8. Segmentation 119 oider histogram [1151. The edge detection methods are well developed and are in large based on derivative operators or best fit models. The more successful algorithms model the computational aspects of the human visual system and lead to operators that extract edges at various scales or resolutions. Many such operators can be approximated by the difference of Gaussians [1161. Texture segmentation uses variations of split and merge algorithms. Also the concept of fractal dimension [117] is used to segment textured images. These algorithms may be applied to segment masses from a mammogram. They are not however of much use in the case of microcalcifications that are only a few pixels in size. Shape analysis employs a description of the object boundaries and may involve morpho logical processing such as thinning etc. [118]. These processes are followed by extraction of numerical features which are then employed to classify the observed abnormalities. 8.2 Detection and Segmentation of Microcalcification Clusters The first step in automated mammographic image analysis is the detection of calcifications, if present. The criteria commonly used by radiologists for recognition of calcifications must first be quantified. This was done based on a large number of cases in collaboration with an experienced radiologist. Features for recognition of calcifications include intensity level, local contrast measured by offset or ratio of the pixel to the average within a window, shape and size tests, gradient tests, etc. In this chapter we describe three methods for the detection and segmentation of microcal cification clusters. Mammograms will typically include identification marks. To exclude these from the analysis, a simple routine was applied to separate breast from non-breast areas. 8.3 Local Histogram Thresholding Algorithm Thresholding of image grey level histogram is commonly used for object segmentation. We implemented a modified version of a locally adaptive method described in [36]. Each image is Chapter 8. Segmentation 120 first subdivided into 100 square regions. In each region the grey level histogram is computed and the pixels belonging to the non breast areas are automatically excluded. The grey scale histogram is repeatedly smoothed using a simple three tap Finite Impulse Response (FIR) averaging filter until it has less than three modes. If the resulting smoothed histogram is bimodal and a significant valley is still found then a threshold value is chosen. The thresholds for regions whose histograms are unimodal or do not exhibit a clear valley are interpolated from neighbouring regions. In this way local thresholds are found, that vary from region to region but are constant within each region. The local threshold values that are found may be used in a preliminary segmentation of the image. The input image is a 32mm x 32mm section of a mammogram containing microcalci fications, shown in Figure 8.la and the processed output is given in Figure 8.lb. Clearly the choice of sub-region size and boundary affects the results and a large portion of the background breast structure is misclassified as calcifications. To correct for the effect of boundaries, a grid, as suggested by Davies et al. 1351, was placed in five overlapping positions as shown in Figure 8.2. The subregion of concern is the shaded area. In this way there are five independently determined local threshold values for each sub region. Figure 8.lc shows the result of accepting as calcification any pixel that has a gray level higher than any of its associated thresholds. Figures 8.ld, 8.le and 8.lf are the results of progressively raising the requirements of pixel gray level to be higher than 2, 3 or all five of threshold values obtained from the five grids. The amount of debris is clearly decreasing but few pixels which, when assessed visually, may be potential candidates are also eliminated. The resulting segmented objects are subjected to further tests involving size of the area and gradient values. The size test was imposed using the heuristic knowledge that significant microcalcifications are rarely greater than 1.6mm x 1.6mm in area. The introduction of this constraint removes most of the remaining artifacts as shown in Figure 8.3e. Figures 8.3a-d are reproduced from Figure 8.1 for ease of comparison. To further reduce the false positive rate, the following gradient test is applied. The Roberts Chapter 8. Segmentation 121 (a) (b) (c) (d) (e) (f) Figure 8.1: Segmentation using the Local Histogram Thresholding algorithm 122 Chapter 8. Segmentation Figure 8.2: The subimage grid displacement Chapter 8. Segmentation 123 (a) (b) (c) (d) (e) (f) Figure 8.3: Effects of area and gradient tests on segmentation Chapter 8. Segmentation 124 gradient operator for one edge of the pixel F(i,j) is given as: R=R+R R= IF(i,j)—F(i—l,j— 1)1 = F(i — 1,j) — F(i,j — 1)1 For computational simplicity we use a modified operator R = IRI + IRI. The average edge strengths of each pixel in the perimeter of each object is measured and weak edges are removed. Figure 8.3f shows that the application of this edge constraint removes several scattered de bris. It cannot be known whether the remaining border line cases are truly microcalcifications or not without a correlated pathologic study. Figure 8.4 shows the results of processing two other mammograms. It can be seen that the algorithm produces satisfactory results when calcifica tions are of good contrast. However when calcifications appear superimposed on other densities, the algorithm segments the entire mass with its associated calcifications. The algorithm fails to correctly detect or segment low contrast microcalcifications. 8.3.1 Object Labeling The above algorithm is computationally expensive due to the repeated shifting of the subregion forming grid. To reduce the computational effort, an alternative method which uses a two pass approach is taken. In the first pass the grid is fixed and the local histogram is used to determine a threshold for each pixel as described above. Potential microcalcification pixels are marked. Each object is then labeled using 6-connected neighbours then object boundaries are marked but they are not segmented from the background. These object boundaries will involve discontinuities at the edges of each subregion since the thresholds change at these locations. This is corrected for in the second pass as follows. The mean grey level of each object is computed and a unique threshold is allocated to each object. This threshold value lies halfway Chapter 8. Segmentation 125 (a) (b) (c) (d) (e) (f) Figure 8.4: Examples of segmentation using the Locai Histogram Thresholding algorithm Chapter 8. Segmentation 126 between the object mean grey level and the background and is used for extracting the object. This method produces results that are comparable with the previous method, but with lower computational effort. 8.4 Edge Detection and Region Growing Algorithm The above algorithm fails to correctly detect microcalcifications of low contrast. Thus another algorithm which identifies the pixels that may potentially belong to microcaicifications, using edge detection technique, instead of obtaining thresholds from the local histograms, is employed. Eight different edge detecting algorithms were implemented for the gradient test and their per formance compared. These included the following four gradient-based edge detectors, namely the Roberts, Sobel, Kirsh, and Prewit operators. These operators are simple to implement, and they have short execution times. However they are all sensitive to noise with comparable performances. Of these four convolutional masks, the Roberts operator gives the thinnest, i.e. the best localized edges, while the Kirsh operator gives the thickest edges with about 4 pixels marked for a one-pixel-wide edge. Two approximations to the second order derivative Laplacian operator were also imple mented. The edges are located by identifying the zero-crossings of the output from the Lapla cian operators. These operators are quite sensitive to noise with many spurious edges detected. To overcome this problem, we implemented the Marr-Hildreth edge detector which is the Lapla clan of the Gaussian operator with variable width of the Gaussian smoother. We found this edge detector to be essentially useless for this application. The reason is that microcalcifica tions are small objects with pronounced edges. The Gaussian smoother tends to eliminate the smaller and fainter microcalcifications and lead to many false negatives. The last edge detector that we tested was based on morphological operations. Essentially if an eroded version of the image is subtracted from the original scene, the resulting output may be thresholded to give an edge map. We found this morphological edge detector to be quite useful for larger microcalcifications, however single pixel calcifications are once again eliminated Chapter 8. Segmentation 127 by the erosion operation despite their high contrast. For these reasons we selected the Roberts operator as the edge detector of choice for the initial detection of edges. Spurious outputs from this operator were then eliminated in the subsequent steps of the algorithm. After the identification of boundary pixels, region growing techniques were employed to grow the calcifications. A local threshold was calculated based on the grey levels of each edge pixel, its neighbouring pixels, and the direction of the steepest gradient. This threshold was used in growing additional pixels at the boundaries of each object. If the seed pixel was not a true edge of a connected object, the object normally grows to be unrealistically large. This criteria was used to eliminate false detections. Finally, the resulting segmented objects were subjected to tests involving shape, size and gradient at the boundaries. Figure 8.5 shows sections of three mammograms containing clusters of microcalcifications and the output of the various segmentation routines. The left column is the original digitized mammograms. The top image is a degenerating fibroadenoma with benign microcalcifications. The two lower images are malignant formations. The central column is the output of the local area thresholding algorithm (algorithm 1), and the right column is the output of the edge detection and region growing algorithm (algorithm 2). It can be seen that in all the three examples the second algorithm outperforms the first one. 8.5 Neural Networks Techniques Various neural network architectures have been used in object detection and image segmentation problems. The basic idea is to treat image segmentation as a pixel classification problem. If m distinct regions exist in an image, then each pixel can have one of m different labels. The network dynamics, in a supervised, or unsupervised fashion, should then lead the network to a stable state such that each pixel has the desired label. In [1191 a constraint satisfaction neural network (CSNN) is developed and applied to segmentation of CT, PET, and MRI images. In [31] an algorithmic approach is taken towards detection of clusters of microcalcifications in digitized mammograms. A multilayered Perceptron is then used to reduce the false positive Chapter 8. Segmentation Figure 8.5: Comparison of two segmentation algorithms 128 Chapter 8. Segmentation 129 rate and increase the specificity of the detection algorithm with only moderate decrease in its sensitivity. In [32] the performance of three different neural networks are compared in reduction of false positive detection of microcalcification clusters. The input to the network consisted of seven features extracted from the segmented objects, and the network output would indicate if they represent a true cluster. In this section we examine the ability of two neural networks to operate directly on the raw image and segment the microcalcifications by assigning one of two labels to each pixel. 8.5.1 The Modified Hopfield Network: The first neural network considered is a modified version of a Hopfield network. This network has self organizing properties that favor the formation of compact regions and therefore a training set of images is not required. The derivation of this network is given in [120] in detail. We provide here a summary description of it. One neuron is assigned to each pixel in a digitized mammogram. Each neuron is connected only locally to the neurons of the pixels in a pre-defined neighbourhood. This neighbourhood may be defined arbitrarily as either the four connected adjacent pixels the eight connected pixels ( referred to as N 3 ). ( referred to as N 2 ) or The choice of neighbourhood is treated as a design parameter and will have an impact on the performance of the system. The input U, to each neuron i, comes from two sources, namely the local neighbourhood, and a bias input. Thus the input U is such that (8.1) where W, are the synaptic connections between neuron i and neuron neuron j, j, 3 is the output of V and I is the bias input to neuron i. We let the bias input be a normalized version of the pixel grey levels in the raw input mammogram, so that I=2—1 (8.2) Chapter 8. Segmentation 130 where g is the grey level of the pixel i in the observed image, and L is the maximum number of grey levels. We will assume symmetric weights, i.e. W, We constrain W, to be in the range { 0 , 1 = }. W, (8.3) The choice of values of W dictate the influence of the neighbourhood. For example if we set Wj = 0 for all i and j, all neurons are decoupled from each other and no segmentation takes place beyond the initial estimate. If W the closest 8 neighbours, and W = 1 for 0 for all other neurons, then these neighbouring pixels will have a strong effect on the classification of the central pixel under consideration while all others will only have an indirect effect. This is similar to a morphological operation that fills holes in the segmentation mask. The energy function for this network also has two components as follows: If a pixel belongs to the background (and not microcalcifications) then there is a high probability that its neighbours also belong to the background. Therefore if a pair of adjacent pixels has similar values then the energy contribution of this pair should be small. One of the terms in the energy function therefore is Also if the grey level of a pixel is very close to that of the background parenchymal pattern (i.e. relatively darker), or to that of the object (i.e. relatively brighter than its neighbourhood), then it is highly likely that in the stable state the corresponding neuron will be labeled as background or object respectively. Since we have chosen the input bias to be proportional to the image grey levels, then the product I,V. should contribute less towards the total energy value. The second part of the energy expression therefore is Finally, the network energy is given by E = - I V 1 . (8.4) Chapter 8. Segmentation 131 The Discrete Model: In this formulation we use bi-state neurons. Each neuron i, has two possible states, correspond ing to the two labels of its associated pixel. The neural output V, may take the values of +1, or 4, corresponding to the ON or OFF states of the neuron. The output of each neuron is derived from its input by a simple threshold logic 1+1 if U>O i_i if U<O, vt= where S. is a local estimate of threshold. A change in the output of the neuron i, tiV will result in a change in the energy of the network LE such that — - (>z wjvj (8.5) + In the above equation it can be seen that /V has the same sign as the expression in the bracket and therefore < 0. Hence the network dynamics tends to obtain the minimum energy. The Continuous Model: In this formulation the neural response is graded and is an approximation of the sigmoidal function. This function is similar to the ‘S’ function used in fuzzy sets and can be written as —1 if U—1 (U+1) — 2 1 if —1<U<0 ) 2 1—(1—U if 0<U<1 +1 if U>1 Once again the network time evolution leads it to a minimum energy state. The dynamics of the network are governed by simultaneous solutions to the differential equations 8EdU 1 dt ÔV. 86 Chapter 8. Segmentation 132 By Euler approximation we have U(t + At) = U(t) + At ( WVj(t) + I — U ( 2 t)) (8.7) It can be shown that the energy function of such a system is monotonically decreasing and therefore the network dynamics lead the system to a stable output representing the required binarized mask of pixels belonging to the microcalcifications. 8.5.2 The Feed Forward Neural Network: The second neural network is a three layer Perceptron. The input layer receives the grey level values of each pixel and the output layer flags the presence of microcalcifications. A hidden layer is added to ensure that a piecewise linear discrimination between the two classes may also be accommodated. The back propagation training a’gorithm is used for the adjustment of the synaptic connections and the training and test set of images are kept separate. It is well known that this network architecture does not have shift-invariant properties. It follows that the network needs to be trained on all the shifted versions of a given input pattern. This requires a prohibitively long training time, and it is unlikely that successful training can be accomplished for images of realistic size and complexity. To overcome this restriction we formulated the problem differently as follows. The detection of microcalcifications is essentially a local operation. A single microcalcifica tion is detected by a human observer by comparison of its features with those of its immediate vicinity. We can therefore examine a window around each pixel to decide whether it should be labeled as a microcalcification or not. In this way the network becomes shift-invariant. Based on the above analysis we chose a three layer perceptron with 81 neurons in the input layer, 4 neurons in the middle layer, and one neuron in the output layer. The input consists of grey levels in a 9 x 9 window centered on a pixel and the output is the label assigned to this pixel. Sections from four mammograms containing numerous microcalcifications of various sizes and shapes were chosen as the training set. training pattern vectors were generated by a moving 9 x 9 window. Chapter 8. Segmentation 133 Automated segmentation Object Background Manual Object 364 161 segmentation Background 371 101504 Table 8.1: Classification matrix for algorithm 1 8.6 Results and Discussions Performance of each of the two neural networks and the two above algorithms were tested using selected images from our database of 68 images. Local Histogram Thresholding algorithm is referred to as algorithm 1, and Edge Detection and Region Growing algorithm is referred to as algorithm 2. We have chosen to compare the results on a pixel by pixel basis instead of on an object by object basis, since both the detection and accurate segmentation of the shape of microcalcifications are clinically significant. Table 8.1 gives a typical classification matrix for the first algorithm. In this particular case a 320 x 320 pixel square window from a mammogram was segmented manually and also using the first algorithm. A total of 525 pixels were labeled manually as belonging to microcaicifications, and the balance of 101875 pixels as background parenchyma. The algorithm correctly labeled 364 pixels as being microcalcifications, missing the other 161 pixels. It also falsely labeled 371 pixels as being microcalcifications. In this image there were 11 individual calcifications, nine of which were detected by the algorithm. The shape of the calcifications however were not segmented accurately leading to pixel misciassifications. We use the standard definitions of Sensitivity TP TP+FN (8.8) Specificity TN TN+FP (8.9) and = and accuracy = TP+TN TP+TN+FP+FN (8.10) Chapter 8. Segmentation Automated segmentation Object Background 134 Manual Object 507 18 segmentation Background 26 101849 Table 8.2: Classification matrix for algorithm 2 where TP, TN, FP, FN are the number of true positive, true negative, false positive and false negative cases respectively. For this typical case Specificity is 0.996, Sensitivity is 0.69, and accuracy is 0.995. The high value of accuracy is due to the very large number of background pixels. In such cases care should be exercised in interpretation of the results. The classification matrix for our second algorithm and the same input image is given in Table 8.2. 507 of the 525 microcalcification pixels were correctly labeled and 26 pixels of the background were misclassified. We note from this table that both Sensitivity and Specificity have increased considerably. This second algorithm has consistently performed better than the first, indicating the superiority of an edge based technique to thresholding of local histogram. The modified Hopfield network was also tested with the same input image. No manual segmentation was performed as this network has self organizing properties. We found that this network can be tuned for a given set of input images to achieve acceptable segmentation. In particular the network dynamics are such that segmented microcalcifications are morpho logically compact, i.e. without holes in them. Unfortunately, the choice of neural response, the various gain factors, and the choice of neighbourhood size and shape are problem depen dent. This limits the application and generalization of this network to other data sets. Further investigation into the proper choice of these parameters is required. The multilayer Perceptron was tested with the same image data base. We first verified the shift invariant property of the network by training it on a subset of images and testing its performance on shifted versions of the same images. A typical classification matrix for this network is given in Table 8.3. We notice a comparable performance with the above algorithmic Chapter 8. Segmentation Automated segmentation Object Background 135 Manual Object 350 175 segmentation Background 87 101788 Table 8.3: Classification matrix for Perceptron approaches. Using the proposed techniques all of the microcalcification clusters in the original set of mammograms were correctly detected. Additionally, other minute calcifications were segmented which were not readily visible on the original mammograms. These results indicate that both algorithmic and neural network based approaches can segment microcalcifications from the background parenchyrnal pattern, and that a combination of approaches may be used for improved sensitivity and specificity. These techniques could prove of assistance for objective classification of some malignant breast lesions. The computational complexity of these routines depends on the contents of each image. In our experience all calculations can be performed within a few seconds on an Apollo 4500 workstation. This time is only slightly longer than image digitization and is sufficiently rapid to make this system usable in a clinical setting to provide a second (objective) opinion to the radiologist’s interpretation of the mammogram. 8.7 Summary Automatic detection and segmentation of microcalcifications may be achieved by application of algorithmic techniques or by use of artificial neural networks. We have developed two al gorithmic approaches using firstly a local area thresholding method and secondly a technique based on edge detection followed by region growing to segment microcalcifications from the background parenchyma. We also selected two neural network architectures and implemented object detection techniques on them. Further In this chapter we reported on the performance Chapter 8. Segmentation 136 and comparative merits of each one of these four systems. Of these four systems the approach based on a modified Hopfield network has been the least successful one in that it requires manual specification of parameters which are dependent on the contents of each image. The perceptron performs similar to the first algorithm. The second algorithm, based on edge detec tion and region growing, outperforms the first algorithm, which was based on local histogram thresholding. The second algorithm was therefore chosen for use in the next chapter. The objective of the work in this chapter and the following chapter is to develop a system for the automatic recognition of microcalcifications which are the most common signs of abnor mality in a mammogram. At present, no practical system is available to provide radiologists means for the computer assisted detection and quantitative analysis of mammographic images which will likely lead to improved diagnosis. Chapter 9 Classification Microcalcifications are non-specific indicators of breast carcinoma in conventional mammo graphic examinations. Currently, radiologists inspect mammograms with a viewing box and based on a subjective appearance of microcalcifications, determine the presence or absence of a suspicious lesion. The objective of the work described in this chapter is to extract various numerical features from microcalcifications in a digitized mammographic image and to use this information to construct a piecewise linear discriminant function for classification of images into two groups of clearly benign or possibly malignant class. Recent advances in digital image acquisition and analysis have made it possible to consider automated diagnosis of mammographic images. Towards this end, a good selection of small number of features would be helpful for a successful discrimination of mammograms [121]. We therefore evaluated over 100 numerical features extracted from individual as well as clusters of microcalcifications in digitized mammographic images. These features quantify the number, size, shape, roughness, and configuration of clusters of microcalcifications. The features are then examined individually and also in various combinations to test their potential in discriminating between benign and malignant microcalcification patterns. Pattern Recognition Methods Pattern recognition methodologies can be divided into two broad categories of supervised and unsupervised classification. Supervised learning may be further divided into statistical decision making and syntactic or structural. 137 Chapter 9. Classification 138 The image represents a point in N-dimensional feature space. In the training phase the feature vector is considered as the input and the desired classification as the output. A dis criminant function is then computed to best separate the feature space into various classes. This function may be linear, quadratic, or generally non-linear. Piecewise linear discriminant functions are the most popular in image pattern recognition since most practical feature vec tor spaces are not linearly separable and the general non-linear discriminant functions are not mathematically tractable. In computing the discriininant function, if the a priori knowledge of the class statistics is unobtainable then a heuristic distribution-free approach is followed. One such decision rule is based on the minimum distance of the observed vector from the class centroid. Alternatively, k-nearest-neighbour classifiers may be used. The statistical pattern recognition methods commonly use Bayes’ rule and require knowledge of a priori probabilities of each class. These may be parametric such as Gaussian distributions or non-parametric. Unsupervised learning techniques such as clustering [122], and fuzzy set reasoning [1231 have also been used in pattern recognition tasks. 9.1 Data Base We have examined over 400 mammograms from patients who have recently undergone breast biopsies at Vancouver Clinic of British Columbia Cancer Agency (BCCA). From this group 68 typical images were selected. These images contain isolated fine and coarse calcifications as well as clusters of microcalcifications that may be clearly benign, clearly malignant or suspicious of malignancy. Participating radiologists from the Clinic have performed a blind diagnosis of the images without access to other relevant clinical or pathological information. The images were then digitized with 100gm sampling interval in both spatial directions using the two dimensional CCD device described in chapter 3. Each image was processed by the “edge detection and region growing” segmentation routine described in chapter 8. Each image was associated with a binary mask representing the pixels that belong to microcalcifications. Each image and its associated mask was checked manually Chapter 9. Classification 139 and if necessary, the segmentation mask was corrected. This step was included to ensure that any errors in automatic segmentation do not affect the pattern classification stage. Quantitative features were then extracted from each image. 9.2 Features My approach was to extract features from microcalcifications individually as well as collectively from all the microcalcifications present in the image. Screen-film mammograms were imaged and discretized to 1280x1024 pixels with 12 bit accuracy using the system as described in chapter 3. Each image was then subjected to a segmentation algorithm based on edge detection and region growing as explained in chapter 8. A binary mask created in this way would tag the pixels that were calcified. In order to eliminate the effect of errors in automated segmentation on the classification results, each binary mask was examined and corrected if necessary. For each microcaicification the following groups of features were computed: 1. Photometric variables: (a) Mean and variance of intensity of the object (microcalcification); (b) Contrast of the object over the background parenchyma; (c) The Optical Density (OD) of each pixel was measured as OD(x,y) = T(x,y) 10 —20log (9.1) where T is the local optical transmittance of the mammogram for each pixel at coordinates (x, y). (d) A histogram of OD was constructed for each object and minimum, maximum, range, mean and variance of OD were calculated. (e) The skewness of the OD histogram is the third moment of its distribution normalized by its second moment; and was calculated from (OD(x, y) ODmean) 3 (n 1)ODa,. — ODskew = — (9.2) Chapter 9. Classification 140 where n is the total number of the pixels in the object, ODva,. is the variance of OD normalized by ODean, and ODskew is a measure of asymmetry in OD distribution. (f) The long-tailed nature of the optical density distribution was measured by its kur tosis, being its fourth moment normalized by the square of its second moment (OD(x, y) ODrnean) 4 ( flOD I var — ODkt = (9.3) — 2. Size variables: (a) Area, being n, the number of connected pixels forming the microcalcification. (b) Perimeter, using 4-connectedness (i.e. only edge adjacent pixels are considered neigh bours) and appropriate corrections for square tessellations P = ni + v’irn 2 + 2.0fl3 (9.4) where n , n 1 , and n 2 3 are the number of edge pixels in the object with 1, 2, and 3 non-object neighbours respectively. (c) Mean radius, being the distance from the object centroid to the object edge. 3. Shape variables: (a) Compactness, calculated as the normalized ratio of area to perimeter squared; C = (9.5) where n is the total number of the pixels in the object. (b) Moment of inertia of the mask considered as a uniform laminar object, normalized by the mask area squared; Inert 2rd ( 2 x,y) where d(x, y) is the distance of the pixel at (x, y) from the object centroid. (9.6) Chapter 9. Classification 41 (c) Variance of the object radii; (d) Sphericity, defined as the ratio of the radii of two concentric circles enclosing the object boundary; (e) Eccentricity, defined as the ratio of the eigenvalues of the matrix of second moments of the object; i.e. let M 2 be such a matrix 2 M lyx where I, = I = — — ) Iyy and I and I, are variances of x-coordinate and y-coordinate values. If eigenvalues of M 2 are denoted as A 1 and A 2 and A 1 2 A then Eccent = 2 / 1 A A (9.7) (f) Elongation, defined as the ratio of major axis to minor axis of the least squares best fit of the object boundary to an effipse. 4. Roughness variables: coarse and fine roughness measured as energy content of se’ected frequency bands in the spectrum of the object radii. The coordinates (x, y) of the object boundary are first expressed in polar coordinates (r, 0). Using Fourier coefficients we can write r(8) = + (a cos(iO) + /Jb sin(iO)) (9.8) then the mean radius is ; (a ,b 1 ) determine the least squares best fit of the object 1 boundary to a circle; and (a ,b 2 ) determine the least squares best fit to an effipse. The 2 major and minor axis of this ellipse are given by 0 + 2a + b a (9.9) Finally the energy content within a given frequency band i 1 to i 2 is calculated from (a + b) (9.10) Chapter 9. Classification 142 We used the parameters (i 1 1 (i = 11 i 2 = = 2 3 i = 10) to measure coarse boundary variation and 31) to estimate fine boundary variation. Features of all individual microcalcifications in an image were combined to derive overall characteristic features for each image. Thus minimum, maximum, range, mean, variance and normalized standard deviation were calculated for each feature of individual microcalcifications. The center of mass of all masks within an image was identified and simple statistical data were computed for distance of each object to the cluster center. To estimate the degree of scatteredness of calcifications, the object-to-cluster-center distance was weighted by object area or object perimeter and similar statistical metrics were calculated. The mean and variance of the distance between pairs of microcalcifications were used as other features to estimate scatteredness. Additionally, a convex polygon enclosing each microcalcification cluster was calculated and several features relating to its size and shape were computed. Other features measured relate to the degree of alignment of calcifications. Radiologists often consider a cluster as being suspicious if the major axis of individual microcalcifications are aligned in a linear or branching fashion. We used two methods of calculating alignment as follows [59]: 1. For each object i, we measure the degree to which the major axis of all other objects in the scene are aligned with it, Alin = (1 — A— -.—)A?, M, . — (9.11) where k is the total number of microcalcifications in the image, Ar is the aspect ratio, A is a unit vector along the major axis of each object, and M. M 2 implies a dot product. . Note that long and thin objects are weighted more than round and compact objects. Simple statistical data from the distribution of alignment values were used to describe the whole cluster. Chapter 9. Classification 143 2. Alternatively, the direction of the major axis of the convex polygon enclosing the cluster was first identified. Orientation of each microcaicification was then compared to this vector to derive alignment features. Many other common features from the literature [124] (e.g. normalized central moments, Fourier coefficients, and other invariant features for shape recognition) were also implemented. A complete list of features is given in Appendix C. 9.3 Analysis A commercially available statistical analysis package (Bio-Medical Data Processing) was used to evaluate the feature vectors. The Fisher statistics for each individual feature was calculated. Individual histograms were also computed for each feature. The features were grouped together based on what attribute of the cluster was being mea sured. From each group of the features one or two representative ones were selected based on their discriminating power. These features were used in calculation of a piece wise linear discriminant function. The jackknife method was used to ensure that the training set and the test cases were kept separate. 9.4 Results Based on a data base of 68 images we have determined the most effective features when used individually as listed in Table 9.1 with their relative discriminating value. Nineteen features showed F-values higher than 4 and were considered significant in classification. It can be seen from Table 9.1 that the most important featllre is the number of microcalcifications in a cluster. This is confirmed by current practice of radiologists who consider the presence of more than 10 microcalcifications per square centimeter as highly suspicious for malignancy. Various features measuring alignment, compactness and inertia were also found to be significant. Combining a few selected features resulted in a two-class discriminant function that correctly classified 64 of the 68 selected mammograms, i.e. with an overall accuracy of 94.1% Chapter 9. Classification 144 Feature Number of Microcalcifications Range of Alignmentl Max. of Alignmentl Var. of Alignmentl Mean of Alignmentl Minimum Inertia Minimum Compactness Range of Alignmentll Max. of Alignmentll Mean of Alignmentll Normalized Std. Dev. of Area Mm. Mean Radius Mm. Perimeter Mm. Distance to Center Var. of Alignmentll Normalized Std. Dev. of Perimeter Mm. Pair Distance Average Mean Radius Average Pair Distance Fisher Value 73.30 46.94 44.57 27.11 26.57 21.95 16.47 11.8 11.52 10.61 7.06 6.92 6.6 5.36 5.28 4.9 4.87 4.35 4.16 Classification Accuracy 89.7 91.2 91.2 88.2 85.3 79.4 79.4 82.4 80.9 80.9 69.1 67.7 64.7 55.9 85.3 66.2 48.5 51.5 66.2 Table 9.1: Discriminant Power of Features (%) Chapter 9. Classification 145 The three “best” features taken as a combination are: X3 = Number of microcalcifications in the cluster X75 = Area of the convex polygon enclosing the cluster X59 = Variance of mutual alignment of microcalcifications X3, the number of microcalcifications in the cluster, is the feature with the most discrim inating power. The next best feature, when taken together with X3, is X75, which measures the compactness of the cluster. The third best feature, when taken together with X3 and X75, is X59, which measures the degree of mutual alignment of the individual microcalcifications within a cluster. Addition of other features did not significantly improve the classification re sults. The automatic selection of these features by the software package for the calculation of discriminant function, confirms and quantifies the experience of radiologists in separating the benign and malignant classes. The resulting discriminant function for the benign and malignant groups are: = 0.57932X3 = 1.78940X3 — — 0.00004X59 0.00028X59 — — 6.10094X75 — 15.05314X75 1.92331 — 11.96746 (9.12) (9.13) The classification matrix for the entire data set is given in Table 9.2. This table shows that even with the use of over 100 features the two classes in the training set can not be separated with 100% accuracy. 53 out of 55 benign cases and 11 of 13 malignant cases were automatically recognized. Given the fact that in current practice only two out of five biopsies prove to be malignant, the sensitivity (85%), specificity(96%), and accuracy (94.1%) figures derived from Table 9.2 represent a major improvement over current clinical practice. Table 9.3 represents the jackknife classification matrix for this data set. The performance accuracy of the algorithms drops to 92.6%. Chapter 9. Classification 146 Algorithmic Benign Malignant [ Benign“Truth” Malignant L 53 2 2 11 Table 9.2: Classification performance of computer algorithms on training images. Algorithmic Benign Malignant “Truth” Benign Malignant 53 2 3 10 Table 9.3: Jackknife classification performance of computer algorithms. 9.5 Summary We evaluated the potential of over 100 features in discriminating benign from malignant mi crocalcification clusters. The results indicate that a few features, when taken in combination, are capable of successfully discriminating the two classes of mammograms, with a success rate of over 92%. The computational complexity of feature calculation routines depends on the contents of each image, and in our experience all calculations can be performed within a few seconds on an Apollo 4500 workstation. This time is only slightly longer than image digitiza tion and is sufficiently rapid to make this system usable in a clinical setting to provide a second (objective) opinion to the radiologist’s interpretation of the mammogram. Chapter 10 Conclusions, and recommendations for future work 10.1 Objectives The objective of this research has been to investigate the following two hypotheses: 1. Using image restoration techniques in digitized mammograms, it is possible to improve the visibility of minute microcalcifications which are common signs of early breast cancer. In some forms of breast cancer these microcalcifications may be present in the vicinity of the primary tumors and are believed to be highly prognostic. 2. Using quantitative image features extracted from microcalcifications it is possible to clas sify a cluster as benign or madignant. 10.2 Summary of the work In order to carry out this investigation a series of hardware and software tools were developed. Chapter 3 describes the development and characterization of a film digitization device based on a two dimensional CCD array sensor. It is shown that this device can provide the required spatial and photometric resolution necessary in digitizing mammograms. We have described the derivations of system parameters and noise characteristics and have also implemented measures to reproduce the original image in the digital form with a high degree of fidelity. The distinguishing features of the newly developed system are: i) a fast method of digitizing mammograms and ii) acquisition of images with a high spatial and photometric resolution. Each 147 Chapter 10. Conclusions, and recommendations for future work 148 frame contains over 1.3 million pixels digitized to 12 bits per pixel. The sampling interval is 6.8 ,um without optical magnification. The fixed pattern and the random noise are minimized using background subtraction and signal averaging techniques. The resulting image appears comparable with that obtained by a laser-scanning microdensitometer obtained at a fraction of the time required. Chapter 4 reports on the clinical evaluation of digitized images and concludes that the use of the conventional magnification mammography procedure may be substantially reduced. It is shown in chapter 5 that application of known restoration techniques to digitized mam mograms can greatly increase the visibility of image details. Locally adaptive filters are used in chapter 6 to improve the results of image restoration in the presence of signal dependent radiographic noise. The results of image restoration algorithms on the data base of 30 images show a marked improvement in detectability of smallest particles of microcalcifications when judged by a human observer. In chapter 7, we have given mathematical derivations to show the benefits of using image extension and circular deconvolution in frequency domain image restoration via the DFT. To test the second hypothesis, we have developed image segmentation routines, including the evaluation of neural network techniques, for the automated detection and extraction of microcalcifications. Automatic detection and segmentation of microcalcifications may be achieved by the ap plication of algorithmic techniques or by the use of artificial neural networks. We selected two neural network architectures and implemented object detection techniques on them. The first neural network considered is a modified version of a Hopfield network. One neuron is assigned to each pixel in a digitized mammogram. Each neuron is connected only ally to a pre-defined neighbourhood. The network dynamics favor the formation of compact regions and lead the system to a stable output representing the required binarized mask of pixels belonging to the microcalcifications. Chapter 10. Conclusions, and recommendations for future work 149 The second neural network is the classical three layer perceptron with error back propagation scheme for training. The input layer receives the gray level values of the normalized image and the output layer flags the presence of microcalcifications. Further we have developed two algorithmic approaches to segment microcalcifications. In the first algorithm, thresholding of local image grey level histogram is used for object segmenta tion. In the first pass, each object is labeled and object boundaries are marked but they are not segmented from the background. In the second pass, the discontinuities due to region bound aries are corrected for by allocating a unique threshold value for each object commensurate with the local background. In an alternative algorithm we employ edge detection to identify the pixels that may po tentially belong to microcalcifications. Region growing techniques are then applied and the resulting segmented objects are subjected to tests involving shape, size and gradient. We have examined over 400 mammograms of patients with biopsy proven benign or malig nant abnormalities. Participating radiologists have performed a blind diagnosis of the images without access to other relevant clinical or pathological information. The images were then digitized at 12 bits with 50km or 100am sampling intervals. Preliminary results indicate that both algorithmic and neural network based approaches can segment microcaicifications from the background parenchymal pattern, and that a combination of approaches may be used for improved sensitivity and specificity. Finally, in chapter 9 we have calculated over 100 photometric and morphological features from different microcalcification patterns in 68 digitized mammographic images. These feature vectors were used in computation of a discriminant function that separates the benign and ma lignant classes with overall accuracies better than can be obtained subjectively by experienced radiologists. Chapter 10. Conclusions, and recommendations for future work 10.3 150 Suggestions for future work The problem of automatic detection and recognition of several abnormalities in mammograms, with clinically acceptable sensitivity and specificity, is as yet unsolved. New film digitization devices should be designed and built with larger sensors so that all or most of a mammogram can be digitized with pixels of 25jm . This is necessary so that image 2 detail up to 201p/mm can be captured. New developments in CCD and associated technologies are making it possible to acquire whole breast digital images directly from the patient without the use of radiographic film. This development will open new avenues for characterization of the physical imaging system. Application of image restoration techniques to this new direct-digital modality needs to be investigated. Enhancement of mammographic images to increase the conspicuity of abnormalities is an other area of on-going research. Differentiation of calcifications from dense background needs to be further investigated. Since only about 50% of malignant diseases of breast manifest microcalcifications, the next step in the analysis of mammograms will be the detection of masses. Soft tissue image pro cessing algorithms will have to be developed for segmenting masses from the normal breast background. A number of segmentation algorithms needs to be evaluated for their effectiveness in discriminating poorly defined masses from normal breast parenchyma. Measurements based on texture will be the major criteria here. Various features of the segmented masses will have to be measured in order to classify them as benign or suspicious for malignancy. The third step in the analysis of mammograms will be the detection of secondary signs of cancer. This includes comparison of images of both breasts of the same subject. Although the parenchymal pattern and size of the two breasts may not be identical, radiologists consider a lack of symmetry between the two images as suspicious. Some measure of asymmetry can be computed even though exact registration of the two images is often not possible. Another secondary sign of cancer is the dilation of a single duct. Linear or non-linear Chapter 10. Conclusions, and recommendations for future work 151 edge detection algorithms may be applied to segment the prominent ductal patterns from the background. The ductal size can then be measured and used as evidence of abnormality. Skin thickness can be measured and any skin retraction can be located by changes in curva ture in unexpected places. Any architectural distortions of the breast will have to be identified. A degree of confidence measure may be assigned to any and all measurements. In this way the combined effects of primary and secondary signs of cancer may be weighed for a final classification of the mammogram as normal or suspicious of malignancy. Bibliography 152 Bibliography 153 Bibliography [1] L. W. Bassett and R. H. G. Cold, “Breast Cancer Detection,” Grune & Stratton, 1987. [2] National Cancer Institute of Canada, “Candian Cancer Statistics,” 1992. [3] Cancer Control Agency of B.C., “B.C. Cancer statistics,” Vancouver: 1988-89. [4] M. Moskowitz, “Mammography to screen asymptomatic women for breast cancer,” Amer ican J. Roentgenology, vol.143, pp.457- 459, Sept. 1984. [5] E. S. Paredes, “Atlas of Film-Screen Mammography,” Baltimore: Urban & Schwarzen berg, 1989. [6] J. Tucker and B. Stenkvist, “Whatever happened to cervical cytology automation?,” Anal. Cell. Pathol., Vol 2. pp.259-266, 1990. [7] J. Chen, M. J. Flynn, and M. Rebner, “Regional contrast enhancement and data com pression for digital mammographic images,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 752-758, 1993. [8] P. C. Tahoces, J. Correa, M. Souto, C. Gonzalez, L. Gomez, and J. J. Vidal “Enhancement of chest and breast radiographs by automatic spatial filtering,” IEEE Trans. Med. Image., Vol. 10, No. 3, pp. 330-335. Sept. 1991. [9] S. Yabashi, M. Hata, K. Kubo, and T. Ishikawa, “Image processing for recognition of tumor on mammography,” in Proc. mt. Symp. on Noise and Clutter Rejection in Radars and Imaging Censors, Japan, Nov. 1989. [101 R. Gordon, and R. M. Rangayyan, “Feature enhancement of film mammograms using fixed and adaptive neighborhoods,” Applied Optics, Vol. 23 No. 4, pp. 560-564, 15 Feb. 1984. [11] A. P. Dhawan, C. Buelloni, and R. Gordon “Enhancement of mammographic features by optimal adaptive neighborhood image processing,” IEEE Trans. Med. Image., Vol. 5, No. 1, pp. 8-15, March 1986. [12] A. P. Dhawan and E. L. Royer, “Mammographic feature enhancement by computerized image processing,” Computer Methods and Programs in Bio-medicine, Vol. 27, pp. 23-35, 1988. [13] W. M. Morrow, R. M. Rangayyan, and J. E. L. Desautels “Feature adaptive enhancement and analysis of high-resolution digitized mammograms,” in Proc. Annual mt. Conf. of IEEE Eng. in Med. and Biol., Vol. 12, No. 1, pp. 165-166, 1990. [14] N. Karssemeijer and L. V. Erning, “Iso-precision scaling of digitized mammograms to facilitate image analysis,” in Proc. SPIE, Image Processing, Vol. 1445, pp. 166-177, 1991. Bibliography 154 [15] W. M. Morrow, R. B. Paranjape, R. M. Rangayyan, and J. E. L. Desautels “Region-based contrast enhancement of mammograms,” IEEE Trans. Med. Image., Vol. 11, No. 3, pp. 392-406, Sep. 1992. [16] W. Qian, L. P. Clarke, M. Callergi, H. D. Li, R. A. Clark, and L. M. Silbiger, “Adaptive order statistic filtering and wavelet transform for feature enhancement in mammography,” in Proc. Annual mt. Conf. of IEEE Eng. in Med. and Biol., Vol. 15, No. 1, pp. 62-63, 1993. [17] W. Qian, L. P. Clarke, M. Callergi, H. D. Li, and R. A. Clark, “Tree-structured nonlinear filters in digital mammography,” IEEE Trans. Med. Image., Vol. 13, No. 1, pp. 25-36, March 1994. [18] A. Lame, S. Song, J. Fan, W. Huda, J. Honeyrnan, and B. Steinbach, “Adaptive mutiscale processing for contrast enhancement,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 521-532, 1993. [19] I. E. Magnin, F. Cluzeau, C. L. Odet, and A. Bremond “Mammographic texture analysis: an evaluation of risk for developing breast cancer,” Optical Engineering, Vol. 25 No. 6, pp. 780-784, June 1986. [20] C. B. Caidwell, S. J. Stapleton, and D. W. Holdsworth “Characterization of mammo graphic parenchymal pattern by fractal dimension,” in Proc. SPIE, Vol. 1092, pp. 10-16, 1989. [21] C. B. Caidwell, S. J. Stapleton, D. W. Holdsworth, R. A. Jong, W. J. Weiser, G. Cooke, and M. J. Yaffe “Characterization of mammographic parenchymal pattern by fractal dimension,” Phys. Med. Biol., Vol. 35, No. 2, pp. 235-247, 1990. [22] M. Kallergi, K. Woods, L. P. Clarke, W. Qian, and R. A. Clark, “Image segmentation in digital mammography: comparison of local thresholding and region growing algorithms,” Computerized Medical Imaging and Graphics, Vol. 16, No. 5, pp. 323-331, 1992. [23] S. Hajnal, P. Taylor, M. H. Dilhuydy, B. Barreau, and J. Fox, “Classifying mammo grams by density: rationale and preliminary results,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 478-489, 1993. [24] H. P. Chan, C. J. Vyborny, H. MacMahon, C. E. Metz, K. Doi, and E. A. Sickles, “Evaluation of digital unsharp-mask filtering for the detection of subtle mammographic microcalcifications,” in Proc. SPIE, Medicine XIV / FAGS IV, Vol. 626, pp. 347-348, 1986. [25] H. P. Chan, K. Doi, C. J. Vyborny, C. E. Metz, H. MacMahon, P. M. Jokich, and S. Galliorta, “Digital mammography: development of a computer-aided system for detection of microcalcifications,” in Proc. SPIE, Medical Imaging, Vol. 767, pp. 367-370, 1987. Bibliography 155 [26] H. P. Chan, K. Doi, S. Caihorta, C. J. Vyborny, H. MacMahon, and P. M. Jokich, “Image feature analysis and computer aided diagnosis in digital radiography: automated detection of microcalcifications in mammography,” Medical Physics, Vol. 14, No. 4, pp. 538-548, July/Aug. 1987. [27] H. P. Chan, K. Doi, K. L. Lam, C. J. Vyborny, R. A. Schmidt, and C. E. Metz, “Digital characterization of clinical mammographic microcalcifications: applications in computeraided detection,” in Proc. SPIE, Medical Imaging II, Vol. 914, pp. 591-593, 1988. [28] H. P. Chan, K. Doi, C. J. Vyborny, K. L. Lam, and R. A. Schmidt, “Computer-aided detection of microcalcifications in mammograms, methodology and preliminary clinical studies,” Invst. Radiology, Vol. 23, pp. 671 , 1988. 664 [29] H. P. Chan, K. Doi, C. J. Vyborny, R. A. Schmidt, C. E. Metz, K. L. Lam, T. Ogura, Y. Wu, and H. Macmahon, “Improvement in radiologists’ detection of clustered microcalci fications on mammograms,” Invst. Radiology, Vol. 25, pp. 1102-1110 , 1990. [30] R. M. Nishikawa, M. L. Giger, K. Doi, F. F. Yin, Y. Wu, Y. Jiang, W. Zhang, Z. Huo, U. Bick, C. J. Vyborny, and R. A. Schmidt, “Computer-aided detection and diagnosis of masses and clustered microcalcifications from digital mammograms,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 422-432, 1993. [31] Y. Wu, K. Doi, M. L. Giger, and K. M. Nishikawa, “Computerized detection of clustered microcalcifications in digital mammograms: application of artificial neural networks,” Medical Physics, Vol. 19, No. 3, pp. 555-560, May/June 1992. [32] K. S. Woods, J. L. Solka, C. E. Priebe, C. C. Doss, K. W. Bowyer, and L. P. Clarke, “Comparative evaluation of pattern recognition techniques for detection of microcalcifi cations,” in Proc. of SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 841-52, 1993. [33] B. W. Fam, S. L. Olson, P. F. Winter, and F. J. Scholz, “The detection of calcification clusters in film screen mammograms: a detailed algorithmic approach,” in Proc. SPIE, Med. Imaging II, Vol. 914, pp. 620-634, 1988. [34] B. W. Fam, S. L. Olson, P. F. Winter, and F. J. Scholz “Algorithm for the detection of fine clustered calcifications on film mammograms,” Radiology, Vol. 169, pp. 333-337, 1988. [35] D. H. Davies, D. K. Dance and C. H. Jones, “Automatic detection of clusters of micro calcifications in digital mammograms using local area thresholding techniques,” in Proc. SPIE, Medical Imaging III: Image Processing, Vol. 1092, pp. 153-157, 1989. [36] D. H. Davies, D. K. Dance and C. H. Jones, “Automatic detection of microcalcifications in digital mammograms,” in Proc. SPIE, Medical Imaging IV: Image Processing, Vol. 1233, pp. 185-191, 1990. Bibliography 156 [37] D. Zhao, M. Shridhar, and D. 0. Daut, “Morphology on detection of calcifications in mammograms,” iii Proc. IEEE mt. Conf. Acous. Speech and Signal Processing, part III, pp. 129-132, 1992. [38] D. Zhao, M. Shridhar, and D. C. Daut, “Rule-based morphological feature extraction of microcalcifications in mammograms,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1.905, pp. 702-715, 1993. [39] A. Strauss, A. Sebbar, S. Desarnaud, and P. Mouillard, “Automatic detection and seg mentation of microcalcifications on digitized mammograms,” in Proc. Annual mt. Conf. of IEEE Eng. in Med. and Biol. Vol. 14, No. 5, pp. 1938-1939, 1992. [40] N. Karssemeijer, “A stochastic model for automated detection of calcifications in mam mograms,” in Proc. XIIth Conf. on Information Processing in Medical Imaging, Wye College, Kent, 1991. [41] N. Karssemeijer, “Recognition of clustered microcalcifications using a random field model,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 776-786, 1993. [42] I. N. Bankman, W. A. Christens-Barry, D. W. Kim, I. N. Weinberg, 0. B. Gatewood, and W. R. Brody, “Automated recognition of microcalcification clusters in mammograms,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 731-738, 1993. [43] I. N. Bankman, W. A. Christens-Barry, D. W. Kim, I. N. Weinberg, 0. B. Gatewood, and W. R. Brody, “Algorithm for detection of microcalcification clusters in mammograms,” in Proc. Annual mt. Conf. of IEEE Eng. in Med. and Biol. Vol. 15, No. 1, pp. 52-53, 1993. [44] J. Dengler, S. Behrens, and J. F. Desaga, “Segmentation of microcalcifications in mam mograms,” IEEE Trans. Med. Imag. Vol. 12, No. 4, pp. 634-42, Dec. 1993. [45] H. Barman and G. Granlund, “Hierarchical feature extraction for computer-aided analysis of mammograms,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visual ization, Vol. 1905, pp. 492-503, 1993. [46] W. Qian, L. P. Clarke, M. Callergi, H. D. Li, R. Veithuizen, R. A. Clark, and M. L. Sil biger, “Tree-structured nonlinear filter and wavelet transform for microcalcification seg mentation in mammography,” in Proc. SPIE, Biomedical Image Processing and Biomed ical Visualization, Vol. 1905, pp. 509-520, 1993. [47] W. B. Richardson Jr., “Wavelet packets applied to mammograms,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 504-508, 1993. [48] S. M. Lai, X. Li, and W. F. Bischof “On techniques for detecting circumscribed masses in mammograms,” IEEE Trans. on Med. Imaging, Vol. 8 No. 4, pp. 377-386, Dec. 1989. Bibliography 157 [49] T. K. Lau and W. Bischof “Automatic detection of breast tumors using the asymmetry approach,” Computers and Biomedical Research, Vol. 24, pp. 273-295, 1991. [50] F. F. Yin, M. L. Giger, K. Doi, C. E. Metz, C. J. Vyborny, and R. Schmidt, “Computerized detection of masses in digital mammograms: analysis of bilateral subtraction images,” Medical Physics, Vol. 18, No. 5, PP. 955-963, Sep/Oct 1991. [51] M. L. Giger, F. F. Yin, K. Doi, C. E. Metz, R. Schmidt, and C. J. Vyborny, “Investigation of methods for the computerized detection and analysis of mammographic masses,” in Proc. SPIE, Medical Imaging IV: Image Processing, Vol. 1233, pp. 183-184, 1990. [521 P. Miller and S. Astley, “Detection of breast asymmetry using anatomical features,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 433-442, 1993. [53] D. Brzakovic, X. M. Luo, and P. Brzakovic, “An approach to automated detection of tumors in mammograms,” IEEE Trans. on Medical Imaging, Vol. 9 No. 3, pp. 233-241, Sept. 1990. [54] D. Brzakovic, P. Brzakovic, and M. Neskovic, “An approach to automated screening of mammograms,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualiza tion, Vol. 1905, pp. 690-701, 1993. [55] M. Terauchi and Y. Takeshita, “An approach toward automatic diagnosis of breast cancer from mammography,” in Proc. IEEE Pacific Rim Conf. on Comm. Computers and Sig. Proc., Victoria, BC, Vol. 2, pp. 594-597, May 1993. [56] W. P. Kegelmeyer Jr., “Computer detection of stellate lesions in mammograms,” in Proc. SPIE, Biomedical Image Processing and Three-Dimensional Microscopy, Vol. 1660, pp. 446-453, 1992. [57] W. P. Kegelmeyer Jr., “Evaluation of stellate lesion detection in a standard mammogram data set,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 787-798, 1993. [58] H. D. Li, M. Kallergi, L. P. Clarke, W. Qian, and R. A. Clark, “Markov random field model for mammogram segmentation,” in Proc. Annual mt. Conf. of IEEE Eng. in Med. and Biol. Vol. 15, No. 1, pp. 54-55, 1993. [59] S. H. Fox, U. M. Pujare, W. G. Wee, M. Moskowitz and R. V. P. Hutter, “A computer analysis of mammographic microcalcifications: global approach,” in Proc. IEEE 5th mt. Conf. on Pattern Recognition, pp. 624-631, 1980. [60] I. E. Magnin, M. ElAlaoui, and A. Bremond, “Automatic microcalcifications pattern recognition from x-ray mamrnographies,” in Proc. SPIE, Science and Engineering of Med ical Imaging, Vol. 1137, pp. 170-175, 1989. Bibliography 158 [61] L. Shen, R. M. Rangayyan, and J. E. Leo Desautels, “An automatic detection and classi fication system for caicifications in mammograms,” in Proc. of SPIE, Biomedical Image Processing and Biomedical Visualization, SPIE Vol. 1905, pp. 799-805, 1993. [62] J. Parker, D. R. Dance, D. H. Davies, L. J. Yeoman, M. J. Michell, and S. Humphreys, “Classification of ductal carcinoma in-situ by image analysis of calcifications from mam mograms,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 832-840, 1993. [63] G. L. Cole and M. D. Fox, “Differentiation of mammographic tumor types using the binary joint transform correlator,” in Proc. Annual mt. Conf. of IEEE Eng. in Med. and Biol., Vol. 11, pp. 358-360, 1989. [64] J. Kilday, F. Palmieri, and M. D. Fox, “Linear discriminant based mammographic tumor classification using shape descriptors,” in Proc. Annual mt. Conf. of IEEE Eng. in Med. and Biol., Vol. 12, No. 1, pp. 434-435, 1990. [65] J. Kilday, F. Palmieri, and M. D. Fox, “Classifying mammographic lesions using comput erized image analysis,” in IEEE Trans. Med. Irnag., Vol. 12, No. 4, pp. 664-669, Dec., 1993. [66] K. M. Allain, and G. L. Cote, “Computer classification of stellate carcinomas and fi broadenomas using digitized mammograms,” in Proc. Annual mt. Conf. of IEEE Eng. in Med. and Biol., Vol. 15, No. 1, pp. 56-57, 1993. [67] A. P. Dhawan, Y. S. Chitre, M. Moskowitz, and E. Gruenstein, “Classification of mographic microcalcification and structural features using an artificial neural network,” in Proc. Annual mt. Conf. of IEEE Eng. in Med. and Biol., Vol. 13, No. 3, pp. 1105-1106, 1991. [68] A. P. Dhawan, Y. S. Chitre, and M. Moskowitz, “Artificial neural network based clas sification of mammographic microcalcifications using image structure features,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 820-831, 1993. [69] Y. Chitre, A. P. Dhawan, and M. Moskowitz, “Artificial neural network based classi fication of mammographic microcalcifications using image structure features,” in Proc. Annual mt. Conf. of IEEE Eng. in Med. and Biol. Vol. 15, No. 1, pp. 50-51, 1993. [70] C. J. Burdett, H. G. Longbotham, M. Desa.i, W. B. Richardson, and J. F. Stoll, “Nonlinear indicators of malignancy,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 853-860, 1993. [71] C. Kimme-Smith, S. Dardashti, and L. W. Bassett, “Glandular tissue contrast in CCD digitized mammograms,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 465-470, 1993. Bibliography 159 [72] E. T. Y. Chen, J. S. J. Lee, and A. C. Nelson, “The comparison of mammogram images using different quantization methods,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, PP. 465-470, 1993. [73] S. Astley, I. Hutt, S. Adamson, P. Miller, P. Rose, C. Boggis, C. Taylor, T. Valentine, J. Davies, and J. Armstrong, “Automation in mammography: computer vision and human perception,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 716-730, 1993. [74] M. L. Giger, R. M. Nishikawa, K. Doi, F. F. Yin, C. J. Vyborny, R. A. Schmidt, C. E. Metz, Y. Wu, H. MacMahon, and H. Yoshimura, “Development of a smart workstation for use in mammography,” in Proc. SPIE, Image Processing, Vol. 1445, pp. 101-103, 1991. [75] D. V. Beard, R. E. Johnston, E. Pisano, B. Hemminger, and S. M. Pizer, “A radiology workstation for mammography: preliminary observations, eyetracker studies, and design,” in Proc. SPIE, Medical Imaging V: FAGS Design and Evaluation, Vol. 1446, pp. 289-296, 1991. [76] M. L. Giger, “Future of breast imaging, computer aided diagnosis,” in Syllabus: Categor ical Course in Physics, Technical Aspects of Breast Imaging ed. A.Haus and M.J.Yaffe, Radiological Society of North America, Oak Brook, IL, pp. 257-270, 1992. [77] A. Cues, A. R. Cowen, and C. J. S. Parkin, “A clinical workstation for digital mammog raphy,” in Proc. SPIE, Biomedical Image Processing and Biomedical Visualization, Vol. 1905, pp. 806-817, 1993. [78] L. T. Cook, M. L. Giger, L. H. Wetzel, M. D. Murphy and S. Batnitzky, “Digitized film radiography,” Invest. RadioL, Vol. 24, No. 11, pp. 910-916, 1981. [79] R. M. Hanson, “Digital radiography,” Phys. Med. Biol., Vol. 33 No. 7, pp. 751-784, 1988. [80] P. Bunch and R. Vanmeter, “Noise characterization and reduction in a scanning microdensitometer,” in Proc. SPIE, Medical Imaging, Vol. 767, pp. 236-249, 1987. [81] J. Vranckx and B. Strul, “Performance evaluation of a laser digitizer,” in Proc. SPIE, Medical Imaging, Vol. 767, pp. 524-528, 1987. [82] B. R. Whitting, E. Muka, T. E. Kocher, and M. J. Flynn, “High resolution, high per formance radiographic scanner,” in Proc. SPIE, Medical Imaging IV: Image Formation, Vol. 1231, pp. 295-305, 1990. [83] K. Yip, T. E. Kocher, E. Muka, B. R. Whitting and A. R. Lubinsky, “A performance analysis of medical x-ray film digitizers,” in Proc. SPIE, Medical Imaging IV: Image Formation, Vol. 1231, pp. 508-525, 1990. [84] 5. Batnitzky, S. J. Rosenthal, E. L. Siegel, L. H. Wetzel, M. D. Murphy, C. C. Cox, J. H. McMillan, A. W. Templeton and S. J. Dwyer III, “Teleradiology: an assessment,” Radiology, Vol. 177, pp. 11-17, 1990. Bibliography 160 [85] B. Jaggi, “Design of a quantitative microscope for image cytometry using a solid state detector in the primary image plane,” M.Sc. Thesis, Simon Fraser University, Vancouver, B.C., Canada, 1989. [86] B. Jaggi, M. J. Deen and B. Palcic “Design of a solid state microscope,” Optical Engi neering, Vol. 28 No. 6, pp. 675-682, June 1989. [87] B. Jaggi, S. S. S. Poon, B. Pontifex, J. J. P. Fengler and B. Palcic “Evaluation of a quantitative microscope for image cytometry,” in Proc. Society for Analytical Cytology Meeting, Ashville, N.C., March 1990. [88] B. Palcic, B. Jaggi and C. MacAulay “The importance of image quality for computing texture features in biomedical specimens,” in Proc. SPIE, OR laser, Los Angeles, 1990. [89] R. Holland, S. H. J. Veling, M. Mravunac and J. H. C. L. Hendriks “Histologic multifo cality of Tis, T1-2 breast carcinomas: implications for clinical trials of breast-conserving surgery,” Cancer, Vol. 56, pp. 979-990, September 1985. [90] G. T. Barnes “The use of bar pattern test objects in assessing the resolution of film/screen systems,” in The Physics of Medical Imaging: Recording System Measurements and Tech niques, A.G.Haus ed., Collected papers of American Association of Physicists in Medicine Summer School, North Carolina, pp. 138-151, July 1979 [91] K. R. Castleman “Digital Image Processing,” p. 181, Prentice Hall, N.J. 1979. [92] B. Jaggi, S. S. S. Poon, C. MacAullay and B. Palcic “Imaging systems for morphometric assessment of absorption or fluorescence in stained cells,” Cytometry, Vol.9, pp. 566-572, 1988. [93] B. Jaggi, B. Pontifex, J. Swanson, and S. S. Poon, “Performance evaluation of a 12-bit, 8 Mpel/s digital camera,” in Proc. SPIE, Cameras, Scanners, and Image Acquisition Systems, Symposium on Electronic Imaging Science and Technology, SanJose, Jan 31Feb 4, 1993. [94] A. K. Jam, “Fundamentals of Digital Image Processing,” Prentice-Hall, NY, 1989. [95] R. Gonzales, “Digital Image Processing,” 1987. [96] A. Rosenfeld and D. Kak “Digital Picture Processing,” 1986. [97] P. C. Bunch “Detective quantum efficiency of selected mammographic screen-film combi nations,” in Proc. SPIE, Medical Imaging III: Image Formation, Vol. 1090, 1989. [98] A. C. Haus, Ed., “The Physics of Medical Imaging,” 1979. [99] D. R. Jacobson and C. R. Wilson, “Focal spot size, magnification and sharpness in mam mography,” Med. Phys., Vol. 19, No. 3, pp. 832, May/June 1992. Bibliography 161 [100] M. Rabbani and R. Shaw, “Detective quantum efficiency of imaging systems with ampli fying and scattering mechanisms,” J. of Optical Soc. of America, Vol. 4, pp. 495, May 1987. [101] J. W. Motz and M. Danos, “Quantum noise-limited images in screen-film systems,” in Proc. SPIE, Application of Optical Instrumentation in Medicine X, Vol. 347, pp. 62-66, 1982. [102] B. B.. Hunt and T. M. Cannon, “Nonstationary assumptions for Gaussian models of images,” IEEE Trans. Syst., Man, Cybern., Vol. 6, pp. 876-881, 1976. [103] D. T. Kuan, A. A. Sawchuk, T. C. Strand, and P. Chavel, “Adaptive noise smoothing filter for images with signal dependent noise,” IEEE Trans. Pattern Anal. Machine Intell., Vol. 7, pp. 165-177, 1985. [104] M. Rabbani, “Baysian filtering of Poisson noise using local statistics,” IEEE Trans. Acous. Speech .4 Sig. Proc., Vol. 36, No. 6, June 1988. [105] A. M. Tekalp and M. I. Sezan, “Quantitative analysis of artifacts in linear space-invariant image restoration,” Multidimensional Systems and Signal Processing, Vol. 1, pp. 143-177, 1990. [106] J. W. Woods, J. Biemond, and A. M. Tekalp, “Boundary value problem in image restora tion,” in Proc. sixth mt. Conf. Acous. Speech and Signal Proc. pp. 18.11.1-18.11.4, 1985. [107] B. White and D. Brzakovic, “Two methods of image extension,” Computer Vision, Graph ics and Image Processing, Vol. 50, pp. 342-352, 1990. [108] H. Lim, K. C. Tan and B. T. G. Tan, “Edge errors in inverse and Wiener filter restoration of motion-blurred images and their windowing treatment,” CVGIP: Graphical Models and Image Processing, Vol. 53, No. 2, pp. 186-195, March 1991. [109] H. C. Tan, H. Lim and B. T. G. Tan, “Windowing techniques for image restoration,” CVGIP: Graphical Models and Image Processing, Vol. 53, No. 5, pp. 491-500, Sept. 1991. [110] A. K. Jam, “Fundamentals of Digital Image Processing,” Prentice Hall, NY, pp.283, 1989. [1111 F. J. Harris, “On the use of windows for harmonic analysis with the discrete Fourier transform,” Proc. of the IEEE, Vol. 66, No. 1, Jan 1978. [112] B. R. Hunt and T. M. Cannon, “Nonstationary assumptions for Gaussian models of images,” IEEE Trans. Acous. Speech é.4 Sig. Proc., Vol. 36, No. 6, June 1988. [113] T. Y. Young, K. S. Pu, “Handbook of Pattern Recognition and Image Processing,” Aca demic Press Inc., 1986. [114] R. Vistnes, “Texture models and image measures for texture discrimination,” Computer Vision, Vol. 3, pp. 313-336, 1989. mt. J. Bibliography 162 [115] S. K. Pal and N. R. Pal, “Segmentation based on measures of contrast, homogenity and region size,” IEEE Trans. Sys. Man. Cyber., Vol. SMC17, No. 5, Sept./Oct. 1987. [116] D. Marr and E. C. Hildreth, “Theory of edge detection,” Proc. of the Royal Society of London, Part B, Vol. 207, pp. 187-217, 1980. [117] B. B. Mandeibrot, “The fractal geometry of nature,” Freeman, San Francisco, Ca., 1982. [118] R. M. Haralick, S. R. Sternberg and X. Zhuang, “Image analysis using mathematical morphology,” IEEE Trans. Pattern Analysis and Machine Intelligence, Vol. PAMI9, No. 4, July 1987. [119] W. C. Lin, E. C. K. Tsao, C. T. Chen, and Y. J. Feng, “Neural networks for medical image segmentation,” in Proc. SPIE, Image Processing, Vol. 1905, pp. 376-85, 1991. [120] A. Ghosh, N. R. Pal, S. K. Pal, “Object background classification using Hopfield type neural network,” mt. J. Pattern Recognition and Artificial Intelligence, Vol. 6, No. 5, pp. 989-1008, 1992. [121] M. Lanyi, “Morphologic analysis of microcalcifications,” in “Early Breast Cancer,” Ed. J. Zander and J. Baltyer, Springer-Verlag, 1.985. [122] C. Chen “Statistical Pattern Recognition,” Spartan, 1973. [123] 5. K. Pal and A. Rosenfeld “Image enhancement and thresholding by optimization of fuzzy compactness,” Pattern Recognition Letters, Vol. 7, pp. 77-86, 1988. [124] R. Y. Wong, and E. L. Hall, “Scene matching with invariant moments,” Computer Graph ics and Image Processing, Vol. 8, pp. 16-24, 1978. [125] Kodak, “Image Insight: geometry of image formation,” 1988. [126] W. Webb, Ed., “The Physics of Medical Imaging,” 1988. [127] Kodak, “Image Insight: screen imaging,” 1988. [128] V. L. Newhouse, “Progress in Medical Imaging,” 1988. [129] A. G. Haus, Ed., “The Physics of Medical Imaging,” 1979. [130] R. M. Nishikawa and M. J. Yaffe, “Signal to noise properties of mammographic film-screen systems,” Med. Phys., Vol. 12, No. 1, June 1985. [131] A. G. Haus, “Evaluation of image blur in medical Photography, vol. 61, numbers 1-2, 1985. imaging,” Medical Radiography and [132] F. S. Alcorn, “The protocol and results of training non- radiologist to scan mammograms,” Radiology, Vol. 99, pp. 523-529, June 1971. Bibliography 163 [133] J. N. Wolfe, “Risk for breast cancer development determined by mammographic parenchy mal pattern,” Cancer, Vol. 37, pp. 2486-2492, 1976. [134] R. G. Payster et al., “Mammographic parenchymal patterns and the prevalence of breast cancer,” Radiology, Vol. 125, pp. 387-391, Nov. 1977. [135] R. L. Egan et al., “Mammographic parenchymal patterns and risk of breast cancer,” Radiology, Vol. 133, pp. 65-70, Oct. 1979. [136] P. M. Krook et al., “Mammographic parenchymal patterns as a risk indicator for prevalent and incident cancer,” Cancer, Vol. 41, pp.1093-1097, 1978. [1371 D. F. Rideout et al., “Patterns of breast parenchyma on mammography,” Journal of the Canadian association of radiologist, Vol. 28, Dec. 77. [138] N. F. Boyd et al., “Observer variation in the interpretation of xero-mammograms,” JNCI Vol. 68, no. 3, March 1982. [1391 A. G. Gale et al.,”Computer aids to inammographic diagnosis,” The British J. of Radi ology, Vol. 60, pp. 887-891, 1979. [140] M. Lanyi, “Morphologic analysis of microcalcifications,” in Early Breast Cancer, Ed. J. Zander and J. Baltyer, Springer- Verlog 1985. [141] E. A. Sickles, “Mammographic features of early breast cancer,” American Journal of Roentgenology, Vol. 143, pp. 461-464, Sept. 1984. [142] E. J. Baines, D. V. McFarlane, A. B. Miller and collaborating radiologists, “Sensitivity and specificity of first screen mammography in 15 NBSS centres,” Journal of Canadian Association of Radiologists, V.39, pp. 273-276, 1988. [143] S. Edeiken,”Mammography and palpable cancer of the breast,” Cancer, Vol. 61, pp. 263265, 1988. [144] P. E. Burns, M. G. A. Grace, A. W. Lees, C. May, “False negative mammograms causing delay in breast cancer diagnosis,” J. Can Assoc. Radiol. Vol. 30, pp. 74-76, 1979. [145] B. D. Mann, A. E. Guiliano, L. W. Bassett, M. S. Barber, W. Hallauer, D. L. Morton, “Delayed diagnosis of breast cancer as a result of normal mammograms,” Arch. Surg. Vol. 118, pp. 23-24, 1983. [146] 0. J. Walker, V. Gebski, A. 0. Langlands, “The misuse of mammography in the man agement of breast cancer revisited,” Med. J. Australia Vol. 151, pp. 509-512, 1989. [147] J. B. Buchanan, J. S. Spratt, L. S. Heuser, “Tumor growth, doubling times and the inability of the radiologist to diagnose certain cancers,” Radiol. Clin. North America Vol. 21, pp. 115-126, 1983. Bibliography 164 [148] J. B. Martin, M. Moskowitz, J. It. Milbrath, “Breast cancer missed by mammography,” American J. Roentgenol. Vol. 132, pp.737-739, 1979. [149] L. Kalisher, “Factors influencing false negative rates in xeromammography,” Radiology, Vol. 133, pp. 297-301, 1979. [150] J. N. Wolfe, K. A. Buck, M. Salane, N. J. Parekh, “Xeroradiography of the breast: overview of 21,057 consecutive cases,” Radiology, Vol. 165, pp. 305-311, 1987. [151] S. A. Feig, G. S. Shaber, A. Patchefsky, C. F. Schwartz, J. Edeiken, H. I. Libshitz, It. Ner linger, It. F. Curley, J. D. Wallace, “Analysis of clinically occult and mammographically occult breast tumors,” Am. J. Roentgenol., Vol. 128, pp. 403-408, 1977. [152] L. L. Fajardo, B. J. Hillman, C. Frey, “Correlation between breast parenchymal patterns and mammographers’ certainty of diagnosis,” Invest. Radiol., Vol. 23, pp. 505-508, 1988. [153] It. Holland, J. H. C. Hendriks, M. Mravunac, “Mammographically occult breast cancer: a pathologic and radiologic study,” Cancer, Vol. 52, pp. 1810-1819, 1983. [154] D. Ikeda, I. Andersson, “Ductal carcinoma in situ: atypical mammographic appearances,” Radiology, Vol. 172, pp. 661-666, 1989. [155] P. H. M. Peeters, A. L. M. Verbeek, J. H. C. L. Hendriks, It. Holland, M. Mravunac, G. P. Vooijs, “The occurrence of interval cancers in the Nijmegen screening programme,” Br. J. Cancer, Vol. 59, pp. 929-932, 1989. [156] B. Moskovic, H. D. Sinnett and C. A. Parsons, “The accuracy of mammographic diagnosis in surgically occult breast lesions,” Clinical Radiology, Vol. 41, pp. 344-346, 1990. [157] Council on Scientific Affairs “Council report, Early Detection of Breast Cancer,” Journal of American Medical Association, Vol. 252, No. 21, Dec. 1984. [158] M. L. Giger and K. Doi “Investigation of basic imaging properties in digital radiography,” Medical Physics, Vol. 11, No. 3 pp. 287-295 May/June 1984. [159] D. Goidgof and It. Archarya, Eds, Biomedical Image Processing and Biomedical Visu alization, Proc. of SPIE Biomedical Image Processing Conf., Vol. 1905, SanJose, Ca., 1993. [160] A. C. Gale, S. M. Astly, D. R. Dance, and A. Y. Cairns, Eds. Digital Mammography, Proc. of the 2nd mt. Workshop on Digital Mammography, York, England, Elsevier:NY, 1994. [161] K. W. Bowyer and S Astley, Eds., State of the Art in Digital Mammographic Image Analysis, World Scientiflc:London, 1994. Appendix A Mammography In this Appendix I will describe the process of mammographic image formation and review the factors affecting primary diagnosis of breast cancer from mammographic images. Mammographic Image Formation System A.1 The components of a typical X-ray imaging system are shown in Figure A.1. The photons emit ted by the X-ray tube enter the patient where they may be scattered, absorbed or transmitted without interaction. The primary photons recorded by the image receptor form the image, but the scattered photons create a background signal which reduces contrast. The receptor is a combination of a fluorescent screen and the radiographic film. The screen is added to increase the detective efficiency of the receptor. After chemical development, the film is illuminated and viewed by a radiologist at a distance and magnification appropriate to the detail in the image. The film is viewed as a transparency to provide a greater dynamic range of intensity. The maximum optical density of a film is over 3 while that of a photographic print is only 2 on a logarithmic scale. In a variation of screen-film mammography, a dry non-silver photographic system is used known as Xeroradiography. In a process similar to photocopying, the photoconductive prop erties of amorphous selenium powder is exploited. The resulting image has better spatial resolution and is edge enhanced. This technique has poor broad area contrast and is slower than screen-film receptors thus producing 5 to 10 times more radiation exposure to the patient. Interpretation of a mammogram directly depends on the quality of the mammographic image. The intensity patterns in a mammogram are composed of the image of the breast 165 __ _____ ____ ____ Appendix A. Mammography 166 X-RAY MAMMOGRAPHY SYSTEM • / Focal Spot Beam Limiter • d, Mass/Density Breast / Penumbra / I - 2 d Screen Film A / / /__ 7 Figure A.1: X-ray mammography system Photocell Appendix A. Mammography 167 degraded by the system transfer function and embedded in noise. The breast image is made up of the image of any abnormality superimposed on the image of the background anatomic structures. These abnormalities must be observed, differentiated and classified as benign or malignant. The perception of the sharpness of the signal is influenced by two factors: contrast and extent of blurring. Both of these factors are affected by noise. I will first consider various components along the imaging path and characterize their overall transfer functions, then consider sources of noise and numerically specify a typical screening mammography unit. Effects of digitization are also considered. The numerical measurements given are for the CGR X-ray mammography unit at the B.C. Breast Screening Centre. A.2 Sources of Image Blur In Mammography Image blur is caused by (a) motion; (b) geometric distortion and(c) the screen-film receptor. A.2.1 Motion Blur Motion blur may be due to patient movement or involuntary organ movement. The degree of motion blur is proportional to exposure time. For a given electron beam voltage setting (in Ky), exposure time cannot be reduced below a minimum without adversely affecting the contrast. This minimum exposure time is a characteristic of screen-film combination and is automatically controlled by a photocell measuring the exposure. The use of a breast compression device reduces motion blur to negligible levels. A.2.2 Geometric Distortion Geometric distortion is a function of (a) effective focal spot size; (b) focal spot to object distance and (c) object to film distance [125]. The effective focal spot size is the actual focal spot size foreshortened by the anode angle. The anode angle cannot be smaller than about 100 due to the so-called “heel” effect. The heel effect is the absorption of x-ray photons by the anode itself. Appendix A. Mammography 168 This absorption is strongest for rays that are almost parallel to the anode surface. The heel effect contributes to the non-uniformity of x-ray, weakening the exposure on the anode side of the film (point C on Figure A.2). This non-uniformity is over and above the non-uniformity due to the inverse square law which causes the central ray (point A on Figure A.2) to be stronger than the peripheral rays (points B and C). For a given Ky setting the size of the focal spot is limited by the heat dissipation capacity of the anode. Rotating anodes achieve much smaller focal sizes [126j. The effective focal spot size for screening mammography is 0.3 mm. For minimum distortion the subject should be imaged parallel to the film and perpendicular to the central ray. Higher focal-spot to subject distances provide a more uniform radiation and reduce the penumbra shadow but require increased exposure time or higher Ky settings. High subject-film distances increase magnification but also increase the blur. In dedicated mammographic units the anode to receptor distance is set to 65 cm and the breast is in contact with the receptor. If the average breast thickness is taken as 4.5 cm, the penumbra shadow will thus be less than 23 1 um. The magnification ratio of an object 4.5 cm from the film will be 1.07. A.2.3 The Radiographic Receptor The receptor is a screen-film combination. The screen is a layer of fluorescent material that absorbs almost 90% of x-ray energy and fluoresces with hundreds of times more photons, usually in the blue range. Use of intensifying screens reduces patient dose by an order of magnitude. Double emulsion films are usually sandwiched between two screens. They reduce exposure requirements but contribute to blur. For this reason only single emulsion films are used in mammography. The main source of blur in x-ray mammography, other than the focal spot size, is the screen [127]. Factors affecting receptor blur are: • Distance of phosphor particles in the screen from the film i.e. layer thickness, and any air gap thinner layers produce lower blur but are also less efficient in X-ray photon Appendix A. Mammography 169 THE HEEL EFFECT Anode C A B Note: A = 100% intensity AB consists of a slight increase over 100% intensity and then a general decrease in intensity as B is approached. AC = consists of a considerable decrease in intensity as C is approached. Figure A.2: The heel effect Appendix A. Mammography 170 absorption and conversion. • Phosphor particle size — • Presence of reflector layer larger particles produce more exposure and more blur. — increases photon absorption efficiency of the film, and the effective distance of fluorescent particles to film. • Presence of light absorbing pigments or dyes in the screen • Cross over ratio of unabsorbed x-ray photons — — reduces scatter. increases blur on double emulsion films. • Film emulsion type, density and grain size. High quality screen-film combinations allow: • Lower patient exposure through lower Ky setting and/or lower exposure time. • Faster response and hence less motion blur. • Longer tube life through lower Ky. Factors affecting exposure time are: • Inherent film speed. • Absorption efficiency of the screen i.e. screen thickness, packing density, phosphor type and size. • Spectral match of the film and screen. • Fluorescent efficiency of the screen i.e. ratio of absorbed x-ray to light produced. • Presence of a reflecting layer in the screen. • Presence of light absorbing dyes in the screen. Appendix A. Mammography 171 The above parameters are used by the manufacturers of screen-film systems to obtain ac ceptable images at the lowest possible exposure to the patient. We will need to consider these factors to understand the process of image formation and subsequently develop algorithms for image restoration of mammograms. A.3 A.3.1 The System Characteristic Parameters The Modulation Transfer Function The overall degradation of images can be measured by test objects [125]. Neglecting the effects of noise and poor contrast we can measure or calculate the system point spread function (PSF), the line spread function (LSF), the square wave response function (SWRF), the edge response function (ERF) or the modulation transfer function (MTF). MTF is the magnitude of optical transfer function (OTF). The MTF of three popular screen-film combinations are given in Figure A.3. In the screening mammography program in B.C. the electron beam energy utilized is usually about 27 Ky and the cathode current is set to 100 to 250 mAs range. The fluorescent screen is a single Kodak Min-R medium screen with an Ortho-M single emulsion film. This arrangement gives a spatial resolution of 16 cycles/mm when MTF is 4% of its peak value. The films are processed on an X-omatic processor on extended cycle. A.3.2 Contrast Contrast is a function of both inherent subject contrast and characteristics of the screen-film combination. The subject contrast, 5, is defined as the ratio of the difference in x-ray fluence incident on the receptor between an image point and an arbitrary reference point to the mean value of the two fluences. Since S ranges from 0 to 2, a percent contrast is often quoted, given by 100(S/2). The subject contrast is affected by the scatter-to-primary ratio R = S’/P, the thickness of the target tissue and the difference in linear attenuation coefficients of the normal tissue and any Appendix A. Mammography 172 MTF CURVES FOR THREE KODAK SCREEN-FILM COMBINATIONS 1.O% 4% MIN-R— OrthoM LAX Medium T-MAT G 4’ 0.5 LEX Regular T-MATG 4% o 4% 4% 4 4 S 4% 4% 4% 0 • 4% 0 I I 5 10 Spatial Frequency (cycles/mm) Figure A.3: MTF curves for three Kodak screen-film combinations [131] 15 Appendix A. Mammography 173 abnormality that may be present in the breast. Figure A.4 gives the variation of contrast with photon energy for two objects of importance in mammography. The upper curve is for a 100 calcium hydroxyapatite and the lower curve is for a 1 mm um glandular tissue. The contrast is relative to normal breast tissue. The film characteristics with and without a screen are given in Figure A.5. This is a plot of optical density against exposure. Optical density is defined as 1 —log o T where T is the transmission ratio. A film with a steep characteristic curve has higher contrast but lower latitude. A major contribution to contrast degradation is scatter radiation discussed later. The dynamic range for screening mammography is over 3 optical densities. The lowest densities are limited by base (the natural tint of the film) and fog (background exposure). Let C(u, v) be the Fourier Transform of the x-ray fluence distribution after passing through the subject. If we consider a point object of contrast S then the x-ray fluence is S5(x,y) and C(u, v) = 5, a constant. If we assume the focal spot to be negligibly small and take the Modulation Transfer Function of the film to be 1, then the MTF of the system is due to that of the screen 3 MTF ( u, v). The resulting image contrast is e MTF 10 (u, v). C(u, v) 3 G log . (A.1) . where G is the gradient of the film characteristic curve at optical density of 1 [130]. In general the point slope of the characteristic curve q(Q) should be used instead of G, since G is a constant and does not reflect the non-linearity of the film, where Q is the number of incident . The contrast transfer function, (CTF), is defined as 2 quanta/mm CTF(u, v, q) = MTF(u, v) q(Q) . (A.2) Detectability of an object, of course, depends not only on the contrast but also on the size. The minimum dose required to visualize an object increases as the inverse fourth power of the size of the object [126]. Appendix A. Mammography 174 VARIATION OF CONTRAST WITH PHOTON ENERGY 1— I — 0.1 (P0 5 Ca 0 3 ) 4 H - C 0.01 - Gland .0.001 I. 10 20 30 I. - 40 Energy (keV) Figure A.4: Variation of contrast with photon energy [126] 50 __0 ________________ Appendix A. Mammography 175 SCREEN-FILM CHARACTERISTIC CURVES A. Characteristics of two screen film combinations Kodak Ortho.G Filna/Lanex Screens Base+FogDensity=O.19 Kodak Oitho-G Film/Lanex Screens Base+FogDensity=O.19 t=2.65 Log Relative Exposure B. Comparison of characteristic curves of film (A) with screen-film (B) 3.0 I 0 0.5 I I I I 1.0 1.5 2.0 2.5 Log (relative exposure) Figure A.5: Screen-film characteristic curves Appendix A. Mammography A.3.3 176 Noise Noise in radiography is due to artifacts and mottle. There are three sources of mottle: (1) film graininess; (2) quantum mottle—random distribution of absorbed x-ray quanta; (3) structure mottle—microscopic non-uniformities of screen. The signal to noise ratio (SNR) is defined as zD /W(u, v) SNR (A.3) where W(u,v) is the Wiener noise power spectrum. In general W(u,v,Q) is a function of exposure Q. Substituting from equation A.1 above and taking the subject contrast to be S we get SNR(u, v) = 10 e s MTF G log (u, v) 3 \/W(u, v) . SNR(u,v) . S. -./NEQ (A.4) (A.5) where NEQ, the noise equivalent quanta is a measure of performance. The effect of the receptor on SNR is given by the detective quantum efficiency D SNRI1 ( A6.) The ideal detector—the one counting the photons—will have DQE of 1 [130]. Figure A.6 gives the quantum noise, film noise and total noise for the Min-R screen Ortho-M film com bination [130]. Clearly, at high spatial frequencies the film noise dominates. The noise power due to screen mottle is only about 0.2% of the total noise power. The ability of the film to record quantum mottle decreases to negligible levels beyond 10 cycles/ mm due to the effects of system MTF. Therefore the detection of small objects is limited by the film granularity. A.3.4 The X-ray Scatter During passage through the subject, x-rays are scattered. The thicker or denser portions of the breast produce more scatter [1281. Use of a compression device reduces the scatter. Figure A.7 Appendix A. Mammography 177 NOISE POWER SPECTRUM r T p I I I r T11111 F r 4 io — — —i.. — ‘I 1— •%.., •I%, I. c) 106 ‘I’ 1’ I 0.1 0.2 I I I I I 0.5 .1 til 1.0 F. 2 5 10 20 40 Spatial frequency (cycles/mm) film noise quantum noise total noise Figure A.6: Noise power spectrum; film noise (solid line); quantum noise (dash line); and total noise [130] Appendix A. Mammography 178 gives a plot of scattered to primary radiation ratio as a function of radiation field and thickness of breast. As the field size increases, typically above 8 cm, the ratio increases slowly because x-rays are absorbed within the breast. The receptor usually has higher efficiencies for scattered radiation, worsening the ratio. A grid is often employed to reduce the incidence of scattered radiation on the screen-film. Grids may be stationary and focused with absorption ratios of 3 to 15. The shadow of this grid may be a disturbing artifact on the film. Moving grids are employed in screening mammography units to reduce the effect of this shadow. The grid has 44 lines per centimeter and achieves a grid ratio of 5:1. A scatter degradation factor may be defined as SDF = P/(P + S’) where P is the primary radiation intensity and S’ is the scattered radiation intensity in a given local neighborhood. For a unit with S’/P ratio of 7, typical of chest radiography, we have SDF = 0.125. If now a 12:1 grid is employed, the measured value of SDF improves to 0.5 [1281. The primary beam is also attenuated by the grid. The primary transmission factor Tp for this grid is 0.62 and hence the overall detective quantum efficiency due to scatter and grid is DQE = Tp.SDF = 0.31. Use of a grid requires an increased exposure of about 2 to 8 and gives a contrast improvement of 1.5 to 3.5 [126]. A.3.5 Digitization The analogue x-ray image is sampled and quantized to form a digital image. The digital image is then processed and displayed. Following the work of Giger et al. [158] I assume a square sampling grid and a finite sampling aperture. In general, the sampling interval may not be equal to the aperture. This is the case for 2-D GCD cameras with fill factor less than unity. When they are equal, we can refer to congruent pixels that completely cover the scene. The analogue image is modulated by the sampling aperture to form the pre-sampling image. This image is then operated on by a two-dimensional comb function. The resulting optical transfer function is: Appendix A. Mammography 179 RATIO OF SCATTERED TO PRIMARY RADIATION AS A FUNCTION OF RADIATION FIELD C 1.0 0.8 Phantom 6cm thick - 0.60.4 Phantom 8cm thick - 0.2. 0 2 I- I I 1 4 6 8 10 12 I -I-- 14 16 Diameter of circular radiation field cm Figure A.7: Ratio of scattered to primary radiation as a function of radiation field [128] Appendix A. Mammography OTF(u,v) 180 {[OTFA(u v) OTFs(u, v)]* (u, v, —, —) }• OTFF(u, v) OTFD(U, v) where A, S, F and D refer to analogue image, sampling aperture, linear filter and the display unit respectively and * denotes convolution. If the image preprocessing filter is non-linear the analysis will be more complex. If the aperture for sampling and display are squares of width Ws and WD, then MTF(’u,v) = {[MTFA(u,v).sinc(-irwsu).sinc(lrWsvyj* fl MTFF(u, v) .sinc(7rWDu) sinc(irWjjv) . where 111 In our case Lx = = Ws = 0.1 mm, The display aperture on a 640 x 480 VGA monitor is 0.375 mm x 0.4375 mm. The digital MTF shows a ‘false’ increase over the analogue MTF due to aliasing. To minimize the effects of aliasing, the sampling and display intervals should be decreased. The Nyquist rate for a resolution of 16 cycles/mm (4% MTF of analogue screenfilm) corresponds to pixel sizes of 31 1 um. The digitization process also affects the noise. Assuming square sampling and display pixels of width Ws and WD respectively: Total noise = {[WSA(u,v)sinc2(7rWsu)sinc2(IrWsv)] . * itt (,,_-_)} MT]4(u, v) 2 sinc ( ’JrWDu) 2 sinc ( ’IrWDv) + WSE(U, v) . where Ws is the Wiener noise power spectrum, and E refers to electronic noise. The effect of digitization on signal to noise ratio and threshold contrast can thus be determined. Appendix A. Mammography A.4 181 Mammographic Image Interpretation Reading of mammograms is a highly skilled task requiring long training periods [132]. In screening applications, the films are classified into two groups of “normal” or “suspicious” requiring further investigation. For diagnostic purposes the type and extent of abnormality is also determined leading to recommendation for treatment modality. For biopsy, excision or other surgery, the location of abnormality should also be specified. It has been claimed that the mammographic parenchymal pattern is an indicator of the risk of breast cancer. Wolfe [133] demonstrated a strong relationship between density patterns and the risk of cancer. In current radiological practice, mammograms are classified into one of four “Wolfe grades” and closer observations in shorter time intervals are recommended for the high risk group. These grades in order of increasing risk are: • Ni: Breasts composed primarily of fat • P1: Prominent ducts in the subareolar region involving approximately one third of the breast • P2: Prominent duct pattern involving the major portion of the breast • DY: Considerable amount of collagen or dysplasia with or without identified ducts. Although several researchers [134, 135, 136, 137] have reported results contradictory to Wolfe’s findings, the Wolfe classification scheme is still used. There is a good deal of variability in classifying mammograms. The variations are not only inter-observer but also for the same observer at different time intervals [138]. Individual radiol ogists in general rely on their experience and do not employ quantitative measurements in their line of reasoning for classification. Several studies have attempted to quantify these mammo graphic features and relate them either to risk of cancer or to pathologically proven cases. In an early attempt a computer aided package called MAMMCAD [139] has been developed as an expert system to aid the radiologist. Appendix A. Mammography 182 The common signs of abnormality can be divided into three categories namely, calcifications, masses and secondary signs of cancer. Clusters of microcalcifications are present in over 50% of malignant diseases of the breast. Because of their radio opacity, they constitute the best clue for early detection and over half of clinically occult breast malignancies are discovered on the basis of microcalcifications alone [140]. Microcalcifications are typically smaller than 0.5 mm each and a group of at least 3 within 1 cm 2 constitute a cluster [141]. 75%-80% of clustered calcifications prove to be within benign lesions, but biopsy is performed in most cases in order to detect those 20%-25% that are cancerous. Malignant clustered calcifications have thin linear, curvilinear and branching shapes and round or oval shapes indicate benign lesions. Consistent use of shape for differential diagnosis is difficult mainly due to the small size of microcalcifications. Malignant masses are invasive of the surrounding tissue and consequently are often spiculated and ill defined. The absence of a well-defined edge makes detection difficult especially in dense fibroglandular tissue. It also makes differential diagnosis of benign and malignant masses difficult. Table A.1 gives common features of images of benign and malignant lesions. The most difficult of early cancers to find by mammography are those that contain no calcification and are also surrounded by isodense tissue, impairing delineation of their tumor masses. These can only be detected by secondary signs of cancer which include a single dilated duct, architectural distortion, asymmetry between the right and left breasts, and a developing density as compared to a previous mammogram. A.5 Radiographically Occult Cancers Perception of an object in a radiographic image depends on its size, the amount of illumination, the sharpness of object boundaries, the degree of contrast and the presence of noise. The mammalian visual system has evolved to be responsive to very dwindling fluxes of photons. Since a nerve pulse involves the movement of millions of atoms or ions and since the energy of a single photon is only able to disturb a single atom or molecule, the visual system is a highly Appendix A. Mammography 183 Benign Malignant Relative density Slight to marked increase Character of density Shape Homogeneous Borders Well-circumscribed, regular and smooth; thin layer of surrounding fat Not invaded; displaced trabeculae pushed aside smoothly; no increased vascularity Coarse, isolated, few and countable, not punctate, more apt to have polarity, widely scattered, similar in density may be in penphery of lesion Same size or larger than clinical measurement Always definite increase denser than benign Non-homogeneous, centre densest Tentacled, ragged, spicu lated, variable; spicules heavier toward nipple Poorly circumscribed, irregular, fuzzy or halo Surrounding tissues Calcifications Relative size Round, oval, lobulated Infiltrated; trabeculae retracted irregularly and thickened; increased vasdu larity Numerous, punctate, uncount able, variable density, confined to a measurable area, less polarity, diffuse in lesion, more central Smaller than clinical measurement, often by a factor of 2-4 Table A.1: Radiographic characteristics of breast masses Appendix A. Mammography 184 efficient photomultiplier. The central region of the human retina, the fovea, which at the focal plane of the lens subtends an angle of about 1.5° is lined almost entirely with closely packed cones. Within the fovea, the spacing of cones is sufficiently close (about 10 sum) to enable grating resolution of about 60 cycles per degree. If a film is viewed at a distance of 30 cm, one degree of fovea images about 5 mm of the film. This corresponds to a limiting spatial resolution of 12 cycles per mm. This limit is attainable only for sufficiently high levels of illumination and low levels of image noise. Very low spatial frequencies below 0.2 cycles per degree are also difficult to perceive as can be seen from figure A.8 [126]. The spatial contrast sensitivity of the eye is tuned to sharp edges of about 1 to 3 cycles per degree and is also a function of luminance. In general, there is an inverse relationship between the size of the object and the minimum contrast required for visibility. These limitations of the human visual system contribute to misclassification of mammograms. If it is assumed that a “true” answer exists to the question whether an abnormality is present, for example from subsequent pathological examination of excised tissue, then the radiological decision can be classified into four categories of true positive (TP), false positive (FP), true negative (TN) and false negative (FN). Proven malignancy Radiologic decision yes no yes TP FP no FN TN The following three standard definitions are used in the literature: Sensitivity = Specificity = Accuracy = number of correct positive assessments number of truly positive cases . TP TP + FN number of correct negative assessments TN number of truly negative cases TN + FP number of correct assessments TP + TN total number of cases TP + TN + FN + FP The sensitivity measures the ability of a radiologist to catch the cancers (positives) and specificity measures his/her ability to reject the negative cases. Appendix A. Mammography 185 CONTRAST SENSITIVITY OF THE EYE 1000 1 i s0o 2 SOOcd m .— - •— C,, Q cf_I JU C’, C 0 1OO.0Scdm rJ) 1 01 I_I OS I 1_i S 10 Spatial frequency (cycles/degree) Figure A.8: Contrast sensitivity of the eye [126] 50100 Appendix A. Mammography 186 The Canadian National Breast Screening Study reported an overall sensitivity across 15 screening centers of 75% [142], or 25% false negative mammography. With symptomatic women, false negative rates of as high as 44% among women aged less than 51 years had been re ported [143]. False negative mammography, by providing false reassurance, may cause further delay in diagnoses and treatment. Burns et at. [144] found that 42% of 50 patients with delayed diagnosis due to negative mammograms had metastatic tumor involvement in their axillary lymph nodes at diagnosis. Mann et at. [1451 also observed that 58% of their patients with false negative mammograms have positive axillary nodes at delayed diagnosis. Walker et at. [146] attributed more advanced stage at presentation in women with initially negative mammograms to the delay in their definitive treatment from the first falsely negative mammographic report. Approximately 5% of the false negative mammograms are attributable to poor mammo graphic technique, and an additional 30% are thought to result from observer oversight due to rushed interpretation, heavy caseload, and eye fatigue [147]. Both of these factors are po tentially correctable. However, an unavoidable cause of false negative mammograms is the radiographic density of the breast [148, 149, 150]. A significant correlation between decreas ing diagnostic certainty and increasing complexity of the mammographic breast parenchyma pattern has been shown in two independent studies [151, 152]. Tumor-related phenomena, such as tumor growth pattern and lack of tumor calcifications, may also hamper the visualization of a clinically evident lesion on the mammogram. Roland et at. [153] observed that 5 of the 15 mammographically occult cases were invasive lobular carci nomas, and four of these were situated in dense to very dense breasts. Intraductal carcinomas are also easy to miss especially when they are noncalcified and have only subtle radiological signs [154]. In a Nijmegen breast screening program, twenty-four radiologically occult cancers were located in dense breast parenchyma and approximately half of these were either lobular invasive or ductal non-invasive cancers [155]. Studies on the subject of false negative mammography have generally not contrasted false Appendix A. Mammography 187 negative mammograms with those that were truly positive. In addition mammograms have generally been read for the purposes of review with knowledge that they contained breast cancer. There are biological explanations for the non-visualization of certain breast cancers. Biolog ical characteristics such as extensive radiological density (the predominant component of sheet like, nodular or linear densities), histologic types such as infiltrative lobular or noninfiltrative cancers, as well as gross tumor size less than 2 cm., have been found to associate significantly and independently with false negative mammograms. False positive diagnosis is also undesirable, although it is less damaging. In one study [156] of 69 surgically occult but mammographically apparent cases, 75% of breast biopsies performed proved to be benign. The problem in this case is not detection but recognition (of features) and classification. A high rate of false positive is obviously wasteful of resources, in addition to its undesirable effects on patients. Appendix B Questionnaire For Clinical Evaluation of Images Observer Number Case Number Parameters (Mark with “X”) 1. Number of Microcalcifications (a) < 10/Cm 2 (b) 10 or > 10/Cm 2 2. Shape of Microcalcifications (a) semi-lunar (tea-cup) (b) linear (c) round (d) Oval / rod shaped (e) curvilinear (f) branching 3. Density of Microcalcifications (a) uniformly dense (b) eggshell (c) poorly defined (smudgy) 188 Appendix B. Questionnaire For Clinical Evaluation of Images 4. Margination of Microcalcifications (a) smooth borders (b) less sharp to irregular 5. Spatial Arrangement (a) scattered (may have polarity) (b) clustered or confined to one area 6. Relation to Mass (a) no associated mass (b) concentrated in core of mass (c) concentrated in periphery of mass (d) scattered evenly throughout mass 7. Overall Clinical Assessment of Microcalcifications (a) benign_ (b) malignant 189 Appendix C Features C. 1 Features from Individual Microcalcifications: • Geometrical: Microcalcification identification (id); coordinates of the enclosing box (xl, yl, x2, y2); coordinates of the center of mass (xm, ym); intensity of the background. • Photometric: Gray level intensity of microcalcifications (mean-intensity, var-intensity); contrast (depth). • Statistics from the Optical Density histogram: mean (odmean); minimum (odmin); max imum (odmax); range (odrange); variance (odvar); skewness (odskew); kurtosis (odkurt). • Morphological: area; perimeter; mean-radius; var-radius; compactness; inertia; elonga tion; coarse and fine measures of boundary roughness (bdycrc, bdycrc2); distance-tocluster-center; aspect ratio; angle-of-major-axis; alignmentl; alignmentll; alignment-withcluster-axis; C.2 Features from Microcalcification Clusters: Photometric Features mm-contrast, max-contrast, range-contrast, mean-contrast, var-contrast Morphological Features: • total area of microcalcifications (total-area); area of the smallest microcalcification (mm area); area of the largest microcalcification (max-area); range, average, variance and 190 Appendix C. Features 191 normalized standard deviation of area (range-area, mean-area, var-area, nrstdv-area); • Simple statistics from the perimeter (mm- perimeter, max-perimeter, range-perimeter, mean-perimeter, var-perimeter, nrstdv-perimeter) and the average radii (mm-mean-radius, max-mean-radius, range-mean-radius, mean-mean-radius, varmean-radius) of each calci ficat ion. • Descriptions of variance of radii of each microcalcification (mm-var-radius, max-varradius, range-var-radius, mean-var-radius, var-var-radius). • Statistics of compactness (mm-compact, max-compact, range-compact, mean-compact, var-compact); and inertia (mm-inertia, max-inertia, range-inertia, mean-inertia, var-inertia). • Coordinates of the cluster center (cluster-center-x, cluster-center-y). • Distance of each microcalcification from the cluster center: mm-distance-to-center, maxdistance-to-center, range-distance-to-center, mean-distance-to-center, var-distance-to-center • Distance of each microcalcification from the cluster center weighted by the area of the mi crocalcification: mm- area-x- distance, max-area-x-dist ance, range-area-x-distance, mean area-x-distance, var-area-x-distance. • Distance of each micro calcification from the cluster center weighted by the circumference of the microcalcification: min-perimeter-x- distance, max-perimeter-x-dist ance, range-perimeter x- distance, mean-perimeter-x-dist ance, var-perimeter-x- distance • Distance between each pair of microcalcifications: mm-pair-distance, max-pair-distance, range-pair-distance, mean-pair- distan Ce, var-pair-distance. • Descriptions of the convex polygon enclosing the cluster: polygon-area, polygon-perimeter, polygon-compactness, polygon-aspect-ratio, mm-polygon-radii, max- polygon radii, range- polygon-radii, mean-polygon-radii, var-polygon-radii Appendix C. Features 192 • Ratio of the area of the convex polygon filled by microcalcifications: polygon-caic-area ratio. • Alignment measures: alignment-field-intensity, magnetic-potential mm- alignmentl, max alignmentl, range- alignmentl, mean- alignmentl, var- alignmentl mm- alignmentll, max alignmentll, range- alignmentll, mean-alignmentll, var-alignmentll, mean-alignment-withcluster-axis. • Moments: moment2O, momentO2, momentli, moment3O, momentO3, moment2l, mo mentl2. • Normalized shape invariant moments: invariant- shape 1, invariant-shape2, invariant-shape3, invariant-shape4, invariant-shape5, invariant-shape6.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Digitization and analysis of mammographic images for...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Digitization and analysis of mammographic images for early detection of breast cancer Aghdasi, Farzin 1994
pdf
Page Metadata
Item Metadata
Title | Digitization and analysis of mammographic images for early detection of breast cancer |
Creator |
Aghdasi, Farzin |
Date Issued | 1994 |
Description | X-ray mammography is the proven method for early detection of breast cancer. Digital processing and analysis of mammographic images can potentially assist in improved performance of radiologists in earlier detection and recognition of abnormalities. In this work a novel image acquisition system based on an area scanning CCD array has been developed for the digitization of mammograms at high spatial and photometric resolutions. The system characteristic parameters were measured. The quality of the resulting images in terms of sharpness and noise content is comparable with that obtained by the more expensive and slower drum laser-scanning microdensitometer. The clinical application of soft-copy display of digitized images are evaluated. To further improve the quality of the images, restoration algorithms were applied to restore the images from the degrading effects of the system’s blur and noise. Performance of three filtering techniques was compared. A new method for the reduction of boundary truncation artifacts in image restoration was suggested and studied. The process of radiographic image formation was modeled and two locally adaptive smoothing filters were employed to counter signal-dependent radiographic noise before application of restoration filters. The results of the restored images show a marked improvement in detectability of smallest particles of microcalcifications when judged by a human observer. Image segmentation routines were developed to separate microcalcifications from the background parenchymal pattern. Performances of two algorithmic approaches to segmentation and two artificial neural networks were compared. Over 100 numerical features were automatically extracted from the clusters of microcalcifications. These features were evaluated for their ability to separate the benign and malignant formations. Using a database of 68 digitized mammograms a discriminant function was calculated. The sensitivity and specificity of this approach in recognition of malignant microcalcification clusters is shown to be comparable to that of trained radiologists. |
Extent | 4485566 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-06-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0065327 |
URI | http://hdl.handle.net/2429/8738 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1995-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1995-982454.pdf [ 4.28MB ]
- Metadata
- JSON: 831-1.0065327.json
- JSON-LD: 831-1.0065327-ld.json
- RDF/XML (Pretty): 831-1.0065327-rdf.xml
- RDF/JSON: 831-1.0065327-rdf.json
- Turtle: 831-1.0065327-turtle.txt
- N-Triples: 831-1.0065327-rdf-ntriples.txt
- Original Record: 831-1.0065327-source.json
- Full Text
- 831-1.0065327-fulltext.txt
- Citation
- 831-1.0065327.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065327/manifest