SPATIAL DATA AND THE VORONOI TESSELLATION
This article contains the following executables: SPATIAL.ZIP COASTS.ZIP
Hrvoje Lukatela and John Russell
Hrvoje Lukatela holds a Masters Degree in Geodetic Engineering from the University of Zagreb and has practiced as a survey engineer, software engineer, and database designer for more than 20 years. He is the author of the Hipparchus Geopositioning Model. John Russell began work as a programmer for IBM in Toronto in 1957 and was IBM Canada's first senior systems engineer. More recently, John was vice president, technology for the 1988 Calgary Olympic Winter Games. They can be contacted at Geodyssey Ltd., 300 815 8th Ave. SW, Calgary, Alberta, Canada T2P 3P2, Fax: 403-266-7117.
In the early days of computing, the data we worked with consisted of integers, real numbers, and characters. Later, we moved on to time and money data. Today, as we increasingly deal with environmental and other geographic information, we need new ways of looking at spatial data.
For millennia, cartographers have attempted to represent the round Earth on flat maps. The first four decades of geographic information systems (GIS) have attempted to automate this process, typically using a "flat Earth" paradigm of map sheets and two-dimensional coordinates. The result has been an unwieldy collection of complex math, preset views, and location-dependent precision.
An alternative is to model the Earth using a "round Earth" paradigm. In this way, we can roam freely with our geographic applications, modeling surface features without restriction, and calculating spatial relationships with uniform high precision.
In this article we'll demonstrate an approach to representing the location, storage, retrieval, and manipulation of data in terms of its spatial relationships. We'll use elementary trigonometry and three-dimensional vector algebra to develop programs that demonstrate the key ideas. Then we'll build on these concepts to show how you can develop a complete GIS that has unprecedented speed and precision, without the use of a conventional GIS solution.
To illustrate these concepts, let's build a simple geographical atlas that lets you roam anywhere on the globe, viewing surface features at varying scales. In the general case, we would model our geographic features of interest as points, lines, areas, or volumes.
- Points might represent cities or survey monuments.
- Lines might represent roads or flight paths.
- Areas might represent islands or properties.
- Volumes might represent geological formations or controlled airspaces.
For simplicity, this application will deal only with line objects. The geographic location of a line object can be given by an ordered set of vertex coordinates. Figure 1 illustrates some sample application objects. Listing One (page 96) provides their numeric specification in the familiar terms of latitude and longitude--the angles that give the location of geographic features relative to the equator and a prime meridian. The frame of reference is geocentric, meaning that the angles are measured from the center of the Earth; see Figure 2. Latitude is labeled 
and longitude is labeled 
.
While early scientists thought of the planet as a perfect sphere, we now know it is somewhat flattened at the poles, an "ellipsoid of rotation." However, since the eccentricity of the Earth is not great (less than a third of one percent), we'll assume for the moment that the Earth is indeed a perfect sphere.
Since latitudes and longitudes are angles, when we work with them we must be prepared to calculate sines, cosines, tangents, arc tangents, and the like. Even with today's math coprocessors, this can get messy. For instance, have you ever tried to find the tangent of 90 degrees? You will if your application deals with objects in the polar regions. Generally, such calculations lack a geographically uniform distribution of precision. Luckily, a point's location on the Earth's surface can be represented in other ways.
Consider a 3-D geocentric space having three orthogonal axes projecting through the equator and the poles. Call these axes X, Y, and Z. Now we can locate a point on the surface with the three coordinates x,y,z; see Figure 3. The X axis projects through the Atlantic ocean just off West Africa, the Y axis projects through the Indian ocean just west of Sumatra, and the Z axis projects through the North Pole. The pictured surface point P(x,y,z) might be somewhere in northern Afghanistan.
Given the 3-D space just described, there's another way to describe the location of a surface point. Instead of referring its coordinates directly, we could describe the vector perpendicular to the surface at that point. For a perfectly spherical Earth, this normal would pass through the center. It has unit length, and its direction is defined by the angles formed between it and the X, Y, and Z axes. These angles are called direction angles.
We'll be working with the cosines of the direction angles--direction cosines--labeled di, dj, and dk, respectively. The point in Afghanistan can now be referred to as P(di,dj,dk); the point off West Africa has the coordinates (1, 0, 0); the point in the Indian ocean is at (0,1,0); the North Pole point is at (0,0,1); and the South Pole point is at (0,0,-1); see Figure 4.
Recording direction cosines as double types in C typically provides sub-millimetric precision. This usually surpasses the precision of your very best field data. Listing Two (page 96) shows some geometric vector-algebra functions and their supporting structures and constants.
Most developers are familiar with latitude and longitude. In addition, there are "flat-Earth" coordinate systems such as UTMs and State Plane that are used by surveyors and map makers. Few, however, are familiar with direction cosines. Consequently, if our new system is to be of any use, we'll need an efficient method of converting between direction cosines and these other coordinates. For simplicity, we'll restrict input to latitudes and longitudes. We begin by converting a file of geographic data to direction cosines.
At first glance, using direction cosines to locate a point on the Earth's surface seems inefficient, since we're trading two items (latitude and longitude) for three. But in modeling geographic objects, we often have multiple locations associated with specific objects. For example, a line object such as a coastline, river, or road is usually modeled as an ordered sequence of connected vertices. In such instances, we might select a single, "central" location and relate all the associated vertices to it. But will this "differential" position encoding be effective?
In developing planar projections, map makers look for the recognizability of shapes (conformity) and the uniformity of scale in all directions (isometry). One of their best efforts is the stereographic projection which, over moderate distances, produces a view of the Earth that's both conformal and isometric. (Despite its name, this projection of 3-D onto 2-D provides no depth perception.)
If an object is restricted in size, it can be represented in the plane of a stereographic projection without significant distortion. This means we can use a specific scheme of differential location recording in which each vertex of the line is encoded as a stereographic planar displacement from some central position. As such, this differential value will have just two components, say dx and dy.
Using only short int types for dx and dy, resolution of better than a meter can be maintained for surface objects as large as ten kilometers in extent. For better resolution, we can use float or long int types; for poorer, we can use signed char.
So, typically, we'll have traded in three doubles for two short ints, a significant reduction in storage requirements. We refer to these differentially encoded coordinates as local coordinates. The full, three-element direction cosine global coordinates can easily be reconstructed at any time, using only elementary vector algebra.
Since we're creating an application to select and display terrestrial "objects," it makes sense to store the data externally under some kind of object-oriented scheme. But how should the objects be indexed?
Conventional wisdom suggests that we index our data on the basis of decomposition (or hashing) of the object's coordinates. For this application, however, let's try something different.
First, let's establish a file as the general repository for the local coordinates of all the vertices of all of the objects modeled. We'll provide access to the individual parts of this file using file pointers.
Next, let's set up an index file of object headers. Each header will hold the object's identifier, the global coordinates of its center, a file pointer to the local coordinates of its vertices, and the vertices count. The object's identifier can serve as a link to its other attributes (if any). The header will also contain an estimate of the object's geographic extent, described shortly.
When we load the database with an object, we can determine its "center" by calculating the "vector mean" of the direction vertices' cosines. We can then use this center to differentially encode coordinates for the vertices.
Because this application will let you zoom in and out through a wide range of scales, we're providing two classes of line objects: those required for close-ups (dense) and those needed only for wide-area presentations (sparse). Since the application is to be interactive, we'll want to reduce unnecessary data retrieval and processing time (especially if we don't have floating-point hardware).
To determine if objects are "onscreen" or not, the application will need to know their geographic extents. Using vector algebra, we can calculate these as surface distances, for which we'll need arc (or great circle) distances. While we're loading objects into the database, it will prove useful to calculate and store the maximum great circle distance that any vertex is displaced from the object's center; see Figure 5.
Listing Three (page 98) and the called functions in Listing Two provide code to read and convert location data to direction cosines, differentially encode them, calculate distances, and build a location-dependent database.
Our application provides a window on the world, so to speak, by displaying objects that come within a field of view you select. A view is defined in terms of location and scale. The location of the display's center can be expressed as a latitude and longitude. Scale can be expressed as the ratio between a distance on the screen and a distance on the ground. Figure 6 illustrates such a window, while Listing Four (page 99) and the called functions in Listing Two show the code needed to establish an initial view and scale.
Now that field of view is defined, we can locate objects that might come into that field using distance calculations. If you think of the display as circular rather than rectangular, then you can calculate a maximum radius for the display. You can go to the database and select those objects that might be displayed. (The graphics-library clipping function will fine-tune the selection later.)
The header for each object contains the maximum distance of any vertex from the object's center. This was calculated and stored when we loaded the database. So, to determine if the object might be in the field of view, simply: 1. find the distance between its center point and that of the display; 2. subtract the maximum radius for the object; and, 3. subtract the maximum radius for the display.
If the result is negative, you'll want to retrieve the object from the database for further processing; otherwise, ignore it. Figure 7 illustrates both of these conditions. Listings Four and Two provide code to make the selection and bring the selected objects into memory.
Next we need to project each object's vertices into the plane of the display (projection plane), which is generally not the same as that of the display. For simplicity, we'll go back to the sphere and reproject the object's vertices using, for this example, the stereographic projection. Other projections--gnomonic, orthographic, Mercator, and the like -- might also be used. Gnomonic is the easiest, but stereographic looks better and is worth the effort. Listings Four and Two give the code to perform the projections and draw the objects. (For more about map projections, see "Map Projections Used By The U.S. Geological Survey, Sec. Ed.," Geological Survey Bulletin 1532, Department of the Interior, U.S. Government Printing Office, Washington, DC, 1984.)
Suppose you want to change the scale or view of the display. Simply modify these items and repeat the previous operations. A simple outside loop that changes the scale or map center point will work. Listing Four shows code to accept changes via the sign and arrow keys.
That completes our simple atlas application. Even with the slowest PC, you can now inspect the world's coastlines without preselection of view or scale. To more fully exercise the system, raw world-coastline data (in ASCII form) is available electronically, as is a prebuilt world-coastline database (in binary); an executable View program in DOS real mode, compiled for VGA with math coprocessor emulation; and ASCII source code for the programs; see "Availability" on page 5.
Since the Earth is closer to an ellipsoid of rotation than a sphere, we'll need to extend our vector algebra. The required quadratic vector algebra has been fully implemented in the Hipparchus Geopositioning Model with significant improvements in speed and precision over conventional geodesy methods. (See Geodesy, by Henry D. Bomford, Oxford University Press, 1973.)
For this sample application, we calculated a local center point for each object and then used this to select objects from the database. We also used these center points to encode the large number of vertex coordinates associated with our objects.
Suppose we could precalculate a set of center points that would serve the same purposes for all the objects in the database. Ideally, such a set of center points would provide both fast spatial indexing and a flexible association with objects. In such a spatial index, each indexed database "bucket" would hold some prescribed maximum number of object-defining coordinates. Then we could have geographically large cells for surface regions where we've little or no data and geographically small cells where we have a lot of data.
The Hipparchus Geopositioning Model implements just such a scheme using a flexible partitioning system called a "Voronoi cell structure." Figure 8 shows one such tessellation of the Earth. The structure illustrated would be suitable for indexing population-related data objects. Voronoi cell structures are always global, even if the application is localized. A cell structure is defined by its cell center points. For each cell, the structure includes a unique cell identifier, the global coordinates of the cell's center point, the cell's maximum radius, and an ordered list of neighbor-cell identifiers. The boundaries between cells exist only mathematically.
The special property of the Voronoi cell structure is that any surface point can be classified unambiguously as belonging to one cell or another on the basis of surface distance. A point is always closer to the center point of its "owner" cell than to the center point of any other cell. For a discussion of the Voronoi tessellation of the plane, see Algorithms, Second Edition, by Robert Sedgewick (Addison-Wesley, 1988). For a description of its adaptation to the surface of the ellipsoid, see "Hipparchus Geopositioning Model: An Overview," by H. Lukatela in Proceedings of Auto Carto 8 (American Society for Photogrammetry and Remote Sensing, 1987), and the Hipparchus Tutorial, by Ron V. Gilmore (Geodyssey, 1992).
In the context of a Voronoi cell structure, an object's vertices are associated with their closest-cell center points as well as an object header. Objects can then be defined without geographic size restriction of any kind. Objects can consist of sets of points, lines, or regions spanning any number of cells. Regions can be nonsimply connected: An island group can be modeled as a single object, islands can have interior lakes with islands, and so on. Volumes can be modeled as regions having elevation or depth attributes. Figure 9 shows the intersection of two overlapping region objects in the Voronoi context.
Cell center points rather than object centers are used for the differential encoding of coordinates. Lists are maintained for each case:
- For point set objects, lists of cells occupied by their points.
- For line set objects, lists of cells traversed by their line segments.
- For region objects, lists of their interior cells as well as the cells traversed by their boundary rings.
For more about these data structures, refer again to "Hipparchus Geopositioning Model: An Overview" and the Hipparchus Tutorial.
When used as an index to objects stored externally, the Voronoi cell structure proves remarkably effective in reducing unnecessary disk accesses. Not only are all the cells containing object data known to the application program, but cells associated with open windows are known as well. As you pan and zoom the window, precise retrieval instructions can be fed to the database.
References to random locations are traced to their owner cells by a geographically direct search route. Ownership of a point by a particular cell is confirmed when a comparison of distances with the cell's immediate neighbor cell center points shows them to be more distant.
In this application, we had to search our entire index to determine which data was to be selected. This was because we knew of no way to map directly from the 3-D ordered domain of our real-world objects into the linearly ordered domain of the computer. But when we associate these objects with a Voronoi cell structure, the situation changes.
The unambiguous classification of object vertices into a specific, linearly ordered structure of cells makes possible the use of hierarchical searches for the data, resulting in significant efficiencies.
Since the order of cell identifiers in a cell structure is irrelevant to its algorithmic operation, cells can be arranged in any order. Therefore, data-access bias can be arbitrarily imposed without affecting the logic of the application.
The demand for efficient handling of crushing volumes of spatial data has arrived. Round-Earth vector algebra and the Voronoi tessellation can be combined to provide unrestricted modeling and efficient manipulation of terrestrial objects. Precise spatial indexing can be provided on the basis of distance calculations rather than coordinate decomposition. Monolithic geographic information systems may soon be history.
_SPATIAL DATA AND THE VORONOI TESSELLATION_
by Hrvoje Lukatela and John Russell
[LISTING ONE]
* Townsite
50.42 -100.13
50.41 -100.15
50.40 -100.16
50.39 -100.18
50.38 -100.20
50.37 -100.20
50.37 -100.28
50.20 -100.28
50.20 -100.16
50.20 -100.00
50.42 -100.00
50.42 -100.13
* Highway
50.45 -99.39
50.31 -100.04
50.17 -100.20
49.56 -100.48
49.42 -101.26
* Highway
50.31 -100.04
50.31 -101.28
* Flight Path
50.00 -101.28
49.85 -99.40
* River
49.37 -99.81
49.38 -99.83
49.42 -99.86
49.43 -99.88
49.44 -99.90
49.46 -99.92
49.47 -99.94
49.47 -99.96
49.48 -99.98
49.49 -100.00
49.49 -100.02
49.49 -100.04
49.48 -100.06
49.48 -100.08
49.47 -100.11
49.47 -100.13
49.47 -100.16
49.46 -100.18
49.46 -100.20
49.45 -100.22
49.44 -100.24
49.43 -100.24
49.42 -100.21
49.41 -100.19
49.41 -100.18
49.41 -100.16
49.41 -100.14
49.41 -100.12
49.41 -100.10
49.40 -100.08
49.38 -100.08
49.37 -100.10
49.36 -100.12
49.35 -100.15
49.35 -100.17
49.36 -100.19
49.37 -100.21
49.38 -100.23
49.38 -100.26
49.39 -100.28
49.39 -100.30
49.39 -100.33
49.38 -100.35
49.38 -100.38
49.38 -100.41
49.37 -100.43
49.36 -100.45
49.35 -100.47
49.34 -100.49
49.33 -100.51
49.32 -100.52
49.31 -100.53
49.29 -100.53
*
[LISTING TWO]
/* ---------------------------------------------------------------- *
* ALGEBRA functions: A sampling of geometronical vector algebra *
* functions and their supporting manifest constants and structure *
* declarations. *
* The following code is derived from similar functions which are *
* a small part of the Hipparchus Library. For simplicity, it lacks *
* the "fuzz control" and other programming elements of practical *
* numerical significance. *
* Programmer: Hrvoje Lukatela, September 1992. *
* Geodyssey Limited, Calgary - (403) 234 9848, fax: (403) 266 7117 *
------------------------------------------------------------------ */
#include <math.h>
#define PI 3.14159265358979324
#define DEG2RAD (PI / 180.0) /* degrees to radians... */
#define RAD2DEG (180.0 / PI) /* ... and vice versa */
#define LC_SCALE 32000.0 /* local coordinate scale factor */
struct plpt { /* point in a Cartesian projection plane */
double est;
double nrt;
};
struct lclpt { /* local (object) coordinates */
short est;
short nrt;
};
struct dpxl { /* display screen coordinates */
short x;
short y;
};
struct ltln { /* point latitude-longitude, radians */
double lat;
double lng;
};
struct vct3 { /* 3-D vector; x,y,z direction cosines */
double di;
double dj;
double dk;
};
struct vct2 { /* as above, in plane, internal use */
double di;
double dj;
};
struct indexRec { /* line segment database index record */
struct vct3 center; /* nominal object center point */
double radius; /* in radian arc measure */
long fileOffset; /* offset in the coordinate data file */
short vertexCount; /* count of coordinate vertices */
short segmentId; /* for possible application use? */
};
/* ----- Transform Latitude and Longitude Angles to Direction Cosines ----- */
void LatLongToDcos3(const struct ltln *pa, struct vct3 *pe) {
double cosphi;
cosphi = cos(pa->lat);
pe->di = cosphi * cos(pa->lng);
pe->dj = cosphi * sin(pa->lng);
pe->dk = sin(pa->lat);
return;
}
/* ---- Transform Direction Cosines to Latitude and Longitude Angles ----- */
void Dcos3ToLatLong(const struct vct3 *pe, struct ltln *pa) {
pa->lat = atan2(pe->dk, sqrt(pe->di * pe->di + pe->dj * pe->dj));
pa->lng = atan2(pe->dj, pe->di);
return;
}
/* ---- Normalize a 3-D Direction Cosine Vector ---- */
void NormalizeDcos3(struct vct3 *vcc) {
double d;
d = 1.0 / sqrt(vcc->di * vcc->di +
vcc->dj * vcc->dj +
vcc->dk * vcc->dk);
vcc->di *= d;
vcc->dj *= d;
vcc->dk *= d;
return;
}
/* ----- Normalize a 2-D Direction Cosine Vector ---- */
void NormalizeDcos2(struct vct2 *vcc) {
double d;
d = 1.0 / sqrt(vcc->di * vcc->di + vcc->dj * vcc->dj);
vcc->di *= d;
vcc->dj *= d;
return;
}
/* ----- Spherical Arc (Great Circle Distance) - First Approximation ----- */
double ArcDist(const struct vct3 *pea, const struct vct3 *peb) {
double chord, sqChord;
sqChord = (peb->di - pea->di) * (peb->di - pea->di) +
(peb->dj - pea->dj) * (peb->dj - pea->dj) +
(peb->dk - pea->dk) * (peb->dk - pea->dk);
chord = sqrt(sqChord);
return(chord + ((sqChord * chord) / 24));
}
/* ----- Direct Stereographic Projection (Map, Sphere to Plane) ----- */
void MapStereo(const struct vct3 *p0,
const struct vct3 *pe, struct plpt *pw) {
struct vct3 prln;
double t, am, bm, ap, bp, cp, xi, yi, zi;
/* ---------------------------------------------------------------- */
/* Find tangency point relative values. */
cp = sqrt(p0->di * p0->di + p0->dj * p0->dj);
am = -(p0->dj / cp);
bm = p0->di / cp;
ap = -(p0->dk * bm);
bp = p0->dk * am;
/* Intersection of the projection line and the intersection plane. */
prln.di = -(p0->di + pe->di);
prln.dj = -(p0->dj + pe->dj);
prln.dk = -(p0->dk + pe->dk);
NormalizeDcos3(&prln);
t = -((p0->di * pe->di + p0->dj * pe->dj + p0->dk * pe->dk - 1.0) /
(p0->di * prln.di + p0->dj * prln.dj + p0->dk * prln.dk));
xi = pe->di + prln.di * t;
yi = pe->dj + prln.dj * t;
zi = pe->dk + prln.dk * t;
/* Stereographic plane coordinates are the oriented distances from
the intersection point to the meridian and prime vertical plane. */
pw->est = am * xi + bm * yi;
pw->nrt = ap * xi + bp * yi + cp * zi;
return;
}
/* ----- Inverse Stereographic Projection (Un-Map, Plane to Sphere) ---- */
void UnMapStereo(const struct vct3 *p0,
const struct plpt *pw, struct vct3 *pe) {
struct vct3 prln;
double gcx, am, bm, ap, bp, cp, cpsq;
double xe, ye, ze, xc, yc, zc, lymx, lxmy, root, t;
/* ---------------------------------------------------------------- */
/* Find the sphere/plane tangency point values: ap, bp, cp are
components of the "North" vector, and am, bm of "East" vector
in this point. The "East" vector has no "Z" axis component. */
gcx = sqrt(pw->est * pw->est + pw->nrt * pw->nrt);
cpsq = p0->di * p0->di + p0->dj * p0->dj;
cp = sqrt(cpsq);
am = -(p0->dj / cp);
bm = p0->di / cp;
ap = -(p0->dk * bm);
bp = p0->dk * am;
/* Find Cartesian coordinates of the projection center (xc,yc,zc)
and the projected point in the plane of projection (xe,ye,ze). */
xc = -p0->di;
yc = -p0->dj;
zc = -p0->dk;
xe = -xc + ap * pw->nrt + am * pw->est;
ye = -yc + bp * pw->nrt + bm * pw->est;
ze = -zc + cp * pw->nrt;
/* Find the intersection of ptc-pte line and the sphere.
Solution requires solving a quadratic in t, the line parameter. */
prln.di = -gcx;
prln.dj = -2.0;
NormalizeDcos2((struct vct2 *)&prln);
lymx = prln.dj * gcx - prln.di;
lxmy = -(prln.di * gcx + prln.dj);
root = sqrt(1.0 - (lymx * lymx));
t = lxmy - root; /* Find the closer of the two quadratic roots. */
/* Substitute the parameter in the parametric line equations. */
prln.di = xc - xe;
prln.dj = yc - ye;
prln.dk = zc - ze;
NormalizeDcos3(&prln);
pe->di = xe + t * prln.di;
pe->dj = ye + t * prln.dj;
pe->dk = ze + t * prln.dk;
NormalizeDcos3(pe);
return;
}
/* ----- Initialize Projection Plane / Display Translation and Scaling ----- */
void SetPlaneDisplay(double *xfmArray,
const struct plpt *w1, const struct plpt *w2,
const struct dpxl *d1, const struct dpxl *d2) {
double dx, dy, du, dv;
dx = (w2->est) - (w1->est);
dy = (w2->nrt) - (w1->nrt);
du = (double)((d2->x) - (d1->x));
dv = (double)((d2->y) - (d1->y));
xfmArray[0] = dx / du;
xfmArray[1] = dy / dv;
xfmArray[2] = du / dx;
xfmArray[3] = dv / dy;
xfmArray[4] = w1->est - xfmArray[0] * ((double)d1->x + 0.5);
xfmArray[5] = w1->nrt - xfmArray[1] * ((double)d1->y + 0.5);
xfmArray[6] = ((double)d1->x + 0.5) - xfmArray[2] * w1->est;
xfmArray[7] = ((double)d1->y + 0.5) - xfmArray[3] * w1->nrt;
return;
}
/* ----- Translate/Scale Point from a Projection Plane to the Display ----- */
void PlaneToDisplay(const double *xfmArray,
const struct plpt *w, struct dpxl *d) {
d->x = (short)(xfmArray[6] + xfmArray[2] * w->est);
d->y = (short)(xfmArray[7] + xfmArray[3] * w->nrt);
return;
}
/* ----- Translate/Scale Point from the Display to a Projection Plane ----- */
void DisplayToPlane(const double *xfmArray,
const struct dpxl *d, struct plpt *w) {
w->est = xfmArray[4] + xfmArray[0] * ((double)d->x + 0.5);
w->nrt = xfmArray[5] + xfmArray[1] * ((double)d->y + 0.5);
return;
}
[LISTING THREE]
/* ---------------------------------------------------------------- *
* BUILD Program: Construct a seamless global database of line *
* segments from text files defining the location of line vertices. *
* Two levels of detail are provided: dense and sparse. The input *
* consists of two flat ASCII text files; one each for dense and *
* sparse coastline segment vertices. The segment is begun with a *
* line containing an asterisk marker (*) as its first character. *
* The remainder of such line is ignored, and may be used to name *
* the segment. End-of-file is signalled by a similar marker line. *
* Line segment vertex coordinates follow each marker line, one *
* latitude/longitude pair per line. Coordinates are given in *
* degrees, as free-format, white-space or comma delimited decimal *
* fraction character strings. Numbers are signed according to *
* the international geographic coordinates sign convention: *
* westerly longitudes and southerly latitudes are negative. *
* ---------------------------------------------------------------- */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "algebra.c"
#define LINE_LNGTH 128
#define MAX_VRTX 1024
void buildFiles(FILE *, FILE *, FILE *);
void main(void) {
FILE *fpTextLines;
FILE *fpIndex;
FILE *fpCoords;
char *fnTextLines0 = "coast0.lns"; /* dense input segments */
char *fnTextLines1 = "coast1.lns"; /* sparse input segments */
char *fnIndex0 = "coast0.idx"; /* dense database index */
char *fnIndex1 = "coast1.idx"; /* sparse database index */
char *fnCoords = "coast.dat"; /* composite coordinate file */
/* ---------------------------------------------------------------- */
if ((fpTextLines = fopen(fnTextLines1, "rt")) == NULL) {
fprintf(stderr,"Input file (%s) open failed.\n", fnTextLines1);
exit(1);
}
if ((fpIndex = fopen(fnIndex1, "wb")) == NULL) {
fprintf(stderr,"Index file (%s) open failed.\n", fnIndex1);
exit(1);
}
if ((fpCoords = fopen(fnCoords, "wb")) == NULL) {
fprintf(stderr,"Data file (%s) open failed.\n", fnCoords);
exit(1);
}
fprintf(stderr, "\nSparse line segment input file...\n");
buildFiles(fpTextLines, fpIndex, fpCoords);
fclose(fpTextLines);
fclose(fpIndex); /* Note that the coordinate file stays open! */
if ((fpTextLines = fopen(fnTextLines0, "rt")) == NULL) {
fprintf(stderr,"Input file (%s) open failed.\n", fnTextLines0);
exit(1);
}
if ((fpIndex = fopen(fnIndex0, "wb")) == NULL) {
fprintf(stderr,"Index file (%s) open failed.\n", fnIndex0);
exit(1);
}
fprintf(stderr, "\nDense line segment input file...\n");
buildFiles(fpTextLines, fpIndex, fpCoords);
fclose(fpTextLines);
fclose(fpIndex);
fclose(fpCoords);
fprintf(stderr, "\nTwo-level coastline database created:\n");
fprintf(stderr, "Index files: %s, %s\n", fnIndex0, fnIndex1);
fprintf(stderr, "Coordinate data file: %s\n", fnCoords);
}
/* ---- Read Line Segment File, Write Index File and Vertex Coordinates ---- */
void buildFiles(FILE *fpTextLines, FILE *fpIndex, FILE *fpCoords) {
static struct vct3 vertex[MAX_VRTX];
struct indexRec indexRec;
struct plpt stereoPlaneVertex;
struct ltln inVertex;
struct lclpt shortVertex;
double s, d;
int i, fileLine, lineCount;
long totalVertexCount;
char inLine[LINE_LNGTH + 1];
/* ---------------------------------------------------------------- */
fileLine = lineCount = indexRec.vertexCount = 0;
indexRec.segmentId = 0;
totalVertexCount = 0L;
while (fgets(inLine, LINE_LNGTH, fpTextLines)) {
fileLine++;
if (inLine[0] == '*') { /* line segment header, end-of-file? */
if (indexRec.vertexCount) { /* process accumulated segment */
fprintf(stderr,"line:%d vertices:%d \r",
lineCount, indexRec.vertexCount);
indexRec.center.di = 0.0; /* find object "center" */
indexRec.center.dj = 0.0;
indexRec.center.dk = 0.0;
for (i = 0; i < indexRec.vertexCount; i++) {
indexRec.center.di += vertex[i].di;
indexRec.center.dj += vertex[i].dj;
indexRec.center.dk += vertex[i].dk;
}
NormalizeDcos3(&indexRec.center);
indexRec.radius = 0.0; /* center-far-vertex distance */
for (i = 0; i < indexRec.vertexCount; i++) {
d = ArcDist(&indexRec.center, vertex + i);
if (d > indexRec.radius) indexRec.radius = d;
}
if (indexRec.radius < 1.0e-10) {
indexRec.radius = 0.0;
s = 0.0;
}
else s = LC_SCALE / indexRec.radius;
indexRec.fileOffset = ftell(fpCoords);
fwrite(&indexRec, sizeof(struct indexRec), 1, fpIndex);
for (i = 0; i < indexRec.vertexCount; i++) {
MapStereo(&indexRec.center, vertex + i,
&stereoPlaneVertex);
shortVertex.est = (short)(s * stereoPlaneVertex.est);
shortVertex.nrt = (short)(s * stereoPlaneVertex.nrt);
fwrite(&shortVertex, sizeof(struct lclpt),1, fpCoords);
}
totalVertexCount += indexRec.vertexCount;
indexRec.vertexCount = 0;
lineCount++;
}
}
else { /* next in a series of line segment vertices */
inVertex.lat = DEG2RAD * atof(strtok(inLine, " ,\t\n"));
inVertex.lng = DEG2RAD * atof(strtok(NULL, " ,\t\n"));
if (((fabs(inVertex.lat) < 1.0e-10)
&& (fabs(inVertex.lng) < 1.0e-10)) /* 0.0 lat, 0.0 long? */
|| (fabs(inVertex.lat) > 0.5 * PI)
|| (fabs(inVertex.lng) > PI)) { /* lat/long out of range? */
fprintf(stderr,"\nBad data, file line %d", fileLine);
exit(1);
}
if (indexRec.vertexCount == MAX_VRTX) {
fprintf(stderr,"\nSegment vertex buffer overflow...");
exit(1);
}
LatLongToDcos3(&inVertex, vertex + indexRec.vertexCount++);
}
}
fprintf(stderr,"...processed, line segments:%d, vertices:%ld\n",
lineCount, totalVertexCount);
}
[LISTING FOUR]
/* --------------------------------------------------------------- *
* VIEW Program: View the line segments from a global database. *
* Display areas of the Earth in stereographic projection at dif- *
* ferent scales and with different map center (tangency) points. *
* The user interface is simple: it provides for change of scale *
* and shift of map center point using the sign and arrow keys. *
* For each map scene, calculate its radius of view. Then retrieve *
* from the coordinate file only those coastline segments that *
* might come into view. Draw the line segments, relaying on the *
* display-graphics to clip the parts still outside the window. *
* --------------------------------------------------------------- */
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <string.h>
#include <graph.h> /* Using MS C 6.0 graphics... */
#include "algebra.c"
#define DISPL_WIDE 0.24 /* screen width, meters */
#define DISPL_HIGH 0.18 /* screen height, meters */
#define EARTH_RAD 6.371e6 /* approximate Earth's radius, meters */
#define MERIDIAN 1 /* grid drawing selectors */
#define PARALLEL 0
#define COLOR_FRAME 7
#define COLOR_GRID 1
#define COLOR_COAST 3
#define COLOR_SCALE 7
#define COLOR_PROMPT_TEXT 5
#define COLOR_PROMPT_KEYS 7
void drawDataLines(FILE *, FILE *);
void drawGrid(int, int);
void drawGridSegment(const struct ltln *, double, double, int);
static double xfmArray[8]; /* Plane/display transfomation */
static double maxDispDist; /* Radius of view */
static struct vct3 displayCenter; /* Map center, spherical */
void main(void) {
/* Initial map scale and center (projection plane tangency) point: */
double worldScale = 100.0e6;
struct ltln llStart = {DEG2RAD * 50.0, DEG2RAD * -100.0};
struct plpt pUpperLeft, pLowerRight, pNewCntr;
struct dpxl dUpperLeft, dUpperMid, dUpperRight,
dLowerLeft, dLowerMid, dLowerRight,
dLeftMid, dRightMid, dNewCntr, dCntr;
double worldWide, worldHigh;
struct videoconfig vcnfg;
struct vct3 sphVx;
int ich;
char outStr[32];
char *fnIndex0 = "coast0.idx";
char *fnIndex1 = "coast1.idx";
char *fnCoordinates = "coast.dat";
FILE *fpIndex0;
FILE *fpIndex1;
FILE *fpCoordinates;
/* ---------------------------------------------------------------- */
LatLongToDcos3(&llStart, &displayCenter); /* initial map center */
if ((fpIndex0 = fopen(fnIndex0, "rb")) == NULL) {
fprintf(stderr,"Index file (%s) open failed.\n", fnIndex0);
exit(1);
}
if ((fpIndex1 = fopen(fnIndex1, "rb")) == NULL) {
fprintf(stderr,"Index file (%s) open failed.\n", fnIndex1);
exit(1);
}
if ((fpCoordinates = fopen(fnCoordinates, "rb")) == NULL) {
fprintf(stderr,"Data file (%s) open failed.\n", fnCoordinates);
exit(1);
}
if (_setvideomode(_VRES16COLOR) == 0) { /* assume VGA graphics */
fprintf(stderr, "Graphics mode set failed.\n");
exit(1);
}
_getvideoconfig(&vcnfg);
dUpperLeft.x = dUpperLeft.y = 0;
dLowerRight.x = vcnfg.numxpixels - 1;
dLowerRight.y = vcnfg.numypixels - 1 - 20;
_setcliprgn(dUpperLeft.x, dUpperLeft.y,
dLowerRight.x, dLowerRight.y);
dCntr.x = (dUpperLeft.x + dLowerRight.x)/2;
dCntr.y = (dUpperLeft.y + dLowerRight.y)/2;
dLowerLeft.x = dLeftMid.x = dUpperLeft.x;
dUpperRight.x = dRightMid.x = dLowerRight.x;
dUpperMid.x = dLowerMid.x = dCntr.x;
dUpperRight.y = dUpperMid.y = dUpperLeft.y;
dLowerLeft.y = dLowerMid.y = dLowerRight.y;
dLeftMid.y = dRightMid.y = dCntr.y;
for (;;) {
worldWide = (worldScale * DISPL_WIDE) / EARTH_RAD;
worldHigh = (worldScale * DISPL_HIGH) / EARTH_RAD;
pUpperLeft.est = -worldWide * 0.5;
pUpperLeft.nrt = worldHigh * 0.5;
pLowerRight.est = worldWide * 0.5;
pLowerRight.nrt = -worldHigh * 0.5;
SetPlaneDisplay(xfmArray,
&pUpperLeft, &pLowerRight, &dUpperLeft, &dLowerRight);
UnMapStereo(&displayCenter, &pUpperLeft, &sphVx);
maxDispDist = ArcDist(&displayCenter, &sphVx);
_clearscreen(_GCLEARSCREEN);
_setcolor(COLOR_GRID);
if (worldScale > 10.0e6) drawGrid( 1, 1);
else if (worldScale > 3.0e6) drawGrid( 2, 3);
else drawGrid(10, 15);
_settextcolor(COLOR_PROMPT_TEXT);
_settextposition(vcnfg.numtextrows, 20);
_outtext("Press space bar to interrupt this scene...");
_setcolor(COLOR_COAST);
if (worldScale > 20.0e6) drawDataLines(fpIndex1, fpCoordinates);
else drawDataLines(fpIndex0, fpCoordinates);
sprintf(outStr, worldScale > 20.0e6 ?
"1 : %.0lfM" : "1 : %.1lfM", worldScale / 1.0e6);
_settextposition(vcnfg.numtextrows - 2, 37);
_settextcolor(COLOR_SCALE);
_outtext(outStr);
_setcolor(COLOR_FRAME);
_rectangle(_GBORDER, dUpperLeft.x, dUpperLeft.y,
dLowerRight.x, dLowerRight.y);
_settextcolor(COLOR_PROMPT_TEXT);
_settextposition(vcnfg.numtextrows, 1);
_outtext(" Press: (+)|(-) to change scale, (\x1b)|(\x18)|"
"(\x19)|(\x1a) to move center, (Esc) to quit.");
_settextcolor(COLOR_PROMPT_KEYS); /* highlight key characters */
_settextposition(vcnfg.numtextrows, 10); _outtext("+");
_settextposition(vcnfg.numtextrows, 14); _outtext("-");
_settextposition(vcnfg.numtextrows, 35); _outtext("\x1b");
_settextposition(vcnfg.numtextrows, 39); _outtext("\x18");
_settextposition(vcnfg.numtextrows, 43); _outtext("\x19");
_settextposition(vcnfg.numtextrows, 47); _outtext("\x1a");
_settextposition(vcnfg.numtextrows, 67); _outtext("Esc");
do {
if (ich = getch()) { /* non-0 scan code, ACSII character */
switch (ich) {
case 45: worldScale *= 2.0; break; /* - */
case 43: worldScale /= 2.0; break; /* + */
case 27: _setvideomode(_DEFAULTMODE); exit(0);/* Esc */
default: putch('\a'); ich = 0; /* invalid key */
}
if (ich) { /* OK, scale changed, enforce limits */
if (worldScale < 1.5625e6) worldScale = 1.5625e6;
if (worldScale > 200.0e6) worldScale = 200.0e6;
}
}
else { /* arrow or diagonal key, move map tangency point */
switch (ich = getch()) { /* get "extended scan code" */
case 73: dNewCntr = dUpperRight; break; /* up/right */
case 72: dNewCntr = dUpperMid; break; /* up */
case 71: dNewCntr = dUpperLeft; break; /* up/left */
case 75: dNewCntr = dLeftMid; break; /* left */
case 79: dNewCntr = dLowerLeft; break; /* down/left */
case 76: dNewCntr = dCntr; break; /* 5, center */
case 80: dNewCntr = dLowerMid; break; /* down */
case 81: dNewCntr = dLowerRight; break;/* down/right */
case 77: dNewCntr = dRightMid; break; /* right */
default: putch('\a'); ich = 0; /* invalid key */
}
if (ich) { /* OK, center was re-positioned */
DisplayToPlane(xfmArray, &dNewCntr, &pNewCntr);
UnMapStereo(&displayCenter, &pNewCntr, &displayCenter);
}
}
} while (ich == 0); /* i.e. until valid input was obtained */
}
}
/* ---- Traverse the Line Index and Display Close Line Segments ---- */
void drawDataLines(FILE *fpIndex, FILE *fpCoordinates) {
struct indexRec indexRec;
struct lclpt shortVx;
struct vct3 sphereVx;
struct plpt stereoPlaneVx;
struct dpxl displayVx;
double s;
int i;
/* ---------------------------------------------------------------- */
rewind(fpIndex); /* whole index will be searched sequentially */
while (fread(&indexRec, sizeof(struct indexRec), 1, fpIndex)) {
if (kbhit()) if (getch() == 32) break; /* space bar, break */
/* Skip line segments which are so far away from the display
that they can't have any vertices in it... */
if (ArcDist(&displayCenter, &indexRec.center)
> maxDispDist + indexRec.radius) continue;
s = indexRec.radius / LC_SCALE;
fseek(fpCoordinates, indexRec.fileOffset, SEEK_SET);
for (i = 0; i < indexRec.vertexCount; i++) {
fread(&shortVx, sizeof(struct lclpt), 1, fpCoordinates);
stereoPlaneVx.est = s * (double)shortVx.est;
stereoPlaneVx.nrt = s * (double)shortVx.nrt;
UnMapStereo(&indexRec.center, &stereoPlaneVx, &sphereVx);
MapStereo(&displayCenter, &sphereVx, &stereoPlaneVx);
PlaneToDisplay(xfmArray, &stereoPlaneVx, &displayVx);
if (i == 0) _moveto(displayVx.x, displayVx.y);
else _lineto(displayVx.x, displayVx.y);
}
}
}
/* ---- Display Latitude/Longitude "Rectangles" in the Display ---- */
#define RCT_LAT_DEG 10 /* Rectangle extent in latitude... */
#define RCT_LNG_DEG 15 /* ... and longitude, in degrees */
#define RCT_HALFDIAG_DEG 9 /* Rect. center-to-corner distance */
#define RCT_LAT_RAD (DEG2RAD * RCT_LAT_DEG) /* same in radians */
#define RCT_LNG_RAD (DEG2RAD * RCT_LNG_DEG)
void drawGrid(int densityLat, int densityLng) {
struct ltln llVx;
struct vct3 sphereVx;
double sLat, sLng;
int i, j, k;
/* ---------------------------------------------------------------- */
for (i = -80; i <= 80; i += RCT_LAT_DEG) {
for (j = -180; j < 180; j += RCT_LNG_DEG) {
llVx.lat = i * DEG2RAD + 0.5 * RCT_LAT_RAD;
llVx.lng = j * DEG2RAD + 0.5 * RCT_LNG_RAD;
LatLongToDcos3(&llVx, &sphereVx); /* grid rectangle center */
if (ArcDist(&displayCenter, &sphereVx)
> maxDispDist + DEG2RAD * RCT_HALFDIAG_DEG) continue;
sLat = (RCT_LAT_RAD / densityLat);
sLng = (RCT_LNG_RAD / densityLng);
for (k = 0; k < densityLat; k++) {
llVx.lat = i * DEG2RAD + k * sLat;
llVx.lng = j * DEG2RAD;
drawGridSegment(&llVx, sLng / 4, RCT_LNG_RAD, PARALLEL);
}
if (i == 80) continue;
for (k = 0; k < densityLng; k++) {
llVx.lat = i * DEG2RAD;
llVx.lng = j * DEG2RAD + k * sLng;
drawGridSegment(&llVx, sLat / 4, RCT_LAT_RAD, MERIDIAN);
}
}
}
}
/* ---- Display a Segment of a Meridian or a Parallel in Short Steps ---- */
void drawGridSegment(const struct ltln *llVxStart,
double step, double maxDist, int m) {
struct ltln llVx;
struct vct3 sphereVx;
struct plpt stereoPlaneVx;
struct dpxl displayVx;
double d = 0.0;
/* ---------------------------------------------------------------- */
LatLongToDcos3(llVxStart, &sphereVx);
MapStereo(&displayCenter, &sphereVx, &stereoPlaneVx);
PlaneToDisplay(xfmArray, &stereoPlaneVx, &displayVx);
_moveto(displayVx.x, displayVx.y);
do {
d += step;
if (d + 1.0e-10 > maxDist) d = maxDist;
llVx.lat = llVxStart->lat + (m * d); /* <m> is 0 or 1 */
llVx.lng = llVxStart->lng + ((1 - m) * d);
LatLongToDcos3(&llVx, &sphereVx);
MapStereo(&displayCenter, &sphereVx, &stereoPlaneVx);
PlaneToDisplay(xfmArray, &stereoPlaneVx, &displayVx);
_lineto(displayVx.x, displayVx.y);
} while (d != maxDist);
}
Copyright © 1992, Dr. Dobb's Journal