Colour Mathematical Morphology For Neural Image Analysis

Real-Time Imaging 8, 455–465 (2002) doi:10.1006/rtim.2002.0288, available online at http://www.idealibrary.com on Colour Mathematical Morphology For Neural Image Analysis T his paper presents an algorithm for automatic neural image analysis in immunostained vertebrate retinas. We present a useful tool for cell quantification avoiding the losst of information of traditional binary techniques in automatic recognition of images. The application is based on the extension of the mathematical morphology to colour images. In qthe paper, we define the basics and more complex morphological operations to vectorial image processing. We propose and demonstrate a colour image reconstruction by geodesic transformations. In addition, we adapt the morphological segmentation of greyscale image to the segmentation of multispectral images of retinas of monkeys. # 2002 Elsevier Science Ltd. All rights reserved. F. Ortiz1, F. Torres1, E. De Juan2 and N. Cuenca3 1 Department of Physics, Systems Engineering and Signal Theory University of Alicante, P.O. Box 99, 03080 Alicante, Spain E-mail: fortiz@disc.ua.es 2 Department of Physiology, Genetics and Microbiology University of Alicante, P.O. Box 99, 03080 Alicante, Spain E-mail: fortiz@disc.ua.es 3 Department of Biotechnology, University of Alicante P.O. Box 99, 03080 Alicante, Spain Introduction For years, analysis of biomedical images has been increasingly important in medical research and diagnosis generation. The use of image processing tools is very essential when images to analyse are very numerous or these do not have good quality. The medical researcher only interprets the results and offers his final diagnosis with the automatically processed image. In this paper, we present a new algorithm for segmentation and classification of neural images of monkeys. Number and diameter quantification of cells is important in cell biology, but it is a hard and tedious work. The commercial and available software for automatic image analysis is very useful in this sense, 1077-2014/02/$35.00 but all the images need to be binarized in this software. The binary process transforms a greyscale or colour image to a binary image with a consequent loss of information: the programs confuse two cells that are connected like a single particle and cells with background. In the present study, we use images from normal retinas of monkeys to identify the number of dopaminergic and calretinin cells. Both cell types are important in the retinal visual information processing and are involved in retinal disease. In Parkinson’s disease, retinal dopaminergic and calretinin cells are altered and produce visual abnormalities. The quantification of these cell types allows us to determine the degree of damage in this disease. In this paper, we present an r 2002Elsevier Science Ltd. Alll rights reserved. 456 F. ORTIZ ETAL. algorithm for neural image analysis based on geodesic transformations for the reconstruction of colour images. We will extend the classical morphological operations to colour images. Morphological image processing is a nonlinear image processing developed by Matheron and Serra [1,2]. This processing technique has proved to be a powerful tool for many computer-vision tasks in binary and greyscale images, such as edge detection, noise suppression, image enhancement, skeletonization, segmentation, pattern recognition, etc. The extension of mathematical morphology to colour image is not straightforward [3]. Mathematical morphology is based on the set theory (complete lattice) where the notion of order is very important. In binary and greyscale images, the pixels are ordered by their value but, in colour images, each pixel is vector-valued (RGB, HSI, HSV, etc.). There is, therefore, no natural order for vectors, and as such, the extension of morphological operations to colour images requires a specific study of the order in multivariate data. In the following section, we comment the state of the art in colour morphology. Later, the method chosen for colour image processing will be shown. Next, we will explain the application of colour morphology in the processing of images from retinas of monkeys. Later on, we develop the algorithm used here and we show the results obtained. Colour Morphology and Lattice Theory In colour images, the pixels are represented by vectorial values in which vector element is a greyscale image P(x, y)=[P1(x, y), P2(x, y), P3 (x, y)]T [4]. Trahanias and Venetsanopoulos [5] summarized several techniques for ordering multivariate data. The two main approaches to processing are marginal ordering and vectorial ordering. In [3] Comer and Delp commented on the differences between marginal and vectorial processing. With marginal ordering, each component P1, P2 or P3 is ordered independently and the operations are applied in each colour channel of the image (Figure 1). The use of marginal ordering in colour image processing is the most straightforward approach. Nevertheless, this method may introduce visual changes in colour and may be unacceptable in applications that use colour for object recognition (as in our case). A vectorial P1 Processing P1 P2 Processing P2 P3 Processing P3 Figure 1. Marginal approach of processing. method for morphological processing is more advisable for avoiding the above-mentioned disadvantages. In vectorial ordering only one processing is done on threedimensional (3D) data (Figure 2). In vectorial data, there are several ways of establishing the order: K K K K Ordering by one component. Canonical ordering. Ordering by distance. Lexicographical order In ordering by one component, the ordering is decided by just one element. In canonical ordering, all three components of a colour space must have either higher or lower values than another vectorial colour. In ordering by distance, a distance function is used as an order measure, etc In [6,7] these methods are discussed in greater detail. The lattice description of morphology allows morphological theorems and techniques to be applied to images other than binary or greyscale [8–10]. A lattice is a partially ordered set in which any two elements possess a least upper bound (called supremum) and a greatest lower bound (infimum). The supremum and the infimum are represented by the symbols _ and ^, respectively. P1 P2 P3 P1 Processing P2 P3 Figure 2. Vectorial approach of processing. The colour of each pixel is a vector. ALGORITHM FOR NEURAL IMAGE ANALYSIS A lattice is complete if every subset of the lattice has a single supremum and a single infimum. The morphological operations must fulfil this latter condition. In order to calculate a dilation or an erosion, the notion of supremum and infimum is very important. Colour Spaces and Vectorial Processing Several coordinate systems are available for representing colour images [4]. The most common one is the RGB colour system (red, green and blue components). Nevertheless, for image processing it is more advisable to use HSI, HSV or HLS colour spaces. These spaces are more closely akin to the human interpretation of colours. The components of these colour models are the human perceptual attributes of colour: hue, saturation and luminance or intensity. The luminance or intensity is the best attribute for image processing and the hue map contains all the spectral colours. As such, we use the HSV colour space for processing. In this colour system, an image can be represented by f : Z2 ! R3 : X ¼ ðx; yÞ ! P ¼ ðH; S; VÞ ð1Þ where 0  S  1, 0  V  1 and 01  H  3601. In a discrete lattice, these values are scaled to integers in the range 0–255. Figure 3 illustrates the geometrical interpretation of the HSV model. The transformation from RGB to the HSV model has been carried out with the Foley and van Dam formulas [11]. We use the lexicographical method to form a complete lattice structure in HSV. This algorithm requires an internal order in each of the components H, S, V, and another order or preference between the components [12,13]. The lexicographical method is the order in which words are arranged in dictionaries: first, the order is decided by a component, followed by a second, and finally by a third value. There are a number of ways of ordering the H,S,V components, relative to one another. The preference or disposition of the components depends on the application and the properties of the images. Ordering with luminance in the first place is the best way of preserving the contours of the objects in the image. In situations in which the objects of interest are highly coloured or in which only objects of a specific colour are of interest, the operations with hue in the first position are the best. Due to the specific shape of the HSV space, a problem arises with the lexicographical method. Saturation and value are totally ordered sets, but hue is not. It is angle valued, H(x,y) 2 [0, 2p). Hue is also module coordinate: a hue angle y=y+2p. In addition, one cannot order hues from the lowest to the highest value. It does not make any perceptual sense to say, for example, that blue is greater than red. To order hues, Hanbury [13] and Peters [14] used a hue-valued structuring function. Hues are ordered according to the absolute value of a distance function between the image hue and a reference hue. The hue circles are partially ordered through the magnitude of the distance: dðHi ; Href Þ ¼ (


