Edge detection in images using Marr-Hildreth filtering techniques

Journal of Neuroscience Methods, 26 (1988) 75-82 Elsevier

75

NSM 00867

Edge detection in images using Marr-Hildreth filtering techniques T . G . Smith Jr. 1, W.B. M a r k s 2 G . D . L a n g e 1 W . H . Sheriff Jr. 3 a n d E.A. N e a l e 4 Laboratories of I Neurophysiology and 2 Neural Control, ~ Instrumentation and Computer Section, N1NCDS, N1H, Bethesda, MD 20892 (U.S.A.), and 4 Laboratory of Developmental Neurobiologk', NICHD, N1H, Bethesda, M D 20892 ( U. S. A. ) (Received 2 February 1988) (Revised 14 April 1988) (Accepted 5 June 1988)

Key words: Edge detection of neuronal border; Marr-Hildreth Laplacian of Gaussian filter: Binary neuronal silhouette Details of the morphology of light microscope images of horseradish peroxidase labeled mammalian neurons in cell culture were investigated. A modified Marr-Hildreth edge-detecting algorithm was used in an image processor to obtain a continuous border of the labeled neurons. The interior of the border was filled to obtain isolated binary silhouettes of the neurons. These silhouettes can be used for further quantitative studies.

Introduction

Edge detection is frequently used in image enhancement, processing and analysis (Gonzalez and Wintz, 1977; Castleman, 1979; Dudgeon and Mersereau, 1984). There are many strategies for edge detection and most of them involve differentiation or gradient measurements. In our research with cultured mammalian neurons, we have a need to generate binary (black = 0; white = 1) silhouettes from grey scale images of neurons that have been labeled with horseradish peroxidase (Neale et al., 1978). These neurons are surrounded by other structures with edges and are superimposed on background noise. We adapted available image processing software to detect the edges of a neuron and fill the interior of its continuous border in order to generate the binary silhouette. Most of the common edge detection techniques tried produced boundaries with gaps. The one technique

that usually produced a continuous boundary was based on an algorithm developed by Marr and Hildreth (1980).

Materials and Methods

The Marr Hildreth algorithm An essential feature of the Marr-Hildreth method is the assumption that the edges of an object in an image should be associated with regions of maximum gradient in intensity. In one dimension the loci of maxima and minima of slope are found at inflection points. These are points where the first derivative of the function is non-zero while the second derivative vanishes. In two dimensions the equivalent statement is that the Laplacian of the function vanishes. The definition of the Laplacian operator and the previous assertion are summarized in the equation: G( x, y ) =~q~F( x, y ) = 0 2 F / 0 x 2 + 02F/Oy 2 = 0

Correspondence." T.G. Smith Jr., Laboratory of Neurophysiology, NINCDS, Building 36, Room 2C02, NIH, Bethesda, MD 20802, U.S.A.

As can be seen the Laplacian operator (50) produces a non-directional second derivative [G(x,y)] of a two-dimensional image [F(x,y)].

76 The locus of zeros in the resulting function is the locus of maxima or minima in the gradient (the two-dimensional equivalent of slope) or else the gradient itself is zero. The loci of zero crossings, that is, the places where the change in the gradient goes from positive to negative or vice versa are the loci of inflections in the image. These inflections occur at the edges of objects. To implement this operation in a computer we use digital filtering, which is the replacement of each point in an image with a weighted sum of its neighbors. The array of weights is called a 'kernel'. To implement the Laplacian operation we digitally filter with a 3 x 3 kernel whose central row and column has the values - 1 4 - 1, the corners are zero and the sum of the kernel elements is zero. This replaces each picture element (pixel) by the second difference of the image at that point in the x and y directions. Marr and Hildreth recognized that differentiation often generates spurious edges due to the amplification of high frequency noise. They, therefore, smoothed the image. In a computer smoothing is implemented by filtering with a kernel which is often larger than 3 x 3: it must be as large as the wavelength of the lowest spatial frequency that is considered noise. For theoretical reasons Marr and Hildreth chose to smooth with a kernel having a Gaussian profile. They adjusted the amount of smoothing by the size of the array, which is proportional to the standard deviation (bandwidth) of the Gaussian. Since the second difference and smoothing operations are both linear Marr and Hildreth were able to combine them into a single kernel, the Laplacian (second difference) of a Gaussian, which is shaped like a 'Mexican hat' (Fig. 2A). Marr and Hildreth pointed out that the positive central disk surrounded by a negative opposing annulus resembles the excitatory-center inhibitory-surround (lateral inhibition) receptive fields found in m a n y sensory systems. Marr and Hildreth argued that the filtering kernels should span a range of sizes so that the full range of scales and spatial frequencies of the image are sampled, as in the visual system. For our purposes it was not necessary to use their complete algorithm since our light microscope images do not contain the large range of scales

