Computer simulation of cortical polymaps: A proto-column algorithm

Neural Networks, Vol. 5, pp. 187-206, 1992 Printed in the USA. All rights reserved.

0893-6080/92 $5.00 + .00 Copyright © 1992 Pergamon Prcss Ltd.

INVITED AR TICLE

Computer Simulation of Cortical Polymaps: A Proto-Column Algorithm PIERRE LANDAU AND ERIC SCHWARTZ* NYU MedicalCenter and Courant Institute of MathematicalSciences (Received and accepted 6 September 1991) Abstract--Neo-cortical sensory areas of the vertebrate brain are organized in terms of sensory (topographic) maps of peripheral sense organs. Known for over 40 years, cortical topography has been generally modeled in terms of a continuous map of a peripheral sensory surface onto a cortical surface. However, the details of cortical architecture do not conform to this concept: most, if not all cortical areas are represented by an interlacing of multiple featural or peripheral maps in the form of columns consisting of distinct classes of neural input to the cortex. Since such an architecture is at best piecewise continuous, it is useful to introduce a new term which focuses attention on the many-to-one interlaced nature of neo-cortex. We use the term "polymap" to refer to a cortical area which consists of more than one input system, interlaced in a globally topographic, but locally columnar fashion. The best known example of a cortical polymap is provided by the ocular dominance column system in layer IV of primate striate cortex, but the puff~extra-puff and orientation systems of surrounding layers also illustrate this concept, as do the thick-thin-interstripe columns of V-2, and the direction columns of MT. Since polymap architecture seems to be a common architectural pattern in neo-cortex, we have addressed the problem of computational modeling of polymap systems. In the present paper, we present an algorithm, based on the computational geometric constructs of Generalized Voronoi Polygon and Medial Axis, along with related image warps, which is capable of providing a general method for simulating polymap systems. We illustrate this idea with the ocular dominance column system of V-l, and show computer simulations of the structure of binocular stimuli, as they would appear at the level of layer IV in V-1. Comparison of these results to recent binocular 2DG experiments is presented, and methods of generalizing these techniques to other common polymap cortical areas are outlined. This work provides a demonstration of computational neuroscience. The basic anatomical structure of many neo-cortical areas cannot be understood even on a qualitative level, and certainly cannot be modeled, without the development of a detailed set of geometric structures and algorithms which make clear the eelationship of an experimentally observed polymap to its afferent inputs.

Keywords--Vision, Cortex, Ocular dominance, Stereo, Computational neuroscience, Column, Polymap.

1. INTRODUCTION AND PROBLEM SPECIFICATION--AN IMAGE PARADIGM FOR CORTICAL ARCHITECTURE

amples include the ocular dominance, cytochrome oxidase puff/extra-puff, and orientation column systems of primate striate cortex, the thick-thin-interstripe columns of V-2, and the direction columns of MT. As these architectures represent the merging of information from several histological and/or featural laminae, we use the term "polymap" to describe them, justifying this usage because the conventional terms "mapping" and "topographic mapping" are not sufficiently expressive of the complexity and structure of neo-cortical architecture. The structure of polymaps raise difficult conceptual, experimental, and computational issues. The demonstration of an algorithm for simulating them is of fundamental significance to the understanding,

In this paper, an algorithm is demonstrated for modeling polymap architectures of the cerebral neocortex, where the term "polymap" emphasizes the joint occurence of topographic mapping of multiple submodalities, interlaced in the form of macroscopic patches ("columns") into a single cortical layer. Ex-

*Supported by NIMH#MH45969-01A1 and AFOSR88-0275. Requests for reprints should be sent to Pierre Landau,Computational Neuroscience Laboratory,NYU MedicalCenter, 550 First Ave., New York, NY 10016. 187

188

description, and modeling of neo-cortical functional architectures. It is necessary to introduce two new data structures, motivated by work in computational geometry, to construct a detailed algorithm for simulating polymaps: • Cortical proto-column: We introduce the use of Voronoi polygons and Generalized Voronoi polygons as a means of defining the "preimages" of observed columns, which we term "proto-columns." A preliminary description of these ideas has appeared previously (Schwartz, Merker, Wolfson, & Shaw, 1988). The proto-column of a given column establishes the region from which to map into the observed column. • Medial and dual-Medial contours: The medial axis or "skeleton" of a polygon is obtained from repeatedly applying morphological thinning operations to the polygon (Blum, 1967). We define the medial contour to be the bounding contour of the medial axis. The dual medial axis is defined in an analogous way for the "outside" of the polygon. These constructs allow us to deal with the irregular (nonconvex) nature of proto-columns and columns. By generalizing the well known medial axis construction they provide a means for constructing an image warp of the proto-column domain into the column range. The proto-column algorithm has been constructed to allow for a fully general columnar topology. Multiply connected columnar patterns are handled with no extra difficulty ~. A technical appendix provides a detailed discussion of the issues related to generalized Voronoi polygons, medial and dual-medial contours, which are the principal geometric constructions involved in this work. The main body of the text presents a discursive treatment of them, along with image simulations of high resolution scenes mapped into "binocular" visual cortex, which are compared to a recent binocular 2-deoxyglucose (2DG) experiment in primate V-1. We also have measured, via computer simulation, the mean and variance of the "compression" of the within-column map function of V1, and show that the considerable range of variance is consistent with available measurements. After presenting our algorithm, and its application to the case of the ocular dominance column system, applications of it to other cortical columnar systems

~An example of nonsimply connected columnar architecture is the cytochrome oxidase puff/extra-puff system of V1. The extrapuff regions, which contain "holes" where the "puffs" reside, are multiply connected. We have also observed isolated examples of multiply connected ocular dominance columns in our own data.

P. Landau and E. Schwartz

are outlined: V-1 orientation columns and "puffs," the thick-thin-interstripe system of V-2, and MT direction columns. 2. CORTICAL POLYMAP AS NEURAL IMAGE Implicit in this work is the notion of "neural image." Our goal is to create a computer graphic simulation of polymaps at a level of scale at which these architectures may be viewed in a quasi-continuous sense. One way to justify this is to think in terms of a 2DG or cytochrome oxidase experiment: in both cases, the results are an image. Another way is to think in terms of the average local pulse density of firing, in response to an external patterned stimulus. This approach to cortical architecture is motivated by the following observations: • Most current experimental knowledge of cortical architectures is expressed, both descriptively and quantitatively, in terms of "images." The 2-deoxyglucose technique, which has been of major importance in elucidating the architecture of visual cortex, can produce an autoradiographic "image" of a cortical map (Schwartz, Christman, & Wolf, 1981; Schwartz, C h r i s t m a n , & Wolf, 1983; Schwartz, Munsif, & Albright, 1989a, Tootell, Silverman, Switkes, & de Valois, 1982), or of the columnar representation of a binocular stimulus (Tootell, Hamilton Silverman, & Switkes, 1988a). Similarly, the cytochrome oxidase technique (Wong-Riley, 1979) provides a detailed image of the ocular dominance column pattern. Recently, our lab has constructed computer flattened models of macaque ocular dominance columns using this technique (Schwartz et al., 1988). This data was used to construct an image simulation of the ocular dominance column pattern (Rojer & Schwartz, 1990a), which is used in the present paper. • Pragmatically, it is not feasible (at present) to model an entire cortical structure except in "image" terms. Just as the notion of "neural image" is fundamental to our experimental understanding of cortical architectures, it is also fundamental to large scale computer modeling of the cortex. A macroscopic amount of visual cortex contains hundreds of millions of neurons, and it is hardly likely that current computer resources are adequate to modeling such structures at the (single neuron) network level. Although some aspects of columnar architecture (such as multi-column connections to single neurons (Katz, Gilbert, & Wiesel, 1989)) cannot easily be represented within this paradigm, there is a large amount of experimental data which is fully consistent with it and justifies its use in the present work.

Computer Simulation of Cortical Polymaps • Computationally, viewing cortical architectures in terms of quasi-continuous, or map representations, has been a fruitful source of ideas in the context of novel architectures for machine vision. The notion that "cortical images" were of computational significance originated with the Gestalt psychologists. Kohler (1969) concretized this notion, expressing it as a distribution of electric field. In modern terms, his interpretation represents a misplaced concreteness which is best replaced with the notion that topographic mapping in the brain itself is sufficient to provide an "image like" representation with potential computational function (Schwartz, 1977b). This idea has been elaborated in terms of the notion of "computational anatomy" (Schwartz, 1980), and in more recent work along similar lines (Mallot, 1987; Mallot, von Seelen, & Gianannakopoulos, 1990). From a purely technological point of view several research groups are currently in the process of fabricating VLSI sensors based on the complex logarithmic map model for V-1 (van der Spiegel et al., 1989; Wallace, Bederson, Ong, & Schwartz, in press). Other work, in which a simplified model of the ocular dominance column pattern was used as a basis for an efficient stereo algorithm in machine vision (Yeshurun & Schwartz, 1989a) has recently been implemented in a real time vision system with excellent results (Ballard, Becker, & Brown, 1988). Approaches to vision based on space-variant, active vision (Balo och & Waxman, 1991; Burt, 1988; Dickmanns, 1989; Rojer & Schwartz, 1990b; Yeshurun & Schwartz, 1989b) are becoming increasingly popular in machine vision. These observations suggest that a better understanding of the combination of columns and maps in the brains of vertebrates, from the point of view of workers in real time machine vision, is worthy of additional attention.

189

