home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
OS/2 Shareware BBS: Science
/
Science.zip
/
gmt_os2.zip
/
src
/
triangulate.c
< prev
next >
Wrap
C/C++ Source or Header
|
1995-06-09
|
12KB
|
357 lines
/*--------------------------------------------------------------------
* The GMT-system: @(#)triangulate.c 2.16 09 Jun 1995
*
* Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
* See README file for copying and redistribution conditions.
*--------------------------------------------------------------------*/
/*
* triangulate reads one or more files (or stdin) with x,y[,whatever] and
* outputs the indeces of the vertices of the optimal Delaunay triangulation
* using the method by Watson, D. F., ACORD: Automatic contouring of raw data,
* Computers & Geosciences, 8, 97-101, 1982. Optionally, the output may take
* the form of (1) a multi-segment file with the vertice coordinates needed to
* draw the triangles, or (2) a grdfile based on gridding the plane estimates
*
* Author: Paul Wessel
* Date: 1-JAN-1993
* Version: 2.0
*
*/
#include "gmt.h"
struct EDGE {
int begin, end;
} *edge;
main (argc, argv)
int argc;
char **argv; {
int i, j, ij, ij1, ij2, ij3, np, k, fno, n = 0, n_alloc = 0, n_files = 0, n_edge;
int n_expected_fields = 2, x, y, n_args, *link, n_fields, compare_edge();
int one_or_zero, i_min, i_max, j_min, j_max, p, dummy[4];
BOOLEAN error = FALSE, map_them = FALSE, multi_segments = FALSE, single_precision = FALSE;
BOOLEAN nofile = TRUE, done = FALSE, b_in = FALSE, b_out = FALSE, do_grid = FALSE, set_empty = FALSE, more;
double west = 0.0, east = 0.0, south = 0.0, north = 0.0, empty = 0.0, zj, zk, zl, zlj, zkj, *in;
double *xx, *yy, *zz, vx[4], vy[4], xinc2, yinc2, idx, idy, xp, yp, a, b, c, f;
double xkj, xlj, ykj, ylj;
float *grd;
char line[512], format[512], EOL_flag = '>', *outfile = 0;
FILE *fp = NULL;
struct GRD_HEADER header;
argc = gmt_begin (argc, argv);
for (i = 1; i < argc; i++) {
if (argv[i][0] == '-') {
switch (argv[i][1]) {
/* Common parameters */
case 'H':
case 'J':
case 'R':
case 'V':
case ':':
case '\0':
error += get_common_args (argv[i], &west, &east, &south, &north);
break;
/* Supplemental parameters */
case 'b':
single_precision = !(argv[i][strlen(argv[i])-1] == 'd');
if (argv[i][2] == 'i')
b_in = TRUE;
else if (argv[i][2] == 'o')
b_out = TRUE;
else
b_in = b_out = TRUE;
break;
case 'E':
empty = atof (&argv[i][2]);
set_empty = TRUE;
break;
case 'F':
header.node_offset = TRUE;
break;
case 'G':
outfile = &argv[i][2];
break;
case 'I':
gmt_getinc (&argv[i][2], &header.x_inc, &header.y_inc);
break;
case 'M':
multi_segments = TRUE;
if (argv[i][2]) EOL_flag = argv[i][2];
break;
case 'Z': /* Z-column present */
n_expected_fields = 3;
break;
default:
error = TRUE;
gmt_default_error (argv[i][1]);
break;
}
}
else
n_files++;
}
if (argc == 1 || gmt_quick) {
fprintf (stderr, "triangulate %s - Delaunay triangulation\n\n", GMT_VERSION);
fprintf (stderr, "usage: triangulate <infiles> [-F] [-G<grdfile> [-H] [-I<dx>[m|c][/<dy>[m|c]]] [-J<parameters>]\n");
fprintf (stderr, " [-M[<flag>]] [-R<west/east/south/north>] [-V] [-:] [-Z] [-b[i|o][d]]\n\n");
if (gmt_quick) exit (-1);
fprintf (stderr, " infiles (in ASCII) has 2 or more columns. If no file(s) is given, standard input is read.\n");
fprintf (stderr, "\n\tOPTIONS:\n");
fprintf(stderr, " -F Force pixel registration (only with -G) [Default is gridline registration]\n");
fprintf(stderr, " -G Grid data. Give name of output gridfile and specify -R -I\n");
explain_option ('H');
fprintf(stderr, " -I sets the grid spacing for the grid. Append m for minutes, c for seconds\n");
explain_option ('J');
fprintf (stderr, " -M output triangle edges as multiple segments separated by a record\n");
fprintf (stderr, " whose first character is <flag> ['>']\n");
fprintf (stderr, " [Default is to output the indeces of vertices for each Delaunay triangle]\n");
explain_option ('R');
explain_option ('V');
fprintf (stderr, "\t-Z means a third column (z) is present in binary input (see -b)\n");
explain_option (':');
fprintf (stderr, "\t-b means binary i/o. Append i or o if only input OR output is binary\n");
fprintf (stderr, "\t Append d for double precision [Default is single precision].\n"); explain_option ('.');
fprintf (stderr, "\t -bo is ignored if -M is selected [Default is ASCII i/o]\n");
explain_option ('.');
exit (-1);
}
if (b_in && gmtdefs.io_header) {
fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data cannot have header -H\n", gmt_program);
error++;
}
do_grid = (outfile && header.x_inc > 0.0 && header.y_inc > 0.0);
if (do_grid && (west == east || south == north)) {
fprintf (stderr, "%s: GMT SYNTAX ERROR. Must specify -R, -I, -G for gridding\n", gmt_program);
error++;
}
if (do_grid && b_in && n_expected_fields == 2) {
fprintf (stderr, "%s: GMT SYNTAX ERROR. Must specify -Z for gridding of binary xyz data\n", gmt_program);
error++;
}
if (error) exit (-1);
if (b_in) gmt_input = (single_precision) ? bin_float_input : bin_double_input;
if (b_out) gmt_output = (single_precision) ? bin_float_output : bin_double_output;
if (west != east && !do_grid) {
map_them = TRUE;
map_setup (west, east, south, north);
}
/* Now we are ready to take on some input values */
x = (gmtdefs.xy_toggle) ? 1 : 0; y = 1 - x; /* Set up which columns have x and y */
if (n_files > 0)
nofile = FALSE;
else
n_files = 1;
n_args = (argc > 1) ? argc : 2;
n_alloc = GMT_CHUNK;
xx = (double *) memory (CNULL, n_alloc, sizeof (double), "triangulate");
yy = (double *) memory (CNULL, n_alloc, sizeof (double), "triangulate");
if (do_grid) {
zz = (double *) memory (CNULL, n_alloc, sizeof (double), "triangulate");
n_expected_fields = 3;
}
n = 0;
for (fno = 1; !done && fno < n_args; fno++) { /* Loop over input files, if any */
if (!nofile && argv[fno][0] == '-') continue;
if (nofile) { /* Just read standard input */
fp = stdin;
done = TRUE;
}
else if ((fp = fopen (argv[fno], "r")) == NULL) {
fprintf (stderr, "triangulate: Cannot open file %s\n", argv[fno]);
continue;
}
if (!nofile && gmtdefs.verbose) fprintf (stderr, "triangulate: Reading file %s\n", argv[fno]);
if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (line, 512, fp);
more = ((n_fields = gmt_input (fp, n_expected_fields, &in)) == n_expected_fields);
while (more) {
if (map_them)
geo_to_xy (in[x], in[y], &xx[n], &yy[n]);
else {
xx[n] = in[x];
yy[n] = in[y];
}
if (do_grid) zz[n] = in[2];
n++;
if (n == n_alloc) { /* Get more memory */
n_alloc += GMT_CHUNK;
xx = (double *) memory ((char *)xx, n_alloc, sizeof (double), "triangulate");
yy = (double *) memory ((char *)yy, n_alloc, sizeof (double), "triangulate");
if (do_grid) zz = (double *) memory ((char *)zz, n_alloc, sizeof (double), "triangulate");
}
more = ((n_fields = gmt_input (fp, n_expected_fields, &in)) == n_expected_fields);
}
if (fp != stdin) fclose (fp);
if (n_fields == -1) n_fields = 0; /* Ok to get EOF */
if (n_fields%n_expected_fields) { /* Got garbage or multiple segment subheader */
fprintf (stderr, "%s: Cannot read line %d, aborts\n", gmt_program, n);
exit (-1);
}
}
xx = (double *) memory ((char *)xx, n, sizeof (double), "triangulate");
yy = (double *) memory ((char *)yy, n, sizeof (double), "triangulate");
if (do_grid) zz = (double *) memory ((char *)zz, n, sizeof (double), "triangulate");
if (gmtdefs.verbose) fprintf (stderr, "triangulate: Do Delaunay optimal triangulation\n");
np = delaunay (xx, yy, n, &link);
if (gmtdefs.verbose) fprintf (stderr, "triangulate: %d Delaunay triangles found\n", np);
if (do_grid) {
grd_init (&header, argc, argv, TRUE);
header.x_min = west; header.x_max = east;
header.y_min = south; header.y_max = north;
if (header.node_offset) {
one_or_zero = 0;
xinc2 = 0.5 * header.x_inc;
yinc2 = 0.5 * header.y_inc;
}
else {
one_or_zero = 1;
xinc2 = yinc2 = 0.0;
}
idx = 1.0 / header.x_inc;
idy = 1.0 / header.y_inc;
header.nx = rint ( (header.x_max - header.x_min) * idx) + one_or_zero;
header.ny = rint ( (header.y_max - header.y_min) * idy) + one_or_zero;
grd = (float *) memory (CNULL, (int)(header.nx * header.ny), sizeof (float), "triangulate");
if (!set_empty) empty = gmt_NaN;
for (i = 0; i < header.nx * header.ny; i++) grd[i] = empty; /* initialize grid */
for (k = ij = 0; k < np; k++) {
/* Find equation for the plane as z = ax + by + c */
vx[0] = vx[3] = xx[link[ij]]; vy[0] = vy[3] = yy[link[ij]]; zj = zz[link[ij++]];
vx[1] = xx[link[ij]]; vy[1] = yy[link[ij]]; zk = zz[link[ij++]];
vx[2] = xx[link[ij]]; vy[2] = yy[link[ij]]; zl = zz[link[ij++]];
xkj = vx[1] - vx[0]; ykj = vy[1] - vy[0];
zkj = zk - zj; xlj = vx[2] - vx[0];
ylj = vy[2] - vy[0]; zlj = zl - zj;
f = 1.0 / (xkj * ylj - ykj * xlj);
a = -f * (ykj * zlj - zkj * ylj);
b = -f * (zkj * xlj - xkj * zlj);
c = -a * vx[1] - b * vy[1] + zk;
i_min = floor ((MIN (MIN (vx[0], vx[1]), vx[2]) - header.x_min + xinc2) * idx);
i_max = ceil ((MAX (MAX (vx[0], vx[1]), vx[2]) - header.x_min + xinc2) * idx);
j_min = floor ((header.y_max - MAX (MAX (vy[0], vy[1]), vy[2]) + yinc2) * idy);
j_max = ceil ((header.y_max - MIN (MIN (vy[0], vy[1]), vy[2]) + yinc2) * idy);
for (j = j_min; j <= j_max; j++) {
yp = header.y_max - j * header.y_inc - yinc2;
p = j * header.nx + i_min;
for (i = i_min; i <= i_max; i++, p++) {
xp = header.x_min + i * header.x_inc + xinc2;
if (!non_zero_winding (xp, yp, vx, vy, 4)) continue; /* Outside */
grd[p] = a * xp + b * yp + c;
}
}
}
dummy[0] = dummy[1] = dummy[2] = dummy[3] = 0;
if (write_grd (outfile, &header, grd, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
fprintf (stderr, "triangulate: Error writing file %s\n", outfile);
exit (-1);
}
}
if (multi_segments) { /* Must find unique edges to output only once */
n_edge = 3 * np;
edge = (struct EDGE *) memory (CNULL, n_edge, sizeof (struct EDGE), "triangulate");
for (i = ij1 = 0, ij2 = 1, ij3 = 2; i < np; i++, ij1 += 3, ij2 += 3, ij3 += 3) {
edge[ij1].begin = link[ij1]; edge[ij1].end = link[ij2];
edge[ij2].begin = link[ij2]; edge[ij2].end = link[ij3];
edge[ij3].begin = link[ij1]; edge[ij3].end = link[ij3];
}
qsort ((char *)edge, n_edge, sizeof (struct EDGE), compare_edge);
for (i = 1, j = 0; i < n_edge; i++) {
if (edge[i].begin != edge[j].begin || edge[i].end != edge[j].end) j++;
edge[j] = edge[i];
}
n_edge = j + 1;
if (gmtdefs.verbose) fprintf (stderr, "triangulate: %d unique triangle edges\n", n_edge);
sprintf (format, "%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format);
for (i = 0; i < n_edge; i++) {
printf ("%c Edge %d-%d\n", EOL_flag, edge[i].begin, edge[i].end);
printf (format, xx[edge[i].begin], yy[edge[i].begin]);
printf (format, xx[edge[i].end], yy[edge[i].end]);
}
free ((char *)edge);
}
else if (b_out)
fwrite ((char *)link, sizeof (int), 3*np, stdout);
else
for (i = ij = 0; i < np; i++, ij += 3) printf ("%d\t%d\t%d\n", link[ij], link[ij+1], link[ij+2]);
free ((char *) xx);
free ((char *) yy);
free ((char *) link);
if (gmtdefs.verbose) fprintf (stderr, "triangulate: Done!\n");
gmt_end (argc, argv);
}
int compare_edge (a, b)
struct EDGE *a, *b; {
if (a->begin < b->begin)
return (-1);
else if (a->begin > b->begin)
return (1);
else {
if (a->end < b->end)
return (-1);
else if (a->end > b->end)
return (1);
else
return (0);
}
}