Nearly Automatic Liver Contour Segmentation

Nearly Automatic Liver Contour Segmentation A thesis submitted in fulfillment of the requirements for the degree of Master of Science By Ofer Eliassa...
Author: Gregory Cross
3 downloads 0 Views 636KB Size
Nearly Automatic Liver Contour Segmentation A thesis submitted in fulfillment of the requirements for the degree of Master of Science By

Ofer Eliassaf Supervised By Prof. Leo Joskowicz The Selim and Rachel Benin School of Computer Science and Engineering The Hebrew University of Jerusalem Jerusalem, Israel May 27, 2009

Contents 1 Introduction 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 1

1.2 1.3 1.4

Computed tomography (CT) . . . . . . . . . . . . . . . . . . . . . . Liver anatomy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Method overview . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 3 5

1.5 1.6

Novel aspects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Thesis organization . . . . . . . . . . . . . . . . . . . . . . . . . . .

6 8

2 Literature Overview 2.1 Medical image segmentation overview . . . . . . . . . . . . . . . . .

9 9

2.1.1 2.1.2

Clustering methods . . . . . . . . . . . . . . . . . . . . . . . 10 Thresholding . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.1.3 2.1.4 2.1.5

Region growing . . . . . . . . . . . . . . . . . . . . . . . . . 11 Graph partitioning . . . . . . . . . . . . . . . . . . . . . . . 12 Semi automatic techniques . . . . . . . . . . . . . . . . . . . 12

2.2

2.1.6 Level Set methods . . . . . . . . . . . . . . . . . . . . . . . 12 Maximum likelihood and Gaussian density estimation . . . . . . . . 13

2.3 2.4 2.5

Bayesian classification . . . . . . . . . . . . . . . . . . . . . . . . . 14 Anisotropic diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Active contours and level set methods . . . . . . . . . . . . . . . . . 17

2.6

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3 The Method 19 3.1 Method considerations . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2

Method overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 i

3.3 3.4 3.5

Intensity model generation . . . . . . . . . . . . . . . . . . . . . . . 22 Voxel classification . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Morphological operations . . . . . . . . . . . . . . . . . . . . . . . . 25 3.5.1 3.5.2 3.5.3

3.6 3.7

Connected component . . . . . . . . . . . . . . . . . . . . . 26 Remove holes . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Smart opening . . . . . . . . . . . . . . . . . . . . . . . . . 27

Geodesic active contour refinement . . . . . . . . . . . . . . . . . . 29 Algorithm iterations and resolutions . . . . . . . . . . . . . . . . . . 32

4 Experimental Results 4.1 4.2 4.3

34

Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 Hadassah Ein-Karem datasets . . . . . . . . . . . . . . . . . . . . . 34 MICCAI data sets . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

5 Conclusions 5.1 5.2

45

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

ii

List of Figures 1.1

A multi slice CT scanner . . . . . . . . . . . . . . . . . . . . . . . .

2

1.2 1.3 1.4

A CT slice containing a liver (marked in red) . . . . . . . . . . . . . Liver structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Method output: a 3d, axial, sagital and coronal views of a segmen-

4 5

tation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.1

Region growing core algorithm . . . . . . . . . . . . . . . . . . . . . 11

3.1 3.2 3.3

Liver, liver tumor and liver vessels gray level probability functions . 21 Algorithm core . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 Anisotropic diffusion smoothing . . . . . . . . . . . . . . . . . . . . 26

3.4 3.5

Connected component operation . . . . . . . . . . . . . . . . . . . . 27 Remove holes operation . . . . . . . . . . . . . . . . . . . . . . . . 28

3.6 3.7 3.8

Smart opening operation . . . . . . . . . . . . . . . . . . . . . . . . 30 Sigmoid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Geodesic active contours . . . . . . . . . . . . . . . . . . . . . . . . 32

4.1 4.2

A segmentation result of MICCAI liver-seg013 data set . . . . . . . 38 Segmentation results of some MICCAI test cases. . . . . . . . . . . 40

iii

List of Tables 4.1

Example of the liver class iteration results . . . . . . . . . . . . . . 35

4.2 4.3 4.4

Hadassah Ein-Karem datasets results . . . . . . . . . . . . . . . . . 37 MICCAI train sets results for robustness check . . . . . . . . . . . . 43 MICCAI test sets results as reported by MICCAI grand challenge organizers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

iv

Acknowledgements I would like to thank the people that helped me during the thesis research. First, I would like to thank my thesis supervisor Prof. Leo Joskowicz for his support, advice and guidance during my thesis. Moti Freiman, a Phd. candidate at Prof. Joskowicz lab, for his good advises, ideas and his effort in writing the articles. Next, Dr. Sosna of Haddassah Medical Center. Dr. Sosna provided us the CTA data sets of livers, as well as some quantitative information. Dr. Sosna also explained to us all we needed to know in CTA and liver analysis. I wish to express my warm gratitude to my fellow student and good friend Yoav Taieb, for his support during the thesis and for his friendship. Yoav headed liver tumors segmentation research which was a combine effort with my research in Prof. Joskowicz lab. Most importantly I want to dedicate this thesis and to thank my wife Alexandra for her love, support and understanding during the days and nights I spent on the research and the thesis writing.

v

Abstract Accurate liver segmentation from abdominal computed tomography (CT) images is an important step for computer aided diagnosis (CAD) for liver CT. It can be used in many clinical applications such as assessment of hepatomegaly and liver cirrhosis, liver regeneration after hepatectomy, hepatic transplantation planning for size compatibility, and monitoring of patients with liver metastases, where the disease is related to an enlargement of the liver. To be of practical clinical use, the liver segmentation and volume estimation must be accurate, robust, fast, and nearly automatic, so that the radiologist can perform it without the help of a technician. This thesis presents a new algorithm for a nearly automatic liver segmentation and volume estimation in abdominal CT datasets and its clinical evaluation and validation. The algorithm uses a multi resolution iterative scheme which repeatedly applies smoothed Bayesian classification followed by set of morphological operators and active contours refinement. It uses multi-class and voxel neighborhood information to compute an accurate intensity distribution function for each class. Its advantages are that it requires only one single user-defined seed voxel inside the liver, with no manual adjustment of internal parameters. We conducted experiments on two clinically validated datasets totaling 56 CTAs. The first dataset included 26 CT scans taken from the “Hadssah-EinKarem” hospital. We compared our method to a ground truth manual segmentation performed by a radiologist and to two commercial methods. The second dataset was taken from the Segmentation of the Liver Competition 2007 (SLIVER07) [15]. This competition was started as part of the workshop 3D Segmentation in the Clinic: A Grand Challenge, on October 29, 2007 in conjunction with MICCAI 2007. The algorithm was tested on new unseen datasets by using 5 different scoring systems which was combined to a total score. The results show that the method is on par with the state of the art in liver segmentation. Our method is accurate, robust, efficient, and easy to use. It is an effective tool for hepatic volume estimation in a clinical setup.

Chapter 1 Introduction This thesis presents a new method for a nearly automatic liver segmentation. This chapter introduces the topics of the thesis. Section 1.1 discusses some background topics. Section 1.2 presents a brief explanation on the Computed Tomography technology. Section 1.3 presents some background on the liver anatomy. Section 1.4 provides a brief overview of the method. Section 1.5 discusses the novel aspects of the method. Section Section 1.6 presents the organization of the thesis.

1.1

Background

Liver segmentation and volume estimation from abdominal CT scans is a key task in many clinical applications. They include the assessment of hepatomegaly and liver cirrhosis, liver regeneration after hepatectomy, hepatic transplantation planning for size compatibility, and monitoring of patients with liver metastases, where the disease is related to an enlargement of the liver, among others. To be of practical clinical use, the liver segmentation and volume estimation must be accurate, robust, fast, and nearly automatic, so that the radiologist can perform it without the help of a technician. Nearly automatic liver segmentation, necessary for hepatic volumetry, is know to be a very challenging task due to the ambiguity of liver boundaries, the complexity of the liver surface, and the presence of surrounding organs. Many previous works propose automatic or semi-automatic methods for liver segmentation. These methods can be classified into four main categories:

1

Figure 1.1: A multi slice CT scanner

1) atlas-guided methods [24] 2) region growing methods [18, 29] 3) deformable models and level sets based methods [2] 4) classification based methods [23]. In many cases, the methods produce segmentations with errors, as they do not make full use of the intensity and geometric information present in the original images. To circumvent this problem, several works [12] propose a hybrid two-phase approach consisting of a preliminary segmentation (usually atlas- or intensity-based) followed by a deformable model refinement. However, when the preliminary segmentation is off (e.g. due to a local minima), it is very unlikely that the refinement step will manage to significantly improve it.

