Region based techniques for segmentation of volumetric histo-pathological images

Computer Methods and Programs in Biomedicine 61 (2000) 23 – 47 www.elsevier.com/locate/cmpb

Region based techniques for segmentation of volumetric histo-pathological images P.S. Umesh Adiga, B.B. Chaudhuri * Computer Vision and Pattern Recognition Unit, Indian Statistical Institute, 203, B.T. Road, Calcutta 700 035, India Received 17 April 1998; received in revised form 15 March 1999; accepted 25 March 1999

Abstract In this article we have presented the application of three region based segmentation techniques namely, seeded volume growing, constrained erosion–dilation techniques and 3-D watershed algorithm. The algorithms are suitably extended to apply on 3-D histo-pathological images. Suitable modifications and extension for each algorithm is done to obtain better segmentation. A quantitative as well as qualitative comparison of the three methods is presented. Modifications to these algorithms for obtaining better results are discussed. The modifications include, (1) design of adaptive similarity measures to control the seeded volume growing and (2) rule-based merging of the over-segmented cells in the case of the 3-D watershed algorithm. Some results and quantitative study is also presented. © 2000 Elsevier Science Ireland Ltd. All rights reserved. Keywords: Histo-pathological; Volume growing; Erosion; Dilation; 3-D watershed

1. Introduction Automatic image analysis and decision making based on histological images is an interesting and non-trivial problem in histo-pathology. The difficulty in analysis of histo-pathological images arises due to closely clustered cells which are touching or overlapping each other. Proper segmentation of these touching cells is the key to automation of histo-pathological image processing. There exists as many segmentation techniques, as there are segmentation problems. * Corresponding author. Tel.: + 91 33 5778086 ext. 2852; fax: + 91 33 5776680. E-mail address: [email protected] (B.B. Chaudhuri)

Considering what the pathologists try to visualize or quantify in an image, one has to design the segmentation techniques. The tissue architecture and image acquisition noise artifacts gives rise to many problems in segmentation. These have to be solved by semi-automatic processing with user fed information. The imaging artifacts such as nonuniform voxel intensity, large background intensity variation, low and uneven gradient magnitude depicting the cell features such as cell boundary, large amount of noise content, etc. has to be corrected by noise removal and enhancement methods. The fluorescent dye properties such as saturation time, etc. may produce artifacts or misrepresents the features in the image and have to be corrected. Apart from these, even under the

0169-2607/00/$ - see front matter © 2000 Elsevier Science Ireland Ltd. All rights reserved. PII: S 0 1 6 9 - 2 6 0 7 ( 9 9 ) 0 0 0 2 6 - 7

24

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

same imaging conditions, the microscope images may have different artifacts and noise levels depending upon the specimen properties. One of the problems with the conventional 2-D imaging of the thin tissue specimen is that the object features in the image may not represent the real feature values. This is because 2-D images can only partially show the tissue architecture and features of the cells in the tissue specimen. Preparation of extremely thin ideal tissue specimens for 2-D imaging is a relatively difficult task. To reduce the error, maximum number of specimens have to be studied from the same patient during each diagnosis session. This is an extremely time consuming task. Most of the wrong measurement of the cell features in the specimen image and the related problems can be reduced by use of volumetric images. Use of thick tissue specimens and the 3-D imaging results in complete and detailed representation of the cells. The study of spatial distribution of the cells, tissue architecture, tumor grading, structural and geometric feature measurements, counting of fluorescence in situ hybridization signals, etc. can be done more accurately by 3-D image analysis. Among recent developments in the visualization of the histopathological images, confocal laser scanning microscopy (CLSM) is one of the most exciting new developments. The confocal microscope can be considered as a 3-D imaging instrument for collecting data from spatial structures, especially biological ones. The CLSM technique with its capacity for optical sectioning, high resolution reconstruction and most importantly 3-D imaging of the fluorescence labeled or reflecting structures, allows accurate and detailed information to be gathered on the architecture, directionality, size and shape of the cells in the tissue specimen. Relatively few works with practical biological results from such images have been presented. The use of quantitative study is even more recent [1 – 3]. The main hurdle in developing completely automatic image analysis tools is the isolation of clustered cells in the histo-pathological images. The low-level edge based methods may not segment the cells completely due to low and uneven gradient magnitude where the cells touch one

another or overlap. Simple region isolation techniques based on similarity of the voxels belonging to the same region fails to provide the accurate segmentation due to fine textured nature of the cytoplasm and near similar image characteristics of the cytoplasm of all the cells present in the image. More adaptive voxel similarity measures have to be incorporated to determine the membership of a particular voxel to its corresponding region. We have tried to solve many of these problems in the region based segmentation techniques described in this paper. The paper is presented in few sections. Section 2 gives a brief background of the related work. Section 3 gives the details of material and image acquisition. Section 3 also includes some of the image enhancement techniques which has to be considered before designing the system. Section 4 explains segmentation methods such as constrained seeded volume growing technique for the segmentation, constrained erosion–dilation technique, the 3-D watershed algorithm applied to volumetric histo-pathological images and an automatic merging technique for reducing the oversegmentation of the cells. Section 5 deals with experimental results and a comparative study of the three region based segmentation techniques.

2. Background A region is defined as a connected set of voxels verifying some homogeneity criterion [4]. Similarity of different properties of the voxels belonging to the same cell holds the key to the success of this algorithm. The voxel properties include, intensity, gradient magnitude, global and local morphological criteria, etc. The most simple and primitive region based technique is thresholding the image based on a local histogram of the voxel intensity and then labeling the isolated regions by connected component labeling algorithm [5,6]. Histogram based thresholding techniques does not exploit the fact that the points from the same object are generally, spatially close due to the surface coherence [7]. Some of the generally established region based segmentation algorithms are split-and-merge algorithm [8], region growing [9],

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

multiple thresholding [10], etc. As stated earlier, due to the inherent characteristics of the 3-D histo-pathological images, direct application of these images need not result in proper segmentation result. There are several variations to the region growing techniques such as seeded region growing [11], successive erosion and dilation [12], or successive pealing and thickening, region growing based on simple distance map, the watershed algorithms [13], etc. In some of the earlier experiments [14,15], the region growing technique were tested successfully on the simple 2-D images where a single cell is present in the image. In such cases the segmentation becomes relatively less complicated. In case of complex 3-D images of a tissue specimen, where several cells are compactly arranged, and the separation between the touching objects is not characterized by proper uniform gradient peaks, more morphological constraints such as size and shape constraints are to be used to obtain optimum results. A more adaptive similarity measure is also necessary to take care of textured cell chromatin. In many cases interactive correction and use of some heuristic methods besides controlled volume growing is necessary for the complete isolation of the cells. In this paper we have tried to present the region based methods which could solve some of these problems.

