Interactive Liver Tumor Segmentation Using Graph-cuts and Watershed

. We present in this paper an application of minimal surfaces and Markov random ﬁelds to the segmentation of liver tumors. The orig-inality of the work consists in applying these models to the region adjacency graph of a watershed transform. We detail the assumptions and the approximations introduced in these models by using a region graph instead of a pixel graph. This strategy leads to an interactive method used to delineate tumors in 3D CT images. We detail our strategy to achieve relevant segmentations of these structures and compare our results to hand made segmentations done by experienced radiologists. This paper summarizes our participation to the MICCAI 2008 3 workshop called: ”3D segmentation in the clinic : A Grand Challenge II”.


Introduction
Possible applications related to liver tumor segmentation are mainly radiotherapy and surgery planning.In both scenarios, the knowledge of the exact location and volume of the tumors is a key problem.Liver tumors present low contrasted boundaries and exhibit a large variability of shapes, sizes and locations in the liver.Due to these multiple difficulties, automatic or model based approaches seem to be inadequate.A few semi-automatic approaches have already been proposed in the literature.Recent approaches are mainly based on active contours [12], level-sets [14], as well as machine learning [10].These methods have the interesting property to allow an interaction with its user through land marks positioning or interactive refinement of the segmentation.The interactive approach of medical images segmentation is, to our mind, the most reliable method to provide robust results.We propose in this paper an interactive segmentation method based on watershed and graph cuts to extract liver tumors boundaries in 3D CT images.
The studied CT images 4 were acquired on one 64-slice and two 40-slice CT scanners using a standard four-phase contrast enhanced imaging protocol.The resulting images have a slice thickness of 1mm or 1.5mm and an in-plane resolution of 0.6-0.9mm.Depending on several parameters such as the patient size and cardiac rhythm, the images present different contrasts in the liver.Several cases are illustrated in figure 1.The liver segmentation presents some difficulties that have to be taken into account to design a relevant segmentation protocol.Depending on the acquisition time, tumors can be in positive or negative contrast with the liver.Secondly, the tumors boundaries are not well defined and perceptual properties have to be used to define the exact contours of the tumors.

Images segmentation strategy
The first step of our methodology is the manual definition of a sub-volume containing one or more tumors that need to be segmented in order to reduce the computation time of the segmentation algorithm.Our strategy is then entirely based on the use of the region adjacency graph of an unsupervised watershed segmentation [2], as originally proposed by Li et al. [11].From this first low-level segmentation, a region adjacency graph is extracted and used for the following optimization steps.Our protocol is motivated by the following simple observation: the liver presents two kinds of tissues: tumoral and healthy tissues.The classification of liver pixels in one of these two classes provides the tumor boundaries.To achieve this classification we model the liver pixels as a Markov random field and the classification is performed through the maximum a posteriori estimation.This classification step is supervised by user defined markers that specify both tumoral and healthy tissues.The markers are used to locate the tumors and to estimate the grey levels characteristics of these structures.However, the liver pixels classification needs also that the liver boundaries are extracted.This task is realized by computing a minimal surface based on user defined markers.The user has finally to specify normal liver tissues, tumoral tissues, and external tissues surrounding the liver.Note that the liver boundaries extraction is only necessary if the sub-volume considered for the segmentation contains also nonliver tissues.In any other cases, the liver boundaries are not computed, and only the classification step is performed.The different energy minimization strategies (minimal surfaces and Markov random field) are based on the computation of a minimal graph cut using Boykov and Kolmogorov algorithm described in [5] and implemented in the Boost Graph Library5 .

The watershed transform
An unsupervised watershed transform of the morphological gradient of the original CT image is used in our work to produce a region adjacency graph.Both minima of the gradient and the watershed transform are computed using the 6-neighborhood adjacency system.The watershed transform [2] allows to obtain a partition of an image composed of small and numerous homogeneous regions.Moreover important contours of the image are preserved during the segmentation and regions of the partition are mostly composed of homogenous pixels (pixels of similar grey values).The quality of this first unsupervised segmentation is important to guarantee a minimal loss of information, an ideal situation would be that important information (contours and/or homogenous regions) about the original image is accessible from this first segmentation.These observations about the watershed transform are not theoretically guaranteed but are verified when working on natural images.Another important point is that the watershed transform algorithm based on hierarchical queues exhibits a linear complexity [13].The time needed by the watershed transform is in practice negligible compared to a standard graph cuts algorithm [5].

