立命館学術成果リポジトリ

Doctoral Thesis reviewed

by Ritsumeikan University

R

OBUST AND

A

CCURATE

M

ULTI-

O

RGAN

S

EGMENTATION FOR

M

EDICAL

I

MAGE

A

NALYSIS

( 医用画像解析のための頑健で正確な複数臓器セグ

メンテーション法 )

March 2016

2016 年 3 月

Doctoral Program in Advanced Information Science and Engineering

Graduate School of Information Science and Engineering

Ritsumeikan University

立命館大学大学院情報理工学研究科

情報理工学専攻博士課程後期課程

DONG Chunhua

ドン チュンファ

Supervisor: Professor CHEN Yen-Wei

研究指導教員: 陳 延偉 教授

Abstract

Organ segmentation from CT volumes is an important prerequisite for surgery planning, computer aided surgery, computer assisted intervention and image guid-ed surgery. Accurate segmentation of organs from CT images is considerguid-ed a challenging task: large variations in shape make accurate segmentation difficult and existing lesions (e.g. tumors) exhibit considerable variations for the organ’s anatomical structure. In order to segment the abdominal CT images into differ-ent tissues, various approaches have been proposed. Recdiffer-ently, there are two main trends for organ segmentation: (1) Graphical models based interactive methods, such as graph cuts and random walks; and (2) Anatomical models based meth-ods, such as probabilistic atlas and statistical shape model. In this thesis, our research focuses on improving the classical random walks-based and probabilis-tic atlas-based methods for organ segmentation and their applications. The main contributions of this dissertation are summarized as below:

1. An efficient interactive method based on random walks is proposed for med-ical image segmentation. The conventional random walks (RW) based 3D segmen-tation method faced the following two problems: heavy compusegmen-tational burdens and inaccurate segmentation due to a large-scale graph. In order to overcome the above problems, we use a knowledge-based segmentation framework based on the random walks for volumetric medical image in a slice-by-slice manner. The proposed method employs the previous segmented slice as prior knowledge (the shape and intensity constraints) for automatic segmentation of the adjacent slice, which can reduce the graph scale and significantly speed up the optimization pro-cedure for the graph. With a small number of user-defined seeds, we can obtain the segmentation results of the start slice in the volume, which can be used as prior knowledge of the segmented organ in the following slices. According to this prior knowledge, the object/background seeds can be dynamically updated for the adjacent slice. The average Dice coefficients of the liver and spleen are 0.950 and 0.952, respectively.

2. An automatic algorithm based on an iterative probabilistic atlas is proposed for multi-organ segmentations. A probabilistic atlas based on the human anatom-ical structure has been widely used for organ segmentation in a Bayes framework.

How to register the probabilistic atlas to the patient volume is the main challenge. Additionally, there is the disadvantage that the conventional probabilistic atlas may cause bias toward the specific patient study because of the single reference. As an improvement over the previous registration scheme, we use a template matching strategy based on an iterative probabilistic atlas for fully automatic multi-organ (liver and spleen) segmentation. We firstly detect a bounding box for the tar-get organ based on human anatomical localization. Then, the probabilistic atlas is used as a template to find the organ in this bounding box by using template matching. The average Dice coefficients of the liver and spleen are 0.930 and 0.922, respectively.

3. As an application, a system is develop to assess locoregional therapy (LT) of hepatocellular carcinoma (HCC) by applying our proposed segmentation and registration techniques. In this system, the proposed knowledge-based of the ran-dom walks segmentation method can be used to segment the livers, tumor and treated region from CT image. The proposed anatomical structure constraints based the non-rigid registration method can be used to automatically register the segmented treated region after LT to the segmented tumor before LT. The utility of our system for real clinical use was evaluated by doctors.

Contents

1 Introduction 1

1.1 Medical Imaging Technology . . . 1

1.2 Medical Image Analysis . . . 6

1.3 Medical Image Enhancement . . . 7

1.4 Medical Image Segmentation . . . 9

1.5 Medical Image Registration . . . 10

1.6 Medical Image Visualization . . . 11

1.7 Contributions of the dissertation . . . 12

2 Image Segmentation Methodology 21 2.1 Thresholding Method . . . 21

2.2 Region Growing Method . . . 23

2.3 Watershed Method . . . 24

2.4 Clustering Method . . . 25

2.5 Snakes Method . . . 26

2.6 Level Set Method . . . 26

2.7 Geodesic Active Contour Method . . . 28

2.8 Statistical Shape/Appearance Model Method . . . 29

2.9 Probabilistic Atlas Method . . . 30

2.10 Intelligent Scissors Method . . . 31

2.11 Normalized Cut Method . . . 31

2.12 Graph Cut Method . . . 33

2.13 Random Walks Method . . . 34

3 Efficient Interactive Segmentation Algorithm of Multi-Organ based on Random Walks 41 3.1 Introduction . . . 41

3.2 Random Walkers for Image Segmentation . . . 43

3.3 Efficient and Accurate Approach for Volumetric Medical Image . . . 46

3.4 Extension to Multi-Organ Segmentation . . . 53

3.6 Discussion . . . 65

3.7 Conclusion . . . 69

4 Automatic Segmentation Algorithm of Multi-Organ based on Com-putational Anatomical Models 73 4.1 Introduction . . . 73

4.2 Methods . . . 74

4.3 Results . . . 85

4.4 Discussion . . . 98

4.5 Conclusion . . . 101

5 Evaluation System for Locoregional Therapy of Hepatocellular Carcinoma 107 5.1 Introduction . . . 107 5.2 System Design . . . 108 5.3 Segmentation Module . . . 108 5.4 Registration Module . . . 110 5.5 Visualization Module . . . 115 5.6 Conclusion . . . 117 6 Conclusion 123 Publication List 125 Acknowledgement 129

List of Figures

1.1 Examples of NC CT images, ART CT images, phase-PV CT images and phase-DL CT images. . . 3 1.2 Examples of T1-weighted MR images, T2-weighted MR images and

PD-weighted MR images. . . 5 1.3 The techniques on analyzing medial image. . . 7 1.4 Diagram of the typical intensity-based registration algorithms. . . 10 1.5 The contributions of this dissertation. . . 12 2.1 Segmentation results of one slice from an abdominal CT scan using

the thresholding method. (a) The input image; (b) The segmenta-tion result with the threshold [100, 200]. . . 22 2.2 Segmentation process based on the region growing method. (a)Start

by choosing an arbitrary seed pixel and compare it with neighbour-ing pixels; (b) Region is grown from the seed pixel by addneighbour-ing in neighbouring pixels that are similar, increasing the size of the region. 23 2.3 Segmentation results of one slice from an abdominal CT scan using

various region growing method (threshold [130, 200]). (a) The input image and one seed point; (b) Segmentation result of the connect-ed threshold method; (c) Segmentation result of the neighborhood connected method. . . 24 2.4 Segmentation results of one slice from an abdominal CT scan

us-ing watershed method. (a) The input image; (b) The watershed results with parameters: Threshhold=0.006, Level=0.02; (c) The watershed results with parameters: Threshhold=0.005, Level=0.01. 25 2.5 Illustration of zero set in a level set [26]. . . 27 2.6 Segmentation process based on the shape detection level set method.

(a) The input image and one seed point; (b) Gradient magnitude of the smoothed image; (c) Sigmoid of the gradient magnitude, where the sigmoid is used to compute the speed term for the front propa-gation; (d) The segmentation result. . . 28 2.7 Segmentation process based on the GAC method. . . 29

3.1 Illustration of the random walker for segmentation. Three seed points (S1, S2, S3) belong to three different classes (Class 1, Class 2 and Class 3), alternately fix the probability of each class to unity and the remaining nodes to zero. (a) Initial seed points. (b) Probability of the node belong to Class 1. (c) Probability of the node belong to Class 2. (d) Probability of the node belong to Class 3. (e) Final segmentation resulting from assigning each node the class that corresponds to its greatest probability. The probabilities at each node sum to unity. . . 44 3.2 The basic ideas of our proposed method. . . 47 3.3 The whole procedure of our proposed method. . . 48 3.4 Extraction of the abdominal region from the raw CT data. (a) One