3.1. Regular Mappings of 2-D Surfaces into 2-D Surfaces

general approach to the structure of 2-D mapping, in the context of cortical topography, has been discussed in earlier work (Schwartz, 1984), based on the properties of the Jacobian of a 2-D mapping. This earlier discussion will be summarized here, since it provides a general analysis of 2-D mappings. The appropriate mathematical concept in this case is that of "regular map," that is a map (R2 --~ R2) whose Jacobian (J here) is nonzero and bounded. Regular maps are one-to-one, continuous, and invertible. As an example, the topographic map of the retina onto the cortex (ignoring any columnar substructure, for the moment) is usually conceptualized by anatomists as a regular map. This is a more formal way of stating the "image" paradigm for cortical maps mentioned earlier in this paper. One unfortunate confusion in the anatomical and physiological literature is that the "magnification factor" of a neural map has universally been substituted for the Jacobian of the map, as a means of characterizing its differential structure. Typically, an experimenter will move an electrode, and plot the distance in the cortex which corresponds to a small distance at the surface of the retina. The ratio of these two quantities is called the "magnification factor" (Daniel & Whitteridge, 1961). It is clear that this quantity is (ideally) a differential measurement, and is a scalar. However, the differential structure of a regular map (of R ~ ~ R 2) is represented by a 2 x 2 matrix of partial derivatives, the Jacobian matrix. An earlier analysis of this material (Schwartz, 1984) is summarized in Appendix A, which demonstrates that the differential structure of a regular mapping is the sum of three geometric invariant components: a rotation, a dilation, and a shear. The shear component is a compression along one axis and an expansion of equal magnitude along another axis, with no change in area. In the case that the mapping has no shear, the differential structure of the mapping is characterized by only two numbers: the rotation and the dilation. The dilation is the "magnification factor," and this is the single number which usually emerges from physiological or anatomical experimental work. In the case that the shear component of the map is zero, the map is conformal (equivalently isotropic): its differential expansion is the same in all directions (Ahlfors, 1966). There is good reason to believe that the V-1 map of primates is globally conformal (Dow, Snyder, Vautin, & Bauer, 1985; Schwartz, 1977a; Tootell et al., 1982; van Essen et al., 1984), although several authors have found small anisotropies (Too-

The experimental literature of cortical topography has been dominated by the concept of "magnification factor" (Daniel & Whitteridge, 1961). This scalar descriptor is not adequate to characterize an arbitrary 2-D mapping, unless it is conformaF. A more

2Conformal mapping is in fact a good approximation to the global map for primate cortical area V-1 (Schwartz, 1980; Schwartz et al., 1989a).

The ability to construct realistic and detailed simulations of cortical polymaps thus provides basic insight into both the nature of cortical architecture in a variety of neural systems, and into design principles for high performance machine vision systems.

3. COMPUTATIONAL DESCRIPTION OF POLYMAP SYSTEMS

190

tell et al., 1982; van Essen, Newsome, & Maunseli, 1984). At the present time, it is believed that the maps of cat V-I and V-2 and primate V-2, have considerable shear (i.e., are nonconformal) (Mailot, 1987: Mallot et al., 1990). As seen in Appendix A, the "compression" or anisotropy of a map is expressed as a simple function of the dilation and shear components of the Jacobian. It is the "compression" which has been used in the experimental literature. Thus, (LeVay, Hubel, & Weisel, 1975) measured a 2:1 compression (perpendicular: parallel to the local column boundary) within a single ocular dominance column. This is the data that is modeled in detail in this paper. More recently, (Tootell et al., 1982; Tootell, Switkes, Silverman, & Hamilton, 1988b) have argued for a "compression" less than 2:1. Later in this paper, a modification of the present algorithm is presented, which deals with compressions different, on the average, from 2:1. In the event that the shear component of a map is approximately zero, powerful numerical tools can be brought to bear in order to model the mapping. For example, it can be shown that for a conformal mapping, the details of the entire map are determined by the boundary shape (the boundary of V-1 in a flattened model) S, the boundary of the (hemispherical) retina, and a single internal point. This somewhat counterintuitive idea is a restatement of the Riemann mapping theorem, and has been used to computationally reconstruct the full cortical map function from a cytochrome oxidase experiment which provided only the cortical boundary and the location of the representation of the blind spot in the cortex (Weinshall & Schwartz, 1987). This conformal approximation provided excellent agreement with independent measurements of cortical magnification factor(see (Frederick & Schwartz, 1990) for a complete discussion of the algorithms, implementations, and simulations asociated with this work, and see Figure 6). The assumption of a regular mapping to represent cortical topography is justified by the small size of the ocular dominance columns (and other related columnar systems: puffs, orientation columns). These column systems are on the order of 1 mm -~, and the surface area of V-1 in monkeys is about 1000 mm-'. The individual column structure is only evident on a scale that is about 0.1% of the global map scale. Using conventional electrophysiological techniques, it is very difficult to observe the columnar structure at all. Classical topographic map measurements in el-

'In other work (Schwartz et al., 1988: Schwartz, Shaw, & Wolfson, 1989b), it is shown that it is possible to llatten monkey V-I with negligeable error. For technical reasons it is much easier to work with flattened cortical models than with three dimensional ones. and it is shown that negligible error is introduced in the (digital) flattening process.

P. Landau and E. Schwartz

fect average out the small scale "grain" of the columnar systems, and thus justify the continuum model of a "regular mapping." For many purposes, it is sufficient to conceptualize the global map architecture of cortex independently from the local columnar architecture. But this becomes impossible if one's goal is to model binocular "image" experiments in cortex, or, in fact, if one wishes to attain a deeper understanding of the nature of cortical mappings. In the case of polymaps, the columns pertaining to a single receptor type form a piecewise regular map; the set of associated proto-columns form a regular map. The transformations described below allow us to characterize the piecewise nature of the columns' regular mapping.

3.2. Specific Map Model The map model used in this paper is a numerical conformal mapping, constructed following the description in (Frederick & Schwartz, 1990). This mapping is qualitatively similar to earlier complex log approximations (Dow et al., 1981; Dow et al., 1985; Schwartz, 1977b: Schwartz, 1980), but is based directly on a series of 2-DG topography experiments conducted in our laboratory (Schwartz et al., 1989a). Numerical conformal mappings of this kind provide a good model of primate V-1 topography. Figure 6 shows an example of a high resolution visual scene mapped to the surface of a digitally flattened striate cortex, constructed using these methods (Frederick & Schwartz, 1990). It is important to emphasize that the remainder of the present paper is not dependent on the details of topographic mapping assumed, since the principal focus of the present work is oriented towards to modeling the effects of multiple columnar maps, constrained by a given observed columnar pattern. It is necessary to start with some regular map, as described above, and we have chosen to use a numerical conformal mapping. Any other regular mapping function could be substituted in the following discussion without any loss of generality4. 'Although some care has been taken to develop a general map model, it should be emphasized that the results of this paper are independent of the details of the particular regular mapping chosen. In the case of primate V-l, there are relatively few choices for a map function consistent with current data which are sufficiently well characterized mathematically for computational modeling purposes. For example, (Tootell et al., 1982) do not provide any two dimensional map function. (van Essen et al., 1984) do provide two dimensional reconstructions of V-I topographic data, but do not state the details of the mapping that was used, (Dow et al., 1985) uses a complex logarithmic mapping which is essentially identical to the one used here. Despite a large amount of experimental work in this area, there are surprisingly few viable two-dimensional map functions that have so far been published in a computationally usable form other than the complex log mapping and the numerical conformal mapping generalization of it that is used in this paper.

Computer Simulation of Cortical Polymaps

191

3.3. Proto-columns and Protomaps We will now introduce a columnar aspect to cortical architecture, in the form of the ocular dominance column pattern. Naturally, this adds a considerable degree of complexity to the mathematical characterization of visual cortical architecture. A simple notion of regular mapping of R 2 --9 R 2 is no longer adequate. In order to proceed with the more complex notion of a mapping which interlaces in the form of a columnar pattern, two new concepts are needed: protomap and proto-column. Following the previous discussion, assume access to some regular mapping function which provides the topographic representations of the left and right retinas. These two mappings must be interlaced, somehow, into the observed pattern of columns, which, for simplicity, are taken to exist in layer IV of primary visual cortex. In addition, it is assumed that a particular columnar data set is available to constrain this problem. Figure 1 shows a recent, digitally flattened reconstruction of the full ocular dominance column pattern of macaque V-1 obtained from a oneeyed monkey (one month following enucleation) (Schwartz et al, 1988). A purely synthetic column model can be constructed instead, using a band-pass filtered white noise technique recently described by (Rojer & Schwartz, 1990a). The synthetic column pattern is virtually identical to the " r e a l " column pattern (see Figure 2). For reasons of simplicity, in the following analysis synthetic ocular dominance columns are used. Once again, as in the case of a choice of map functions, it is not significant at the algorithmic level which particular columnar pattern is used. From this point on, we consider a "columnar pattern" to be a binary array of pixels, with black pixels corresponding to the representation of one eye, and white pixels to the other eye -~. No further assumptions are made, concerning the topology of these " c o l u m n s . " They need not be convex; they need not be simply connected (i.e., they may have arbitrary Euler number°). Our only assumption is that a single