Approximate minimal surfaces
We detail now how to extract the liver boundaries by an approximate minimal surface using a region adjacency graph [16].The combination of graph-cuts with a watershed low-level segmentation provides us with an explicit and efficient way to compute approximate minimal surfaces.Our basic assumption is that the minimal surface to be computed is embedded in the watershed low-level segmentation contours.We propose thus to solve the following combinatorial problem: finding a surface composed of a finite union of watershed contours such that the surface minimizes a given geometric functional.We solve this problem by using graph-cuts optimization on a region adjacency graph.
Following the formulation of Caselles et al. [7], we want to find a surface S defined by a finite union of watershed contours that minimizes the following energy function: where g is a positive and strictly decreasing function and ∇I(x, y) is the modulus of the gradient of the image I (image contrast) along the surface S.
Note that Cauchy-Crofton formulaes can be used to minimize the energy function E(S) by computing a minimal graph cut as described by Boykov et al. in [4].
Let us consider G = (V, E, W ) as the pixel graph of an image I. Classically V is the set of nodes and represents the pixels of I, E is the set of edges representing neighborhood relations between pixels and W is a positive weight assigned to each edge of E. In our terminology, an edge linking two nodes i and j is written e i,j and the corresponding edge weight is denoted by w i,j .From the pixel graph, we define the region adjacency graph G R = (V R , E R , W R ) of the watershed transform where V R is the set of nodes (i.e the regions of the watershed transform).E R is the set of edges (i.e the neighborhood relation between regions) and W R is the weights of the edges.Let us define F (ri,rj ) as the set of edges of the pixel graph connecting two regions r i and r j of the low-level watershed segmentation: Note that the set F (ri,rj ) depends on the adjacency system of the pixel graph G.The set of edges of the pixel graph describes also implicitly a set of surfaces between the regions r i and r j as illustrated in figure 2. Let S (ri,rj ) denote the set of surfaces that could cross the edges of F (ri,rj ) .Following Cauchy-Crofton formulaes with the V 6 adjacency system, the energy function E(S (ri,rj ) ) can be approximated by: where ∇I(m) and ∇I(n) are the gradient magnitudes of the end points of e m,n .In the following, we consider the strictly positive and decreasing function g: The parameter k ∈ R + is a free parameter that can be used as a smoothing term as shown by Allène et al. in [1].In our application this parameter was set to k = 2.The function g works as an edge indicator of the image I and takes a small value if neighbors pixels m and n take different grey values p m and p n .
The energy E(S (ri,rj ) ) of the boundary between two regions is simply obtained by summing the local contrasts along the boundaries between two regions.The edge weights of the region adjacency graph are then set such that the weight of a graph cut equals the energy function of the surface it implicitly defines: The liver boundaries are finally extracted by computing a minimal graph cut of the region adjacency graph with weights given by equation 5.The minimal cut is computed on the region adjacency graph with two additional nodes s and t, respectively connected to the markers of the liver and the markers of the external tissues.In the following, we denote the markers that specify the liver tissues as the set of regions M t and M h , respectively for tumoral and healthy tissues.The markers specifying the tissues surrounding the liver are denoted by M ext .The edge weights of the graph are summarized in table 1. 1. Edge weights for approximate minimal surfaces.