slice of the raw CT data; (b) Intensity variations; (c) The bounding box of the abdomen; (d) Extraction of the abdominal region. . . . 49 3.5 Steps of the RWNBT method. (a) The segmented liver (red) of

the previous slice; (b) The current slice; (c) Candidate pixel by applying a threshold to the GMM; (d) The rough object (red) and background (green) seed points; (e) The fine seed points using the NBT method; (f) The initial segmentation result by RWNBT; (g) Smoothing the boundary by Fourier transform; (h) Comparison of final(red) and groundtruth (blue) liver border; (k) Visualisation of the segmented liver volume. . . 51 3.6 Comparison of the manual segmentation (blue) with the

segmenta-tion results of our method (red). The first row is the segmentasegmenta-tion result in Case 9. The second row is the segmentation result in Case 16. The third row is the segmentation result of pathological case with unusual liver shape in Case 21. . . 57 3.7 Our technique performed on the Group-I dataset with Dice

mea-surement. The first 20 data is normal cases and the remaining 5 data is pathological cases. . . 58 3.8 Comparison of the runtime and accuracy of the classical RW3D

(blue) and our proposed method RWNBT (red) for different sized 3D image. . . 60 3.9 Comparison of liver segmentation results with RWNBT method and

RW3D method in Case 4. . . 60 3.10 Comparison of liver segmentation results with RWNBT method and

3.11 Comparison of manual segmentation (blue) and the segmentation results of our method (red) for the Group II dataset (High-resolution). The first row is Case 1 of Group-II dataset. The second row is Case 3 of Group-II dataset. (a) The seed points on one start slice; (b)-(d) Different slices of the segmented volume; (e) Visualization of

the segmented volume. . . 62

3.12 The segmentation accuracy for the Group II dataset (High-resolution). 63 3.13 Comparison of the segmentation results of our method with the Group III dataset (Low-resolution). The first row is Case 11 of Group-III dataset. The second row is Case 13 of Group-III dataset. (a) The seed points on one start slice; (b)-(d) Different slices of the segmented volume. . . 64

3.14 The liver segmentation accuracy for the Group III dataset (Low-resolution). . . 64

3.15 The spleen segmentation accuracy for the Group III dataset (Low-resolution). . . 65

3.16 Comparison of the segmentation results with RWNBT method and RW3D method in Case 6 of Group-III dataset. . . 66

3.17 Effect of resolution on segmentation accuracy in Case 1 of Group-I dataset. . . 67

3.18 Effect of noise on segmentation accuracy for one 2D slice. (a) The original image and seed points; (b) The tumor segmentation results using RWNBT, Dice=0.974; (c) The segmentation result Dice=0.935 at Gaussian noise with standard deviations sigma = 0.3; (d) Dice similarity coefficient between the segmentations at varying levels of noise. . . 68

4.1 The framework of our proposed segmentation method. . . 75

4.2 Schematic of organ bounding box construction. . . 76

4.3 The flowchart of the bone segmentation method. . . 77

4.4 Intensity distribution of one CT volume (blue line). The yellow line is the assumed Gaussian distribution (G(µ, σ)). . . 78

4.5 (a) Bone segmentation accuracy with different thresholds (P); (b) P = 1%; (c) P = 2%; (d) P = 5%. . . 79

4.6 Steps of automatic extraction of a seed point in the third slice of CT volume. (a) One CT slice; (b) The extracted abdominal region; (c) The central zone of the abdomen; (d) The input image is threshold-ed using the optimum threshold [Tlow, Thigh]; (e) The central point (Green point) of the largest object; (f) The final seed point. . . 79

4.7 Four different reference bones. (a) Reference 1 (shoulder to legs); (b) Reference 2 (Liver Region to legs); (c) Reference 3 (Liver Region to pelvis); (d) Reference 4 (Liver Region). . . 80 4.8 The bounding box for the training data. (a) Examples of bone

segmentation; (b) Registration of bone (Gray-reference bone, Blue-moving bone before registration, Green-the Blue-moving bone after reg-istration); (c) The bounding box of the liver (the ensemble of six registered livers); (d) The bounding box of the spleen (the ensemble of six registered spleens) . . . 81 4.9 Building the spleen (L = 2) likelihood map. . . 82 4.10 The diagram of the iterative atlas construction. . . 83 4.11 The bone segmentation results with various contrast phases. The

first row is the phase-ART images (2 cases); The second row is the phase-PV images (2 cases); The third row is the phase-DL images (2 cases). . . 87 4.12 Effect of the segmented bones on the liver/spleen segmentation

re-sults with the phase-PV CT images. The first row is the bone segmentation result only contains the bone parts; The second row is the bone segmentation result contains the bone parts and trivial non-bone parts. . . 87 4.13 Registration results of four typical data with different bone

ref-erence (Gray-refref-erence bone, Blue-moving bone before registration, Green-moving bone after registration). The first row is the registra-tion results with Reference 1 and the second row is the registraregistra-tion results with our four references ( Reference 1, 2, 3 and 4). . . 89 4.14 The estimated bounding box for the liver and spleen on Case 22. . . 90 4.15 The estimated bounding box for the liver and spleen on Case 5. . . 91 4.16 Comparison of different liver segmentation results on Case 11. (a)

The original CT image; (b) The manual segmentation; (c) Like-lihood map; (d) ConSeg-ConAtlas method; (e) ConSeg-IterAtlas method; (f) OurSeg-IterAtlas method. . . 93 4.17 The ROC of liver with different thresholds on Case 11. . . 93 4.18 The final liver segmentation result performed on two cases. The

first row is Case 11 and the second row is Case 1. . . 93 4.19 Qualitative comparison of different segmentation methods for the

liver. The first 6 data is ART-phase CT cases; the next 12 data is PV-phase CT cases; the remaining 4 data is DL-phase CT cases. . . 94 4.20 The segmentation accuracy can be improved by increasing the

4.21 Comparison of different spleen segmentation results on Case 11. (a) The original CT image; (b) The manual segmentation; (c) Like-lihood map; (d) ConSeg-ConAtlas method; (e) ConSeg-IterAtlas method; (f) OurSeg-IterAtlas method. . . 96 4.22 The ROC of spleen with different thresholds on Case 11. . . 97 4.23 The final spleen segmentation result performed on two cases. The

first row is Case 11 and the second row is Case 1. . . 97 4.24 Qualitative comparison of different segmentation methods for the

spleen. The first 6 data is ART-phase CT cases; the next 12 data is PV-phase CT cases; the remaining 4 data is DL-phase CT cases. 98 4.25 Three examples of segmentations for the pathological liver cases

(the tumor, enlarged livers with unusual liver shapes). . . 99 4.26 Examples from three different patients of segmentations for the

pathological spleen cases (the tumor, enlarged spleens with unusual spleen shapes). . . 100 5.1 The diagram of our system. . . 109 5.2 Steps of the segmentation module. . . 110 5.3 Illustration of landmark positions before LT and after LT. (a) the

branch of the right hepatic (the first one is before LT and last one is after LT); (b) the umbilical portion; (c) the first branch of the portal vein; (d) the second branch of the portal vein. . . 113 5.4 Steps of the registration module. . . 116 5.5 3D visualization of registration results with the same constrained

weights (λ=0.0001) using different landmarks (from 1 to 3 land-marks). (a) The original segmented data (Purple represent the fixed landmark before LT, Green represent the moving landmark after LT); (b) The classical non-rigid registration result; (c) The proposed non-rigid registration result with one landmark; (d) The proposed non-rigid registration result with two landmarks; (e) The proposed non-rigid registration result with three landmarks. . . 117

List of Tables

3.1 Segmentation accuracy obtained by the state-of-the-art methods for the liver on Group-I dataset. . . 61 3.2 Effect of resolution on segmentation accuracy. . . 63 3.3 Segmentation accuracy obtained by the state-of-the-art methods on

the Group-III dataset. . . 66 4.1 Segmentation accuracy obtained by different methods. . . 94 4.2 Evaluation of the MICCAI testing dataset. . . 95 4.3 Comparative statistics p values between different segmentation