planar domain (in this case V-l, layer IV) has been tesselated by a pattern of " c o l u m n s . " It is implicit in this analysis that there is no overlap of "columns," and that all pixels in the planar area belong to either one or the other of the "black" or "white" systems. In other words, consider a model in which each "neuron" (or pixel) belongs to one column. For the example of ocular dominance columns, this column membership index is binary. In the discussion section, we indicate how this may be generalized to more than two types of column, and also discuss generalizations to allow for "fuzzy" membership in a column. At this point, the problem being addressed is of the class of "image warping" (Wolberg, 1990). The goal of simulting a binocular cortex, given a binocular input "retinal" image requires a "warp" of the incoming image data from the left and right eyes into the regions defined by their respective columns ("black" and "white" pixels). We wish to have a name to attach to the "incoming" data from the left and right eyes, prior to its interlacing in columnar cortex. In earlier work these images have been called " p r o t o m a p s " : they are the topographic representation of the data prior to its mixing in the columns of visual cortex (Schwartz et al., 1988). This is a slightly difficult concept, since

qn the case where there are more than two classes of the columnar system being modeled, such as with orientation columns, consider black to represent the current class, and white the union of all other classes. "Some definitions: a connected component of a set is a subset of maximal size such that any two of its points can be joined by a connected curve lying entirely within the subset. In a discrete (2D) set. it is important to use 4-connectivity for figure and 8connectivity for ground (or vice-versa) to avoid ambiguity (Duda & Hart, 1973). The number of holes in a figure is the number of connected components in the complement of the figure minus one. Letting C be the number of connected components in a figure and H be the number of holes, the Euler number E is defined to be C - H. The Euler number is a measure of connectivity: simply connected figures have Euler number 1; an annulus has Euler number 0, and so on.

FIGURE 2. Synthetic column patterns: (a) represents a section from a real ocular dominance pattern. From its power spectrum, we derive the parameters of an anisotropic filter which is applied to white noise and thresholded to obtain (b), a synthetic column image (from Rojer & Schwartz, 1990a).

FIGURE 1. Pattern of cortical ocular dominance columns, hand-traced from a reconstructed image.

192

these "protomaps" only exist in an implicit form. There are no real protomaps in primate cortex, only the columnar mappings of layer IV which are already interlaced. However, some biological motivation for this concept is provided by the developmental studies of (Rakic, 1976; Rakic, 1988). Rakic showed that prior to three weeks before birth, populations of neurons preferring one eye or the other coexist in an intermixed fashion on the developing visual cortex. Shortly before birth, these populations segregate, and ocular dominance columns form. The biological existence of "protomaps," in a prenatal stage of development, has no direct importance to the present algorithm. They provide a regular map from which the actual columnar map of V-1 can be constructed. Now consider two protomaps in superposition over an actual cortical columnar pattern such as the one shown in Figure 1. Our final problem is to effect an image warp from each protomap into its respective columnar targets. This is to be accomplished in such a way that every pixel in both protomaps is represented in either the black or white columnar target regions, such that there is no overlap, and such that every pixel in a columnar region has a source in one of the two protomaps. In order to proceed, the regions in each protomap that correspond to each column in the actual cortex must be constructed (see Figure 3). In other words, each ocular dominance column has a "preimage" in its corresponding protomap. This "preimage" is a columnar region that is somewhat larger than the real column, and which maps into the real column in such a way as to form the observed columnar cortex (see Figure 5). Once again, proto-columns do not necessarily have an independent existence, either in prenatal development, or in the mature cortex. However, they are the crucial conceptual link which allows us to connect two simple topographic mappings (the protomaps) into a single columnar cortex.

FIGURE 3. Pattern of proto-columns generated by the ocular dominance columns In Figure 1. These correspond to the Generalized Voronoi polygons of each column.

P. Landau and E. Schwartz

One way of clarifying the concept of a protocolumn is to describe a Gedanken experiment which might be used to measure them. This description is based on a pilot experiment that was performed several years ago in collaboration with Mikhail Pankrotov of the Retina Foundation, Boston, M A 7. In that work, a high powered opthalmic laser was used to create small photolesions in the retina of a Macaque monkey. These lesions are highly controllable with respect to retina position and size, as the fundus of the retina is directly visualized as the laser is aimed and fired, and the laser power determines the size of the lesion. A pattern of 50 micron lesions were created in the retina, and about one month later processed the brain with the cytochrome oxidase technique. In addition, a retinal whole mount of the eye which had the laser lesions was made. The pattern of spots from the laser was clearly visible in the retinal whole mount, and the corresponding degeneration was visible in the cortex as a localized region of depressed cytochrome oxidase activity. Now, imagine that a series of spots was made in the retina, and each spot localized in the cortex. If this were done with very small spots in a high density grid, and the corresponding spots located in cortex, then the locus of points surrounding a single cortical ocular dominance column would correspond to a retinal area which is termed the retinal proto-column. A simple global topographic mappng (e.g., the complex log, or a numerical conformal mapping) of this measured retinal proto-column would mark out a locus of points in the cortex which surrounded the original cortical column and was perhaps twice the "width" of it. This is what we term the cortical proto-column: it is the area of "recruitment" of an actual cortical ocular dominance column from a hypothetical monocular cortical "protomap." Until such time as this type of experiment is pursued on a large scale, it would seem that experimental data on the details of local polymap architecture in the cortex will be very hard to obtain. Therefore, for the purposes of the present paper, an algorithm has been constructed which produces a plausible set of protocolumns, based on a notion of minimal effort: it is assumed that the relationship is such that a minimal distance condition exists between a proto-column and its column, with the added constraint that the entire protomap is covered, as a jigsaw puzzle, by a tesselation into proto-columns. Later in this paper, simple ways are suggested to allow overlap or "fuzzy" columns. For now, we focus attention on a precise mapping of the entire (set) of proto-maps, into the observed columnar map. A convenient way of constructing such a set of proto-

7Merker, Pankrotov, and Schwartz, 1984, (unpublished).

Computer Simulation of Cortical Polymaps columns which satisfies these constraints is to consider the Generalized Voronoi Diagram of the observed column system. Voranoi Diagram of a Point Set (Preparata & Shamos, 1985). The Voronoi diagram of a finite set S of points in the plane is a partition of the plane such that each element of the partition is the locus of points which are closer to one member of S than to any other member. Each such region is the Voronoi polygon of the corresponding point (see Figure 9). The Voronoi Diagram consists of the boundaries of the Voronoi polygons; it is formed by the intersection of bisectors of the line segments connecting the points in the given set. A practical O(n log n) algorithm for finding the Voronoi diagram of a point set has been recently developed by (Fortune, 1987), although in our current implementation, a simpler "brute force" O(n2) Voronoi algorithm, which is adequately fast for the present purposes, is used. Generalized Voronoi Polygon (Preparata & Shamos, 1985). The Generalized Voronoi polygon V(T) of a subset T of points is the locus of points p such that there exists a point in T which is nearer to p than is any point not in T. A single observed cortical column is associated with T (i.e., corresponds to the pixels of a given column). V(T) is then a larger region which contains T, and has the property that it is the locus points which are closest to T (see Figs. 3 and 4). V(T) is called the proto-column corresponding to T. Two points should be emphasized: 1. The proto-column of a given column is a function of the neighboring columns, not just its corresponding column. 2. The union of the proto-columns of all of the columns provides a perfect tesselation of the protomap surface: every point in an observed column belongs to one and only one proto-column, the proto-columns are "nearest" to their corresponding columns, and the union of the proto-columns provides a perfect covering (jigsaw-puzzle) of the protomap.

193 Given these definitions, we can now provide a definition for the proto-column set of a given pattern of columns: it is the set of Generalized Voronoi Polygons generated by the columns. This definition provides a complete tesselation of the protomap, contingent on the details of the observed columns, in such a way that for each point in a proto-column, the point is closer to the contour of the corresponding column than to any other column contour. Figure 3 shows the proto-columns corresponding to the left eye columns of primate V-1. The Voronoi Diagram and its generalizations are fundamental to the computational geometry of patterns. They have received alternate names in the literature of different fields: in mathematics, they are called Dirichlet regions, in crystallography they are called Wigner-Seitz cells. They provide an elegant and computationally powerful means of classifying the regions associated with cortical columns. One additional property that they provide the medial axis, which is needed in order to construct an image warp of the highly irregular (generally nonconvex and possibly multiply connected) proto-column onto its equally complex column. The Voronoi Diagram of a point set is closely related to the problem of constructing the "skeleton," or medial axis (Blum, 1967) of the point set. The medial axis of a polygon P is the locus or of its internal points such that each p e or is equidistant from at least two line segments on the boundary of P (Preparata & Shamos, 1985). T h e medial axis is frequently visualized as the result of a "grass-fire" set on the boundary of the polygon, which propagates inwards. The medial axis or skeleton of the polygon is the locus of points where the front of the fire meets. This construction provides a simple way of characterizing the symmetry axis of an irregular polygon. For the purposes of the present work, it is necessary to generalize the medial axis in two ways. First, we define a medial contour as the boundary of the medial axis. This provides a closed contour which just surrounds the medial axis, and is necessary

FIGURE 4. A synthetic pattern of ocular dominance columns. (a) their corresponding generalized Voronol polygons, each shaded differently (b). (c) portrays the medial and dual-medial contours El. In this case each column consists of a single annulus, as there are no embedded holes.

194

FIGURE 5. Texture mapping using the proto-columns of the previous figure.

