home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Power-Programmierung
/
CD1.mdf
/
magazine
/
drdobbs
/
1992
/
06
/
contour.asc
< prev
next >
Wrap
Text File
|
1992-05-18
|
19KB
|
694 lines
_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 ();
}