meth-ods for the liver and spleen. (Method1:ConSeg-ConAtlas; Method2:ConSeg-IterAtlas; OurMethod:OurSeg-IterAtlas) . . . 100

Chapter 1

Introduction

With the rapid development of medical imaging technologies, high-resolution med-ical images have been achieved to assist scientists and physicians in clinmed-ical diagno-sis and treatment planning. Research efforts have been focused on medical image processing and analysis to extract clinically meaningful information about anatom-ical structures. Medanatom-ical image segmentation and registration are fundamental pro-cesses in most clinical applications that assist medical diagnosis, surgical planning and disease treatments. In this dissertation, we focus on the segmentation and registration of medical images and their applications. Thus, we provide a brief overview of medical imaging technologies used to obtain different medical image modalities and describe basic processing and analysis technologies.

1.1

Medical Imaging Technology

Medical imaging techniques are used to produce an image that reveals the interior of a body for clinical diagnosis and surgical planning. At present, a large number of imaging modalities [1, 2] suitable for clinical application are available. In this chapter, we describe the major types of medical imaging modalities.

1.1.1

X-Ray Radiography

X-ray radiography is the oldest and most widely used medical imaging technologies. X-rays were discovered in 1895 by Wilhelm Conrad Roentgen. X-rays are a type of radiation called electromagnetic waves. By passing X-rays, which are a form of ionizing radiation through the body, we can generate images of internal structures, because different tissues absorb different amounts of radiation. Calcium in bones absorbs X-rays the most; therefore, bones look white. Fat and other soft tissues absorb less radiation and appear gray. X-rays were originally used to image bones,

which were easily distinguishable from soft tissues. At present, X-ray imaging is also used in specialized examinations of areas such as the lungs, breasts, and cavities. The main features of X-ray scans can be summarized as:

• Use an X-ray and observe shadows (ionizing) • Different rates of absorption (tissues)

• Relatively high resolution • Low cost

1.1.2

Computer Tomography (CT)

Computed tomography (CT) [3] is a medical imaging technology that produces detailed cross-sectional images of specific areas inside the body using multiple X-ray projections at different angles. In a CT system, digital image processing is used to generate three-dimensional (3D) images that can represent the internal anatomical structure of an object. CT images can be used to generate 3D images of soft tissues, the pelvis, blood vessels, lungs, the brain, the heart, the abdomen, and bones. CT scans are often the preferred method for detecting many cancers such as liver, lung, and pancreatic cancers. The main features of CT scans can be summarized as:

• Use shaped X-ray beams and observe shadows (ionizing)

• Different radiation absorption of tissues on a series of cross sections • High resolution and anatomical detail

• Intermediate cost

CT is used as a diagnostic tool and a guide for interventional procedures. In many instances, contrast materials, such as intravenous iodinated agents, are painlessly injected into a peripheral vein to generate detailed images showing both blood vessels and tissues. This is useful to differentiate structures, such as blood vessels, from their surroundings. Moreover, contrast materials can help obtain functional information about tissues. Contrast-enhanced CT [4] has recently be-come a popular imaging modality to detect pathologies by enhancing the contrast between a lesion and the normal surrounding structures. The common four en-hancement phases are described as follows:

(1). Precontrast phase (NC). Precontrast scans are acquired before intra-venous contrast administration. Images are used to detect calcification, visualize

Figure 1.1: Examples of phase-NC CT images, phase-ART CT images, phase-PV CT images and phase-DL CT images.

hemorrhage from trauma, and identify lesions by comparing the surrounding nor-mal parenchyma.

(2). Arterial (ART) phase. Scanning in the ART phase is performed nearly 30 seconds after the iodinated contrast injection is inserted into the vein. During this phase, the main portal veins and artery show early enhancement, but the hepatic veins will not yet be enhanced, as shown in the second column of Fig 1.1. (3). Portal venous (PV) phase. Scanning in the ART phase is performed approximately 30 s after the iodinated contrast is injected. During this phase, the main portal veins and artery show early enhancement; however, the hepatic veins will not be enhanced, as shown in the third column of Fig 1.1.

(4). Delayed (DL) phase. Scanning in the DL phase is performed 5-10 mins post contrast. The DL phase is useful to characterize lesions with persistent or late enhancement, such as cavernous hemangiomas. Portal and hepatic veins are somewhat enhanced; however, they are not enhanced as clearly as in the PV phase.

In this dissertation, we focus on CT images that are stored in the Digital Imaging and Communications in Medicine (DICOM) image format. In addition to pixel data, DICOM files contain medical information in the file header, such as patient ID, name, etc.

1.1.3

Magnetic Resonance Imaging (MRI)

Magnetic resonance imaging (MRI) [5] creates 3D cross-sectional images of organs and internal structures using strong magnetic fields and radio waves. An MRI scan

can be used to differentiate normal and diseased soft tissues. The main features of MRI scans can be summarized as:

• Use magnetic fields and detects radiation from an organ (non-ionizing) • Different relaxation time or nuclei density

• Relatively high resolution and differentiation of tissues with functional dif-ferences

• Intermediate cost

Compared to an X-ray, ultrasound scan, or CT scan, an MRI scan can reveal abnormal tissue (such as tumors, bleeding, injury, and blood vessel diseases) more clearly with a contrast material. MRI is the investigative tool of choice for neu-rological cancers, as it has better resolution and visualization than CT. In the MRI system, repetition time (TR) and echo time (TE) are basic pulse sequence parameters. These control times determine the tissue contrast in MRI scans, i.e., T1-weighted, PD-weighted (Proton or spin density), and T2-weighted images, as shown in Fig 1.2. Short TR and TE will generate a T1-weighted image, long TR and short TE will generate a PD-weighted image, and long TR and TE will generate a T2-weighted image.

1.1.4

Positron Emission Tomography (PET)

Positron emission tomography (PET) is a nuclear medicine technique. The detailed images of specific areas inside the body generated by PET are produced using a scanner and a small amount of radiotracer (radioactive substance), which is injected into the patient’s vein. PET scans can provide physicians with information about how tissues and organs are functioning. PET is often used to evaluate neurological diseases. The main features of PET scans can be summarized as:

• Administer isotope and observe the radiation (ionizing) • Track the radiation intensity and location

• Relatively low resolution and functional imaging • Low cost

In addition to patient feedback, a PET scan can potentially help identify disease in its earliest stages. Compared to CT or MRI scans, PET scans are more sensitive, i.e., PET scans can identify changes at a cellular level. PET may detect the early

Figure 1.2: Examples of T1-weighted MR images, T2-weighted MR images and PD-weighted MR images.

onset of disease before it is evident in CT or MRI scans. Combined PET and CT (PETCT) scans have emerged as a new medical imaging technique. This technique acquires images from both devices sequentially and then combines image pairs into a single superposed image. PETCT is a primary tool for the delineation of tumor volumes, staging, and the preparation of patient treatment plans. The combination has been shown to improve oncologic care by positively affecting active treatment decisions, disease recurrence monitoring, and patient outcomes, such as disease-free progression.

1.1.5

Ultrasound Imaging

Diagnostic ultrasound, also known as medical sonography or ultrasonography, can construct images of the internal structures of the body using high frequency sound waves. An ultrasound machine sends sound waves into the body and converts the reflected sound echoes into a picture on the horizontal axis. Ultrasonic images (sonography) are effective for imaging pregnancy, symptoms of pain, swelling, and infection. The main features of Ultrasound scans can be summarized as:

• Different sound-wave propagation velocity and tissue density • Poor resolution and high noise

• Low cost

Compared to the other medical imaging technologies that do not use ionizing radiation, ultrasound examinations are extremely safe. Ultrasound provides real-time imaging of the structure and movement of internal organs.

1.1.6

Other Types of Medical Imaging

Other types of medical imaging techniques, such as elastography, thermography, biomarker imaging, functional MRI (fMRI), and single-photon emission computed tomography, have specific uses in the diagnosis of disease.

1.2

Medical Image Analysis