1.2

Computed tomography (CT)

A CT (Computerized Tomography) scanner is a special kind of X-ray machine. In standard X-rays, a beam of energy is aimed at the body part being studied. A plate behind the body part captures the variations of the energy beam after it passes through skin, bone, muscle, and other tissue. While much information can be obtained from a regular X-ray, a lot of detail about internal organs and other structures is not available. In Computed Tomography, the X-ray beam moves in a circle around the body. This allows many different views of the same organ or structure, and provides much greater detail. The X-ray information is sent to a computer that interprets the

2

X-ray data and displays it in 2-dimensional form on a monitor. Newer technology and computer software makes three-dimensional (3-D) images possible. CT scans may be done with or without contrast.A CT image taken with contrast is called CTA. ”Contrast” refers to a substance taken by mouth or injected into an intravenous (IV) line that causes the particular organ or tissue under study to be seen more clearly. Contrast examinations may require you to fast for a certain period of time before the procedure. The CT image units are called Hounsfield Units (HU) which measures radiodensity, which is the relative transparency of a substance to X-ray radiation. The technique of CT scanning was developed by the British inventor Sir Godfrey Hounsfield, who was awarded the Nobel Prize for his work. CT scanners first began to be installed in 1974. Because of advances in computer technology, CT scanners have vastly improved patient comfort because they are now much faster. These improvements have also led to higher-resolution images, which improve the diagnostic capabilities of the test. CT scans have already allowed doctors to inspect the inside of the body without having to operate or perform unpleasant examinations. CT scanning has also proven invaluable in pinpointing tumours and planning treatment with radiotherapy.

1.3

Liver anatomy

The liver is located in the upper right-hand portion of the abdominal cavity, beneath the diaphragm, and on top of the stomach, right kidney, and intestines. Shaped like a cone, the liver is a dark reddish-brown organ that weighs about 3 pounds [1]. There are two distinct sources that supply blood to the liver: oxygenated blood flows in from the hepatic artery and nutrient-rich blood flows in from the hepatic portal vein. The liver holds about one pint (13 percent) of the body’s blood supply at any given moment. The liver consists of two main lobes, both of which are made up of thousands of lobules. These lobules are connected to small ducts that connect with larger ducts to ultimately form the hepatic duct. The hepatic duct transports 3

Figure 1.2: A CT slice containing a liver (marked in red)

the bile produced by the liver cells to the gallbladder and duodenum (the first part of the small intestine). The liver regulates most chemical levels in the blood and excretes a product called bile, which helps carry away waste products from the liver. All the blood leaving the stomach and intestines passes through the liver. The liver processes this blood and breaks down the nutrients and drugs into forms that are easier to use for the rest of the body. More than 500 vital functions have been identified with the liver. Some of the more well-known functions include the following: • production of bile, which helps carry away waste and break down fats in the small intestine during digestion • production of certain proteins for blood plasma • production of cholesterol and special proteins to help carry fats through the body • conversion of excess glucose into glycogen for storage (glycogen can later be converted back to glucose for energy)

4

Figure 1.3: Liver structure

• regulation of blood levels of amino acids, which form the building blocks of proteins • processing of hemoglobin for use of its iron content (the liver stores iron) • conversion of poisonous ammonia to urea (urea is an end product of protein metabolism and is excreted in the urine) • clearing the blood of drugs and other poisonous substances • regulating blood clotting • resisting infections by producing immune factors and removing bacteria from the bloodstream When the liver has broken down harmful substances, its by-products are excreted into the bile or blood. Bile by-products enter the intestine and ultimately leave the body in the form of feces. Blood by-products are filtered out by the kidneys, and leave the body in the form of urine.

1.4

Method overview

The segmentation begins when the user selects a seed point inside the liver. This is the only interaction of the user with the method. After the seed selection the 5

method is activated. In each iteration the first stage is building a model which will later be used by a Bayesian classifier. In the first iteration the model is learned by looking in the seed point selection. In other iterations the model is learned by looking at the previous iterations. In the third step we apply some morphological operations for fixing the classification errors and for solving some “edge cases” problems. In the fourth step we apply an active contours algorithm for fine tuning the result. The four steps described above are performed in sequence iteratively until convergence. The iterations are necessary for fine-tuning the intensity model to improve the classification and to minimize the importance of the initial voxel seed selection. After each iteration, the internal parameters of the multi-class intensity model are updated before they are used in the smoothed Bayesian classification step. This coupling is designed to overcome a biased classification due to ambiguous liver boundaries and biased seed selection. To speed up the segmentation and make it more robust and accurate, we use a multi resolution approach. The first few iterations are performed on a down sampled CTA dataset to obtain a rough contour segmentation. Subsequent iterations are performed on the original CTA dataset until no further improvements can be made.

1.5

Novel aspects

Up until now there is no automatic algorithm for segmenting the liver that can actually be used in a clinical routine. Most of the algorithms that are being used clinically are semi-automatic and they require many user interaction (such as taking seed points and parameters adjustments). Most of the previous methods for segmenting the liver are using statistical algorithms. The main drawbacks of such methods are that they require building some statistical model, which requires physiological research of its own. They are very sensitive to the process by which the statistics was collected. If some special case will be presented to those methods they will likely fail. Expanding those methods to other organs require additional research per organ. Our method suggests different concept. To segment the liver we do not define only the liver properties, we also define the background properties. This concept 6

Figure 1.4: Method output: a 3d, axial, sagital and coronal views of a segmentation

is not used by any of the state of the art algorithms in the SLIVER competition. The key advantage of our method is that it segments the liver without the need of a large number of datasets. The method fits itself automatically to a new case even if it is special. We use only one user-selected seed which makes the method almost (or nearly) automatic. The segmentation takes about 6 minutes on average. Experimental results from two datasets totaling 56 CTA images show that the our method is accurate, efficient, and robust to seed selection when compared to manually generated ground truth segmentation. The work was already published in [5] and in [6].

7

1.6

Thesis organization

This thesis consists of five chapters. Chapter 2 is a survey of previous work as well as some literature background used by the method. Chapter 3 provides an overview of our method. Chapter 4 describes our experimental validation and discusses its results. Chapter 5 concludes with a summary of our contributions and presents topics for future work.

8

Chapter 2 Literature Overview This chapter reviews some popular work on medical images segmentation techniques. Section 2.1 provides a brief summary of segmentation techniques. The rest of the chapter goes into details in the methods that we actually used to segment the liver. Section 2.2 describes the maximum likelihood principle and how to estimate Gaussian probability parameters according to the maximum likelihood principle. Section 2.3 describes the Bayesian classification algorithm we used. Section 2.4 describes the anisotropic diffusion smoothing algorithm and section 2.5 describes the Active contours algorithm.

2.1

Medical image segmentation overview

Image segmentation is defined as the partitioning of an image into non-overlapping, constituent regions which are homogeneous with respect to some characteristic such as intensity or texture [4]. If the domain of the image is given by I, then the segmentation problem is to determine the sets Sk ⊂ I whose union is the entire image I. Thus, the sets that make up a segmentation must satisfy I=