0 0 0-1-1-1-1-1-1-1-1-1 0 0-1-1-2-2-3.3-3-2-20-1 - 1- 2- 3- 3- 3- 3- 3-1- 1-2-3-3-3-2-2-2-3-3-3-2-1-2-3-3-4 0 2 3 2 -1-2-3-3 0 3 910 9 -1-3-3-2 2 9172017 - 1- 3- 3- 2 3 1020242010 -1-3-3-2 2 9172017 -1-2-3-3 0 3 910 9 - 1-2-3-3-4 0 2 3 2 -1-1 - 2- 3- 3- 3-22- 20-1 -1-2-3-3-3-3-3-3-3-2-1-1 0 0 -1-1 -2-2-3-3-3-2-2-1 0 0 0-1 -1-1 -1-1 -1-1

0 0 0 1-1 0 0 3 - 3 - 2,. 1 - 1 0 1- 1 0 -4-3-3-2-! 3 0-3-3-2-1 9 2-2-3-3-1 3 -23- 3- 1 9 2-2-3-3-1 3 0-3-3-2-1 0 -4-3-3-2-1 3- 3- 3 - 2- 1 -1 0 -1 0 0 -1-1 0 0 0

Fig. 1. A 15×15 Marr-Hildreth, "Mexican hat' filter kernel. The filter kernel is the second difference of a Gaussian. See text. found in the 'real world' images. We find that by combining the application of a single 9 x 9 to 17 x 17 filter kernel with a threshold operation* for detecting the zero-crossings, we can construct the edges of a neuronal image. However, this method requires h u m a n intervention to select the correct filter and to establish the position of the interior of the cell. Marr and Hildreth were inter° ested in edge detection in the context of human (and perhaps robotic) vision and could make fewer assumptions about the nature of objects to be detected than we. Our large filtering kernel (Fig. 1) is a wide band-pass filter which captures most of the largescale features of the original image. The zerocrossings, which represent the edges in the image, are recovered by, first, a threshold operation on a suitable 'zero' level (e.g. the median value). This produces black-white borders at the original edges. We used a 3 x 3 second difference (Laplacian) filtering operation to detect this border rather than a first difference because the former is omni-directional. However, it produces a biphasic

* The term' threshold operation', in image processingparlance, has a very specific meaning, viz., if the intensity of a pixel has a value equal to or less than some threshold value, its new value is set to a logical zero (black); otherwise, it is set to a logical one (white). Moreover, image processers have described the operation as 'thresholding' or 'threshold, -ed' (Gonzalez and Wintz, 1977; Castleman, 1979; Dudgeon and Mersereau, 1984).

77 positive-negative impulse at the edges. The positive-going one-pixel-wide edge is recovered by a second thresholding step. The interior of the borders are then filled with a constant intensity value, which is greater than that of the edges, to produce a silhouette of the cell distinct from unfilled b a c k g r o u n d structures. The filled cell is isolated by another threshold operation on the filled region. A final 3 x 3 second difference and thresholding step is used when a border-only image is required, The L a p l a c i a n - G a u s s i a n filter kernels are square arrays with from 3 to 17 columns and rows. The values in the kernels are generated from the radially symmetric function:

