lEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL 37, NO 6. J U N E 1990
647
Communications Unsupervised Classification of Cell Images Using Pyramid Node Linking FARSHID ARMAN
AND
JOHN A. PEARCE
Abstract-In this communication we describe a segmentation technique which combines two properties in an iterative and hierarchial matter to correctly segment and classify the given cell images. The technique is applied to digital images taken from microscope slides of cultured rat liver cells, and the goal is to classify these cells into one of three possible classes. The first class cells (I) are morphologically normal and stain the darkest. The second class cells (11) are slightly damaged showing both nuclear and cytoplasmic swelling with resultant lessening of staining affinity. The third class cells (111) are markedly damaged as demonstrated by the presence of cytoplasmic vacuolization, or are completely disintegrated. First class cells are classified by taking advantage of their staining affinity; the original gray level image is segmented into four gray levels. The darkest is then classified as type I. Type 111 cells are classified by using high busyness as a characteristic; the standard deviation of the original image is segmented into four busyness levels. The highest level is classified as type 111 cell. Assuming only the three cell types are present in any given image, the remaining non-background unclassified pixels are determined to belong to type I1 cells.
Fig. I . A typical input image of cultured, damaged rat hepatocytes. below, and from the neighboring nodes o n the same level. This results in several trees with roots at the upper most levcl and leaves on bottom level producing smooth simply-connected regions with well-defined boundaries on the original image. The pyramid node linking approach t o segmentation has many advantages over other approaches. First, in most segmentation schemes one faces the problem of defining a window size over which the property is t o be calculated; however, due t o the multiresolution property of pyramid node linking, all portions of the image will b e processed at the appropriate resolution [ 5 ] .Second, due to the iterative nature of the procedure, most decisions are reversible and are not final until stabilization has occurred. Furthermore, pyramid-based algorithms are easily implementable in parallel [6]-[8]. In our classification scheme, first class cells are classified by taking advantage of their staining affinity; pyramid node linking is applied to the original gray level image segmenting it into four gray levels. The darkest is then classified as type I. Type I11 cells are classified using high busyness [9] a s a characteristic; pyramid node linking is applied to the standard deviation of the original image segmenting it into four busyness levels. The highest level is classified as type 111 cells. Assuming only the three types are present in any given image, the remaining unclassified pixels, not belonging to the background, are determined to belong to type I1 cells. This paper is organized a s follows: Section I1 describes the segmentation procedure, Section 111 illustrates the process of classifying the segmentation results, Section IV reports and analyzes experimental results, and the last section summarizes the process and presents our concluding remarks.
I. INTRODUCTION Given a digital image of a microscope slide (Fig. 1) the task is to segment and classify the cells in the image into three separate classes. The cells under investigation are cultured rat hepatocytes classified based on texture and gray level. The three classes are defined a s follows: type I cells are morphologically normal and stain the darkest; type I1 cells are slightly damaged, showing both nuclear and cytoplasmic swelling with resultant lessening of staining affinity; and type 111 cells are either markedly damaged, as demonstrated by the presence of cytoplasmic vacuolization, or completely disintegrated. In the current manual classification method [1]-[3], a square grid, with each square being approximately the size of a cell, is superimposed on the slide. Each point of the grid is then classified as one of the three types depending on the cell type it is covering. The frequency of occurrence of each of the three cell types is recorded and divided by total number of cells examined. Using this method, in [ 11 10 000 cells and in [2] 2 0 000 cells were examined. It is apparent that existing computer a n d image processing technology could b e utilized to speed up the process significantly, decrease the vulnerability to the human error, and eliminate subjectivity that exists in the process of cell examination by humans. In o u r research, an existing method of image segmentation, called pyramid node linking [4], [5] has been applied with a few modifications t o cell segmentation. Briefly, pyramids are hierarchical data structures built by successively reducing the spatial resolution of the image to be segmented from one level to next. T h e reduction is in the form of averaging some, typically square, neighborhood to get one node for the next higher level. The process is repeated until there are four nodes remaining on the top-most level. Then, in a bottom-to-top iterative process, the nodes from level t o level are linked using information from the level above, the level
11. SEGMENTATION In calculating different properties for segmentation, one faces the problem of defining the area over which the property is calculated. If the window size is too small, then the algorithm might segment small microstructures a s different segments, an undesirable outcome. O n the other hand, if the window size is too large, the algorithm might overlook the desired segments, artificially coalescing adjacent regions. Therefore, it is necessary t o accurately pick the window sizes so that they are large compared to the pattern elements contributing to the property computation, but small enough so that the segments d o not overlap, thus preserving the
Manuscript received March 30, 1989; revised November 16, 1989. This work was supported in part by the National Science Foundation Grant EET 845 1 123. The authors are with Department of Electrical and Computer Engineering, The University of Texas at Austin, Austin, TX 78712. IEEE Log Number 9034761,
0018-9294/90/0600-0647$01.OO O 1990 IEEE
1
64 8
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 37, NO. 6. JUNE 1990
boundaries between segments. By using pyramid node linking [4], [ 5 ] we have overcome this problem. Briefly, pyramids are hierarchical data structures with a rectangular array of nodes in a top-to-bottom sequence of levels [5]. Image resolution decreases a s one moves from the bottom level (finest) to the top most (coarsest) level. The values of the nodes in each level are built by averaging the values of the nodes, in some neighborhood, of the level below (each pixel in the image represents a node in the bottom most layer). T h e node calculated in this manner is referred to a s the father of the nodes below. The values of the nodes can be the gray level intensity, local standard deviation, or edge pixels, among others. This averaging process is repeated until there are four nodes remaining on the top-most level. Next, in a bottom-to-top iterative process, the links between the nodes from level t o level are assigned weights using information from the level above, from the level below, and from the neighbor nodes on the same level; the weights are then used to recalculate the value of each node a s before. This two-step process is iterated until the weights have stabilized (the stabilization always occurs [lo]), each node has chosen one unique father, and has been assigned the value of the chosen father. This process results in smooth simply-connected regions with well-defined boundaries on the bottom level, the original image. Pyramid node linking is described in more detail in the appendix. A . Background Segmentation
Pyramid node linking segments the image, including the background, into several regions. However, in most cell images the background, which is the medium the cells are embedded in, can be segmented using much simpler and faster methods. Segmentation of the background in cell images prior to pyramid node linking has several advantages. First, since the background pixels are not included in the process and there are fewer pixels t o consider, node linking will be much faster. Second, because we only treat the cell pixels, the stabilization of the weights could b e achieved earlier, resulting in further speedup in execution time of the algorithm. In other words, no iterations are spent in segmenting cell-background boundary pixels. Third, we are assured that cell pixels will be segmented into maximum number of possible regions (see Appendix). In addition, since slide backgrounds may vary due to processing variations resulted during slide preparations, pyramid node linking could segment the background into several subregions. Therefore, the background is segmented prior to Pyramid Node Linking using a standard thresholding technique [9]. Vacuoles, which have the same gray-level intensity a s the background are consequently merged with the cell pixel using information such a s region area and the neighbors. After finding the background areas, the nodes belonging to it are flagged as undefined during the pyramid node linking process. The pyramid is built and linked a s described earlier applying the equations to the defined nodes only (compare Figs. 2 and 3). 111. CLASSIFICATION METHOD
Referring to the original definitions of the three different cell types, it is apparent that gray level as well as a measure of busyness have to be used in order to achieve the correct classification. Thus, using these two measures, two pyramids, a gray level pyramid and a standard deviation pyramid, were constructed. Type I cells are defined as morphologically normal and not affected by the pharmacological experimentations and they have the highest staining affinity of the three types. In order to classify these cells, a gray level pyramid was initialized [Fig. 3(a)] and the input image was segmented into four gray levels [Fig. 3(b)] and the pixels with darkest of the four gray levels were classified a s type I cell pixels. Type 111 cells are defined a s severely damaged cells, evident by heavy vacuolization of the cytoplasm. They could also b e disintegrated in some cases. These cells are characterized a s having a higher variability o r busyness [9] than the other two cell types. In order to measure this property, the standard deviation of the input image was calculated using a 3 x 3 window, (any other measure of busyness could have been used here). Consequently, a pyramid
Fig. 2. (a) Initial gray level pyramid of eight levels with the background included. (b) The same pyramid after 12 iterations. Notice only two gray levels at the top most level are linked to cell areas of level zero, the other two have been used to segment the background.
Fig. 3 . (a) Initial gray level pyramid of eight levels with the background excluded, notice the expansion of cell nodes as the levels increase. (b) The same pyramid after 12 iterations, the top most level now holds four distinct gray levels which are used only to segment the cell areas in level zero.
was initialized using the resultant image as the input image to level zero, and node linking was performed to segment the standard deviation image into four levels of busyness. Assuming all the nonbackground pixels in the image belong to one of the three cell types, type I1 cells were classified a s pixels not belonging to either of class I or class 111; the final result is shown in Fig. 4. Another advantage of pyramid node linking is that due to the linking process, the averages are not across regions boundaries [4];i.e., since each father uses its sons to recalculate its value, the boundaries across segments are well preserved. This results in high resolution of boundaries in the final segmentation (Fig. 5). IV. EXPERIMENTAL RESULTS W e have classified over 150 cells using the method described. The algorithm used is completely automatic and needs no human intervention. The area threshold, used for merging of the small regions after thresholding is set for microscope images taken under
649
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 37. NO. 6. JUNE 1990
area. Second, clusters of mitochondria tend to "fool" the algorithm; when the clusters occur in the type I1 cells, the algorithm classifies them as type I . This is due to the high staining affinity of the mitochondria. W e believe when used in great volumes, as is needed in the actual pharmacological studies, the errors mentioned above will average to zero, providing a much greater accuracy than is possible by manual classification.
Fig 4 Final output of the example shown in figure I , dark represents type I cells, gray repreaents type I1 cells, and white represents type 111, the detected background is shown in black
IV. CONCLUSION Through the hierarchical and iterative process of pyramid node linking we have combined t w o properties, namely gray level and standard deviation, to segment and classify digital images of microscope slides of cultured rat liver cells. This algorithm depends on no inputs from the user. The output of this algorithm is an image in which all pixels have a label classifying them into one of the three cell classes or the background. This image is in turn used to calculate the percent area of each cell type with respect to total cell area. The advantages of this algorithm versus the manual method include significant speed u p of cell classification, substantial decrease in vulnerability to human error, and the elimination of subjectivity existent in the process of cell examination by humans. In short, the major contribution of this communication has been to demonstrate the use of two local properties, namely gray level and standard deviation, in a hierarchical fashion to reliably segment and classify cell images with minimal input from the user. APPENDIX
PYRAMIDNODE LINKING Starting with the original 2" X 2" image in level zero, each node at level I , 0 < l 5 n - I , is calculated a s an average of 2c X 2c nodes at level 1 - 1 ( t h e integer c 2 1 is called the span factor [ I l l ) . Mathematically: let a node at location (i,j ) in level 1 be assigned the location [ i, j, I ] , then if c = 2, there exists 16, (2c)', nodes located between
and
[2i - 1, 2j
-
1, 1
-
11,
[2i - 1, 2j
+ 2, 1
-
11,
[2i + 2, 2; - 1, 1 - 1 1 ,
[2i + 2 , 2j
+ 2, 1
-
11
(1)
on the level below called the sons of the node [ i, j , 1 1. The value of the node [ i, j , I ] , g [ i, j , I ] , is then initialized using the following equation:
g [ i , j, I ] Fig. 5. Boundaries of detected type I cells are shown in white overlapped on the original image.
lOOOx magnification and will require adjustment for a different magnification. Although the best judgment regarding the results of this algorithm could be made by a statistical comparison with the results obtained by manual classification, such comparison could be in error for two reasons. The goal of the classification is to obtain the percentage of each cell type area with compared to total cell area. However, hand classification (where a grid is superimposed and used a s the sampling point guide) results in cell (number) classification whereas our algorithm results in pixel (area) classification. In addition, hand classification, also known a s the relative volume percent (RVP) method [3], is accurate only when thousands of microscope slides are examined. Applying pyramid node linking t o as many cells as in [l] and [2] is not feasible due to the great volume of cell images that would have t o b e processed. The 150 cells classified have been compared visually with the actual definitions provided and the results are reasonably good. There are several sources of error in the classification process. First, cell-background edges as well as mitochondria-cytoplasm edges are both classified a s type I11 cells due to their high busyness characteristics. This, however, is a very small percentage of the image
=
1 C 7 4c k =all sons g [ i i , j;, I
-
11,
for 1 5 I i n
-
I.
(2) The node formed by this method is called thefather of the nodes below. Each node has c 2 fathers; in other words, it has been used in calculating the value of c z nodes at the level above. If c = 2 the four fathers of a node at [ i, j, I ] are located at
[ L(i [ L(i +
1 ) / 2 J , L(j
-
1)/2J, l
+ 11,
11/21> L(j
-
1)/2J I
+
[ L(i [ L(i
L(j + W J L(; +
11,
and -
+
1)/2J9
1)/2J, l + I ] W 2 J 1 + 11
(3)
where LxJ denotes the largest integer less than x . Next, in the node linking step, each node chooses its best father [4], [5].This choosing process is based on a closeness, in property value, between a node and its kth father, evaluated using ~5~ C5k = g[i, j, I ] - g[i;, j[, /
for k = I , 2,
+
+ 11
. . . , c 2 and 0
5
I
5
n
-
2
(4)
where g[i{, j{, I 1 1 is the kth father of g[i. j , I ] . Then, the node-father weight is
650
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 37, NO. 6. JUNE 1990
w(n.X,I
)
=
116: h = 311 fathers
f o r k = I , 2,
1,6~,
. . . , c’
[ I I ] W. 1. Grosky and R. Jain, “A pyramid-based approach to segmen-
(5)
where n is the node at level I , a n d h is its kth father. Then, using w , each node property g is recalculated:
tation applied to region matching,” IEEE Trans. Part. Anal. Mach. Intell., vol. PAMI-8, no. 5, pp. 639-650, Sept. 1986.
Real-Time Filtering for the Estimation of SteadyState Visual Evoked Brain Potentials THOMAS F. COLLURA
w h e r e g [ i l , j j , 1 - 11 i s t h e k t h s o n o f g [ i , j , l ] a n d w h e r e w ( n L , f , 1 - 1 ) is the weight between the node g [ i, j , I 1, f,and its k t h son nk at level 1 - 1. As the process iterates in a bottom-to-top fashion the weight function w approaches either zero or one, with one being representative of a node most likely to be the father and zero being representative of a node least likely t o be the father. In cases where 6, = 0 in (4), the kth node-father link gets the weight of one in ( 5 ) , and all the other c’ node-father links will get the weight of zero. It is possible that a node may not be picked a s a father by any of the nodes on the lower level (when the denominator of (6) is zero). In such an event, we let the node g [ i, j , I ] be the average of its eight brothers, and the iteration continues. During the final phase, the tree generation phase, the segment values s are assigned top-down starting at a chosen level L 5 n 1 to level 1 = 0. At level L
where [ i”, J ” , I 1 is the father of [ i, j , I ] chosen during node Iinking. At the end of this stage, s [ i , j , 01 will have the final segmentation of the image into 22“1-L’ homogeneous regions where the original image is 2“ X 2“; for L 5 3 there will be a maximum of 256 homogeneous regions. ACKNOWLEDGMENT The authors would like to thank Dr. E. M. B. Sorensen for providing the microscope slides used in this research. REFERENCES [ I ] E. M. B. Sorensen and D. Acostd, “Erythromycin estolate-induced
toxicity in cultured rat hepatocytes,” Toxicol. Lett., no. 27, pp. 7382, 1985. [2] -, “Protective effects of calcium,” Alternative Methods in Toxicology, Vol. 3: In Vitro Toxicology. New York: Liebert, 1985, pp. IO I - 139. 131 M. A. Hayat, Principles and Techniques of Electron Microscopy: Biological Application. New York: Van Nostrand Reinhold, 1970. ch. 6 , pp. 239-295. 141 T . H. Hong, K. A. Narayanan. S . Peleg, and A . Rosenfeld, “Image smoothing and segmentation by multiresolution pixel linking: Further experiments and extensions,” IEEE Trans. Sysr., M a n , Cybern., vol. SMC-12, no. 5, pp. 611-622, Sept./Oct. 1982. [SI P. J . Burt, T . H. Hong and A . Rosenfeld, “Segmentation and estimation of image region properties through cooperative hierarchial computation,” IEEE Trans. Syst., Man, Cybern., vol. SMC-I I , no. 12, pp. 802-809, Dec. 1981. 161 N. Ahuja and S . Swamy, “Multiprocessor pyramid architectures for bottom-up image analysis,” Multiresolution Image Procrssing and Analysis, A. Rosenfeld, Ed. New York: Springer-Verlag. 1984, ch. 3, pp. 38-58. [7] C. Lengauer, B. Sabata, and F. Arman, “A mechanically derived systolic implementation of pyramid initialization,” in Proc. Workshop on Hardware Specification Verijcution, Synthesis: Mathematical Aspects? Ithaca, NY. July 1989. [8] C. R. Dyer, “A VLSI pyramid machine for hierarchial parallel image processing,” in Proc. PRIP ’81, Aug. 1981, pp. 381-386. [9] A. Rosenfeld and A. Kak, Digital Picture Processing. New York:
Academic, 1982. 1101 S . Kasif and A . Rosenfeld, “Pyramid Linking is a special case of ISODATA,” IEEE Trans. Syst., Man, Cybern., vol. SMC-13. vol. I , pp. 84-85, Jan./Feb. 1983.
Abstract-It is shown that EEG visual evoked potentials elicited by repetitive stimuli in the range of 2 to 20 per second can be readily estimated in real time using a simple filtering approach. This measurement takes advantage of the fact that a comb filter will pass the important Fourier harmonics of the signal to provide an estimate of the evoked activity, plus track time-variations in the signal. Results on human subjects demonstrate the effectiveness of the approach.
INTRODUCTION Repetitive visual stimuli such a s flashing lights and reversing checkerboards are commonly used in the production of brainevoked responses for research and clinical studies. The brain response is buried in significant amounts of background EEG activity, which necessitates the use of some type of averaging or filtering to increase the signal-to-noise ratio. Various methods have been studied, and there is increasing emphasis on adaptive techniques and techniques which can estimate single responses [1]-[4]. One reason for this emphasis is the realization that successive evoked potentials are not identical, and that short-term changes in brain responsivity create significant variation between them [ 5 ] . A major problem in improving signal-to-noise ratio is the fact that there may be no “prototypical” evoked potential; the “average” evoked potential is a fiction, much a s the “typical man” who is characterized by a certain average height and weight, but may never be represented by a particular individual. For this reason, there is value in techniques which d o not make assumptions about the evoked potential other than the fact that it is time-locked to the stimulus. While averaging follows this assumption, amplitude and latency variations can reduce its accuracy in developing a valid estimate based upon many responses. Also, averaging of dozens of responses may be too slow to produce a good estimate before shortterm changes degrade the measurement. Digital techniques which overcome these limitations may not be suited t o low cost or realtime operation. This technique is founded o n the recognition that variation between successive evoked potentials not only degrades the accuracy of averaging, but is itself a valuable part of the information in the signal. This work provides a method which allows the investigator to track short-term changes in evoked responses, and to make them available as part of the measurement paradigm. METHOD The method is based on the strategy of increasing signal-to-noise ratio before the system has time to change significantly. In practical terms, it is necessary to produce an estimate of the evoked responses within 5 to I O s and to keep that estimate up to date, thus tracking changes in brain state. Using stimulus rates between 2 and 2 0 per second, successive cortical evoked responses lead to a periodic wave. Simple overlap (linear superposition) is a dominant mechanism, particularly at the lower rates [6]. This is shown diagrammatically in Fig. 1. Above four responses per second, these are commonly referred to a s “steady-state’’ evoked potentials. This periodic evoked wave can Manuscript received April 20, 1989; revised September 5 , 1989. This work was supported by NlH Training Grant GM-01090-14. The author is with the Section of Epilepsy and Clinical Neurophyisology, The Cleveland Clinic Foundation, Cleveland, OH 44195. IEEE Log Number 9035447.
00 I8-9294/90/0600-0650$0 1 .OO 0 1990 IEEE
647
Communications Unsupervised Classification of Cell Images Using Pyramid Node Linking FARSHID ARMAN
AND
JOHN A. PEARCE
Abstract-In this communication we describe a segmentation technique which combines two properties in an iterative and hierarchial matter to correctly segment and classify the given cell images. The technique is applied to digital images taken from microscope slides of cultured rat liver cells, and the goal is to classify these cells into one of three possible classes. The first class cells (I) are morphologically normal and stain the darkest. The second class cells (11) are slightly damaged showing both nuclear and cytoplasmic swelling with resultant lessening of staining affinity. The third class cells (111) are markedly damaged as demonstrated by the presence of cytoplasmic vacuolization, or are completely disintegrated. First class cells are classified by taking advantage of their staining affinity; the original gray level image is segmented into four gray levels. The darkest is then classified as type I. Type 111 cells are classified by using high busyness as a characteristic; the standard deviation of the original image is segmented into four busyness levels. The highest level is classified as type 111 cell. Assuming only the three cell types are present in any given image, the remaining non-background unclassified pixels are determined to belong to type I1 cells.
Fig. I . A typical input image of cultured, damaged rat hepatocytes. below, and from the neighboring nodes o n the same level. This results in several trees with roots at the upper most levcl and leaves on bottom level producing smooth simply-connected regions with well-defined boundaries on the original image. The pyramid node linking approach t o segmentation has many advantages over other approaches. First, in most segmentation schemes one faces the problem of defining a window size over which the property is t o be calculated; however, due t o the multiresolution property of pyramid node linking, all portions of the image will b e processed at the appropriate resolution [ 5 ] .Second, due to the iterative nature of the procedure, most decisions are reversible and are not final until stabilization has occurred. Furthermore, pyramid-based algorithms are easily implementable in parallel [6]-[8]. In our classification scheme, first class cells are classified by taking advantage of their staining affinity; pyramid node linking is applied to the original gray level image segmenting it into four gray levels. The darkest is then classified as type I. Type I11 cells are classified using high busyness [9] a s a characteristic; pyramid node linking is applied to the standard deviation of the original image segmenting it into four busyness levels. The highest level is classified as type 111 cells. Assuming only the three types are present in any given image, the remaining unclassified pixels, not belonging to the background, are determined to belong to type I1 cells. This paper is organized a s follows: Section I1 describes the segmentation procedure, Section 111 illustrates the process of classifying the segmentation results, Section IV reports and analyzes experimental results, and the last section summarizes the process and presents our concluding remarks.
I. INTRODUCTION Given a digital image of a microscope slide (Fig. 1) the task is to segment and classify the cells in the image into three separate classes. The cells under investigation are cultured rat hepatocytes classified based on texture and gray level. The three classes are defined a s follows: type I cells are morphologically normal and stain the darkest; type I1 cells are slightly damaged, showing both nuclear and cytoplasmic swelling with resultant lessening of staining affinity; and type 111 cells are either markedly damaged, as demonstrated by the presence of cytoplasmic vacuolization, or completely disintegrated. In the current manual classification method [1]-[3], a square grid, with each square being approximately the size of a cell, is superimposed on the slide. Each point of the grid is then classified as one of the three types depending on the cell type it is covering. The frequency of occurrence of each of the three cell types is recorded and divided by total number of cells examined. Using this method, in [ 11 10 000 cells and in [2] 2 0 000 cells were examined. It is apparent that existing computer a n d image processing technology could b e utilized to speed up the process significantly, decrease the vulnerability to the human error, and eliminate subjectivity that exists in the process of cell examination by humans. In o u r research, an existing method of image segmentation, called pyramid node linking [4], [5] has been applied with a few modifications t o cell segmentation. Briefly, pyramids are hierarchical data structures built by successively reducing the spatial resolution of the image to be segmented from one level to next. T h e reduction is in the form of averaging some, typically square, neighborhood to get one node for the next higher level. The process is repeated until there are four nodes remaining on the top-most level. Then, in a bottom-to-top iterative process, the nodes from level t o level are linked using information from the level above, the level
11. SEGMENTATION In calculating different properties for segmentation, one faces the problem of defining the area over which the property is calculated. If the window size is too small, then the algorithm might segment small microstructures a s different segments, an undesirable outcome. O n the other hand, if the window size is too large, the algorithm might overlook the desired segments, artificially coalescing adjacent regions. Therefore, it is necessary t o accurately pick the window sizes so that they are large compared to the pattern elements contributing to the property computation, but small enough so that the segments d o not overlap, thus preserving the
Manuscript received March 30, 1989; revised November 16, 1989. This work was supported in part by the National Science Foundation Grant EET 845 1 123. The authors are with Department of Electrical and Computer Engineering, The University of Texas at Austin, Austin, TX 78712. IEEE Log Number 9034761,
0018-9294/90/0600-0647$01.OO O 1990 IEEE
1
64 8
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 37, NO. 6. JUNE 1990
boundaries between segments. By using pyramid node linking [4], [ 5 ] we have overcome this problem. Briefly, pyramids are hierarchical data structures with a rectangular array of nodes in a top-to-bottom sequence of levels [5]. Image resolution decreases a s one moves from the bottom level (finest) to the top most (coarsest) level. The values of the nodes in each level are built by averaging the values of the nodes, in some neighborhood, of the level below (each pixel in the image represents a node in the bottom most layer). T h e node calculated in this manner is referred to a s the father of the nodes below. The values of the nodes can be the gray level intensity, local standard deviation, or edge pixels, among others. This averaging process is repeated until there are four nodes remaining on the top-most level. Next, in a bottom-to-top iterative process, the links between the nodes from level t o level are assigned weights using information from the level above, from the level below, and from the neighbor nodes on the same level; the weights are then used to recalculate the value of each node a s before. This two-step process is iterated until the weights have stabilized (the stabilization always occurs [lo]), each node has chosen one unique father, and has been assigned the value of the chosen father. This process results in smooth simply-connected regions with well-defined boundaries on the bottom level, the original image. Pyramid node linking is described in more detail in the appendix. A . Background Segmentation
Pyramid node linking segments the image, including the background, into several regions. However, in most cell images the background, which is the medium the cells are embedded in, can be segmented using much simpler and faster methods. Segmentation of the background in cell images prior to pyramid node linking has several advantages. First, since the background pixels are not included in the process and there are fewer pixels t o consider, node linking will be much faster. Second, because we only treat the cell pixels, the stabilization of the weights could b e achieved earlier, resulting in further speedup in execution time of the algorithm. In other words, no iterations are spent in segmenting cell-background boundary pixels. Third, we are assured that cell pixels will be segmented into maximum number of possible regions (see Appendix). In addition, since slide backgrounds may vary due to processing variations resulted during slide preparations, pyramid node linking could segment the background into several subregions. Therefore, the background is segmented prior to Pyramid Node Linking using a standard thresholding technique [9]. Vacuoles, which have the same gray-level intensity a s the background are consequently merged with the cell pixel using information such a s region area and the neighbors. After finding the background areas, the nodes belonging to it are flagged as undefined during the pyramid node linking process. The pyramid is built and linked a s described earlier applying the equations to the defined nodes only (compare Figs. 2 and 3). 111. CLASSIFICATION METHOD
Referring to the original definitions of the three different cell types, it is apparent that gray level as well as a measure of busyness have to be used in order to achieve the correct classification. Thus, using these two measures, two pyramids, a gray level pyramid and a standard deviation pyramid, were constructed. Type I cells are defined as morphologically normal and not affected by the pharmacological experimentations and they have the highest staining affinity of the three types. In order to classify these cells, a gray level pyramid was initialized [Fig. 3(a)] and the input image was segmented into four gray levels [Fig. 3(b)] and the pixels with darkest of the four gray levels were classified a s type I cell pixels. Type 111 cells are defined a s severely damaged cells, evident by heavy vacuolization of the cytoplasm. They could also b e disintegrated in some cases. These cells are characterized a s having a higher variability o r busyness [9] than the other two cell types. In order to measure this property, the standard deviation of the input image was calculated using a 3 x 3 window, (any other measure of busyness could have been used here). Consequently, a pyramid
Fig. 2. (a) Initial gray level pyramid of eight levels with the background included. (b) The same pyramid after 12 iterations. Notice only two gray levels at the top most level are linked to cell areas of level zero, the other two have been used to segment the background.
Fig. 3 . (a) Initial gray level pyramid of eight levels with the background excluded, notice the expansion of cell nodes as the levels increase. (b) The same pyramid after 12 iterations, the top most level now holds four distinct gray levels which are used only to segment the cell areas in level zero.
was initialized using the resultant image as the input image to level zero, and node linking was performed to segment the standard deviation image into four levels of busyness. Assuming all the nonbackground pixels in the image belong to one of the three cell types, type I1 cells were classified a s pixels not belonging to either of class I or class 111; the final result is shown in Fig. 4. Another advantage of pyramid node linking is that due to the linking process, the averages are not across regions boundaries [4];i.e., since each father uses its sons to recalculate its value, the boundaries across segments are well preserved. This results in high resolution of boundaries in the final segmentation (Fig. 5). IV. EXPERIMENTAL RESULTS W e have classified over 150 cells using the method described. The algorithm used is completely automatic and needs no human intervention. The area threshold, used for merging of the small regions after thresholding is set for microscope images taken under
649
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 37. NO. 6. JUNE 1990
area. Second, clusters of mitochondria tend to "fool" the algorithm; when the clusters occur in the type I1 cells, the algorithm classifies them as type I . This is due to the high staining affinity of the mitochondria. W e believe when used in great volumes, as is needed in the actual pharmacological studies, the errors mentioned above will average to zero, providing a much greater accuracy than is possible by manual classification.
Fig 4 Final output of the example shown in figure I , dark represents type I cells, gray repreaents type I1 cells, and white represents type 111, the detected background is shown in black
IV. CONCLUSION Through the hierarchical and iterative process of pyramid node linking we have combined t w o properties, namely gray level and standard deviation, to segment and classify digital images of microscope slides of cultured rat liver cells. This algorithm depends on no inputs from the user. The output of this algorithm is an image in which all pixels have a label classifying them into one of the three cell classes or the background. This image is in turn used to calculate the percent area of each cell type with respect to total cell area. The advantages of this algorithm versus the manual method include significant speed u p of cell classification, substantial decrease in vulnerability to human error, and the elimination of subjectivity existent in the process of cell examination by humans. In short, the major contribution of this communication has been to demonstrate the use of two local properties, namely gray level and standard deviation, in a hierarchical fashion to reliably segment and classify cell images with minimal input from the user. APPENDIX
PYRAMIDNODE LINKING Starting with the original 2" X 2" image in level zero, each node at level I , 0 < l 5 n - I , is calculated a s an average of 2c X 2c nodes at level 1 - 1 ( t h e integer c 2 1 is called the span factor [ I l l ) . Mathematically: let a node at location (i,j ) in level 1 be assigned the location [ i, j, I ] , then if c = 2, there exists 16, (2c)', nodes located between
and
[2i - 1, 2j
-
1, 1
-
11,
[2i - 1, 2j
+ 2, 1
-
11,
[2i + 2, 2; - 1, 1 - 1 1 ,
[2i + 2 , 2j
+ 2, 1
-
11
(1)
on the level below called the sons of the node [ i, j , 1 1. The value of the node [ i, j , I ] , g [ i, j , I ] , is then initialized using the following equation:
g [ i , j, I ] Fig. 5. Boundaries of detected type I cells are shown in white overlapped on the original image.
lOOOx magnification and will require adjustment for a different magnification. Although the best judgment regarding the results of this algorithm could be made by a statistical comparison with the results obtained by manual classification, such comparison could be in error for two reasons. The goal of the classification is to obtain the percentage of each cell type area with compared to total cell area. However, hand classification (where a grid is superimposed and used a s the sampling point guide) results in cell (number) classification whereas our algorithm results in pixel (area) classification. In addition, hand classification, also known a s the relative volume percent (RVP) method [3], is accurate only when thousands of microscope slides are examined. Applying pyramid node linking t o as many cells as in [l] and [2] is not feasible due to the great volume of cell images that would have t o b e processed. The 150 cells classified have been compared visually with the actual definitions provided and the results are reasonably good. There are several sources of error in the classification process. First, cell-background edges as well as mitochondria-cytoplasm edges are both classified a s type I11 cells due to their high busyness characteristics. This, however, is a very small percentage of the image
=
1 C 7 4c k =all sons g [ i i , j;, I
-
11,
for 1 5 I i n
-
I.
(2) The node formed by this method is called thefather of the nodes below. Each node has c 2 fathers; in other words, it has been used in calculating the value of c z nodes at the level above. If c = 2 the four fathers of a node at [ i, j, I ] are located at
[ L(i [ L(i +
1 ) / 2 J , L(j
-
1)/2J, l
+ 11,
11/21> L(j
-
1)/2J I
+
[ L(i [ L(i
L(j + W J L(; +
11,
and -
+
1)/2J9
1)/2J, l + I ] W 2 J 1 + 11
(3)
where LxJ denotes the largest integer less than x . Next, in the node linking step, each node chooses its best father [4], [5].This choosing process is based on a closeness, in property value, between a node and its kth father, evaluated using ~5~ C5k = g[i, j, I ] - g[i;, j[, /
for k = I , 2,
+
+ 11
. . . , c 2 and 0
5
I
5
n
-
2
(4)
where g[i{, j{, I 1 1 is the kth father of g[i. j , I ] . Then, the node-father weight is
650
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 37, NO. 6. JUNE 1990
w(n.X,I
)
=
116: h = 311 fathers
f o r k = I , 2,
1,6~,
. . . , c’
[ I I ] W. 1. Grosky and R. Jain, “A pyramid-based approach to segmen-
(5)
where n is the node at level I , a n d h is its kth father. Then, using w , each node property g is recalculated:
tation applied to region matching,” IEEE Trans. Part. Anal. Mach. Intell., vol. PAMI-8, no. 5, pp. 639-650, Sept. 1986.
Real-Time Filtering for the Estimation of SteadyState Visual Evoked Brain Potentials THOMAS F. COLLURA
w h e r e g [ i l , j j , 1 - 11 i s t h e k t h s o n o f g [ i , j , l ] a n d w h e r e w ( n L , f , 1 - 1 ) is the weight between the node g [ i, j , I 1, f,and its k t h son nk at level 1 - 1. As the process iterates in a bottom-to-top fashion the weight function w approaches either zero or one, with one being representative of a node most likely to be the father and zero being representative of a node least likely t o be the father. In cases where 6, = 0 in (4), the kth node-father link gets the weight of one in ( 5 ) , and all the other c’ node-father links will get the weight of zero. It is possible that a node may not be picked a s a father by any of the nodes on the lower level (when the denominator of (6) is zero). In such an event, we let the node g [ i, j , I ] be the average of its eight brothers, and the iteration continues. During the final phase, the tree generation phase, the segment values s are assigned top-down starting at a chosen level L 5 n 1 to level 1 = 0. At level L
where [ i”, J ” , I 1 is the father of [ i, j , I ] chosen during node Iinking. At the end of this stage, s [ i , j , 01 will have the final segmentation of the image into 22“1-L’ homogeneous regions where the original image is 2“ X 2“; for L 5 3 there will be a maximum of 256 homogeneous regions. ACKNOWLEDGMENT The authors would like to thank Dr. E. M. B. Sorensen for providing the microscope slides used in this research. REFERENCES [ I ] E. M. B. Sorensen and D. Acostd, “Erythromycin estolate-induced
toxicity in cultured rat hepatocytes,” Toxicol. Lett., no. 27, pp. 7382, 1985. [2] -, “Protective effects of calcium,” Alternative Methods in Toxicology, Vol. 3: In Vitro Toxicology. New York: Liebert, 1985, pp. IO I - 139. 131 M. A. Hayat, Principles and Techniques of Electron Microscopy: Biological Application. New York: Van Nostrand Reinhold, 1970. ch. 6 , pp. 239-295. 141 T . H. Hong, K. A. Narayanan. S . Peleg, and A . Rosenfeld, “Image smoothing and segmentation by multiresolution pixel linking: Further experiments and extensions,” IEEE Trans. Sysr., M a n , Cybern., vol. SMC-12, no. 5, pp. 611-622, Sept./Oct. 1982. [SI P. J . Burt, T . H. Hong and A . Rosenfeld, “Segmentation and estimation of image region properties through cooperative hierarchial computation,” IEEE Trans. Syst., Man, Cybern., vol. SMC-I I , no. 12, pp. 802-809, Dec. 1981. 161 N. Ahuja and S . Swamy, “Multiprocessor pyramid architectures for bottom-up image analysis,” Multiresolution Image Procrssing and Analysis, A. Rosenfeld, Ed. New York: Springer-Verlag. 1984, ch. 3, pp. 38-58. [7] C. Lengauer, B. Sabata, and F. Arman, “A mechanically derived systolic implementation of pyramid initialization,” in Proc. Workshop on Hardware Specification Verijcution, Synthesis: Mathematical Aspects? Ithaca, NY. July 1989. [8] C. R. Dyer, “A VLSI pyramid machine for hierarchial parallel image processing,” in Proc. PRIP ’81, Aug. 1981, pp. 381-386. [9] A. Rosenfeld and A. Kak, Digital Picture Processing. New York:
Academic, 1982. 1101 S . Kasif and A . Rosenfeld, “Pyramid Linking is a special case of ISODATA,” IEEE Trans. Syst., Man, Cybern., vol. SMC-13. vol. I , pp. 84-85, Jan./Feb. 1983.
Abstract-It is shown that EEG visual evoked potentials elicited by repetitive stimuli in the range of 2 to 20 per second can be readily estimated in real time using a simple filtering approach. This measurement takes advantage of the fact that a comb filter will pass the important Fourier harmonics of the signal to provide an estimate of the evoked activity, plus track time-variations in the signal. Results on human subjects demonstrate the effectiveness of the approach.
INTRODUCTION Repetitive visual stimuli such a s flashing lights and reversing checkerboards are commonly used in the production of brainevoked responses for research and clinical studies. The brain response is buried in significant amounts of background EEG activity, which necessitates the use of some type of averaging or filtering to increase the signal-to-noise ratio. Various methods have been studied, and there is increasing emphasis on adaptive techniques and techniques which can estimate single responses [1]-[4]. One reason for this emphasis is the realization that successive evoked potentials are not identical, and that short-term changes in brain responsivity create significant variation between them [ 5 ] . A major problem in improving signal-to-noise ratio is the fact that there may be no “prototypical” evoked potential; the “average” evoked potential is a fiction, much a s the “typical man” who is characterized by a certain average height and weight, but may never be represented by a particular individual. For this reason, there is value in techniques which d o not make assumptions about the evoked potential other than the fact that it is time-locked to the stimulus. While averaging follows this assumption, amplitude and latency variations can reduce its accuracy in developing a valid estimate based upon many responses. Also, averaging of dozens of responses may be too slow to produce a good estimate before shortterm changes degrade the measurement. Digital techniques which overcome these limitations may not be suited t o low cost or realtime operation. This technique is founded o n the recognition that variation between successive evoked potentials not only degrades the accuracy of averaging, but is itself a valuable part of the information in the signal. This work provides a method which allows the investigator to track short-term changes in evoked responses, and to make them available as part of the measurement paradigm. METHOD The method is based on the strategy of increasing signal-to-noise ratio before the system has time to change significantly. In practical terms, it is necessary to produce an estimate of the evoked responses within 5 to I O s and to keep that estimate up to date, thus tracking changes in brain state. Using stimulus rates between 2 and 2 0 per second, successive cortical evoked responses lead to a periodic wave. Simple overlap (linear superposition) is a dominant mechanism, particularly at the lower rates [6]. This is shown diagrammatically in Fig. 1. Above four responses per second, these are commonly referred to a s “steady-state’’ evoked potentials. This periodic evoked wave can Manuscript received April 20, 1989; revised September 5 , 1989. This work was supported by NlH Training Grant GM-01090-14. The author is with the Section of Epilepsy and Clinical Neurophyisology, The Cleveland Clinic Foundation, Cleveland, OH 44195. IEEE Log Number 9035447.
00 I8-9294/90/0600-0650$0 1 .OO 0 1990 IEEE