K [

Sk

(2.1)

k=1

T where Sk Sj = φ for k 6= j and each Sk is connected. Ideally, a segmentation method finds those sets that correspond to distinct anatomical structures or regions 9

of interest in the image. Image segmentation plays a crucial role in many medical imaging applications by automating or facilitating the delineation of anatomical structures and other regions of interest. In this section, I will briefly describe several common approaches that have appeared in the recent literature on medical image segmentation.

2.1.1

Clustering methods

Those methods are based on the well known K-means algorithm originally introduced by [14]. This is an iterative algorithm that is used to partition an image into k clusters based on their means as follows: 1. Pick k cluster centers, either randomly or based on some heuristic 2. Assign each pixel in the image to the cluster that minimizes the variance between the pixel and the cluster center 3. Re-compute the cluster centers by averaging all of the pixels in the cluster 4. Repeat steps 2 and 3 until convergence is attained (e.g. no pixels change clusters) In this case, variance is the squared or absolute difference between a pixel and a cluster center. The difference is typically based on pixel color, intensity, texture, and location, or a weighted combination of these factors. k can be selected manually, randomly, or by a heuristic. This algorithm is guaranteed to converge, but it may not return the optimal solution. The quality of the solution depends on the initial set of clusters and the value of k [27, 21].

2.1.2

Thresholding

The input to a thresholding operation is typically a gray scale or color image [7]. In the simplest implementation, the output is a binary image representing the segmentation. Black pixels correspond to background and white pixels correspond to foreground (or vice versa). In simple implementations, the segmentation is determined by a single parameter known as the intensity threshold. In a single pass, each pixel in the image is compared with this threshold. If the pixel’s intensity is higher than the threshold, the pixel is set to, say, white in the output. If it is less than the threshold, it is set to black. In more sophisticated implementations, multiple thresholds can be specified, so that a band of intensity values can be set to white while everything else is set to black. For color or multi-spectral images, it 10

Input Seed point inside the object Output Segmentation of the object Procedure main-class.empty main-class.push(seed) while (main-class.has-neighbors()) do a) x = pop-neighbor(main-class) b) if (F (X)) main-class.push(x) Figure 2.1: Region growing core algorithm F (X) is a function represents “similarity” may be possible to set different thresholds for each color channel, and so select just those pixels within a specified cuboid in RGB space. Another common variant is to set to black all those pixels corresponding to background, but leave foreground pixels at their original color/intensity (as opposed to forcing them to white), so that that information is not lost. Whereas the conventional thresholding operator uses a global threshold for all pixels, adaptive thresholding changes the threshold dynamically over the image. This more sophisticated version of thresholding can accommodate changing lighting conditions in the image, e.g. those occurring as a result of a strong illumination gradient or shadows.

2.1.3

Region growing

The main idea behind the “Region Growing” family of segmentation methods is that all the pixels belonging to the segmented object are connected and are “similar” [26]. The generic Region Growing algorithm is given in Fig. 2.1. The difference between different types of region growing methods is in the predicate “F”. Where the most commonly used region growing criteria is based on the observation that an objects gray values are usually within some range around a mean value.

11

2.1.4

Graph partitioning

Graphs can often used in image segmentation. The pixels (or groups of pixels) are considered as vertices and the edges define the (dis)similarity among the neighborhood pixels. Some popular algorithms of this category are minimum mean cut, minimum spanning tree-based algorithm, normalized cut, etc. The normalized cuts method was proposed by Shi and Malik [22]. In this method, the image being segmented is modelled as a weighted, undirected graph. Each pixel/voxels is a node in the graph, and an edge is formed between every pair of pixels [27]. The weight of an edge is determined by the similarity between the pixels. The image is partitioned into disjoint sets (segments) by removing the edges connecting the segments. The optimal partitioning of the graph is the one that minimizes the weights of the edges that were removed (the cut). This algorithm seeks to minimize the normalized cut, which is the ratio of the cut to all of the edges in the set.

2.1.5

Semi automatic techniques

In this kind of segmentation, the user segments the image with the mouse in a rough manner and algorithms are applied so that the path that best fits the edge of the image is shown. Techniques like Livewire or Intelligent Scissors are used in this kind of segmentation [27].

2.1.6

Level Set methods

Those types of segmentation are based on curve propagation. The main idea behind those techniques is to evolve a curve towards the lowest potential of a cost function, where its definition reflects the task to be addressed and imposes certain smoothness constraints [27]. While those techniques can be very efficient, they suffer from various limitations like deciding on the sampling strategy, estimating the internal geometric properties of the curve, changing its topology, addressing problems in higher dimensions, etc.. In section 2.5 we provide details of the active contours algorithm.

12

2.2

Maximum likelihood and Gaussian density estimation

First we will define the principle of Maximum Likelihood (ML). Let x1 , . . . , xm be random variables from a random sample from a discrete distribution whose joint probability distribution is P (x|θ) where x = (x1 , . . . , xm ) is a vector in the sample and θ is a parameter from some parameter space (which could be a discrete set of values say class membership) [20]. When P (x|θ) is considered as a function of θ it is called the likelihood function. The ML principle is to select the value of θ that maximizes the likelihood function over the observations (training set) x1, . . . ., xm . If the observations are sampled independent and identically distributed (a common, not always valid, assumption), then the ML principle is to maximize: θ ∗ = argmax θ

m Y i=1

m m Y X P (xi |θ) = argmax log( P (xi |θ)) = argmax log(P (xi |θ)) θ

θ

i=1

(2.2)

i=1

since it is more convenient to maximize the log likelihood. This gives us a way to estimate the Gaussian distribution P (x1 , . . . , xm ) ∼ N(µ, E) with mean vector µ and covariance matrix E. The parameters of the density function are denoted by θ = (µ, E) and for every vector x ∈ Rn we have: P (x|θ) =

1

1

n 2

2Π |E|

1 2

T E −1 (x−µ)

e− 2 (x−µ)

(2.3)

where |E| is the determinant of the covariance matrix. Now, if we assume that the variables x1 , . . . , xm are drawn i.i.d. The parameter estimation would be recovered by taking derivatives of: L(θ) = argmax θ

n X

log(P (xi |θ))

i=1

with respect to θ, i.e., ∇θ L = 0. We have: m

m

Xn X1 1 L(θ) = − log|E| − log(2π) − (xi − µi )T E −1 (xi − µi ) 2 2 2 i=1 i=1

13

Notice that E is a full rank symmetric matrix. By taking partial derivatives with respect to µ we get: m X ∂L −1 =E (µ − xi ) = 0 ∂µ i=1 P Since E is a full rank symmetric matrix we obtain µ = m1 m i=1 Xi . By taking partial derivatives with respect to E we get: m

1 X ∂L = 0 => E = (xi − µi )(xi − µi )T . ∂E m i=1 The results we obtained suggests that the mean and variance which maximizes the likelihood are just the empirical mean and variance.

2.3

Bayesian classification

The Bayesian Classifier technique is based on the so-called Bayesian theorem. Despite its simplicity, This algorithm can often outperform more sophisticated classification methods. The classifier uses empirical data to evaluate/approximate the conditional probability distributions that arise from Bayes’ theorem. This allow one to estimate quantities (probabilities, averages, etc.) about an individual member of a population by combining information from empirical measurements on the individual and on the entire population. Formally, let c1 . . . .ck be the output classes and let x1 . . . .xn be the input variables [19]. The purpose of the lerarner is to classify to which class ci belongs the input variable xj . This can be done if given an input variable xj we calculate P (ci |xj ) to each class cj , and by choosing the class with the maximum probability as the output class. According to Bayes’ theorem we know that: P (ci |xj ) =

P (xj |ci )P (ci) P (xi )

(2.4)

where P (xj |ci ) is called “likelihood”, P (ci) is called “prior”, P (xj ) is called “evidence” and P (ci |xj ) is called “posterior”. By translating this equation to English

14

we get:

likelihood · prior (2.5) evidence Formally, the purpose is to solve the maximization problem (which is the decision posterior =

rule of the classification): Cf inal = argmax ci

P (xj |ci )P (ci ) P xj

(2.6)

It is easy to see that the evidence is just a normalization factor and as such is not affecting the argmax result, so there is no reason to calculate it. In cases where the classes are equaly distributed or unknown, it is possible to as that the prior is the same to all of the classes. This leaves us with the follwing maximization problem: ML : Cf inal = argmax P (xj |ci )

(2.7)

ci

which is known as Maximum Likelihood (ML). When the prior is not uniform we have: MAP : Cf inal = argmax P (xj |ci )P (ci)

(2.8)

ci

When classifying voxels in an image the random variables xj ’s are usually some function of the voxels, such as their gray values, and the classes random variables represent some objects to classify inside the image. The Bayes calculation that is applied creates k probability images (k is the number of classes), and the decision rule should select the best match. It is sometimes useful to smooth the probability maps before applying the ML/MAP decision rules. This can lead to the effect of noise reduction since it is not likely that single voxels will be classified different than their surroundings. When smoothing the probability maps its important not to smooth the boundaries of the objects so that the segmentation near the boundaries will not be hurt.

2.4

Anisotropic diffusion

Many image processing applications are characterized by noisy or corrupted image data. There are some simple filters used to smooth images, where the Gaussian 15

filter is one of the most popular among them. Kondernik [11] showed that the Gaussian filter is equivalent to the soultion of some diffusion equation. The problem is that the Gaussian filter smoothes the edges inside the image. Perona an Malik [17] then presented the Anisotropic diffusion technique. Anisotropic diffusion has been formulated as a process that will mainly take place in the interior of regions, and it will not affect the region boundaries. The basic idea is to evolve a family of smoothed images It from an initial image I0 using the following partial differential equation: It = div(g(||∇I||)∇I) = g(||∇I||)∆I + ∇g · ∇I, I(x, y, 0) = I0 (x, y)