Medical imaging has been a prominent role in the diagnosis and treatment of disease. Research efforts have been devoted to medical image processing and ana-lyzing for extracting the clinically meaningful information about anatomic struc-tures. The fundamental technologies in analyzing medical image [6, 7] are shown in Fig 1.3, with special emphasis on medical image enhancement, medical image segmentation , medical image registration and medical image visualization.

Medical imaging has played a prominent role in the diagnosis and treatment of disease. Research efforts have been devoted to medical image processing and analysis to extract clinically meaningful information about anatomic structures. The fundamental technologies used to analyze a medical image [6, 7] are shown in Fig 1.3, with special emphasis on medical image enhancement, segmentation, registration, and visualization.

(1). Medical image enhancement: reduces high noise and increases low contrast to enhance the quality of medical images.

(2). Medical image segmentation: separates a medical image into specific objects, such as the liver, spleen, and heart.

(3). Medical image registration: registers two medical images by finding an optimal geometric transformation.

(4). Medical image visualization: displays quantitative information about anatomic tissues and functions that are affected by disease, including surface and volume renderings.

Figure 1.3: The techniques on analyzing medial image.

In our research, we use and extend the enhancement, segmentation, and regis-tration algorithms provided by the Insight Toolkit (ITK) and use the Visualization Toolkit (VTK) to display anatomic tissues.

The ITK is an open-source, cross-platform system that provides developers with an extensive suite of software tools for medical image analysis. The ITK employs numerous leading-edge algorithms to enhance, register, and segment mul-tidimensional data. The VTK is an open-source, freely available software system for 3D computer graphics, image processing, and visualization. The VTK supports a wide variety of visualization algorithms. Developers can use, debug, maintain, and extend the software because ITK and VTK are open-source projects.

1.3

Medical Image Enhancement

With the rapid development of advanced medical equipment, surgeons desire high-resolution medical images to assist diagnosis and interpretation. However, data acquisition devices and illumination conditions often degrade the quality of medical images; thus, development of methods to improve quality is a major focus for medical image processing. The primary purpose of medical image enhancement is to solve the problems associated with low contrast or high levels of noise in medical images. Various medical image enhancement technologies, which mainly concentrate on the following filters, have been proposed.

(1). Gradient Magnitude Filter. The magnitude of the image gradient is used extensively in medical image enhancement, primarily to help determine object contours and separate homogeneous regions. ITK filters attempt to reduce this ambiguity by including the magnitude term when appropriate. The ITK pro-vides itk::GradientMagnitudeImageFilter to compute the magnitude of the image gradient at each pixel location using a simple finite differences approach.

finding the statistical mean of the neighborhood of the corresponding input pixel. Note that this algorithm is sensitive to outliers in the neighborhood. The ITK providesitk::MeanImageFilter to remove noise from images.

(3). Median Filter. This filter computes the value of each output pixel as the statistical median of the neighborhood around the corresponding input pixel. This filter is particularly efficient against salt-and-pepper noise. In other words, it is robust to gray level outliers. The ITK provides itk::MedianImageFilter to reduce the strength of noise at each pixel.

(4). Gaussian Filter. Blurring is the traditional approach for denoising im-ages. The Gaussian filter is the most commonly used blurring method. This filter computes the convolution of the input image using a Gaussian kernel. Large Gaus-sian variances will produce large convolution kernels and require correspondingly longer computation times. The ITK provides itk::DiscreteGaussianImageFilterto smooth the image by implementing a convolution with a kernel.

(5). Gradient Anisotropic Diffusion Filter. The drawback of image denoising (smoothing) is that it tends to blur sharp boundaries in the image. Anisotropic diffusion attempts to reduce noise and texture while preserving (en-hancing) edges. This filter implements an N-dimensional version of the classic PeronaMalik anisotropic diffusion equation [11]. The ITK provides itk::Gradient AnisotropicDiffusionImageFilterto smooth homogeneous regions and preserve edges.

(6). Curvature Anisotropic Diffusion Filter. Theitk::Curvature Anisotrop-ic DiffusionImageFilterperforms anisotropic diffusion on an image using a modified curvature diffusion equation (MCDE) [12]. MCDE does not exhibit the edge en-hancing properties of classic anisotropic diffusion. It is less sensitive to contrast than classic Perona-Malik style diffusion and preserves finer detailed structures in images.

(7). Fourier Transform Filter. One of the most common image processing operations performed in the Fourier domain is masking the spectrum to eliminate a range of spatial frequencies from the input image. The ITK example FFTIm-ageFilterFourierDomainFiltering.cxx is typically performed by taking the input image, computing its Fourier transform using a Fast Fourier Transform (FFT) fil-ter, masking the resulting image in the Fourier domain with a mask and finally taking the result of the masking and computing its inverse Fourier transform.

These enhancement operations are performed to remove noise, modify image brightness/contrast, or the distribution of grey levels. Image enhancement is com-monly a preprocessing step for image segmentation and registration.

1.4

Medical Image Segmentation

Image segmentation [13] separates an image into several ”homogeneous” regions according to pixel attributes. In the medical domain, elements with similar pixel attributes belong to the same tissue type or organ. Accurate segmentation of anatomical structures is a key enabling technology for medical applications, such as diagnostics, planning, and guidance. Since most medical image studies rely on large amounts of data, manual segmentation of organs by a medical specialist is time consuming. To segment organs accurately, various sophisticated segmentation methods [14] have been proposed. Such methods can generally be classified as automatic and interactive (semi-supervised) methods.

Automatic segmentation methods [15] rely on gradient/intensity analysis or prior knowledge (e.g., location, shape, and appearance) to detect specific objects in images. There is a great deal of automatic image segmentation literature. Well-known algorithms include thresholding [16, 17, 18], region growing [19, 20, 21], clustering-based [22, 23, 24], watershed [25, 26, 27], active contour [28, 29], and level set [30, 31] methods. However, these methods are often susceptible to imaging artifacts, ambiguous intensities, and low contrast, low signal-to-noise ratio imag-ing modalities because no prior knowledge of shape or pose is assumed. As image segmentation problems can easily be translated into optimization problems, some algorithms based on minimizing cost functions have been proposed. The most widely used criteria include expectation maximum (EM) [32, 33], Markov ran-dom field (MRF)[34, 35], and maximum a posteriori estimation (MAP) [36, 37], in which cost functions are minimized by different optimization strategies. Modeling the problem as an energy minimization process, these methods, such as a proba-bilistic atlas [38, 39, 40], statistical shape model [41, 42], and active appearance model [43, 44], aim to produce desirable segmentation by achieving global opti-mization. Anatomical shape priors are very useful for medical image segmentation because they handle shadows and weak edges effectively, and they have been used successfully for different anatomical structures segmentation, such as heart, liver, breast, and kidney. These priors can be encoded as statistical shape models, usual-ly derived from large data set segmentations; thus, they are often computationalusual-ly expensive. Unfortunately, they are of limited use for pathological cases due to high variation in structure, shape, size, and localization of lesions. Despite this drawback, they are used extensively due to their robust and reliable performance in clinical applications.

Interactive segmentation methods, also called semi-supervised methods, require the user to mark some initial seed points or the boundary edge of the region of interest manually. Interactive segmentation has become an effective and popular technique in recent years because the combination of human expertise and ma-chine intelligence can improve the segmentation performance with minimal user

intervention. Numerous interactive segmentation methods based on intelligent s-cissors [45, 46], normalized cut [47, 48], graph cut [49, 50], and random walks (RW) [51, 52, 53] have been used in clinical applications. The main advantage of these segmentation methods is to achieve high accuracy segmentation results without requiring prior knowledge from a large amount of labeled data. This tech-nique utilizes expert knowledge to produce accurate segmentation of anatomical structures, which benefits the measurement and diagnosis of diseases.

Our research focuses on the probabilistic atlas (anatomical model) based auto-matic segmentation method and random walks (graphical model) based interactive segmentation method.

1.5

Medical Image Registration

