_CONTOURING DATA FIELDS_ by Bruce (Bear) Giles [LISTING ONE] #include #include #include #include #if defined (NEVER) #include #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 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 #include 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 (); }