|
Scientific data is often in gridded data format, and contouring is a way to view it.
June 01, 1992
URL:http://www.drdobbs.com/database/contouring-data-fields/184408780
Copyright © 1992, Dr. Dobb's Journal
Copyright © 1992, Dr. Dobb's Journal
Copyright © 1992, Dr. Dobb's Journal
Copyright © 1992, Dr. Dobb's Journal
Bear, a programmer for National Systems & Research, can be reached at 325 Broadway, MC 425, R/E/FS4, Boulder, CO 80303 or [email protected].
A great deal of scientific and engineering data is available in a gridded (or regularly sampled) data format. For example, the Nested Grid Model (NGM) used by the National Weather Service consists of a series of atmospheric measurements at every 1.25 degrees of latitude and 2.5 degrees of longitude across the continental United States. This data is generated by supercomputers in Suitland, Maryland every 12 hours.
The most common way to view such data is through contouring. An example of contoured data familiar to many people is the topological maps used by hikers and skiers. Such maps are much easier to use than books specifying the elevation of thousands of different locations.
Contour lines of gridded data have several properties: They never cross, never split or join other contour lines, and never stop, except at the edge of the data field.
While it is easy to imagine situations where the physical contour lines violate one or more of these rules, such situations never arise in sampled data. The closest possibility is the unlikely event of adjoining data points having precisely the same value as the contour level. This could split the contour line, but because this causes a division-by-zero error, in these cases the contour line is simply stopped. In practice this rarely occurs.
We are using a linear interpolation between data points, so no more than one contour line ever crosses any grid cell-face for any contour level. (This is not guaranteed if more powerful interpolation techniques are used.)
Allowing for rotations, it is easy to show that there are only four possible orientations of a contour line within a single grid cell; see Figure 1. The final possibility has a symmetrical alternative, as in Figure 2. Either of these last two possibilities may be used as long as a consistent choice is made for each contour level.
As mentioned, a simple linear interpolation provides the precise entry and exit points of the grid cell. A straight line is then drawn between these points to form the contour. While it may be tempting to apply curve-smoothing functions (such as cubic splines) to the contour lines, this should be avoided. If curve smoothing is done, there is an excellent chance of erroneous results (such as crossed contour lines).
The contouring algorithm is quite simple: Ignoring boundary conditions for the moment, we simply find any cell with a contour line within it and successively apply the preceding rules. The contour line can only cross a cell-face once, so maintaining a simple bitmap of all crossed cell-faces can prevent endless processing of closed contours.
Of course, in the real world contours often cross the edge of the data field. To handle this, each border is scanned for inward contour lines, and these contours are followed until another border is encountered. The aforementioned bitmap of cell-faces is not cleared as these contour lines are generated. All interior cell-faces are then scanned for unmarked contour lines. These form closed contour lines and are followed until the contour returns to its initial point.
After all interior points have been checked, the bitmap is cleared and the next level (if any) is contoured.
Listing One (page 91) shows the contouring routine. Listing Two (page 95) contains a simple PostScript file generator illustrating the requirements of Polyline() and Text(). Finally, Listing Three (page 95) contains a sample application that generates a data field and then contours it.
The contour routine is called with four arguments: a pointer to the data field itself, its width and height, and the contouring interval. In the accompanying code floats are used to save space, although when contouring a nearly flat field it may be necessary to promote all floats to doubles for accuracy.
The data field is assumed to be stored in row-major order, starting at the upper-left corner. All data points must contain valid data, although the extensions to allow for missing data are straightforward.
The application writer must provide two extra routines, called by the contouring module. The first is Polyline(), which is called with a list of points on the contour line. All coordinates are within the range 0 to 1.
The second routine is Text(), which is used to label the contour lines. It is called with a text string and two coordinates, once again within the range 0 to 1.
As mentioned before, modifying the code to allow for missing data and/or irregular edges is relatively straightforward. The primary headache is modifying startEdge() to follow a jagged border.
In many cases, nonrectangular grids are trivial to support. As long as each contour cell has four corners, the precise topology of the data points is irrelevant; simply remap the coordinates within Polyline() and Text(). For triangular cells the algorithm is even simpler than for rectangular cells, as there are fewer cases or rotations to deal with.
Filled contours, similar to the weather map in USA Today, are straightforward in concept. Perform a polygon fill (as described in Michael Abrash's February 1991 "Graphics Programming" column) on the polygon formed by the current contour line connected to all interior contour lines. (For example, contour lines A, B, and C would form the larger polygon a..ab..bac..ca where x..x is all points within contour line X.) Unfortunately, it requires substantial analysis to determine which contour lines lie "within" another contour line; see Figure 3 and Figure 4.
Determining the area within a contour line can be accomplished by appropriately modifying the preceding modification. Simply count the number of pixels which would have been set had the fill been performed, and divide by the total number of pixels within the display.
_CONTOURING DATA FIELDS_ by Bruce (Bear) Giles[LISTING ONE]
#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <values.h>
#if defined (NEVER)
#include <ieeefp.h>
#else
#define NaN 0xFFFFFFFF
#define isnanf(x) ((x) == NaN)
#endif
typedef unsigned short ushort;
typedef unsigned char uchar;
#define DEFAULT_LEVELS 16
/* Mnemonics for contour line bearings */
#define EAST 0
#define NORTH 1
#define WEST 2
#define SOUTH 3
/* Mnemonics for relative data point positions */
#define SAME 0
#define NEXT 1
#define OPPOSITE 2
#define ADJACENT 3
/* Bit-mapped information in 'map' field. */
#define EW_MAP 0x01
#define NS_MAP 0x02
typedef struct
{
float x, y;
} LIST;
typedef struct
{
short dim_x; /* dimensions of grid array... */
short dim_y;
float max_value; /* statistics on the data... */
float min_value;
float mean;
float std;
short contour_mode; /* control variable */
float first_level; /* first (and subsequent) */
float step; /* contour level */
char format[20]; /* format of contour labels */
float *data; /* pointer to grid data */
char *map; /* pointer to "in-use" map */
LIST *list; /* used by 'Polyline()' */
ushort count;
} GRID;
typedef struct
{
short x;
short y;
uchar bearing;
} POINT;
#define MXY_to_L(g,x,y) ((ushort) (y) * (g)->dim_x + (ushort) (x) + 1)
#define XY_to_L(g,x,y) ((ushort) (y) * (g)->dim_x + (ushort) (x))
extern void Text();
extern void Polyline();
/* Contour generation. */
void Contour();
int scaleData();
void startLine();
void startInterior();
void drawLine();
void markInUse();
uchar faceInUse();
float getDataPoint();
void initPoint();
void lastPoint();
uchar savePoint();
/* inc > 0 increment to use between contour levels.
* < 0 number of contour levels to generate [abs(inc)].
* = 0 generate default number of contour levels.
*/
void Contour (data, dim_x, dim_y, inc)
float *data;
int dim_x;
int dim_y;
double inc;
{
GRID grid;
grid.data = data;
grid.dim_x = dim_x;
grid.dim_y = dim_y;
/* Allocate buffers used to contain contour information */
if ((grid.map = malloc ((dim_x + 1) * dim_y)) == NULL)
{
fprintf (stderr, "Contour(): unable to allocate buffer! (%d bytes)\n",
(dim_x + 1) * dim_y * sizeof (LIST));
free ((char *) grid.map);
return;
}
grid.list = (LIST *) malloc (2 * dim_x * dim_y * sizeof (LIST));
if (grid.list == (LIST *) NULL)
{
fprintf (stderr, "Contour(): unable to allocate buffer! (%d bytes)\n",
2 * dim_x * dim_y * sizeof (LIST));
free ((char *) grid.map);
return;
}
/* Generate contours, if not a uniform field. */
if (scaleData (&grid, inc))
startLine (&grid);
/* Release data structures. */
free ((char *) grid.map);
free ((char *) grid.list);
}
/* scaleData--Determine necessary statistics for contouring data set: global
* maximum & minimum, etc. Then initialize items used by rest of module. */
int scaleData (grid, inc)
GRID *grid;
double inc;
{
ushort i;
float step, level;
float sum, sum2, count;
float p, *u, *v, r;
char *s;
short n1, n2;
int first, n;
long x;
sum = sum2 = count = 0.0;
first = 1;
s = grid->map;
u = grid->data;
v = u + grid->dim_x * grid->dim_y;
for (i = 0; i < grid->dim_x * grid->dim_y; i++, u++, v++, s++)
{
r = *u;
sum += r;
sum2 += r * r;
count += 1.0;
if (first)
{
grid->max_value = grid->min_value = r;
first = 0;
}
else if (grid->max_value < r)
grid->max_value = r;
else if (grid->min_value > r)
grid->min_value = r;
}
grid->mean = sum / count;
if (grid->min_value == grid->max_value)
return 0;
grid->std = sqrt ((sum2 - sum * sum /count) / (count - 1.0));
if (inc > 0.0)
{
/* Use specified increment */
step = inc;
n = (int) (grid->max_value - grid->min_value) / step + 1;
while (n > 40)
{
step *= 2.0;
n = (int) (grid->max_value - grid->min_value) / step + 1;
}
}
else
{
/* Choose [specified|reasonable] number of levels and normalize
* increment to a reasonable value. */
n = (inc == 0.0) ? DEFAULT_LEVELS : (short) fabs (inc);
step = 4.0 * grid->std / (float) n;
p = pow (10.0, floor (log10 ((double) step)));
step = p * floor ((step + p / 2.0) / p);
}
n1 = (int) floor (log10 (fabs (grid->max_value)));
n2 = -((int) floor (log10 (step)));
if (n2 > 0)
sprintf (grid->format, "%%%d.%df", n1 + n2 + 2, n2);
else
sprintf (grid->format, "%%%d.0f", n1 + 1);
if (grid->max_value * grid->min_value < 0.0)
level = step * floor (grid->mean / step); /* odd */
else
level = step * floor (grid->min_value / step);
level -= step * floor ((float) (n - 1)/ 2);
/* Back up to include add'l levels, if necessary */
while (level - step > grid->min_value)
level -= step;
grid->first_level = level;
grid->step = step;
return 1;
}
/* startLine -- Locate first point of contour lines by checking edges of
gridded data set, then interior points, for each contour level. */
static
void startLine (grid)
GRID *grid;
{
ushort idx, i, edge;
double level;
for (idx = 0, level = grid->first_level; level < grid->max_value;
level += grid->step, idx++)
{
/* Clear flags */
grid->contour_mode = (level >= grid->mean);
memset (grid->map, 0, grid->dim_x * grid->dim_y);
/* Check edges */
for (edge = 0; edge < 4; edge++)
startEdge (grid, level, edge);
/* Check interior points */
startInterior (grid, level);
}
}
/* startEdge -- For a specified contour level and edge of gridded data set,
* check for (properly directed) contour line. */
static
void startEdge (grid, level, bearing)
GRID *grid;
float level;
uchar bearing;
{
POINT point1, point2;
float last, next;
short i, ds;
switch (point1.bearing = bearing)
{
case EAST:
point1.x = 0;
point1.y = 0;
ds = 1;
break;
case NORTH:
point1.x = 0;
point1.y = grid->dim_y - 2;
ds = 1;
break;
case WEST:
point1.x = grid->dim_x - 2;
point1.y = grid->dim_y - 2;
ds = -1;
break;
case SOUTH:
point1.x = grid->dim_x - 2;
point1.y = 0;
ds = -1;
break;
}
switch (point1.bearing)
{
/* Find first point with valid data. */
case EAST:
case WEST:
next = getDataPoint (grid, &point1, SAME);
memcpy ((char *) &point2, (char *) &point1, sizeof (POINT));
point2.x -= ds;
for (i = 1; i < grid->dim_y; i++, point1.y = point2.y += ds)
{
last = next;
next = getDataPoint (grid, &point1, NEXT);
if (last >= level && level > next)
{
drawLine (grid, &point1, level);
memcpy ((char *) &point1, (char *) &point2, sizeof (POINT));
point1.x = point2.x + ds;
}
}
break;
/* Find first point with valid data. */
case SOUTH:
case NORTH:
next = getDataPoint (grid, &point1, SAME);
memcpy ((char *) &point2, (char *) &point1, sizeof (POINT));
point2.y += ds;
for (i = 1; i < grid->dim_x; i++, point1.x = point2.x += ds)
{
last = next;
next = getDataPoint (grid, &point1, NEXT);
if (last >= level && level > next)
{
drawLine (grid, &point1, level);
memcpy ((char *) &point1, (char *) &point2, sizeof (POINT));
point1.y = point2.y - ds;
}
}
break;
}
}
/* startInterior -- For a specified contour level, check for (properly
directed) contour line for all interior data points. Do _not_ follow contour
lines detected by the 'startEdge' routine. */
static
void startInterior (grid, level)
GRID *grid;
float level;
{
POINT point;
ushort x, y;
float next, last;
for (x = 1; x < grid->dim_x - 1; x++)
{
point.x = x;
point.y = 0;
point.bearing = EAST;
next = getDataPoint (grid, &point, SAME);
for (y = point.y; y < grid->dim_y; y++, point.y++)
{
last = next;
next = getDataPoint (grid, &point, NEXT);
if (last >= level && level > next)
{
if (!faceInUse (grid, &point, WEST))
{
drawLine (grid, &point, level);
point.x = x;
point.y = y;
point.bearing = EAST;
}
}
}
}
}
/* drawLine -- Given in initial contour point by either 'startEdge' or
'startInterior', follow contour line until it encounters either an edge
or previously contoured cell. */
static
void drawLine (grid, point, level)
GRID *grid;
POINT *point;
float level;
{
uchar exit_bearing;
uchar adj, opp;
float fadj, fopp;
initPoint (grid);
for ( ;; )
{
/* Add current point to vector list. If either of the points is
* missing, return immediately (open contour). */
if (!savePoint (grid, point, level))
{
lastPoint (grid);
return;
}
/* Has this face of this cell been marked in use? If so, then this is
* a closed contour. */
if (faceInUse (grid, point, WEST))
{
lastPoint (grid);
return;
}
/* Examine adjacent and opposite corners of cell; determine
* appropriate action. */
markInUse (grid, point, WEST);
fadj = getDataPoint (grid, point, ADJACENT);
fopp = getDataPoint (grid, point, OPPOSITE);
/* If either point is missing, return immediately (open contour). */
if (isnanf (fadj) || isnanf (fopp))
{
lastPoint (grid);
return;
}
adj = (fadj <= level) ? 2 : 0;
opp = (fopp >= level) ? 1 : 0;
switch (adj + opp)
{
/* Exit EAST face. */
case 0:
markInUse (grid, point, NORTH);
markInUse (grid, point, SOUTH);
exit_bearing = EAST;
break;
/* Exit SOUTH face. */
case 1:
markInUse (grid, point, NORTH);
markInUse (grid, point, EAST);
exit_bearing = SOUTH;
break;
/* Exit NORTH face. */
case 2:
markInUse (grid, point, EAST);
markInUse (grid, point, SOUTH);
exit_bearing = NORTH;
break;
/* Exit NORTH or SOUTH face, depending upon contour level. */
case 3:
exit_bearing = (grid->contour_mode) ? NORTH : SOUTH;
break;
}
/* Update face number, coordinate of defining corner. */
point->bearing = (point->bearing + exit_bearing) % 4;
switch (point->bearing)
{
case EAST : point->x++; break;
case NORTH: point->y--; break;
case WEST : point->x--; break;
case SOUTH: point->y++; break;
}
}
}
/* markInUse -- Mark the specified cell face as contoured. This is necessary to
* prevent infinite processing of closed contours. see also: faceInUse */
static
void markInUse (grid, point, face)
GRID *grid;
POINT *point;
uchar face;
{
face = (point->bearing + face) % 4;
switch (face)
{
case NORTH:
case SOUTH:
grid->map[MXY_to_L (grid,
point->x, point->y + (face == SOUTH ? 1 : 0))] |= NS_MAP;
break;
case EAST:
case WEST:
grid->map[MXY_to_L (grid,
point->x + (face == EAST ? 1 : 0), point->y)] |= EW_MAP;
break;
}
}
/* faceInUse -- Determine if the specified cell face has been marked as
* contoured. This is necessary to prevent infinite processing of closed
* contours. see also: markInUse */
static
uchar faceInUse (grid, point, face)
GRID *grid;
POINT *point;
uchar face;
{
uchar r;
face = (point->bearing + face) % 4;
switch (face)
{
case NORTH:
case SOUTH:
r = grid->map[MXY_to_L (grid,
point->x, point->y + (face == SOUTH ? 1 : 0))] & NS_MAP;
break;
case EAST:
case WEST:
r = grid->map[MXY_to_L (grid,
point->x + (face == EAST ? 1 : 0), point->y)] & EW_MAP;
break;
}
return r;
}
/* initPoint -- Initialize the contour point list.
* see also: savePoint, lastPoint */
static
void initPoint (grid)
GRID *grid;
{
grid->count = 0;
}
/* lastPoint -- Generate the actual contour line from the contour point list.
* see also: savePoint, initPoint */
static
void lastPoint (grid)
GRID *grid;
{
if (grid->count)
Polyline (grid->count, grid->list);
}
/* savePoints -- Add specified point to the contour point list.
* see also: initPoint, lastPoint */
static
uchar savePoint (grid, point, level)
GRID *grid;
POINT *point;
float level;
{
float last, next;
float x, y, ds;
char s[80];
static int cnt = 0;
last = getDataPoint (grid, point, SAME);
next = getDataPoint (grid, point, NEXT);
/* Are the points the same value? */
if (last == next)
{
fprintf (stderr, "(%2d, %2d, %d) ", point->x, point->y,
point->bearing);
fprintf (stderr, "%8g %8g ", last, next);
fprintf (stderr, "potential divide-by-zero!\n");
return 0;
}
x = (float) point->x;
y = (float) point->y;
ds = (float) ((last - level) / (last - next));
switch (point->bearing)
{
case EAST : y += ds; break;
case NORTH: x += ds; y += 1.0; break;
case WEST : x += 1.0; y += 1.0 - ds; break;
case SOUTH: x += 1.0 - ds; break;
}
/* Update to contour point list */
grid->list[grid->count].x = x / (float) (grid->dim_x - 1);
grid->list[grid->count].y = y / (float) (grid->dim_y - 1);
/* Add text label to contour line. */
if (!(cnt++ % 11))
{
sprintf (s, grid->format, level);
Text (s, grid->list[grid->count].x, grid->list[grid->count].y);
}
/* Update counter */
grid->count++;
return 1;
}
/* getDataPoint -- Return the value of the data point in specified corner of
* specified cell (the 'point' parameter contains the address of the
* top-left corner of this cell). */
static
float getDataPoint (grid, point, corner)
GRID *grid;
POINT *point;
uchar corner;
{
ushort dx, dy;
ushort offset;
switch ((point->bearing + corner) % 4)
{
case SAME : dx = 0; dy = 0; break;
case NEXT : dx = 0; dy = 1; break;
case OPPOSITE: dx = 1; dy = 1; break;
case ADJACENT: dx = 1; dy = 0; break;
}
offset = XY_to_L (grid, point->x + dx, point->y + dy);
if ((short) (point->x + dx) >= grid->dim_x ||
(short) (point->y + dy) >= grid->dim_y ||
(short) (point->x + dx) < 0 ||
(short) (point->y + dy) < 0)
{
return NaN;
}
else
return grid->data[offset];
}
[LISTING TWO]
#include <stdio.h>
typedef struct
{
float x, y;
} LIST;
void Polyline (n, list)
int n;
LIST *list;
{
short x0, x1, y0, y1;
if (n < 2)
return;
printf ("newpath\n");
printf ("%.6f %.6f moveto\n", list->x, 1.0 - list->y);
list++;
while (--n)
{
printf ("%.6f %.6f lineto\n", list->x, 1.0 - list->y);
list++;
}
printf ("stroke\n\n");
}
void Text (s, x, y)
char *s;
float x, y;
{
printf ("%.6f %.6f moveto (%s) show\n", x, 1.0 - y, s);
}
void psOpen ()
{
printf ("%%!\n");
printf ("save\n");
printf ("\n");
printf ("/Helvetica findfont 0.015 scalefont setfont\n");
printf ("\n");
printf ("72 252 translate\n");
printf ("468 468 scale\n");
printf ("0.001 setlinewidth\n");
printf ("\n");
printf ("newpath\n");
printf (" 0 0 moveto\n");
printf (" 0 1 lineto\n");
printf (" 1 1 lineto\n");
printf (" 1 0 lineto\n");
printf (" closepath\n");
printf ("stroke\n");
printf ("clippath\n");
printf ("\n");
printf ("0.00001 setlinewidth\n");
printf ("\n");
}
void psClose ()
{
printf ("restore\n");
printf ("showpage\n");
[LISTING THREE]
#include <stdio.h>
#include <math.h>
float data[400];
void main (argc, argv)
int argc;
char **argv;
{
float *s;
int i, j;
double x, y, r1, r2;
s = data;
for (i = 0, s = data; i < 20; i++)
for (j = 0; j < 20; j++)
{
x = ((double) i - 9.5) / 4.0;
y = ((double) j - 9.5) / 4.0;
*s++ = (float) (10.0 * cos(x) * cos(y));
}
psOpen ();
Contour (data, 20, 20, 2.0);
psClose ();
}
Copyright © 1992, Dr. Dobb's Journal
Terms of Service | Privacy Statement | Copyright © 2024 UBM Tech, All rights reserved.