to perform the image warps which are required to map the proto-columns into their columns. In addition, to deal with the case of multiply connected columns, we define the "dual" medial axis, which, in a fairly obvious way, is the "outside" counterpart of the medial axis contour (see Figs. 12 and 13). Having determined the medial and dual-medial contours corresponding to a column/proto-column pair, the final step is to perform an image warp. This warp is effected by a piece-wise bilinear mapping, using the medial contour to provide symmetry. This is illustrated in Figure 5. Performing this for every column/proto-column pair in the cortex produces a binocular image model of that cortex. Thus, several constructions have been used to provide a means of warping the proto-column into each corresponding column in the simplest and most symmetric possible way, given the observed column shape. There is a considerable amount of mathematical baggage associated with this work, but it is unavoidable to define these constructs in order to achieve a fully general and correct proto-column algorithm. The reader interested in implementing this algorithm, or in fully understanding it, is advised to read Appendix B.

3.4. Image Warp The definitions of Voronoi diagram, Generalized Voronoi diagram, medial and dual-medial contour provide a means of constructing a proto-column to column mapping which represents a minimal disruption of the protomap(s), and which produces a symmetric image warp of the resulting proto-columns into columns. We believe that these represent a parsimonious set of assumptions about the detailed nature of columnar topography, given the current paucity of relevant experimental data. In fact, the only attempt to directly assess the nature of this image warp that we are aware consists of a remark in a paper of (LeVay et al., 1975), in which it was claimed that the magnification within an ocular dominance column is twice as large in the direction parallel to the column boundary as in the perpendicular

P. Landau and E. Schwartz

direction. Our algorithm produces such a two-fold compression. This seems to be a reasonable observation, since it suggests that the left and right eye representations, which occupy twice the area of the actual observed cortex, must be "squeezed" by a factor of two in order to "fit" into the columnar cortex. In the discussion section of this paper, we will address other suggestions that have been made concerning the ."compression" of the cortical map within a single ocular dominance column (Tootel, Silverman, & DeValois, 1981; Tootell et al., 1988b), and outline a modification of the proto-column algorithm to deal with this possibility.

4. EXAMPLE OF PROTO-COLUMN ALGORITHM USING THE V-I OCULAR DOMINANCE COLUMN SYSTEM To model the cortical image produced on layer IV of V-1 by a binocular input: • plot the contour of the visual cortex (Obtained from digital or physical flattening of anatomical data) • construct a conformal mapping from the retinal half-disk to cortex. (Or use some other choice of regular map function. In this paper, the algorithm described by (Frederick & Schwartz, 1990) is used). • for each column on on the cortex, contruct a generalized voronoi polygon (proto-column) of that column. • for each contour point of each column, construct a matched pair of quadrilaterals to be used for (locally) bilinear image warping. In the first phase, the retinal image is mapped to the cortical image using a numerical conformal mapping technique. This gives the protomap for one eye, and similarly for the input from the other eye. In the second phase, the protomap of each eye's input is warped into the ocular dominance columns proper to that eye. This gives a stepwise continuous image in the overall shape of the cortex. Repeating the process for the right eye input "fills in" the portions of this image which had no left-eye input. Each binary image of a single eye's ocular dominance columns consists of a number of connected components, one for each column. The digitized data which we start with is a sampled representation of the "continuum" cortical model that we have analyzed here. A coarse sampling function can result in aliased representations of columns and their contours, which would confuse our algorithm. Pragmatically, we have found that it is sufficient that the diagonal intersample distance be smaller than half of any intercontour distance (see Appendix B for brief discussion of aliasing in this context).

Computer Simulation of Cortical Polymaps In the following discussion, we refer to both continuous and appropriately sampled discrete column contours interchangeably. The graphical examples are provided for illustrative purposes only, and represent a much coarser sampling than would be used in practice. Multiply-connected columns are decomposed into simpler components (annuli, described in Appendix B), one for each contour in the original column s. Columns with no holes have one component, while those with one hole have two components, one contained within the other (see Figure 14). Each contour of each column has a generalized voronoi polygon delimited by a medial and a dualmedial contour; these two contours impose a structure on the region which will facilitate the warping step. The warping takes advantage of the decomposition of columns into annuli. Each annulus can be considered to have two addressing indices: a radial index, representing distance from its inner contour, and a circumferential index (Wolberg, 1990). To perform the warp, the contour is sampled along its circumferential index, and the normal at each such sample point is used to define a pair of quadrilaterals (see Figure 16). Warping an image from one quadrilateral to another is easily performed by a bilinear transformation, which is computed and applied separately for each quadrilateral pair. Our algorithm proceeds in stages:

Stage 1: Resample the binary image to ensure that our aliasing constraints are satisfied. Trace the contours of the individual columns, generating a list of contour points. In cases where a column is multiply connected, it will generate more than one contour. Stage 2: Collect all the contour points found by the above procedure and compute the Voronoi diagram for these points. Stage 3: For each contour, construct the medial and dual-medial contours. The region bounded by these two is the protoannulus. The portion of the protoannulus interior to the column itself is the annulus. At the end of this stage, one annulus and one protoannulus exists for each contour in the original binary image.

~We have found several multiply connected columns in our ocular dominance column data. It is our impression that other researchers have not found multiply connected ocular dominance columns, but they have also generally used manual, or mixed manual-computerizedmethods of reconstruction. Any manual intervention is highly likely to result in "editing" out any multiply connected columns which might exist. The need to deal with multiply connected columns is clearly required in V-l, since the cytochrome oxidase puff/extra-puff columnar system is multiply connected: the extra-puff regions are all (at least) doubly connected.

195 Stage 4: For each contour point on each contour, find the normal to the contour and its intersections with the median and dual-medial contours, and use them to construct a pair of quadrilaterals, one contained in the protoannulus, the other in the annulus. The resulting set of matched quadrilaterals is used to warp any protoimage into the areas specified by the columns in our original binary image.

4.1. Comparison of Proto-column Model to Existing Data Tootell et al., 1988a, 1988b) shows an example of the 2-DG image produced in layer IV of V-1 by a binocular stimulus. We have simulated this experiment, using the methods described above. A figure (a crosshair) is presented to each eye (see Figure 7). One eye is assumed to be on axis, while the other is rotated by 2° around a point 1° off axis. For each eye, the projection onto the protocortex is computed using the logarithmic conformal map w = log(z + 0.3), which is in good agreement with several recent topographic map studies (Dow et al., 1981; Dow et al., 1985; Schwartz et al., 1989a). To produce this simulation, a synthetic set of ocular dominance columns was generated using the method of (Rojer & Schwartz, 1990a). The appropriate set of proto-columns were then determined. Finally, we texture mapped the image of Tootell et al.'s (1988) visual stimulus into the protomaps and then into the columns of the simulated cortex. The result obtained (see Figure 7) is qualitatively similar to the data presented by (Tootell et al., 1988b), which represents the only available data, to our knowledge, on the global pattern of binocular stimulation in V-1. The ocular dominance column system introduces a characteristic modulation, or ripple, onto the binocular image. The amplitude of this ripple is determined by the disparity, and the spatial frequency is determined by the column width. Figure 6 shows a good example of this in the letters of the eye chart.

DISCUSSION

5.1. Compression of CMF by Columns In attempting to measure the cortical magnification factor (CMF), researchers have noted a different measurement parallel to the local column borders from that perpendicular to them. Le Vay et al., (1975) observed a 2:1 compression of magnification factor parallel to the long axis of an ocular dominance column, while others (Tootell et al., 1982) suggest that the compression is closer to 1.5:1. More re-

196

P. Landau and E. Schwartz

'

/,-4

@

,It

g FIGURE 6. (a) and (d) depict the retinal Images on the left and right eyes. There is a diagonal disparity of 10 ° of arc between the two eyes, which is too small to see in the retinal images above. (b) and (e) depict the corresponding cortical images obtained by mapping the right hemifield of each eye into the left cortex; (c) and (f) the corresponding column images, and (g) the binocular image.

cently, these same authors claim that a compression of 1:1 is correct (Tootell et al., 1988b). One possible reason for discrepancies in measurement is variance in the compression itself, due to variance in the width of the ocular dominance columns. This can be estimated with our algorithm, providing a lower bound on the experimentally observed range of compression in V1. A set of synthetic columns was created (Rojer & Schwartz, 1990a), and a checkerboard pattern was warped into them. These synthetic columns are based on spectral estimates of real data, and have a variance in width which is comparable to the real data. This allowed us to make numerous individual measurements of the local differential structure of our mapping. For each check (representing an infinitessimal square), the components of the Jacobian af/Ox, Og/ay, Of/Oy, ag/Ox (see Appendix A) were measured, using the fact that the positions of the corners of the checks in the mapped image were easy to determine. For each of 50 checks the compression ratio (see Appendix A) was calculated, and determined to have a mean 2 = 2.06 and a standard deviation a = 0.41. We would expect a similar result from a comparable physiological ex-

periment which measured the differential structure of the mapping at a large number of points. In other words, it should be expected that experimental measurements of compression should be in the range of 2:1, but with a considerable variance about this value, as seems to be the case. An additional source of this varying estimate is based on the fact that the only direct measurement of within-column compression is that of (LeVay et al., 1975), which reported 2:1 compression, in agreement with the simulation of the present paper. The evidence for 1.5:1 or 1:1 compression, as suggested at various times by Tootell and his collaborators, is based on indirect arguments concerning the fact that the vertical and horizontal meridians are of different lengths, and that the representation of data along equal eccentricity segments are different. It has been pointed out that this reasoning is incorrect (Schwartz, 1985), since the different lengths of the vertical and horizontal meridians are expected, based on the conformal structure of the cortical map. Moreover, if the cortical compression were 1:1, there would be a 100% anisotropy in the measured cortical magnification factor across more than one column. This is

