Training Material

Image Segmentation

By Vincenzo Positano, PhD

First appeared in: Casciaro S., Distante A., Samset E. (editors) (2005) "Augmented Reality in Surgery", Text Book of the 2005 ARISER Summer School

adapted for web by Eigil Samset

Abstract

Image segmentation, which associates different tissues imaged by an acquisition device to different logical classes, is often a pre-requisite in 3D visualization (i.e. surface rendering) and image registration and fusion. Moreover, image segmentation is a key issue in surgical planning and training, where a model of mechanical properties of organs is introduced. In this paper, we introduce the problem of image segmentation referring to three examples of image processing algorithms that should be representative of the main possible approaches to the problem.

Introduction

Image processing is a large research field that can be classified in three parts following the traditional computer vision approach. Early vision includes image restoration and basic image segmentation, as well as the image formation problem. Mid-Level vision is about advanced image segmentation, image registration and volumetric image visualization. Finally, high-level vision covers image recognition and understanding, linking the image processing field with the machine intelligence one. In the application of image processing to virtual reality in surgery, about all aspects of computer vision are important, but the main topics of interest are perhaps image segmentation, image registration and 3D visualization.

Three-D visualization represents perhaps the most important connection between the VR device and the user involved in the simulation. Volume rendering works directly on image data, but requires more computational power. Surface rendering approach uses instead image segmentation to extract the relevant regions to visualize. Surface rendering is used in many VR applications, because 3D visualization must be performed in real time or in near real time. So that, Image segmentation is often a needed pre-processing step in VR application. Moreover, it is a main key in surgical planning and training. In fact, surgery is heavily dependent on patient data. In surgery planning, surgeons interact with models of patient anatomy. In surgery training, they operate on a model that is built from patient data. This obviously requires that the models be as accurate as possible for the available data. Image segmentation associates different tissues imaged by an acquisition device to different logical classes. Of course, image segmentation often requires as prerequisite image restoration in order to correct for noise and artefacts. As showed in the following, segmentation algorithms must be adapted to the data to be processed in order to achieve the best performances. An incredible number of image segmentation algorithms was developed in the last years; in this paper, we refer to three example of image processing algorithms, that should be representative of the main possible approaches to the problem.

Medical Image Segmentation

Medical image segmentation continues to be a growing field [1,2,3,4], with an impressing number of papers published for year. Issues related to both general and object-specific segmentation continued to be studied, with both boundary finding and intensity-based region growing being considered as alternatives. More recently, deformable models were discovered and then introduced into the field [5]. The use of deformable models has been pioneered by Kass [6] and generalized for tracking medical 3D data [7,8,9,10]. Active contour models or snakes was improved by the Gradient Vector Flow (GVF) concept [11,12], where a new class of external force, named GVF field, that are vector fields derived from images by minimizing an energy functional in a variational framework. On the region growing side, original ideas aims at utilizing the multiparametric data coming from different data sources, as MRI images acquired with different pulse sequences [13]. Basic decision tree (minimum spanning tree) and statistical clustering strategies were employed to combine the information coming from different channels. Another interesting concept applies forms of scale space theory to the segmentation problem, addressing the idea that relevant medical image features show up at a variety of scales, even within the same dataset.

The key difficulty with segmentation of medical images is perhaps the presence of gray scale inhomogeneities. In MRI, for instance, these are caused by radiofrequency (RF) wave attenuation by tissue, inhomogeneities in the static magnetic field, non uniform RF coil transmission or sensitivity, magnetic susceptibility of tissue. These grayscale inhomogeneities are known to change over time and with different acquisition parameters and to appear in images as local changes in tissue mean and variance. Many different approaches to correcting these inhomogeneities have been reported using statistical methods [14,15] and anisotropic smoothing filtering [16,17,18,19]. In particular, Gerig et al. [17] demonstrated the efficient noise reduction and sharpening of object boundaries by applying non linear anisotropic filtering to MR data of the brain. In this paper, we describe first a method for image filtering, the anisotropic filtering approach, often used as a pre-processing step in medical image processing. The anisotropic filter is also an example of a scale-space method. The second section is about gray level-based segmentation, an in particular about fuzzy C-mean algorithm. Finally, we describe the GVF-snake algorithm, one of the most know approach to contour-based segmentation.

Image selective smoothing

Smoothing is often used as a pre-processing step in medical image analysis in order to reduce image noise. However, the smoothing operation has a blurring effect on image region boundaries, reducing the effectiveness of the following segmentation procedure. A possible solution is the use of an anisotropic filter, that is a filter that works in different way on different parts of the image.  

