FULLY AUTOMATED CYST SEGMENTATION IN ULTRASOUND IMAGES OF KIDNEY Abouzar Eslami Department of Electrical Eng. Sharif universty of tech. Azadi St., Tehran, Iran email: [email protected]

Mehran Jahed Department of Electrical Eng. Sharif universty of tech. Azadi St., Tehran, Iran email: [email protected]

ABSTRACT Cyst is one of the most common lesions in kidney and ultrasound imaging is appropriate tool for detecting these lesions. This study develops an automated approach for cyst segmentation in kidney’s ultrasound images. The approach includes three steps: initially, ultrasound image is transformed under a special function derived from Gibbs joint probability function. This transform suppresses noise and discriminates cyst and other tissues. Next, transformed image is decomposed to its low resolution component. Segmentation, morphological operations and coarse boundary detection (performed in low resolution) determines the initial contour employed in the final step. In last step, precise edge detection is performed in unity resolution using active contours model. Proposed approach is designed such that it overcomes noise, imaging artifacts and handles multi cyst cases. Coarse segmentation and then fine boundary extraction is an efficient scheme since segmentation is performed in low resolution (where SNR is relatively high) and lesion boundary is extracted precisely in high resolution (where details available).

hereditary cysts can affect kidney development [3],[4]. There are different proposed approaches for automatic segmentation. In the case of breast lesions [1], Madabhushi et. al. employed seed point extraction and region growing. Although this technique is a classic scheme but it is impractical in presence of imaging artifacts and multi cyst cases. The seed point is defined as the global point which has most probability to belong to the object, while artifacts resemble cysts in histogram and texture and they can have even bigger cystis probability than actual cysts. Also in the case of multi cysts, the approach just extracts one of the cysts. Presence of artifacts aborts employing level–sets and geometric deformable models unless extra operations be employed for discriminating artifacts and cysts. Shen et. al. in segmentation of prostate [2], premised that the prostate has a definite location in image. This is true since their images were taken through rectum. In these images the prostate is located in front of probe and mostly in the center of image. So in their approach there is a constant initial contour (identical for all images) and segmentation is performed during active shape model deformation. In the current study, B-scan images are taken from abdomen and as it will be shown, presence of cyst is probable over a vast region of image. Therefore it is impossible to define a predefined initial contour. Although active shape model is common in medical image segmentation (specially in MR images taken from brain) but since cyst doesn’t have a particular shape, this technique is quite inappropriate in the case of renal cysts. In this work, initially a transform is applied to the ultrasound image. This transform assigns to each pixel a probability value which determines how much it is probable to belong to a cyst. This transform is derived from Gibbs joint probability equation and employs pixel’s intensity, coordinates and texture characteristics. Transformed image is decomposed to low resolution component. In this resolution a course segmentation extracts cysts from other tissues. Exploiting geometric features of objects, artifacts are eliminated and remaining objects (cysts) form the initial contour needed for snake–based boundary detection. Since boundary of each object is determined in unity resolution in the absence of other objects and the approach doesn’t use maximum probability for extracting the seed point, it does not encounter problems in regard to multi objects (cyst).

KEY WORDS Snake, Modified GGVF, Wavelet, Gibbs probability function.

1 Introduction and Related Works Employing digital image processing techniques in medical image understanding has been increased recently. Computer aided diagnosis (CAD) refers to the usage of computers in helping doctors to recognize the abnormal areas in medical images [1]. Although renal cysts (fluid filled objects) are extremely common in adults but they are not mentioned adequately in CAD yet. Different edge and boundary detectors are applied on ultrasound images, while automatic segmentation of lesions in ultrasound, except for few applications (e.g. automated breast lesion segmentation [1] and automated prostate segmentation [2]), is seldom studied. One of these overlooked areas is cyst segmentation in ultrasound images of kidney while simple renal cysts are commonly present in

 of subjects over the age of 50. Approximately  of patients treated by long term hemodialysis develop renal cysts and even some kinds of

458-098

Minoo Naroienejad Department of Radiology Iran Medical Science University Hemmat highway, Tehran, Iran email: [email protected]

127

0.1 0.1 0.1 0.1 0.1