Image registration [54, 55] is used to find the optimal transformation to transfor-m structures of interest in the satransfor-me coordinate systetransfor-m. A great deal of research has been devoted to image registration. Applications of image registration in the medical field include (1) fusion of anatomical images with functional images, e.g., registering CT/MRI images in PET/fMRI images; (2) intervention and treatment planning; (3) computer-aided diagnosis and follow-up; (4) atlas construction; (5) assisted/guided surgery; (6) anatomy segmentation; and (7) computational mod-eling.

Fig 1.4 illustrates intensity-based image registration methodologies. The main idea is to find the optimal geometric transformation by an optimization process in which a similarity measure is minimized or maximized (also known as the cost function). The similarity measure is used to calculate distinctive intensity in the same region of the input images. The optimizer defines the search strategy, which is a function. An interpolator resamples voxel intensity into the new coordinate system based on the geometric transformation.

The choice of the geometric transformation is crucial to a registration algorith-m. According to the types of geometric transformation, the registration method can be classified as rigid and non-rigid registration. Rigid transformation provides rigid registration; otherwise, non-rigid registration is obtained. Rigid registration describes the global motion of a medical image. The simplest choice is a rigid transformation with six degrees of freedom (parameters), three rotation angles, and three translation vectors on each axis. Rigid registration is used when the subject can be considered a rigid object, e.g., human bone. Non-rigid registration can capture local motion. Non-rigid transformations have more degrees of freedom, such as B-spline, thin-plate spline, elastic, and fluid transformations. Non-rigid registration is a powerful tool for deforming objects, e.g., the liver, spleen, or heart.

1.6

Medical Image Visualization

Medical image visualization [56, 57] is the most fundamental process to diagnose a wide range of medical diseases. Using visualization techniques, we can characterize the structure of organs and tissues in medical data to obtain accurate measure-ment of changes in organ volume. Recently, the most common algorithms for 3D visualization of medical data are surface and volume rendering.

In surface rendering, an organ represented by medical data can be visualized by its surface estimated from the image. A surface fitting technique, such as marching cubes, can be used to convert the volume data into sets of polygons that represent the anatomical surface of interest. A linear approximation leads to triangular surfaces for visualization. This algorithm is applied to place surface patches at each contour point with hidden surface removal and shading to visualize the surface. In volume rendering, rather than visualizing its surface, the tissues in every voxel are visualized through calculations based on the opacities and colors assigned to different tissues. A transfer function, based on optical behavior, maps each intensity value to a color and opacity, ranging from fully opaque to fully transparent. Surface rendering is achieved by a preprocessing step that extracts the surface. Volume rendering is achieved through a classification process in which opacities and colors are assigned for each voxel based on the original intensity values.

Through visualization, we can characterize the quality of medical images. For 3D visualization of medical objects in this dissertation, we apply the VTK to display anatomic tissues.

1.7

Contributions of the dissertation

The contributions of this dissertation primarily include the three aspects summa-rized in Fig 1.5.

Figure 1.5: The contributions of this dissertation.

First, we propose an efficient interactive random walks-based algorithm for multi-organ segmentation. This algorithm is described in Chapter 3. The conven-tional RW-based 3D segmentation method imposes a heavy computaconven-tional burden and yields inaccurate segmentation due to a large-scale graph. To overcome the above problems, we use a knowledge-based segmentation framework for an RW-based volumetric medical image in a slice-by-slice manner. The proposed method employs the previous segmented slice as prior knowledge (i.e., shape and intensity constraints) for automatic segmentation of the adjacent slice, which can reduce the graph scale and significantly speed up the optimization procedure for the graph. With a small number of user-defined seeds, a high-quality image segmentation re-sult can be achieved automatically by combining a narrow band threshold method and a combinational RW algorithm. Experimental results show that the average DICE similarity measure of the liver and spleen are 0.950 and 0.952, respectively. The segmentation results of the proposed method were more accurate than state of the art interactive segmentation methods (p < 0.001).

Second, we propose an automatic iterative probabilistic atlas-based algorithm for multi-organ segmentation. This algorithm is described in Chapter 4. Proba-bilistic atlas based on human anatomical structure has been widely used for organ

segmentation in a Bayes framework. The method to register the probabilistic at-las to the patient volume is the main challenge. In addition, the conventional probabilistic atlas may cause bias in a specific patient study due to the single ref-erence. To improve the previous registration scheme, we use a template matching strategy based on an iterative probabilistic atlas for fully automatic multi-organ (liver and spleen) segmentation. We first detect a bounding box for the target organ based on human anatomical localization. Then, the probabilistic atlas is used as a template to find the organ in this bounding box using template match-ing. The average DICE coefficients of the liver and spleen are 0.930 and 0.922, respectively. In addition, our results show improved segmentation accuracy for multiple organs (p < 0.00001) compared to conventional and recently developed atlas-based methods.

Finally, as an application, we develop a system to evaluate locoregional therapy (LT) of hepatocellular carcinoma. This system is described in Chapter 5. In this system, the proposed knowledge-based RW segmentation method is applied to segment the livers, tumors, and treated regions from CT images. Meanwhile, the proposed non-rigid registration method with anatomical structure constraints is applied to automatically register the treated region after LT to the tumor before LT. Medical specialists have assessed our system and indicated that this work could be used in clinical applications in future.

Bibliography

[1] H. Ammari: Biomedical Imaging Modalities. An Introduction to Mathemat-ics of Emerging Biomedical Imaging, (2008)(Chapter 1) 3-13, Springer-Verlag, Berlin.

[2] M.E. Clouse: Current diagnostic imaging modalities of the liver. Surgical Clin-ics of North America, 69(2)(1989) 193-234.

[3] D.J. Brenner, E.J. Hall: Computed Tomography-An Increasing Source of Ra-diation Exposure. New England Journal of Medicine, 357(22)(2007) 2277-2284. [4] H.M. Cynthia, N.P. Andrew, B. Natalie, K. James, L.F. Yu, C. Jodie: Strate-gies for Reducing Radiation Dose in CT. Radiologic Clinics of North America, 47(1)(2009) 27-40.

[5] W. Robert, Brown, Y.-C.N. Cheng, E.M. Haacke, M.R. Thompson, R. Venkate-san: Magnetic Resonance Imaging: Physical Principles and Sequence Design, 1999, Wiley-Lis, New York.

[6] J.L. Prince, J.M. Links: Medical imaging signals and system. Pearson Educa-tion, 2006.

[7] A. Macovski: Medical imaging systems. Prentice-Hall, 1983.

[8] X. Zhou, T. Kitagawa, T. Hara, H. Fujita, X. Zhang, R. Yokoyama, H. Kondo, M. Kanematsu, H. Hoshi: Constructing a Probabilistic Model for Automated Liver Region Segmentation Using Non-contrast X-Ray Torso CT images. Proc. IEEE International Conference on International Conference for Medical Image Computing and Computer-Assisted Intervention (MICCAI 2006), (2006) 856-863.

[9] Z.G. Pan, J.F. Lu: A Bayes-Based Region-Growing Algorithm for Medical Im-age Segmentation. Journal of Computing in Science and Engineering, 9(4)(2004) 32-38.

[10] H. Park, P.H. Bland, C.R. Meyer: Construction of an abdominal probabilis-tic atlas and its application in segmentation. IEEE Transactions on Medical Imaging, 22(2003) 483-492.

[11] P. Perona, J. Malik: Scale-space and edge detection using anisotropic dif-fusion. IEEE Transactions on Pattern Analysis Machine Intelligence, 12(1990) 629639.

[12] R. Whitaker, X. Xue: Variable-Conductance, Level-Set Curvature for Image Denoising. Proc. International Conference on Image Processing (ICIP 2001), 3(2001) 142-145.

[13] N.R. Pal, S.H. Pal: A review on image segmentation techniques. Journal of Pattern Recognition, 26 (1993) 1277-1294.

[14] D.L. Pham, C. Xu, J.L. Prince: Current methods in medical image segmen-tation. Annual Review of Biomedical Engineering, 2(2000) 315-337.

[15] N. Sharma, L.M. Aggarwal: Automated medical image segmentation tech-niques. Journal of Medical Physics, 35(1)(2010) 3-14.