25

col of Alers [16] was used. Briefly after deparafination and rehydration, the samples were digested in a pronase E solution (0.05% in phosphate buffered saline, PBS, pH 7.0) for 2 h at 37°C, with vigorous vortexing every 15 min. The cell suspension was then centrifuged, washed twice in cold PBS before filtering through a nylon mesh (pore size 30 mm) to remove any remaining tissue fragments. The cells were resuspended in a cold Carony’s fixative (methanol:acetic acid= 3:1) and dropped on to clean slides for the purpose of image acquisition.

3.1. Image acquisition The working principle of the CLSM can be briefly stated as follows. A sharp laser beam is focused on a particular point of the focal plane in a thick tissue specimen. The reflected laser beam being modulated by the specimen property at the respective focal point, carries information about the specimen at that point. The focal point is moved in a lateral plane in a raster pattern for scanning and the information about the specimen in a corresponding focal plane is obtained. Then the focal plane is varied in the axial direction and next image slice is obtained. The lateral and axial resolutions are controlled with the help of a computer interfaced to the CLSM. In this way a stack of optical sections depicting the 3-D information about the tissue specimen is obtained.

3. The material and design considerations We have processed thick specimens of the prostate carcinoma. The specimen preparation, fixing and image acquisition has been done with the help of pathologists, physicians and physicists of Bio-Medical Image Analysis Division, Institute of Pathology, GSF, Munich, Germany. Routinely processed formalin-fixed and paraffin embedded tissue specimens from radical-prostatectomies of several patients with prostatic adenocarcinoma were used to acquire the images. Several 30 mm thick sections were prepared and stained with methylene blue. From these sections the corresponding tumor areas were selectively cut out with a fine scalpel under microscopic control. For desegregation of tumor nuclei the modified proto-

3.1.1. Set-up features for image acquisition Fluorescence images are scanned using a confocal laser scanning microscope (CLSM) Zeiss LSM 410. Essential setup features for the acquisition of FISH signals in tissue sections are as follows: Lens Zeiss PNF 100x, 1.3, zoom =2, realized by scanning unit). The scanned field of 62.5× 62.5 mm is sampled to 256×256 pixels giving a pixel size of 0.25 mm2 in the x and y directions. Excitation laser lines are selected according to the fluorochromes used. For propidium iodide (PI) used as DNA counter stain and fluorescein isothiocyanate (FITC) labeled signals both are excited by the Argon line 488 nm. Lateral resolution of ca. 0.5mm is obtained by the system. The lateral resolution is related to the wavelength of the

26

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

excitation laser as resellat =0.8l/2(NA) where NA is the numerical aperture which is 1.3 in the present case. The emission of FITC and PI fluorescent materials is measured simultaneously in two separate channels using a bandpass 515 – 565 nm for FITC and a lowpass LP 590 for the PI channel, respectively. The axial distance selected between two subsequent confocal images depends on further evaluation strategy. If a spatial isotropic representation of the 3-D data is intended, the axial distance is identical to the lateral distance of the pixels, i.e. 0.25 mm. For a 16 m thick section this results in a sequence of 64 scanned images. If only the sampling theorem should be met, roughly a 0.5 mm distance is enough, noting that the confocal axial resolution of a lens with numerical aperture NA = 1.3 is given by: reselaxi = 1.4lh/(NA)2, where h is the refractive index of the object. Axial resolution is approximately equal to 1 mm. The FITC emission is associated with the green channel and the PI emission with the red channel of a RGB color image. The image data are stored on disk in TIFF format.

3.2. Design considerations Several quantitative distortions can be introduced into 3-D stacks of optical sections during data collection. Once the sources of errors are identified, some can be prevented by modifying the microscope setup, while others are to be corrected either at the time of data acquisition or afterwards. Most of the errors that can be corrected by image processing are, low axial resolution, the uneven illumination, photobleaching effects, out of focus information or blur, various kinds of noise artifacts, etc. Confocal microscope has an inherent property of minimizing the outof-focus information in an image due to its unique pinhole point detector. This helps in doing away with the computationally expensive deconvolution process. Image enhancement is necessary to reduce the textural nature of cell chromatin, reinforce the boundary where the cells touch one another so that smooth and proper region growing can be implemented.

One of the important aspects before image enhancement, is to separate the background and foreground regions. Window slicing which is also known as amplitude thresholding converts all the pixels or voxels which are below a threshold Ts to the lowest gray value that is zero (indicating dark). The voxels having intensity value above the threshold is kept undisturbed. This gives us an image of uniform background while the foreground consists of different gray levels. We have observed that the histogram of our image data sets exhibits unimodal property. Hence the valley point indicating likely threshold, may not be found explicitly in the histogram of the images. In such cases it is often possible to define a good threshold at the shoulder of a histogram [17]. Since both the valleys and the shoulders correspond to the concavity of the histogram, a threshold can be chosen by analyzing the concavity of the histogram [18]. There are several methods to choose the shoulder point of the histogram. We have used an empirical formula, t= (m9 ks) where s is the standard deviation of the gray value and k is a constant determined experimentally. Here, t=(m+ ks), when the mode point of the histogram lies close to the minimum gray value in the image and t= (m−ks), when the mode point lies close to the maximum gray value in the image. Sometimes, window slicing results in the creation of holes in the cells and small noisy island like object structures in the background. This is due to the presence of dense non cellular matters in the background and dark intracellular objects within the cell. We have developed a size and shape heuristic based filter to reduce this effect. The size of all the isolated objects is calculated. If the size is below a predefined size threshold then such objects are considered as artifacts and its voxel intensity are changed to background intensity. Similarly, the holes are identified. If the holes are smaller than the predefined threshold, then the original gray values of the voxels belonging to the holes are restored. Fig. 1b and Fig. 1c show the result of window slicing and size filtering on a representative image slice of CLSM image

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

27

stack. After window-slicing and size filtering of the image stack, the background noise is reduced and the background gray value becomes spatially uniform. Another problem with the confocal image stack is the attenuation of light along the depth of the specimen. The uneven illumination along the depth of the specimen results in the spatial variation of light intensity in the image volume. Besides the optical problems, the photobleaching of the specimen contributes to the degradation of

Fig. 2. Graph of average foreground intensity against the depth of the image stack, (a) before intensity restoration, (b) after intensity restoration.