Hi Href


2p
Hi Href
θ Geometrical interpretation of HSV colour space. ð2Þ 0º Reference hue Infimun Supremun Figure 3.
Href
 p
Href
4p π/2 π S
if
Hi
if
Hi The infimum in the hue set is the reference hue and the supremum is ((infimum + p) mod 2p). Figure 4 shows the lattice of the hue circle with red (y=01) as the reference hue (infimum). H V 457 3π/4 Figure 4. Hue circle of HSV colour space. Reference hue (infimum) located in 01. 458 F. ORTIZ ETAL. Dilation of hue is defined as the selection of the image hue value with the greatest absolute distance. Likewise, erosion of hue is defined as the selection of the hue pixel value that generates the least measurement of hue difference. Peters does not address how a pixel is chosen if two distinct colour vectors have identical distance measurements. Hanbury uses an angle criterion. In this situation, we define infimum or supremum as the first to be chosen by the positioning of the pixels on the structuring element. Thus, we impose a total ordering on the hue component. Immunostaining Technique in Wholemount Retinas Figure 5. Now, we explain the process to obtain images of retinas of monkeys with immunocytochemistry method. Application for Retinas of Monkeys Retinas of monkeys were fixed in 4% paraformaldehyde in 0.1 M phosphate buffer at pH 7.4 for 1 h and then washed in 0.1 M PB and cryoprotected in 15% sucrose for 1/2 h, 20% sucrose for 1 h and 30% sucrose overnight at 41C. The following day, wholemount retinas were put through a freeze–thaw procedure. The retinas were incubated in 10% normal donkey serum in 0.1 M PB 0.5% Triton X-100 for 1 h at 41C. Without washing, the retinas were then incubated in a mixture of two primary antibodies, goat anti-calretinin (1:500) and rabbit anti-tyrosine hydroxylase (1:500) in 0.1 M PB, 0.5% Triton X-100 for 4 days at 41C under agitation. After further washes in PB, secondary incubation took place with donkey antigoat IgG Rhodamine and donkey anti-rabbit (FITC) (Jackson Immuno-reagent) for 1 h at 1:100 dilution, with agitation at room temperature. The retinas were washed, mounted in Vectashield mounting medium (Vector) and coverslipped for viewing by confocal microscopy (Zeiss 510 microscope). The Pinhole diameter was 77 mm giving an optical slice thickness of less than 0.9 mm. Images were viewed and captured at single focus planes. Control sections were obtained by omitting the primary antibody. Double immunostained wholemount labeled retinas showed two populations of amacrine cells. The dopaminergic amacrine cells labeled in green and the calretinin amacrine cells (type AII amacrine cells) labeled in red. Figure 5 shows the resulting image. The image is highly saturated and the objects to be analysed present high luminance. Original image. Monkey retina. The objective of our algorithm is the segmentation and identification of the two cell types present in the original image in Figure 5. For the analysis of each cell type, it is necessary to establish its pattern of colouring, size and form beforehand. Figures 6 (a) and (b) show, in detail, the two cell types present in the original image. There is a great amount of calretinin amacrine cells. These cells can be connected or can overlap each other. The colour of these cells is red (hue angle about 01). The dopaminergic amacrine cells (labeled in green) are present in a lesser quantity. Their size is greater than that of the red cells. In addition, these cells have ramifications and a hole inside. For the vectorial processing, we use a lexicographical order (denoted by olex) in which the first ordering decision is the hue distance, followed by value, and then saturation (3). For an image f and two vector pixels f(xi,yi)=Pi and f(xj,yj)=Pj: Pi oPj 8 dðHi ; Href ÞodðHj ; Href Þ > > > > or > > < dðHi ; Href Þ ¼ dðHj ; Href Þ if and Vi oVj > > or > > dðHi ; Href Þ ¼ dðHj ; Href Þ > > : and Vi ¼ Vj and Si oSj ð3Þ We can define the basic vectorial morphological operations (erosion ev and dilation dv by using the order olex. The erosion of an image f by a flat structuring element ALGORITHM FOR NEURAL IMAGE ANALYSIS 459 Figure 6. Detail of monkey retina: (a) The calretinin amacrine cells (red cells), and (b) the dopaminergic amacrine cells (green cells). K can be written as follows: evK ðf Þðx; yÞ ¼ minolex f ðx þ s; y þ tÞ ð4Þ ðs;tÞ2K We also define the vectorial dilation as dvK ðf Þðx; yÞ ¼ maxolex f ðx s; y tÞ ð5Þ ðs;tÞ2K The new definition of minolex and maxolex is important because they are now vectorial operations. The minolex and maxolex are equivalent to vectorial set operators of infimum ^v and supremum _v , respectively. The infimum and supremum are calculated by the order of olex defined in Eqn. (3): 8 < Pi if Pi  Pj ð6Þ Pi ^v Pj ¼ : P if Pi 4Pj j Pi _v Pj ¼ 8 < Pi : P j if Pi We develop a colour vectorial reconstruction by dilation, in which a marker image is infinitely dilated (until stability or idempotency) by a squared structuring element K. The selection of the marker image is very important for obtaining good results in the reconstruction. Reference hue Infimun Pj π/2 ð7Þ if Pi oPj Algorithm The algorithm we propose for segmentation of colour images is divided into four modules: 1. 2. 3. 4. and green cells. For this purpose, we use an extension of geodesic transformations to colour images. Geodesic dilations and erosions are seldom used in practice. However, when they are iterated until stability, they allow the definition of powerful morphological reconstruction algorithms. Morphological reconstruction techniques are at the basis of numerous transformations of high level for removing all the unwanted objects from an image. Reconstruction module. Filtering module. Cell-separation module. Gradient and statistical module. The aim of the first module (M1) is the classification of the original image into two classes of objects: red cells θ π 3π/4 0º Supremun Figure 7. Hue circle of HSV colour space. Reference hue (infimum) is placed in 901. 460 F. ORTIZ ETAL. In our algorithm, the marker image g is an eroded image f by a squared structuring element K of size 5 5: g ¼ evK ðf Þ ð8Þ We define the extension of classical geodesic dilation [3,14] to vectorial data as follows: dvf ðgÞ ¼ ^v ðdv ðgÞ; f Þ ð9Þ The new reconstruction by dilation for colour images is defined as the geodesic dilation of g with respect to f until stability: h ð1Þ ðn dðnÞ ðgÞ ¼ d f dvf v vf 1Þ ðgÞ ðnþ1Þ where n is such that dðnÞ ðgÞ. vf ðgÞ ¼ dvf i ð10Þ Figure 8. Detail of geodesic reconstruction of colour images: (a) a section of original image, (b) marker image for red objects, (c) marker image for green objects, (d) red object reconstruction and (e) green object reconstruction. ALGORITHM FOR NEURAL IMAGE ANALYSIS 461 The image marker is different for the reconstruction of each cell type. For the red cells, the erosion is calculated by a reference hue of 01 (Figure 4). For the green cells the reference hue is defined at about 901 (Figure 7). It is important to remember that the infimum (^v ) and supremum (_v ) operators select pixels according to the order olex. Figure 8 (a) shows a more detailed section of the original image. Figures 8 (b) and (c) illustrate the marker images for the red and the green objects, respectively. The images processed by our extension of geodesic transformation to colour images are presented in Figures 8 (d) and (e). Contrary to the morphological opening, the reconstruction by dilation preserves the shape of the components that are not removed by the erosion: all image features that cannot contain the structuring element are removed, while the others are unaltered. In addition, our algorithm allows the identification of objects of any colour through an easy change in the reference hue angle, contrary to RGB morphological processing, in which only three colours (red, green and blue) can be perfectly identified. The following step in our algorithm is the image filtering (M2). Each resulting type of image in the first module is now separately processed as a greyscale image. We now need to eliminate the objects that, due to their form and size, cannot be cells. A simple morphological opening filter removes all features that are not big enough to contain the structuring element. Nevertheless, the opening filter can modify the shape of the objects. This way, we use a geodesic filter of area opening which removes the features that have an area (number of pixels) of less than the selected threshold (30 pixels in our images). Finally, in this module, a closing filter removes internal holes from the cells. The results can be seen in Figure 9(a) and (b). The third module is only executed in the processing of the calretinin amacrine cells (red cells). The aim of this module is to separate cells that are connected or overlap with another cell [15]. The marker-controlled segmentation provides us with a powerful tool for solving this problem [2,16]. The algorithm calculates the watershed transformation of the complement of the image of euclidean distances. The distances are calculated in each cell from its edge or contour line. The regional maxima of the distance function is used as a marker to prevent over-segmentation. All of the steps of the methodology used are illustrated in Figure 10. We can see that the Figure 9. Image filtering of red objects: (a) threshold of red objects, and (b) result of geodesic area opening and closing filters. cells that were connected are now perfectly disconnected. Finally, in the fourth module (M4), the morphological gradient is used to determine an edge map of cells. As such, each cell is identified (Figure 11(a) and (b)). Different statistics can be calculated with the identifica- Figure 10. Separation of cells by the intersection of watersheds and connected cells. 462 F. ORTIZ ETAL. M1 Color reconstruction M2 Threshold Area opening Closing M3 Distance Regional max Complement Watershed Intersection M4 Gradient Stats Figure 11. Detail of cells identification: (a) red cells and (b) green cells. Figure 12. Modules of the algorithm. Four steps of processing: M1 (reconstruction), M2 (filtering), M3 (cell separation) and M4 (gradient and statistics). tion of the cells, for example, the amount of red and green cells, area, etc. A labeling (connected component) of the binary images of the gradient is used for particle counting. The block diagram of modules is illustrated in Figure 12. All modules have been implemented in the software developed by the authors and executed in the Windows system. At present, there is no software for automatic neural image analysis that carries out all of these operations. The first module is the one that requires the longest execution time, since the reconstruction is an operation that employs many morphological colour operations. Results In this section we show the results of our investigation: the output images of the different steps of our algorithm. Figure 13. Vectorial colour reconstruction by dilation of calretinin amacrine cells (red cells). ALGORITHM FOR NEURAL IMAGE ANALYSIS Figure 14. Filtering of the reconstructed image of red cells. Firstly, we show the results for the analysis of the calretinin amacrine cells (Figures 13–17). Secondly, Figure 18–20 illustrate the successive steps for obtaining the identification and the statistics of the dopaminergic amacrine cells. Figure 13 is the result of the application of the vectorial geodesic reconstruction of the original image. The green objects have been dimmed. Figure 14 shows the output of module 2. The possible cells have been selected by an area opening and a closing filter. The geometrical features in the cells have been preserved. Figure 15. Image of distances. 463 Figure 16. Intersection of watershed transformation and filtered image of figure 12. Separation of overlapping cells. The processing results of module 3 are illustrated in Figures 15 and 16. Figures 15 shows the distance function of the cells. The influence zones of each cell are determined by the watershed lines (Figure 16). The intersection between the watershed transformation and the filtered image of Figure 14 separates the connected cells. Finally, the algorithm calculates the gradient of the separate cells. The identification of the red cells on the original image is shown in Figure 17. The labeling used for counting cells gives a value of 179. Figure 17. Identification of calretinin amacrine cells (red cells) in original image. 464 F. ORTIZ ETAL. Figure 18. Vectorial colour reconstruction by dilation of dopaminergic amacrine cells (green cells). For the other type of objects in the original image (dopaminergic amacrine cells), module 1 of our algorithm also calculates a vectorial geodesic reconstruction. The output image of the colour reconstruction step is illustrated in Figure 18. The value of the green cells is greater than that of the red cells. The filtering of module 2 (Figure 19) provides a simple binary image in which only the possible cells are present. The edge detection of the cells by morphological gradient offers the identification of the green cells (Figure 20). The amount of green cells present is 2. Conclusions In this paper, we have presented an algorithm for analysing immunostained retinas. An extension of the Figure 20. Identification of dopaminergic amacrine cells in original image. classical morphological operations to segmentation of multispectral images has been introduced. A good method for vectorial mathematical morphology in HSV colour space has been detailed and the colour geodesic transformations have been presented. Our software incorporates colour morphology to avoid the losst of information in cell quantification present in the thresholding or binarization of images. Thus, we improve the automatic recognition of images in cell biology. In future studies we shall continue the research and extension of numerical geodesy to colour images. Acknowledgements We would like to thank the members of the Center of Mathematical Morphology (ENSMP) for their useful comments on the extension of morphological tools to colour images. This work has been partially supported by DGESIC PB98-0972. References Figure 19. Filtering of the reconstructed image of the green cells. 1. Serra, J. (1988) Mathematical Morphology, Theoretical Advances, Vol. 2. New York: Academic Press. 2. Soille, P. (1999) Morphological Image Analysis. Principles and Applications. Berlin: Springer-Verlag. 3. Comer, M. & Delp, E. (1999) Morphological Operations for colour image processing. Journal of Electronic Imaging, 8: 279–289. ALGORITHM FOR NEURAL IMAGE ANALYSIS 4. Sangwine, S. & Horne, R. (1998) The Colour Image Processing Handbook. London, UK: Chapman & Hall. 5. Trahanias, P. & Venetsanopoulos A. (1992) Colour edge detectors based on multivariate ordering. In: Proceedings of SPIE, Visual Communications and Image Processing’92, Petros Maragos (Ed) University of Toronto, Vol. 1818, 1396–1407. 6. Chanussot, J. (1998) Approches Vectorielles ou Marginales pour le Traitement d’Images Multicomposantes. PhD thesis, E.I.S.A., University of Savoie. 7. Chanussot, J. & Lambert, P. (1998) Total ordering based on space filling curves for multivalued morphology. In: Proceedings of the 4th International Symposium on Mathematical Morphology and its Applications, Amsterdam, The Netherlands, pp. 51–58. 8. Serra, J. (1993) Anamorphoses and function lattices (multivalued morphology). In: Dougherty, E. (Ed), Mathematical Morphology in Image Processing, pp. 483–521. 9. Ronse, C. (1990) Why mathematical morphology needs complete lattices. Signal Processing, 21: 129–154. 10. Sartor, L. & Weeks, A. (2001) Morphological operations on colour images. Journal of Electronic Imaging, 10: 548-559. 465 11. Foley, J. & van Dam, A. (1982) Fundamentals of Interactive Computer Graphics, The Systems Programming Series. Reading, MA: Addison-Wesley. 12. Ortiz, F., Torres, F., Puente, S., Candelas, F. & Gil, P. (2000) Use of the hue/saturation/intensity colour spaces to the morphological processing of colour images. In: Proceedings of the 1st International Conference in Graphics and Image Processing, Saint-Etienne, France, pp. 219–224. 13. Hanbury, A. (2001) Lexicographical order in the HLS colour space. Technical Report N-04/01/MM, Centre de Morphologie Mathématique. École des Mines de Paris. 14. Peters II, A. (1997) Mathematical morphology for anglevalued images. In: Proceedings of SPIE, Non-Linear Image Processing VIII, Vol. 3026, pp. 84–94. 15. Lantuéjoul, C. (1982) Geodesic segmentation. Multicomputers and Image Processing: Algorithms and Programs. New York: Academic Press. 16. Meyer, F. (1992) Colour image segmentation. In: Fourth Int. Conf. on Image Proc. and its Appl., Maastricht, The Netherlands, 303–304.