Computer Simulation of Cortical Polymaps

197

a

b

d

c

e

f

FIGURE 7. Stimulus (after Tootell, 1987), and Its mappings. (a) The stimulus extends over the central 2° of the visual field, which is otherwise blank. (b) mapping of the halfplane by the complex function w = Iog(z + 0.3), assuming a central fixation point. (c) as (b), but the eye is rotated by 2° around a point 1° from the optic axis along one of the rays. (d) Mapping of the left eye's view onto the ocular dominance columns, (e) Mapping the binocular views onto the corresponding ocular dominance column patterns. (f) detail from (e) which resembles Tootell's original figure (g).

in disagreement with all microelectrode experiments which have been performed, which find that the primate V-1 map is approximately isotropic at this scale. However, it is a simple matter to generalize the algorithm of this paper to account for within-column compressions which are different than 2:1. The Voronoi construction can be "weighted" by modifying the definition of proximity. While this has little impact on the construction of the proto-columns, it introduces a computationally difficult packing problem. Consider a line segment connecting two points. The set of points on the line segment can be partioned into two classes, based on proximity to one or the other endpoint. This construction gives rise to the traditional Voronoi diagram. Now suppose instead that there is a central region of unassigned points, which is temporarily ignored, and assign to each endpoint only those points which are unequivocally closer to it than to the other endpoint: let AB be the line segment. The Voronoi region for point A is given by VA = {p ~ AB, d(p, A)/d(p, B)< 2}, where 2 = 1 for a Voronoi diagram, but could be some other number 0 < 2 - 1 d{p.A}

d{p.B) r

VA

unassigned

VB

Constructing a generalized Voronoi diagram of the columns, using this new definition of distance, gives rise to a set of proto-columns which no longer tesselate the plane, as there remain regions of points which are unassigned to any columns. This presents difficulties for a two-stage mapping which assumes the first mapping is continuous. The process of fitting the proto-columns obtained back together to form a continuous protomap is in general impossible, but can be approximated by a polygon packing problem, which is solved implicitly if the 2 = 1 Voronoi construction is used, but must be solved computationally in other cases. One final remark should be made concerning the suggestion of within-column compressions that are systematically different from 2: 1. As the previous discussion suggests, this would require a protomap which differed markedly in size and shape from the actual cortex. For the 2 : 1 cortex, the protomap and the cortex are both conformai, and are similar in size and shape. For the 1 : 1 cortex, the protomap is conformal, but the cortex is highly anisotropic, and the cortex and protomap have markedly different sizes and shapes. In developmental terms, this would require a major reorganization of the cortical tissue as the ocular dominance columns form. While not impossible, it seems much more parsimonious for the cortex to follow the 2:1 paradigm, since then, the

198

protomap and cortex are essentially identical, and both have essentially identical global map functions. The 1:1 cortex is very awkward to simulate, and would seem to be very awkward for the developing cortex to follow. Moreover, a 1 : 1 cortex would have large anisotropies in magnification factor everywhere, which is in contradiction to all recent experimental evidence. It remains for further experimental work to settle the detailed nature of cortical anisotropies. For the simple types of anisotropies which are currently believed to exist, the present algorithms are sufficient.

5.2. Generalizing the Algorithm to Other Column Systems The main thrust of the present paper was to simulate the ocular dominance column system of macaque striate cortex. This is a system which is believed to be topographically organized at the columnar level of scale, and which consists of two distinct inputs (left and right eyes) which have a very "sharp" border or demarcation between their respective column systems. There are other polymap systems of interest: the orientation column system and puff/extrapuff systems of V-1 are obvious candidates, as well as systems like the cat ocular dominance column system which has "fuzzy" boundaries between columns. We will now briefly outline methods for generalizing the present work to simulate these types of polymap systems.

5.3. Multiple Column Systems Consider the orientation column system of V-1 (Macaque). Orientation columns may or may not be discrete, but certainly have a multiplicity greater than two. It is common to experimentally determine perhaps 5 or 10 ranges of orientation tuning, and "discrete" boundaries can be constructed (if they do not in fact occur biologically) by thresholding. Corresponding to each orientation column is a protomap. If one considered there to be 10 distinct orientation columns, then 10 protomaps would be constructed, and each protomap would represent retinal information that is spatially "filtered" with the corresponding oriented filter. The proto-columns of each orientation column would be determined as in the previous discussion, and by the nature of that construction, each (orientation) protomap would contain a smooth and complete tesselation of its surface. Similarly, the image warp from the proto-columns of each protomap to the corresponding orientation columns would occur precisely as before. Proceeding in this way, we would have the capability to simulate the "neural" image of the oriented laminae of striate cortex.

P. Landau and E. Schwartz

5.4. Fuzzy Column Membership In our above discussion, our use of Voronoi regions has been motivated by the desire for the proto-columns to tesselate the cortex in a nonoverlapping way. This assumption simplifies the analysis, but is neither physiologically plausible for all column systems, nor necessary for our algorithm. Consider the case of a column with fuzzy edges; populations of neurons in some transition region between one column and the other would intermix. This behavior can be readily modeled by this algorithm by assuming that there exists a threshold value below which a receptor population does not belong to a given column. We call "column" the area enclosed by this isopopulation contour. This gives us a binary image, and proceed as before. The only difference is that now some regions in the output image have multiple predecessors in the protomaps 9. The mixing effect is easily modeled by making the value at each point a linear combination of the values projected to it. Weights could be assigned as a linear function of distance from the column midline. A modification of this sort would be expected to provide a good model of the "fuzzy" boundaries of, for example, the cat ocular dominance column system.

5.5. V-1 Cytochrome Oxidase Puff/Extra-Puff Systems Cytochrome oxidase "puffs" are areas most visible in layer III of the cortex which are noticeable after staining normal cortex with cytochrome oxidase. These areas are believed to contain populations of nonoriented color-opponent cells, again grouped by type (red-green, blue-yellow) (Livingstone & Hubel, 1987b). The extra-puff system, consisting of the cortex with the puff regions deleted, is believed to consist of orientation tuned, largely monochromatic cells and is multiply connected. Following the preceding discussion, we would construct two protomaps, one for nonoriented chromatically sensitive filters, and one map for oriented, nonchromatically sensitive filters. The oriented aspect of this latter map could be handled by constructing a set of oriented protomaps, as outlined above. The corresponding proto-columns would then be constructed in each protomap, now making good use of the fact that our algorithm is topologically general, and the final image warp would then provide a model of the puff/extra-puff polymap.

5.6. V-2 Thick-Thin-lnterstripe Primate V-2 (cortical area 18) is divided into three broad classes of stripes labelled "thick," "thin," and 'q'his, of course, makes the columnmapping noninvertible.

Computer Simulation o f Cortical Polymaps

"pale," (interstripe) according to appearance under cytochrome oxidase staining (Curcio & Harting, 1978; Hubel & Livingstone, 1987; Livingstone & Hubel, 1987a; Tootell, Silverman, DeValois, & Jacobs, 1983). Each of these populations have distinct afferent and efferent projections. The thick stripes, which correspond to the magnocellular subdivision of the lateral geniculate nucleus (LGN), and to layer 4B of V-1 project to the motion-detection area MT. Receptors in this area are sensitive to movement and binocular disparity. The thin stripes, which are connected to the blob systems of V-l, contain cortical neurons concerned with color, but are relatively insensitive to movement, stereopsis, and orientation. The pale stripes correspond to the parvoceilular subdivision of the LGN and the interblob regions of V-1. Receptors in this region have high spatial resolution, orientation specificity, and end-stopping (detect ends of line segments). While it appears difficult, with the current level of knowledge about V-2, to imagine the computational utility of juxtaposing representations of such apparently different measures in a receptotopic manner, it is not difficult to imagine simulating this system in an analogous manner to those we have described above. This could either be done directly from the retinal representation, or from some intermediate representation such as our model of V-1. Once again, access to an observation of the final columnar pattern, coupled with some understanding of the stimulus response properties of the respective protomap subsystems, allows a simulation to be performed. 5.7. MT Direction Columns

One final example is provided by the Middle Temporal cerebral cortex (MT), which appears to be concerned with the processing of motion (Gatass & Gross, 1981; Maunsell & Newsome, 1987; Maunsell & van Essen, 1983a; Maunsell & van Essen, 1983b). This area appears to be organized into columns of receptors sensitive to particular orientations of motion (Albright & Desimone, 1987; Albright, Desimone & Gross, 1984). In direct analogy to our discussion of V-1 orientation columns, this system would be simulated by constructing multiple motion orientation protomaps and proto-columns. 6. SUMMARY

In summary, the algorithms described in the present paper provides a generic means of simulating cortical polymaps, using the data structures of protomap, proto-columns, and (dual) medial contour. Due to the relative lack of current experimental data on the fine details of any one of the cortical systems of in-

199