H(x,y) = (1

A

FILTER KERNEL

g

IMAGE

C

FILTERED IMAGE

-~.. •..

•

.

/I,~o

. ...: ,;" ~

"~..-"

"V'.

;.,,f...

f •. :

g

FILTIMAG >0

\.:

N

E ~P~CE>O Fig. 2. One-dimensional illustration of our technique for edge detection with a modified Marr-Hildreth algorithm. A: 'Mexican hat' filter. B: intensity profile through part of cell. C: profile of filtered image to produce zero-crossings of B. D: thresholding C on median value to isolate zero-crossings. E: results of 3 x 3 second difference of D to leave only borders.

X2-q--Y2

In the second term, e = 2.718 . . . . and e to the p o w e r - (x 2 + y 2 ) / 2 o 2 is the Gaussian function in two dimensions. This becomes multiplied by the first term when operated on by the Laplacian filter. A n N × N kernel consisted of a square array H(x,y) of values of this function c o m p u t e d for N values of x and y given by x, y = - M . . . . . - 1 , 0 , 1 . . . . . M with

M = (N - 1)/2

and

o = N/(2.1-2.3)

where

N = o d d integer ( 3 - 1 7 )

To minimize c o m p u t e r r o u n d o f f error, this array was multiplied b y 25. This was the m a x i m u m magnitude for the array that did not cause c o m p u tational overflow. These values were then r o u n d e d to integers for insertion into the image processor. Finally, the kernel was modified by integral increments to symmetrically placed elements so that the sum of all the elements was zero. The 15 × 15 kernel is shown in Fig. 1. Fig. 2 is a one-dimensional schematic illustration of the algorithm. Fig. 2 B - E diagrams intensity level profiles along, say, a basal dendrite of the neuron. Fig. 2A shows our filter kernel. Fig. 2B is the profile of the original image. In Fig. 2C, obtained with the Laplacian operation, the zerocrossings are seen to occur at the inflections of the original image. The edges are recovered by the two

subsequent Boolean operations (Fig. 2D, E), which isolate the zero crossings and generate edges-only images, respectively.

Application to microscopic images The source of the images was a conventional inverted light microscope (Nikon Diaphot) operated in the brightfield mode. The preparation consisted of m a m m a l i a n central neurons grown in a m o n o l a y e r cell culture ( R a n s o m et al., 1977). Individual neurons were injected by intracellular iontophoresis with horseradish peroxidase ( H R P ) via a micropipette whose tip penetrated the cell soma. Histochemical processing with hydrogen peroxide and diaminobenzidine (DAB) caused the injected cells to be filled with an opaque reaction p r o d u c t (Neale et al., 1978). The culture was placed on the microscope stage and the labeled cell visualized *. The cine output of the microscope was imaged onto the sensor of a video camera ( R C A C C D T-2800). A video amplifier was used to adjust analog dark levels and contrast. The video signal was digitized by the

* Cell-cultured neurons grow as a monolayer on a non-neuronal cellular substrate and have a total depth of about 50/~ When viewed with a 10 (N.A. 0.3) or 20 (N.A. 0.44) power microscope objective that has a large depth of focus, essentially all of the processes are in focus and, if a cell has been injected with HRP, it appears as a two-dimensional structure with few if any depth features.

78 A / D converter of an R C I Trappix 5500 real-time image processor, which was controlled by a D E C P D P 1 1 / 8 6 minicomputer. The general purpose image processing software employed was written by the T A U C o r p o r a t i o n specifically for the Trappix 5500. Images were stored temporarily in volatile 5500 m e m o r y or permanently on a Winchester disc. The images illustrated in this paper were p h o t o g r a p h e d with a Polaroid camera from the face of a TV monitor. We have measured the resolution of the combined optical, video and image processor c o m p o nents of our system. A slide containing a black and white grating whose centers were separated by 6.25/~ was viewed on a 13 inch m o n o c h r o m e TV monitor. U n d e r these conditions we could easily resolve 2 ft.

Results A n image of a cultured neuron viewed with a brightfield light microscope, at a magnification of about 150 x on the TV sensor face, is shown in Fig. 3A. The labeled n e u r o n consists of an ovoid soma from which a n u m b e r of b r a n c h e d dendrites project. Also in the field of view are a n u m b e r of other cells, some of which are unstained neurons.

The goal of the experiment is to extract a binary profile or silhouette of the labeled neuron and eliminate the other objects in the field. We prefer to work with a negative of the original image, as shown in Fig. 3B (diagrammed in Fig. 2B). The image filtered to produce the zero crossings is shown in Fig. 4A (compare Fig. 2C), and its intensity histogram in Fig. 4B. For each neuron, filter kernels of several widths were tried. The 'correct' (usually a 9 x 9 to 17 x 17 kernel) filter was the one that (1) produced the most continuous borders, and (2) gave the least distortion of the relative sizes of the different parts of the n e u r o n (e.g. soma vs dendrite diameters). This selection is the most subjective and time consuming part of the operation. For example, a single run for one kernel took 3-~5 min. In a few cases the borders were not completely continuous, but this could be easily corrected by an image processing procedure which closed any gaps by drawing a line across the gaps. Thresholding the filtered image at its median intensity value, to isolate the zero crossings*, pro* Note that the zero crossings do not have an intensity value of zero, but have a value near the median level of the image (see Fig. 4B). Images are displayed with an intensity value of 0-255 for an 8-bit image.

I k

.

.

.

.

.

.

Fig. 3. A: bright-field light microscope linage of a mouse spinal cord neuron grown in cell culture and labeled by horseradish peroxidase injection. Magnification approximately x 150. B: negative image of A.

79

Fig. 4. A: zero-crossings image of Fig. 3B produced by 'Mexican hat' filter of Fig. 1. B: intensity histogram of A, with intensity on the abscissa (black on left, white on right) and number of pixels of each intensity on the ordinate. d u c e d the i m a g e in F i g . 5 A ( c o m p a r e with Fig. 2D). F i l t e r i n g of that image with a 3 × 3 L a p l a cian kernel, followed b y a n o t h e r t h r e s h o l d o p e r ation, p r o d u c e d Fig. 5B ( c o m p a r e with Fig. 2E). A s can be seen, the n e u r o n is d e t e c t a b l e in these images b u t there is a great deal of noise p r e s e n t after the t h r e s h o l d p r o c e d u r e . I n part, this is a c o n s e q u e n c e of the fact that, after filtering, the grey levels of all the structures fell into the s a m e

range; that is, the h i s t o g r a m of the filtered image was n a r r o w (Fig. 3B). N o n e t h e l e s s , the b o r d e r s of the n e u r o n a l i m a g e are c o n t i n u o u s a n d w i t h o u t breaks, since after filling within the b o r d e r of the n e u r o n * , it c a n be e x t r a c t e d with a t h r e s h o l d p r o c e d u r e , as can be seen in Fig. 6. This is b e c a u s e

* The operator identifies the inside of the cell with a cursor.

Fig. 5. A: image of thresholding Fig. 4A at the median value of Fig. 4A. B: edges image of A, produced by edge detection with 3 x 3 Laplacian and threshold operation.

80

l Fig. 6. Silhouette image of neuron produced by filling within the continuous border of neuron and thresholding on the intensity value of the filled neuron.

the pixel intensity level of the filled neuron has a value greater than the other parts of the image. Thus, it can be extracted from the noisy multiedged image (Fig. 5B). These results were reproducible with several dozen other neurons examined. In most cases the neurons did not fill the field of view of the microscope-TV camera combination and a single image was sufficient for cell silhouette extraction. A few neurons, however, overfilled the field. Their entire extent was captured as subimages and gathered into a single composite image with the image processor. The resulting reduction may have produced some loss of resolution and perhaps of some of the finest dendritic processes.

Discussion

The extraction of binary images of neurons for analytic purposes has a long history. Most of these techniques which involve computers employ some type of camera lucida system a n d a digitizing pad (Uylings et al., 1986). In this paper, on the other hand, we have developed methods for extracting binary silhouettes or borders of neurons, injected

with HRP, from grey scale images containing unstained neurons and other background material with the use of digital image processing procedures. One of our goals has been to develop a system which achieves the same result as previous methods (Uylings et al., 1986), but to do it in a way that is less time-consuming, involves less investigator judgement and is more automatic. Using image processing techniques, it might be imagined that the contrast between the stained neuron and the other parts of our images was sufficiently large that the neuron could be extracted by a simple threshold procedure without the additional filtering steps. Such was not the case, however, as the intensity values of the fine processes of the neurons were often in the same range as the unwanted objects. After some experimentauon, it became clear that a procedure was required that would capitalize on the continuity of the neuron's various parts~ viz. an edge-detecting scheme. We tried a variety of such schemes and they succeeded in varying degree. While not exhaustive they include the Roberts and Sobel operators, non-maximal suppression, simple Laplacian and various directional derivatives. Except for Marr-Hildreth, however, they all failed in one important regard. They did not produce continuous, unbroken borders necessary for our filling routines. The precise reason(s) why the other techniques failed and M a r r - H i l d r e t h succeeded is (are) unclear. Perhaps it is related to the greater smoothing range of this filter, which would produce comparatively less noise. One of the fundamental problems of image processing is that it has no general theory for assessing performance of image manipulation. The field is very much one of 'cut and try'. Other schemes might work quite well where unbroken borders are not required. It is obvious that the effectiveness of these procedures is due in large part to the lack of any visible internal structure within the neuron. This a consequence of the opacity of the H R P - D A B reaction product. How much of a problem visible internal structures will present with other preparations depends on a number of variables. However, most image processing systems have an 'erase' feature that could remove such visible internal structures.

81

The effectiveness of these procedures in isolating the injected neuron from other cells, even those with processes that crossed the injected neuron, is due to the fact that the density gradients at the border of that neuron were much greater than other structures and because of the continuity of that border. Although there is nothing novel about the basic objective, namely edge detection, the approach is novel in that it was inspired by a computational theory of vision. The Marr-Hildreth 'Mexican hat' filter simulated the lateral inhibition of sensory systems and is realized as the second difference of a Gaussian bandpass filter, where the bandwidth is varied to optimize the trade-off between resolution and noise. By and large, the image processing procedures used by biologists are those derived from applications in physics, astronomy, space and engineering. These are usually adequate but occasionally, as in the present case, an image processing procedure arises from the biological sciences that has been useful to other branches of science (Lunscher and Beddoes, 1986).

Acknowledgements The authors are grateful to Drs. Phillip Nelson, Robert Macdonald and Raymond Pun for perfor-

ruing the intracellular horseradish peroxidase injections of the spinal cord neurons used for this study.

References Castleman, K.R. (1979) Digital Image Processing, PrenticeHall, Englewood, N J, U.S.A. Dudgeon, D.E. and Mersereau, R.M. (1984) Multidimensional Digital Signal Processing, Prentice-Hall, Inc., Englewood Cliffs, N J, U.S.A. Gonzalez, R. and Wintz, P. (1977) Digital hnage Processing. Addison-Wesley, Reading, MA, U.S.A. kunscher, W.H.H.J. and Beddoes, M.P. (1986) Optimal edge detector design l: parameter selection and noise effects. IEEE Transactions on Pattern Analysis and Machine Intelligence, PAM1 8: 164-177. Marr, D. and Hildreth, E. (1980) Theory of Edge Detection, Proc. R. Soc. Lond. Ser. B, 207:187 217. Neale, E.A., Macdonald, R.L. and Nelson, P.G. (1978), lntracellular horseradish peroxidase injection for correlation of light and electron microscopic anatomy with synaptic physiology of cultured mouse spinal cord neurons, Brain Res., 152: 265-282. Ranson, B.R., Neale, E., Henkart, M., Bullock, P.N. and Nelson, P.G. (1977) Mouse spinal cord in cell culture. I. Morphology and intrinsic neuronal electrophysiologic properties, J. Neurophysiol,, 40:1132-1150. Uylings, H.B.M., Ruiz-Marcos, A. and Van Pelt, J. (1986) The metric analysis of three-dimensional dendritic tree patterns: a methodological review, J Neurosci. Methods. 18: 127 151.