(2.9)

where ∇ indicates the gradient operator with respect to the image, div indicates the divergence operator, t is the iteration step, ||∇I|| is the gradient magnitude, ∆ indicates the Laplacian operator with respect to image and the function g is an “edge-sttoping” function called “diffusion function”. This function is chosen to satisfy g(x) → 0 when x → ∞ so that the diffusion is stopped near edges, for example, g(||∇I||) = e−( or similarly: g(||∇I||) =

|∇I| 2 ) k

1 1 + ( |∇I| )2 k

(2.10)

(2.11)

The discrete version of the diffusion equation can be calculated iteratively by using the following equation, presented by Black et al [3]: Ist+1 = Ist +

λ X g(∇Is,p)∇Is,p |ηs | p∈η

(2.12)

s

where Ist is a discretely samples image, s denotes the pixel direction position in a discrete grid and t denotes discrete time stamp (iteration). The constant λ ∈ R+ is a scalar that determines the rate of diffusion, ηs represents the spatial neighbors of pixel s, and |ηs | is the number of neighbors. they approximated the image gradient magnitude in a particular direction as: ∇Is,p = Ip − Ist , p ∈ ηs 16

2.5

Active contours and level set methods

Active contours, also known as “snakes”, are parametric curves which one tries to fit to an image, usually to the edges within an image [28]. The 3D version is often known as deformable models or active surfaces in the literature. The original model is due to Kaas et al. [10], but many modifications have been proposed in the literature. Let C : [0, 1] → Ω2 denote a curve in the 2D image domain Ω2 . The original energy functional is then given by: E(C) = α |

Z

1 ′

0

2

|C (q)| dq + β {z Eint

Z

1 ′′

0

2

|C (q)| dq − λ } |

Z

0

1

|∇f (C(q))| dq {z }

(2.13)

Eext

The goal is to solve the optimization problem: E(C) = 0. where Eint is the internal energy which is meant to enforce the smoothness of the curve and Eext is the external energy which moves the contour towards object edges. In fact, the active contour model is a regularized gradient edge detector. The coefficients α and β in Eq. 2.13, which must be positive, govern the restoring forces associated with the elasticity and stiffness of the snake respectively. Either of these coefficients may be allowed to vary with q, along the snake. These methods main drawbacks are their inability to move into concavities of an objects boundary, their inability to find the borders when it is initialized too far distant from the actual border location and the model cannot change topology: A curve must always stay closed and is implicitly not allowed to cross.

2.6

Summary

The image segmentation problem has no one good solution since the problem is extremely difficult. I have presented some of the most popular techniques which attempt to deal with this problem. Although they produce reasonable segmentations in many situations, at some point local ambiguities and errors introduced by the segmentation process can only be resolved by application specific knowledge. The liver segmentation is very complicated so we had to make use of some of the techniques presented here, however, none of the methods worked well “out of the 17

box” so we had to use many of the methods and we had to insert some heuristics that without them the results were poor.

18

Chapter 3 The Method This chapter describes a method for segmenting the liver inside a CT scan. The content of this chapter is as follows. Section 3.1 describes the considerations that led to the given method. Section 3.2 describes the work flow of the method. Section 3.3 provides an overview of the intensity model generation. Section 3.4 provides an overview on the Bayesian voxels classification. Section 3.5 provides an overview on the morphological operations. Section 3.6 provides overview on the active contour refinement. Section 3.7 explains the need in working with multiple resolutions and with multiple iterations.

3.1

Method considerations

Nearly automatic CT-based liver segmentation is known to be a very challenging task. The main difficulties include the ambiguity of the liver boundary, the complexity of the liver surface, the presence of surrounding organs, the contrast variability between liver parenchyma and vessels, the different tumor sizes and shapes, and the presence of many small metastases. In addition, liver segmentation often requires the simultaneous identification and visualization of liver contours, blood vessels, and tumors, each with their own features and characteristics. Many statistical methods for organ segmentation exists and some can show reasonable results in liver segmentation. However, those statistical techniques will not work in the same success when segmenting the liver with the presence of tumors

19

since they are characterized by different shape and sizes. This lead us to choose a different type of algorithm. When analyzing the edges of the liver inside a CT scan it is sometimes impossible to determine where exactly the edges are by just observing the derivatives value of the voxels on the edges. There are places where the derivatives are very strong (near the lung, for example) and there are places where the derivatives do not exist at all (if the liver intersects with the spleen). This means that it is impossible to base the liver segmentation entirely on edges based algorithms, but it is possible to use edge based segmentation in addition to other methods. When analyzing the gray level histograms of the liver it seems that they are distributed normally with different parameters, which depends on the CT machine, contrast matter and more. figure 3.1 shows an example of such histogram of a given liver. Since we want our algorithm to work on a wide range of CT scans and on a wide range of liver examples in a nearly automatic way, the algorithm should automatically adjust itself to a given example. This lead us to base our algorithm on gray level classification algorithm which will be based on the normal distribution feature, but because of the variance of the normal distribution parameters it is important to adjust this classification to a given example. The active contours family of algorithms may be in use in liver segmentation. Those algorithms will not work correctly if the initial guess for them will be too far from the actual shape of a given liver, since the shape of the liver is very complex and since the boundaries are not always clear. However, if a good initial guess will be given as inputs to this family of algorithms, there is a big added value to them, since they can find the contour in a sub-pixel accuracy over the edge of the real boundary of the liver.

3.2

Method overview

The segmentation begins when the user takes a seed point inside the liver. This is the only interaction of the user with the method. After the seed selection the method is activated. In each iteration the first stage is for learning and adjusting our model which is used by the next step - to fit the current liver example. The second stage is for using the gray level distribution feature with a Bayesian classi20

Figure 3.1: Liver, liver tumor and liver vessels gray level probability functions The x axis represents gray values (Hunsfield units) and the y axis the probability for the gray value Note that the gray levels are distributed normally. fication, according to the model produced from the previous step. The third step is for fixing the classification errors which are caused because of the mathematics do not always fit exactly to the real life. The fourth step is for fine tuning the contour of the liver in a sub-pixel segmentation. The four steps described above are performed in sequence iteratively until convergence. The iterations are necessary for fine-tuning the intensity model to improve the classification and to minimize the importance of the initial voxel seed selection. After each iteration, the internal parameters of the multi-class intensity model are updated before they are used in the smoothed Bayesian classification step. This coupling is designed to overcome a biased classification due to ambiguous liver boundaries and biased seed selection. To speed up the segmentation and make it more robust and accurate, we use a multi resolution approach. The first few iterations are performed on a down sampled CTA dataset to obtain a rough contour segmentation. Subsequent iterations are performed on the original CTA 21

Input CT scan with liver and one seed point inside the liver Output Segmentation of the liver Procedure For each level in the pyramid (resolution) from highest to lowest do: Repeat until convergence: 1. multi-class intensity model generation 2. smoothed Bayesian classification 3. Morphological operations 4. geodestic active contour-based refinement Figure 3.2: Algorithm core dataset until no further improvements can be made.

3.3

Intensity model generation

The purpose of this stage is producing a model to the Bayesian classification. The main challenge in Bayesian classification is to define the probability to a given voxel to be part of the object and to the background. In order to do so it is necessary to define the background. The first step is the construction of an intensity model that differentiates between the liver , and the remaining organs and tissues, which are interpreted as background. The main challenges are to overcome the ambiguous boundary between the liver parenchyma and the outside organs such as the kidney or spleen. To overcome these difficulties we use the observation that the gray levels of the liver are distributed normally with a given mean and variance parameters, which depends on the CT machine, contrast matter and more. This observation led us to propose to first classify the voxels into a main class which represents the liver, according to its gray level values. For this main class, a five class model is then built. Its central class is the main class (e.g. the liver class) and the remaining four represent other organs and tissues. We model each class with a normal distribution defined by its mean intensity value and variance. In the first iteration, the mean

22

and the variance of a rectangular neighborhood around the manually selected seed inside the liver is computed and interpreted as the parameters of the liver class. We then define 4 new thresholds according to the main class model as follows: µnear−low = µcentral − knear × σcentral µnear−low = µcentral − knear × σcentral µnear−high = µcentral + knear × σcentral µf ar−low = µcentral − kf ar × σcentral µf ar−high = µcentral + kf ar × σcentral

(3.1)

