Segmentation of MRI and CT Data Using Locally Excitatory
Globally Inhibitory Oscillator Networks

Naeem Shareef

Department of Computer and Information Science
The Ohio State University, Columbus, Ohio 43210, USA
shareef@cis.ohio-state.edu

Introduction

Theoretical considerations state that objects may be represented by neural groups that are bound together by temporal correlation of their firing activities. Recent experiments provide support for this view with evidence of oscillations. One important work shows oscillations in the range 40 - 60Hz and synchronization (phase-locking with zero phase shift) between spatially separate neural groups in the visual cortex [3]. Many models have been proposed that try to simulate the experimental findings but few have addressed the difficult engineering problem of image segmentation. Recently an oscillator network model called LEGION [6] was shown to rapidly synchronize oscillator groups that represent a feature of an object and desynchronize oscillator groups representing different features. As suggested by biological considerations, the model is stimulus-driven and uses only local connections. A global inhibitor is used to desynchronize oscillator groups. The model has been shown to successfully segment multiple objects in network simulations [5][6].

Segmentation is an important preprocessing step for automated processing of commonly used medical image data such as Computed Tomography (CT) and Magnetic Resonance Imaging (MRI). Time consuming manual and/or semi-automatic methods in use today have many drawbacks. In this paper we identify and address the main difficulties in segmenting this kind of data which current methods have not been able to overcome. The first is that the acquisition process produces data with widely varying intensity values, ill-defined or blurry boundaries, and sampling artifacts such as speckled noise and undersampled areas. The second difficulty is the enormous amount of data produced by sampling devices. Data is a stack of 2D images, called volume data, representing material in a 3D volume of space. Typical low resolution data dimensions of 128 x 128 x 128 and 256 x 256 x 256 translate to over 2 - 16 MB of datapoints.

In this paper we describe and simulate the four key aspects of LEGION's computational dynamics with a computer algorithm. Our simple approach produces promising results on CT and MRI of real medical data.

Dynamics of LEGION

LEGION is a collection of relaxation oscillators each constructed by coupling an excitatory unit x and an inhibitory unit y in a feedback loop (Figure 1). Unit x receives input from an input stimulus (I) and neighboring oscillators (S). The behavior of an oscillator is defined mathematically with differential equations with respect to time. An analysis can be found in [5][6]. In this paper we describe LEGION with an XY plot, called a phase diagram, that shows the temporal behavior of oscillators in the network. The pair of values at units x and y at some time is represented by a location (x, y) in the graph. A trail or path in the graph illustrates value changes of an oscillator unit over time. The graph itself is partitioned by a cubic function called the x-nullcline and a sigmoid function called the y-nullcline as shown in Figure 2. The values of input stimulus (I) and local coupling (S) have the effect of vertically translating the x-nullcline. The position of the x-nullcline defines two types of behavior. Figure 2a shows the case where an oscillator produces oscillations. Here a single fixed point is created where the two nullclines intersect. The oscillator travels on a counter-clockwise circular path called a limit cycle that alternates between the left leg of the x-nullcline, called the silent phase, and the right leg, called the active phase. The value at unit x increases and decreases at unit y in the silent phase while the opposite occurs in the active phase. Transitions between phases occur at the two humps of the x-nullcline. An oscillator jumps up at the left knee (LK) from the silent to the active phase and jumps down at the right knee (RK) from the active back to the silent phase. Jumps are instantaneous and are denoted by horizontal lines in Figure 2a. Figure 2b shows that when the the x-nullcline shifts downward two additional fixed points are created on either side of LK. The fixed point left of LK is a stable fixed point and thus prevents oscillators in the silent phase from reaching the LK.