Fig. 1. Result of image enhancement, (a) a representative image slice of CLSM image stack, (b) after window slicing, (c) after size filtering, (d) after surface crispening, (e) result of spatial averaging, (f) after directional Gaussian smoothing.

intensity. Photobleaching can be modeled as a first order decay process and hence is computationally corrected [19]. When we plot the average image intensity of each image slice against depth of the stack, we have observed that the illumination degradation is not linear. It is not possible to define the variation of image intensity as a linear decay as is evident from Fig. 2a. There are earlier works by Aslund et al. [20], Rigaut et al. [21], etc. which give sufficient restoration of image intensity along the depth of the image stack. We have implemented a simple method of restoration of intensity of the foreground voxels by comparing it with the highest intensity image slice in the stack. Let Ii be the image slice having maximum average image intensity, i.e. Ii = max{I1, I2,…., In } where I1, I2,…., In are the intensities of 1, 2,…., nth image slice in the stack. We consider this image slice i as the standard image slice and increase the average image intensity of remaining image slices in the stack to be on par with average image intensity of image slice i. Increasing the average image intensity of the whole image slice increases the background intensity too, which is undesirable. Hence we have to consider only those voxels that belong to a particular region of interest where the intensity compensation is necessary. The foreground (highlighted

28

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

by locally high gray value or image intensity due to specific fluorescent material) of the image is considered as the region of interest in the present case. Let the mean intensities of the image slices 1, 2,…., n, in the image stack I1, I2,…., In. The variation in the average image intensity of the foreground is plotted as shown in Fig. 2 The maximum average intensity of foreground of the image slices is considered as a reference image slice, i.e. if Ii is the average intensity of reference image slice, such that Ii ]Ij for all j values then, image slice i is considered as the reference image slice and Ii as the standard image intensity value. Let b= Ii − Ik be the difference of average intensity of the foreground in kth image slice and the reference image slices i. Then for the kth image slice the gray level of each pixel of the foreground is enhanced by a factor cbk, i.e. if I(x, y, k) is the intensity of voxel at (x, y, k), then the enhanced intensity is given as I(x, y, k) =I(x, y, k) + cbk where c is an experimentally chosen constant. Ideally c should be 1. This simple addition of the average value to the image intensity results in the loss of pixel sensitivity. A trade off optimizing the requirements of light intensity and loss of sensitivity is useful. This tradeoff is also imaging and application dependent. As the confocal microscopy images do not give clear details of the intracellular structures as well as our interest being limited to measure the quantitative features of cells and the tissue, we have ignored the pixel sensitivity issue. Fig. 3 shows one of the original image slice and after the intensity restoration.

Fig. 3. Result of image intensity restoration shown over a single image slice, (a) original image slice, (b) after intensity restoration.

Fig. 4. Diagrammatic representation of directional Gaussian filter.

The next step in image enhancement is to crispen the cell surface. We have used the conventional process of adding a high pass signal to the image for this purpose. If u(x, y, z) is the input image and g(x, y, z) is the gradient magnitude image, then the crispening process can be written as 6(x, y, z)=u(x, y, z)+lg(x, y, z) where l] 1. We have used l=1 as the scale factor. Fig. 1d shows the result of surface crispening. Surface crispening process also enhances the noisy peaks in the foreground, which is undesirable. One of the simple ways of reducing these noisy peak signal is to subject the image stack to smoothing. Simple spatial averaging blurs the cell boundary while smoothing. This is an undesirable effect. Symmetrical Gaussian smoothing also blurs the edges while smoothing, though the influence of far off voxels on the smoothing is reduced. To reduce the blurring effect on the edges, we have used directional Gaussian filtering. A 5× 5 Gaussian filter is subdivided into six directional windows as shown in Fig. 4 The Gaussian convolved values of the pixels in each directional window is calculated as shown. The maximum of the convolution values in all the directional windows gives the desired result. If 6(x, y, u) are calculated in several directions as,

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

6(x,y,u)=

1 SS u(x,y)G(x −m,y −n) Nu(m,n)  Wu

!

"

(x 2 +y 2) is a Gaussian 2s 2 filter, Nu is the total number of pixels present in the directional window Wu as shown in Fig. 4. A direction u* is found such that u(x, y) − 6(x, y, u*) is minimum. Then the output image 6(x, y) < 6(x, y, u*) is the desired result. Directional Gaussian smoothing technique is applied separately to all the image slices in the image stack. Fig. 1f shows the result of directional Gaussian smoothing in comparison with a simple spatial averaging (Fig. 1e) to reduce the noisy peaks in the image while protecting the edges. One of the important enhancement steps to be carried out is the improvement of axial resolution of the image stack. Poor resolution of the image stack leads to error in segmentation. Due to anisotropy in voxel lattice, the direct 3-D processing algorithms fail to properly make use of spatial neighborhood relations. To avoid these errors, and to enhance the qualitative and quantitative accuracy of visualization and analysis, a suitable interpolation process has to be used to increase the axial resolution of the stack. We have done this by virtually inserting the interpolated image slice in between the actual image slices in the stack. Interpolation is done based on morphing. Consider the two tone version of two neighboring image slices (source images) j and ( j+ 1) where the object is represented by gray level 1 and background by 0 as shown in Fig. 5a. For creating a image slice in between j and ( j+1), we have to choose some control point in the interpolated image slice towards which the source images are distorted for the purpose of morphing. It can be argued that if by some means we can get the overall boundary of the objects in the interpolated image slice, we can use these edge pixels as control points and fill the gray level within these boundary points using morphing or weighted averaging technique. The overall boundary of the objects in the interpolated image slice can be found by using simple logical operations as follows. The two tone version of the neighboring image slices j and ( j+ 1) are subjected to pixel-topixel X-OR operation. The resulting image conwhere G(x,y)= exp

29