Perona and Malik [16] formulated smoothing as a diffusive process that is suppressed or stopped at boundaries by selecting proper spatial diffusion strengths. In particular, depending on the values assumed by diffusion strength, the filter is able to realize intraregion smoothing in preference to smoothing across boundaries. In other words, the non linear anisotropic diffusion equation is:

                                                    (1)

The diffusion strength is controlled by c(x,t). The vector x represents the spatial coordinate, while the variable t in our discrete implementation corresponds to iteration step n. The function I(x,t) is the image intensity. The diffusion function c(x,t) assumes a constant value for linear isotropic diffusion. In that case the diffused image derives from application of the Laplacian operator to the image isotropically. It was shown [20] that the price for eliminating the noise with linear diffusion is blurring of edges. This causes their detection and localization to be difficult.

In order to preserve edges, the diffusion must be reduced or even blocked when close to a discontinuity. We choose to be a function of  gradient magnitude evaluated on image intensity I(x,t):

                                                                     (2)

The diffusion coefficient c(x,t) monotonical decreases with increasing gradient ÑI. The parameter K is the diffusion constant and it is chosen in order to preserve edge strength at object boundary and to reduce noise contribute. In order to exploit the relation between K parameter and image gradient I, it is preferable to discuss the flow function f(I) defined as the product c ·I. The flow strength is dependent on the relationship between K and I. The maximum flow is produced at image location with I = K. When I is below K, the flow function reduces to zero because in almost homogeneous regions the flow is minimal. For I larger than K the flow function again decreases to zero, halting diffusion at locations of high gradients. Therefore, a proper choice of the diffusion function not only preserves, but also enhances object edges usually situated at high gradient values. Obviously, it requires setting the K parameter correctly. In [17] Gerig et al. demonstrated by simulation procedure the dependence of the edge slope on K and N parameters, and the effectiveness of anisotropic diffusion filtering to reduce noise and to enhance object boundaries.

The Figure 1 shows the effect of different levels of anisotropic filter on a MRI cardiac image. Should be noted that the noise in the image is reduced proportionally to the filtering level but the transitions between different tissues remain sharp.

Fuzzy C-Mean Cluster Segmentation

The Fuzzy c-mean (FCM) approach [21,22] is able to make unsupervised classification of data in a number of clusters, by identifying different tissues in an image without the use of an explicit threshold. The FCM algorithm performs a classification of image data by computing a measure of membership, called fuzzy membership, at each pixel for a specified number of classes. The fuzzy membership function, constrained to be between zero and one, reflects the level of similarity between the image pixel at that location and the prototypical data value or centroid of its class. Thus, a membership value near unity means that the image pixel is close to the centroid for that particular class. FCM is formulated as the minimization of the following objective function with respect to the membership function u and centroids v:

d                                                          (3)

where W represents the pixel location in image domain, q is a parameter greater that one which determines the amount of fuzziness of the classification (typically q=2is used), is the membership value at location j for class k, is the intensity value at j-th location, is the centroid of the class k, and C is the number of classes.

When the above objective function is minimized, the value of approaches one only if the pixel intensity at j-th location is close to the centroids of class k. Similarly, the value of  approaches zero only if the pixel intensity at j-th location is far from the centroids of class k. Also, the pixels with the same intensity value are grouped into the same groups with the same probability.

The minimization of  is based on suitably selecting u and v by using an iterative process through the following equations:

                          (4)

                                         (5)

The algorithm stops when the value of  converges.

The results of the algorithm are the intensity values that characterize the tissue classes (), and the  masks that describe the distribution of the classified tissues along the processed image.