The thresholds help us to produce the remaining 4 classes by learning the average and the variance of the next 4 classes: • Voxels with gray values in the range [µnear−low , µnear−high ] are ignored, since they are in the main class range. • Class far-low - Voxels with gray values lower than µf ar−low . • Class near-low - Voxels with gray values in the range [µnear−low , µf ar−low ]. • Class far-high - Voxels with gray values bigger than µf ar−high . • Class near-high - Voxels with gray values in the range [µnear−high , µf ar−high ]. In subsequent iterations, the segmented region from the previous iteration is used to compute the mean and the variance of the liver class. The remaining 4 classes are generated like in the first iteration. Formally, the liver class is defined as: 2 Xliver ∼ N(µliver , σliver )

(3.2)

The other four classes, Xi are modeled as: Xi ∼ N(µi , σi2 )

23

(3.3)

for each i ∈ {near − low, near − high, f ar − low, f ar − high}, The means and variance are defined as in the first iteration. where the factors Knear and Kf ar are determined from Chebyshev’s inequality: Pr(|Xi − µliver | ≥ Kσliver ) ≤

1 k2

(3.4)

This inequality ensures that at least (1 − k12 ) × 100% of the values are within k standard deviations from the mean. By setting Knear = 2.2 and Kf ar = 4, we ensure that at least 80% of the main class voxels will classified as belonging to the liver, and at least 70% of the voxels that belong to the near classes will classified as near, (regardless of the normal distribution assumption). This ensures that each voxel from the ambiguous class boundary has a high probability of being in both the liver and near classes. We have set the exact factors by running exhaustive search on those parameters. Fig 3.1 describes the Gaussian model of the liver, tumor and vessels classes and their ambiguous boundaries and present a classified liver where the liver is in blue, vessels are in red, and tumors are in green.

3.4

Voxel classification

This step produces a unique classification of each voxel. The classification relies on the voxels intensity values and their neighboring voxels. Neighborhood information is important since voxel intensity values are correlated. Our classification follows the method in Teo et al [16]. First, we compute the probability C(j,i) of a voxel Vj with intensity value vj to belong to class ci where i ∈ {near − low, near − high, f ar − low, f ar − high} using Bayes rule: Pr(vj |cj ) ∝ Pr(cj |vj ) Pr(cj )

(3.5)

where Pr(cj |vj ) is obtained from the intensity model that was generated as described in the previous section 3.3 by using the Gaussian probability density function to each class as follows: P r(cj |vj ) = p

1 2Πσi2

24



e

(vj −µi )2 2σ 2 i

(3.6)

Pr(ci ) is assumed to be from a uniform distribution. The resulting five maps quantify the membership probability of each voxel to each class. Next, we incorporate the neighborhood information to the classification process by smoothing the probability maps C(j,i) for each class i separately using the anisotropic diffusion equation proposed by Perona and Malik [17]: k∇Ik 2 ∂I = ∇ · ∇I, with conduct a = e−( b ) ∂t

(3.7)

where ∇I is the image gradient, a is the conductivity coefficient of the local image edges smoothing, which inversely depends on the image gradient magnitude and on scaling factor b that helps to differentiate between noise and real image gradients. The anisotropic smoothing process smoothes small peeks in the membership probability maps considered as mis-classifications while preserving sharp edges between different objects. We repeat the smoothing operation iteratively k times for getting better results. Using exhaustive search on the MICCAI training set we set the parameter k to 3 for getting the best results. Finally, we apply the Maximum Likelihood (ML) rule: C(j,f inal) = argmax P r(vj |cj )

(3.8)

ci

The decision rule sets the final membership of each voxel to the class with the highest probability that the voxel belongs to it. The binary segmentation image is generated by selecting only the voxels that belong to the liver class.

3.5

Morphological operations

The third step of the algorithm iterations is the identification of the liver region. The Bayesian classification yields most of the liver region, but several problems may occur after that classification: • There are holes inside the liver that can be caused by the presence of vessels (with brighter intensity values) and the tumors (with darker intensity values). • Other disconnected components outside the liver which correspond to differ25

(a)

(b)

Figure 3.3: Anisotropic diffusion smoothing This figure shows the liver class probability image before (a) and after (b) the smoothing. ent organs may be classified as the liver, since they have the same intensity values. • Some organs like the kidney, heart or the spleen may be connected to the liver. We deal with those situations by using some morphological operations.

3.5.1

Connected component

We use this morphological operation to remove the organs which are not connected to the liver. After this operation we remain with one connected component - the one which is classified as the liver class and contains the user seed point. Another usage to this operation is in the “smart opening” operation - after we disconnect the other organs components from the liver component using the opening operation, we remove them by using connected component as will be described later.

3.5.2

Remove holes

This operation is used in order to remove the holes inside the liver that was caused by the presence of tumors and vessels inside the liver. Tumors and vessels have 26

(a)

(b)

Figure 3.4: Connected component operation Figure (a) shows the Bayesian classification result. Each class is labeled in a different color. The transparent color represents the liver class. Figure (b) shows the liver class (in red) after the connected component operation. As can be seen many other organs were removed by this operation. different values of gray levels. The “Remove Holes” operation is a two dimensional operation that we activate on each slice. We paint the background from the top right corner (which is assumed to be a background voxel in all slices). After the paint is done all the voxels with a background class that are not painted are considered as holes. We mark those holes as the liver class. After this operation we have the liver as one connected component with no holes. Other organs such as the spleen, kidney or the heart may still be connected to it.

3.5.3

Smart opening

This operation is used in order to separate some other organs such as spleen or heart from the liver class. The operation is based on the observation that the other organs are being connected to the liver with a small interface radius. In this operation we use the following steps: • Opening operation. The opening is done with a radius big enough to remove the other organs and leave us with the liver. This radius is being determined 27

(a)

(b)

Figure 3.5: Remove holes operation This figure shows “remove hole” operation effect before (a) and after (b) the operation. according to the liver volume since there is a direct correlation between the liver volume and the interface radius with the other organs. This operation makes other organs to be in a different connected component. • Connected component. We paint the liver connected component around the user seed point and we remove the rest of the component. In rare cases it is possible that the user seed is on a very small connected component, if for example the user selected a seed on a vessel or a small tumor. For this reason we paint a neighboring minimal radius around the seed point so that the liver will definitely remain after this stage. • Closing operation. To expand the liver back to its true dimensions we start painting back pixels that were removed by the opening operation and are closer than a second radius to the remaining voxels of the liver. This radius is larger than the first radius. This way the other organs were removed from the liver. Some voxels of the other organs which lie in the interface area still remain labeled as the liver, however those voxels are negligible and their area is shrinking between the algorithm iteration. 28

• Remove holes. The opening operation may caused some holes inside the liver so now is the time to close them so that we can have a liver without holes. As can be seen there is a trade off between two types of kernels 1. Using a big opening kernel which will make other organs be removed for sure, but will also remove many liver voxels. 2. Using small opening kernel which will make all/most of the liver voxels remain but can accidentally make other organs remain connected. This operation uses opening and closing radius (kernels). There is a different spacing between two different CT machines. This lead us to the understating that we need to normalize the radiuses in the spacing of the CT or we won’t be able to set one radius to all the data sets. Another issue is the difference between human beings anatomical structure and sizes. A small radius will work on a child CT but won’t work on and adult and vise versa - a big radius will work fine on an adult but not on a child. These consideration lead us to determine the size of the kernel according to a formula which combines the liver volume with the spacing. radius =

liverV olume 1, 000, 000 × capacityF actor × Spacing

(3.9)

where liverCapacity is the liver volume in the image coordinates, capacity factor is a factor to determine the size of the kernel and spacing is the spacing from the image information (depends on the CT machine). The capacity factor is learned by an exhaustive search process.

3.6

Geodesic active contour refinement

The segmentation result after the morphological operators may not good enough. The “smart opening” operation can cause the object contour to shrink a little so there may be a need for a minor correction of the segmentation. In addition, the previous stages didn’t produce a sub-pixel accuration of the contour. To fix this problem we use active contours segmentation. When given a good initial guesses and a good energy functions the active contours yield very good results. Before 29

(a)

(b)

Figure 3.6: Smart opening operation Figure (a) shows the classified liver before the “smart opening” operation. As can be seen the heart (the big component from the right) is classified as the liver since it is connected in the 3d. 2 other small components are also classified as the liver. Figure (b) show the result of the operation - the heart and the other components are is no longer connected to the liver. applying the active contour we create a feature image out of the original image in the following way: We first activate a threshold function as follows:    µliver + σliver if I(x) > µliver − σliver ′ I (x) = I(x) if µliver − σliver ≤ I(x) ≤ µliver + σliver   µliver − σliver if I(x) > µliver − σliver