sists of the portion of the object which is common to the objects in both j and ( j+ 1) image slices. This is shown diagrammatically in Fig. 5a. The medial axis of this image gives the boundary of the objects in the interpolated image slice. If the objects in j and ( j +1) image slices are laterally shifted or there is a strong variation in shape of the objects between two source images then the resulting medial axis need not be continuous. To join the disconnected boundary, we have applied the following simple operation. The contours of the objects in the two source images are obtained over a two tone version of the source images j and ( j+1). The common portion of the boundary of the objects in j and ( j +1) image slices are obtained by pixel-to-pixel AND operation of the two boundary map of the source images. The result of the AND operation is added to the medial axis of the interpolated image slice by logical OR operation. This results in linking the broken contour. The process is shown diagrammatically in Fig. 5b. To fill the intensity information within the boundary contour of interpolated image slice, the window-sliced version of source images j and ( j+ 1) are distorted towards the position of the contour based control points in the interpolated image slice. Then, the two deformed images are blended with simple weighted averaging to generate the gray value in interpolated image slice. Let Gj = Ij 0, Ij 1, Ij 2,…, Ij n be the gray value of the object pixels in source image j and Gj + 1 = Ij 0 + 1, Ij 1 + 1, Ij 2 + 1,…, Ij n + 1 are the gray value of the object pixels in source image ( j+ 1). Then, the gray value of corresponding object pixels in interpolated image slice is given as (w1G1 + w2G2)/2w where w1, w2 and w are the constants which are calculated experimentally. In our experiment we have used all the weights as 1 that simplifies the method to simple averaging. For the pixels that fall outside the boundary contour in the interpolated image slice, the gray value is given as zero. As this method incorporates the high level knowledge, i.e. contour information, in the interpolation process, it is superior to the other intensity based interpolation processes. Fig. 6 shows the result of interpolation. Fig. 6a and Fig. 6b are the two source images while Fig. 6c is the result of interpolation by morphing.

30

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

Fig. 5. Schematic diagram explaining interpolation by morphing, (a) interpolation when the cells are not laterally shifted or deformed, (b) interpolation when the cells are laterally shifted or deformed.

Since most conspicuous features to human eyes in images occur at the places with higher gradient magnitudes, which are also where the

contour or boundary point usually locates, using contours as control lines for morphing are very effective.

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

31

Fig. 6. Result of interpolation by morphing, (a) source image 1, (b) source image 2, (c) interpolated image slice.

4. Segmentation system The segmentation of histo-pathological images can be considered good if the result shows the following properties. 1. Complete and detailed isolation of the clustered cells. 2. Minimum over-segmentation and under-segmentation of the cells. 3. Each visually perceivable 3-D cell has been given a unique label. 4. The segmentation process should be less interactive and the results should be comparable to manual segmentation. 5. The process should be fast and efficient. We have implemented the modified 3-D versions of volume growing based on voxel similarity properties, Successive pealing and thickening technique and 3-D watershed with a rule based merge technique to reduce over-segmentation.

4.1. Constrained seeded 6olume growing This method is based on the well established region growing techniques for similar gray values inside the object. All the region growing techniques need a starting point from where the region starts to grow. We call these starting points as seeds. Fig. 7, shows the flow diagram of the seeded volume growing technique.

4.1.1. Seed marking Seeds can be marked interactively by clicking the mouse over the cell(s) of interest and considering a small group of connected voxels around the

selected voxel to calculate the seed properties. Seeds can also be chosen randomly all over the image. The number of seeds chosen in random methods must be high enough so that all the cells to be segmented will have at least one seed in it. We have used a more reliable and algorithmic approach for automatic seed marking. 1. A 3× 3× 3 six connected structuring element (say K) with unit magnitude is defined. 2. Erode the two tone version of the noise filtered and enhanced image volume by K. 3. Measure the size of each eroded individual object in the image. If the object size is below a predefined size-threshold, mark that object as seed and no further erosion of the seed is allowed. The similarity measures, calculated on the basis of a few seed voxels may lead to wrong results in case of noisy images. This is because the seed point may fall on a noisy voxel resulting in erroneous seed properties. To avoid this we have set a minimum seed size consisting of few (ca. 35) connected voxels. When an eroded object reaches the minimum size threshold, it will not be eroded further but marked as a seed with a unique label. When few iterations of erosion gives unique seed to every cell, eroding further is not needed. Hence we have put up a criteria that after first three iterations, for every few (two) successive erosions if the number of individual seeds does not change then the algorithm should understand that all the cells are given a unique seed and stops further erosion, i.e. if n(t) is the number of marked seed points at time t then the erosion can be stopped if n(t)− n(t− 2) = 0. The seeds are then labeled with a 3-D component

32

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

labeling algorithm. Fig. 8c shows the marked seeds (marked regions) inside the cells.

4.1.2. Implementation For a voxel to be considered as belonging to the cell of interest, the following properties should hold good.

1. Voxel gray value should be within the similarity measure threshold limits. 2. Voxel gradient magnitude should be within the predefined limits. 3. Global morphological condition such as size, compactness, elongation, etc. should also satisfy the predefined conditions, if any.

Fig. 7. Flow diagram of seeded volume growing technique.

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

33

Fig. 8. Result of seeded volume growing for segmentation of cells, (a) original image slices, (b) result of local thresholding and component labeling, (c) after seeds are superposed on the cells, (d) result of seeded volume growing shown over few image slices.

4.1.3. Similarity measure based on 6oxel gray 6alue The general assumption we made is that the voxels belonging to the same cell are approximately similar in their gray level. If agray is the gray intensity similarity measure, then

agray = m−uq





where u represents the image voxel, 1 m = ni= 1ui is the mean gray value of the marked n region, n is the number of voxels in the already

34

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

marked region and uq is the voxel belonging to the neighborhood of the marked regions which are placed in an ordered queue for testing the similarity. If agray is within a predefined threshold range, then the voxel satisfies the first and a major condition. In case of highly texture images, the similarity measure needs to be normalized as follows. agray =

m − uq s

where s is the standard deviation of the gray level of the marked regions. The standard deviation of the gray values is given by s=





1 ni= 1 (ui − m)2 n− 1

1/2

.

Use of normalized similarity measure considerably slows down the process of volume growing and should be avoided wherever is possible. The edge

Fig. 8. (Continued)

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

preserving smoothing techniques such as directional Gaussian, etc. are the better options.

4.1.4. Similarity measure based on the gradient magnitude This measure is defined based on the assumption that the variation in voxel intensity within a cell is smooth causing similar gradient magnitude and at the surface of the cell the voxel intensity varies abruptly showing a local peak in the gradient map. This measure helps in marking the cell surface there by stopping the seed growth beyond the cell surface into neighboring cells. The gradient similarity measure is defined as agrad =

)



1 n % 9u − ( 9uq ) ni = 1 i

)