Furthermore exploiting geometric features makes the approach powerful in eliminating the artifacts. The boundary extracted in low resolution is just the initial contour and the total boundary is extracted by snake technique, so the total boundary has the same resolution as image. Proposed algorithm is applied to an actual population of patients and experimental results are presented in the final sections of the manuscript.




   

0.1 0.1 0.1 0.1 0.1

0.1 0.1 1 0.1 0.1

0.1 0.1 1 0.1 0.1

0.1 0.1 1 0.1 0.1

0.1 0.1 1 0.1 0.1

0.1 0.1 1 0.1 0.1

(b)

image. Since convex probes are used in image capturing, polar representation of locations ( ) is more useful. Generally       and    . The upper and lower limits can vary from one apparatus to another but are identical in different images of one apparatus. Considering the training cases, the following function is proposed as a part of the zero order clique potential.

  



  



   

   

(4)

   
        

2.2 intensity This part of the zero order clique potential rises from the fact that cysts have lower density than normal tissues [4]. Our experimental findings imply that the probability distribution of the cyst’s intensity can be discriminated from that of normal tissue’s intensity, with acceptable accuracy. Equation (5) determines the role of intensity in the zero order potential.

(1)

  

  

Assuming that the random field is Gibbs random field (a simplifying assumption, premised frequently in ultrasound image processing both in segmentation [1] and speckle reduction [5], [6]), the joint probability function is defined as:

  N    
  N (2)




0.1 0.1 0.1 0.1 0.1

Figure 1. Line detection kernels; (a) vertical, ker  ; (b) horizontal, ker .

Major parts of our proposed approach (segmentation and boundary extraction) are performed in transform domain instead of gray scale domain. In the first step of our approach, the image is transformed by a function derived from Gibbs joint probability function. Actually, transformed image is a probability map in which the value of each pixel represents its probability to belong to a cyst. Although we do not seek the seed point (pixel associating with maximum cystis probability), but processing this probability map has interesting features in comparison with the crude image manipulation. As the transformed image has less noise and better contrast therefore cysts become more recognizable, fig. 2. Assume that, the undergoing image has the size of
  , therefore N (the set of known random variables) is determined by N  
  
  
  
 . Where  is the intensity of pixel  and   is its location. The cystis probability is a function of N defined by the conditional probability of:

  N 

1 1 1 1 1

(a)

2 Transform Domain

  
       N    N

0.1 0.1 0.1 0.1 0.1

(5)

Where () and ( ) denote the mean and the variance of undergoing image and 
is a statistical value, estimated using training data.

2.3 texture

(3)

Here four different features are extracted to represent texture. Two of these features are the vertical and horizontal lines in a    neighborhood and are calculated using 2D correlation of neighborhood with vertical and horizontal kernels shown in figure 1. The other two are the expectation value of the difference between the central pixel and those on the perimeter of two different size windows. Suppose N  is the rectangular window of size  around (then N  has 25 elements), also N   and N  are the pixels locating on the perimeter of    and    rectangular windows, respectively (certainly around ). So N   has 24 and N  , has 32 elements. By this notation following potentials are proposed as higher order potentials:

Where k is a normalizing constant,   is the potential associating with clique   (each clique is a subset of random variables) and  , its associating weight. It can be shown that, non-negative exponent  guarantees that Gibbs random field is Markov as well and therefore    can be composed just by utilizing a neighborhood around the pixel instead of all the pixels. More details on GRF can be found in [7]. There are three ingredients involved in construction of energy function :

2.1 position Existence of skin, fat tissue, and glandular tissue, directly below the probe, and dept of kidney, causes cystis probability, not to have a unique distribution over the whole of

128

  N   ker
   

(6)

  N   ker
   

(7)

(a)

these images has frequency content of initial image in special region of frequency space and can be down sampled by 2 to reach in half resolution components of  . Iterating this decomposition through a dyadic wavelet transform for 3 times with spline biorthogonal mother functions 6 and 8 leads to
 resolution approximation which has sufficient signal to noise ratio for reliable segmentation. In proposed approach, symmetry is an essential characteristic for filters, because it makes them linear phase and prevents dephasing during analysis stages. It is vital to know the amount of delay caused by filtering in each level to determine the boundary in unity resolution associating with the boundary in
 resolution. Since it is not feasible to have both orthogonality and symmetry we employed biorthogonal wavelet functions [9]. More details on wavelet transform and multi scale processing can be found in [8],[10].

