Automatic segmentation of the thoracic organs for image registration and RT planning
author: Alexandra Koulouri supervisors: Prof. D. Hawkes and Dr. J. McClelland
This work is submitted as part requirement for the MSc degree in Medical Image Computing at University College London. It is substantially the result of my own work except where explicitly indicated in the text.
CMIC Department of Medical Physics and Bioengineering University College London University of London 2009
Abstract We present a novel method for the automatic segmentation of the thoracic cavity and the detection of human lungs and the major thoracic organs, as a necessary pre-processing step for a subsequent deformable registration scheme. Our method is divided into two parts. In the first part, a coarse separation of the body into two subregions, one with the skeleton, fat and skin and another with the lungs, the heart and the rest soft tissues is achieved by approximating the rib cage with a c-spline curve. The second part employs curve fitting methods in order to detect the boundaries of the lungs and the heart with better accuracy, based on the estimation of the first part. In particular, the proposed method uses an automatic threshold technique for the rib cage definition, followed by a cubic interpolation scheme for a rough approximation of the cross-sectional inner thoracic area in CT transversal slice. A robust coupled phase Level Set framework is then applied to efficiently separate the pulmonary cavities and the heart from the rest biological structures. Level Sets initialization embeds information from the first part. The algorithms and the mathematical tools developed during the project prove feasibility for an accurate and reliable method for the segmentation of thoracic cavity from CT data, that in the future could be used to benefit patients with lung or thoracic diseases and more precisely as a first step in oncological applications and radiotherapy treatment planning.
Acknowledgements Firstly, I would like to thank Professor D. Hawkes and Dr. J. McClelland, for their invaluable guidance and assistance throughout my year at UCL. Moreover, I would like to thank the Latsis Public Benefit Foundation and the EPRC for their financial support. What is more, I would not get to the end of my project work without the invaluable advice and support of friends within and outside the College.
i
Preface With the advances in computer technology over the past decades, the amount of research in the field of medical image has grown dramatically. Computed Tomography(CT), magnetic imaging (MRI) and other image modalities provide an effective non-invasive way of mapping many anatomic structures and allow doctors and physicians to virtually interact with body structures and learn life saving information. Nowadays, with the well established field of medical image processing, we are able to go from a simple visualization of anatomical structures to patient diagnosis, advanced surgical planning and human organ simulations. Although current volume visualization techniques can potentially provide accurate and high quality 3D view of anatomical structures, their utilization for reliable and efficient analysis is still under research. These limitations derive from the complexity and the variability of the anatomic organs. The problem is even more difficult if we consider the shortcomings of imaging modality, such as sampling artifacts, noise, low contrast etc. which may cause the boundaries of anatomical structures to be indistinct and disconnected. In order to tackle these difficulties there are techniques which can efficiently delineate and separate out biological structures e.g. heart, kidney, or region of interest and are known as Image Segmentation Techniques. Image Segmentation Techniques benefit either from the structural characteristics or the statistical and stochastic features of a human organ morphology. Today there is a wide variety of segmentation approaches depending on the specific application, the imaging modality (CT, MRI, etc.) and other factors. There is currently no segmentation algorithm that provides acceptii
able results for every medical data-set and exactly this variability is what makes segmentation a very challenging problem. The present dissertation is another challenge in the field of medical image processing, focusing on the segmentation of the thoracic cavity from CT images. This dissertation is a part of a more general effort of the CMIC group to contribute with state-of-the-art methods to an efficient lung radiotherapy planning.
iii
Contents 1 Introduction 1.1 Project Aim and Objectives . . . . . . . . . . . . . . . . . . . 1.2 Lung Cancer . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Brief Disease Description . . . . . . . . . . . . . . . . . 1.2.2 Image Acquisition . . 1.2.3 Current Treatments . . 1.3 Data Acquisition Information 1.4 Motivation . . . . . . . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
2 3 4 4
. . . .
4 4 6 7
1.5 Text Structure . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2 Previously Applied Methods 10 2.1 Region Growing Methods . . . . . . . . . . . . . . . . . . . . . 11 2.2 2.3 2.4 2.5
Active Shape Models Deformable Models . Hybrid Methods . . . Summary . . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
11 12 13 13
3 Mathematical Background 3.1 Curve Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Least Squares Method . . . . . . . . . . . . . . . . . . 3.1.2 Cubic Splines . . . . . . . . . . . . . . . . . . . . . . .
14 14 14 16
3.1.3 Splines vs Polynomials . 3.2 Optical Flow . . . . . . . . . . 3.2.1 Mathematical Approach 3.2.2 Initial Value Problem . . iv
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
19 20 21 22
3.2.3
Numerical Calculation . . . . . . . . . . . . . . . . . . 23
3.3 Euler Lagrange Equation . . . . . . . . . . . . . . 3.4 Naive Bayes Classifier . . . . . . . . . . . . . . . 3.5 Level Sets Methods . . . . . . . . . . . . . . . . . 3.5.1 The Chan-Vese Level Set Method . . . . . 3.5.2 Front propagation without Re-initialization:
. . . . . . . . . . . . . . . . . . . . . . . . Level Sets
. . . .
24 26 27 29
Approach . . . . . . . . . . . . . . . . . . . . . . . . . 35 4 Proposed Method 38 4.1 Pre-Processing . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.1.1 Median Filter . . . . . 4.1.2 Contrast enhancement 4.2 Rib Cage Approximation . . . 4.3 Heart Position . . . . . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
40 40 41 44
4.4 Body Volume Definition . . . . . . . . . . 4.5 Thoracic Cavity Extraction . . . . . . . . 4.5.1 Levels Sets Functions Initialization 4.5.2 2-Phase Level Sets Implementation
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
45 47 48 49
5 Experimental Results 55 5.1 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.1.1 Results and Discussion . . . . . . . . . . . . . . . . . . 61 6 Conclusions and Future Work
v
70
List of Figures 1.1 Cross-section CT image
. . . . . . . . . . . . . . . . . . . . .
5
1.2 Boundary Condition in Lungs Cavity . . . . . . . . . . . . . . 1.3 Sliding Motion . . . . . . . . . . . . . . . . . . . . . . . . . .
8 8
3.1 Splines vs Polynomials . . . . . . . . . . . . . . . . . . . . . . 20 3.2 Signed distance function φ(~x) defining the regions Ω+ and Ω− as well as the boundary ∂Ω. . . . . . . . . . . . . . . . . . . . 29 3.3 All possible cases in the initial position of the curve A . . . . . 31 4.1 Flow Chart . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.2 re-scaling graph . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.3 (A)P (X ≥ T hresh) and (B) Segmented Image . . . . . . . . . 42 4.4 (A) Estimation of minimum distance in each angular direction and (B) pixels with the minimum distance . . . . . . . . . . . 43 4.5 (A)Result of curve fitting and (B) coarse separation of thoracic region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.6 4.7 4.8 4.9
Filling process, 4-connected region . . . . . . . . . (A) Mask: result of filling process and (B) the CT (A)CT slice, (B) outer part and (C) inner part. . Graph of initial zero level set functions . . . . . .
. . . slice . . . . . .
. . . .
. . . .
. . . .
. . . .
46 46 47 49
4.10 (A)φ2 level set function and (B) φ1 level set function . . . . . 49 4.11 evolution of φ2 (A) without rib cage stopping criterion and (B) with stopping criterion. . . . . . . . . . . . . . . . . . . . 51 4.12 evolution of φ2 (A) without (1 − h(φ1 )h(φ2 )) term and (B) with (1 − h(φ1 )h(φ2 )) term. . . . . . . . . . . . . . . . . . . . 51
vi
4.13 Initial zero level sets, (B) final zero level sets (500 iterations) and (C) segmented inner body part . . . . . . . . . . . . . . . 54 5.1 (A) Initial Image and (B) Enhanced Image, Median filtering 255 (P SN R = 20 log10 ( RM ) where RMS is the mean least square S error between the initial and the filtered image). . . . . . . . . 56 5.2 (A) Thresholded Image and (B) Coarse Separation. . . . . . . 56 5.3 Displacement estimation of rigid body parts . . . . . . . . . . 57 5.4 (A) Heart position according to Classification and (B)Classification Map: magenta is the heart regions, light violet the rest gray value region and in light blue the background . . . . . . . . . 58 5.5 Body mask . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.6 Top left: Initial zero levels, top right: Zero level (iter=600) , bottom left: Final contour, bottom right: Segmented regions . 60 5.7 1st patient, Final rib cage curves (A) slice 84, (B) slice 102, (C) slice 161 and (D) 174. . . . . . . . . . . . . . . . . . . . . 62 5.8 1st patient, Colomns:(A)Initial contours, (B)Final contours. . . 63 5.9 1st patient, Colomns:(A) Lung segmentation using proposed method (green) vs Region Growing (Red) and (B) Heart segmentation, automatic(blue)-manual(magenta). . . . . . . . . . 64 5.10 2nd patient: specificity, sensitivity and DICE for the segments segment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.11 2nd patient: specificity, sensitivity and DICE for the lungs segment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.12 2nd patient, slice 80-no satisfactory rib cage approximation. . . 67 5.13 3rd patient, Colomns:(A)Initial contours, (B)Final contours and (C) Final rib cage curve. . . . . . . . . . . . . . . . . . . 68
vii
List of Tables 5.1 Results of Heart Segmentation for the 1st patient . . . . . . . 65 5.2 Results of Lungs Segmentation for the 1nd patient . . . . . . . 65 5.3 Mean value and standard dev. of Specificity, Sensitivity and DICE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.4 Results of Heart Segmentation for the 3rd patient . . . . . . . 69 5.5 Results of Lungs Segmentation for the 3rd patient . . . . . . . 69
1
Chapter 1 Introduction Accurate and automatic segmentation of medical image is a field of great interest, as it allows the fast extraction of image features and other regions of interest for further medical and clinical examination. Furthermore, image segmentation is indispensable for proper 3D visualization of consecutive stacked images such as Computer Tomographies or MR Images. Recently, 3D reconstruction and visualization of biological structures has been extensively used in medical applications. The virtual inspection of human organs is very helpful for advanced surgical planning or patient diagnosis. It is even possible to simulate the functionality of an organ or a tissue as it is in reality. The present project intends to design a useful segmentation tool for further research in the field of lung radiotherapy, improving treatment planning or providing essential information to the clinical doctors. More precisely, the proposed method focuses on the automatic segmentation of the whole thoracic cavity (lungs, heart) from CT images. The segmentation scheme, followed by the registration process, enables the estimation of the motion of the organs during the breath cycle and thus provides useful information to the radiotherapist about the position of the organs and the specification of the dose distribution that should be used between the tumor and the organs at risk.
2
1.1
Project Aim and Objectives
The main goal of this thesis is to design an efficient, accurate and automatic method for the segmentation of the whole thoracic cavity from CT images. The method is divided into two parts. In the first part a “rough” separation of the cross sectional slice into two subregions is performed, while in the second an accurate extraction of the pulmonary cavity and the major organs (i.e. Heart) is achieved using knowledge of the first part. In particular, the coarse separation of the body into two subregions, one with the skeleton, fat and skin and another with the lungs, the heart and the rest soft tissues is achieved approximating the rib cage with a c-spline curve. This curve roughly segments the body and is the starting point for the segmentation of the thoracic cavity. For the refined segmentation of the cavity, a coupled level set scheme is used for the simultaneous extraction of the heart and the lungs. For this purpose, the level set framework employs two curves. The curves’ initialization is based on the definition of an initial position for the heart and the lungs. Naive Bayes classifier is applied for the initialization of the curve for the heart detection and the rib cage curve of the first part is used for the lungs initialization. The method was applied in 3 computer tomographies of approximately 250 images each, giving promising results for the automatic detection of inner thoracic organs and the extraction of the thoracic cavity above the diaphragm. Upper abdomen area which includes diaphragm and liver remain a difficult task due to the lack of ribs and the intensity similarities between the organs and the body fat. Hence, only lungs can be segmented efficiently in this area. Although the knowledge of the position of the organs and the rib cage in the upper thoracic region could help to formulate a future scheme for estimations below the diaphragm.
3
1.2
Lung Cancer
1.2.1
Brief Disease Description
Lung cancer is the malignant lesions in the lung tissues, usually in the cells lining air passages. There are two main types of lung cancer, the small cell cancer and non-small cell lung cancer. Lung cancer has one of the lowest survival outcomes of any cancer because over two-thirds of patients are diagnosed at a late stage when curative treatment is not possible. Earlier diagnosis and referral to specialist teams would make a significant difference to survival rates. According to latest surveys in England and Wales (source“National Cancer Institute”) around the 25% of all lung cancer patients are alive one year after diagnosis falling to 7% at five years.
1.2.2
Image Acquisition
There are 2 main tests to diagnose and to examine the lung cancer progress, the spiral CT scans (helical CT scans) and the Bronchoscopy [35]. • the spiral CT scans (helical CT scans) rotate rapidly around the body, taking more than one hundred pictures in sequence. These scans can detect smaller lung tumors than a conventional CT scan and the examination takes only a few minutes. • With bronchoscopy, a careful examination of the inside of the lung airways is made, taking samples (biopsies) of the cells. Normally a thin, flexible tube called a bronchoscope is used and the test is carried out under local anaesthetic.
1.2.3
Current Treatments
There are different ways and medical options of an efficient and appropriate treatment planning. Surgery, radiotherapy and chemotherapy may be used separately or together to treat people with lung cancer [35].
4
Figure 1.1: Cross-section CT image A multidisciplinary team who specialize in treating lung cancer and in giving information and support including surgeons, who are experienced in lung surgery, oncologists, radiologists, pathologists counselors and psychologists have to take under consideration a number of factors before any treatment decision is made. In particular some of the basic factors are • the patient general health • the type of lung cancer (small cell or non-small cell) • the size and position of the tumour • whether it has spread beyond the lung (the stage of the cancer). Today the most common treatment methods considering the above factors are 1. Surgery 2. Radiotherapy 3. CHART Radiotherapy 4. Chemotherapy 5. Radio-frequency ablation 6. Cancer growth inhibitors
5
7. Controlling symptoms 8. Cryosurgery. As the present project is a part of a radiotherapy planning treatment, we focus only on radiotherapy method, giving some further detail in the following paragraph. Radiotherapy Radiotherapy treats cancer by using high-energy x-rays to destroy the cancer cells, while doing as little harm as possible to healthy cells. The treatment is given in the hospital radiotherapy department. The number of treatments, and the length of time they take, will depend on the type of cancer as well as its size and position. To make sure that the radiotherapy works as well as possible, it has to be carefully planned. Before the external treatment, the therapy radiographers takes pictures with a special CT scanner. CT imaging is extremely important for the radiotherapy planning as it visualizes the tumor motion throughout the breath cycle and enables the collection of measurements about the position, the size of the tumor and the surrounding healthy tissues with manual or semi-automatic processes. All these measurements are used to define the x-ray beam direction and the dose distribution. Frequently marks may be drawn on the skin to help the radiographer, to position the patient accurately and to show where the rays have to be directed.
1.3
Data Acquisition Information
In the current project a Siemens Somaton Sensation Multi-slice CT scanner standard diagnostic scan with voxel dimensions of 0.97 x 0.97 x 1 mm which covers whole lungs was used for the data acquisition. The data was acquired using the Inhalation and Exhalation breath-hold CT scheme [11], thus two volumes are obtained. The main disadvantages of this process are the difficulty of the patient to hold his breath reproducibly 6
and that muscles and tissues can be different for breath-hold and free breathing cases, giving inaccurate information. Today, a 4D-CT data acquisition scheme is applied which tries to capture thoracic motion during free breathing with main drawback the increased patient irradiation. Although, as a first approach, the Inhale-Exhale data is useful for the validation of the proposed segmentation method.
1.4
Motivation
For a fast and efficient radiotherapy treatment planning of lung cancer, a successful registration procedure, which aligns the CT data and tracks accurately the thoracic volume during breathing cycle, is extremely essential for the radiologists. Breathing motion changes dramatically the positions of the thoracic biological structures. During radiotherapy the position of the soft tissues and the tumor have to be defined as the x-rays beam direction has to be accurate. Registration procedure intends to capture the tissues deformations and displacements during breath cycle and to give information about the magnitude and the direction of these changes which will be used and evaluated from the clinician during planning process. Thus, the application of an efficient registration algorithm is extremely important in radiotherapy planning. Currently deformable registration formulations [9], [38] assume a smooth and continuous deformation vector field (displacement field). By definition, these algorithms have poor performance in the cases where there are discontinuities in motion field, giving inaccurate results. Sliding motion of the lung against chest wall (figure 1.2) during breathing creates a discontinuous vector field which is a challenging problem for image registration procedures. Actually, the existence of pleura which is a transparent membrane which covers the lungs and lines the inside of the chest walls, provokes this sliding and subsequently the displacement vectors on the two sides of pleura are discontinuous. Thus, in many cases lung deformable registration algorithms fail to correctly register the CT volume, as it is referred by McClelland et al. [10]. 7
As a solution to this problem Wu et al. [38],[10] proposed to separate the Pleura Fat/Skin Link Lung/Pleura
Link pleura/Fat tissue
Lungs
Link pleura/rib
Figure 1.2: Boundary Condition in Lungs Cavity
Figure 1.3: Sliding Motion thoracic region into two subregions and to apply the registration algorithm in each subregion separately. The manual segmentation was made concerning moving and less moving subregions. Hence, the subregion with the highest motion properties includes the inner thoracic organs (i.e. the lungs, the heart and the liver) and the less moving subregion consists of the skeleton (ribs, spine) and the fat-skin. Registration in separately subregions [38] revealed promising results, capturing the discontinuities of the deformation vector field. Till now there is no automatic segmentation tool for the subregions separation, making the process extremely time consuming. In the present project, we intend to develop the algorithms and the mathematical tools for the segmentation of these two subregions from CT data. The segmentation result 8
can be used for the initialization of a subsequent non-linear registration procedure.
1.5
Text Structure
The rest of the text consists of four chapters. The 2st chapter is introductory and defines a brief description of the previously related works which mostly concerns lung segmentation techniques. The 3rd chapter gives an overview of the mathematical background which is used in the proposed method. C-spline theory for the rib cage approximation as well as the mathematical definition of level sets which is used for the lungs and heart segmentation are explained. In the 4th chapter a step by step description of the proposed method is reported. Finally in the 5th chapter experimental results are presented while in chapter six future work is discussed.
9
Chapter 2 Previously Applied Methods A plethora of different approaches have been proposed in the area of thoracic and abdominal area segmentation in general [18],[24]. However, relatively few methods deal with the segmentation of the thoracic cavity as a preprocessing step for the initialization of the registration procedure. Most of them propose a semi-automatic scheme and usually focus only on the lungs. The automatic segmentation of both the lungs and the soft tissue (heart etc) that lie inside the thoracic cavity is a big challenge because organs and lungs have different image intensity properties while the fat tissue surrounding the rib cage has similar luminance with the heart or other tissues inside the rib cage. Considering that the heart and the fat tissue have noticeable contact regions (like the front side of the thorax) their separation is extremely difficult. Today, the most commonly used methods of lung segmentation are based on region growing techniques [33],[8], which rely on the large gray value contrasts between the lungs and the surrounding tissues. However, these methods usually fail on detecting dense pathologies and most of the time manual intervention is necessary.
10
2.1
Region Growing Methods
Region growing is one of the simplest region-based image segmentation algorithms. Region-growing exploits the fundamental fact that pixels which are close together have similar gray values. Starting from single points (“seed point”), examining the neighboring pixels of the “seed points” and adding them if they are similar to the seeds, this technique defines in the end a smooth connected region. Hu et al. [33] proposed a method which employs the conventional region growing algorithm for lung segmentation in CT data. In many cases this approach was used as the initial step of later methods, as for example the recent work of Rikxoot [8]. In particular the Hu method consists of three major steps: the extraction of the lung region from the CT images by region growing using automatic threshold selection, the separation of the left and right lung, and the final smoothing of the lung boundaries. Automatic threshold selection is based on the Otsu’s method [29]. After the threshold selection, the seed point for the region growing process is chosen by finding low intensity value pixels inside the lungs. For this purpose, the image background is extracted and thus only the body volume with the low intensity regions (lungs) and the higher intensity regions (tissues) are kept. The Hu method has been mainly tested in healthy CT volumes giving quite promising results. However, the density of the lungs is influenced by a variety of factors such as tissue volume, air volume, image acquisition protocol and physical material properties of the lung parenchyma. These factors make the automatic selection of a gray-level threshold difficult and in many times this method fails to segment the lungs accurately.
2.2
Active Shape Models
In the cases where gray value information is not sufficient to distinguish lungs from no-lungs area, active shape modeling could be useful. However, neither the 3D shape model design of a lung, nor the acquisition of a large amount of the required training data are trivial tasks. 11
In [4], a coarse segmentation of the thoracic main biological structure is presented, based on the use of anatomical models and hierarchical model matching techniques. Lelieveldt et al. [4] proposed the construction of complex anatomical models which describe the biological modalities of the thoracic volume, applying volume-density and surface-based approaches. These models are used for the “rough” extraction of the main thoracic organs. The approach of Lelieveldt et al. can be the first step for defining multi-modality training data sets. However, till now, there is not any approach which combines the information from all the thoracic structure.
2.3
Deformable Models
Deformable models (Active Contours) attempt to keep desirable shape characteristics, e.g. smoothness, and simultaneously preserve strong edges in the image data. Deformable model methods can be divided into two categories: parametric deformable model methods, i.e. snakes, and geometric deformable model methods, i.e. level sets. Parametric deformable model methods define object contours by using parametric curves that are deformed under internal and external forces. Level set methods extend the segmentation problem to one dimension higher and represent the curve/surface as the zero level set of this higher-dimensional function, which evolves according to a given form of a speed function. Active contours for lung segmentation are actually used as a second step in most of previous methods [24],[21] in order to refine the segmented region. This approach is necessary since it is important to have prior knowledge of the location and other properties of the object we want to segment. Actually, deformable models have not been used extensively in lung segmentation due to their mathematical complexity and their time-consuming implementation. However, their flexibility and adaptation in different problems enable them to solve difficult tasks.
12
2.4
Hybrid Methods
There are many reported methods for the lung segmentation. The most interesting of them is the segmentation by registration scheme [14]. In this method, a scan with normal lungs is elastically registered to a scan containing pathology. After that the resulting transformation is applied to a mask of normal lungs, and a segmentation of the pathologic lung is found. For the mask definition, a set of healthy lungs is used. According to [14] this scheme gives promising results for the definition of pathologic lungs. However as it is referred, this method is extremely time consuming and with high computational cost. Other methods which , however, do not provide any rigorous validation are presented in [2], [1] and [25].
2.5
Summary
Studies of previous work indicate that segmentation of thoracic cavity as seen in CT images is a problem that has received relatively little attention in the past. Much of the work devoted to create automatic segmentation schemes of lung from CT images is still experimental. Moreover it becomes clear that segmentation of thoracic organs is a complex and difficult task, and a scheme to perform such a segmentation automatically with minimal initialization has not been developed yet.
13
Chapter 3 Mathematical Background In the context of the present project, various mathematical tools and methods from the fields of image processing and numerical analysis have been employed. In this section we describe the background necessary for understanding and implementing the proposed approach.
3.1
Curve Fitting
In this section, mathematical approach of the most common curve fitting methods for the approximation of a the rib cage are presented. More precisely details of 2nd order polynomial approach and cubic spline methods are given as well as a comparison between these two approaches.
3.1.1
Least Squares Method
If the curve is not represented by a known function then an approximation of the curve in small intervals [(xi−N , yi−N ), ...., (xi+N , yi+N )] can be performed with a 2nd order polynomial (curve fitting)
φi (x) = c0 + c1 x + c2 x2
14
(3.1)
According to least Squares Method, coefficients c0 , c1 and c2 are defined calculating the minimum square error of Πi =
N X
(yi+j − φi (xi+j ))2
j=−N
. Hence, N X ∂Πi =0 ⇒2 (yi+j − (c0 + c1 xi,j + c2 x2i+j )) = 0 ∂c0 j=−N N X ∂Πi =0 ⇒2 xi+j (yi+j − (c0 + c1 xi,j + c2 x2i+j )) = 0 ∂c1 j=−N N X ∂Πi =0 ⇒2 x2i+j (yi+j − (c0 + c1 xi,j + c2 x2i+j )) = 0 ∂c2 j=−N
And finally we have the following linear system N X
N X
yi+j = 2N c0 + c1
j=−N N X
j=−N
xi+j yi+j = c0
j=−N N X
xi+j + c2
N X
xi+j + c1
j=−N
x2i+j yi+j
= c0
j=−N
N X
N X j=−N
N X
x2i+j
+ c2
j=−N
x2i+j
+ c1
j=−N
x2i+j
N X
N X
x3i+j
j=−N
x3i+j
+ c2
j=−N
N X
x4i+j
j=−N
The linear system can be solved easily using determinants ¯ PN PN 2 ¯ j=−N xi+j j=−N xi+j ¯ P 2N P P ¯ N N 3 2 D=¯ N j=−N xi+j j=−N xi+j j=−N xi+j ¯ PN PN PN ¯ x4 x3 x2 j=−N
i+j
j=−N
i+j
j=−N
i+j
¯ ¯ ¯ ¯ ¯ ¯ ¯
¯ PN PN PN 2 ¯ j=−N xi+j j=−N xi+j ¯ P j=−N yi+j P P ¯ N N 3 2 Dc0 = ¯ N j=−N xi+j j=−N xi+j j=−N yi+j xi+j ¯ PN PN PN ¯ x4 x3 x2 yi+j j=−N
i+j
j=−N
15
i+j
j=−N
i+j
¯ ¯ ¯ ¯ ¯ ¯ ¯
¯ PN PN 2 ¯ j=−N yi+j j=−N xi+j ¯ P 2N P P ¯ N N 3 Dc1 = ¯ N j=−N xi+j j=−N yi+j xi+j j=−N xi+j ¯ PN P P N N ¯ x2 x2 yi+j x4
¯ ¯ ¯ ¯ ¯ ¯ ¯
¯ PN PN ¯ yi+j j=−N xi+j ¯ P 2N PN PN j=−N ¯ N 2 Dc2 = ¯ j=−N xi+j j=−N xi+j j=−N xi+j yi+j ¯ PN PN PN 2 3 ¯ x x x2 yi+j
¯ ¯ ¯ ¯ ¯ ¯ ¯
j=−N
j=−N
So, c0 =
i+j
i+j
j=−N
j=−N
D c0 Dc , c1 = 1 D D
i+j
i+j
j=−N
j=−N
and c2 =
i+j
i+j
D c2 D
. In the rib cage approximation, as we have only a few data points and the intervals between them are big, we try to construct joined 2nd order polynomials. Hence, we construct joined polynomials to fit 3 data points each.
3.1.2
Cubic Splines
In this section we will describe an efficient way of constructing joined cubic splines through known data points. Briefly we could say that Cubic spline interpolation [31] is based on 4 basic conditions which are 1. Final curve pass through all specified data points (C0 continuity). 2. 1st derivative continuity(C1 ) at interior points. 3. 2nd derivative continuity(C2 ) at interior points. 4. Boundary conditions are specified at the free ends (for open curves only). In particular, considering the problem of constructing two cubic splines to fit three data points p1 = (x1 , y1 ), p2 = (x2 , y2 ) and p3 = (x3 , y3 ) with equations
16
1st Spline x1 ≤ x ≤ x2 f1 (x) = a1 (x − x1 )3 + b1 (x − x1 )2 + c1 (x − x1 ) + d1 f10 (x) = 3a1 (x − x1 )2 + 2b1 (x − x1 ) + c1 f100 (x) = 6a1 (x − x1 ) + 2b1 2nd Spline x2 ≤ x ≤ x3 f2 (x) = a2 (x − x2 )3 + b2 (x − x2 )2 + c2 (x − x2 ) + d2 f20 (x) = 3a2 (x − x2 )2 + 2b2 (x − x2 ) + c2 f200 (x) = 6a2 (x − x2 ) + 2b2 For the 1st Spline the 1st condition C0 is satisfied when y1 = f1 (x1 ) and y2 = f1 (x2 ) ⇒ y1 = d1
(3.2)
and y2 = a1 (x2 − x1 )3 + b1 (x2 − x1 )2 + c1 (x2 − x1 ) + d1 if we substitute x2 − x1 with h1 then y2 = a1 h31 + b1 h21 + c1 h1 + d1
(3.3)
In addition according to C2 condition y100 = f100 (x1 ) and y200 = f100 (x2 ) ⇒ y100 = 2b1 ⇒ b1 = y100 /2 and 17
(3.4)
y200 = 6a1 h1 + y100
(3.5)
⇒ a1 =
y200 − y100 6h1
(3.6)
An the second derivative becomes f100 (x) =
(x − x1 )(y200 − y100 ) + y100 h1
(3.7)
From 3.3 and 3.6,3.4 we have that y2 = h31
y200 − y100 + y100 h21 /2 + c1 h1 + y1 6h1
⇒ (y2 − y1 )/h1 = y200 /6h1 + y100 h1 /3 + c1 ⇒ c1 = (y2 − y1 )/h1 − y200 h1 /6 − y100 h1 /3
(3.8)
Finally, applying the C1 condition for the point p2 = (x2 , y2 ) such as y20 = f10 (x2 ) = f20 (x2 ) ⇒ 3a1 h21 + 2b1 h1 + c1 = c2 Substituting c1 , b1 , a1 and c2 = (y3 − y2 )/h2 − y300 h2 /6 − y200 h2 /3 where h2 = x3 − x2 we lead to
h1 y 00 y1 + 2(h1 + h2 )y200 + h2y300 = 6[(y3 − y2 )/h2 − (y2 − y1 )/h1 ] Hence in general for n number of points we have the system
18
(3.9)
1
0
0
··· ...
h1 2(h1 + h2 ) h2 0 h2 2(h2 + h3 )
h3 .. . .. . ...
0 0 0 0 0
0
hn−2 0
0
y 00 n−1 yn00
0 0
2(hn−2 + hn−1 ) hn−1 0 1 =
(y3 −y2 ) h2 (y4 −y3 ) h3
(yn −yn−1 ) hn−1
y100 y200 y30 .. . .. . .. .
0 − − .. . .. . .. .
−
=
(y2 −y1 ) h1 (y3 −y2 ) h2
(yn−1 −yn−2 ) hn−2
0 where hi = xi+1 − xi . Thus, solving the above system and estimating the nd
2 order derivatives of each point, we can calculate the joined cubic spline coefficients ai ,bi , ci and di .
3.1.3
Splines vs Polynomials
Polynomials are the approximating functions of choice when a smooth function is to be approximated locally. In general a polynomial provides a satisfactory approximation when the interval between the given points is sufficient small . But if a function is to be approximated on a larger intervals, the degree n of the approximating polynomial may have to be chosen unacceptably large. The alternative, as it was described earlier, is to subdivide the intervals into sufficiently smaller intervals so that, on each such interval, a polynomial of relatively low degree can provide a good approximation (joined polyno19
mial process). Although if the number of the given points set is low then polynomials fail to do an accurate curve fitting. Moreover, the discontinuities between the joined polynomials approximation give a finally no-smooth curve. For our rib cage approximating curve, c-spines is a good choice for curve fitting because of ease of data interpolation ,integration, differentiation and for smooth curve interpretation as well as we have only a small number of points (figure 3.1). 1.4
1.2
1
0.8
0.6
0.4 2nd order polynomial data points 2nd order polynomial Cubic Spline
0.2
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Figure 3.1: Splines vs Polynomials .
3.2
Optical Flow
Optical flow is a measure of the apparent motion of objects in a sequence of images. In the current project, optical flow is used for estimating the 3D displacement (volume motion) between inhale and exhale breath-hold images. This information can be extremely useful for the recognition of the tissues in the image because different biological tissues have distinctive motion features (magnitude and phase) during breath-cycle. For example, the back and the 20
spine are usually more static and motionless compared to the liver or other organs.
3.2.1
Mathematical Approach
Consider two volumes I1 (x, y, z) and I2 (x, y, z) defined on <3 which represent the two breath-hold stages. Under the assumption that corresponding pixels have equal grey values, the optical flow from I1 to I2 is defined as a displacement function h(¯ p) = [u(¯ p), v(¯ p), w(¯ p)] such that I1 (¯ p) = I2 (¯ p + h(¯ p)) ⇒ I1 (x, y, z) = I2 (x + u(x, y, z), y + v(x, y, z), z + w(x, y, z)) where p¯ = [x, y, z] ∈ <3 . The ordinary optical flow procedure [16, 3], linearizes the problem and approximates the solution using Taylor series expansion as I1 (x, y, z) = I2 (x, y, z) + u(x, y, z)
∂I2 ∂I2 ∂I2 + v(x, y, z) + w(x, y, z) +ξ ⇒ ∂x ∂y ∂z
I1 (x, y, z) − I2 (x, y, z) = h∇I2 (x, y, z), h(x, y, z)i This can be used only for small displacement h(x, y, z) and that the image (volume) is slowly varying in space. Otherwise this approximation is no longer valid. For any possible displacement and a well-posed definition of the problem, Alvarez et al. [16] considered the following minimization problem Z (I1 (x, y, z) − I2 (x + u(x, y, z), y + v(x, y, z), z + w(x, y, z)))2 dxdydz
EOF = <3
Z
|∇h|2 dxdydz.
+γ <3
(3.10) where the first term is the optical flow constraint and the second term smoothes the displacement field. Thus, large γ values lead to smoother flow 21
field. The estimation of h(x, y, z) which leads to the energy minimization EOF is based on the solution of the Euler-Lagrange equation. Thus, the associated Euler-Lagrange equations are given by the PDE system ∂I2 =0 ∂x ∂I2 γ4v + [I1 (x, y, z) − I2 (x + u(x, y, z), y + v(x, y, z), z + w(x, y, z))] =0 ∂y ∂I2 γ4w + [I1 (x, y, z) − I2 (x + u(x, y, z), y + v(x, y, z), z + w(x, y, z))] =0 ∂z γ4u + [I1 (x, y, z) − I2 (x + u(x, y, z), y + v(x, y, z), z + w(x, y, z))]
where 4 is the Laplace operator. The numerical solutions are obtained by calculating the asymptotic state t −→ ∞ of the parabolic system. ∂u ∂I2 = γ4u + (I1 (¯ x) − I2 (¯ x + h(¯ x)) ∂t ∂x ∂v ∂I2 = γ4v + (I1 (¯ x) − I2 (¯ x + h(¯ x)) ∂t ∂y ∂w ∂I2 = γ4w + (I1 (¯ x) − I2 (¯ x + h(¯ x)) ∂t ∂z where x¯ = (u(x, y, z), v(x, y, z), w(x, y, z)). The uniqueness and existence of the above minimization problem is proved in [17].
3.2.2
Initial Value Problem
Differential equations have multiple solutions (family of solutions) depending on the initial value choice. Initial data play an important role for the numerical calculation of a solution of a differential equation. Typically, the convergence is faster when the initial data is closer to the asymptotic state t −→ ∞. Therefore, for small displacements, u ≡ v ≡ w ≡ 0 comprises a rational choice. For large displacements however, a more sophisticated initialization is necessary. According to linear scale space theory [19], the introduction of a scale factor such as the gaussian filter Gσ can be used in the case of relatively large displacements. Hence the PDE system becomes
22
∂uσ ∂(Gσ ∗ I2 ) = γ4uσ + (Gσ ∗ I1 (¯ x) − Gσ ∗ I2 (¯ x + h(¯ x)) ∂t ∂x ∂vσ ∂(Gσ ∗ I2 ) = γ4vσ + (Gσ ∗ I1 (¯ x) − Gσ ∗ I2 (¯ x + h(¯ x)) ∂t ∂y ∂wσ ∂(Gσ ∗ I2 ) = γ4wσ + (Gσ ∗ I1 (¯ x) − Gσ ∗ I2 (¯ x + h(¯ x)) ∂t ∂z
(3.11) (3.12) (3.13)
where Gσ ∗ I is the convolution of I with a Gaussian filter of standard derivation σ. The convolution with the Gaussian filter blends the information in images (volumes) and provides a connection between the objects in I1 and I2 . The process starts with a large σ0 value in order to capture the largest expected displacements. Thus, the displacement (uσ0 , vσ0 , wσ0 ) is computed using the above PDE system with initial data u ≡ v ≡ w ≡ 0. Next, a number of scales σn < σn−1 < ... < σ0 which satisfy the σi = η i σ0 is chosen. For each scale σi the displacement is estimated applying the above PDE, with initial data (uσi−1 , vσi−1 , wσi−1 ).
3.2.3
Numerical Calculation
Discretization of the above parabolic system with finite differences gives
k uk+1 i,j,l − ui,j,l = γ0.25(uki+1,j,l + uki,j+1,l + uki−1,j,l + uki,j−1,l + uki,j,l+1 + uki,j,l−1 ) τ + [I1,σ (¯ xki,j,l ) − I2,σ (¯ xki,j,l + h(¯ xki,j,l ))]I2,x,σ (¯ xki,j,l + h(¯ xki,j,l ))
where I2,x,σ (¯ xki,j,l + h(¯ xki,j,l )) = 0.5[I2,σ (¯ xki+1,j,l + h(¯ xki+1,j,l )) − I2,σ (¯ xki−1,j,l + h(¯ xki−1,j,l ))]
23
and time step τ=
0.5 4γ +
maxi,j,l (|I1,x,σ (¯ xi,j,l )|2 , |I1,x,σ (¯ xi,j,l )|2 , |I1,z,σ (¯ xi,j,l )|2 )
The same equation is respectively for the v and w.
3.3
Euler Lagrange Equation
The definite integral Z I=
f (x, φ, ∇φ)dx
(3.14)
Ω
is a well defined quantity when Ω is a bounded region, when the integrand f is given as a function of arguments x, φ and ∇φ and φ is given as a function of x. The ”main” problem of calculus of variation is the estimation of the various values of integral I when different choices of φ(x) are substituted into the f . Actually what calculus of variations theory seeks, is the particular function φ which gives to I a minimum or a maximum value. In order to define the φ which gives an extremal to I, the most common method [15] is to create a one-parameter family functions and to try to determine the function which gives minimum. In particular, this family is expressed as u(x) = φ(x) + εh(x) where φ(x) is the extremizing function and h(x) is an arbitrarily selected continues function with continues first order derivatives and which is equal to zero at Ω boundaries. Thus substituting the u(x) to 3.14 we have that Z I(ε) =
f (x, φ(x) + εh(x), ∇φ + ε∇h(x))dx Ω
.
24
(3.15)
Then the necessary condition for I minimization is dI |ε=0 = 0 dε ⇒ dI = dε
Z µ Ω
¶ ∂f ∂f h+ ∇h dx ∂(φ + εh) ∂(∇φ + ε∇h)
For ε = 0 we have that dI |ε=0 = dε
Z µ Ω
¶ ∂f ∂f h+ ∇h dx = 0 ∂φ ∂∇φ
(3.16)
⇒ Z Ω
· ¸ µ ¶ µ ¶ Z Z ∂f d ∂f ∂f ∂f d ∇hdx = h − h dx = − h dx ∂∇φ ∂∇φ Ω ∂∇φ ∂∇φ Ω dx Ω dx
as h(x) = 0 in Ω boundaries. So 3.16⇒ µ
Z h Ω
∂f d ∂f − ∂φ dx ∂∇φ
¶ dx = 0
(3.17)
According to fundamental lemma of the calculus of variations, if a function K(x) is continues in Ω ⊂
where h(x) is an arbitrarily selected function which is continues and Ntimes continuously differentiable and equal to zero at the Ω boundaries then K(x) = 0 everywhere in the open subspace Ω. Thus, 3.17 ⇒ d ∂f ∂f − =0 ∂φ dx ∂∇φ
(3.18)
which is the well known Euler-Lagrange equation. For the case where 25
φ(x, y) : Ω ⊂ <2 → < is 2D Euler-Lagrange equations becomes ∂f ∂ ∂f ∂ ∂f − − =0 ∂φ ∂x ∂φx ∂y ∂φy
(3.19)
or in general m+1
∂f X ∂ ∂f − =0 ∂φ ∂x ∂φ i x i i=1
3.4
(3.20)
Naive Bayes Classifier
Bayes Classifier is based on the estimation of the posterior probabilities [28] P (xi |F1 , ..., Fn )
(3.21)
where xi are classes or outcomes of the dependant variable X, conditional on feature variables F1 ,...,Fn . Actually Bayes classifier determines the xi which corresponds to the maximum of the set of posterior probabilities 3.21. This decision criterion is called Maximum a posteriori probability or MAP. Using Bayes’ Theorem posterior probabilities can be expressed as P (xi |F1 , ..., Fn ) =
P (xi )f (F1 , ...Fn |xi ) f (F1 , ...Fn )
(3.22)
where f (F1 , ...Fn |xi ) is the conditional PDF of features vector given xi . Observing the equation 3.22 we realize that the computation of the posterior probability P (xi |F1 , ..., Fn ) requires only the a priori probability P (xi ) and the f (F1 , ...Fn |xi ). Furthermore, under the ”naive” conditional independence assumptions of the classifier where each feature Fi is conditionally independent of every other feature, equation 3.22 becomes P (xi |F1 , ..., Fn ) = P (xi )
n Y
f (Fj |xi )
j=1
As P (xi )f (F1 , ..., Fn |xi ) = f (xi , F1 , ..., Fn ) 26
(3.23)
and f (xi , F1 , ..., Fn ) = P (xi )f (F1 |xi )f (F2 |xi , F2 )...f (Fn |xi , F1 , F2 , ..., Fn1 ) = P (xi )f (F1 |xi )f (F2 |xi )...f (Fn |xi ) . So finally the decision rule or the MAP criterion can be expressed as argmaxX {P (X = xi )
n Y
f (Fj = fj |X = xi )}
(3.24)
j=1
Or in the case where P (X = xi ) are equiprobable then we have the Maximum Likelihood criterion(ML) n Y argmaxX { f (Fj = fj |X = xi )}
(3.25)
j=1
3.5
Level Sets Methods
In this section we present a brief description of the main level sets concept as well as some mathematical background as it was described by Osher and Sethian [27, 34, 26]. Generally, level set methods intend to fit closed evolving curves or surfaces to specific object boundaries. The level sets framework employs an implicit expression of the curve which enables the adaptation of the curve to various topological changes and even the detection of objects with sharp corners. It is also able to split or merge the curve and thus it is much more efficient and adaptive than the ordinary parametric contours or snakes [26, 23]. In detail, for a given open region Ω+ with smooth boundary we assume the existence of a 2D level set function φ(x, y, t), which is Lipschitz continuous, satisfying φ(x, y, t) > 0 (x, y) ∈ Ω+ (3.26) φ(x, y, t) = 0 (x, y) ∈ ∂Ω φ(x, y, t) < 0 (x, y) ∈ Ω− and the evolving curve defined as A(t) = {(x, y) ∈ Ω|φ(x, y, t) = 0}. A(t) is 27
called the zero level set or moving front. According to [27, 34], φ(x, y, t) evolves based on the Hamilton-Jacobi equation ∂φ + F |∇φ| = 0 (3.27) ∂t where F is defined as a speed function and depends on the image data and the level set function. An edge indicator is also usually applied to H-J evolutionary PDE in order to stop the evolution when the object boundaries are reached. Moreover, an external force v is added for increasing the speed of the numerical convergence. Thus, an edge detector g can be defined as a positive and decreasing function, depending on the image I(x, y) gradient and the evolution equation 3.27 ∂φ ∇φ = g(|∇I|)|∇φ|(div( ) + v) ∂t |∇φ|
(3.28)
in (0, ∞) × <2 where v > 0 is constant. During level set evolution, φ(x, y, t) can develop sharp edges or flat shape, which may lead to further inaccurate computations. For a smooth curve evolution, φ(x, y, t) is initialized to approximate a signed distance function d(~x). A distance function d(~x) : Ω → < is defined as d(~x) = min(k~x − x~I k)
(3.29)
for x~I ∈ ∂Ω. This implies that d(~x) = 0 on the boundary where ~x ∈ ∂Ω otherwise d(~x) = k~x − x~I k. Since the level set function must approximate a signed distance function, we have that |φ(~x)| = d(~x) for all ~x or φ(x) = d(~x) = 0 for ~x ∈ ∂Ω, φ(~x) = d(~x) for ~x ∈ Ω+ and φ(~x) = −d(~x) for ~x ∈ Ω− ( figure 3.5). Moreover according to 3.29 equation it is deduced that a signed distance function satisfies the |∇d| = 1 and thus |∇φ| = 1. Periodically during the evolution, φ(x, y, t) is re-initialized to approximate a signed distance function. In fact the re-initialization procedure for ordinary 28
?Ω Φ=0
ΩΦ<0
Ω+ Φ>0
Figure 3.2: Signed distance function φ(~x) defining the regions Ω+ and Ω− as well as the boundary ∂Ω. level set is a difficult, computationally expensive task as it is based on the solution of the ∂φ = sign(φ0 )(1 − |∇φ|) ∂t where φ0 is the function to be re-initialized and sign() is the sign function. However the re-initialization can be avoided solving a calculus of variation problem [5] as it will be described later in the “Front propagation without re-initialization” section.
3.5.1
The Chan-Vese Level Set Method
Main Idea The Chan-Vese (or without edges) level set method [6] intends to separate and segment the image into distinctive piecewise constant or piecewise smooth intensity regions, using the Mumford-Shah minimization functional [22]. The simplest case of this method, as we employ it in the present project, concerns the segmentation of the thoracic body volume into two regions with approximately piecewise constant intensity. This idea derives from the fact that the thoracic area consists of two intensity regions, one with low and close to zero values which correspond to the lungs and the other one with higher intensity values that includes the main body with the fat, the soft tissues and the ribs cage. In particular, assuming that an image I(x, y) is formed by two approxi-
29
mately piecewise constant intensity regions with intensity values c01 and c02 respectively and that the boundary between them is defined by a curve A0 , then I ≈ c01 inside A0 and I ≈ c02 outside A0 . Mathematically this can be expressed as the minimization of the Mumford-Shah functional Z
Z 2
inf (F1 (A) + F2 (A)) = inf ( A
A
|I(x, y) − c2 |2 dxdy)
|I(x, y) − c1 | dxdy + inside(A)
outside(A)
(3.30) where A is any other variable curve and the constants c1 and c2 are the average values of I inside and outside of A. Obviously, the minimization of 3.30 leads to the estimation of the desired A0 which separates the two regions. inf (F1 (A) + F2 (A)) ≈ 0 ≈ F1 (A0 ) + F2 (A0 ) A
. Thus, if the curve A is outside the object (i.e. outside the lungs), then F1 (A) > 0 and F2 (A) ≈ 0. If A is inside the object, then F1 (A) ≈ 0 and F2 (A) > 0. In the case that A is both inside and outside the object then both terms are greater than zero. Finally, the equation 3.30 is minimized when A = A0 (figure 3.3). Level Set Formulation of the Problem Chen and Vese [6] adapted the Mumford-Shah to the variational level sets, as it was described by Caselles et al. [36]. In particular, Caselles et al. replaced the unknown curve A(t) of the ordinary snakes approach [20] with the implicit function φ and solved the variational minimization problem leading to equation 3.28. This way, Caselles et al. managed to connect the classical “snakes” approach based on energy minimization [36] with the geometric active contour model based on the theory of curve evolution. Chen and Vese continued on this idea by adding the Mumford-Shah term to the level sets energy functional. Thus, the energy minimization functional was formulated according to 30
F1(A)>0, F2(A)~0
F1(A)~0, F2(A)>0
F1(A)>0, F2(A)>0
F1(A)~0, F2(A)~0
Figure 3.3: All possible cases in the initial position of the curve A .
Z inf {E(c1 , c2 , φ)} = inf {µ δ(φ(x, y))|∇φ(x, y)|dxdy c1 ,c2 ,φ Ω Z + v H(φ(x, y))dxdy ZΩ + λ1 |I(x, y) − c1 |2 H(φ(x, y))dxdy ZΩ + λ2 |I(x, y) − c2 |2 (1 − H(φ(x, y)dxdy)}
c1 ,c2 ,φ
(3.31)
Ω
where the first and the second term derives from the level set formulation [36] (only zero level set evolution) and the last two terms from Mumford31
Shah energy functional. More precisely, for the implicit function φ : Ω → < and for evoving only of the zero level set (i.e. φ = 0), the first two terms are defined as
Z
Z
Length[φ = 0] =
|∇H(φ(x, y))|dxdy = Ω
δ(φ(x, y))|∇φ(x, y)|dxdy Ω
Z Area[φ > 0] =
H(φ(x, y))dxdy Ω
where the first term intends to minimize the curve length and it actually smooths its surface while the second term minimizes the φ > 0 surface and works as an initial “force” which pushes the curve inwards or outwards depending on the sign of v. Moreover, the third and fourth Mumford-Shah terms can be expressed using the heaviside function as following Z
Z 2
|I(x, y) − c1 |2 H(φ(x, y))dxdy
|I(x, y) − c1 | dxdy = φ>0
Ω
and Z
Z 2
|I(x, y) − c2 |2 (1 − H(φ(x, y)))dxdy
|I(x, y) − c2 | dxdy = φ<0
Ω
where H(φ(x, y)) is the heaviside function ( H(φ) =
0 φ<0 1 φ>0
and δ(φ(x, y)) is the dirac where ( δε (φ) =
0 φ 6= 0 1 φ=0
32
Energy Minimization Keeping φ fixed and minimizing the functional E(c1 , c2 , φ) (equation 3.31) with respect to the constants c1 and c2 we have R c1 (φ) = if
R Ω
(3.32)
I(x, y)(1 − H(φ(x, y)))dxdy R (1 − H(φ(x, y)))dxdy Ω
(3.33)
H(φ(x, y))dxdy > 0 and R c2 (φ) =
if
I(x, y)H(φ(x, y))dxdy R H(φ(x, y))dxdy Ω
Ω
R
Ω
(1 − H(φ(x, y)))dxdy > 0. The minimization of E(c1 , c2 , φ) with respect to φ leads to the EulerLagrange equation estimation based on the calculus of variation (see Euler Ω
Lagrange section). Thus
δ(φ)[µdiv(
∂E =0⇒ ∂φ
(3.34)
∇φ ) − v − λ1 (I(x, y) − c1 )2 + λ2 (I(x, y) − c2 )2 ] = 0 |∇φ|
(3.35)
∀(x, y) ∈ Ω. Algorithm Steps The major steps of the Chan-Vese algorithm are • Initialize φ0 . • Compute c1 (φ) and c2 (φ) by 3.32 and 3.33. • Evolve the initial surface using ∂φ ∇φ = δ(φ)[µdiv( )−v −λ1 (I(x, y)−c1 )2 +λ2 (I(x, y)−c2 )2 ] (3.36) ∂t |∇φ| • Re-initialize φ locally to the signed distance function.
33
• Check if the solution is stationary. If it is not then repeat the procedure until a steady state ( ∂φ = 0) is reached, which is then an extremal of ∂t E.
34
3.5.2
Front propagation without Re-initialization: Level Sets Approach
Front propagation via level sets [32, 26] can be used for accurately extracting objects with well defined boundaries from 2D images. The efficient estimation of object boundaries using evolving fronts (A(t) or φ(x, y, t) = 0) comes with some basic limitations such that the initial curve should be inside or outside of the desired object and that we have only inwards or outwards curve evolution. In this section, we present a variational scheme of a modified zero level set or front propagation as it was described by Li [5]. In particular, φ evolution is based on the minimization of the following functional E(φ) = ηP (φ) + Em (φ)
(3.37)
where η > 0 is a parameter controlling P (φ) term which is defined as Z P (φ) = Ω
1 (|∇φ| − 1)2 dxdy 2
(3.38)
where Ω ⊂ <2 . P (φ) is actually a penalizing term, which intends to keep φ close to a signed distance function form, as a distance function should satisfy |∇φ| = 1. Em (φ) is called external energy and aims to “drive” the motion of the zero level set towards the object boundary. Em (φ) is defined as Z Z Eλ,κ (φ) = λ g|∇φ|δφdxdy + κ gH(−φ)dxdy Ω
Ω
where • δ is the Dirac function and is approximated by δε =
1 πx [1 + cos ( )], |x| ≤ ε 2ε ε
and H is the Heaviside function • g is the edge indicator function and is given by
35
(3.39)
g=
1 1 + |Gσ ∇I|2
where I is the image and Gσ =
1 e 2π
x2 +y 2 2σ 2
is the Gaussian mask.
• λ > 0 and κ constant. The coefficient κ can be positive or negative depending on the position of initial contour relatively to the object boundaries. If the initial contour is inside the object then κ < 0 otherwise κ > 0. According to the Steepest Descent method [5, 13, 12], the numerical minimization of the functional 3.37 is given by ∂φ ∂E =− ∂t ∂φ
(3.40)
with
∂E ∇φ ∇φ = −η[4φ − div( )] − λδ(φ)div(g ) − κgδ(φ) ∂φ |∇φ| |∇φ|
(3.41)
where 4 is the Laplacian operator and div is the divergence operator. Discrete forms of equations 3.41 and 3.40 can be applied for the front propagation implementation giving accurate results without the need of any reinitialization of φ as the P (φ) minimizes the divergence from a signed distance function. Moreover, the initial φ is no longer required to be a signed distance function as the P (φ) term forces it to approximate a signed distance function during the evolution. Implementation The approximation of equation 3.40 is the first step for the implementation of the front propagation. Equation 3.40 can be written in a difference scheme as k k φk+1 i,j = φi,j + τ L(φi,j )
36
(3.42)
where L(φki,j ) is the discretization of the right hand side of equation 3.41 and τ is the time step which must satisfy τ η < 1/4. Moreover we have to assign an initial level function φ0 for the iteration procedure 3.42. The φ0 can be determined by the circle with center (ai , bi ) and radius Ri . If an initial level set front is defined as A(0) = {(x, y)|(x − ai )2 + (y − bi )2 = Ri2 } then the initial level set function can be described by −ρ0 (x, y) ∈ Outside of A(0) φ0 (x, y) = 0 (x, y) ∈ A(0) ρ0 (x, y) ∈ Inside of A(0)
(3.43)
The coefficient κ of equation 3.39 can be positive or negative depending on the relative position of initial contour to the object boundaries. In our method we will assume that κ < 0 and then the initial contour must be inside the object for proper level set evolution. Applying equation 3.42 for k = 0 : N will reach the boundaries when −1 N φi,j − φN = 0 is satisfied. i,j
37
Chapter 4 Proposed Method The presented method consists of a pre-processing step, where we apply standard image processing techniques in order to improve the quality of the data, and two major parts as shown in figure 4.1.
CT slice
Pre-prossesing
1st part Rib cage Approximation
2nd part
Lungs & Heart “rough” estimation
refined estimation
Level Set Method
Segmented image
Figure 4.1: Flow Chart In the first part we have a coarse distinction of the body into two sub38
regions. The one subregion consists the outer thoracic part with the bones and the fat tissues and the other subregion (the inner thoracic part) with the lungs, the heart and some other tissues and structures like the aorta. In detail, the “rough” separation into the two subregions is performed by estimating a curve which approximates the rib cage. A simple threshold selection technique is applied to segment the rib cage and c-spline interpolation is used to define the curve. In the second part, we apply a coupled level set procedure in order to compensate for the “rough” separation of the first part and accurately extract the whole thoracic cavity. The level set framework employs a two phase scheme for the capture of the lungs and soft tissues (i.e. the heart and some other small tissues). The two-phase level set method intends to segment two objects (the lungs, the heart) and thus two zero level set functions should be defined, one for the lung and one for the heart. Initialization of the two zero level sets is achieved using the rib cage approximating curve. For the initialization of the zero level function which will extract the heart we also have to define an initial small region inside the rib cage where the heart is positioned. For this purpose, a naive Bayes classifier is used to define which of the grey value pixels inside the estimated rib cage are most probable to belong to the heart. The pixels with the highest probability are used for the heart zero level function initialization. In the case of the lungs, the zero level function is initialized using the rib cage curve and the heart zero level set function. In the following sections, we describe in detail each part of the method.
39
4.1
Pre-Processing
The pre-prossing stage employs a median filter followed by a contrast enhancement scheme. Pre-processing is very important as it reduces noise and smooths homogenous parts of the CT images.
4.1.1
Median Filter
The 2D median filter is defined as y(i, j) = med{x(i + r, j + s)|(r, s) ∈ A, (i, j) ∈ Z 2 } where med{} expresses the median value of the filter window and Z 2 = Z × Z determines the discrete image plane. A ⊂ Z determines the size of the filter window. Median filter is preferable in our case as it reduces speckle noise and salt and pepper noise. Moreover its edge-preserving property makes it useful in cases where edge blurring is undesirable.
4.1.2
Contrast enhancement
Contrast enhancement [37, 30] of the image is applied before the filtering. In particular two threshold intensities (ilower , iupper ) are selected. All pixels with intensities below ilower are set to the minimum intensity imin (0 in our case), and all pixels above the iupper are set to the maximum intensity imax (255 in our case). The rest of the pixels get intensity values according to ipixN ewV al =
ipixel − ilower imax iupper − ilower
where ipixel is the intensity of the pixel to be recalculated. The upper and lower threshold intensities, iupper and ilower , are given by ilower = tlower imax and iupper = tupper imax where tlower and tupper are values between 0 and 1 [30]. Adjusting contrast generally helps to discern a specific area of the image.
40
Figure 4.2: re-scaling graph
4.2
Rib Cage Approximation
The fist step of the proposed method is the automatic “rough” separation of the inner body part (i.e. soft tissues like heart and lungs) from the skeleton and the skin. Because ribs are visible (higher intensity value) in all CT testing data, and define a natural border between the inner thoracic area and the outer body part, we can achieve a coarse separation using a simple threshold selection for the ribs segmentation, followed by a curve fitting process. In the following section, we describe step by step the process for the calculation of the approximated rib cage curve. • Image Thresholding: The selection of an appropriate threshold which segments accurately the ribs and the spinal column is based on the fact that bones are the biological structures with the highest luminance in the thoracic CT volumes. Moreover, we can assume that the bones correspond to approximately the 2 ∼ 3% of the total number of grey values of the CT volume and thus we can express it in probabilistic terms as nbones P (X ≥ T hresh) = = 0.03 ntotal where X is the grey values variable.
41
It is known that Z
T hresh
P (X ≥ T hresh) = 1 −
fX (u)du
(4.1)
0
where fX () is the probability density function (PDF) of the grey values. The above PDF can be derived from the normalized volume histogram. Therefore, estimating the histogram of the CT volume and dividing it with the total number of volume pixels, we obtain the PDF fX (x). Finally, solving equation 4.1 we determine the threshold.
(A)
(B)
Figure 4.3: (A)P (X ≥ T hresh) and (B) Segmented Image . In the discrete case, equation 4.1 is approximated by
P (X ≥ T hresh) = 1 −
TX hresh i=0
ni ntotal
= 0.03
where ni is the number of pixel for each grey intensity value. • Application of “open” morphology operation in order to reduce the number of single bright pixels and calculation of the barycenter of the 42
rib cage in each cross-sectional slice according to PN R=
i=1
ri
N
where ri are the coordinates of the bright pixels and N the total number of the bright pixels.
φ
(A)
(B)
Figure 4.4: (A) Estimation of minimum distance in each angular direction and (B) pixels with the minimum distance . • Calculation of the minimum distance between the barycenter and the ribs (bright pixels of the segmented image) in each angular direction (figure 4.4(A)). • Curve fitting process using cubic spline interpolation (figure 4.5).
43
(A)
(B)
Figure 4.5: (A)Result of curve fitting and (B) coarse separation of thoracic region .
4.3
Heart Position
In the previous section, the rib cage approximation steps defined a curve which roughly separates the thoracic volume into two subregions. Indeed, as the separation is coarse, the one subregion (outer subregion) includes the ribs, skin, fat a probably small lungs parts and the other consists of the rest of the lungs, the heart, some other tissues (i.e aorta etc) and possibly some other small area which should belong to the outer subregion. The definition of the lungs area in the inner subregion is easy as lungs are characterized by low luminance. However, the heart which has approximately the same intensity as other tissues is more difficult to be distinguished. For this purpose, a classification scheme is used for the definition of the heart position among the other gray value tissues which are inside the estimated rib cage. According to the theoretical background of section 3.4, we have the dependent variable X with two outcomes X = {x1 = “is heart00 , x2 =
44
“is not heart00 } and the feature vector F¯ = {F1 , F2 , F3 , F4 } ,where F1 is the pixel intensity and F2 , F3 , F4 are the displacements towards x, y and z respectively between the inhale and exhale volume as they calculated by the optical flow of section 3.2. In addition, assuming that the P (X = xi ) are equiprobable and that the conditional PDFs f (Fj = fj |X = xi ) follow the gaussian distribution and thus they can be expressed as
f (Fj = fj |X = xi ) =
1 exp 2πσFj |xi
−
(fj −µF |x )2 j i 2σ 2 Fj |xi
(4.2)
we end up to the ML criterion. Since in our case we have two classes, the ML criterion can be expressed as Q4 j=1
Q4
f (Fj = fj |X = x1 )
j=1 f (Fj = fj |X = x2 )
≷xx12 1
(4.3)
For the estimation of the mean values and the standard deviations of each feature distribution, we segmented manually the heart and we defined the rest non-zero pixels which are inside the rib cage, then we calculated the means and the standard deviations of the intensity values and the displacement for both regions. The training data was extracted from the first patient. Actually, from the whole volume [392 × 297 × 237] we used only the cross sectional slice multiples of 5 between [50 : 220].
4.4
Body Volume Definition
In the rest steps of the method all the processing is made only in the body volume. For the elimination of the background air (black area outside the volume), we apply a filling algorithm. Filling algorithms define 4 or 8-connected regions. After rescaling, low intensities become zero and thus the background is now zero. Using a filling algorithm, we define 4-connected regions starting from the boundary (seed points) of the image and thus we determine the background area (figure 4.7) and respectively the image subspace where
45
body is located.
Figure 4.6: Filling process, 4-connected region
(A)
(B)
Figure 4.7: (A) Mask: result of filling process and (B) the CT slice In the rest of the text Ω will denote the image subspace (body) defined by this process. Background elimination decreases the computational cost and gives better results with the level set method.
46
4.5
Thoracic Cavity Extraction
For the extraction of the inner thoracic area (figure4.8-C), a level set framework based on the combination of the Chan-Vese level set method and the Front Propagation scheme is applied.
A
B
C
Figure 4.8: (A)CT slice, (B) outer part and (C) inner part. In particular as the subregion inside the thorax consists of biological structures with different properties (i.e intensity values, structure morphologies, size etc) a simple level set function φ : Ω → < as defined in section 3.5 cannot capture and separate this area from the outer body “shell” (skeleton-fat). The definition of a single evolutionary partial differential equation which includes criteria and penalizing factors for a complex region such as the thorax does not exist and it is impossible to be designed at least till now. However, considering that the thoracic organ can be actually separated into a low intensity subregion, i.e. the lungs, and a higher intensity subregion, i.e. the heart and some other soft tissues like the aorta, a level set approach employing two level set functions can be designed. Each level set function will have different properties and will evolve with different PDE. Hence, this scheme can accurately segment these different structures (i.e. Lung, Heart etc.). As it was described earlier in section 3.5, the Chan-Vese (CV) level set equation enables the separation of a region into two piecewise constant intensity value subregions. Thus, this method can separate the lung area from 47
the rest volume. On the other hand, the Front Propagation method is able to segment an object using an initial curve which stops only when boundaries are reached. So, the definition of an initial curve inside the heart and its evolution according to 3.41 equation can give the heart segment. In the following section, we will explain the implementation of the level set method, the initialization of the level sets functions and the additional necessary penalizing factors.
4.5.1
Levels Sets Functions Initialization
The level sets function initialization is an extremely important task and most of the previous work relies on manual initialization. In our scheme, the initialization is automatic and is based on the estimation of the rib cage and the position of heart. Let assume that the rib cage curve is defined by the equation C1 (x, y) = 1 where (x, y) ∈ Ω and that the image point with the highest probability to belong in heart is the (xi , yi ) (classification result). Using the point (xi , yi ) we can draw a small curve C2 = {x, y ∈ Ω|(x−xi )2 +(y −yi ) = r} with r w 2−3 pixels radius. Then we can define the initial level sets functions φ1 (x, y, 0) and φ2 (x, y, 0), which correspond to the lungs and the heart (possibly including some other soft tissues too) respectively, using C1 and C2 . In particular, we have −ρ0 (x, y) ∈ Outside of C2 φ2 (x, y, 0) = 0 (x, y) ∈ C2 ρ0 (x, y) ∈ Inside of C2
(4.4)
Determining implicitly the C1 curve as −ρ0 (x, y) ∈ Outside of C1 R(x, y, 0) = 0 (x, y) ∈ C1 ρ0 (x, y) ∈ Inside of C1
(4.5)
we can define the φ1 (x, y, 0) as φ1 (x, y, 0) = −sign(R(x, y)φ2 (x, y, 0))ρ0 (figures 4.9 and 4.10) where sign() is the sign function. 48
Figure 4.9: Graph of initial zero level set functions .
(A)
(B)
Figure 4.10: (A)φ2 level set function and (B) φ1 level set function .
4.5.2
2-Phase Level Sets Implementation
In the present section we define how we are going to combine the Chan Vese and Front Propagation Level Sets schemes of section 3.5 in order to estimate and separate the inner body part from the rib cage and the fat-skin.
49
As we have already defined the initial level set functions φ1 and φ2 in the previous section, we have to write the evolutionary differential equation as well as any restriction applied. The two evolutionary equation, see section 3.5, are
∂φ1 ∇φ1 ∇φ1 = δ(φ1 )[µdiv( )−v−λ1 (I(x, y)−c1 )2 +λ2 (I(x, y)−c2 )2 ]−η[4φ1 −div( )] ∂t |∇φ1 | |∇φ1 | (4.6) for the lungs extraction and
∂φ2 ∇φ2 ∇φ2 = −η[4φ2 − div( )] − λδ(φ2 )div(g ) − κgδ(φ2 ) ∂t |∇φ2 | |∇φ2 |
(4.7)
for the heart and soft tissue where g is the edge indicator as it is described in Front Propagation. In order to prevent any “leakage” of φ2 (figure 4.11-A) during the evolution to the surrounding tissues which are outside the rib cage and they have quite similar intensities with the heart, we introduce a criterion to keep the curve A2 (t) = {x, y|φ2 (x, y, t) = 0} inside the rib cage. This criterion ensures that 4.7 equation evolves only inside the rib cage area and can be expressed as ∂φ1 ∂φ1 = h(R(x, y)) ∂t ∂t where R(x, y) the implicity expression of the C1 (rib cage curve) according to equation 4.5 and h() is the heaviside function ( h(φ) =
0 φ<0 1 φ>0
. In addition, we multiply each equation with the (1 − h(φ1 )h(φ2 )) term in order to reduce any possible overlapping between the two evolutionary function .
50
Zero Level Set
Zero Level Set
(A)
(B)
Figure 4.11: evolution of φ2 (A) without rib cage stopping criterion and (B) with stopping criterion. Zero Level Set
Zero Level Set
(A)
(B)
Figure 4.12: evolution of φ2 (A) without (1 − h(φ1 )h(φ2 )) term and (B) with (1 − h(φ1 )h(φ2 )) term. Numerical Approximation In this section we present the numerical approximations of the first, second derivatives as well as the dirac and heaviside functions of the evolutionary 51
PDEs. The 2D first order partial derivatives based on forward difference scheme are
∂φ φi+1,j − φi−1,j = ∂x ∆x ∂φ φi,j+1 − φi,j−1 = ∂y ∆y
(4.8) (4.9)
while the 2D second order derivatives are ∂φ2 φi+1,j − φi,j + φi−1,j = 2 ∂x ∆x2
(4.10)
∂φ2 φi,j+1 − φi,j + φi,j−1 = 2 ∂y ∆y 2
(4.11)
φi+1,j+j − φi+1,j−1 + φi−1,j−1 − φi−1,j+1 ∂φ2 = ∂x∂y ∆x∆y
(4.12)
The Heaviside and the dirac approximations are respectively 0 Hε (φ) =
1 2
φ < −ε +
φ 2ε
+
1 2π
sin
1
0 δε (φ) =
1 2ε
+
1 2ε
cos
0
πφ ε
πφ ε
−ε ≤ φ ≤ ε φ>ε
φ < −ε −ε ≤ φ ≤ ε φ>ε
(4.13)
(4.14)
∇φ And the curvature kφ or div( |∇φ| ) is
φyy φ2x − 2φxy φxφy + φxx φ2y kφ = |∇φ|3/2
(4.15)
Hence the numerical approximation of the PDEs can be based on these approximation. Particularly, for the equation 4.6 we can write
φk =φk−1 + τ [1 − hε (φk−1 )hε (φk−1 )]{µδε (φ)kφk−1 −vδε (φ) − λ1 δε (φ)(I(x, y) − c1 )2 + λ2 δε (φ)(I(x, y) − c2 )2 k−1 − η(φk−1 xx + φyy − kφk−1 )}
52
and for the equation 4.7
k−1 φk =φk−1 + τ [1 − hε (φk−1 )hε (φk−1 )]hε (R){−η(φk−1 xx + φyy − kφk−1 ) φx φy −λδε (φk−1 )(gkφk−1 + gx 2 + gy 2 ) − κgδε (φk−1 )} 2 1/2 (φx + φy ) (φx + φ2y )1/2
. Algorithm Overview • Initialize φ01 and φ02 . • Compute c1 (φ01 ) and c2 (φ01 ) from 3.32 and 3.33 and the heaviside functions of φ01 , φ02 and R(x, y). • Evolve the initial surfaces applying the
∂φ1 ∇φ1 = [1 − h(φ1 )h(φ2 )]{µδ(φ1 )div( ) − vδ(φ1 )− ∂t |∇φ1 | λ1 δ(φ1 )(I(x, y) − c1 )2 + λ2 δ(φ1 )(I(x, y) − c2 )2 − η(4φ1 − div(
∇φ1 ))} |∇φ1 |
and ∂φ2 ∇φ2 ∇φ2 = (1 − h(φ1 )h(φ2 ))h(R){−η[4φ2 − div( )] − λδ(φ2 )div(g ) − κgδ(φ2 )} ∂t |∇φ2 | |∇φ2 | . • Check if the solutions are stationary. If they are not then repeat the 1 2 procedure until a steady state of ( ∂φ = 0) and ( ∂φ = 0) is reached. ∂t ∂t When the iterative process ends, the thoracic cavity is defined as the common area of the two expanding curves. According to the definition of levels sets (see section 3.5) since φ1 < 0 and φ2 < 0 outside the area of interest and positive inside the lungs and the heart, we can describe the
53
Zero Level Set
Zero Level Set
A
B
C
Figure 4.13: Initial zero level sets, (B) final zero level sets (500 iterations) and (C) segmented inner body part thoracic cavity as the function ( φ(x, y) =
0 φ1 φ2 < 0 1 φ1 φ2 ≥ 1
(4.16)
In order to ensure that there is not any pixel missed out between the heart and the lungs, we apply a filling holes algorithm. Hence, we eliminate any single or any small group of zero values which are inside the φ > 0. Subsequently, the boundary curve of the segmented cavity is defined according to C = |∇φ| .
54
Chapter 5 Experimental Results To verify the theoretical preview, a description of the proposed method is presented in this chapter giving precise results and mentioning further details about the choice of the constants, initial conditions, filter size and any other constraint which was considered.
5.1
Experiments
Pre-processing In the pre-processing stage the images are enhanced by applying an image contrast adjustment according to the re-scaling scheme (fig.4.2) with tlower = 0.2 and tupper = 0.87 and a 2D [5 × 5] median filter. Figure 5.1 depicts the result of the preprocessing. Smoothing is very important for our method, especially when we try to accurately extract the lungs from the whole volume using the idea of ChanVese model, as we need to have two piecewise smooth subregions. Rib Cage Approximation The estimation of the threshold is made using the re-scaled volume and not the smooth result of the median filtering in order to keep the bones intensities in high values and make easier the threshold selection. In figure 5.2 the threshold was 969 for the unit16 image. 55
(A)
(B)
Figure 5.1: (A) Initial Image and (B) Enhanced Image, Median filtering 255 ) where RMS is the mean least square error between (P SN R = 20 log10 ( RM S the initial and the filtered image).
(A)
(B)
Figure 5.2: (A) Thresholded Image and (B) Coarse Separation. After the threshold selection and the definition of the points which are closer to the center of mass in each direction, we “draw” a c-spline that fits to the estimated points and we this way have a first separation (figure 5.2-B). 56
Features estimation and Classification Result For the classification of the grey pixels which are inside the rib cage curve the first step is to estimate the features values (intensity and 3D displacement). The displacement values are estimated using the optical flow of section 3.2. Figure 5.6 shows the result of optical flow in the transverse plane with the smoothing constant (section 3.11) set to γ = 0.1. The result of the classification scheme is displayed in figure 5.4.
Figure 5.3: Displacement estimation of rigid body parts The classification is important for the initialization of the φ2 level set function (function which segments heart) and for the rest process. Incorrect initialization of the heart obstructs the accurate level set evolution. In fact, there are cases where the position of the heart is not in middle of the rib cage and some other times where there are small nodes in the lungs area which can be confused with the heart tissue and can lead to a wrong initialization. In addition, the proposed classification scheme is the first step of a more complete classifier which will includes the liver and other organs 57
(A)
(B)
Figure 5.4: (A) Heart position according to Classification and (B)Classification Map: magenta is the heart regions, light violet the rest gray value region and in light blue the background and the decision of the classifier will determine the future evolutionary level set scheme. As we saw in the previous chapter, there is a penalizing factor which stops the evolution of the heart level sets function inside the rib cage, although a different factor should be used for the case where we are in liver. Definition of the Body Volume Using a filling algorithm, we create a binary mask (Background=1, Body=0) in order to perform further calculations only on the area of the body. Since the background is a connected region with zero intensity, we can easily define the background starting from the image boundaries and stopping when the body is reached (e.g. the intensity in greater than zero). An example of the extracted body mash is shown in figure 5.5. Level Sets Results The Level Set method is applied in every image for the accurate extraction of the inner thoracic part. The initial Level Set functions φ01 and φ02 are constructed according to the functions 4.4, 4.5 and φ1 (x, y, 0) = −sign(R(x, y)φ2 (x, y, 0))ρ0 with ρ0 = 5. The initial zero level of φ02 is a 58
Figure 5.5: Body mask circle with center defined by the pixel with the maximum posterior probability to be heart and a small radius. The zero level of φ01 is determined by the rib cage curve. The initial level function is not necessary to approximate ∇φ ) added in the a distance function due to the penalizing factor 4φ − div( |∇φ| evolutionary PDE (for further details read the front propagation subsection 3.5). The evolution of φ1 is based on the equation 4.6 and the φ2 on the 4.7. In our experiments, we generally choose the parameters as follows • For the equation 4.6 the parameters are λ1 = λ2 = 1, v = 5 and µ = 16. λ1 and λ2 factors define the contribution of the Mumford-Shah term to the final result. v works as an initial “force” to push the initially curve inwards. µ parameter has a scaling roll in PDE, for small values can capture small objects and for higher values large objects or groups of objects, In our experiments, we choose this value to be small enough to capture structure details and big enough to neglect nodes or small gray value fluctuations inside the lungs. • For the equation 4.7 λ = 5 and κ = −3.5. We choose κ < 0 as the initial zero level is inside the object boundaries and as we have referred 59
in section 3.5 we want outwards expansion. The Gaussian operator Gσ for the estimation of g (edge indicator function) has σ = 0.9 and dimension [5 × 5], which are quite good values for image smoothing without loosing big image details. • For the regularized functions δε (x) and hε we set ε = 1.2, η = 0.19/τ and the time step τ = 5 which ensures stability and fast convergence. Figure 5.6 shows the initial contour on the CT slice and the result after approximately 600 iterations. Zero Level Set
Figure 5.6: Top left: Initial zero levels, top right: Zero level (iter=600) , bottom left: Final contour, bottom right: Segmented regions
60
5.1.1
Results and Discussion
The proposed method has been tested using CT images of real patients. For the validation of our result we consider • True Positive (TP): the pixels of the slice which belong to the desired region (Lung or Heart) and were segmented automatically by our method. • False Positive (FP): the pixels of the slice which do not belong to the region but they were wrongly segmented by our method. • False Negative (FN): the pixels of the slice which belong to the desired region but they were not segmented. • True Negative (TN): the pixels of the slice which do not belong to the region and correctly were not segmented by our method. And we analyze our results using the following performance measures Sensitivity =
num of T P num of T P + num of F N
Specif icity =
num of T N num of T N + num of F P
where P = num of T P + num of F N is the number of pixels segmented manually for the heart and applying the regions growing with manual threshold choice for the lungs. and T A B DICE = 2 |A| + |B| where A is the group of the automatic segmented pixels and B the manually segmented for the heart and using region growing for the lungs.
61
Tables and figures show the results for the three patients with volumes sizes [392 × 297 × 237], [387 × 295 × 225] and [417 × 274 × 355] for the 1st , 2nd and 3rd patient respectively. The selected slices are between the diaphragm and the upper part of the heart and the slices’ numbers start from the lower body part to the upper part. The validation for the heart was made by comparing it with manually segmented hearts and for the lungs by comparing with region growing results. First patient The tables 5.1 and 5.2 show the performance of the method for slices of the first patient whereas figures 5.8, 5.9 and 5.7 the respective estimated segmentation results.
Figure 5.7: 1st patient, Final rib cage curves (A) slice 84, (B) slice 102, (C) slice 161 and (D) 174.
62
Figure 5.8: 1st patient, Colomns:(A)Initial contours, (B)Final contours.
63
Figure 5.9: 1st patient, Colomns:(A) Lung segmentation using proposed method (green) vs Region Growing (Red) and (B) Heart segmentation, automatic(blue)-manual(magenta).
64
Table 5.1: Results of Heart Segmentation for the 1st patient 1st Patient slice number Sensitivity Specificity DICE iterations time in sec. 84 0.9805 0.967 0.7983 900 58 102 0.98 0.9837 0.9114 700 44 161 0.9702 0.9907 0.86 700 42 172 0.9491 0.9918 0.8525 700 41
Table 5.2: Results of Lungs Segmentation for 1nd Patient slice number Sensitivity Specificity 84 0.9582 0.9941 102 0.9762 0.9927 161 0.9658 0.9940 172 0.9760 0.9955
the 1nd patient DICE 0.9631 0.9631 0.9672 0.9742
In the above figures and tables, we see the performance of the proposed method for slices of different parts of the thoracic cavity. Figure 5.8 depicts the initial zero level sets and the final contours. Figure 5.9-A compares actually the region growing segmentation with the modified chan-vese level set result and the figure 5.9-B the manual segmentation with the proposed automatic method. Finally, figure 5.7 presents the final curves.
65
Second patient For the second patient we estimated the sensitivity, specificity and DICE for the slices [55 : 128] between the diaphragm and the upper part of the thorax . The validation measurements for the second patient are presented in the following figures: 1
0.98
0.96
0.94
0.92
0.9
0.88
0.86
0.84 Sensitivity DICE Specificity
0.82
0.8 50
60
70
80
90
100
110
120
130
Number of slice
Figure 5.10: 2nd patient: specificity, sensitivity and DICE for the segments segment.
Table 5.3: Mean value and standard dev. of Specificity, Sensitivity and DICE Heart Lungs Mean Std Mean Standard deviation Sensitivity 0.972 0.0294 Sensitivity 0.978 0.0021 Specificity 0.9716 0.0155 Specificity 0.986 0.0042 DICE 0.904 0.1149 DICE 0.964 0.0061 According to figures 5.10, 5.11 and tables 5.3 and 5.3, the sensitivity, the specificity and the DICE are close to 1 and ensure that the segmentation is quite close to manual result (for the heart) and to region growing approximation (for the lungs). Although we can observe that there are some slice with no so good results. This happens (like in slice 80 figure 5.12) because the rib cage approximation is not accurate and has as a result to obstruct 66
1
0.99
0.98
0.97
0.96
0.95
0.94
0.93
0.92
0.91 50
Specificity Sensitivity DICE 60
70
80
90
100
110
120
130
Number of slice
Figure 5.11: 2nd patient: specificity, sensitivity and DICE for the lungs segment. the level set evolution outside the rib cage and to extract the whole heart and thus the measurements are not satisfactory.
Figure 5.12: 2nd patient, slice 80-no satisfactory rib cage approximation. Here, we have to mention that the results of figures 5.10 and 5.11 depend on how reliable are the region growing and the manual segmentation. The manual segmentation of the heart especially in the upper thoracic part, where there is the trachea and some protrusions around the heart volume (aortic tubes or lungs airways) is quite difficult, therefore the comparison of our method with the manual segmentation is quite hard.
67
Third patient For the third patient we chose some images from different parts of the thorax. As we can observe, the results show (figure 5.13) that the segmentation process is better in the middle of the thorax while in the lower part (slice 50) is more inaccurate. In the lower part the rib cage approximation is not very good as there is not rib in the front side (chest wall). Moreover, in the second and third row of the figure 5.13-A we see that the initial rib cage curve (red line) may not include the whole lungs (slice 85) or may includes parts of the outer thoracic body (slice 215), something which not influence the level set evolution.
Figure 5.13: 3rd patient, Colomns:(A)Initial contours, (B)Final contours and (C) Final rib cage curve. 68
Table 5.4: Results of Heart Segmentation for the 3rd patient 3rd Patient slice number Sensitivity Specificity DICE iterations time in sec. 64 0.98 0.98 0.94 1100 61 85 0.98 0.99 0.964 850 45 215 0.985 0.988 0.89 800 42
Table 5.5: Results of Lungs Segmentation for 3rd Patient slice number Sensitivity Specificity 64 0.971 0.99 85 0.964 0.99 215 0.939 0.974
the 3rd patient DICE 0.97 0.973 0.952
In conclusion, we can say that the results are quite promising since the sensitivity and the specificity are very close to 1. Furthermore, it is obvious from table 5.1 that the algorithm’s computational time is directly related to the number of the level set iterations. The different number of iterations in every data set depends on the different size of the inner tissue. For instance, in the first line of 5.8, where we are close to the diaphragm, the level set needs more iterations to accurately converge to these boundaries. Moreover, in the absence of ribs, such as in slice 65 of the third patient in figure 5.13, the initial wrong rib cage approximation prevents the evolution of the φ2 (green line) and the whole soft tissue extraction. Finally, we have to mention that the boundary of the curves are not very smooth, this happens as the Chan-Vese algorithm tries to separate two regions with piecewise constant intensity (lungs and the rest tissues). However, lungs and heart are not two piecewise constant, and thus the algorithm can not fit accurately. There are a couple of ways to solve this problem. The more sophisticated is to define the problem for two piecewise smooth regions [7] although a simple smoothing filtering can enhance the boundaries or the addition of an edge indicator to the Chan-Vese PDE may improve the performance.
69
Chapter 6 Conclusions and Future Work Current applications of segmentation require the design of algorithms with great accuracy and minimum computational cost. The majority of the previous algorithms imposes a significant manual interaction. The proposed method eliminates the manual initialization or interaction without the need of additional memory requirements or any performance restrictions. This way the inner thoracic region can be automatically extracted from the CT images for further examination or for a non-linear registration processing. Moreover, the achieved individual organs segmentation (i.e. lungs, heart) could be the first step for a different radiotherapy planning treatment where each organ can be “tracked” separately. Several improvements can be made to the segmentation scheme some of which require only more time to implement, while others depend on further investigation in this field. In particular: • The robustness of the method can be improved by further testing it on more data sets and adjusting the various parameters to better fit a wider variety of input. • Further improvement of the classification scheme. Use of more features and criteria which can provide more efficient results. Implementation of more sophisticated classification schemes like SVMs. More training data and more classes for classifying not only the heart but also other biological structures. 70
• Semi-implicit level set implementation (much Faster) and Shape knowledge level set method for liver estimation as the boundaries between liver and body fat are obscure. • A limitation of this study is the 2D implementation. So, 3D Initial meshes and 3D level set implementation is the next step. • Extending the current method and combining it with the deformable registration scheme as proposed by [10],[9]. The main goal of this project was to perform segmentation of the inner thoracic part as seen in CT images, without any kind of manual intervention, using the provided mathematical tools and the Matlab-mex software. The results show great potential and the method is dramatically faster than manual segmentation while still giving reliable results. In all it is a significant step in the direction of automatic thoracic segmentation and can provide the basis for future implementations of medical imaging tools for the extraction of essential clinical information.
71
Bibliography [1] Automatic Segmentation and Registration of Lung Surfaces in Temporal Chest CT Scans, volume 523-2005 of Pattern Recognition and Image Analysis. Springer Berlin - Heidelberg, 2005. [2] A. S. El-Baz A. M. Ali and A. A. Farag. A novel framework for accurate lung segmentation using graph cuts. Biomedical Imaging: From Nano to Macro, 2007. ISBI 2007. 4th IEEE International Symposium, 2007. [3] G. Aubert and P. Kornprobst. Mathematical Problems in Image Processing-Partial Differential Equations and Calculus of variations. 2000. [4] M. Ramze Rezaee J. G. Bosch B. P. F. Lelieveldt, R. J. van der Geest and J. H. C. Reiber. Anatomical model matching with fuzzy implicit surfaces for segmentation of thoracic volume scans. IEEE Transactions on Medical Imaging, 18(3):218–229, March 1999. [5] C. Gui C. Li, C. Xu and M. D. Fox. Level set evolution without reinitialization: A new variational formulation. IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2005. [6] T. F. Chan and L. A. Vese. Active contours without edges. IEEE transactions on image processing, vol. 10, February 2001. [7] Mumford D and J. Shah. Optimal approximations by piecewise smooth functions and variational problems,. Communications on Pure and Applied Mathematics, (5):577–685, 1988.
72
[8] S. van de Vorst M. Prokop E. M. van Rikxoort, B. de Hoop and Bram van Ginneken. Automatic segmentation of pulmonary segments from volumetric chest ct scans. IEEE Transactions on Medical Imaging, 28(4):621–630, April 2009. [9] D. Rueckert et al. Nonrigid registration using free-form deformations: application to breast mr images. IEEE Transactions on Medical Imaging, 18(8):712 – 721, August 1999. [10] J.R. McClelland et al. A continuous 4d motion model from multiple respiratory cycles for use in lung radiotherapy. Medical Physics, Vol. 33(2):3348–3358, August 2006. [11] P. J. Keallb et al. The management of respiratory motion in radiation oncology report of aapm task group 76. Medical Physics, Vol. 33(10)::3874–3900, August 2006. [12] L. Evans. Partial differential equations. American Mathematical Society, 1998. [13] M. H. Hayes. Statistical Digital Signal Processing and Modeling. 1996. [14] M. Prokop I. Sluimer and B. van Ginneken. Toward automated segmentation of the pathological lung in ct. IEEE Transactions on Medical Imaging, 24(8):1025–1038, August 2005. [15] S. Fomin I.Gerfand. Calculus of Variations. 2000. [16] J. Weickert L. Alvarez and J. Snchez. A scale-space approach to nonlocal optical flow calculations. In ScaleSpace 99, 1999. [17] M. Lefebure L. Alvarez, J. Esclarin and J. Snchez. A pde model for computing optical flow. Proc. XVI Congreso de Ecuaciones Difrenciales y Aplicaciones, 1999. [18] Varshney Lav R. Abdominal organ segmentation in ct scan images: A survey. August 2002. 73
[19] T. Lindeberg. Scale-space theory in computer vision. Kluwer, 1994. [20] A. Witkin M. Kass and D. Terzopoulos. Snakes: Active contours model. Int. J. Computer Vision, vol. 1, 1988. [21] J. Nascimento M. Silveira and J. Marques. Automatic segmentation of the lungs using robust level sets. In Proceedings of the 29th Annual International Conference of the IEEE EMBS, Cit Internationale, Lyon, France 2007. [22] D. Mumford and J. Shah. Optimal approximation by piecewise smooth functions and associated variational problems. Commun. Pure Appl. Math, vol. 42, 1989. [23] M. Nixon and Al. Aguado. Feature Extraction and Image Processing (second edition). 2002. [24] O. Colliot O. Camara and I. Bloch. Computational modeling of thoracic and abdominal anatomy using spatial relationships for image segmentation. Real-Time Imaging, Vol. 10(4):263 – 273, 2004. Special issue on imaging in bioinformatics: Part III. [25] J. Alirezaie O. Talakoub and P. Babyn. Lung segmentation in pulmonary ct images using wavelet transform. Acoustics, Speech and Signal Processing, 2007. ICASSP 2007. IEEE International Conference, 2007. [26] S. Osher and J. A. Sethian. Fronts propagating with curvature dependent speed: Algorithms based on hamilton-jacobi formulation. Comput. Physics, vol. 79, 1988. [27] Stanley Osher and R. Fedkiw. Level Set Methods and Dynamic Implicit Surfaces. Springer, 2002. [28] A. Papoulis. Probability, Random Variables and Stochastic Processes. 1991. [29] M. Petrou and P. Bosdogianni. Image Processing: The fundamentals. 1999. 74
[30] William K. Pratt. Digital Image Processing: PIKS Inside. John Wiley and Son, 3rd edition, 2001. [31] J. C. Beatty R. H. Bartels and B. A. Barsky. An Introduction to Splines for Use in Computer Graphics and Geometric Modeling. 1987. [32] James A. Sethian R. Malladi and B. C. Vemuri. Shape modeling with front propagation:a level set approach. IEEE transactions on pattern analysis and machine intelligence, vol. 2, 1995. [33] E. A. Hoffman S. Hu and J. M. Reinhardt. Automatic lung segmentation for accurate quantitation of volumetric x-ray ct images. IEEE Transactions on Medical Imaging, 20(6):490–498, June 2001. [34] J.A. Sethian. Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science. 1999. [35] MacMillan Cancer Support. http://www.cancerbackup.org.uk. [36] R. Kimmel V. Caselles and G. Sapino. On geodesic active contours. Int. J. Comput. Vis., vol. 22, 1997. [37] J. Weszka. A survey of threshold selection techniques. Computer Greaphics and Image Processing, 1978. [38] V. Boldea D. Sarrut L. Brard Z. Wu, E. Rietzel and G. C. Sharp. Evaluation of deformable registration of patient lung 4dct with subanatomical region segmentations. Medical Physics, Vol. 35(2):775–781, January 2008.
75