where 9u gives the gradient magnitude of the voxel intensity. Presence of intracellular matters and fine textured nature of the cell chromatin adversely effects the assumption of uniform gradient magnitude within the cell. This can be rectified to certain extent by adopting the normalized similarity measure similar to the one explained regarding agray or by locally smoothing the image before calculating the gradient magnitude. The smoothing of the image with suitable edge preserving filter should be preferred to avoid the slowing down of the volume growing process. We have implemented the volume growing method by a simple queuing algorithm. All the seeds are labeled. The immediate spatial neighbors (18 connected) of each seed are put into a corresponding labeled queue. If s1, s2, s3,…, sn are the seeds unique for the cells 1, 2, 3,…., n and q1, q2, q3,….qn are the corresponding queues for each seed, then queues contain those unmarked voxels which have at least one of their 18-connected spatial neighbor as already marked voxel or seed voxel. Till the queues become empty by converting all the members in to marked voxel or discarding them for not being similar, the queues are not updated. For each voxel in the queue, the intensity and gradient similarity is calculated and if they are within the threshold they are added to the corresponding marked region. If the voxels are not similar then they are checked for being a noisy voxel or not. This is done by checking for

35

its isolation property. If it is an isolated voxel, showing a dissimilar intensity and/or gradient, it is considered as a noisy voxel within the region of interest and added into the marked region otherwise it is discarded. When the gradient magnitude of the voxel shows high gradient magnitude, such a voxel is separately marked as a surface voxel. After every three iterations the morphological criteria for size and shape are calculated to check the over growth. Each iteration consists of emptying all the voxels in the queue. Generally cells are convex in shape. In a healthy tissue the size, shape and other geometrical features of the cells show less deviation. So the shape features can be effectively used to control the growth of the seeds. It can be said that the global morphological criteria aggregate the voxels until fitting the cell under a convexity assumption. After every few (three) iterations of volume growing, the size and elongation (Rmax/Rmin) are checked. If any abnormal deviations are found, then the volume growing is stopped. The general assumption that the cells are convex may not hold good always. In the higher grades of carcinomic stages, the cell starts growing arbitrarily, losing their convex shapes. One has to relax the convexity criteria in such cases. When the image volume is of malignant tumor, the assumption of the convexity of the cell as well as the similar size and shape may not hold, leading to major errors in controlling the volume growing methods. In such cases interactive corrections are necessary to get acceptable segmentation. The process of volume growing is stopped when all the voxels in the foreground are marked as separate cell regions separated by the surface voxels or when the morphological criteria fails to hold good. The resulting marked regions are subject to morphological opening [22]. This cuts off the thin and frail connections between the cells caused due to low gradient magnitude where the cells touch one another. Then the cells are labeled using 3-D component labeling algorithms. Fig. 8 shows the result of seeded volume growing as a sequence of image slices constituting the volumetric data set. Fig. 8a shows the original image slices. Fig. 8b shows the result of simple thresholding and labeling the enhanced image.

36

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

Fig. 8c shows the marked seeds inside each cell. Fig. 8d is the result of volume growing. It is evident from the results that the seeded volume growing gives reasonably good result when the cell surface is marked by a uniform gradient peak. Post-processing and interactive corrections are sometime necessary for complete isolation of touching or overlapping cells.

4.2. Constrained erosion– dilation technique This is one of the simple growing methods which yield acceptable results in many cases. This is useful in automatically cutting off the touching and overlapping cells. This is also known as erosion– dilation technique [12]. Segmentation of cells in a 3-D image by successive pealing and thickening technique can be explained in few steps. In the first step, the objects in the foreground of the image are subject to a erosion process. Erosion with a suitable structuring element, is an iterative process of removing the voxels of the foreground which are having at least one background voxel in its immediate neighborhood. The process continues till a unique signature is obtained for each cell present in the foreground. Care should be taken that no cell signature completely disappears from the spatial image domain. This causes the unrecoverable loss of information. To avoid the loss of information, we have set the minimum size of the object that can not be subject to further pealing as approximately 100 voxels similar to the process of seed marking explained in Section 4.1. In the second step, the cell signatures are subject to controlled dilation under the condition that no two dilated signatures would overlap or touch each other. The number of iterations for dilation process should be greater than the number of iterations done during erosion of the cell signatures. In the third step, the dilated image is subject to one-to-one logical AND operation with the two-tone version of the original image volume. This marks the bifurcation between the cells in the cluster.

Fig. 9. Diagrammatic representation of process of constrained erosion – dilation.

In the fourth step, the 3-D component labeling algorithm is applied and all the isolated cells in the image are uniquely labeled. After labeling the segmentation process can be said as complete. Fig. 9 diagrammatically shows the process of constrained erosion–dilation technique for separation of cells in a cluster. Fig. 10 shows the result of this method over a CLSM image data set as a sequence of image slices. Another method of volume growing is based on defining the zone of influence of each seed or regional minima of the individual cells on the basis of distance measure. Then the seeds are extended to all the voxels in its zone of influ-

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

ence. The zone of influence is marked by watershed line. This method is known as watershed algorithm. The classical watershed method is ex-

37

tended to 3-D and we have improved the watershed algorithm by merge technique to reduce over-segmented objects.

Fig. 10. Result of constrained erosion and dilation, (a) original image slice, (b) after erosion, (c) after successive dilation, (d) result after AND operation.

38

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

Fig. 10. (Continued)

4.3. 3 -D watershed Segmentation based on the use of watershed lines has been originally developed in the framework of mathematical morphology by Lantuejoul

[23] and later on improved jointly with Beucher [24]. Considering gray level image domain as analogous to watershed can be explained by mapping the gray scale image as a topographic relief. The gray level of each voxel stands for the eleva-

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

tions of the corresponding watershed surface. Though mapping of gray levels as a topographical relief in 3-D image is bit fuzzy to understand, the algorithm is a simple extension of watershed in to 2-D to 3-D. Watershed algorithm is an iterative region growing process like the earlier two described methods. The region grows only within the zone of influence (ZOI) of each regional minimum or regional marker. The regional minima Ri of a 3-D

39

object i is defined as a connected group of voxels at the approximate center of the object characterized by the following property. If the shortest distance of this group of connected voxels from the object boundary is di, then it is not possible to traverse to another region of distance longer than di without traversing through the voxels of distance value shorter than di. In simple words, a regional minimum is a group of voxels or a single voxel, belonging to the object and has a maximum

Fig. 11. Schematic diagram explaining the concept of geodesic distance, SKIZ, ZOI, regional minima, etc.

40

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

