home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / triangulate.c < prev    next >
C/C++ Source or Header  |  1995-06-09  |  12KB  |  357 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:   @(#)triangulate.c    2.16  09 Jun 1995
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  * triangulate reads one or more files (or stdin) with x,y[,whatever] and
  9.  * outputs the indeces of the vertices of the optimal Delaunay triangulation
  10.  * using the method by Watson, D. F., ACORD: Automatic contouring of raw data,
  11.  * Computers & Geosciences, 8, 97-101, 1982.  Optionally, the output may take
  12.  * the form of (1) a multi-segment file with the vertice coordinates needed to
  13.  * draw the triangles, or (2) a grdfile based on gridding the plane estimates
  14.  *
  15.  * Author:      Paul Wessel
  16.  * Date:        1-JAN-1993
  17.  * Version:     2.0
  18.  *
  19.  */
  20.  
  21. #include "gmt.h"
  22.  
  23. struct EDGE {
  24.     int begin, end;
  25. } *edge;
  26.     
  27. main (argc, argv)
  28. int argc;
  29. char **argv; {
  30.     int i, j, ij, ij1, ij2, ij3, np, k, fno, n = 0, n_alloc = 0, n_files = 0, n_edge;
  31.     int n_expected_fields = 2, x, y, n_args, *link, n_fields, compare_edge();
  32.     int one_or_zero, i_min, i_max, j_min, j_max, p, dummy[4];
  33.     
  34.     BOOLEAN error = FALSE, map_them = FALSE, multi_segments = FALSE, single_precision = FALSE;
  35.     BOOLEAN nofile = TRUE, done = FALSE, b_in = FALSE, b_out = FALSE, do_grid = FALSE, set_empty = FALSE, more;
  36.     
  37.     double west = 0.0, east = 0.0, south = 0.0, north = 0.0, empty = 0.0, zj, zk, zl, zlj, zkj, *in;
  38.     double *xx, *yy, *zz, vx[4], vy[4], xinc2, yinc2, idx, idy, xp, yp, a, b, c, f;
  39.     double xkj, xlj, ykj, ylj;
  40.     
  41.     float *grd;
  42.     
  43.     char line[512], format[512], EOL_flag = '>', *outfile = 0;
  44.     
  45.     FILE *fp = NULL;
  46.     
  47.     struct GRD_HEADER header;
  48.  
  49.     argc = gmt_begin (argc, argv);
  50.  
  51.     for (i = 1; i < argc; i++) {
  52.         if (argv[i][0] == '-') {
  53.             switch (argv[i][1]) {
  54.         
  55.                 /* Common parameters */
  56.         
  57.                 case 'H':
  58.                 case 'J':
  59.                 case 'R':
  60.                 case 'V':
  61.                 case ':':
  62.                 case '\0':
  63.                     error += get_common_args (argv[i], &west, &east, &south, &north);
  64.                     break;
  65.         
  66.                 /* Supplemental parameters */
  67.         
  68.                 case 'b':
  69.                     single_precision = !(argv[i][strlen(argv[i])-1] == 'd');
  70.                     if (argv[i][2] == 'i')
  71.                         b_in = TRUE;
  72.                     else if (argv[i][2] == 'o')
  73.                         b_out = TRUE;
  74.                     else
  75.                         b_in = b_out = TRUE;
  76.                     break;
  77.                 case 'E':
  78.                     empty = atof (&argv[i][2]);
  79.                     set_empty = TRUE;
  80.                     break;
  81.                 case 'F':
  82.                     header.node_offset = TRUE;
  83.                     break;
  84.                 case 'G':
  85.                     outfile = &argv[i][2];
  86.                     break;
  87.                 case 'I':
  88.                     gmt_getinc (&argv[i][2], &header.x_inc, &header.y_inc);
  89.                     break;
  90.                 case 'M':
  91.                     multi_segments = TRUE;
  92.                     if (argv[i][2]) EOL_flag = argv[i][2];
  93.                     break;
  94.                 case 'Z':    /* Z-column present */
  95.                     n_expected_fields = 3;
  96.                     break;
  97.                 default:
  98.                     error = TRUE;
  99.                     gmt_default_error (argv[i][1]);
  100.                     break;
  101.             }
  102.         }
  103.         else
  104.             n_files++;
  105.     }
  106.  
  107.     if (argc == 1 || gmt_quick) {
  108.         fprintf (stderr, "triangulate %s - Delaunay triangulation\n\n", GMT_VERSION);
  109.         fprintf (stderr, "usage: triangulate <infiles> [-F] [-G<grdfile> [-H] [-I<dx>[m|c][/<dy>[m|c]]] [-J<parameters>]\n");
  110.         fprintf (stderr, "    [-M[<flag>]] [-R<west/east/south/north>] [-V] [-:] [-Z] [-b[i|o][d]]\n\n");
  111.                 
  112.         if (gmt_quick) exit (-1);
  113.                 
  114.         fprintf (stderr, "      infiles (in ASCII) has 2 or more columns.  If no file(s) is given, standard input is read.\n");
  115.         fprintf (stderr, "\n\tOPTIONS:\n");
  116.         fprintf(stderr, "    -F Force pixel registration (only with -G) [Default is gridline registration]\n");
  117.         fprintf(stderr, "    -G Grid data. Give name of output gridfile and specify -R -I\n");
  118.         explain_option ('H');
  119.         fprintf(stderr, "    -I sets the grid spacing for the grid.  Append m for minutes, c for seconds\n");
  120.         explain_option ('J');   
  121.         fprintf (stderr, "      -M output triangle edges as multiple segments separated by a record\n");
  122.         fprintf (stderr, "         whose first character is <flag> ['>']\n");
  123.         fprintf (stderr, "         [Default is to output the indeces of vertices for each Delaunay triangle]\n");
  124.         explain_option ('R');
  125.         explain_option ('V');
  126.         fprintf (stderr, "\t-Z means a third column (z) is present in binary input (see -b)\n");
  127.         explain_option (':');
  128.         fprintf (stderr, "\t-b means binary i/o.  Append i or o if only input OR output is binary\n");
  129.         fprintf (stderr, "\t   Append d for double precision [Default is single precision].\n");        explain_option ('.');
  130.         fprintf (stderr, "\t   -bo is ignored if -M is selected [Default is ASCII i/o]\n");
  131.         explain_option ('.');
  132.         exit (-1);
  133.     }
  134.     
  135.     if (b_in && gmtdefs.io_header) {
  136.         fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", gmt_program);
  137.         error++;
  138.     }
  139.     
  140.     do_grid = (outfile && header.x_inc > 0.0 && header.y_inc > 0.0);
  141.     
  142.     if (do_grid && (west == east || south == north)) {
  143.         fprintf (stderr, "%s: GMT SYNTAX ERROR.  Must specify -R, -I, -G for gridding\n", gmt_program);
  144.         error++;
  145.     }
  146.     
  147.     if (do_grid && b_in && n_expected_fields == 2) {
  148.         fprintf (stderr, "%s: GMT SYNTAX ERROR.  Must specify -Z for gridding of binary xyz data\n", gmt_program);
  149.         error++;
  150.     }
  151.     
  152.     if (error) exit (-1);
  153.  
  154.     if (b_in) gmt_input = (single_precision) ? bin_float_input : bin_double_input;
  155.     if (b_out) gmt_output = (single_precision) ? bin_float_output : bin_double_output;
  156.  
  157.     if (west != east && !do_grid) {
  158.         map_them = TRUE;
  159.         map_setup (west, east, south, north);
  160.     }
  161.  
  162.     /* Now we are ready to take on some input values */
  163.     
  164.     x = (gmtdefs.xy_toggle) ? 1 : 0;    y = 1 - x;        /* Set up which columns have x and y */
  165.     
  166.     if (n_files > 0)
  167.         nofile = FALSE;
  168.     else
  169.         n_files = 1;
  170.     n_args = (argc > 1) ? argc : 2;
  171.     
  172.     n_alloc = GMT_CHUNK;
  173.     xx = (double *) memory (CNULL, n_alloc, sizeof (double), "triangulate");
  174.     yy = (double *) memory (CNULL, n_alloc, sizeof (double), "triangulate");
  175.     if (do_grid) {
  176.         zz = (double *) memory (CNULL, n_alloc, sizeof (double), "triangulate");
  177.         n_expected_fields = 3;
  178.     }
  179.         
  180.     n = 0;
  181.     for (fno = 1; !done && fno < n_args; fno++) {    /* Loop over input files, if any */
  182.         if (!nofile && argv[fno][0] == '-') continue;
  183.         
  184.         if (nofile) {    /* Just read standard input */
  185.             fp = stdin;
  186.             done = TRUE;
  187.         }
  188.         else if ((fp = fopen (argv[fno], "r")) == NULL) {
  189.             fprintf (stderr, "triangulate: Cannot open file %s\n", argv[fno]);
  190.             continue;
  191.         }
  192.  
  193.         if (!nofile && gmtdefs.verbose) fprintf (stderr, "triangulate: Reading file %s\n", argv[fno]);
  194.         
  195.         if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (line, 512, fp);
  196.  
  197.         more = ((n_fields = gmt_input (fp,  n_expected_fields, &in)) == n_expected_fields);
  198.             
  199.         while (more) {
  200.                     
  201.             if (map_them)
  202.                 geo_to_xy (in[x], in[y], &xx[n], &yy[n]);
  203.             else {
  204.                 xx[n] = in[x];
  205.                 yy[n] = in[y];
  206.             }
  207.             if (do_grid) zz[n] = in[2];
  208.             n++;
  209.             
  210.             if (n == n_alloc) {    /* Get more memory */
  211.                 n_alloc += GMT_CHUNK;
  212.                 xx = (double *) memory ((char *)xx, n_alloc, sizeof (double), "triangulate");
  213.                 yy = (double *) memory ((char *)yy, n_alloc, sizeof (double), "triangulate");
  214.                 if (do_grid) zz = (double *) memory ((char *)zz, n_alloc, sizeof (double), "triangulate");
  215.             }
  216.             
  217.             more = ((n_fields = gmt_input (fp,  n_expected_fields, &in)) == n_expected_fields);
  218.         }
  219.         
  220.         if (fp != stdin) fclose (fp);
  221.         if (n_fields == -1) n_fields = 0;    /* Ok to get EOF */
  222.         if (n_fields%n_expected_fields) {    /* Got garbage or multiple segment subheader */
  223.             fprintf (stderr, "%s:  Cannot read line %d, aborts\n", gmt_program, n);
  224.             exit (-1);
  225.         }
  226.     }
  227.     
  228.     xx = (double *) memory ((char *)xx, n, sizeof (double), "triangulate");
  229.     yy = (double *) memory ((char *)yy, n, sizeof (double), "triangulate");
  230.     if (do_grid) zz = (double *) memory ((char *)zz, n, sizeof (double), "triangulate");
  231.  
  232.     if (gmtdefs.verbose) fprintf (stderr, "triangulate: Do Delaunay optimal triangulation\n");
  233.  
  234.     np = delaunay (xx, yy, n, &link);
  235.     
  236.     if (gmtdefs.verbose) fprintf (stderr, "triangulate: %d Delaunay triangles found\n", np);
  237.  
  238.     if (do_grid) {
  239.     
  240.         grd_init (&header, argc, argv, TRUE);
  241.  
  242.         header.x_min = west;    header.x_max = east;
  243.         header.y_min = south;    header.y_max = north;
  244.         if (header.node_offset) {
  245.             one_or_zero = 0;
  246.             xinc2 = 0.5 * header.x_inc;
  247.             yinc2 = 0.5 * header.y_inc;
  248.         }
  249.         else {
  250.             one_or_zero = 1;
  251.             xinc2 = yinc2 = 0.0;
  252.         }
  253.         idx = 1.0 / header.x_inc;
  254.         idy = 1.0 / header.y_inc;
  255.         header.nx = rint ( (header.x_max - header.x_min) * idx) + one_or_zero;
  256.         header.ny = rint ( (header.y_max - header.y_min) * idy) + one_or_zero;
  257.         grd = (float *) memory (CNULL, (int)(header.nx * header.ny), sizeof (float), "triangulate");
  258.         if (!set_empty) empty = gmt_NaN;
  259.         for (i = 0; i < header.nx * header.ny; i++) grd[i] = empty;    /* initialize grid */
  260.         
  261.         for (k = ij = 0; k < np; k++) {
  262.         
  263.             /* Find equation for the plane as z = ax + by + c */
  264.             
  265.             vx[0] = vx[3] = xx[link[ij]];    vy[0] = vy[3] = yy[link[ij]];    zj = zz[link[ij++]];
  266.             vx[1] = xx[link[ij]];    vy[1] = yy[link[ij]];    zk = zz[link[ij++]];
  267.             vx[2] = xx[link[ij]];    vy[2] = yy[link[ij]];    zl = zz[link[ij++]];
  268.             
  269.             xkj = vx[1] - vx[0];    ykj = vy[1] - vy[0];
  270.             zkj = zk - zj;    xlj = vx[2] - vx[0];
  271.             ylj = vy[2] - vy[0];    zlj = zl - zj;
  272.             
  273.             f = 1.0 / (xkj * ylj - ykj * xlj);
  274.             a = -f * (ykj * zlj - zkj * ylj);
  275.             b = -f * (zkj * xlj - xkj * zlj);
  276.             c = -a * vx[1] - b * vy[1] + zk;
  277.             
  278.             i_min = floor ((MIN (MIN (vx[0], vx[1]), vx[2]) - header.x_min + xinc2) * idx);
  279.             i_max = ceil ((MAX (MAX (vx[0], vx[1]), vx[2]) - header.x_min + xinc2) * idx);
  280.             j_min = floor ((header.y_max - MAX (MAX (vy[0], vy[1]), vy[2]) + yinc2) * idy);
  281.             j_max = ceil ((header.y_max - MIN (MIN (vy[0], vy[1]), vy[2]) + yinc2) * idy);
  282.             
  283.             for (j = j_min; j <= j_max; j++) {
  284.                 yp = header.y_max - j * header.y_inc - yinc2;
  285.                 p = j * header.nx + i_min;
  286.                 for (i = i_min; i <= i_max; i++, p++) {
  287.                     xp = header.x_min + i * header.x_inc + xinc2;
  288.                     
  289.                     if (!non_zero_winding (xp, yp, vx, vy, 4)) continue;    /* Outside */
  290.                     
  291.                     grd[p] = a * xp + b * yp + c;
  292.                 }
  293.             }
  294.         }
  295.         dummy[0] = dummy[1] = dummy[2] = dummy[3] = 0;
  296.         if (write_grd (outfile, &header, grd, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  297.             fprintf (stderr, "triangulate: Error writing file %s\n", outfile);
  298.             exit (-1);
  299.         }
  300.                     
  301.     }
  302.     if (multi_segments) {    /* Must find unique edges to output only once */
  303.         n_edge = 3 * np;
  304.         edge = (struct EDGE *) memory (CNULL, n_edge, sizeof (struct EDGE), "triangulate");
  305.         for (i = ij1 = 0, ij2 = 1, ij3 = 2; i < np; i++, ij1 += 3, ij2 += 3, ij3 += 3) {
  306.             edge[ij1].begin = link[ij1];    edge[ij1].end = link[ij2];
  307.             edge[ij2].begin = link[ij2];    edge[ij2].end = link[ij3];
  308.             edge[ij3].begin = link[ij1];    edge[ij3].end = link[ij3];
  309.         }
  310.         
  311.         qsort ((char *)edge, n_edge, sizeof (struct EDGE), compare_edge);
  312.         for (i = 1, j = 0; i < n_edge; i++) {
  313.             if (edge[i].begin != edge[j].begin || edge[i].end != edge[j].end) j++;
  314.             edge[j] = edge[i];
  315.         }
  316.         n_edge = j + 1;
  317.     
  318.         if (gmtdefs.verbose) fprintf (stderr, "triangulate: %d unique triangle edges\n", n_edge);
  319.         
  320.         sprintf (format, "%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format);
  321.         for (i = 0; i < n_edge; i++) {
  322.             printf ("%c Edge %d-%d\n", EOL_flag, edge[i].begin, edge[i].end);
  323.             printf (format, xx[edge[i].begin], yy[edge[i].begin]);
  324.             printf (format, xx[edge[i].end], yy[edge[i].end]);
  325.         }
  326.         free ((char *)edge);
  327.     }
  328.     else if (b_out)
  329.         fwrite ((char *)link, sizeof (int), 3*np, stdout);
  330.     else
  331.         for (i = ij = 0; i < np; i++, ij += 3) printf ("%d\t%d\t%d\n", link[ij], link[ij+1], link[ij+2]);
  332.     
  333.     free ((char *) xx);
  334.     free ((char *) yy);
  335.     free ((char *) link);
  336.     
  337.     if (gmtdefs.verbose) fprintf (stderr, "triangulate: Done!\n");
  338.     
  339.     gmt_end (argc, argv);
  340. }
  341.  
  342. int compare_edge (a, b)
  343. struct EDGE *a, *b; {
  344.     if (a->begin < b->begin)
  345.         return (-1);
  346.     else if (a->begin > b->begin)
  347.         return (1);
  348.     else {
  349.         if (a->end < b->end)
  350.             return (-1);
  351.         else if (a->end > b->end)
  352.             return (1);
  353.         else
  354.             return (0);
  355.     }
  356. }
  357.