[16] Y.S. Frank, S.X. Cheng: Automatic seeded region growing for color image segmentation. Image and vision computing Journal, 23(10) 877-886.

[17] N. Otsu: A threshold selection method from grey level histogram. IEEE Trans-actions on Systems, Man and Cybernetics, 9(1)(1979) 62-66.

[18] M. Sezgin, B. Sankur: Survey over image thresholding techniques and quanti-tative performance evaluation. Journal of Electronic Imaging, 13(1)(2004) 146-165.

[19] R. Adams, L. Bischof: Seeded region growing. IEEE Transaction on Pattern Analysis and Machine Intelligence, 16(6)(1994) 641-647.

[20] S. Kamdi, R.K. Krishna: Image Segmentation and Region Growing Algorith-m. International Journal of Computer Technology and Electronics Engineering (IJCTEE), 2(1)(2012) 103-107.

[21] R. Pohle, K. Toennies: Segmentation of medical images using adaptive region growing. Proc. SPIE on Medical Imaging 2001: Image Processing, 1337, (2001). [22] N. Sharma, A.K. Ray, S. Sharma, K.K. Shukla, S. Pradhan, L.M. Aggarwal: Segmentation of medical images using simulated annealing based fuzzy C Mean-s algorithm. International Journal of Biomedical Engineering and Technology, 2(2009) 260-278.

[23] M.N. Ahmed, S.M. Yamany: A modified fuzzy C-means algorithm for bias field estimation and segmentation of MRI data. IEEE Transactions on Medical Imaging, 21(2002) 193-199.

[24] A.H. Foruzan, Y.W. Chen, R.A. Zoroofi, A. Furukawa, Y. Sato, M. Hori, N. Momiyama: Segmentation of Liver in Low-contrast Images Using K-Means Clustering and Geodesic Active Contour Algorithms. IEICE Transactions on Information and Systems, E96-D(2013) 798-807.

[25] V. Grau, A.U. Mewes, M. Alcaniz, R. Kikinis, S.K. Warfield: Improved water-shed transform for medical image segmentation using prior information. IEEE Transactions on Medical Imaging. 23(4) (2004) 447-458.

[26] J. B. Roerdink, A. Meijster: The Watershed Transform: Definitions, Algo-rithms and Parallelization Strategies. Fundamenta Informaticae, 41(2001) 187-228.

[27] S. Wegner, T. Harms, H. Oswald, E. Fleck: The watershed transformation on graphs for the segmentation of CT images. Proc. 13th International Conference on Pattern Recognition (ICPR 1996), (1996) 498-502.

[28] V. Caselles, R. Kimmel, G. Sapiro: Geodesic active contours. International Journal on Computer Vision, 22(1) (1997) 61-97.

[29] P. Karasev, I. Kolesov, K. Fritscher, P. Vela, P. Mitchell, A. Tannenbaum: Interactive Medical Image Segmentation Using PDE Control of Active Contours. IEEE Transactions on Medical Imaging, 32(11)(2013) 2127-2139.

[30] R. Malladi, J. A. Sethian, B. C. Vemuri: Shape modeling with front propaga-tion: A level set approach. IEEE Transactions on Pattern Analysis and Machine Intelligence, 17(2)(1995) 158-175.

[31] M. Droske, B. Meyer, M. Rumpf, K. Schaller: An adaptive level set method for medical image segmentation. Proc. of the Annual Symposium on Information Processing in Medical Imaging (IPMI 2001), (2001) 416-422.

[32] A.P. Dempster, N.M. Laird, D.B. Rubin: Maximum likelihood from incom-plete data via the EM algorithm. Journal of the Royal Statistical Society: Series B, 39(1)(1997) 1-38.

[33] O. Ronen, J.R. Rohlicek, M. Ostendorf: Parameter estimation of dependence tree models using the EM algorithm. IEEE Signal Processing Letters, 2(8)(1995) 157-159.

[34] H. Derin, H. Elliott, R. Cristi, D. Geman: Bayes smoothing algorithms for segmentation of binary images modeled by Markov random fields. IEEE Trans-actions on Pattern Analysis and Machine Intelligence, 6(6)(1984) 707-719. [35] C. A. Bouman, M. Shapiro: A multiscale random field model for Bayesian

image segmentation. IEEE Transactions on Image Processing, 3(2)(1994) 162-177.

[36] S.F. Chen, L.L. Cao, Y.M. Wang, J.Z. Liu, X.O. Tang. Image Segmentation by MAP-ML Estimations, IEEE Transactions on Processing, 19(9)(2010) 2254-2264.

[37] H. Park, P.H. Bland, C.R. Meyer: Construction of Abdominal Probabilistic Atlases and Their Value in Segmentation of Normal Organs in Abdominal CT Scans. IEICE Transactions on Information and Systems, E93-D(2010) 2291-2301.

[38] X. Zhou, T. Kitagawa, K. Okuo, T. Hara, H. Fujita, R. Yokoyama, M. Kane-matsu, H. Hoshi: Construction of a probabilistic atlas for automated liver seg-mentation in non-contrast torso CT images. Proc. International Congress Series, 1281(2005) 1169-1174.

[39] X. Zhou, T. Kitagawa, T. Hara, H. Fujita, X. Zhang, R. Yokoyama, H. Kondo, M. Kanematsu, H. Hoshi: Constructing a Probabilistic Model for Automated Liver Region Segmentation Using Non-contrast X-Ray Torso CT images. Proc. IEEE International Conference on International Conference for Medical Image Computing and Computer-Assisted Intervention (MICCAI 2006), (2006) 856-863.

[40] M.G. Linguraru, J.K. Sandberg, Z. Li, F. Shah, R.M. Summers: Automated segmentation and quantification of liver and spleen from CT images using nor-malized probabilistic atlases and enhancement estimation. International Journal of Medical Physics, 37(2)(2010) 771-783.

[41] T. Okada, R. Shimada, M. Hori, M. Nakamoto, Y.-W. Chen, H. Nakamau-ra, Y. Sato: Automated segmentation of the liver from 3D CT images using probabilistic atlas and multi-level statistical shape model. Journal of Academic Radiology, 63(2008) 1390-1403.

[42] M. Lee, W. Cho, S. Kim, S. Park, J.H. Kim: Segmentation of interest region in medical volume images using geometric deformable model. Journal of Computers in Biology and Medicine, 42(5)(2012) 523-537.

[43] T.F. Cootes, A. Hill, C.J. Taylor, J.Haslam: The use of active shape mod-els for locating structures in medical images. Image and Vision Computing, 12(6)(1994) 355-366.

[44] T.F. Cootes, G.J. Edwards, C.J. Taylor: Active appearance models. Proc. European Conference on Computer Vision (ECCV 1998), 2(1998) 484-498. [45] A. Mishra, A. Wong, W. Zhang, D. Clausi, P. Fieguth: Improved interactive

medical image segmentation using Enhanced Intelligent Scissors (EIS). Proc. IEEE Conference on Engineering in Medicine and Biology Society (EMBC2008), (2008) 3083-3086.

[46] A. Blake, C. Rother, M. Brown, P. Perez, P. Torr: Interactive Image Seg-mentation using an adaptive GMMRF model. Proc. 8th European Conference in Computer Vision (ECCV 2004), (2004) 1-14.

[47] J. Shi, J. Malik: Normalized cuts and image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(8)(2000): 888-905.

[48] S.X. Yu, J. Shi: Segmentation given partial grouping constraints. IEEE Trans-actions on Pattern Analysis and Machine Intelligence, 26(2) 173183, 2004. [49] Y. Boykov: Graph Cuts and Efficient N-D Image Segmentation. International

Journal of Computer Vision, 70(2)(2006) 109-131.

[50] Y. Boykov, M.P. Jolly: Interactive graph cuts for optimal boundary and region segmentation of objects in N-D images. Proc. IEEE International Conference on Computer Vision (ICCV 2001), 1(2001) 105-112.

[51] L. Grady: Random Walks for Image Segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 28(11) (2006) 1768-1783.

