home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / pscontour.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-08-09  |  26.2 KB  |  819 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:   @(#)pscontour.c    2.22  09 Aug 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. /* pscontour will read a file of points in the plane, performs the
  8.  * Delaunay triangulation, and contours these triangles.  As an option
  9.  * the user may provide a file with indeces of which vertices constitute
  10.  * the triangles.
  11.  *
  12.  * Author:    Paul Wessel
  13.  * Date:    1-JAN-1995
  14.  * Version:    3.0
  15.  *
  16.  */
  17.  
  18. #include "gmt.h"
  19.  
  20. /* Macros used below */
  21. int delta, sum;
  22. /* Returns plus/minus 1 if node2 is after/before node1 */
  23. #define get_direction(node_1, node_2) (((delta = (node_2) - (node_1)) == 1 || delta == -2) ? 1 : -1)
  24. /* Returns the id of the node common to the two edges */
  25. #define get_node_index(edge_1, edge_2) (((sum = (edge_1) + (edge_2)) == 1) ? 1 : ((sum == 2) ? 0 : 2))
  26.  
  27. struct PSCONTOUR {
  28.     int n_alloc, nl;
  29.     double val;
  30.     double angle;
  31.     char type;
  32.     struct PSLINE *L;
  33. } *cont;
  34.  
  35. struct PSLINE {    /* Beginning and end of straight contour segment */
  36.     double x0, y0;
  37.     double x1, y1;
  38. };
  39.  
  40. struct CHAIN {
  41.     struct PT *begin;
  42.     struct PT *end;
  43.     struct CHAIN *next;
  44. };
  45.  
  46. struct PT {
  47.     double x, y;
  48.     struct PT *next;
  49. };
  50.  
  51. struct LABEL {
  52.     double x, y;
  53.     double angle;
  54.     char label[32];
  55.     struct LABEL *next_label, *prev_label;
  56. } *anchor, *old_label;
  57.     
  58.  
  59. main(argc, argv)
  60. int argc;
  61. char *argv[]; {
  62.     int n, np, nx, i, j, ij, k, kk, k2, k3, way1, way2, node, n1, n2, c, r, g, b;
  63.     int n_alloc, box = 1, section = 0, ix, iy, *ind, *vert, *cind, n_fields, n_contours;
  64.     int add, close, br, bg, bb, label_font_size = 9, bad;
  65.     
  66.     BOOLEAN error = FALSE, image_them = FALSE, draw_contours = FALSE, bin = FALSE, dump = FALSE, clip = TRUE;
  67.     BOOLEAN fix_angle = FALSE, more, t_set = FALSE, single_precision = FALSE, multi_segments = FALSE;
  68.     
  69.     double xx[3], yy[3], zz[3], xout[4], yout[4], *in, west, east, south, north, z_min;
  70.     double *xc, *yc, *zc, *x, *y, *z, *xp, *yp, anot_dist, label_angle = 0.0;
  71.     
  72.     char line[512], *cpt_file = CNULL, *t_file = CNULL, dfile[512], EOL_flag = '>';
  73.     
  74.     FILE *fp = NULL, *fp_d = NULL;
  75.     
  76.     struct PEN pen;
  77.  
  78.     argc = gmt_begin (argc, argv);
  79.     
  80.     gmt_init_pen (&pen, 1);
  81.     anot_dist = 4.0; /* Distance in inches between anotations */
  82.     if (!gmtdefs.measure_unit) anot_dist = 10.0;    /* change to cm */
  83.     dfile[0] = 0;
  84.     br = gmtdefs.page_rgb[0];    /* Default box color is page color */
  85.     bg = gmtdefs.page_rgb[1];
  86.     bb = gmtdefs.page_rgb[2];
  87.     
  88.     /* Check and interpret the command line arguments */
  89.     
  90.     for (i = 1; i < argc; i++) {
  91.         if (argv[i][0] == '-') {
  92.             switch(argv[i][1]) {
  93.         
  94.                 /* Common parameters */
  95.             
  96.                 case 'B':
  97.                 case 'H':
  98.                 case 'J':
  99.                 case 'K':
  100.                 case 'O':
  101.                 case 'P':
  102.                 case 'R':
  103.                 case 'U':
  104.                 case 'V':
  105.                 case 'X':
  106.                 case 'x':
  107.                 case 'Y':
  108.                 case 'y':
  109.                 case 'c':
  110.                 case ':':
  111.                 case '\0':
  112.                     error += get_common_args (argv[i], &west, &east, &south, &north);
  113.                     break;
  114.                 
  115.                 /* Supplemental parameters */
  116.             
  117.                 case 'A':
  118.                     for (j = 2, bad = 0; argv[i][j] && argv[i][j] != 'f'; j++);
  119.                     if (argv[i][j])    { /* Found font size option */
  120.                         label_font_size = atoi (&argv[i][j+1]);
  121.                         if (label_font_size <= 0) bad++;
  122.                     }
  123.                         
  124.                     for (j = 2; argv[i][j] && argv[i][j] != 'a'; j++);
  125.                     if (argv[i][j])    { /* Found fixed angle option */
  126.                         label_angle = atof (&argv[i][j+1]);
  127.                         fix_angle = TRUE;
  128.                         if (label_angle <= -90.0 || label_angle > 180.0) bad++;
  129.                     }
  130.                         
  131.                     for (j = 2; argv[i][j] && argv[i][j] != '/'; j++);
  132.                     if (argv[i][j] && gmt_getrgb (&argv[i][j+1], &br, &bg, &bb)) bad++;
  133.                     if (strchr (argv[i], 'o')) box = 2;
  134.                     if (bad) {
  135.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -A option.  Correct syntax:\n", gmt_program);
  136.                         fprintf (stderr, "\t-A[f<fontsize>][a<fixedangle>][/<r/g/b>][o]\n");
  137.                         error += bad;
  138.                     }
  139.                     break;
  140.  
  141.                 case 'b':
  142.                     single_precision = !(argv[i][strlen(argv[i])-1] == 'd');
  143.                     bin = TRUE;
  144.                     break;
  145.                 case 'C':
  146.                     cpt_file = &argv[i][2];
  147.                     break;
  148.                 case 'D':
  149.                     dump = TRUE;
  150.                     strcpy (dfile, &argv[i][2]);
  151.                     break;
  152.                 case 'E':
  153.                     sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
  154.                     break;
  155.                 case 'F':
  156.                     if (gmt_getrgb (&argv[i][2], &gmtdefs.basemap_frame_rgb[0], &gmtdefs.basemap_frame_rgb[1], &gmtdefs.basemap_frame_rgb[2])) {
  157.                         gmt_pen_syntax ('F');
  158.                         error++;
  159.                     }
  160.                 case 'G':
  161.                     anot_dist = atof (&argv[i][2]);
  162.                     break;
  163.                 case 'I':
  164.                     image_them = TRUE;
  165.                     break;
  166.                 case 'M':    /* with -D, create one multiple line segments */
  167.                     multi_segments = 2;
  168.                     if (argv[i][2]) EOL_flag = argv[i][2];
  169.                     break;
  170.                 case 'N':
  171.                     clip = FALSE;
  172.                     break;
  173.                 case 'T':
  174.                     t_file = &argv[i][2];
  175.                     t_set = TRUE;
  176.                     break;
  177.                 case 'W':
  178.                     if (gmt_getpen (&argv[i][2], &pen)) {
  179.                         gmt_pen_syntax ('W');
  180.                         error++;
  181.                     }
  182.                     else
  183.                         draw_contours = TRUE;
  184.                     break;
  185.                     
  186.                 /* Options not recognized */
  187.                         
  188.                 default:
  189.                     error = TRUE;
  190.                     gmt_default_error (argv[i][1]);
  191.                     break;
  192.             }
  193.         }
  194.         else
  195.             fp = fopen (argv[i], "r");
  196.     }
  197.     
  198.     if (gmt_quick || argc == 1) {    /* Display usage */
  199.         fprintf (stderr,"pscontour %s - Contour xyz-data by triangulation\n\n", GMT_VERSION);
  200.         fprintf(stderr,"usage: pscontour <xyzfile> -C<cpt_file> -J<params> -R<west>/<east>/<south>/<north>\n");
  201.         fprintf (stderr, "    [-A[f<fontsize>][/r/g/b][a<angle>][o]] [-B<tickinfo>] [-D<dumpfile>] [-E<az>/<el>]\n");
  202.         fprintf (stderr, "    [-F<r>/<g>/<b>] [-H] [-I] [-K] [-M[<flag>]] [-N] [-O] [-P] [-T<indexfile>]\n");
  203.         fprintf (stderr, "    [-U] [-V] [-W<pen>] [-X<x_shift>] [-Y<y_shift>] [-c<ncopies>] [-:] [-b[d]]\n\n");
  204.         
  205.         if (gmt_quick) exit(-1);
  206.         
  207.         fprintf(stderr,"    -C Color palette table\n");
  208.         explain_option ('j');
  209.         explain_option ('R');
  210.         fprintf (stderr, "\n\tOPTIONS:\n");
  211.         fprintf (stderr, "    -A Annotation format information.\n");
  212.         fprintf (stderr, "       Append f followed by desired font size in points [Default is 9].\n");
  213.         fprintf (stderr, "       Append /r/g/b to change color of text box [Default is %d/%d/%d]\n", br, bg, bb);
  214.         fprintf (stderr, "       Append o to draw outline of text box [Default is no outline]\n");
  215.         fprintf (stderr, "       Append a<angle> to force anotations at this fixed angle [Default follows contour]\n");
  216.         explain_option ('b');
  217.         fprintf (stderr, "    -D to Dump contour lines to individual files (but see -M)\n");
  218.         fprintf (stderr, "    -E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n");
  219.         fprintf (stderr, "    -F Set color used for Frame and anotation [%d/%d/%d]\n",
  220.             gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  221.         fprintf (stderr, "    -G Gap between anotations in %s [Default = %lg]\n", gmt_unit_names[gmtdefs.measure_unit], anot_dist);
  222.         explain_option ('H');
  223.         fprintf (stderr, "    -I Color triangles using the cpt file\n");
  224.         explain_option ('K');
  225.         fprintf (stderr, "      -M Used with -D.   Create a single multiple segment file where contours are separated by a record\n");
  226.         fprintf (stderr, "         whose first character is <flag> ['>'].  This header also has the contour level value\n");
  227.         fprintf(stderr,"    -N do NOT clip contours/image at the border [Default clips]\n");
  228.         explain_option ('O');
  229.         explain_option ('P');
  230.         fprintf(stderr,"    -T file with triplets of point indeces for each triangle\n");
  231.         fprintf(stderr,"       [Default performs the Delauney triangulation on xyz-data]\n");
  232.         explain_option ('U');
  233.         explain_option ('V');
  234.         fprintf (stderr, "    -W selects contouring and sets contour pen attributes\n");
  235.         explain_option ('X');
  236.         explain_option ('c');
  237.         explain_option (':');
  238.         fprintf (stderr, "\t-b means binary input for xyz and node ids (ints)) [ASCII]\n");
  239.         fprintf (stderr, "\t   Append d for double precision [Default is single precision].\n");        explain_option ('.');
  240.         explain_option ('.');
  241.         exit(-1);
  242.     }
  243.  
  244.     /* Check that the options selected are mutually consistant */
  245.     
  246.     if (!project_info.region_supplied) {
  247.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", gmt_program);
  248.         error++;
  249.     }
  250.     if (!(draw_contours || image_them)) {
  251.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify one of -W or -I\n", gmt_program);
  252.         error++;
  253.     }
  254.     if (!cpt_file) {
  255.         fprintf (stderr, "%s: GMT SYNTAX ERROR -C option:  Must specify a color palette table\n", gmt_program);
  256.         error++;
  257.     }
  258.     if (t_set && !t_file) {
  259.         fprintf (stderr, "%s: GMT SYNTAX ERROR -T option:  Must specify an index file\n", gmt_program);
  260.         error++;
  261.     }
  262.     if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) {
  263.         fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  Enter azimuth in 0-360 range, elevation in 0-90 range\n", gmt_program);
  264.         error++;
  265.     }
  266.     if (bin && gmtdefs.io_header) {
  267.         fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", gmt_program);
  268.         error++;
  269.     }
  270.     
  271.     if (error) exit (-1);
  272.  
  273.     if (bin) gmt_input = (single_precision) ? bin_float_input : bin_double_input;
  274.  
  275.     if (t_set) {
  276.         if ((fp_d = fopen (t_file, "r")) == NULL) {
  277.         fprintf (stderr, "pscontour: Could not open index file %s\n", t_file);
  278.             exit (-1);
  279.         }
  280.     }
  281.     
  282.     if (dump && dfile[0] == 0) {
  283.         fprintf (stderr, "pscontour: contours wil be written to file contour\n");
  284.         strcpy (dfile,"contour");
  285.     }
  286.  
  287.     read_cpt (cpt_file);
  288.     if (image_them && gmt_continuous) {
  289.         fprintf (stderr, "pscontour: -I option requires constant color between contours!\n");
  290.         exit (-1);
  291.     }
  292.     
  293.     map_setup (west, east, south, north);
  294.     
  295.     ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
  296.         gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies, gmtdefs.dpi,
  297.         gmtdefs.measure_unit, gmtdefs.paper_width, gmtdefs.page_rgb, gmt_epsinfo (argv[0]));
  298.         
  299.     echo_command (argc, argv);
  300.     if (gmtdefs.unix_time) timestamp (argc, argv);
  301.     if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0);
  302.     
  303.         if (clip) map_clip_on (-1, -1, -1, 3);
  304.  
  305.     if (fp == NULL) fp = stdin;
  306.     
  307.     n_alloc = GMT_CHUNK;
  308.     x = (double *) memory (CNULL, n_alloc, sizeof (double), "pscontour");
  309.     y = (double *) memory (CNULL, n_alloc, sizeof (double), "pscontour");
  310.     z = (double *) memory (CNULL, n_alloc, sizeof (double), "pscontour");
  311.     
  312.     if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (line, 512, fp);
  313.     
  314.     ix = (gmtdefs.xy_toggle) ? 1 : 0;        iy = 1 - ix;              /* Set up which columns have x and y */
  315.  
  316.     n = 0;
  317.     more = ((n_fields = gmt_input (fp,  3, &in)) == 3);
  318.         
  319.     while (more) {
  320.         if (bad_float (in[2])) continue;    /* Skip when z = NaN */
  321.         
  322.         x[n] = in[ix];
  323.         y[n] = in[iy];
  324.         z[n] = in[2];
  325.         n++;
  326.         
  327.         if (n == n_alloc) {
  328.             n_alloc += GMT_CHUNK;
  329.             x = (double *) memory ((char *)x, n_alloc, sizeof (double), "pscontour");
  330.             y = (double *) memory ((char *)y, n_alloc, sizeof (double), "pscontour");
  331.             z = (double *) memory ((char *)z, n_alloc, sizeof (double), "pscontour");
  332.         }
  333.             
  334.         more = ((n_fields = gmt_input (fp,  3, &in)) == 3);
  335.     }
  336.     if (fp != stdin) fclose (fp);
  337.     if (n_fields == -1) n_fields = 0;    /* Ok to get EOF */
  338.     if (n_fields%3) {    /* Got garbage or multiple segment subheader */
  339.         fprintf (stderr, "%s:  Cannot read X Y Z line %d, aborts\n", gmt_program, n);
  340.         exit (-1);
  341.     }
  342.     
  343.     x = (double *) memory ((char *)x, n, sizeof (double), "pscontour");
  344.     y = (double *) memory ((char *)y, n, sizeof (double), "pscontour");
  345.     z = (double *) memory ((char *)z, n, sizeof (double), "pscontour");
  346.     
  347.     /* Map transform */
  348.     
  349.     for (i = 0; i < n; i++) geo_to_xy (x[i], y[i], &x[i], &y[i]);
  350.     
  351.     if (fp_d) {
  352.         n_alloc = 3 * GMT_CHUNK;
  353.         ind = (int *) memory (CNULL, n_alloc, sizeof (int), "pscontour");
  354.  
  355.         ij = np = 0;
  356.         
  357.         if (bin)    /* Binary input */
  358.             more = (fread ((char *)ind, sizeof (int), 3, fp_d) == 3);
  359.         else        /* ascii input */
  360.             more = (fgets (line, 512, fp_d) != CNULL);
  361.             
  362.         while (more) {
  363.             if (!bin && sscanf (line, "%d %d %d", &ind[ij], &ind[ij+1], &ind[ij+2]) != 3) continue;
  364.             ij += 3;
  365.             np++;
  366.             if (ij == n_alloc) {
  367.                 n_alloc += (3 * GMT_CHUNK);
  368.                 ind = (int *) memory ((char *)ind, n_alloc, sizeof (int), "pscontour");
  369.             }
  370.             if (bin)    /* Binary input */
  371.                 more = (fread ((char *)&ind[ij], sizeof (int), 3, fp_d) == 3);
  372.             else        /* ascii input */
  373.                 more = (fgets (line, 512, fp_d) != CNULL);
  374.         }
  375.         ind = (int *) memory ((char *)ind, ij, sizeof (int), "pscontour");
  376.         fclose (fp_d);
  377.     }
  378.     else    /* Do Delauney triangulation */
  379.             
  380.         np = delaunay (x, y, n, &ind);
  381.     
  382.     /* Get PSCONTOUR structs */
  383.     
  384.     n_contours = gmt_n_colors + 1;
  385.     cont = (struct PSCONTOUR *) memory (CNULL, n_contours, sizeof (struct PSCONTOUR), "pscontour");
  386.     
  387.     for (i = 0; i < gmt_n_colors; i++) {
  388.         cont[i].val = gmt_lut[i].z_low;
  389.         cont[i].type = (gmt_lut[i].anot) ? 'A' : 'C';
  390.         cont[i].angle = (fix_angle) ? label_angle : gmt_NaN;
  391.     }
  392.     cont[gmt_n_colors].val = gmt_lut[gmt_n_colors-1].z_high;
  393.     cont[gmt_n_colors].type = (gmt_lut[gmt_n_colors-1].anot & 2) ? 'A' : 'C';
  394.     cont[gmt_n_colors].angle = (fix_angle) ? label_angle : gmt_NaN;
  395.     for (i = 0; i < n_contours; i++) {
  396.         cont[i].n_alloc = GMT_SMALL_CHUNK;
  397.         cont[i].L = (struct PSLINE *)  memory (CNULL, GMT_SMALL_CHUNK, sizeof (struct PSLINE), "pscontour");
  398.     }
  399.     
  400.     ps_setpaint (pen.r, pen.g, pen.b);
  401.     ps_setline (pen.width);
  402.     if (pen.texture) ps_setdash (pen.texture, pen.offset);
  403.  
  404.     for (i = ij = 0; i < np; i++, ij += 3) {    /* For all triangles */
  405.     
  406.         k = ij;
  407.         xx[0] = x[ind[k]];    yy[0] = y[ind[k]];    zz[0] = z[ind[k++]];
  408.         xx[1] = x[ind[k]];    yy[1] = y[ind[k]];    zz[1] = z[ind[k++]];
  409.         xx[2] = x[ind[k]];    yy[2] = y[ind[k]];    zz[2] = z[ind[k]];
  410.         
  411.         nx = get_triangle_crossings (x, y, z, &ind[ij], &xc, &yc, &zc, &vert, &cind);
  412.         
  413.         if (image_them) {    /* Must color the triangle slices according to cpt file */
  414.         
  415.             /* First paint background color for entire triangle */
  416.             
  417.             z_min = MIN (zz[0], MIN (zz[1], zz[2]));
  418.             get_rgb24 (z_min, &r, &g, &b);
  419.             if (project_info.three_D) {
  420.                 for (k = 0; k < 3; k++) xy_do_z_to_xy (xx[k], yy[k], project_info.z_level, &xout[k], &yout[k]);
  421.                 ps_patch (xout, yout, 3, r, g, b, FALSE);
  422.             }
  423.             else
  424.                 ps_patch (xx, yy, 3, r, g, b, FALSE);
  425.                 
  426.             /* Then loop over contours and paint the part that is higher up in the color scheme */
  427.             
  428.             for (k = k2 = 0, k3 = 1; k < nx; k++, k2 += 2, k3 += 2) {
  429.                 xout[0] = xc[k2];    yout[0] = yc[k2];
  430.                 xout[1] = xc[k3];    yout[1] = yc[k3];
  431.                 node = get_node_index (vert[k2], vert[k3]);
  432.                 
  433.                 /* We either add one vertix and paint that triangle or we must add the other
  434.                    two vertices and paint a polygon.  The if-test will tell us which one */
  435.                    
  436.                 if (zz[node] > zc[k2]) {    /* Add this single vertex to path */
  437.                     xout[2] = xx[node];
  438.                     yout[2] = yy[node];
  439.                     get_rgb24 (zc[k2], &r, &g, &b);
  440.                     if (project_info.three_D)
  441.                         for (kk = 0; kk < 3; kk++) xy_do_z_to_xy (xout[kk], yout[kk], project_info.z_level, &xout[kk], &yout[kk]);
  442.                     ps_patch (xout, yout, 3, r, g, b, FALSE);
  443.                 }
  444.                 else {    /* Must add the other two vertices */
  445.                     n1 = (node + 1) % 3;
  446.                     n2 = (node + 2) % 3;
  447.                     way1 = get_direction (vert[k2], vert[k3]);
  448.                     way2 = get_direction (n1, n2);
  449.                     if (way1 * way2 < 0) i_swap (n1, n2);
  450.                     xout[2] = xx[n1];    yout[2] = yy[n1];
  451.                     xout[3] = xx[n2];    yout[3] = yy[n2];
  452.                     get_rgb24 (zc[k2], &r, &g, &b);
  453.                     if (project_info.three_D)
  454.                         for (kk = 0; kk < 4; kk++) xy_do_z_to_xy (xout[kk], yout[kk], project_info.z_level, &xout[kk], &yout[kk]);
  455.                     ps_patch (xout, yout, 4, r, g, b, FALSE);
  456.                 }
  457.             }
  458.         }
  459.         
  460.         if (draw_contours && nx > 0) {    /* Save contour lines for later */
  461.         
  462.             if (project_info.three_D)
  463.                 for (k = 0; k < 2*nx; k++) xy_do_z_to_xy (xc[k], yc[k], project_info.z_level, &xc[k], &yc[k]);
  464.                 
  465.             for (k = k2 = 0; k < nx; k++) {
  466.                 c = cind[k];
  467.                 n = cont[c].nl;
  468.                 cont[c].L[n].x0 = xc[k2];
  469.                 cont[c].L[n].y0 = yc[k2++];
  470.                 cont[c].L[n].x1 = xc[k2];
  471.                 cont[c].L[n].y1 = yc[k2++];
  472.                 n++;
  473.                 if (n >= cont[c].n_alloc) {
  474.                     cont[c].n_alloc += GMT_SMALL_CHUNK;
  475.                     cont[c].L = (struct PSLINE *)  memory ((char *)cont[c].L, cont[c].n_alloc, sizeof (struct PSLINE), "pscontour");
  476.                 }
  477.                 cont[c].nl = n;
  478.             }
  479.                 
  480.             /* for (k = 0; k < nx; k++) ps_line (&xc[2*k], &yc[2*k], 2, 3, FALSE, FALSE); */
  481.         }
  482.         
  483.         if (nx > 0) {
  484.             free ((char *)xc);
  485.             free ((char *)yc);
  486.             free ((char *)zc);
  487.             free ((char *)vert);
  488.             free ((char *)cind);
  489.         }
  490.     }
  491.     
  492.     /* Draw contours */
  493.     
  494.     if (draw_contours) {
  495.         
  496.         struct CHAIN *head_c, *last_c, *this_c;
  497.         struct PT *p, *q;
  498.         
  499.         anchor = old_label = (struct LABEL *) memory (CNULL, 1, sizeof (struct LABEL), "pscontour");
  500.  
  501.         for (c = 0; c < n_contours; c++) {
  502.         
  503.             if (cont[c].nl == 0) {
  504.                 free ((char *)cont[c].L);
  505.                 continue;
  506.             }
  507.             
  508.             head_c = last_c = (struct CHAIN *) memory (CNULL, 1, sizeof (struct CHAIN), "pscontour");
  509.  
  510.             while (cont[c].nl) {
  511.                 this_c = last_c->next = (struct CHAIN *) memory (CNULL, 1, sizeof (struct CHAIN), "pscontour");
  512.                 k = 0;
  513.                 this_c->begin = (struct PT *) memory (CNULL, 1, sizeof (struct PT), "pscontour");
  514.                 this_c->end = (struct PT *) memory (CNULL, 1, sizeof (struct PT), "pscontour");
  515.                 this_c->begin->x = cont[c].L[k].x0;
  516.                 this_c->begin->y = cont[c].L[k].y0;
  517.                 this_c->end->x = cont[c].L[k].x1;
  518.                 this_c->end->y = cont[c].L[k].y1;
  519.                 this_c->begin->next = this_c->end;
  520.                 cont[c].nl--;
  521.                 cont[c].L[k] = cont[c].L[cont[c].nl];
  522.                 while (k < cont[c].nl) {
  523.                     add = 0;
  524.                     if (fabs(cont[c].L[k].x0 - this_c->begin->x) < SMALL && fabs(cont[c].L[k].y0 - this_c->begin->y) < SMALL) {
  525.                         p = (struct PT *) memory (CNULL, 1, sizeof (struct PT), "pscontour");
  526.                         p->x = cont[c].L[k].x1;
  527.                         p->y = cont[c].L[k].y1;
  528.                         p->next = this_c->begin;
  529.                         add = -1;
  530.                     }
  531.                     else if (fabs(cont[c].L[k].x1 - this_c->begin->x) < SMALL && fabs(cont[c].L[k].y1 - this_c->begin->y) < SMALL) {
  532.                         p = (struct PT *) memory (CNULL, 1, sizeof (struct PT), "pscontour");
  533.                         p->x = cont[c].L[k].x0;
  534.                         p->y = cont[c].L[k].y0;
  535.                         p->next = this_c->begin;
  536.                         add = -1;
  537.                     }
  538.                     else if (fabs(cont[c].L[k].x0 - this_c->end->x) < SMALL && fabs(cont[c].L[k].y0 - this_c->end->y) < SMALL) {
  539.                         p = (struct PT *) memory (CNULL, 1, sizeof (struct PT), "pscontour");
  540.                         p->x = cont[c].L[k].x1;
  541.                         p->y = cont[c].L[k].y1;
  542.                         this_c->end->next = p;
  543.                         add = 1;
  544.                     }
  545.                     else if (fabs(cont[c].L[k].x1 - this_c->end->x) < SMALL && fabs(cont[c].L[k].y1 - this_c->end->y) < SMALL) {
  546.                         p = (struct PT *) memory (CNULL, 1, sizeof (struct PT), "pscontour");
  547.                         p->x = cont[c].L[k].x0;
  548.                         p->y = cont[c].L[k].y0;
  549.                         this_c->end->next = p;
  550.                         add = 1;
  551.                     }
  552.                     if (add) {    /* Got one */
  553.                         if (add == -1)
  554.                             this_c->begin = p;
  555.                         else if (add == 1)
  556.                             this_c->end = p;
  557.                         cont[c].nl--;
  558.                         cont[c].L[k] = cont[c].L[cont[c].nl];
  559.                         k = 0;
  560.                     }
  561.                     else
  562.                         k++;
  563.                 }
  564.                 last_c = this_c;
  565.             }
  566.             free ((char *)cont[c].L);
  567.             
  568.             this_c = head_c->next;
  569.             while (this_c) {
  570.                 xp = (double *) memory (CNULL, GMT_SMALL_CHUNK, sizeof (double), "pscontour");
  571.                 yp = (double *) memory (CNULL, GMT_SMALL_CHUNK, sizeof (double), "pscontour");
  572.                 n_alloc = GMT_SMALL_CHUNK;
  573.                 p = this_c->begin;
  574.                 n = 0;
  575.                 while (p) {
  576.                     xp[n] = p->x;
  577.                     yp[n++] = p->y;
  578.                     q = p;
  579.                     p = p->next;
  580.                     free ((char *)q);
  581.                     if (n == n_alloc) {
  582.                         n_alloc += GMT_SMALL_CHUNK;
  583.                         xp = (double *) memory ((char *)xp, n_alloc, sizeof (double), "pscontour");
  584.                         yp = (double *) memory ((char *)yp, n_alloc, sizeof (double), "pscontour");
  585.                     }
  586.                 }
  587.                 last_c = this_c;
  588.                 this_c = this_c->next;
  589.                 free ((char *)last_c);
  590.                 
  591.                 close = (xp[0] == xp[n-1] && yp[0] == yp[n-1]);
  592.                 draw_contour (xp, yp, n, cont[c].val, cont[c].type, cont[c].angle, close, anot_dist);
  593.                 if (dump) dump_contour (xp, yp, n, cont[c].val, section++, close, dfile, &multi_segments, EOL_flag);
  594.  
  595.                 free ((char *)xp);
  596.                 free ((char *)yp);
  597.             }
  598.                     
  599.         }
  600.         plot_labels (label_font_size, box, br, bg, bb);
  601.     }
  602.                             
  603.     if (pen.texture) ps_setdash (CNULL, 0);
  604.     
  605.         if (clip) map_clip_off ();
  606.  
  607.     if (frame_info.plot) {
  608.         ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  609.         map_basemap ();
  610.     }
  611.     
  612.     ps_setpaint (0, 0, 0);
  613.         
  614.     if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);
  615.  
  616.     ps_plotend (gmtdefs.last_page);
  617.     
  618.     free ((char *)x);
  619.     free ((char *)y);
  620.     free ((char *)z);
  621.     free ((char *)ind);
  622.     
  623.     gmt_end (argc, argv);
  624. }
  625.  
  626.  
  627. int get_triangle_crossings (x, y, z, ind, xc, yc, zc, v, cindex)
  628. double x[], y[], z[];
  629. int ind[];
  630. double *xc[], *yc[], *zc[];
  631. int *v[], *cindex[]; {
  632.     int i, j, k, k2, i1, nx, n_alloc, *vout, *cind;
  633.     double xx[3], yy[3], zz[3], zmin, zmax, dz, frac, *xout, *yout, *zout;
  634.     
  635.     xx[0] = x[ind[0]];    yy[0] = y[ind[0]];    zz[0] = z[ind[0]];
  636.     xx[1] = x[ind[1]];    yy[1] = y[ind[1]];    zz[1] = z[ind[1]];
  637.     xx[2] = x[ind[2]];    yy[2] = y[ind[2]];    zz[2] = z[ind[2]];
  638.  
  639.     zmin = MIN (zz[0], MIN (zz[1], zz[2]));
  640.     zmax = MAX (zz[0], MAX (zz[1], zz[2]));
  641.     
  642.     i = 0;    j = gmt_n_colors - 1;
  643.     while (gmt_lut[i].z_low <= zmin && i < gmt_n_colors) i++;
  644.     while (gmt_lut[j].z_high >= zmax && j > 0) j--;
  645.     
  646.     nx = j - i + 2;
  647.     
  648.     if (nx == 0) return (0);
  649.     
  650.     n_alloc = 2 * nx;
  651.     xout = (double *) memory (CNULL, n_alloc, sizeof (double), "pscontour");
  652.     yout = (double *) memory (CNULL, n_alloc, sizeof (double), "pscontour");
  653.     zout = (double *) memory (CNULL, n_alloc, sizeof (double), "pscontour");
  654.     vout = (int *) memory (CNULL, n_alloc, sizeof (int), "pscontour");
  655.     cind = (int *) memory (CNULL, nx, sizeof (int), "pscontour");
  656.     
  657.     k = k2 = 0;
  658.     while (i <= j) {
  659.         zout[k2] = zout[k2+1] = gmt_lut[i].z_low;
  660.         cind[k++] = i;
  661.         k2 += 2;
  662.         i++;
  663.     }
  664.     zout[k2] = zout[k2+1] = gmt_lut[j].z_high;
  665.     cind[k] = j + 1;
  666.         
  667.     for (k = k2 = j = 0; k < nx; k++, k2 += 2) {
  668.         for (i = 0; i < 3; i++) {
  669.             i1 = (i == 2) ? 0 : i + 1;
  670.             if ((zout[k2] >= zz[i] && zout[k2] < zz[i1]) || (zout[k2] <= zz[i] && zout[k2] > zz[i1])) {
  671.                 dz = zz[i1] - zz[i];
  672.                 if (dz == 0.0) {    /* Contour goes along ende */
  673.                     xout[j] = xx[i];    yout[j] = yy[i];
  674.                 }
  675.                 else {
  676.                     frac = (zout[k2] - zz[i]) / dz;
  677.                     xout[j] = xx[i] + frac * (xx[i1] - xx[i]);
  678.                     yout[j] = yy[i] + frac * (yy[i1] - yy[i]);
  679.                 }
  680.                 vout[j++] = i;
  681.             }
  682.         }
  683.         if (j%2) j--;    /* Contour went through a single vertice only, skip this */
  684.     }
  685.     
  686.     *xc = xout;
  687.     *yc = yout;
  688.     *zc = zout;
  689.     *v = vout;
  690.     *cindex = cind;
  691.     
  692.     return (j/2);
  693. }
  694.  
  695. int draw_contour (xx, yy, nn, cval, ctype, cangle, closed, gap)
  696. double *xx, *yy, cval, cangle, gap;
  697. char ctype;
  698. int nn, closed; {
  699.     int i, j;
  700.     double dist, angle, dx, dy, width, one = 1.0;
  701.     char label[100], format[50];
  702.     struct LABEL *new_label;
  703.     
  704.     if (nn < 2) return;
  705.         
  706.     sprintf (format, "%lg contour\0", cval);
  707.     ps_comment (format);
  708.     
  709.     ps_line (xx, yy, nn, 3, closed, TRUE);
  710.     
  711.     if (ctype == 'A' || ctype == 'a') {    /* Annotated contours */
  712.         if (!gmtdefs.measure_unit) one = 2.5;
  713.         get_format (cval, format);
  714.         sprintf (label, format, cval);
  715.         dist = (closed) ? gap - one : 0.0;    /* Label closed contours longer than 1 inch */
  716.         for (i = 1; i < nn; i++) {
  717.  
  718.             dx = xx[i] - xx[i-1];
  719.             if (fabs (dx) > (width = gmt_half_map_width (yy[i-1]))) {
  720.                 width *= 2.0;
  721.                 dx = copysign (width - fabs (dx), -dx);
  722.                 if (xx[i] < width)
  723.                     xx[i-1] -= width;
  724.                 else
  725.                     xx[i-1] += width;
  726.             }
  727.             dy = yy[i] - yy[i-1];
  728.             dist += hypot (dx, dy);
  729.             if (dist > gap) {    /* Time for label */
  730.                 new_label = (struct LABEL *) memory (CNULL, 1, sizeof (struct LABEL), "pscontour");
  731.                 new_label->x = 0.5 * (xx[i-1] + xx[i]);
  732.                 new_label->y = 0.5 * (yy[i-1] + yy[i]);
  733.                 strcpy (new_label->label, label);
  734.                 if (bad_float (cangle)) {    /* Must calculate label angle */
  735.                     angle = d_atan2 (dy, dx) * R2D;
  736.                     if (angle < 0.0) angle += 360.0;
  737.                     if (angle > 90.0 && angle < 270) angle -= 180.0;
  738.                 }
  739.                 else
  740.                     angle = cangle;
  741.                 new_label->angle = angle;
  742.                 new_label->prev_label = old_label;
  743.                 dist = 0.0;
  744.                 old_label->next_label = new_label;
  745.                 old_label = new_label;
  746.             }
  747.         }
  748.     }
  749. }
  750.  
  751. int plot_labels (size, box, r, g, b)
  752. int size, box, r, g, b; {
  753.     /* box = 1: white box only, box = 2: white box + draw outline */
  754.     double dx, dy;
  755.     struct LABEL *this, *old;
  756.     
  757.     dx = 0.5 * gmtdefs.anot_offset;
  758.     dy = 0.05 * gmtdefs.anot_offset;
  759.     box--;
  760.     ps_comment ("Contour annotations:");
  761.     
  762.     ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  763.     
  764.     for (old = anchor; old->next_label; old = old->next_label) {    /* First draw boxes if not 3-D*/
  765.         this = old->next_label;
  766.         gmt_textbox3d (this->x, this->y, project_info.z_level, size, gmtdefs.anot_font, this->label, this->angle, 6, box, dx, dy, r, g, b);
  767.     }
  768.     for (old = anchor; old->next_label; old = old->next_label) { /* Then labels */
  769.         this = old->next_label;
  770.         gmt_text3d (this->x, this->y, project_info.z_level, size, gmtdefs.anot_font, this->label, this->angle, 6, 0);
  771.     }
  772.         
  773.     ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  774.     
  775.     this = anchor;
  776.     while (this) {    /* Free memory */
  777.         old = this;
  778.         this = old->next_label;
  779.         free ((char *)old);
  780.     }
  781. }
  782.  
  783. int dump_contour (xx, yy, nn, cval, id, interior, file, m_segment, flag)
  784. double *xx, *yy, cval;
  785. char *file, flag;
  786. int nn, id;
  787. BOOLEAN interior, *m_segment; {
  788.     int i;
  789.     double lon, lat;
  790.     char fname[512], format[80];
  791.     FILE *fp;
  792.     
  793.     if (nn < 2) return;
  794.     
  795.     sprintf (format, "%s\t%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
  796.     if (*m_segment) {
  797.         if (*m_segment == 2) {    /* Must create file the first time around */
  798.             fp = fopen (file, "w");
  799.             *m_segment = TRUE;
  800.         }
  801.         else    /* Later we append to it */
  802.             fp = fopen (file, "a+");
  803.         fprintf (fp, "%c %lg contour\n", flag, cval);
  804.     }
  805.     else {
  806.         if (interior)
  807.             sprintf (fname, "%s_%lg_%d_i.xyz\0", file, cval, id);
  808.         else
  809.             sprintf (fname, "%s_%lg_%d.xyz\0", file, cval, id);
  810.         fp = fopen (fname, "w");
  811.     }
  812.     for (i = 0; i < nn; i++) {
  813.         xy_to_geo (&lon, &lat, xx[i], yy[i]);
  814.         fprintf (fp, format, lon, lat, cval);
  815.     }
  816.     fclose (fp);
  817. }
  818.  
  819.