The dynamics of a connected network of oscillators in LEGION is determined by four key components: input stimulus, local coupling, the global inhibitor, and oscillator potential. First, positive input stimulus (I>0) produces the case in Figure 2a, i.e. oscillations, and the negative input stimulus (I<=0) produces the case in Figure 2b, i.e. no oscillations. Thus oscillations are stimulus-driven. Second, two coupled oscillating units synchronize if they are close enough on the limit cycle. Synchronization is defined as two or more oscillators that simultaneously jump up and down. Terman and Wang [5] show that networks of two or more oscillators are able to synchronize through local coupling at an exponential rate. Third, the global inhibitor desynchronizes separate groups using a gating mechanism that allows only one synchronized group into the active phase at a time. When a synchronized group enters the active phase the global inhibitor is excited and in turn sends inhibition to the entire network through S. This causes the x-nullclines of all oscillators to shift down to the case in Figure 2b. Groups in the silent phase are prevented from jumping up. The group in the active phase counteracts the inhibition with local coupling. Inhibition stops when this group jumps down to allow the next group to jump up. Thus groups are distinguished temporally with a phase shift.

The fourth component is the oscillator potential. It was introduced into the model to solve the problem of fragmentation because many tiny regions, e.g. noise, could hinder the networks ability to synchronize and desynchronize groups. It functions as an ON/OFF switch to I and so determines oscillations. From analysis and computer simulations in [7] it was determined that synchronized groups containing at least one oscillator (called a Leader) with enough synchronized neighbors continue to oscillate synchronously in the steady state while all other groups cease oscillations. Essentially, leaders in a group are able to reinforce their potentials in the active phase at each cycle. They are also able to preserve oscillations and synchronization of the group via S. All oscillators in nonleader groups experience a reduction in potential values and will eventually enter the case in Figure2b and thus cease oscillations in the steady state.

Region Growing Algorithm

We apply LEGION for segmentation of two- and three-dimensional grey level image data. LEGION can be applied to image data in a straightforward one-to-one mapping between oscillators and pixels. Coupling between oscillators preserves pixel adjacency. Neighborhood topologies can be varied by simply adding or removing neighborhood connections. I and S are appropriately set in order to extract desired features. For example, to segment binary images (black and white) oscillators which map to black pixels receive positive I and thus will produce oscillations. Assigning positive connection weights to neighboring oscillators such that pixel adjacency is preserved in the image as well as the network produces groups of synchronized oscillators representing connected black pixel regions. Thus, groups of oscillators representing connected black pixel regions will oscillate due to I, group and synchronize due to S, and desynchronize from disconnected regions due to the global inhibitor. Applying this scheme to grey level images is equivalent to global thresholding because I binarizes the grey pixel values. Instead, we apply positive I to all oscillators and set the connection weight between coupled units to the reciprocal of the difference between corresponding pixel values, called the pixel difference test. Thus, two coupled oscillators will group and synchronize if the difference between their pixel intensity values is small.

Computer simulation of LEGION is expensive for even modest sized volume datasets. Instead, we construct an algorithm by mimicking the four model dynamics presented in the previous section. Applying positive I to all oscillators in the network translates to incorporation of all pixels in an image into regions. Grouping via local coupling implies that pairs of neighboring pixels are grouped into the same region if the difference in pixel intensity is small. The global inhibitor implements a tagging mechanism to each synchronized group. The oscillator potential implies that only regions containing at least one leader will remain in the final segmentation.

We construct a new region growing algorithm that contains features from seeded region growing and merging techniques[1][2][4]. The algorithm has three steps: Leader Identification, Recruiting, and Background Collection. Initially, all pixels are unlabeled. Leader pixels are identified as those which have a majority of neighbor pixels that pass the pixel difference test. From a leader, Recruiting grows a region by initially assigning the leader a unique label and then iteratively adding unlabeled boundary pixels that pass the pixel difference test with a neighboring pixel already in the region. After all leaders are grown, the Background Collection step labels all remaining unlabeled pixels as background. Note that during Recruiting a region may engulf other leaders and that the background region may be disjoint. Note that Recruiting ensures connected regions. Thus, regions are iteratively grown from "homogeneous" parts of the image using local intensity information.

Results

We applied our algorithm to CT and MRI data of the human head. First we segmented the CT image shown in Figure 3 (left) of an axial slice of the jaw area. Bone appears as highly contrasted bright white and thus is easily segmented with simple segmentation methods. The result of our algorithm is shown in Figure 3 (right). Each segmented region is color coded and the background region is shown in black. Due to presentation difficulties two spatially separate regions that appear to have the same color are in fact two different segmented regions. The algorithm accurately segments bone and even difficult soft tissue areas.