[52] L. Grady, T. Schiwietz, S. Aharon, R. Westermann: Random walks for in-teractive organ segmentation in two and three dimensions: Implementation and validation. Proc. IEEE International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI 2005), 3750(2005) 773-780. [53] L. Grady, G. Funka-Lea: Multi-Label Image Segmentation for Medical

Ap-plications Based on Graph-Theoretic Electrical Potentials. Proc. IEEE Interna-tional Conference on Computer Vision and Mathematical Methods in Medical and Biomedical Image Analysis (ECCV 2004), 3117(2004) 230-245.

[54] F.P. Oliveira, J.M. Tavares: Medical image registration: a review. Computer Methods in Biomechanics and Biomedical Engineering, 17(2)(2014) 73-93.

[55] V.R.S Mani, Dr.S. rivazhagan: Survey of Medical Image Registration. Journal of Biomedical Engineering and Technology, 1(2)(2013) 8-25.

[56] K.U. Jayaram, H.M. Hung, K.S. Chuang:Surface and Volume Rendering in Three-Dimensional Imaging: A Comparison. Journal of Digital lmaging, 4(3)(1991) 159-168.

[57] R. Shahidi:Surface rendering versus volume rendering in medical imaging: techniques and applications. Proc. the 7th conference on Visualization (VIS’96), (1996) 439-440.

Chapter 2

Image Segmentation Methodology

Organ segmentation from medical images plays a prominent role in analyzing and processing medical images. Organ segmentation identifies the boundaries of ob-jects, such as organs or abnormal regions, in images. The extracted anatomical structure makes it possible to analyze organ shape, detect volume changes, and determine precise intervention and treatment plans. In this dissertation, we focus on organ segmentation in CT images.

Image segmentation means to separate an image into several ”homogeneous” regions according to pixel attributes (intensity, texture or color). If we treat each pixel as a sample and use gray level as the feature, image segmentation can be treated as a pixel classification task. The pixel classification task can be defined as:

Cl = arg max

c

(pc(x)) (2.1)

where c = 1, 2, 3, . . . , C and C is the number of the class. x is the feature vector for the pixel. pc(x) is the probability of the pixel belong to the c−th class.

Finally, for each pixel, the label of the pixel belongs to the Cl-th class by comparing the probabilities. Hence, the image segmentation problem can be regarded as the problem how to calculate the probability of the pixel. We review the common and most important techniques available for medical image segmentation in the following.

2.1

Thresholding Method

Thresholding [2, 3, 4] technology is one of the most basic intensity-based segmen-tation methods. The simplest approach identifies objects using a single threshold value T , such that image intensity values (I (x, y)) above and below the threshold

are object pixels (set to value 1) and background pixels (set to value 0) respectively. The segmentation process using a single threshold value is defined as follows:

B (x, y) = 

1 I (x, y) ≥ T

0 otherwise (2.2)

In addition, the user can also define two thresholds: Lower and Upper. For each voxel, if the voxel value belongs to the intensity range [Lower, U pper], it is assigned the object; otherwise, it is assigned to the background, as shown in the Equation 2.3. The ITK source code for this method can be implemented by BinaryThresholdImageFilter.cxx. In this code, the SetLowerThreshold() and SetUpperThreshold() methods define the threshold range of the input image inten-sities. Fig 2.1 illustrates the segmentation results of this method on an abdominal CT image.

B (x, y) = 1 Tlower ≤ I (x, y) ≤ Tupper

0 otherwise (2.3)

Figure 2.1: Segmentation results of one slice from an abdominal CT scan using the thresholding method. (a) The input image; (b) The segmentation result with the threshold [100, 200].

In the thresholding method, the segmentation problem can be viewed as a selection process for the appropriate threshold value. One common way to select the threshold is to analyze the peaks and valleys of the image histogram [5, 6]. The ideal case is that the histogram presents only two dominant modes. In this case, the value of T is selected as the valley point between the two dominant modes. However, in medical images, the intensity distribution is typically not bimodal. This may lead to inaccurate partitioning of the images into various anatomical structures. On the other hand, the thresholding technique does not utilize the spatial information of images. Therefore, this method can only be applied to well-circumscribed objects. When applied to medical images, the thresholding technique is only used as a preprocessing operation. As a result, expert intervention is required to achieve a refined result.

2.2

Region Growing Method

Region-based segmentation algorithms are recursive operations that connect neigh-boring pixels with similar features (gray value or texture). Region growing [7, 8, 9, 10] is a simple region-based image segmentation method, in which regions are connected based on a predefined criterion and the initial seed points. The se-lection process of the seed points can be performed automatically or manually. Each seed point is considered an initial region. Then, the regions are grown from these seed points to adjacent points with similar image features. The criterion of similarity/homogeneity is crucial for the segmentation result. For example, an appropriate range of similarity threshold value helps improve segmentation results. The whole procedure is iterated until all pixels are partitioned into special regions. The region growing approach is opposite to a split and merge approach. Region growing begins by choosing an initial seed pixel and comparing it to neighboring pixels (Fig 2.2).

Figure 2.2: Segmentation process based on the region growing method. (a)Start by choosing an arbitrary seed pixel and compare it with neighbouring pixels; (b) Region is grown from the seed pixel by adding in neighbouring pixels that are similar, increasing the size of the region.

Several implementations of region growing are available in ITK, such as Con-nectedThresholdImageFilter.cxx, OtsuThresholdImageFilter.cxx, Neighborhood-ConnectedImageFilter.cxx and IsolatedConnectedImageFilter.cxx. Most source code must provide the seed points and user-specified lower/upper thresholds. Fig 2.3 illustrates the segmentation results of the connected threshold method in an abdominal CT image.

Since the regions are grown primarily based on the image’s intensity or edge in-formation, the region growing method is sensitive to noise. It can achieve expected segmentation if the edges of the images are well presented. However, for medical data, images usually contain significant amounts of noise or intensity variation, which may result in holes or over-segmentation.

Figure 2.3: Segmentation results of one slice from an abdominal CT scan using various region growing method (threshold [130, 200]). (a) The input image and one seed point; (b) Segmentation result of the connected threshold method; (c) Segmentation result of the neighborhood connected method.

2.3

Watershed Method

The watershed algorithm [11, 12, 13] is a special type of region-based segmentation. Rather than using image intensity or edge information, the watershed approach assigns pixels to corresponding regions using gradient information on image fea-tures and analyzing weak points along region boundaries. The intuitive idea of this work originates from a concept taken from geology; rain over a landscape or topographic relief flows into lower basins due to gravity. With increased rainfall, the water in catchment basins (regions) increases gradually. Then, the water may spill into another basin. Finally, small basins can be merged to form a larger basin. In an image, pixel intensity is considered the height of the landscape. The progress of the water with overflow continues until the water level reaches the highest peak in the landscape. Thus, regions are established by associating pixels in the image based on local geometric structures, such as gradient or curvature magnitude. The ITK example WatershedSegmentation2.cxx illustrates how to preprocess and segment images using the itk::WatershedImageFilter. There are two parameters: SetLevel() controls watershed depth and SetThreshold() controls the lower thresholding of the input. Fig 2.4 shows the output of an abdominal CT image with different labels, where a label denotes membership of a pixel in a particular segmented region.

Compared to classical region-growing methods, the watershed technique is less sensitive to user-defined thresholds. However, similar to the drawback of the region-growing method, the watershed approach is also sensitive to noise. For structures with a low signal-to-noise ratio, the segmentation result may not be desirable. Moreover, it can result in over-segmentation due to local minimums, which, in practice, give too many regions.

Figure 2.4: Segmentation results of one slice from an abdominal CT scan using watershed method. (a) The input image; (b) The watershed results with parame-ters: Threshhold=0.006, Level=0.02; (c) The watershed results with parameparame-ters: Threshhold=0.005, Level=0.01.

2.4

Clustering Method