(b)

Figure 2. (a) ultrasound image of a kidney cyst; (b) associating cystis probability map.

(a)

The first process performed on the low resolution approximation is hard thresholding. Notice that since we used probability map instead of crude image and also as we decomposed it to low resolution component, there is no need to employ complicated segmentation approaches. The threshold value is set to    
 experimentally (approximated from training set), where  and  are mean and variance of transformed image’s histogram. Consequent binary image is filtered by horizontal    and vertical    median filters. This eliminates objects whose length or width is equal to one. After that, morphological operations (closing and opening) eliminate holes and scratches of the objects and gives them smooth corners. Regarding to the anatomic definition of cyst (which defines them circular objects), two different size disks are used as structure element. Fig. 4 shows the effect of applying size constraint (median filtering) and morphological operations in eliminating undesirable patterns. After modification of deformed objects, number of objects are counted and then each object is treated individually. Therefore whatever is the number of objects, the approach can extract them and determine their boundary. Here no assumption is made on the number of objects while in seed point–based algorithms the number of objects is assumed as one. Furthermore, in proposed approach geometric characteristics of each object can be extracted which is an important tool for eliminating artifacts. Although the objects such as some imaging artifacts resemble cysts in texture, but they have different shapes and their shape model presents efficient discriminating features. Based on medical definition of cysts, these object are usually circular and have a simple shape while the others do not have circular shapes and also their basic direction tends to vertical axis. so discriminating features can be constructed exploiting principle component analysis (PCA).

(b)

Figure 3. ambiguity reduction in canny edge map of the image in fig. 2 applied to (a) unity resolution, (b) associating approximation in
 resolution.

   N 
  

(8)

   N 
  

(9)

where  denotes the correlation operator. Exploiting these texture features gives the probability map, low-frequency inherent and reduces the spike-like noise (notice the large size windows employed). Such characteristic is extremely desirable since we deal with a noisy image (ultrasound images have high amount of noise and speckle, and we just exerted a    median filter for noise reduction).

3 Low Resolution Processes Ultrasound images have poor quality compared to MRI, CT and X-Ray images. Although exploiting probability map reduces speckles but it is still inappropriate for segmentation. Segmentation in unity resolution can lead to intersection and splitting the cyst to two or more objects. In order to reduce the risk of segmentation, the low resolution approximation is used. Indeed most of the processes are performed in low resolution. Decomposing to low resolution has desirable effect on enhancing the signal to noise ratio. Fig. 3 shows the result of canny edge detector applied to the transformed image (shown in fig. 2) with unity resolution and its associating
 resolution approximation. As it can be inferred, employing low resolution approximation reduces the ambiguity resulted from noise and artifacts. Based on wavelet theory, image  can be decomposed to 4 different images        . Each one of

Let 
and  , (
  ) be the eigen values, V
and V be the eigen vectors associating with covariance matrix calculated from object pixels’ coordinates. An object is not classified as a cyst if both of the following conditions be

129

  that moves through the special domain of image to minimize the energy function 

(a)

(b)

(c)

(d)




 
V
  

"x¼    # x¼¼       x $ (11)


where " and # are weighting parameters that control the snakes tension and rigidity, and x ¼   and x¼¼   denote the first and second derivatives of x  with respect to . These terms construct the internal energy. The external energy function   is derived from the image so that it takes its smaller values at the features of interest, such as boundaries. Given a gray-level image  !, viewed as a function of continuous position variables  !, a typical external energy designed to lead an active contour toward step edges is defined as the gradient of gaussian of image. In order to have bigger capture region and also to overcome concave boundaries, Xu et al proposed generalized gradient vector flow (GGVF). Here, since transformed images include speckles and edges with different strength in the edge map, a number of modifications are performed on conventional GGVF. In the modified GGVF, the external force field is replaced by vector field v     resulted from minimizing the equation (12).