Figure 2 shows an example of the use of FCM method on MRI abdominal images [23]. In this kind of problem, three tissue classes can be identified: background/air signal, fat signal and signal related to other tissues such as muscle, blood, etc. Consequently, we imposed that C=3 in the FCM algorithm. Three masks were extracted from each image, each related to a tissue distribution and to a representative value  for the extracted tissues. Figure 1 shows a typical example of masks extracted from an MR image. In each mask, white pixels represent pixels with a high value of (i.e. values near 1.0), dark pixels are related to low values of (i.e. values near 0.0). In other words, white pixels represent the background/air area in Figure 1.b, non-fat tissues in Figure 1.c, and fat tissues in Figure 1.d. Only background/air ( and fat () maps are used in our methodology. The representative intensity values are  (air), (non-fat tissue) and  (fat tissue).

GVF snake

A deformable model (also known as "snake" is a surface (a curve in the 2D space) that moves through the spatial domain of the image itself. Snakes or active contours can move under the influence of internal forces coming from within the curve itself and external forces computed from the image data. The internal and external forces are defined so that the snake will conform to an object boundary or other desired features within an image. In general, the standard snake algorithm suffers two limitations: firstly, it is very sensitive to the initialisation curve that must be as near as possible to the detecting contour; secondly, active contours have difficulties to progress into boundary concavities. The introduction of a new external force in the active contour model, as described in, overcomes previous drawbacks. It has been defined as Gradient Vector Flow snake (GVF-snake) and it is computed starting from the gradient vectors of a grey-level or binary edge (gradient) map derived from the image. The use of GVF-snake, brings to more stable and accurate solution to contour detection, also with rough curve initialisation. 

From a mathematical point of view, a 3D deformable model can be defined as a surface that moves through the spatial domain of an image to minimize the following energy functional:

                                               (6)

where:  and s e [0,1], a and b are weight parameters that control the mechanical properties of the deformable model, i.e. tension and rigidity respectively, x' and x''’ denote the first and the second derivatives of x(s) with respect to s, and  is the potential associated to the external forces.  is derived from the image so that it takes on its smaller values at the features of interest. If we want the deformable model to be attracted to edge points, the potential should depend on the image gradient: . Thus the curve is under control of two types of forces: the internal forces (the first two terms in eq(6)) and the image force (the potential term). The choice of external force has been shown to have a deep impact on the deformable model behavior in practical situations [10]. The GVF field is a static type of external force. Given a gradient mapdefined on a data volume, then 3D GVF can be defined to be a vector field v(x) that minimizes the following functional:

                                                                                                     (7)

This formulation forces the field to vary slowly in homogeneous regions and to keep v nearly equal to the gradient map when it is large. The m parameter governs the tradeoff between the first and the second term in the integrand. In fig. 3b, a GVF field representation of the MR image in Fig 3.a is shown. As expected, the GVF vectors are pointing out into the boundary cavity, thus confirming the peculiarity of this operator to enter into boundary concavities.

By applying the calculus of variations [24], it can be shown that 3D GVF field can be found by solving the following Euler equations:

(8)

where:  v(x) = (u(x),v(x),q(x)), . In equations (8) we can treat v(x) as function of time as follows:

                           (9)

where v_t  is the first derivative of v with respect to time. Therefore, given a 3D edge map, 3D GVF is computed using (9). Therefore, the external field v(x) resulting from this calculus of variations is used in eq.(6) as potential force E_ext(x), yelding:

    (10)

where  and  are respectively the second and the fourth derivative of x(s) with respect to s. We call the parametric deformable surface solving eq.(8) the 3D GVF snake.

In the figure 4, the GVF snake is applied to the detection of the endocardial border in a MRI image of the human heart [25]. In (a), the original MRI image and the corresponding GVF field is shown. In (b), the original placement of the snake is depicted, on both the MRI image and the GVF map. In  (c), the snake converged to maximum of the GVF field, that corresponds to the epicardial border in the MRI image.

 

Discussion and Conclusions

Some conclusions can be inferred from the examples showed in the paper. First, image segmentation use in explicit or explicit way some modeling of the images involved in the processing. Optimal parameters setting in the formulation of anisotropic filtering and GVF-snake depends from the nature of the image under examination. In the fuzzy C-mean algorithm, the number of tissue classes to be used depends from the image nature as well. Nowadays, the approach to image segmentation explicitly includes the algorithm dependence from data, and the so called model-based segmentation, that incorporates a model of the structure to be identified in the algorithm, is extensively used. Another popular approach is the Bayesian one, where the prior knowledge of the object to be recognized is included as a –priori hypothesis. The use of machine learning approach to learn the optimal algorithm setting from data is also currently explored. Neural networks, Hidden Markov Models and Support kernel machines was used together with other methods in medical image segmentation. All the three algorithms herein examined are iterative, and the number of iteration needed to reach the algorithm convergence is not a-priori known.  This means that these algorithms are difficult to use where real time constraints are introduced, as in virtual reality applications. Generally speaking, complex image processing approaches require a lot of computational power, and are difficult to implement in real time. Possible solutions are the reduction of the data size, by image downsampling in the space or time domain, and the use of peculiar computer architectures, as parallel computers or dedicated hardware. Finally, the performances of image segmentation algorithms strongly depend from the image quality. Consequently, the optimization of the acquisition procedure is the first an perhaps the most important step in order to achieve an effective image segmentation procedure.

References

1. Duncan JS, Ayache N. Medical image analysis: progress over two decades and the challenges ahead. IEEE Transactions on Pattern Analysis and Machine Intelligence; 22(1) 2000;85–106.
2. Masero V, Leon-Rojas JM, Moreno J.Volume reconstruction for health care: a survey of computational methods. Ann N Y Acad Sci. 2002 Dec;980:198-211
3. Pham DL, Xu C, Prince JL.Current methods in medical image segmentation. Annu Rev Biomed Eng. 2000;2:315-37
4. Udupa JK. Three-dimensional visualization and analysis methodologies: a current perspective.Radiographics. 1999 May-Jun;19(3):783-806.
5. Singh A, Goldgof D, Terzopoulos D. Deformable Models in Medical Image Analysis. 1998. IEEE Press.
6. Kass M., Witkin A., Terzopoulos D. Active contours models. Int. J. Computer Vision, vol. 1: pp.321-331, 1987.
7. Ayache, I. Cohen, and I. Herlin. "Medical image tracking". In A. Blake and A. Yuille, eds. Active Vision, Chapt 17, MIT-Press, 1992.
8. E. Bardinet, L. Cohen, and N. Ayache. Tracking medical 3D data with a deformable parametric model". Proc. Eur. Conf. of Computer Vision (ECCV’96), April 14-18, 1996, Cambridge, England.
9. L.D. Cohen. On active contour models and ballons". Computer Vision, Graphics and Image Processing: Image Understanding, vol. 53, n. 2, pp. 211-218, 1991.
10. L.D. Cohen, and I. Cohen. "Finite element methods for active contour models and balloons for 2D and 3D images". IEEE Trans On Pattern Analysis and Machine Intelligence, 15, ppp. 1993.
11. C. Xu and J.L. Prince. "Gradient Vector Flow: A New External Force for Snakes". IEEE Proc. Conf. On Comp. Vis. Patt. Recog. CVPR’97, pp. 66-71, 1997.
12.     C. Xu and J.L. Prince. "Snakes, Shapes, and Gradient Vector Flow". IEEE Trans. On Image Process., pp.  359-369, 1998.
13. M. O'Donnell, J. Goie, and W. Adams. Towards an Automatic Algorithm for NMR imaging II: Initial Segmentation algorithm. Medical Physics,
1986, vol. 13, pp. 293-247.
14. Dawant, Zijdenbos, Margolin. "Correction of intensity variations in MR images for computer-aided tissue classification", IEEE Trans. Med. Imaging, 12(4): 770-781, 1993.
15. Condon et al. "Image non-uniformity in magnetic resonance imaging: its magnitude and methods for its correction". Brit. J. Radiol., 60: 83-87, 1987.
16. P. Perona and J. Malik. "Scale space and edge detection using anisotropic diffusion". IEEE Trans. On Pattern Analysis and Machine Intelligence, vol. 12, n. 7, pp629-639, 1990.
17. G. Gerig, O. Kubler, R. Kikins, and F.A. Jolesz. "Nonlinear anisotropic filtering of MRI data". IEEE Trans. On Medical Imaging, vol. 11, n. 2, pp.221-232, 1992.N.
18. N. Nordstrom. "Biased anisotropic diffusion: a unified generalization and diffusion approach to edge detection". Image Vision Comput., 8(4): pp. 318-327. 1990.
19. M. S. Atkins, B.T. Mackiewich. "Fully automatic segmentation of the brain in MRI". IEEE Trans. on Medical Imaging, 17(1), pp. 98-107, 1998.
20. L. Alvarez, P.L. Lions, and J.M. Morel. "Image selective smoothing and edge detection by nonlinear diffusion. SIAM J. Number. Anal., 29(3): pp. 845-866, 1992.
21. Udupa JK, Samarasekera S. Fuzzy connectedness and object definition: Theory, algorithm and application in image segmentation. Graphical Models and Image Processing, 1996; 58: 246-261
22. Bezdek J, Hall L, Clarke L. Review of MR image segmentation using pattern recognition. Med Phys. 1993;20:1033-1948.
23. Positano V, Gastaldelli A, Sironi AM, Santarelli MF, Lombardi M, Landini L. An accurate and robust method for unsupervised assessment of abdominal fat by MRI. J Magn Reson Imaging. 2004;20(4):684-689.
24. R. Courant, and D. Hilbert. "Methods of mathematical Physics". Vol. 1. Interscience, New York, 1953.
25. Santarelli M.F., Positano V., Michelassi C., Lombardi M., Landini L.. (2003). Automated cardiac MR image segmentation: theory and measurement evaluation. Medical Engineering and Physics, 25(2):149-159.