terest, we find it necessary at the present time to make assumptions of parsimony in order to proceed with this work. In the case of primate ocular dominance columns, we have demonstrated that the assumptions that we have made (minimal distance of proto-columns to columns, simple local bilinear image warp, etc.) are well justified: we have in fact produced a simulation which strongly resembles the only available data in this area. We have also made a quantitative prediction concerning the behavior of the cortical magnification factor in the presence of columns. For the case of the other columnar systems outlined here, we intend to produce the corresponding proto-column simulations in the near future. These simulations will provide clear experimental predictions which are accessible to modern histological techniques (2DG, optical mapping, etc.). This work has provided a link between the observable column patterns of the neo-cortex, and the neural maps which they represent, an important first step in understanding the detailed structure of cortical polymaps. REFERENCES Ahlfors, L. (1966). Complex analysis. New York: McGraw-Hill. Albright, T. D., & Desimone, R. (1987). Local precision of visuotopic organization in the middle temporal area (MT) of the macaque. Experimental Brain Research, 65,582-592. Albright, T. D., Desimone, R., & Gross, C. G. (1984). Columnar organization of directionally selective cells in visual area MT of the macaque. Journal of Neurophysiology, 51, 16-31. Baloch, A. A., & Waxman, A. M. (1991). Visual learning: adaptive expections, and behavioural conditioning of the mobile robot mavin. Neural Networks, 4, 271-302. Ballard, D., Becker, T., & Brown, C. (1988). The rochester robot (Tech. Rep. No.). University of Rochester, Department of Computer Science. Blum, H. (1967). A transformation for extracting new descriptors of shape. In W. Wathen-Dunn (Ed.), Models for the perception of speech and visual form. Cambridge, MA: MIT Press. Burt, P. J. (1988). Algorithms and architectures for smart sensing. Proceedings DARPA Image Understanding Workshop, 139153. Curcio, C. A., & Harting, J. K. (1978). Organization of pulvinar afferents to area 18 in the squirrel monkey: evidence for stripes. Brain Research, 143, 155-161. Daniel, M., & Whitteridge, D. (1961). The representation of the visual field on the cerebral cortex in monkeys. Journal of Physiology, 159, 203-221. Dickmanns, E. D. (1989). Subject-object discrimination in 4Ddynamic scene interpretation for machine vision. Proceedings, Workshop on Visual Motion, 298-304. Dow, B. M., Snyder, A. Z., Vautin, R. G., & Bauer, R. (1981). Magnification factor and receptive field size in foveal striate cortex of monkey. Experimental Brain Research, 44,213-228. Dow, B., Vautin, R. G., & Bauer, R. (1985). The mapping of visual space onto foveal striate cortex in the macaque monkey. Journal of Neuroscience, 5, 890-902. Duda, R. O., & Hart, P. E. (1973). Pattern classification andscene analysis. New York: Wiley. Fortune, S. (1987). A sweepline algorithm for voronoi diagrams.

ACM Symposium on Computational Geometry. Frederick, C., & Schwartz, E. L. (1990, March). Conformal image

200 warping. IEEE Computer Graphics and Applications, pp. 5461. Gatass, K., & Gross, C. G. (1981). Visual topography of striate projection zone (MT) in posterior superior temporal sulcus of the macaque. Journal of Neurophysiology. 46, 621. Hubel, D. H., & Livingstone, M. S. (1987). Segregation of form, color and stereopsis in primate area 18. Journal of Neuroscience, 7, 3378-3415. Katz, L. C., Gilbert, C. D., & Wiesel, T. N, (1989). Local circuits and ocular dominance columns in monkey striate cortex. Journal of Neuroscience, 9, 1389-1399. Kohler, W. (1969). The task of Gestalt psychology. Princeton, N J: Princeton University Press. LeVay, S., Hubel, D. H., & Wiesel, T. N. (1975). The pattern of ocular dominance dominance columns in macaque visual cortex revealed by a reduced silver stain. Journal Comparative Neurology, 159, 559-576. Livingstone, M. S., & Hubel, D. H. (1987a). Connections between layer 4B of area 17 and the thick cytochrome oxidase stripes of area 18 in the squirrel monkey. Journal of Neuroscience, 7, 3371-3377. Livingstone, M. S., & Hubel, D. H. (1987b). Psychophysical evidence for separate channels for the perception of form, color, movement and depth. Journal of Neuroscience. 7, 3416-3468. Mallot, H. A. (1987). Point images, receptive fields, and retinotopic mapping. Trends in Neuroscience, 10(8), 310-311. Mallot, H. A., von Seelen, W., & Giannakopoulos, F. (1990). Neural mapping and space-variant image processing. Neural Networks, 3(3), 245-264. Maunsell, J. H. R., & Newsome, W. T. (1987). Visual processing in monkey extrastriate cortex. Annual Reviews Neuroscience, 10, 363-401. Maunsell, J. H. R., & van Essen, D. C. (1983a). Functional properties of neurons in middle temporal visual area of the macaque monkey. I. selectivity for stimulus direction, speed and orientation. Journal of Neurophysiology, 49, 1127-1147. Maunsell, J. H. R., & van Essen, D. C. (1983b). Functional properties of neurons in middle temporal visual area of the macaque monkey. II. binocular interactions and sensitivity to binocular disparity. Journal of Neurophysiology, 49, 11481167. Prcparata, E P., & Shamos, M. I. (1985). Computational geometry: An introduction. New York: Springer-Verlag. Rakic, P. (1976), Prenatal genesis of connections subserving ocular dominance in the rhesus monkey. Nature, 261,467-471. Rakic, P. (1988). Specification of cerebral cortical areas. Science, 241, 170. Rojer, A., & Schwartz, E. L. (1990a). Cat and monkey cortical columnar patterns modeled by bandpass-filtered 2D white noise. Biological Cybernetics, 62, 381-391. Rojer, A. S., & Schwartz, E. L. (1990b). Design considerations for a space-variant visual sensor with complex-logarithmic geometry, lOth International Conference on Pattern Recognition, 2,278-285. Schwartz, E. L. (1977a). The development of specific visual projections in the monkey and the goldfish: Outline of a geometric theory of receptotopic structure. Journal of Theoretical Biology, 69, 655-685. Schwartz, E. L. (1977b), Spatial mapping in primate sensory projection: analytic structure and relevance to perception. Biological Cybernetics, 25, 181-194. Schwartz, E. L. (1980). Computational anatomy and functional architecture of striate cortex: A spatial mapping approach to perceptual coding. Vision Research, 20, 645-669. Schwartz, E. L., Christman, D. R., & Wolf, A. P. (1981). Positron emission tomography studies of human visual cortex. Society for Neuroscience Abstracts, 7,367. Schwartz, E. L. (1984). Anatomical and physiological correlates

P. Landau and E. Schwartz of human visual perception. IEEE Transactions on Systems, Man and C.vbernetics, SMC-14, 257-271. Schwartz, E. L. (1985). On the mathematical structure of the retinotopic mapping of primate striate cortex. Science, 227, 1066. Schwartz, E. L,, Christman, D. R., & Wolf, A. P. (1983). Human primary visual cortex topography imaged via positron-emission tomography. Brain Research, 104, 104-112. Schwartz, E. L., Merker, B., Wolfson, E., & Shaw, A. (1988). Computational neuroscience: Applications of computer graphics and image processing to two and three dimensional modeling of the functional architecture of visual cortex. IEEE Computer Graphics and Applications, 8(4), 13-28. Schwartz, E. L., Munsif, A., & Albright, T. D. (1989a). The topographic map of macaque VI measured via 3D computer reconstruction of serial sections, numerical flattening of cortex, and conformal image modeling. Investigative Opthalmology Supplement, 298. Schwartz, E. L., Shaw, A., & Wolfson, E. (1989b). A numerical solution to the generalized mapmaker's problem. IEEE Transactions on Pattern Analysis and Machine Intelligence, 11, 10051008. Tootell, R. B., Silverman, M., Switkes, E., & deValois, R. (1982). Deoxyglucose analysis of retinotopic organization in primate striate cortex. Science, 218,902-904. Tootell, R. B,, Silverman, M. S., & DeValois, R. L. (1981). Spatial frequency columns in primary visual cortex. Science, 214, 813-815. Tootell, R. B., Silverman, M. S., DeValois, R. L., & Jacobs, G. H. (1983). Functional organization of the second cortical visual area in primates. Science, 220, 737-739. Tootell, R. B. H., Hamilton, S. L., Silverman, M. S., & Switkes, E. (1988a), Functional anatomy of macaque striate cortex. I. ocular dominance, binocular interactions, and baseline conditions. Journal of Neuroscience, 8, 1531-1568. Tootell, R. B. H., Switkes, R., Silverman, M. S., & Hamilton, S. L. (1988b). Functional anatomy of macaque striate cortex. ii, retinotopic organization. Journal of Neuroscience, 8. van der Spiegel, J., Kreider, F., Claiys, C., Debusschere, I., Sandini, G., Dario, P., Fantini, F., Belluti, P., & Soncini, G. (1989). A foveated retina-like sensor using CCD technology. In C. Mead & M. Ismail, (Eds.), Analog VLSI implementations of neural networks. Boston: Kluwer. van Essen, D. C., Newsome, W. T., & Maunsell, J. H. R. (1984). The visual representation in striate cortex of the macaque monkey: Asymetries, anisotropies, and individual variability.

Vision Research, 24,429-448. Wallace, R., Bederson, B., Ong, P.-W., & Schwartz, E. (in press).