Most clustering-based segmentation algorithms [14] do not depend on training data. This method uses an iterative alternation between image segmentation and characterization of the properties of each class. Image segmentation can be treated as a clustering process. First, the pixel calculates the texture feature vector based on local neighborhoods. Then, each pixel is classified into the specific region based on the corresponding texture feature vector. K-means [15, 16, 17], fuzzy C-means [18, 19, 20], and EM [21, 22] are three typical clustering approaches. The K-means algorithm produces a hard segmentation result. It iterates computing the mean intensity for each class, and then the pixel is assigned to a class based on the closest mean. The fuzzy C-means algorithm produces a soft segmentation result. It allows the pixel to have cluster membership, and then the pixel is assigned to a class based on the maximum membership value. The EM algorithm applies the same clustering principles. However, in this procedure, the distribution of the data is assumed to follow a Gaussian mixture model. It iteratively computes between the maximum likelihood estimations (the covariance, mean, and mixing coefficients) and posterior probabilities.

Clustering-based algorithms typically require proper initialization, i.e., initial parameters or segmentation. Among the above methods, the K-means and fuzzy C-means methods are less sensitive to initialization than EM. Since these clustering-based algorithms do not use spatial information, they are sensitive to noise and intensity inhomogeneity. However, because they do not use spatial information, they require less computation time. Moreover, the clustering algorithm highlights a flaw in the stopping criterion for image segmentation.

2.5

Snakes Method

The first deformable model-based segmentation was the Snake method [23, 24, 25]. In this method, the model is deformed iteratively until the explicit curve representation corresponding to the minimum total energy is found to match the object contour in the image. The energy is defined as follows:

E [Cur(s)] = − Z |∇I(Cur(s))|ds + v1 Z |Cur0(s)|2ds + v2 Z |Cur00(s)|2ds (2.4) where Cur(s) is a parametric curve with parameter s, and I is the image. The first term is the external energy and the last two terms are the internal energies of the snake. The external energy is used to drive the curve towards points with high gradient magnitude. The first internal energy is a continuity term to measure the length of the snake (weighted by v1), and the second internal energy is a

smoothness term to measure the stiffness (weighted by v2). The energy function

of the snake is minimized iteratively using an optimization method until good segmentations (curves) are achieved, i.e., ∂E/∂Cur = 0. This equation results in a partial differential equation (PDE), and it is difficult to solve directly in most cases. Instead, an energy optimum is found by evolving the curve in the steepest descent direction of the energy, such that ∂Cur/∂t = −∂E/∂Cur. This is the basic idea behind popular variational methods used in image processing. However, note that any solution will be a local optimum; thus, the result strongly depends on proper initialization or an appropriate starting curve.

Snake methods have some intrinsic problems. The main drawback is the explicit curve representation that cannot handle changes in the topology of the evolving contour. This usually requires complex re-parameterization algorithms. Second, contours cannot split or merge automatically. This makes recognition of multiple objects difficult. Third, since the nodes of the curve are affected by local image features, segmentation is sensitive to proper initialization. This means that the energy will change if the parameterization of the curve changes.

2.6

Level Set Method

The Snakbased deformable model cannot handle topological changes in the e-volving contour directly. Moreover, this is computationally expensive, especially for 3D cases. To address these curve issues, Osher and Sethian [26] devised a level set approach. Rather than manipulating the contour directly, the contour is embedded as the zero level set of a higher dimensional function called the level-set function, ψ (x, t), as shown in Fig 2.5. Over time (t), the level-set function evolves

under the control of a PDE. At some future time, the evolving contour Cur can be obtained by extracting the zero level set Cur = {ψ (x, t) = 0}.

Figure 2.5: Illustration of zero set in a level set [26].

By differentiating with respect to time, the motion of the curve dx/dt can be coupled with evolution of the level set function ψ under the control of PDE, i.e., dψ/dt = −dx/dt∇ψ. When the motion is restricted along the normal direction of the curve, dx/dt = F n, where the normal n = ∇ψ/|∇ψ|. F is referred to as the speed function. The basic idea is to generate a speed function that pulls the curve towards the boundaries of the target object while the curve is regularized with curvature motion. A typical speed function can be based on image gradients, such that it approaches zero when the norm of the image gradient is large (i.e., there is a distinct edge).

The ITK exampleShapeDetectionLevelSetFilter.cxxillustrates how to segment images using the filter itk::ShapeDetectionLevelSetImageFilter. The implementa-tion of this filter in ITK is based on the paper by Malladi et al [27]. This filter requires two inputs: one is an initial level set (or an initial contour) and another is a feature image. Fig 2.6 shows the whole segmentation process of an abdominal CT image using the level set method.

The level set [28, 29] approach is a popular method for organ segmentation. The major advantage of the level set approach is that it allows arbitrary topological changes. As a result, the initial curve does not require the same topology as the segmented object. Unlike the Snake algorithm, which was motivated by energy minimization, the level set approach was motivated by geometric curve evolution. However, this motivation is based on the geometrical aspects of the image and the curve; therefore, it is relatively insensitive to initial conditions but sensitive to noise.

Figure 2.6: Segmentation process based on the shape detection level set method. (a) The input image and one seed point; (b) Gradient magnitude of the smoothed image; (c) Sigmoid of the gradient magnitude, where the sigmoid is used to com-pute the speed term for the front propagation; (d) The segmentation result.

2.7

Geodesic Active Contour Method

The geodesic active contour (GAC) model proposed by Caselles et al. [30, 31] provides an alternative model for edge detection that is derived from the classical Snakes model. This approach can be considered a special case of the classical Snakes model approach because the underlying energy can be interpreted as the geodesic distances (locally shortest paths) of a contour in Riemannian space [32] with a metric induced by the image intensity.

Let us consider a particular class of the Snake model where the rigidity coeffi-cient is set to zero, i.e., v2 = 0. Assume that the Equation 2.4 can be reduced as

follows. E [Cur(s)] = − Z |∇I(Cur(s))|ds + v1 Z |Cur0(s)|2ds (2.5) Observe that, by minimizing the above Equation 2.5, we are attempting to locate the curve at the points of maxima |∇I| (acting as ”an edge detector”) while maintaining a certain smoothness in the curve (the object boundary). Equation 2.5 can be extended by generalizing the edge detector part. Here, let g be a strictly decreasing function. Thus, the gradient − |∇I| can be replaced by an edge detector function g(|∇I|)2. As a result, the energy to be minimized has the following formulation: EGAC[Cur(s)] = Z L(Cur) 0 g(|∇I(Cur(s))|)2ds + v1 Z L(Cur) 0 |Cur0(s)|2ds (2.6) where L(Cur) is the length of Cur and g (|∇I|) is a decreasing edge indicator function, i.e., g (|∇I|) = 1/(1 + |∇I|). Observe that minimizing the energy

func-tion EGAC[Cur(s)] gives a minimal curve, where the length of the curve Cur(s) is

weighted by the function g (|∇I|). The geodesic active contours are derived from Equation 2.6.

The ITK example GeodesicActiveContourImageFilter.cxx illustrates how to segment images using the filter itk::GeodesicActiveContourLevelSetImageFilter. The implementation of this filter in ITK is based on the paper by Caselles [30]. This filter requires two inputs: an initial level set (or an initial contour) and a feature image. Fig 2.7 shows the whole segmentation process of an abdominal CT image using the GAC method. The preprocessing of the GAC method is similar to the level set method.

Figure 2.7: Segmentation process based on the GAC method.

The GAC model [33, 34] does not suffer from some of the limitations of the classical Snake model. Due to the intrinsic nature of this model, changes in the geometry of contours are handled automatically; thus, multiple objects can be recognized simultaneously. Furthermore, this approach uses fewer parameters; therefore, the need for preprocessing the image is reduced. This approach is less prone to variation on the edges, thereby allowing an object with non-ideal edges to be recognized.

2.8

Statistical Shape/Appearance Model Method

Statistical models of shape and appearance are powerful tools for interpreting med-ical images. They provide a standard statistmed-ical system for analyzing the variations in shape and/or texture (intensity) of the same organ. In medical imaging applica-tions, the general shape, location and intensity of an anatomical organ is typically known a priori. This information can be used to aid the segmentation process especially when image contrast is low or when the object boundary is indistinct. The prior knowledge (shape and intensity) can be used to constrain the topology and general shape of the result. When a set of training samples for the same organ is available, standard statistical shape/intensity analysis can be applied. Model