Category : Files from Magazines
Archive : DDJ0692.ZIP
Filename : CONTOUR.ASC
Output of file : CONTOUR.ASC contained in archive : DDJ0692.ZIP
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 ();
}
Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!
This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.
But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/