Space variant vision with an MOS complex log sensor. Weinshall, D., & Schwartz, E. L. (1987). A new method for measuring the visuotopic map function of striate cortex: Validation with macaque data and possible extension to meas u r e m e n t of the h u m a n map. Society for Neuroscience Abstracts, 13, 1291. Wolberg, G. (1990). Digital image warping. Los Alamitos, CA: IEEE Computer Society Press. Wong-Riley, M. T. T. (1979). Changes in the visual system of monocularly sutured or enucleated cats demonstrable with cytochrome oxidase histochemistry. Brain Research, 171, 11-28.

Yeshurun, Y., & Schwartz, E. L. (1989a). Cepstral filtering on a columnar image architecture: a fast algorithm for binocular stereo segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 11(7), 759-767. Yeshurun, Y., & Schwartz, E. L. (1989b). Shape description with a space-variant sensor: Algorithms for scan-path, fusion and convergence over multiple scans. IEEE Transactions on Pattern Analysis and Machine Intelligence, 11, 1217-1222.

C o m p u t e r S i m u l a t i o n o f Cortical P o l y m a p s

201

A P P E N D I X A: R E G U L A R M A P P I N G S O F 2-D S U R F A C E S I N T O 2-D S U R F A C E S

For a regular map F:R-" --~ R-'; F:(x,y) -+ (l'(x, y), g(x, y)), the Jacobian is: j

[~l'/0x

OflOv]

As mentioned in earlier work (Schwartz, 1984). it is instructive to rewrite the Jacobian in terms of its symmetric and antisymmetric parts: 2J = S + A = (J + j r ) + ( j _ j r )

It is further useful to rewrite the symmetric part as the sum of a diagonal and a traceless matrix: S = [

21; I; + g,

I;. + g,] = D + T 2g, = [ tO tO] + [ 2 1 ; - tr Ll; + Y,

9f, + g, ] -g,. - tr J

where tr represents the trace of the Jacobian, i.e.. tr = I; + g,. The traceless component of the Jacobian can be diagonalized by an orthogonal transformation, since it is symmetric and positive definite, which we write explicitly:

FIGURE 8. A circular region (annulus) which "grows" outward, and one which "grows" inward,

APPENDIX

The differential structure of a regular mapping has been dissected into the sum of three geometric invariant components: • A rotation (the antisymmetric part A) • A dilation (the diagonal part D) and • A shear (the traceless part T). The shear component is a compression along one axis and an expansion of equal magnitude along another axis, of magnitude equal to 2. The direction of the principal axes are given by the rotation matrix R. In the case that the mapping has no shear, then T (traceless, symmetric part of the Jacobian) is zero and the differential structure of the mapping is characterized by only two numbers: the rotation and the dilation. The expansion is the "ma D nification factor," and this is the single number which usually emerges from physiological or anatomical experimental work. In the literature, the apparent compression of the CMF perpendicular to the axes of ocular dominance columns is sometimes expressed as a ratio of gain along the two axes. In the above analysis, if we consider only the unrotated, symmetric portion RSR ~ = R D R -~ + R T R ' = D + R T R =[tr+). 0 ] 0 tr - 2 ' and examine the ratio of the two nonzero terms, we obtain the "'compression ratio" of the mapping. 7.1. Example

Suppose we consider a square ("differential patch") with coordinates, (0,0), (0,10), (10,0), (10,10). Suppose that in the mapping, we measure the image of this square as having coordinates (0,0), (1,5), (11,5), (10.0). Then the Jacobian

,__ [0,10 o] 1/10

5 0 '

and the symmetric component and its traceless form, respectively, are

s= [o.l OoO;] • "

'

o.<

-0.25J"

Diagonalizing T gives us R, and

[000, 0 :9 1 and the compression ratio derived from the two diagonal elements of this matrix is 2.03:1.

B: C O L U M N S , P R O T O C U L U M S , ANNULI

8.1. Definitions

In order to be able to define a mapping from a column to a protocolumn, we must guarantee that the proto-column is deformable into (has the same Euler number as) its column. To resolve this difficulty, we introduce the notion of an annulus which topologically has one and only one "'hole" in it. We will show that we can decompose each column and its proto-column into a series of such annuli, and can reduce the general column mapping problem to two cases: that of an annulus which grows inward, and that of an annulus which grows outward (see Figure 8). To map the entire protoimage into the corresponding column image, we decompose both into corresponding annuli, and map each annulus in sequence. The contours in our image are the set of Jordan curves u' which form the boundaries of the columns. These contours have been sampled at some spatial sampling frequency and are represented by their respective sets of sample points. Given a set P of N labelled points p ~ P in the plane, the Voronoi polygon" of each p is the locus of points closer to this point than to any other point in S: vp = {q E S I ~p' E P, p' ¢ p, d(q, p') < d(q, p)} where d(x, y) is the Euclidean distance from x to v. The Voronoi diagram is the set of all such polygons t: (see Figure 9). A Voronoi polygon is convex and has no more than N - 1 sides (Preparata & Shamos, 1985). The medial axis, or skeleton, of a polygon (Blum, 1967; Duda & Hart, 1973), is defined in one of two ways: the result of a morphological, "prairie fire" thinning operation, or equivalently, as the loci of proximity of the sides of the containing polygon (Preparata & Shamos, 1985). To construct the medial axis in this latter way, the Voronoi procedure is modified to compute adjacency sets of edges rather than points, where the distance between a point and an edge is simply the shortest distance to the edge

u% Jordan curve is a closed curve in the plane which separates a simply connected interior region from an unbounded exterior region (Preparata & Shamos, 1985). t~Also known as Voronoi region. ~:Properly speaking, the polygons of points on the convex hull of the set of points P extend to infinity, but they are restricted here to an image which is a subset of the plane.

P. L a n d a u and E. S c h w a r t z

202

%

p"

a,.~

qi\ J

0

/ . -

)-~"

o /

/

/

FIGURE 10. Two types of medial axes can be created, one using the Voronoi tesselation due to the vertices of a polygon, the other due to the edges of the polygon. As the number of vertices Increases, these two are approach each other.

p',

qo''\'

o

Vc is a Voronoi polygon associated with the contour point p of contour c. We associate with oc~ the set of edges Ecp = {e,pp,}which delimit the region. We orient these edges in a clockwise direction. Each edge is labeled by the two Voronoi polygons it separates (see Figure 9). Consider the set of edges associated with the c-th contour

Ec = U Ecp - E" = [..J e,.pp, - E',

FIGURE 9. Portion of Voronoi diagram due to some points on contour C~. We see two adjacent points p and p'. The Voronol polygon Ucpis the area delimited by the points {a, b, c, d, el. The edge set { ab, bc, de, e a / w e represent as E~. = {ec,~, e~.p., e~.-, e~op,,e~.¢}.

d(x, E) = inf(d(x, y)), y E E. The Voronoi diagram VE obtained is in a sense the dual of that obtained from the vertices of the polygon V~.. The medial axis of a polygon is the maximal subset of the edges in the Voronoi diagram such that both endpoints of the edge are contained in the polygon (see Figure 10). By analogy, we define the dual-medial axis as the subset of the edges in the Voronoi diagram which are contained in the proto-column (generalized Voronoi polygon) but have both endpoints exterior to the column itself. Thus edges which bisect adjacent contour points (resp. edges) are not part of the medial or dual-medial axes. In the limit where the size of an edge becomes very small and the number of vertices large, the "medial" portion of both Ve and V~ are identical. This is due to the fact that as an angle becomes very acute, its cosine approaches 1, so that the distance between an interior point and its nearest contour vertex is not appreciably different the distance between it and its nearest contour edge. In this work, we used Ve as it is simpler to compute. The term medial contour is a generalization of the medial axis: It is the locus of points on the bounding contour of the medial axis, ordered in a counterclockwise direction ~-~.The contour will not generally be convex, though it is considered to lie arbitrarily close to the medial axis. Analogously we define the dual-medial contour, the bounding contour of the dual-medial axis. We call the region of the plane delimited by these two contours the protoannulus, and the portion of the protoannulus internal to the contour itself the annulus. 8.2. Construction of a Protoannulus Let 6", represent the set of contour points of contour c. These points have been obtained by sampling the contours, as described above. We distinguish two types of contours: external, bounding a column, and internal, bounding a hole. Let C = t_J,C, represent the union of all such sets (i.e., the union of all column and hole boundary point sets). We call V the Voronoi diagram for C, and V, C V the subset of V corresponding to points in Co. Each oep

~In the case where the original column has one or more holes, an additional medial contour will be generated for each such hole, corresponding to the inner hull of the medial axis which is reachable from the hole.

p

p.p

where E" = {e,,p, E E~.Ip, p ' are adjacent contour points on the same contour}, (see Figure 11). The set Ec consists of two concentric oriented contours which delimit a region of the plane in the form of an annulus, though in many cases, the inner contour is degenerate: the annulus has been shrunk such that the hole is no longer visible. DEFINITIONS. In the following, we assume a discrete, isotropic sampiing o f the plane o f sufficiently high frequency not to introduce aliasing errors ~. Let p be a point in contour C~. Then there is no contour point q ECc. ~ Cc on some other contour such that q is closer to p than either o.f p's neighbors in C~. This follows directly from antialiasing constraints we have imposed on the initial sampling. LEMMA. Let ~ be a sampling o f the plane o f high enough frequency to avoid aliasing. Let C~ C ~ be a closed contour o f a column, either interior or exterior. Let {ecpp,} represent the clockwise-oriented set of edges of the Voronoi polygon associated with point p ~ C~. Each such edge separates the Voronoi polygon of point p from that o f pohzt p'. Then the set E~ = {..J e,pp, - E,'-, p,p"