(3.10)

where I is the original CTA data, I’ is the new image, and µliver and σliver are the liver class parameters as computed in the intensity model. This thresholding in necessary since we do not want that irrelevant edges and intensity values which belongs to other organs will influence on our feature image. After the threshold operation we smooth the results using an anisotropic diffusion to preserve the edges of the liver. Next, the gradients are computed on the resulting 3d volume. After the gradient volume is calculated we enhance it by using a sigmoid function. 30

Figure 3.7: Sigmoid

The sigmoid function, also called the sigmoid curve (von Seggern 2007, p. 148) or logistic function, is the function: 1 1 + e−x

(3.11)

In our implementation we use the function: y = (max − min)

1 + min (1 + e−(x−α) )

(3.12)

where max is the maximal gradient value and min is the minimal gradient value. The parameter α specify the center of the “S” shape and the normalization in max and min is for bounding the result between max and min. This function makes strong edges (bigger than α) in the gradient image be stronger (close to max), weak edges (smaller than α) be weaker (close to min) and it is not changing the edges with gradient intensity of α. The parameter α was set to 15 after running an exhaustive search on the MICCAI database. The sigmoid function result is chosen as the feature map to the active contour function. We use the Geodesic active contour algorithm. The geodesic active contours technique is based on active contours evolving in time according to intrinsic geometric measures of the image. The evolving contours naturally split and merge, allowing the simultaneous detection of several objects and both interior and exte31

(a)

(b)

Figure 3.8: Geodesic active contours Figure (a) shows the classified liver before applying the geodesic active contours algorithm. Figure (b) shows the algorithm result. rior boundaries. The proposed approach is based on the relation between active contours and the computation of geodesics or minimal distance curves. We chose the geodesic active contours algorithm since it allows stable boundary detection when their gradients suffer from large variations including gaps.

3.7

Algorithm iterations and resolutions

Another important aspect of the algorithm is its iteration scheme. The classification is done iteratively in each resolution until the Bayesian classes are converging. This is done for the following reasons: • Robustness. This means that two different seeds should yield a very similar results. • Fine tuning to the algorithm results. • Algorithm results confidence: In some situations, mainly when a “bad” seed is chosen by the user (if the seed surrounding area is not fully inside the 32

liver). The algorithm results can sometimes not converge. This means that the confidence of the result is low in such cases. The only examples we encountered are bad seed selection of the user. Another important aspect of the algorithm is the work in multiple resolutions. This is done since the work in the given CT resolution can take a long time. It is impossible to run many iterations in the highest resolution thus it is necessary to give the highest resolution a very fine-tuned guess, so that only one iteration will be done in this resolution. This way most of the work is done in the lower resolutions and takes less time, and we still get very accurate results.

33

Chapter 4 Experimental Results This chapter will describe the experimental results. Section 4.1 will describe the experimental setup, section 4.2 will describe the experiments and the results on the Hadassah CTA datasets and section 4.3 will describe the experiments and the results on the MICCAI datasets.

4.1

Experimental setup

We implemented our algorithm using the ITK software library [13]. The Bayesian module (implemented also in ITK) was taken from [8]. We used an Intel Core2 Dual Quad 2.8 GHz PC with 2GB of RAM for all the experiments that will be described in this chapter. We evaluated our algorithm on two CTA datasets. The first dataset was taken from the radiology department at Hadassah Ein Kerem hospital from Dr. Sosna’s lab. The second dataset was taken from the MICCAI 2007 workshop as part of a grand challenge competition.

4.2

Hadassah Ein-Karem datasets

This experiment was conducted on 26 CTA datasets. The datasets were acquired from a Brilliance 16, Philips Medical Systems. Its main characteristics are: 16 channel detector scanner with 2 mm slice thickness and 1 mm increment, with 200-250 mAs, 120 kVp. 34

#Iter. 1 2 3 4 5

Near Avg. 64.34 71.40 72.00 72.01 72.01

Low Var. 84.78 100.27 105.94 105.99 105.99

Near Avg. 160.14 140.56 140.22 140.21 140.21

High Var. 108.52 111.62 105.74 105.71 105.77

Far Avg. -671.82 -608.76 -604.63 -604.61 -604.61

Low Var. 194631 203142 204609 204615 204615

Far Avg. 386.17 344.89 343.68 343.65 343.65

High Var. 31450.0 32251.5 32255.5 32255.4 32255.4

Object Avg. Var. 112.38 253.53 109.21 274.10 109.28 268.55 109.26 267.62 109.26 267.50

Table 4.1: Example of the liver class iteration results The datasets were taken from 26 adults patients with an average age of 54 years, all in the range of 32-78. The clinical findings of the datasets included 13 healthy livers, 3 livers with cirrhosis, 4 with metastases, and 6 with primary liver tumors. The datasets were previously used in a study by Sosna et al. [9]. The ground truth of the experiment was obtained from an experienced technologist in a manual segmentation technique of the liver contours. The manual segmentation was performed on Vitrea 2 workstation (Vital Images, Plymouth, MN). We also compared our result to a fully automatic prototype software and to a semi automatic paintbrush tool in a Philips EBW 3.0 workstation, (Philips Medical Systems, Groningen, The Netherlands). The results of those methods were reported in [9]. After we segmented the liver contours we compared our method to the other methods by comparing the following aspects: Liver Volume The liver volumes were compared by using absolute differences and by liner regression to the ground truth. The mean absolute difference between the estimated liver volumes between the manual tracing and our method was 6.48% (std 5.08%) of the liver volume, with Philips automatic method (2) 10.55% (std 5.98%), and with Philips semi-automatic method (3) 2.65% (std 1.45%). This means that our algorithm estimated the liver volume more accurately than Phillips automatic method but less accurately than the semi-automatic method. The calculated correlation coefficient of the liver volume between the manual segmentation and our method measurements was 0.96 (which is an excellent correlation). The correlation between Philips automatic method and the ground truth was 0.97, and between the semi-automatic method and the 35

ground truth was 0.99. Robustness We compared our method to the Philips automatic method by comparing volume differences of 3 different initial seed points. the mean difference between the different initializations using our method was 1.1% (std 2.34) and using Philips automatic method was 2.09% (std 2.85), which indicates that our algorithm is highly robust to seed selection. Table 4.1 shows an example of the resulting classes during the iterations. It can be seen that the algorithm converge very fast, so actually not more than 4 iterations are needed in each resolution. Run time We recorded the run time of each method. The mean segmentation time using the manual tracing was 7:24 minutes (range 4:32 - 10:26 minutes). Our method mean running time was 6:05 minutes. (range 4:22 - 8:13 minutes.) Note that only a few seconds are required from the radiologist for seed selection which makes it even more suitable for clinical use. The Philips automatic method running time was 0:37 mins (range 0:28 - 1:13 minutes). The Philips automatic segmentation was done in a dedicated workstation and not on a standard PC and therefore it is not directly comparable to our method running time. The segmentation time using Philips semi automatic method was 5:14 mins (range 3:47 - 6:41 minutes) but requires much more user interaction time. The conclusions from the measurements are that our method estimated the liver volume much better than the Philips automatic algorithm. The Philips semiautomatic method achieved more accurate results but required 5 minutes interaction time compared to our method, which only required a few seconds for seed selection. Our method also achieved a much better correlation (0.96 vs. 0.79) than that reported by Doi et al [30] on 35 datasets (albeit different ones, but not available to us). Table 4.2 shows the summary of the experimental results. The first column indicates the method used. The second column shows the estimated volume difference (mean and standard deviation in parethesis). The third column shows the correlation coefficient with the ground-truth manual tracing segmentation. The 36

Method Ours Philips semi-automatic Philips automatic

Volume diff % 6.48 (5.08) 2.65 (1.48) 10.55 (5.98)

Correlation Coefficient 0.96 0.99 0.97

Robustness 1.10 (2.34) 2.09 (2.85)

Table 4.2: Hadassah Ein-Karem datasets results last column show the robustness of the algorithm according different seed selections.

4.3

MICCAI data sets