Figure 4. low resolution processes (a) probability map, (b) binary image (result of hard thresholding), (c) result of median filtering and morphological operation, (d) the boundary of object classified as cyst.

met:





(10)



%         v
J $$!

(12) where  is an edge map associating with undergoing image (derived from image) exploiting directional edge detectors such as Sobel and Prewitt. J is a vector field resulted from  whose shrinkage forces are deleted and the weights % and  are proposed as:

The objects who satisfy the above conditions will be dismissed and the boundary of the others will be extracted, exploiting canny edge detector (fig. 4). This boundary is used for precise boundary extraction with snake.

"   #    &      ' (13)   
%  (14)

%   


4 Snake Snakes, since their introduction by Kass et. al. [11], have been widely employed in different applications, including: edge detection, segmentation and motion tracking. Boundary detection in this approach is the result of an energy function minimization. The movement of snake is under the influence of internal forces derived from contour curvature and external forces derived from image. Deformable models are broadly classified as either parametric deformable models or geometric deformable models according to their representation and implementation [13]. Topological flexibility has long been claimed as the major advantage of geometric deformable models over parametric deformable models. Such flexibility is desirable in some applications But it is not always demanded. For instance in the field of medical image processing, object shape is known and is usually used for coarse segmentation in lower resolution. Here high level boundary detectors as active contours, are used just for fine boundary detection and they are not allowed to change the object topology. So in spite of longer history, parametric active contours and snakes are employed here. A traditional snake is a curve x     ! , 

The basic idea rises from the effect of speckle in gradient vector field. The vectors adjacent to speckle turns toward it and so a rotation happens in vector field therefore the curl of associating vector field has big values in the location of speckles. Employing this footprint of speckles, modified weights (% and ) are proposed to control the evolution of gradient vector field. Finally, the snake evolves from its initial location (determined in previous chapter) to the cyst boundary under the influence of conventional internal force with #   and new vector field. Details of implementation of snake evolution can be found in [12].

5 Experimental Results Proposed approach (fig. ??) is applied to ultrasound images of renal cyst, captured by ESAOTE technos mp . In this experiment, 24 healthy cases and 19 patients including females and males aged from 45 to 80, are involved. Of these 43 cases, 17 are used for training and others for test. Images

130

Figure 6. The boundaries delineated by specialist and resulted from the proposed approach.

Table 1. Boundary errors. column 1 and 2: results of proposed approach. column 3: results of multiscale radial segmentation, quoted from [14]. column 4: results of region growing and conventional snake quoted from [1].

Figure 5. The flowchart of proposed procedure.

include cases with more than one cyst and cases with imaging artifacts. The accuracy of proposed approach should be assessed in two fields:



ME FP  FN  TP 

Object extraction and discriminating cyst and similar objects (artifacts and ...).

, 




,- 








(16)

 

(17)

* in [1] 6.687 20.85 24.96 75.04

(18)

where  , + are the boundary and area of object delineated by specialist and   , + denotes the boundary and area of the object resulted by the approach. The two most right columns of table 1 pertains to two other approaches for automatic segmentation, both using seed point extraction for determining the cyst location. One uses coarse to fine radial (multi scale) edge detection for boundary detection [14] and the other utilizes snake approach for boundary detection with initial contour determined by region growing [1]. Although our proposed method has more complexity in comparison, it has better quantified results and even more reliability, desirable features in medical imaging applications. It can also cope with multi cyst cases as well as those including artifacts. On the other hand the referenced approaches provide undesirable results in similar circumstances.

6 Conclusion This study proposed and provided a practical framework for a robust automated segmentation scheme for ultrasound images of renal cysts. As the first step, image is transformed to restrict speckles and to provide better contrast between cysts and other tissues. Since the proposed transform employs texture information it has desirable effects on noise reduction. Hard thresholding in low resolutions leads to better and more reliable segmentation. Each object in consequent binary image is treated individually which

(15)

 



  