where E~ = {e,.,,. @ E~lp. p' are adjacent contour points on the same contour}

is composed of two closed contours, one interior to (termed medial contour), one exterior to (termed dual medial contour), the closed contour C~. These two concentric contours bound a region o f the plane called a protoannulus. Further, none o f the protoannuli generated by different contours in the plane will overlap, and the set o f all protoannuli will form a partition o f the plane.

t4A necessary and sufficient condition is that the intersample distance be smaller than the smallest inter-contour distance divided by 2V~. This ensures that the distance between adjacent contour points on the same contour is always smaller than the distance between adjacent contours. The samplng frequency must also be high enough not to alias representations of an individual contour boundary. If both these conditions are satisfied, there is never a pixel shared among two contours, and contours never have spurious articulation points.

Computer Simulation o f Cortical Polymaps

203

,. .................

qo i

qo

.k.

oo.~,.

,k.

.o- .

!,

,,,,\, .o - . .

\

q,, '.,

\,

'\\

/ !

o

q •~

\\

o

;/

\\\,

/

/

!

o

FIGURE 11. The two adjacent points p and p' have Voronoi polygons with edge sets E~ and E~o.. In constructing our medial and dual-medial contours, we will form the union of these two sets, but will then subtract e~p, and e~,p, which divide the polygons belonging to adjacent contour points on C¢.

Proof. By (modified) induction. Initial case: C~ is a closed contour; the Voronoi polygon for each point along C~ is convex, and shares an edge with both of its neighbors along the contour. This is always true because they are its closest neighbors. When two Voronoi polygons, each delimited by a clockwise contour, are

juxtaposed, and the shared edges are eliminated, there remains a clockwise path delimiting both regions (see Figure 11). Induction: if there exists a set of adjacent contour points of length n constructed by such a procedure, then the property holds true for a set of length n + 1, unless this closes the contour (termination). The regular case is similar to the initial case. For the termination case, when the last point on the contour is added, closing the contour, the two " e n d " segments of the contours are removed, leaving two closed contours. As all shared boundary elements which are eliminated from E, are perpendicular to the contour C, and necessarily present in the original set, the contours formed will be on either side of the intial contour. The annulus generated cannot overlap any other annulus because, by its nature, the Voronoi diagram is a partition of the plane into discrete regions of points, and any partition of the original partition is still a partition of the plane.

COROLLARY. Any set of multiply connected regions in the plane can be partitioned into a set of annuli, as defined above. Proof. We simply consider each interior and exterior contour of a column separately. Each such contour generates a single annulus, a result which follows from above. 8.3. Construction o f an A n n u l u s

An annulus is the subset of the protoannulus contained in the column proper. It is constructed, similarly to E~ above, from regions ~% = {p E o,v,lp E C,}. This corresponds to considering the portion each Voronoi polygon which is internal to the column.

,\ q o

,f....

~c/.

i--t

o

".\ \....~;, \\

o

E

~'~.

o

....... E

/

/!

"\

\

,.

o

>.,~., //:

I

',

\

o o

,

/.f F ~

,.

,\

o

/ \

o

\\__-.J

[]

FIGURE 12. (a) A simple, convex column (shown with two adjacent columns). This column Is delimited by a single contour. For illustrative purposes, we choose six labelled contour points A, B, C, D, E, F, though In our algorithm, the grid Is much denser. (b) the corresponding Voronoi polygons. (c) the set E~of edges which delimit the proto-column, composing two disjoint sets. termed "dual-medial contour" and "medial contour." (d) the set I~, which delimit the column. The medial contour remains unchanged; the column contour is shown. The medial contour Is effectively degenerate, but for our purposes we treat It as a true contour.

o

. . . . . . . . ,/"" ~\\ o

i

I

"~

/

"l"

",,

I

\

!

/

,

/

. . . . . . J.-----i.--,:,

,.,

.,'/\",

/

I

/ L'f

/

~

~,.~. ",,

/"..

!

I

I

~,..'", / \ /

1

°

'\,

..: v

'~

o

~ , \'

"- --'"".... ""i..... ~,---".:]'~'1)

\.

'~'-.. i

\/"

/

//

\ j..-:.........~.

'.,

I

[] a

..~

\

o

o

o F. . . . . . . . . . . . i ~

o

", ......,,"" "~'~

,.. (

o

[]

..*,

o

~""* "",,. ,/"

,,,,°"

" "-

"~, t . , \ \

o

..)".Q-"J . . . .~. . \ " :>.....

°

~

o

.

,,,j.L ..........! , / []

°

o

[]

FIGURE 13. (a) A column which is no longer convex. (b) the corresponding Voronoi polygons. (c) the set E, of edges which delimit the proto-column, and (d) the set E,, which delimit the column. Note that here part of the dual-medial contour is degenerate. o 0

:

0

JR...

~

o

/

\,

[]

0 .'¢"

'\

o

°

0

'

o

o

0

o

,'

[]

;

.. o

o

°o.'~-.

.

',,.(I.

~.\

0

),

o

o

//..,

o

i-

.j o

o ,,,

o

[]

\

..~- . . . . . . . o

/ []

FIGURE 14. (a) A column with a hole In It. This column has two contours; each will generate a annulus. (b) the corresponding Voronol polygons. (c) the set E~ of edges which delimit the two annull. Together, these a n n u l form the proto-column. (d) the set I~, which delimit the column. Note that plxels In the outer annulus will be mapped outward, whereas those In the Inner annulus will be mapped inward.

204

Computer Simulation of Cortical Polymaps

205

As in our definition of the protoannulus, above, we define the

f

set

E~

=

U ,~., -

p,p'

~:', e

where @" = {e • l~,lp, p' are adjacent contour points on the same contour}. By construction, the annulus is entirely contained in the protoannulus. Further, it sometimes represents the inner portion of the protoannulus. In this case, we call the mapping "outward." In the case where the annulus represents the external portion of the protoannulus, we term the mapping "inward" (see Figure 8).

8.4. Examples In light of the above discussion, consider Figure 12, which depicts a column with its two neighboring columns. The column under consideration is convex, which represents the simplest case of our algorithm. We show the Voronoi polygons which would correspond to choosing six sample points along the contour of the column, though in practice, the number of sample points used would be at least an order of magnitude greater. The medial and dual-medial contours are readily seen to corresepond to the edges of the Voronoi polygons. The medial contour bounds a slit which is infinitessimal in width (Ahlfors, 1966).

O

[]

[]

FIGURE 16. A portion of a contour with Its accompanying medial and dual-medial contours. To perform the warping itself, we sample the contour along its circumferential index, and use the normals at those points to construct a set of matched quadrilaterals. In the figure above ~abcd will be warped into babel using a blllnear transformation.

Figure 13 depicts the complications associated with introducing a nonconvex aspect to the column. A portion of the dual-medial contour is slit-like here as well. Figure 14 depicts a column with an embedded hole. This column has two contours, each of which gives rise to a medial and dual-medial contour. The outer contour has a medial contour which corresponds to the inner contour's dual-medial contour. In this case, the medial contour of the inner contour is reduced to an infinitessimal circular region. In Figure 15, the result of embedding a different colunn in the original column is considered. The results are very similar to those in the previous figure, with the exception that the medial contour of the hole boundary is no longer degenerate.

8.5. Decomposition into Quadrilaterals

[]

[]

2 FIGURE 15. (a) A column with another column embedded In It. The column consists of two concentric annull. (b) the corresponding Voronol polygons. (c) the set Ef of edges form four contours which delimit the two annull. (d) the set I~l of edges which delimit the column. This resembles the previous case, except that here the Innermost contour of El Is no longer degenerate. (e) Proto-annull which will represent the mapping from (f) annull.

The principal advantage of the above construction, where we have reduced arbitrary columns to annuli, is that annuli can be addressed by two quasi-polar coordinate indices: a "circumferential" or annular index, and a "radial" index, expressed as distance from the medial contour. A similar approach has been followed by (Wolberg, 1990) in his image warping algorithm. Given such an addressing scheme, it is natural to decompose the problem in a radial fashion. We sample the contour along its circumference (or use the previous sample points), and consider each radial slice in turn. At each sample point along the contour, we construct the normal, and calculate the intersection points of the normal with the medial and dual-medial contours (see Figure 16). Each two normals corresponding to adjacent contour points define two quadrilaterals, formed by connecting pairwise the intersection points with each of the contours (see Figure 16). These quadrilaterals define a polygonal approximation to the contours. To map one quadrilateral to another can be accomplished by a bilinear transformation (Wolberg, 1990). The coefficients of the bilinear transformation found by solving the linear equations determined by the coordinates of the intersection points. To warp the entire image, we warp each quadrilateral area in sequence, each using its own bilinear warping. As in other sampling cases, the above construction is sensitive to aliasing. In this case, if the sampling is too coarse to properly describe the curvature of the contour, the normals may intersect each other before intersecting the medial or dual-medial contour. This is simple to test for and avoid.

206 Algorithm (stated here for the one eye): P=O for each contour C~ E St find the contour points Pep of contour C, add these to P: P = P U P,~ construct the Voronoi diagram V for P

P. Landau and E. Schwartz for each o,p @ V construct the set of edges {ecpf}. for each contour Cc ~E SL using the edges {e,cp}find the protoannulus contour E~ using the v,p find the annulus contour El construct a matched set of quadrilaterals in the annulus and protoannulus using the resulting quadrilateral sets, texture map a full proto image to a full column image.