The Voronoi spatial model : two-and three-dimensional applications in image analysis

Traditional raster-based map or image manipulation techniques are known frequently to produce different results from vector implementations of the same operation. In addition, rasters require resampling and information loss for most resizing and rotation transformations, while vector maps in a GIS require complex and error-prone procedures to define neighbours, adjacency, connectedness, etc. The two views of space appear incompatible. Many GIS and remote sensing users have wished for an integrated spatial model having the best features of, and permitting identical operations on, both data types in the same data structure. In the GIS world, the integration of interpolation questions using any mixture of point, line and polygon data is very desirable, as is a data structure that is locally updateable, and hence dynamic in the sense that the result of object editing or moving is immediately available for subsequent analysis. In remote sensing, raster/vector integration is required, as are consistent techniques for handling spatial or temporal data inhomogeneity or anisotropy. There is also a need for the consistent handling of hierarchical spatial entities, from the pixel up to the classified polygon, together with error estimates of the consequences of this aggregation. A variety of other desired properties of the basic spatial model and consequent data structure could be enumerated. The Voronoi spatial model, as currently under study at Laval, appears to address many of these issues. Spatial adjacency is defined as directly as it is in a raster, and yet it is rotationally invariant like a vector map. Both previous spatial models may be readily incorporated in it. It is dynamically modifiable, it forms the basis for a variety of useful interpolation or resampling operations on point, line or polygon data, together with error estimates, and it can be used in a time-sequential fashion to express the evolution of the image or map. The same basic spatial operations are as applicable to terrain model construction as to polygon or raster maps, and the techniques are not inherently restricted to two dimensions. It therefore appears to have real potential as a general spatial method with applications in remote sensing, GIS and other spatial disciplines. At a conference several years ago, a friend (Geoff Dutton) and the first author were enthusiastically discussing the “intrinsic structure of data” (visualized as a cloud of points in three dimensions) without knowing exactly what we were talking about. What we did mean, however, was that artifacts of scale, coordinate system and algorithm did not affect the underlying spatial adjacency relationships-“topology” if you will-between neighbouring objects in space. We held these truths to be self-evident! The reality, of course, was (and is) completely otherwise. We worked (and still work) under four tyrannies (at least): the tyranny of being forced to think in terms of an often arbitrary global coordinate system (eg, raster); the tyranny of defining relationships on the basis of distance; the tyranny of defining relationships on the basis of the intersections of line segements (eg, vector GIS); and the tyranny 1 Dept des sciences géodésiques et de télédétection, Université Laval, Sainte-Foy, Quebec, Canada GlK 7P4 11 of being forced to rebuild our world in a batch environment whenever some local change is made. As we all know to our cost, coordinates do not transform to connected graphs without toil, blood, sweat and tears. Why is this so? The Newton-Leibnitz controversy has been mentioned elsewhere [5], and how the Newtonian view of space has dominated our thought until very recently. Initially, however, the quickest way to illustrate the point is to take an application that is of interest in both raster and vector worlds: digital elevation models (DEMs). These are usually expressed as either a network of triangles (a TIN) where relationships are expressed between objects and a triangle edge represents a pair of adjacent data points, or as an elevation grid (a DTM) where relationships are implied between square pieces of space. These two views-expressing the relationships between objec ts , or having a space-covering tiling where relations between adjacent tiles are known-are usually considered mutually exclusive, and we resign ourselves to a technologic schizophrenia where map objects in one world dissolve into fragments (pixels) in the other, and apparently equivalent operations become strangely distorted between worlds (eg, overlay, corridor and topology). Must we live also with the heartbreak of schizophrenia? Of course, one aspect that makes the whole problem even more infuriating is the impression that performing basic cartographic operations by hand is not particularly demanding (tedious yes; demanding no). Can we simulate, on the computer, the typical coordinate-free manual construction techniques? “The Problem”, therefore, is to determine whether the limitations resulting from the raster/vector dichotomy, together with those associated with excessive dependence on global coordinates, can be ameliorated by the use of an alternative spatial model. In order to address these questions, let us review the properties of the two tradtional “models of space”. This will be examined first with respect to the previously mentioned interpolation problem for DEMs. INTERPOLATION A PROTOTYPICAL QUESTION IN SPATIAL ANALYSIS An interesting sidelight on our research undertaken in the last few years is its origins in the interpolation problem for arbitrarily distributed data points (originally geologic data). This is approximately mid-way between the raster world of image analysis and the vector world of most GISs. Traditionally, Voronoi spatial model ITC Journal 1992-1 some kind of function must be described about a portion of the map, and elevations estimated at grid nodes for subsequent contouring. The more we looked at this problem, the less satisfied we were with almost any of the underlying assumptions, arbitrary techniques and approximations being used. To name merely a few: there is no need to estimate intermediate values onto a grid, in fact it is counterproductive; polynomial functions in Cartesian coordinates create massive problems in smooth surface continuity; the selection of a set of neighbouring points for local approximation or averaging can never be achieved successfully (ie, without surface discontinuities caused by arbitrary acceptance or rejection of points) using an exclusively metric criterion; and any weighting function used to average these neighbours will not be successful if exclusively metric in nature. There had to be a better, a more “natural” approach! A similar feeling occurred with the examination of GIS problems, where spatial neighbourhood/adjacency issues were combined with database concepts. It is clearly desirable that one object in the database be related to one object in the spatial referencing system, but, in a raster system, space is sliced into little squares which then have to be allocated to objects (with limited resolution). Alternatively, in vector systems, a great deal of overhead is expended on connecting little pieces of lines that (approximately) meet each other-with considerable difficulties in extracting the graph network (topology) from the coordinate information (geometry). The old familiar feeling recurred when working with finite-element or finite-difference models of groundwater flow and similar problems. Either you live with a very coarse grid-based representation of your aquifer or rock formation, or you create a system of equations and spatial elements that even hardened practitioners agree is not intuitively simple. It became apparent that, in interpolation as in other spatial questions, the “adjacency problem” (finding the neighbouring data points to the location where an elevation estimate was desired) was only poorly handled by a metric method (selecting all points within a specified distance, with or without various restricting conditions). Thus spatial adjacency entered as a key issue to be resolved. The set of spatially adjacent points had to be stable under various minor perturbations of sampling location or coordinate system, and data points had to enter or leave the neighbourhood set in a consistent manner. Considerable effort was expended in attempting to define such a system, and it rapidly became clear that the answer did not lie in a more clever computer data structure, but in a consistent concept of adjacency of individual points in space. A triangulation technique was implemented, but this did not in itself say which of the many possible triangulations of a set of data points was appropriate. A few years ago, a reasonable consensus was reached that the Delaunay triangulation was the only one with the property of stability whatever the data input order, and with only local perturbation when a new point was added (eg, [7, 1, 12]). This appeared to settle the problem. Then, however, a new issue arose: how to provide a suitable weighting for each of the data points selected as neighbours to the interrogation point. Further work led to the examination of the basis of the Delaunay triangulation-the underlying Voronoi region set that was its dual. This could readily be envisaged as the “blowing-up of bubbles” around each data point until they all met. Thus the graph-theoretic result (the linkage of adjacent points by a triangular network) was achieved by a metric method (the definition of “zones of influence” about each point such that all points falling within its zone were closer to that point than to any other). The interpolation problem could be reduced to the insertion of a desired interrogation or sampling point into the previouslyconstructed Voronoi diagram of the data points, and examining the changes. This work was described by Sibson [ 13] and Gold [5], and some of the interpolation problems to be resolved