distance value from the nearest background voxel. Fig. 11, shows the concept of regional minima concerning a 3-D object. The most simple way to find the regional minimum is to threshold the distance map such that each cell gets at least one regional marker. The distance map defines the zone of influence of each regional minima. We have used path generated distance transforms (PGDT) as proposed by the Borgefors [25] for finding the zone of influence. The PGDT is briefly explained below. The two-tone version of the image is obtained by thresholding the image at suitable level as explained earlier. All the background voxels are given a distance value 0 while the foreground voxels are given a very high distance value (ideally infinity). In the first scan from the top-left, a new distance value is computed for the foreground voxel. A 3× 3× 3 neighborhood of each foreground voxel is considered along the raster scan. Compute the sums of the values of the already visited neighbors and the corresponding local distances. In a 26 neighborhood there are 13 already visited voxels for each central voxel. Thus the new distance value for the central voxel is the minimum of the 13 sums. This scan computes the distances from the left-up-top. In the second scan, the scan direction is reversed. Again for each voxel, sum of the already visited neighbor distance values and local distances are computed. The only difference is that now the central voxel itself must be included, adding local distance zero. Thus the new, final value is the minimum of the fourteen sums. The second scan computes the distances from right-down-bottom. The geodesic distance di of a voxel i from the regional minima Ri of a 3-D object is the length of a shortest existing path which is included in F and linking voxel i and regional minimum Ri of the object. Fig. 11, presents the concept of geodesic distance. The zone of influence (ZOI) of a regional minima Ri consists of all those voxels whose geodesic distance di from the regional minima Ri is smaller than their geodesic distance to all other regional minima Rj " i where j= 1, 2,....,N N being the total number of regional minima in the image volume. Fig. 11 shows the concept of zone of influence

diagrammatically. The zone of influence of each cell marker (regional minima) approximately covers the voxels belonging to the respective cells. Analogous to the watershed terms, ZOI can be called as catchment basins. The skeleton by zone of influence (SKIZ) is defined as the surface consisting of those voxels in the foreground which do not fall into ZOI of any particular regional minima Rj where j= 1, 2,...., N. This SKIZ surface forms the watershed lines separating the zone of influence of different regional minima. We have implemented the method of region growing in watershed based on the work by Vincent and Soille [15], and Beucher and Meyer [13], and extended the method with few modifications to encompass the 3-D objects. Here we suppose that the image voxels are stored in a simple array and direct access to any of the voxels of the image volume is possible. The regional markers can be obtained by either thresholding the distance map or by super imposing the seeds on the distance map. The seeds obtained as explained in Section 4.1, can also be considered as regional minima. In case of multiple seeds representing the same cell, marking the regional minima by seeds results in over-segmentation. Our experience says that, over-segmentation occurs even if we use the threshold distance map for marking the regional minima. This may be due to noisy holes that are still present in the foreground. After obtaining the regional minima for each cell, they are labeled using component labeling algorithm. If di is the distance value of the regional minima i, then consider all the voxels at distance (di + 1) in an ordered queue. Extend the labeled regional minima into the voxels in the queue having distance value (di + 1). Here we can also check that each voxel in the queue is having at least one voxel in its neighborhood as a voxel of the labeled regional minima. The process is carried out for all regional minima i= 1, 2,…, N where N is the total number of labeled regional markers. The process is repeated till all the regional minima are extended completely into their respective ZOI. The resulting segmented image volume is shown in Fig. 12b. This with all the preprocessing done,

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

41

Fig. 12. Result of 3-D watershed with rule-based merging, (a) original image slices, (b) result of simple 3-D watershed, (c) result after rule-based merging.

42

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

Fig. 12. (Continued)

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

still shows over-segmentation of some of the cells resulting in tiny noisy, objects. This over-segmented result as such may not be acceptable. All the tiny objects which are the parts of the larger cells have to be merged with corresponding cells. This has to be done automatically. The following section explain this merging process.

4.3.1. Rule based merging of the o6er-segmented cells The over-segmentation of the cells by the direct construction of the watershed lines in the image domain, is due to the fact that every regional minimum becomes the center of a catchment basin, i.e. every regional minimum is considered as belonging to a unique cell and its ZOI is calculated. It is to be noted that when we directly reconstruct the watershed lines, all the regional minima may not be of importance. Some of the regional minima may have been produced by noise and others by minor structures in the image such as intracellular matters, textured cell chromatin, dense intercellular fluids, etc. Fig. 12b shows the over-segmentation of a single image slice in the stack, due to direct watershed line construction. From Fig. 12 b, it is evident that the 3-D watershed technique is not foolproof and there are several cases where it over-segments the image volume, though several precautions such as extensive preprocessing is undertaken to avoid the errors in segmentation. Several researchers have tried to address this problem [13,26]. One possible method is make use of hysterisis thresholding to suppress the noisy, weak contours representing the watersheds between small regions. Najman and Schmitt [26], have listed various reasons for ill suitability of hysterisis thresholding in case of watershed. They are: 1. Hysterisis thresholding is suitable for edges and it produces non-closed contours while watershed contours are already closed contours which are obtained as a complimentary to the set of regions. 2. Hysterisis thresholding on watershed segmentation produces barbs, etc. Some of the proposed method to overcome the over-segmentation problem is the hierarchical segmentation [13], hierarchical segmentation

43

using dynamics of contour [25], etc. We have proposed an extension to the classical 3-D watershed which identifies the over-segmentation based on simple size and shape features and merge them with their parent cell. Instead of finding one marker for each object at the initial stages, which is a highly complex process, we allow the watershed to build on the enhanced image volume. Let N be the total number of segmented objects due to the application of a classical 3-D watershed segmentation method. It means there were N regional markers separated by a watershed contour. Let Tsize be the size threshold, i.e. a cell should have a minimum size of Tsize. All the tiny fragments of the cells whose sizes are below Tsize are considered as noisy and need to be merged to corresponding parent cell. We assign a fragment a to a parent cell A, if, the following conditions hold good. 1. ‘A’ and a should be touching neighbors. 2. If a is sharing its boundary with more than one large cell fragment, then the length of the boundary it shares with ‘A’ should be larger than the lengths of the boundary it shares with any other touching large object. If these conditions are satisfied then the fragment a is merged with the parent cell ‘A’. In case when a does not have any large touching neighbor but is a congregation of small fragments, there is a possibility that a single cell has been fragmented to such a level that there is no single fragment of size above Tsize. In such cases, all the tiny fragments which are touching each other and not connected to a larger object, are merged to form a single object. If this merged object is above the threshold Tsize then it is considered as a cell otherwise it is discarded as a noise artifact. Fig. 12c shows the result of automatic merging of the small objects in the over-segmented image volume. There can be still many cases where large individual cells are found to be segmented into two or three different objects. This may be due to the fact that such artifacts do not fall into the size threshold chosen by trial and error methods. This merging based on heuristic conditions gives better result than the result obtained by hysterisis thresholding. Table 1 gives the comparative result of segmentation by classical 3-D wa-