In the this study we used the datasets from the Segmentation of the Liver Competition 2007 (SLIVER07). The goal of this competition is to compare different algorithms to segment the liver from clinical 3D computed tomography (CT) scans. Teams that participated in the liver segmentation contest of this workshop downloaded training and test data and submitted the results of their algorithms on the test data to the competition organizers. Those results can be found in the competition website [15]. We used 30 datasets with ground truth segmentations, segmented by few experienced radiologists [25]. The database is divided into two groups. The first group consists of 20 datasets and is used for training and testing. The second group consists of 10 datasets and is used to evaluate segmentation algorithms. We compared our segmentation results to the ground-truth by using metrics as decribed in [25] and we got the following results for each of the following metric: 1. Volumetric overlap error (in percents). This is the number of voxels in the intersection of segmentation and reference divided by the number of voxels in the union of segmentation and reference, subtracted from 1 and multiplied by 100. This value is 0 for a perfect segmentation and has 100 as the lowest possible value, when there is no overlap at all between segmentation and reference. We got average of 8.55% on the test sets and 7.81% on the train sets. 2. Relative absolute volume difference (in percent). The total volume difference of the segmentation to the reference is divided by the total volume of the reference. 37

Figure 4.1: A segmentation result of MICCAI liver-seg013 data set

The result is multiplied by 100. An undersegmentations will provide negative values and oversegmentation will provide positive values. To compute the score, the absolute value is taken. Note that the perfect value of 0 can also be obtained for a non-perfect segmentation, as long as the volume of that segmentation is equal to the volume of the reference. We got 2.78% on the test sets and 2.63% on the train test. 3. Average symmetric surface distance (in millimeters). The border voxels of segmentation and reference are determined. These are defined as those voxels in the object that have at least one neighbor (of their 18 nearest neighbors) that does not belong to the object. For each voxel along one border, the closest voxel along the other border is determined (using Euclidean distance, not signed, 38

and real world distances, so taking into account the different resolutions in the different scan directions). All these distances are stored, for border voxels from both reference and segmentation. The average of all these distances gives the average symmetric absolute surface distance. This value is 0 for a perfect segmentation. We got 1.46mm on the test sets and 1.28mm on the training sets. 4. Root Mean Square (RMS) symmetric surface distance (in millimeters). This measure is similar to the previous measure, but stores the squared distances between the two sets of border voxels. After averaging the squared values, the root is extracted to give the symmetric RMS surface distance. This value is 0 for a perfect segmentation. We got 2.94mm on the test sets and 2.55mm on the train sets. 5. Maximum symmetric surface distance (in millimeters). This measure is similar to the previous two, but in this case the maximum of all voxel distances is taken instead of the average. This value is 0 for a perfect segmentation. We got 26.72mm on the test sets and 22.77mm on the training sets. Each of the results above was normalized to a score and a total score was computed as described in [25] for both training and test datasets. Our method achieved a total score of 71.8 on the training set and of 67.87 on the test set. This result indicates the the accuracy of the segmentation is very high (almost like a trained human - [25]). It is ranked forth of the semi-automatic method. However, all other semi-automatic methods require significantly more user interaction. The correlation between our volume estimation and the ground truth on the training set was 0.99. Note that since the ground truth is not available for the test set (the organizers of the grand challenge chose to hide it from the competitors), we could not compute the correlation value for it. In addition, we measured the robustness of our method to three different seeds initializations on the 20 training datasets. The absolute volume difference was 0.004%, which is negligible. These results suggest that our method is both accurate and robust to seeds selection.

39

(a) test-001 - sagittal view

(b) test-006 - coronal view

(c) test-007 - axial view Figure 4.2: Segmentation results of some MICCAI test cases. The segmentation outline is given in blue. The reference segmentation is in red.

40

Volume Name

Volumetric Overlap

Relative absolute

Average symmetric

RMS symmetric

Maximum symmetric

Error

volume difference

surface distance

surface distance

surface distance

liver-seg001

9.94

0.33

1.56

2.90

23.12

liver-seg001

9.70

0.10

1.51

2.81

23.85

liver-seg001

10.00

0.43

1.57

2.93

23.26

liver-seg002

8.47

3.16

1.35

2.33

21.48

liver-seg002

8.47

2.91

1.35

2.33

21.48

liver-seg002

8.49

3.30

1.33

2.23

17.82

liver-seg003

8.06

0.59

1.42

2.89

25.23

liver-seg003

7.91

0.31

1.39

2.81

25.23

liver-seg003

8.27

0.78

1.47

3.00

25.09

liver-seg004

6.46

3.10

0.94

1.76

15.66

liver-seg004

6.24

2.92

0.91

1.69

15.39

liver-seg004

6.48

3.17

0.98

1.86

17.59

liver-seg005

6.73

4.35

0.93

1.81

16.67

liver-seg005

6.93

4.25

0.96

1.89

16.67

liver-seg005

6.74

4.40

0.94

1.82

16.67

liver-seg006

7.06

0.54

1.06

2.40

24.72

liver-seg006

7.09

0.72

1.08

2.46

24.88

liver-seg006

7.10

0.71

1.08

2.46

24.88

liver-seg007

6.68

1.48

1.06

2.05

17.46

liver-seg007

6.82

1.88

1.09

2.10

18.24

liver-seg007

6.88

1.32

1.14

2.22

20.11

liver-seg008

8.40

1.13

1.52

2.81

25.28

liver-seg008

8.40

1.29

1.50

2.73

25.38

liver-seg008

8.35

1.43

1.48

2.69

25.38

liver-seg009

6.94

3.35

1.04

1.76

18.35

liver-seg009

6.95

3.41

1.04

1.76

18.35

liver-seg009

6.74

3.07

1.00

1.71

18.92

Continued on Next Page. . . 41

Volume Name

Volumetric Overlap

Relative absolute

Average symmetric

RMS symmetric

Maximum symmetric

Error

surface distance 1.43

surface distance 2.46

surface distance 17.71

liver-seg010

8.89

volume difference 6.58

liver-seg010

9.67

6.55

1.56

2.76

27.74

liver-seg010

9.59

6.45

1.54

2.71

27.74

liver-seg011

7.66

2.97

1.41

3.17

23.50

liver-seg011

8.09

0.80

1.51

3.53

26.55

liver-seg011

8.08

0.19

1.58

3.74

37.90

liver-seg012

6.27

2.07

0.93

2.07

21.22

liver-seg012

6.24

2.01

0.92

2.00

20.63

liver-seg012

6.35

2.31

0.95

2.10

20.61

liver-seg013

11.62

5.16

1.87

3.88

34.96

liver-seg013

11.44

5.43

1.83

3.84

35.72

liver-seg013

12.48

5.55

1.96

3.61

32.04

liver-seg014

7.99

2.78

1.21

2.53

20.84

liver-seg014

8.00

2.75

1.21

2.51

20.84

liver-seg014

8.14

3.08

1.23

2.58

20.84

liver-seg015

4.75

3.06

0.69

1.48

20.74

liver-seg015

4.81

3.20

0.70

1.52

20.74

liver-seg015

4.76

3.09

0.69

1.48

20.74

liver-seg016

9.79

3.48

1.98

3.70

23.93

liver-seg016

9.81

3.35

1.99

3.72

23.93

liver-seg016

10.14

3.96

2.05

3.80

23.82

liver-seg017

7.78

1.17

1.25

2.61

23.41

liver-seg017

7.87

1.44

1.32

2.79

23.26

liver-seg017

7.87

1.80

1.33

2.78

23.26

liver-seg018

7.19

3.24

1.27

2.41

19.74

liver-seg018

6.90

1.09

1.13

2.10

17.91

liver-seg018

6.89

1.08

1.13

2.10

17.52

Continued on Next Page. . . 42

Volume Name

Volumetric Overlap

Relative absolute

Average symmetric

RMS symmetric

Maximum symmetric

Error

surface distance 1.26

surface distance 2.96

surface distance 27.06

liver-seg019

7.08

volume difference 3.76

liver-seg019

7.06

3.75

1.29

3.17

28.73

liver-seg019

7.05

3.75

1.29

3.16

28.65

liver-seg020

7.26

2.20

1.20

2.34

21.19

liver-seg020

7.16

2.40

1.23

2.47

22.77

liver-seg020

7.32

2.04

1.25

2.56

22.71

All (average)

7.81

2.62

1.28

2.55

22.77

Table 4.3: MICCAI train sets results for robustness check

Volume Name

Volumetric Overlap

Relative absolute

Average symmetric

RMS symmetric

Maximum symmetric

Error

volume difference

surface distance

surface distance

surface distance

test-001

7.86

3.34

1.19

2.32

20.19

test-002

7.00

0.47

0.91

1.76

19.55

test-003

8.41

2.79

1.90

3.81

35.57

test-004

12.47

7.52

2.40

5.00