The approach was successful in object extraction and discriminating cyst for all the cases. Presented images in fig. 7 show the boundary extraction accuracy of proposed approach. The first image, is an ultrasound image containing renal cyst, imaging artifact and a dark object similar to cysts and other image is a case with two renal cysts. Table 1presents three quantitative benchmarks commonly used for assessment of the boundary extraction accuracy. Benchmarks include: distance error metrics, mean of error (ME) and two surface error metrics, false positive error (FP, and false negative error (FN). ME is the mean of the distance between two boundaries resulted by approach and delineated by expert, FN is the surface assigned as cyst by expert but omitted by the approach and FP is the surface assigned as cyst by approach while it is not determined by expert. The last benchmark, true positive (TP), represents the surface both delineated by expert and extracted by approach. In concept, ME refers to 1-D error while FP and FN refer to 2-D errors. To reach normalized benchmarks, FP, FN TP are divided by the area of cyst. The best experienced result is ME=2.39, FP=, FN=  and TP=
 and belong to a male patient aged 42.

 

* in [14] 5.860 11.11 31.25 76.12

 . 

Boundary extraction.

(   )* +
+

* 5.916 13.15 19.28 80.71

Best 2.39 5.58 4.18 95.82

131

[5] J. H. Hokland and P. A. Kelly, Markov Model of Specular and Diffuse Scattering in Restoration of Medical Ultrasound Images, IEEE transactions on ultrasonics,feroelectrcs and frequency control, VOL. 43, NO. 4 july 1996. [6] L. Onural, M. Bilge Alp and M. I. Gurelli, Gibbs random field model based weight selection for the 2-D adaptive weighted median filter, IEEE Trans. Pattern Anal. Machine Intell. VOL. 16, NO. 8, Aug. 1994. (a)

(b)

[7] S. Geman and D. Geman, Stochastic Relaxation, Gibbs Distribution and the Baysian Restoration of images, IEEE Trans. Pattern Anal. Machine Intell. PMAI-6(6), Nov. 1984. [8] I. Daubechies, Ten lectures on wavelets, CBMS-NSF conference series in applied mathematics. SIAM Ed, 1992.

(c)

[9] A. Cohen, Ondelettes, analyses multirsolution et traitement numrique du signal, Ph. D. Thesis, University of Paris IX , DAUPHINE, 1992.

(d)

[10] A. N. Akansu and R. A. Haddad Multiresolution Signal Decomposision, Transforms, Subbands Wavelets, second ed. 2001, Academic Press.

Figure 7. Experimental results (right column) associating with different ultrasound images (left column); (a) renal cyst (bottom left) with an imaging artifact (top left) and a cyst resembling object (top right); (c) case with two renal cysts

[11] M. Kass,A. Witkin and D. Terzopoulos, Snakes: Active Contour Models Int. J. comput. Vis. , Vol. 1, pp. 321-331, 1987. [12] S. Lobregt and M. A. Viergeve, A Discrete Dynamic Contour Model, IEEE trans. medical imaging 14(1) pp. 12-24, 1995.

gives the scheme ability to handle multi cyst cases. Inspecting the geometric features of each object is useful for classification and discriminating the cysts from other objects. Processing in low resolution proved to be valuable in discriminating cysts from closely resembling objects. Finally, active contour model (snake) with modified GGVF provided a useful tool in extracting precise boundary detection.

[13] C. Xu, and J. L. Prince, Snakes, Shapes, and Gradient Vector Flow, IEEE trans. Image Process. 7(3) pp. 359369, 1998. [14] A. Eslami, S. Kasaei, M. Jahed, Radial Multiscale Cyst Segmentation in Ultrasound Images of Kidney procedeing of ISSPIT 2004.

References [1] A. Madabhushi and D. N. Mtaxas, Combing low-,highlevel and empirical domain knowledge for automated segmentation of ultrasonic breast lesions, IEEE transactions on medical imaging, 22(2) FEBRUARY 2003. [2] D. Shen, Y. Zhan and C. Davatzikos, Segmentation of Prostate Boundaries From Ultrasound Images Using Statistical Shape Model, IEEE transactions on medical imaging, 22(4) April 2003. [3] H. M. Pollack and B. L. McClennan clinical urigraphy Second ed. Vol. 2 PP. 1245-1250 W. B saunders company. [4] D. Sutton, Text Book of Radiology and Imaging, (Sixth ed. 2003, Vol. 2, PP. 950-952, Churchill Livingstone).

132