44

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

Table 1 Comparison of different methods to reduce the over-segmentation Sl. No.

Actual number of cells present

Number of cells by simple 3- Number of cells by selective D watershed marker technique

Number of cells by rule-based merging technique

01 02 03 04 05 06 07 08 09 10 11 12 13 14 15

23 31 18 42 26 12 05 31 28 37 11 12 16 28 07

42 56 29 67 38 19 05 49 43 56 11 12 21 44 09 % Error= 53%

24 32 21 38 26 12 05 33 31 40 13 12 17 31 08 % Error= 2%

24 31 23 47 29 12 05 34 31 42 12 12 18 30 07 % Error=9%

tershed with suppression of over-segmentation by selective marker techniques and by our rule based merging method.

5. Experimental results and discussion A Silicon Graphics, Indy machine with IRIX 5.3 was used as a platform to run the software. Much of the programming is done using interactive data language (IDL) and C. Most of the experiments needed no human interaction except to load the images. Figs. 8, 10 and 12 show the results of the three region based segmentation methods presented in this paper. For the purpose of qualitative evaluation, we have presented the results on the same data set. Our experience says that seeded region growing is useful mostly when the cells are not compactly arranged like in higher grades of carcinoma. The successive erosion – dilation technique can be used for cutting off the touching cells while watershed technique is useful in isolating cells in cluster and gives acceptable results when used with rule based merging technique to reduce the over-segmentation. For a 256×256×24 size image data, seeded volume growing took 1 min 38 s for segmentation with

80% correct segmentation results. The erosion–dilation technique took 2 min 15 s with 95% correct segmentation. The 3-D watershed technique with rule based merging took 3 min 30 s with 95% correct segmentation. Without rule based merging, results can be obtained within a minute but the segmentation is found to be only 53% correct. These figures are mentioned with respect to a particular data set. For a reliable figures, one has to experiment on a large amount of data before deciding the efficiency of a particular methodology. In the seeded region growing, the marked region gets extended to the neighborhood voxels based on the voxel similarity conditions as well as global morphological criteria. It is difficult to define the zone of influence of each seed as we do not use any distance measure for region growing. In a highly textured image volume, the seeded region growing fails to produce the proper segmentation due to dissimilarities in the voxel properties and uneven as well as low gradient magnitude characterizing the boundary of the cells. The normalization to reduce the effect of noisy peaks in intensity results in slowing down of the whole process. An extensive preprocessing is needed to smooth out the fine texture and influ-

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

ence of noisy voxel values on the volume growing. In the same time the edges which constitutes the cell boundary should be enhanced so that no over growth of the seed in to the region belonging to other cell occurs. The anisotropy in the voxel lattice, uneven illumination, photobleaching, etc. adversely affects the performance of the seeded volume growing technique. In a malignant, tumor tissue specimen images, the cells can not be put under strict convexity and shape constraints, making it very difficult define the global morphological criteria and/or to aggregate the voxels to fit the cell under some convexity assumption. More often seeded volume growing causes the under segmentation due to over growth of the seed of one cell into the touching and/or overlapping neighboring cell. In a highly textured cell, the growth may be too slow and stop before marking the complete cell volume resulting in error. Since many voxels fall outside the similarity measures, they may left unmarked even after the process of volume growing is completed. Most of these factors results in under segmentation. Another cause for error is the presence of holes in a segmented cells due to the dense intracellular matters which do not fall into the similarity criteria. These holes have to be identified and the gray values have to be restored using size and shape filters. The signal to noise ratio is directly proportional to the better performance of the segmentation using volume growing. Presence of noise and other artifacts

45

may result in under growth of the seed or in growing into the region belonging to the other cells. The successive pealing and thickening method appears to be very simple and effective. When the precision of the segmentation is not important as in finding the cell density, testing the clustering property of the cells, etc. this technique can be used for segmentation. It is simple and an effective means to cut off the touching and overlapping cells in a tissue image. One of the problems in this technique is, a cell may break down into more than one fragment while successive pealing is done. Each fragment is then treated as a separate object while thickening. This results in oversegmentation. In case the image contain cells which differ predominantly in their size and shape characters, the segmented image, using successive pealing and thickening technique, may not reflect the variation of shape and size of the cells as in original image. This is because, the separation of the touching cells need not be along the local gradient peaks. Moreover, the precision of the segmentation is limited by the low-level thresholding done on the original intensity image. The separation between the touching and overlapping cells need not be along the local gradient peaks. In the three dimensional watershed algorithm, the regional minima gets extended to the neighborhood voxels based on the geodesic distance of the neighborhood voxel and the zone of influence

Table 2 Comparative results of the number of cells uniquely labeled by different region based segmentation methods Sp. No.

Manual segmentation

Controlled (seeded) volume growing

Successive pealing and thickening

3-D Watershed with heuristic merging

1 2 3 4 5 6 7 8 9 10 11

23 31 18 42 26 12 5 31 28 37 11

21 26 18 33 25 12 5 27 21 33 11

23 35 18 48 31 12 5 30 37 41 11

24 32 21 38 26 12 5 33 31 40 13

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

46

Table 3 Comparative study of the different region based segmentation method corresponding to their ability in measuring simple features a Cell no.

Manual segmentation

Controlled (seeded) volume Successive pealing and growing thickening

3-D Watershed with heuristic merging

Symm. differ. (%)

Shape factor Symm. differ. (%)

Shape factor Symm. differ. (%)

Shape factor Symm. differ. (%)

Shape factor

Sp. No. 1 1 2 3 4 5

0 0 0 0 0

0.656 0.621 0.857 0.633 0.757

+5 −2 −6 +6 +4

0.638 0.577 0.839 0.629 0.739

+7 +3 +4 +7 +4

0.791 0.633 0.932 0.519 0.780

+4 −2 −5 +3 +5

0.645 0.593 0.902 0.647 0.785

Sp. No. 2 1 2 3 4 5 6

0 0 0 0 0 0

0.581 1.221 0.587 0.633 0.589 0.460

+7 +5 +3 −4 −6 +1

0.570 1.416 0.52 0.642 0.552 0.453

+7 +5 +5 +6 +3 +7

0.463 0.971 0.633 0.656 0.598 0.476

+5 +3 +3 −4 −5 −2

0.563 1.201 0.614 0.455 0.575 0.483

Sp. No. 3 1 2 3 4

0 0 0 0

0.988 0.831 0.569 0.683

+3 −8 −3 +4

1.134 0.756 0.578 0.605

+5 +3 +6 +5

1.247 0.873 0.591 0.671