We then applied our algorithm to MRI data which shows soft tissue such as brain and skin(see Figure 4, Figure 5, and Figure 7). It is more difficult to segment MRI data because objects are non-homogeneous with ill-defined boundaries. Figure 4 shows a 2D MRI image showing a sagittal slice. First we changed the neighborhood kernel size to show its effect on segmentation results. Using 4- , 6- , and 24-neighborhood sizes we noticed that smaller kernel sizes produce many small sized regions. The top right image of Figure 4 shows the 4-neighborhood result where the brain is a collection of many tiny regions. The bottom right image shows the 24-neighborhood result where the brain is segmented intact.

Figure 5 shows the algorithm applied to a noisy 2D MRI image showing an axial view. Notice the artifacts introduced by the sampling device such as the considerable amount of speckled noise in the background and undersampling in the nasal area (top portion of the images). The bottom left image shows the nasal area collected into the background region using a 24-neighborhood kernel. We varied the threshold used in the pixel difference test while holding the kernel size at 24-neighborhood. An adaptive thresholding scheme was used where pixel intensities are mapped to a range of threshold values. Various mapping functions may be used. We use a cubic function in all results shown here. A high threshold value allows for grouping between pixels with larger intensity differences. Figure 5 (top right) shows results using a small threshold. Since the change in pixel values is very large, the Recruiting step has trouble filling the entire brain region. The bottom right image shows the brain region filled with a larger threshold. To test robustness to noise, we reduced the image resolution in Figure 5 by half and show the results in Figure 6. Eventhough image quality is reduced, i.e. even more high frequencies, the algorithm produces similar results to Figure 5.

We slice-by-slice hand segmented an entire volume data set for the brain region. Applying the algorithm to volume data simply involves extending the kernel size to 3D neighborhoods. Figure 7 shows results from a single slice from both methods. The holes in the autosegmented results may be easily removed in a postprocessing step using binary operations. The brain regions for the two 3D segmented brain structures were rendered with volume rendering software and results shown in Figure 8 for both methods.

Summary/Conclusion

In this paper we abstracted a computer algorithm from LEGION and applied it to segment large volume datasets from CT and MRI. The novel approach is simple, robust to image artifacts, and easily scalable. It requires setting only two parameters which directly correspond to segmentation results. We are currently investigating threshold and hierarchical network techniques. Here we show segmentation of grey level values, but LEGION may be used as a computational framework upon which other features may be computed. The architecture of LEGION lends itself to implementation on VLSI chip. Thus, our results using a computer algorithm suggest that LEGION may provide a promising computational model for real-time segmentation.

Acknowledgments

This work was supported in part by ONR grant N00014-93-1-0335 and NSF grants IRI-9211419, IRI-9423312, and CDA-9413962. Volume data was provided by ADARC (Alcohol and Drug Abuse Research Center) at McLean Hospital of the Harvard Medical School. We would like to thank DeLiang Wang, Roni Yagel, and Shannon Campbell from the Department of Computer and Information Science at The Ohio State University, Greg Wiet and Bill Mayr from the Department of Otolaryngology at The Ohio State University, and Don Stredney from The Ohio Supercomputer Center, for participation and contribution to this work.

References

[1] R. Adams and L. Bischof, IEEE Transactions On Pattern Analysis and Machine Intelligence, Vol. 16, No. 6., June 1994.

[2] T. Asano and N. Yokoya, Pattern Recognition, Vol. 14, Nos 1 - 6, pp. 267 - 273, 1981.

[3] C. M. Gray and W. Singer, Proceedings of the National Academy of Science USA , Vol. 86, pp. 1698 - 1702, 1989.

[4] R. M. Haralick and L. G. Shapiro, Computer Graphics and Image Processing , Vol. 29, pp. 100 - 132, 1985.

[5] D. Terman and D. L. Wang, Physica D, Vol. 81, pp. 148 - 176, 1995.

[6] D. L. Wang and D. Terman, IEEE Transactions On Neural Networks, Vol. 6, No. 1, pp. 283 - 286, January 1995.

[7] D. L. Wang and D. Terman, Proc. of World Congress on Neural Networks (WCNN-95), pp. II.521 - II.525, 1995.