/*-------------------------------------------------------------------- * $Id: pscontour.c,v 1.4.4.5 2002/02/27 17:58:55 pwessel Exp $ * * Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith * See COPYING file for copying and redistribution conditions. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; version 2 of the License. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * Contact info: gmt.soest.hawaii.edu *--------------------------------------------------------------------*/ /* pscontour will read a file of points in the plane, performs the * Delaunay triangulation, and contours these triangles. As an option * the user may provide a file with indeces of which vertices constitute * the triangles. * * Author: Paul Wessel * Date: 1-JAN-1995 * Version: 3.4.1 */ #include "gmt.h" /* Macros used below */ int delta, sum; /* Returns plus/minus 1 if node2 is after/before node1 */ #define get_direction(node_1, node_2) (((delta = (node_2) - (node_1)) == 1 || delta == -2) ? 1 : -1) /* Returns the id of the node common to the two edges */ #define get_node_index(edge_1, edge_2) (((sum = (edge_1) + (edge_2)) == 1) ? 1 : ((sum == 2) ? 0 : 2)) #define get_other_node(node1, node2) (((sum = (node1 + node2)) == 3) ? 0 : ((sum == 2) ? 1 : 2)) /* The other node needed */ struct PSCONTOUR { int n_alloc, nl; double val; double angle; char type; struct PSLINE *L; } *cont; struct PSLINE { /* Beginning and end of straight contour segment */ double x0, y0; double x1, y1; }; struct CHAIN { struct PT *begin; struct PT *end; struct CHAIN *next; }; struct PT { double x, y; struct PT *next; }; struct LABEL { double x, y; double angle; char label[64]; struct LABEL *next_label, *prev_label; } *anchor, *old_label; int get_triangle_crossings(double *x, double *y, double *z, int *ind, double **xc, double **yc, double **zc, int **v, int **cindex); void draw_contour(double *xx, double *yy, int nn, double cval, char ctype, double cangle, int closed, double gap), plot_labels(int size, int box, int rgb[]); void dump_contour(double *xx, double *yy, int nn, double cval, int id, BOOLEAN interior, char *file); void paint_it (double x[], double y[], int n, double z); main(int argc, char **argv) { int n, np, nx, i, j, ij, k, k2, k3, node1, node2, c, rgb[3]; int n_alloc, box = 1, section = 0, *ind, *vert, *cind, n_fields, n_contours; int add, close, b_rgb[3], label_font_size = 9, bad, n_expected_fields, low, high, last_entry, last_exit; BOOLEAN error = FALSE, image_them = FALSE, draw_contours = FALSE, dump = FALSE, clip = TRUE; BOOLEAN fix_angle = FALSE, more, t_set = FALSE, draw_mesh = FALSE, use_cpt_colors = FALSE, use_cpt_anot = TRUE; double xx[3], yy[3], zz[3], xout[5], yout[5], *in, west, east, south, north; double *xc, *yc, *zc, *x, *y, *z, *xp, *yp, anot_dist, label_angle = 0.0, current_contour = -DBL_MAX; char line[BUFSIZ], *cpt_file = CNULL, *t_file = CNULL, dfile[BUFSIZ]; FILE *fp = NULL, *fp_d = NULL; struct GMT_PEN pen, lpen; argc = GMT_begin (argc, argv); GMT_init_pen (&pen, GMT_PENWIDTH); GMT_init_pen (&lpen, GMT_PENWIDTH); anot_dist = 4.0; /* Distance in inches between anotations */ if (gmtdefs.measure_unit == GMT_CM) anot_dist = 10.0 / 2.54; /* change to 10 cm */ dfile[0] = 0; for (i = 0; i < 3; i++) b_rgb[i] = gmtdefs.page_rgb[i]; /* Default box color is page color */ /* Check and interpret the command line arguments */ for (i = 1; i < argc; i++) { if (argv[i][0] == '-') { switch(argv[i][1]) { /* Common parameters */ case 'B': case 'H': case 'J': case 'K': case 'O': case 'P': case 'R': case 'U': case 'V': case 'X': case 'x': case 'Y': case 'y': case 'c': case ':': case '\0': error += GMT_get_common_args (argv[i], &west, &east, &south, &north); break; /* Supplemental parameters */ case 'A': if (argv[i][2] == '-') use_cpt_anot = FALSE; for (j = 2, bad = 0; argv[i][j] && argv[i][j] != 'f'; j++); if (argv[i][j]) { /* Found font size option */ label_font_size = atoi (&argv[i][j+1]); if (label_font_size <= 0) bad++; } for (j = 2; argv[i][j] && argv[i][j] != 'a'; j++); if (argv[i][j]) { /* Found fixed angle option */ label_angle = atof (&argv[i][j+1]); fix_angle = TRUE; if (label_angle < -90.0 || label_angle > 180.0) bad++; } for (j = 2; argv[i][j] && argv[i][j] != '/'; j++); if (argv[i][j] && GMT_getrgb (&argv[i][j+1], b_rgb)) bad++; if (strchr (argv[i], 'o')) box = 2; if (bad) { fprintf (stderr, "%s: GMT SYNTAX ERROR -A option. Correct syntax:\n", GMT_program); fprintf (stderr, "\t-A[-][f][a][/][o]\n"); error += bad; } break; case 'b': error += GMT_io_selection (&argv[i][2]); break; case 'C': cpt_file = &argv[i][2]; break; case 'D': dump = TRUE; strcpy (dfile, &argv[i][2]); break; case 'E': sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation); break; case 'G': anot_dist = GMT_convert_units (&argv[i][2], GMT_INCH); break; case 'I': image_them = TRUE; break; case 'L': if (GMT_getpen (&argv[i][2], &lpen)) { GMT_pen_syntax ('L'); error++; } else draw_mesh = TRUE; break; case 'M': /* with -D, create one multiple line segments */ GMT_multisegment (&argv[i][2]); GMT_io.multi_segments = 2; break; case 'N': clip = FALSE; break; case 'T': t_file = &argv[i][2]; t_set = TRUE; break; case 'W': k = 2; if (argv[i][k] == '+') use_cpt_colors = TRUE, k++; if (argv[i][k] && GMT_getpen (&argv[i][k], &pen)) { GMT_pen_syntax ('W'); error++; } else draw_contours = TRUE; break; /* Options not recognized */ default: error = TRUE; GMT_default_error (argv[i][1]); break; } } else fp = GMT_fopen (argv[i], GMT_io.r_mode); } if (GMT_quick || argc == 1) { /* Display usage */ fprintf (stderr,"pscontour %s - Contour xyz-data by triangulation\n\n", GMT_VERSION); fprintf (stderr,"usage: pscontour -C -J -R///\n"); fprintf (stderr, "\t[-A[-][f][/r/g/b][a][o]] [-B] [-D] [-E/] [-G]\n"); fprintf (stderr, "\t[-H[]] [-I] [-K] [-L] [-M[]] [-N] [-O] [-P] [-T]\n"); fprintf (stderr, "\t[-U] [-V] [-W[+]] [-X] [-Y] [-c] [-:] [-bi[s][]] [-bo[s][]]\n\n"); if (GMT_quick) exit (EXIT_FAILURE); fprintf (stderr,"\t-C Color palette table\n"); GMT_explain_option ('j'); GMT_explain_option ('R'); fprintf (stderr, "\n\tOPTIONS:\n"); fprintf (stderr, "\t-A Annotation format information.\n"); fprintf (stderr, "\t Append - to disable all contour anotations\n"); fprintf (stderr, "\t Append f followed by desired font size in points [Default is 9].\n"); fprintf (stderr, "\t Append /r/g/b to change color of text box [Default is %d/%d/%d]\n", b_rgb[0], b_rgb[1], b_rgb[2]); fprintf (stderr, "\t Append o to draw outline of text box [Default is no outline]\n"); fprintf (stderr, "\t Append a to force anotations at this fixed angle [Default follows contour]\n"); GMT_explain_option ('b'); fprintf (stderr, "\t-D to Dump contour lines to individual files (but see -M)\n"); fprintf (stderr, "\t-E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n"); fprintf (stderr, "\t-G Gap between anotations in %s [Default = %lg]\n", GMT_unit_names[gmtdefs.measure_unit], anot_dist); GMT_explain_option ('H'); fprintf (stderr, "\t-I Color triangles using the cpt file\n"); GMT_explain_option ('K'); fprintf (stderr, "\t-L draws the triangular mesh with the specified pen\n"); fprintf (stderr, "\t-M Used with -D. Create a single multiple segment file where contours are separated by a record\n"); fprintf (stderr, "\t whose first character is ['>']. This header also has the contour level value\n"); fprintf (stderr,"\t-N do NOT clip contours/image at the border [Default clips]\n"); GMT_explain_option ('O'); GMT_explain_option ('P'); fprintf (stderr,"\t-T file with triplets of point indeces for each triangle\n"); fprintf (stderr,"\t [Default performs the Delauney triangulation on xyz-data]\n"); GMT_explain_option ('U'); GMT_explain_option ('V'); fprintf (stderr, "\t-W selects contouring and sets contour pen attributes\n"); fprintf (stderr, "\t Use + to draw colored contours based on the cpt file\n"); GMT_explain_option ('X'); GMT_explain_option ('c'); GMT_explain_option (':'); GMT_explain_option ('i'); GMT_explain_option ('n'); fprintf (stderr,"\t Default is 3 input columns.\n"); GMT_explain_option ('o'); GMT_explain_option ('.'); exit (EXIT_FAILURE); } /* Check that the options selected are mutually consistant */ if (!project_info.region_supplied) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify -R option\n", GMT_program); error++; } if (!(draw_contours || image_them)) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify one of -W or -I\n", GMT_program); error++; } if (!cpt_file) { fprintf (stderr, "%s: GMT SYNTAX ERROR -C option: Must specify a color palette table\n", GMT_program); error++; } if (t_set && !t_file) { fprintf (stderr, "%s: GMT SYNTAX ERROR -T option: Must specify an index file\n", GMT_program); error++; } if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -E option: Enter azimuth in 0-360 range, elevation in 0-90 range\n", GMT_program); error++; } if (GMT_io.binary[0] && gmtdefs.io_header) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data cannot have header -H\n", GMT_program); error++; } if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = 3; if (GMT_io.binary[0] && GMT_io.ncol[0] < 3) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data (-bi) must have at least 3 columns\n", GMT_program); error++; } if (error) exit (EXIT_FAILURE); GMT_put_history (argc, argv); /* Update .gmtcommands */ if (GMT_io.binary[0] && gmtdefs.verbose) { char *type[2] = {"double", "single"}; fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[0], type[GMT_io.single_precision[0]]); } if (GMT_io.binary[1] && GMT_io.multi_segments) { fprintf (stderr, "%s: GMT Warning. -M for output ignored with -D, -bo\n", GMT_program); error++; } if (t_set) { if ((fp_d = fopen (t_file, GMT_io.r_mode)) == NULL) { fprintf (stderr, "%s: Could not open index file %s\n", GMT_program, t_file); exit (EXIT_FAILURE); } } if (dump && dfile[0] == 0) { fprintf (stderr, "%s: contours will be written to file contour\n", GMT_program); strcpy (dfile,"contour"); } GMT_read_cpt (cpt_file); if (image_them && GMT_continuous) { fprintf (stderr, "%s: -I option requires constant color between contours!\n", GMT_program); exit (EXIT_FAILURE); } GMT_map_setup (west, east, south, north); ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin, gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies, gmtdefs.dpi, GMT_INCH, gmtdefs.paper_width, gmtdefs.page_rgb, GMT_epsinfo (argv[0])); GMT_echo_command (argc, argv); if (gmtdefs.unix_time) GMT_timestamp (argc, argv); if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0); if (clip) GMT_map_clip_on (GMT_no_rgb, 3); if (fp == NULL) { fp = GMT_stdin; #ifdef SET_IO_MODE GMT_setmode (0); #endif } n_alloc = GMT_CHUNK; x = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program); y = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program); z = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program); if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (line, BUFSIZ, fp); n = 0; n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : 3; n_fields = GMT_input (fp, &n_expected_fields, &in); while (! (GMT_io.status & GMT_IO_EOF)) { /* Not yet EOF */ if (GMT_io.status & GMT_IO_MISMATCH) { fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields, n_expected_fields, n); exit (EXIT_FAILURE); } if (!GMT_is_dnan (in[2])) { /* Unless z = NaN */ x[n] = in[0]; y[n] = in[1]; z[n] = in[2]; n++; if (n == n_alloc) { n_alloc += GMT_CHUNK; x = (double *) GMT_memory ((void *)x, (size_t)n_alloc, sizeof (double), GMT_program); y = (double *) GMT_memory ((void *)y, (size_t)n_alloc, sizeof (double), GMT_program); z = (double *) GMT_memory ((void *)z, (size_t)n_alloc, sizeof (double), GMT_program); } } n_fields = GMT_input (fp, &n_expected_fields, &in); } if (fp != GMT_stdin) GMT_fclose (fp); x = (double *) GMT_memory ((void *)x, (size_t)n, sizeof (double), GMT_program); y = (double *) GMT_memory ((void *)y, (size_t)n, sizeof (double), GMT_program); z = (double *) GMT_memory ((void *)z, (size_t)n, sizeof (double), GMT_program); /* Map transform */ for (i = 0; i < n; i++) GMT_geo_to_xy (x[i], y[i], &x[i], &y[i]); if (fp_d) { n_alloc = 3 * GMT_CHUNK; ind = (int *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (int), GMT_program); ij = np = 0; if (GMT_io.binary[0]) /* Binary input */ more = (fread ((void *)ind, sizeof (int), (size_t)3, fp_d) == 3); else /* ascii input */ more = (fgets (line, BUFSIZ, fp_d) != CNULL); while (more) { if (!GMT_io.binary[0] && sscanf (line, "%d %d %d", &ind[ij], &ind[ij+1], &ind[ij+2]) != 3) continue; ij += 3; np++; if (ij == n_alloc) { n_alloc += (3 * GMT_CHUNK); ind = (int *) GMT_memory ((void *)ind, (size_t)n_alloc, sizeof (int), GMT_program); } if (GMT_io.binary[0]) /* Binary input */ more = (fread ((void *)&ind[ij], sizeof (int), (size_t)3, fp_d) == 3); else /* ascii input */ more = (fgets (line, BUFSIZ, fp_d) != CNULL); } ind = (int *) GMT_memory ((void *)ind, (size_t)ij, sizeof (int), GMT_program); fclose (fp_d); } else /* Do Delauney triangulation */ np = GMT_delaunay (x, y, n, &ind); if (draw_mesh) { /* Draw triangular mesh */ GMT_setpen (&lpen); for (i = k = 0; i < np; i++) { /* For all triangles */ xx[0] = x[ind[k]]; yy[0] = y[ind[k++]]; xx[1] = x[ind[k]]; yy[1] = y[ind[k++]]; xx[2] = x[ind[k]]; yy[2] = y[ind[k++]]; ps_line (xx, yy, 3, 3, 1, FALSE); } } /* Get PSCONTOUR structs */ n_contours = GMT_n_colors + 1; cont = (struct PSCONTOUR *) GMT_memory (VNULL, (size_t)n_contours, sizeof (struct PSCONTOUR), GMT_program); for (i = 0; i < GMT_n_colors; i++) { cont[i].val = GMT_lut[i].z_low; cont[i].type = (GMT_lut[i].anot && use_cpt_anot) ? 'A' : 'C'; cont[i].angle = (fix_angle) ? label_angle : GMT_d_NaN; } cont[GMT_n_colors].val = GMT_lut[GMT_n_colors-1].z_high; cont[GMT_n_colors].type = ((GMT_lut[GMT_n_colors-1].anot & 2) && use_cpt_anot) ? 'A' : 'C'; cont[GMT_n_colors].angle = (fix_angle) ? label_angle : GMT_d_NaN; for (i = 0; i < n_contours; i++) { cont[i].n_alloc = GMT_SMALL_CHUNK; cont[i].L = (struct PSLINE *) GMT_memory (VNULL, (size_t)GMT_SMALL_CHUNK, sizeof (struct PSLINE), GMT_program); } GMT_setpen (&pen); for (i = ij = 0; i < np; i++, ij += 3) { /* For all triangles */ k = ij; xx[0] = x[ind[k]]; yy[0] = y[ind[k]]; zz[0] = z[ind[k++]]; xx[1] = x[ind[k]]; yy[1] = y[ind[k]]; zz[1] = z[ind[k++]]; xx[2] = x[ind[k]]; yy[2] = y[ind[k]]; zz[2] = z[ind[k]]; nx = get_triangle_crossings (x, y, z, &ind[ij], &xc, &yc, &zc, &vert, &cind); if (image_them) { /* Must color the triangle slices according to cpt file */ if (nx == 0) { /* No contours go through - easy, but must check for NaNs */ int kzz; double zzz; for (k = kzz = 0, zzz = 0.0; k < 3; k++) { if (GMT_is_dnan (zz[k])) continue; zzz += zz[k]; kzz++; } if (kzz) paint_it (xx, yy, 3, zzz / kzz); } else { /* Must paint all those slices separately */ /* Find vertecies with the lowest and highest values */ for (k = 1, low = high = 0; k < 3; k++) { if (zz[k] < zz[low]) low = k; if (zz[k] > zz[high]) high = k; } /* Paint the piece delimited by the low node and the first contour */ xout[0] = xx[low]; yout[0] = yy[low]; node1 = get_node_index (vert[0], vert[1]); /* Find single vertex opposing this contour segment */ if (node1 == low) { /* Contour and low node make up a triangle */ xout[1] = xc[0]; yout[1] = yc[0]; xout[2] = xc[1]; yout[2] = yc[1]; n = 3; } else { /* Need twe other two vertices to form a 4-sided polygon */ node2 = get_other_node (node1, low); /* The other node needed */ xout[1] = xx[node2]; yout[1] = yy[node2]; if (low == vert[0] || node2 == vert[1]) { /* Add segment in opposite order */ xout[2] = xc[1]; yout[2] = yc[1]; xout[3] = xc[0]; yout[3] = yc[0]; } else { /* Add in regular order */ xout[2] = xc[0]; yout[2] = yc[0]; xout[3] = xc[1]; yout[3] = yc[1]; } n = 4; } paint_it (xout, yout, n, 0.5 * (zz[low] + zc[1])); /* z is contour value */ /* Then loop over contours and paint the part between contours */ for (k = 1, k2 = 2, k3 = 3; k < nx; k++, k2 += 2, k3 += 2) { xout[0] = xc[k2-2]; yout[0] = yc[k2-2]; xout[1] = xc[k3-2]; yout[1] = yc[k3-2]; n = 2; last_entry = vert[k2-2]; last_exit = vert[k3-2]; if (last_exit == vert[k2]) { xout[n] = xc[k2]; yout[n] = yc[k2]; n++; xout[n] = xc[k3]; yout[n] = yc[k3]; n++; if (vert[k3] != last_entry) { /* Need to add an intervening corner */ node1 = get_node_index (last_entry, vert[k3]); /* Find corner id */ xout[n] = xx[node1]; yout[n] = yy[node1]; n++; } } else if (last_exit == vert[k3]) { xout[n] = xc[k3]; yout[n] = yc[k3]; n++; xout[n] = xc[k2]; yout[n] = yc[k2]; n++; if (vert[k2] != last_entry) { /* Need to add an intervening corner */ node1 = get_node_index (last_entry, vert[k2]); /* Find corner id */ xout[n] = xx[node1]; yout[n] = yy[node1]; n++; } } else if (last_entry == vert[k2]) { node1 = get_node_index (last_exit, vert[k3]); /* Find corner id */ xout[n] = xx[node1]; yout[n] = yy[node1]; n++; xout[n] = xc[k3]; yout[n] = yc[k3]; n++; xout[n] = xc[k2]; yout[n] = yc[k2]; n++; } else { node1 = get_node_index (last_exit, vert[k2]); /* Find corner id */ xout[n] = xx[node1]; yout[n] = yy[node1]; n++; xout[n] = xc[k2]; yout[n] = yc[k2]; n++; xout[n] = xc[k3]; yout[n] = yc[k3]; n++; } paint_it (xout, yout, n, 0.5 * (zc[k2]+zc[k2-2])); } /* Add the last piece between last contour and high node */ k2 -= 2; k3 -= 2; xout[0] = xx[high]; yout[0] = yy[high]; node1 = get_node_index (vert[k2], vert[k3]); /* Find corner id */ if (node1 == high) { /* Cut off a triangular piece */ xout[1] = xc[k2]; yout[1] = yc[k2]; xout[2] = xc[k3]; yout[2] = yc[k3]; n = 3; } else { /* Need a 4-sided polygon */ node2 = get_other_node (node1, high); /* The other node needed */ xout[1] = xx[node2]; yout[1] = yy[node2]; if (high == vert[0] || node2 == vert[1]) { /* On same side, start here */ xout[2] = xc[k3]; yout[2] = yc[k3]; xout[3] = xc[k2]; yout[3] = yc[k2]; } else { /* On same side, start here */ xout[2] = xc[k2]; yout[2] = yc[k2]; xout[3] = xc[k3]; yout[3] = yc[k3]; } n = 4; } paint_it (xout, yout, n, 0.5 * (zz[high] + zc[k2])); /* z is contour value */ } } if (draw_contours && nx > 0) { /* Save contour lines for later */ if (project_info.three_D) for (k = 0; k < 2*nx; k++) GMT_xy_do_z_to_xy (xc[k], yc[k], project_info.z_level, &xc[k], &yc[k]); for (k = k2 = 0; k < nx; k++) { c = cind[k]; n = cont[c].nl; cont[c].L[n].x0 = xc[k2]; cont[c].L[n].y0 = yc[k2++]; cont[c].L[n].x1 = xc[k2]; cont[c].L[n].y1 = yc[k2++]; n++; if (n >= cont[c].n_alloc) { cont[c].n_alloc += GMT_SMALL_CHUNK; cont[c].L = (struct PSLINE *) GMT_memory ((void *)cont[c].L, (size_t)cont[c].n_alloc, sizeof (struct PSLINE), GMT_program); } cont[c].nl = n; } /* for (k = 0; k < nx; k++) ps_line (&xc[2*k], &yc[2*k], 2, 3, FALSE, FALSE); */ } if (nx > 0) { GMT_free ((void *)xc); GMT_free ((void *)yc); GMT_free ((void *)zc); GMT_free ((void *)vert); GMT_free ((void *)cind); } } /* Draw contours */ if (draw_contours) { struct CHAIN *head_c, *last_c, *this_c; struct PT *p, *q; anchor = old_label = (struct LABEL *) GMT_memory (VNULL, (size_t)1, sizeof (struct LABEL), GMT_program); for (c = 0; c < n_contours; c++) { if (cont[c].nl == 0) { GMT_free ((void *)cont[c].L); continue; } head_c = last_c = (struct CHAIN *) GMT_memory (VNULL, (size_t)1, sizeof (struct CHAIN), GMT_program); while (cont[c].nl) { this_c = last_c->next = (struct CHAIN *) GMT_memory (VNULL, (size_t)1, sizeof (struct CHAIN), GMT_program); k = 0; this_c->begin = (struct PT *) GMT_memory (VNULL, (size_t)1, sizeof (struct PT), GMT_program); this_c->end = (struct PT *) GMT_memory (VNULL, (size_t)1, sizeof (struct PT), GMT_program); this_c->begin->x = cont[c].L[k].x0; this_c->begin->y = cont[c].L[k].y0; this_c->end->x = cont[c].L[k].x1; this_c->end->y = cont[c].L[k].y1; this_c->begin->next = this_c->end; cont[c].nl--; cont[c].L[k] = cont[c].L[cont[c].nl]; while (k < cont[c].nl) { add = 0; if (fabs(cont[c].L[k].x0 - this_c->begin->x) < SMALL && fabs(cont[c].L[k].y0 - this_c->begin->y) < SMALL) { p = (struct PT *) GMT_memory (VNULL, (size_t)1, sizeof (struct PT), GMT_program); p->x = cont[c].L[k].x1; p->y = cont[c].L[k].y1; p->next = this_c->begin; add = -1; } else if (fabs(cont[c].L[k].x1 - this_c->begin->x) < SMALL && fabs(cont[c].L[k].y1 - this_c->begin->y) < SMALL) { p = (struct PT *) GMT_memory (VNULL, (size_t)1, sizeof (struct PT), GMT_program); p->x = cont[c].L[k].x0; p->y = cont[c].L[k].y0; p->next = this_c->begin; add = -1; } else if (fabs(cont[c].L[k].x0 - this_c->end->x) < SMALL && fabs(cont[c].L[k].y0 - this_c->end->y) < SMALL) { p = (struct PT *) GMT_memory (VNULL, (size_t)1, sizeof (struct PT), GMT_program); p->x = cont[c].L[k].x1; p->y = cont[c].L[k].y1; this_c->end->next = p; add = 1; } else if (fabs(cont[c].L[k].x1 - this_c->end->x) < SMALL && fabs(cont[c].L[k].y1 - this_c->end->y) < SMALL) { p = (struct PT *) GMT_memory (VNULL, (size_t)1, sizeof (struct PT), GMT_program); p->x = cont[c].L[k].x0; p->y = cont[c].L[k].y0; this_c->end->next = p; add = 1; } if (add) { /* Got one */ if (add == -1) this_c->begin = p; else if (add == 1) this_c->end = p; cont[c].nl--; cont[c].L[k] = cont[c].L[cont[c].nl]; k = 0; } else k++; } last_c = this_c; } GMT_free ((void *)cont[c].L); this_c = head_c->next; while (this_c) { xp = (double *) GMT_memory (VNULL, (size_t)GMT_SMALL_CHUNK, sizeof (double), GMT_program); yp = (double *) GMT_memory (VNULL, (size_t)GMT_SMALL_CHUNK, sizeof (double), GMT_program); n_alloc = GMT_SMALL_CHUNK; p = this_c->begin; n = 0; while (p) { xp[n] = p->x; yp[n++] = p->y; q = p; p = p->next; GMT_free ((void *)q); if (n == n_alloc) { n_alloc += GMT_SMALL_CHUNK; xp = (double *) GMT_memory ((void *)xp, (size_t)n_alloc, sizeof (double), GMT_program); yp = (double *) GMT_memory ((void *)yp, (size_t)n_alloc, sizeof (double), GMT_program); } } last_c = this_c; this_c = this_c->next; GMT_free ((void *)last_c); close = (xp[0] == xp[n-1] && yp[0] == yp[n-1]); if (use_cpt_colors) { if (current_contour != cont[c].val) { GMT_get_rgb24 (cont[c].val, rgb); ps_setpaint (rgb); current_contour = cont[c].val; } } draw_contour (xp, yp, n, cont[c].val, cont[c].type, cont[c].angle, close, anot_dist); if (dump) dump_contour (xp, yp, n, cont[c].val, section++, close, dfile); GMT_free ((void *)xp); GMT_free ((void *)yp); } } plot_labels (label_font_size, box, b_rgb); } if (clip) GMT_map_clip_off (); if (frame_info.plot) GMT_map_basemap (); if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0); ps_plotend (gmtdefs.last_page); GMT_free ((void *)x); GMT_free ((void *)y); GMT_free ((void *)z); GMT_free ((void *)ind); GMT_end (argc, argv); } int get_triangle_crossings (double *x, double *y, double *z, int *ind, double **xc, double **yc, double **zc, int **v, int **cindex) { /* This routine finds all the contour crossings for this triangle. Each contour consists of * linesegments made up of two points, with coordinates xc, yc, and contour level zc. */ int i, j, k, k2, i1, nx, n_alloc, *vout, *cind; double xx[3], yy[3], zz[3], zmin, zmax, dz, frac, *xout, *yout, *zout; xx[0] = x[ind[0]]; yy[0] = y[ind[0]]; zz[0] = z[ind[0]]; xx[1] = x[ind[1]]; yy[1] = y[ind[1]]; zz[1] = z[ind[1]]; xx[2] = x[ind[2]]; yy[2] = y[ind[2]]; zz[2] = z[ind[2]]; if (GMT_is_dnan (zz[0]) || GMT_is_dnan (zz[1]) || GMT_is_dnan (zz[2])) return (0); /* Cannot have crossings if NaNs are present */ zmin = MIN (zz[0], MIN (zz[1], zz[2])); /* Min z vertex */ zmax = MAX (zz[0], MAX (zz[1], zz[2])); /* Max z vertex */ i = 0; j = GMT_n_colors - 1; while (GMT_lut[i].z_low <= zmin && i < GMT_n_colors) i++; while (GMT_lut[j].z_high >= zmax && j > 0) j--; nx = j - i + 2; /* Total number of contours */ if (nx <= 0) return (0); n_alloc = 2 * nx; xout = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program); yout = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program); zout = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program); vout = (int *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (int), GMT_program); cind = (int *) GMT_memory (VNULL, (size_t)nx, sizeof (int), GMT_program); /* Fill out array zout which holds the nx contour levels */ k = k2 = 0; while (i <= j) { zout[k2] = zout[k2+1] = GMT_lut[i].z_low; cind[k++] = i; k2 += 2; i++; } zout[k2] = zout[k2+1] = GMT_lut[j].z_high; cind[k] = j + 1; /* Loop over the contour levels and determine the line segments */ for (k = k2 = j = 0; k < nx; k++, k2 += 2) { for (i = 0; i < 3; i++) { /* Try each side in turn 0-1, 1-2, 2-0 */ i1 = (i == 2) ? 0 : i + 1; if ((zout[k2] >= zz[i] && zout[k2] < zz[i1]) || (zout[k2] <= zz[i] && zout[k2] > zz[i1])) { dz = zz[i1] - zz[i]; if (dz == 0.0) { /* Contour goes along ende */ xout[j] = xx[i]; yout[j] = yy[i]; } else { frac = (zout[k2] - zz[i]) / dz; xout[j] = xx[i] + frac * (xx[i1] - xx[i]); yout[j] = yy[i] + frac * (yy[i1] - yy[i]); } vout[j++] = i; /* Keep track of the side number */ } } if (j%2) j--; /* Contour went through a single vertice only, skip this */ } *xc = xout; *yc = yout; *zc = zout; *v = vout; *cindex = cind; return (j/2); } void draw_contour (double *xx, double *yy, int nn, double cval, char ctype, double cangle, int closed, double gap) { int i; double dist, angle, dx, dy, width; char label[128], format[64]; struct LABEL *new_label; if (nn < 2) return; sprintf (format, "%lg contour", cval); ps_comment (format); ps_line (xx, yy, nn, 3, closed, TRUE); if (ctype == 'A' || ctype == 'a') { /* Annotated contours */ GMT_get_format (cval, CNULL, format); sprintf (label, format, cval); dist = (closed) ? gap - 1.0: 0.0; /* Label closed contours longer than 1 inch */ for (i = 1; i < nn; i++) { dx = xx[i] - xx[i-1]; if (fabs (dx) > (width = GMT_half_map_width (yy[i-1]))) { width *= 2.0; dx = copysign (width - fabs (dx), -dx); if (xx[i] < width) xx[i-1] -= width; else xx[i-1] += width; } dy = yy[i] - yy[i-1]; dist += hypot (dx, dy); if (dist > gap) { /* Time for label */ new_label = (struct LABEL *) GMT_memory (VNULL, (size_t)1, sizeof (struct LABEL), GMT_program); new_label->x = 0.5 * (xx[i-1] + xx[i]); new_label->y = 0.5 * (yy[i-1] + yy[i]); strcpy (new_label->label, label); if (GMT_is_dnan (cangle)) { /* Must calculate label angle */ angle = d_atan2 (dy, dx) * R2D; if (angle < 0.0) angle += 360.0; if (angle > 90.0 && angle < 270) angle -= 180.0; } else angle = cangle; new_label->angle = angle; new_label->prev_label = old_label; dist = 0.0; old_label->next_label = new_label; old_label = new_label; } } } } void plot_labels (int size, int box, int rgb[]) { /* box = 1: white box only, box = 2: white box + draw outline */ double dx, dy; struct LABEL *this, *old; dx = 0.5 * gmtdefs.anot_offset; dy = 0.05 * gmtdefs.anot_offset; box--; ps_comment ("Contour annotations:"); ps_setpaint (gmtdefs.basemap_frame_rgb); for (old = anchor; old->next_label; old = old->next_label) { /* First draw boxes if not 3-D*/ this = old->next_label; GMT_textbox3d (this->x, this->y, project_info.z_level, size, gmtdefs.anot_font, this->label, this->angle, 6, box, dx, dy, rgb); } for (old = anchor; old->next_label; old = old->next_label) { /* Then labels */ this = old->next_label; GMT_text3d (this->x, this->y, project_info.z_level, size, gmtdefs.anot_font, this->label, this->angle, 6, 0); } ps_setpaint (gmtdefs.background_rgb); this = anchor; while (this) { /* Free memory */ old = this; this = old->next_label; GMT_free ((void *)old); } } void dump_contour (double *xx, double *yy, int nn, double cval, int id, BOOLEAN interior, char *file) { int i; double out[3]; char fname[BUFSIZ], format[80], suffix[4]; FILE *fp; if (nn < 2) return; out[2] = cval; (GMT_io.binary[1]) ? strcpy (suffix, "b") : strcpy (suffix, "xyz"); sprintf (format, "%s\t%s\t%s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format); if (!GMT_io.binary[1] && GMT_io.multi_segments) { if (GMT_io.multi_segments == 2) { /* Must create file the first time around */ fp = GMT_fopen (file, "w"); GMT_io.multi_segments = TRUE; } else /* Later we append to it */ fp = GMT_fopen (file, "a+"); sprintf (GMT_io.segment_header, "%c %lg contour", GMT_io.EOF_flag, cval); GMT_write_segmentheader (fp, 3); } else { if (interior) sprintf (fname, "%s_%lg_%d_i.%s", file, cval, id, suffix); else sprintf (fname, "%s_%lg_%d.%s", file, cval, id, suffix); fp = GMT_fopen (fname, GMT_io.w_mode); } for (i = 0; i < nn; i++) { GMT_xy_to_geo (&out[0], &out[1], xx[i], yy[i]); GMT_output (fp, 3, out); } GMT_fclose (fp); } void paint_it (double x[], double y[], int n, double z) { int index, k, rgb[3]; struct GMT_FILL *f; if (n < 3) return; /* Need at least 3 points to make a polygon */ index = GMT_get_rgb24 (z, rgb); if (GMT_cpt_skip) return; /* Skip this z-slice */ /* Now we must paint, with colors or patterns */ if (project_info.three_D) for (k = 0; k < n; k++) GMT_xy_do_z_to_xy (x[k], y[k], project_info.z_level, &x[k], &y[k]); if ((index >= 0 && (f = GMT_lut[index].fill)) || (index < 0 && (f = GMT_bfn.fill[index+3]))) { GMT_fill (x, y, n, f, FALSE); /* Contours drawn separately later if desired */ } else ps_patch (x, y, n, rgb, FALSE); /* Contours drawn separately later if desired */ }