+4 −3 −1 +2

1.122 0.783 0.591 0.675

a

Results are shown on a few selected cells of three representative image data.

of each regional marker. Each cell should have only one regional minima/marker to avoid the splitting of a single cell into two different object regions causing the over-segmentation of the image volume. Sometimes an interactive correction of the segmentation results is necessary. The 3D watershed technique is a relatively complex process and the precision of the result is very much dependent on proper selection of regional markers and the marking of zone of influence of each regional minima. Table 2, gives the comparative results based on the ability to mark segmented correct number of cells by each region based method. Table 3, shows the comparative results of seeded volume growing technique, successive pealing and thickening and 3-D watershed with merging with respect to the ability to give the precision in measuring the features such as percentage of symmetric difference in total volume and 3-D shape factor. The results are compared with

manually segmented cells. In the future work, we plan to study the influence of different region based segmentation techniques on the measurement of important cytological and histological features. We also plan to combine these techniques to produce better segmentation results.

Acknowledgements We would like to thank all the scientists, especially K. Rodenacker, Dr Hutzler and Dr Aubele of Biomedical Image Analysis Division, Institute of Pathology, GSF, Munich, Germany for providing the data sets. The work was partly supported under the project Indo–German collaboration of Biomedical research sponsored by International Beuro, Germany and ICMR, India. We thank both organizations for their support.

P.S. Umesh Adiga, B.B. Chaudhuri / Computer Methods and Programs in Biomedicine 61 (2000) 23–47

References [1] R.S. Acharya, C.J. Cogswell, D.B. Goldgof (Eds.), Biomedical Image Processing and 3-D Microscopy, SPIE Proceedings, vol. 16660, 1995. [2] P.S. Umesh Adiga, B.B. Chaudhuri, An efficient cell segmentation tool for confocal microscopy tissue images for quantitative evaluation of FISH signals, J. Microscop. Res. Tech. 44 (1999) 49–68. [3] K. Rodenacker, M. Aubele, P. Hutzler, P.S. Umesh Adiga, Groping for quantitative 3-D image analysis: an approach to quantitative evaluation of fluorescence in situ hybridization in thick tissue sections of prostate carcinoma, Anal. Cell. Pathol. 15 (1997) 19–29. [4] J.M. Chassery, C. Garbay, An iterative segmentation method based on a contextual color and shape criterion, IEEE Trans. on Pattern Anal. Mach. Intell. 6 (1984) 165 – 172. [5] J.K. Udupa, V.G. Ajjanagadde, Boundary and object labeling in three dimensional images, Comput. Vis. Graph. Image Processing: Image Understand. 3 (1990) 355 – 369. [6] L. Thurfjell, E. Bengtsson, B. Nordin, A new 3-D connected component labeling algorithm with simultaneous object feature extraction capability, Int. J. Graph. Models Image Understand. 54 (1992) 357–364. [7] R. Jain, R. Kasturi, B.G. Schunk, Machine Vision, McGraw-Hill, New York, 1995. [8] F. Cheevasuvit, H. Maitre, D. Vidal Madjar, A robust method for picture segmentation based on split-andmerge procedure, Comput. Vis. Graph. Image Processing 34 (1986) 268 – 281. [9] S.W. Zucker, Region growing: childhood and adolescence, Comput. Graph. Image Processing 5 (1976) 382– 399. [10] R. Kohler, A segmentation system based on thresholding, Comput. Graph. Image Processing 15 (1981) 319– 338. [11] R. Adams, L. Bischoff, Seeded region growing, IEEE Trans. Pattern Anal. Mach. Intell. 16 (1994) 641–647. [12] J.C. Russ, The Handbook of Image Processing, CRC Press, London, 1995. [13] S. Beucher, F. Meyer, The morphological approach to segmentation: the watershed transformation, in: E.R. Doughert (Ed.), Mathematical Morphology in Image Processing, Morcel Dekker, NewYork, 1993.

.

47

[14] C. Garbay, Image structure representation and processing: a discussion of some segmentation methods in cytology, IEEE Trans. Pattern Anal. Mach. Intell. 6 (1986) 165 – 172. [15] L. Vincent, P. Soille, Watershed in digital spaces: an efficient algorithm based on immersion simulation, IEEE Trans. Pattern Anal. Mach. Intell. 13 (1987) 583 – 598. [16] J.C. Alers, P.J. Krijtenburg, K.J. Vissers, S.K. Krishnadath, F.T. Bosman, H. Van Dekken, Interphase in situ hybridization to desegregate and intact tissue specimens of prostatic adeno carcinomas, Histochem. Cell Biol. 104 (1995) 479 – 486. [17] P.K. Sahoo, S. Soltani, A.K.C. Wong, Y.C. Chen, A survey of thresholding techniques, Comp. Vis. Graph. Image Processing 41 (1988) 233 – 260. [18] A. Rosenfeld, D. la Torre, Histogram concavity analysis as an aid in threshold selection, IEEE Trans. Syst. Man. Cybern. 13 (1983) 231 – 235. [19] H. Chen, J.R. Swedlow, M. Grote, J.W. Sedat, D.A. Agard, The collection, processing and display of digital 3-D images of biological specimens, in: J.B. Pawley (Ed.), Handbook of Biological Confocal Microscopy, Plenum Press, New York, 1995, pp. 197 – 209. [20] N. Aslund, A. Liljeberg, P.O. Forsgren, S. Wahlesten, 3-D digital microscopy using PHOIBOS scanner, Scanning 9 (1987) 227 – 235. [21] J. Rigaut, S.C Gonzalez, J. Vassey, 3-D Cytometry, in: A. Kriete (Ed.), Visualization in Biomedical Microscopies: 3-D Imaging and Computer Applications, VCH, NewYork, 1992, pp. 205 – 248. [22] R.M. Hralick, L. Shapiro, Computer and Robot Vision, Addison-Wesley Publishing Company, NewYork, 1992. [23] C. Lantuejoul, La Squelettisation et Son Application aux Mesures Topologiques des Mosaiques Polycrystalline, Ph.D. Thesis, Paris School of Mines., 1978. [24] S. Beucher, The watershed transformations applied to image segmentation, Scanning Microscop. 6 (1992) 299 – 314. [25] G. Borgefors, On digital distance transforms in three dimensions, Comput. Vis. Graph. Image Processing 64 (1996) 368 – 376. [26] L. Najman, M. Schmitt, Geodesic saliency of watershed contours and hierarchical segmentation, IEEE Trans. Pattern Anal. Mach. Intell. 18 (1996) 2.