38.36

test-005

9.90

4.17

1.83

3.56

28.22

test-006

9.85

0.16

1.75

3.95

37.86

test-007

6.13

3.94

1.02

2.08

20.60

test-008

6.01

1.81

0.97

1.93

17.21

test-009

5.97

1.96

0.78

1.76

25.45

test-0010

11.91

1.63

1.82

3.24

24.22

All (average)

8.55

1.12

1.46

2.94

26.72

Continued on Next Page. . .

43

Volume Name

Volumetric Overlap

Relative absolute

Error

Average symmetric

RMS symmetric

Maximum symmetric

volume surface surface surface difference distance distance distance Table 4.4: MICCAI test sets results as reported by MICCAI grand challenge organizers.

44

Chapter 5 Conclusions This chapter concludes the thesis. Section 5.1 summarizes the thesis’s background, goals, and contributions. Section 5.2 discusses ideas for improvements and future work.

5.1

Summary

Liver segmentation and volume estimation from abdominal CT scans is a key task in many clinical applications. They include the assessment of hepatomegaly and liver cirrhosis, liver regeneration after hepatectomy, hepatic transplantation planning for size compatibility, and monitoring of patients with liver metastases, where the disease is related to an enlargement of the liver, among others. To be of practical clinical use, the liver segmentation and volume estimation must be accurate, robust, fast, and nearly automatic, so that the radiologist can perform it without the help of a technician. Nearly automatic liver segmentation, necessary for hepatic volumetry, is know to be a very challenging task due to the ambiguity of liver boundaries, the complexity of the liver surface, and the presence of surrounding organs. Many previous works propose automatic or semi-automatic methods for liver segmentation. Up until now there’s no automatic algorithm for segmenting the liver that can actually be used by hospitals. Most of the algorithms that are being used by hospitals are semi-automatic and they require many user interaction (such

45

as taking seed points and parameters adjustments). Most of the previous methods for segmenting the liver are using statistical algorithms. The main drawbacks of such methods are that they require a statistical research. They are very sensitive to the process by which the statistics was collected. If some special case will be presented to those methods they will likely fail. Expanding those methods to other organs require additional research per organ. The presented method suggests an interesting concept - in order to segment the liver we are not only define the liver properties itself, we also define its background properties. This concept is not used by any of the state of the art algorithms we met in the SLIVER competition. Very important advantage of our method is that it segments the liver without intensively learning huge data sets. The method fits itself automatically to a new case even if it’s special. We use only one user-selected seed which makes the method almost (or nearly) automatic. The segmentation takes about 6 minutes on average. We validated our method on 2 datasets. The first dataset was taken from the radiology department at Hadassah Ein Kerem hospital from Dr. Sosna’s lab. The second dataset was taken from the SLIVER07 competition. The results from the two datasets totaling 56 CTA images show that our method is accurate, efficient, and robust to seed selection when compared to manually generated ground truth segmentation. We conclude that the presented method takes a significant step beyond previous alternatives towards the goal of automatic liver segmentation method that can actually be in use by hospitals.

5.2

Future work

The presented method was already used as a basis to an additional research for liver analysis (mainly for liver tumors segmentation). In the future we intend to expand the usage of the method to complete the liver analysis, liver contours, liver tumors contours, liver vessels and liver parts. We intend to write a software for visualization and quantitative analysis of the liver, supporting diagnosis and operation planning to be used in hospitals. In addition to liver analysis we intend to activate the method on other organs 46

such as spleen and kidneys.

47

Bibliography [1] http://www.healthsystem.virginia.edu/uvahealth/adult liver/liver.cfm/. [2] Schenk A, Prause GPM, and Peitgen HO. Efficient semiautomatic segmentation of 3d objects in medical images. [3] Black, G. Sapiro, H. Marimont, and D. Heeger. Anisotropic diffusion. IEEE Transactions on Image Processing, 1998. [4] Joachim M menation.

Buhhmann. available

Learning issues and image segat: http://www.stochastik.math.uni-

goettingen.de/files/kolloquium/buhmann.pdf. [5] M. Freiman, O. Eliassaf, Y. Taieb, L. Joskowicz, and J. Sosna. A Bayesian Approach for Liver Analysis: Algorithm and Validation Study. In Proceedings of the 11th international conference on Medical Image Computing and Computer-Assisted Intervention-Part I, pages 85–92. Springer, 2008. [6] Moti Freiman, Ofer Eliassaf, Yoav Taieb, Leo Joskowicz, Yussef Azraq, and Jacob Sosna. An iterative bayesian approach for nearly automatic liver segmentation: algorithm and validation. Int. Journal for Computer Assisted Radiology and Surgery (IJCARS), 3(5):439–446, November 2008. [7] HIPR. http://homepages.inf.ed.ac.uk/rbf/hipr2/threshld.htm. [8] Melonakos J, Krishnan K, and Tannenbaum. An itk filter for bayesian segmentation: itkbayesianclassifierimagefilter. Insight Journal, 2006. [9] Sosna J, Berman P, Azraq Y, and Libson E. Automated liver segmentation and volume calculation from mdct using a bayesian likelihood maximization 48

technique: Comparison with manual tracing technique. In 92nd Scientific Assembly and Annual Meeting of the Radiological Society of North America, 2006. [10] M. Kaas, A. Witkin, and D. Terzopoulos. Active contour models. Journal of Computer Vision, 1998. [11] J. Koenderink. The structure of images. Biological Cybernetics, 1984. [12] Gao L, Heath DG, Kuszyk BS, and Fishman EK. Automatic liver segmentation technique for three-dimensional visualization of ct data. Radiology, 1996. [13] Ibanez L, Schroeder W, Ng L, and Cates J. The itk software guide. 2005. available at: http://www.itk.org/ItkSoftwareGuide.pdf. [14] J. B. MacQueen. Some methods for classification and analysis of multivariate observations. [15] Segmentation of the Liver Competition. http://sliver07.isi.uu.nl/. [16] Teo P, Sapiro G, and Wandell B. Anisotropic smoothing of posterior probabilities. In Proc. of the 1997 Int. Conf. on Image Processing, 1997. [17] Malik J Perona P. Scale-space and edge detection using anisotropic diffusion. IEEE transactions on Pattern Analysis and Machine Intelligence, 1990. [18] Pohle R and Toennies K. Segmentation of medical images using adaptive region growing. in sonka m, hanson km. SPIE Medical Imaging, 2001. [19] Amnon Shashua. Bayesian decision theory lecture notes. available at: http://www.cs.huji.ac.il/ iml/handouts 2008 Amnon/class2-Bayes.pdf. [20] Amnon Shashua.

Maximum likelihood lecture notes.

at: http://www.cs.huji.ac.il/ MaxEnt.pdf.

available

iml/handouts 2008 Amnon/class3-ML-

49

[21] Amnon Shashua. Spectral analysis - clustering. available at: http://www.cs.huji.ac.il/ iml/handouts 2008 Amnon/iml2006-09clustering.pdf. [22] Jianbo Shi and Jitendra Malik. Normalized cuts and image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22:888– 905, 1997. [23] Lambrou T, Linney A, and Todd-Pokropek A. Wavelet analysis of the liver from ct datasets. [24] Okada T, Shimada R, Sato Y, Hori M, Yokota K, Nakamoto M, Chen Y, Nakamura H, and Tamura S. Automated segmentation of the liver from 3d ct images using probabilistic atlas and multi-level statistical shape model. [25] Bram van Ginneken, Tobias Heimann, and Martin Styner. 3d segmentation in the clinic: A grand challenge. 2007. available at: http://sliver07.isi.uu.nl/. [26] Wikipedia. http://en.wikipedia.org/wiki/region growing. [27] Wikipedia. http://en.wikipedia.org/wiki/segmentation (image processing). [28] O. Wirjadi. Survey of 3d image segmentation methods. Fraunhofer-Institut Technound Wirtschaftsmathematik ITWM, 2007. [29] Nakayama Y, Li Q, Katsuragawa S, Ikeda R, Hiai Y, Awai K, Kusunoki S, Yamashita Y, Okajima H, Inomata Y, and Doi K. Automated hepatic volumetry for living related liver transplantation at multisection ct. Radiology, 2006. [30] Nakayama Y, Li Q, Katsuragawa S, Ikeda R, Hiai Y, Awai K, Kusunoki S, Yamashita Y, Okajima H, Inomata Y, and Doi K. Automated hepatic volumetry for living related liver transplantation at multisection ct. Radiology, 2006.

50

Suggest Documents