Approximate maximum a posteriori estimation of a Markov random field (MRF)
This section details the segmentation method used to detect the tumors in the liver.We are now going to take into account a second assumption about the watershed transform: the unsupervised watershed transform of a natural image is composed of regions of homogenous grey level intensities.This assumption permits us to model the image to be segmented as a Markov random field [6,3], where each random variable corresponds to the mean value inside a region of the watershed transform [17].
Let us consider the pixel graph G = (V, E, W ) of an image, as well as the corresponding region adjacency graph G R = (V R , E R , W R ) of its watershed transform.The binary image restoration problem [8,6] is classically defined as the labeling X of the nodes V : such that X minimizes: where p i ∈ [0, 1] is the grey level of the pixel i, and x i ∈ {0, 1} is a label that has to be assigned to the pixel i. N i is the set of neighbors of the pixel i and u i,j is a positive function.δ(x i = x j ) is the indicator function: δ(x i = x j ) = 1 if and only if x i = x j and δ(x i = x j ) = 0 otherwise.
This equation describes the classical formulation of the maximum a posteriori estimation of a Markov random field [6].The data term, also called the likelihood function P r(p i |x i ), ensures that dark pixels, p i ≈ 0, will be assigned the label x i = 0; and that bright pixels, p i ≈ 1, will be assigned the label x i = 1.This model is regularized with the prior function u i,j which is typically used to guarantee that the resulting segmentation has smooth boundaries.We detail now how each term of this energy function can be approximated by using the region adjacency graph of the watershed transform.
Since we assume that each region of the watershed transform is composed of region of homogenous intensities, we can approximate the likelihood term by: where |r i | is the number of pixels inside the region r i .
Assuming that the data are corrupted with a white gaussian noise, the likelihood function can finally be written as: P r(µ ri |( ) , where µ ri is the mean gray level of the region r i , and µ (xi=0) and µ (xi=1) are the mean values of the pixels expected to take the values x i = 0 and x i = 1.
In our experiments the values of σ 0 and σ 1 were empirically set to σ 0 = 0.25 and σ 1 = 0.2.In our model, σ 0 represents the grey level variance of the tumoral tissues and σ 1 represents the grey level variance of the healthy tissues.We have experimentally found that the variance of tumoral tissues is slightly higher than the grey level variance of the healthy tissues.On the other side, the values of µ xi are estimated from the user defined markers as: where |M h | and |M t | are respectively the number of regions marked as healthy or tumoral.We recall that the markers that specify the liver tissues are denoted as the set of regions M t and M h , respectively for tumoral and healthy tissues.
The markers specifying the tissues surrounding the liver are denoted by M ext .
On the other side, the prior function depends on the boundaries properties of the labeling x.We assume now that all pixels inside a region of the watershed transform are assigned the same label x i , the prior function does thus only depend of the pixels lying in the boundaries between two regions: where N ri is the set of neighbor regions of the region r i and |F (ri,rj ) | is the number of edges of F (ri,rj ) .
The prior function u i,j that we use in our application is a contrast sensitive function: where n is a free parameter describing the strength of the term β * (µ ri −µ rj ) and takes also into account the local contrast.In the following we set the parameter to n = 4.Note that if n is very large, our contrast sensitive model is equal to the classical Ising model of ferromagnetism [9].The prior function takes into account the contrast between two regions and is equal to a constant in the areas where the contrast is low.This function allows to detect correctly the boundaries of high contrasted tumors, whereas low contrasted boundaries are smoothed such that the surface of the object is minimized.
We have finally to minimize the following energy function: As shown in [6], we can minimize this energy function by computing a minimal graph cut.The minimal cut is computed on the region adjacency graph with two additional nodes s and t, respectively connected to the markers of the tumoral tissues and the markers of the healthy liver tissues.The edge weights of the graph are summarized in table 2.

Edge
Weight for 2. Edge weights for approximate maximum a posteriori estimation of a MRF.

Post-processing
An additional post processing step is also proposed to the user after the extraction of the segmented tumors.This step consists in smoothing the segmentation by using a morphological opening of the object representing the tumor.The opening is computed with a structuring element of size 1 using the V6 adjacency system.This additional step permits to obtain slightly smoother tumors boundaries.

Example
Figure 3 illustrates our segmentation strategy on a single slice of a 3D CT image.The obtained results are in good concordance with the expected results obtained by a hand made segmentation.Alternatively, the user can add or delete markers if he is not satisfied with the computed segmentation.

Graphical user interface
We have developed a graphical user interface dedicated to 3D medical image segmentation and visualization.The software is entirely developed in python and C++ based on the VTK6 and our own image analysis library Morph-M7 .
Our software allows the user to explore a highly detailed view of the data-set for easy interpretation.Visualization is particularly important for segmentation validation purposes.Data sets can be explored through 2D orthogonal cuts of the 3D volume and 3D rendering of the whole image.The user can interactively provide the markers needed for the segmentation algorithms by drawing on slices of the 3D image, as illustrated in figure 4. The user can visualize the segmented image as a set of surfaces, one surface for each object and surimpose it on the original image.Combination of all this visualization methods allows an easy and fast interpretation of the segmentation result.

Results
Table 3 summarizes the evaluation scores of our method on a set of 5 CT images presenting 10 tumors with unknown hand made segmentations.The evaluation scores compare our results with the radiologists segmentations.The important point is that these results have been obtained without the knowledge of the hand made segmentations.The mean surface distance between our segmentations and the references is approximatively one and a half millimeter, which represents 2 to 3 voxels of the studied 3D CT images.The volumetric overlap error shows that approximatively 71 % of our segmentation volume is in perfect match with a hand made segmentation.Some evaluation results show that we have misunderstood some structures that had to be extracted.This problem leads to very low scores (see for instance tumor number 5).However our results are promising, considering the low quality of some images of the dataset.In comparison, the mean total score obtained on the training data set with known segmentations was equal to 88.The computation time of our method depends on the image size and the number of refinements of the segmentation.The first step, which consists in extracting a sub-volume, is typically done in one minute.The markers placement requires basically one to two minutes.In our experiments we have drawn markers in three orthogonal slices centered on the tumor in about two minutes.The computation time of the segmentation algorithms is then relatively fast: about 1 second for the watershed transform, 2 seconds for the liver extraction and 2 other seconds for the tumor extraction (for a typical sub-volume of size 100×100×100, computed on a common personal computer).The time needed for the refinement steps is then mainly due to the markers placement.In our experiments, two to five refinement steps were needed to provide the presented results.The typical time spent for these refinements varied from two to five minutes.The total time needed for the segmentation of a tumor is then approximatively equal to five, up to height minutes.

Conclusion
Our method exhibits promising results for the aimed application.However some open problems still remain.First, the segmentation of multiple tumors in the same liver often requires additional user markers to correctly separate the tumors.The developed method merges the tumors in a single object when different tumors are too close.This problem requires that the user adds markers between the merged tumors.This additional interaction speeds down the segmentation protocol.However the used methods (minimal surfaces and Markov random fields on a region adjacency graph) are fast and can be used interactively.Secondly we did not develop any preprocessing step such as filtering of the images.There is thus still some possible improvements of our methodology.Future work will be concentrated on the development of adapted filters to simplify the segmentation and the classification step of our methodology.
The use of a region adjacency graph offers a good trade-off between speed and precision for the computation of minimal surfaces and maximum a posteriori estimation of a MRF.We want also to point out that a region based approach is potentially richer than a pixel approach since a wide class of geometrical and statistical functionals can be computed on each region of the watershed transform and additional constraints, such as curvature smoothing term, could be added to the energy function to minimize.These approaches were already successfully used for various image segmentation problems for medical and material sciences applications [16,15,18].

Fig. 1 .
Fig. 1.Liver 3D CT images.Tumors are indicated by red arrows and the aorta is indicated by a blue arrow.The brightness of the aorta indicates the quantity of contrast agent present in it.

Fig. 2 .
Fig. 2. (a) A region adjacency graph.(b) The set of nodes of the pixel graph considered to compute boundary properties between two regions, with a V4 adjacency system.(c) A curve crossing the edges of the boundary between two regions r1 and r2.

Fig. 3 .Fig. 4 .
Fig. 3. Liver tumor segmentation.(a) User specified markers.In blue, the tumor markers, in red, the liver markers and in green, the external markers.(b) Results of our segmentation strategy.(c) Radiologist hand made segmentation and the liver contours extracted by our method.

Table 3 .
Results of comparison metrics and scores for all ten test tumors.