/*-------------------------------------------------------------------- * $Id: grdview.c,v 1.5.4.2 2002/02/20 18:40:19 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 *--------------------------------------------------------------------*/ /* * grdview will read a topofile and produce a 3-D perspective plot * of the surface z = f(x,y) using PostScript. The surface can * be represented as: * 1) A Mesh plot * 2) A shaded (or colored) surface w/wo contourlines and w/wo * illumination by artificial sun(s). * * grdview calls contours to find the line segments that make up the * contour lines. This allows the user to specify that the contours * should be smoothed before plotting. This will make the resulting * image smoother, especially if nx and ny are relatively small. * * Author: Paul Wessel * Date: 8-DEC-1988 * Modified: 5-JUL-1994 New version 3.0 * Modified: 15-JUL-1998 New version 3.1 * Modified: 15-MAR-1999 New version 3.2 * Modified: 07-OCT-1999 3.3.2 - Bug in image fixed, plus added correct * image mask clipping * Modified: 07-FEB-2000 3.3.4 - Now works if red = -1 (skips that z-slice) * Understands patterns in cpt files (-Qs only) * Version: 3.4.1 */ #include "gmt.h" /* Declarations needed for binning of smooth contours */ struct BIN { struct CONT *first_cont; } *binij; struct CONT { struct POINT *first_point; struct CONT *next_cont; double value; }; struct POINT { double x, y; struct POINT *next_point; }; /* Global variables */ float *grd, *intensity, *topo; int *edge; int bin_inc[4], ij_inc[4]; double x_inc[4], y_inc[4], *x, *y, *z, *v, *xx, *yy; char *c_method[2] = { "colorimage", "colortiles", }; void grdview_init_setup(struct GRD_HEADER *header, float *topo, int two, BOOLEAN draw_plane, double plane_level); int pixel_inside(int ip, int jp, int *ix, int *iy, int bin), quick_idist(int x1, int y1, int x2, int y2); int get_side (double x, double y, double x_left, double y_bottom, double xinc, double yinc, double dx2, double dy2); void copy_points_fw (double x[], double y[], double z[], double v[], double xcont[], double ycont[], double zcont[], double vcont[], int cont, int *n); void copy_points_bw (double x[], double y[], double z[], double v[], double xcont[], double ycont[], double zcont[], double vcont[], int cont, int *n); double get_z_ave (double v[], double next_up, int n); void add_node (double x[], double y[], double z[], double v[], int *k, int node, double X_vert[], double Y_vert[], float topo[], float grd[], int ij, int bin); void paint_it (double x[], double y[], int n, double z, BOOLEAN intens, double intensity); double get_intensity (float *intensity, int k, int nx); struct CONT *get_cont_struct (int bin, double value); struct POINT *get_point (double x, double y); BOOLEAN monochrome = FALSE; main (int argc, char **argv) { int i, j, ij, n_edges, k, k1, n, max, i_bin, j_bin, i_bin_old, j_bin_old; int side, way, nm, nm2, nx_f, ny_f, off, nx, ny; int rgb_facade[3], bin, two, mx, my, sw, se, nw, ne, id, n4, nk; int c, rgb[3], i_start, i_stop, j_start, j_stop, i_inc, j_inc; int smooth = 1, tiling = 0, dpi_i = 100; BOOLEAN mesh = FALSE, draw_plane = FALSE, facade = FALSE, get_contours, bad, set_z = FALSE; BOOLEAN error = FALSE, outline = FALSE, surface = FALSE, pen_not_set, bilinear = FALSE; BOOLEAN first, begin, draw_contours = FALSE, intens = FALSE, drape = FALSE; BOOLEAN saddle, image = FALSE, subset = FALSE, no_nans, q_set = FALSE; double cval, x_left, x_right, y_top, y_bottom, small, z_ave, this_intensity, *xval, *yval; double plane_level = 0.0, dx2, dy2, take_out; double west = 0.0, east = 0.0, south = 0.0, north = 0.0, new_z_level = 0.0, del_off; double data_west, data_east, data_south, data_north, delx, dely, z_val, next_up, xmesh[4], ymesh[4]; double x_pixel_size, y_pixel_size, *x_imask, *y_imask; char *topofile, *intensfile, *drapefile, *cpt_file; struct CONT *start_cont, *this_cont, *last_cont; struct POINT *this_point, *last_point; struct GRD_HEADER header, t_head, i_head, d_head; struct GMT_PEN pen[2]; struct GMT_FILL fill; struct GMT_EDGEINFO edgeinfo; for (i = 0; i < 3; i++) rgb_facade[i] = 0; cpt_file = topofile = intensfile = drapefile = CNULL; argc = GMT_begin (argc, argv); GMT_3D_mode = 1; /* Only do background axis first; do foreground at end */ GMT_init_pen (&pen[0], 3.0 * GMT_PENWIDTH); /* Contour pen */ GMT_init_pen (&pen[1], GMT_PENWIDTH); /* Mesh pen */ GMT_init_fill (&fill, 255 - gmtdefs.basemap_frame_rgb[0], 255 - gmtdefs.basemap_frame_rgb[0], 255 - gmtdefs.basemap_frame_rgb[0]); for (i = 1; i < argc; i++) { if (argv[i][0] == '-') { switch (argv[i][1]) { /* Common parameters */ case 'B': 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 '\0': error += GMT_get_common_args (argv[i], &west, &east, &south, &north); break; /* Supplemental parameters */ case 'C': cpt_file = &argv[i][2]; break; case 'E': sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation); break; case 'G': drapefile = &argv[i][2]; drape = TRUE; break; case 'I': intensfile = &argv[i][2]; intens = TRUE; break; case 'L': if (argv[i][2]) error += GMT_boundcond_parse (&edgeinfo, &argv[i][2]); else bilinear = TRUE; break; case 'N': if (argv[i][2]) { n = sscanf (&argv[i][2], "%lf/%d/%d/%d", &plane_level, &rgb_facade[0], &rgb_facade[1], &rgb_facade[2]); draw_plane = TRUE; if (n == 4) facade = TRUE; if ((n > 1 && n != 4) || (n == 4 && GMT_check_rgb (rgb_facade))) { fprintf (stderr, "%s: GMT SYNTAX ERROR option -N: RGB components must all be in 0-255 range\n", GMT_program); error = TRUE; } } else { fprintf (stderr, "%s: GMT SYNTAX ERROR option -N: Must specify plane level\n", GMT_program); error = TRUE; } break; case 'Q': q_set = TRUE; switch (argv[i][2]) { case 'm': /* Mesh plot */ mesh = TRUE; if (argv[i][3] == '/' && GMT_getfill (&argv[i][4], &fill)) { fprintf (stderr, "%s: GMT SYNTAX ERROR -Qm option: RGB components must all be in 0-255 range\n", GMT_program); error = TRUE; } break; case 's': /* Color wo/ contours */ surface = TRUE; if (argv[i][3] == 'm') outline = TRUE; break; case 'i': /* image w/ clipmask */ case 'I': /* Backward compatibility, gives -Qi */ image = TRUE; if (argv[i][3] && isdigit ((int)argv[i][3])) dpi_i = atoi (&argv[i][3]); break; default: fprintf (stderr, "%s: GMT SYNTAX ERROR: Unrecognized qualifier (%c) for option -%c\n", GMT_program, argv[i][2], argv[i][1]); error = TRUE; break; } monochrome = (argv[i][strlen(argv[i])-1] == 'g'); break; case 'S': smooth = atoi (&argv[i][2]); break; case 'T': tiling = 1; if (argv[i][2] == 's') tiling = 2; break; case 'W': /* Contour and/or mesh pens */ j = (argv[i][2] == 'm' || argv[i][2] == 'c') ? 3 : 2; if (j == 2) { /* Set both */ if (GMT_getpen (&argv[i][2], &pen[0])) { GMT_pen_syntax ('W'); error++; } else pen[1] = pen[0]; draw_contours = TRUE; } else { id = (argv[i][2] == 'm') ? 1 : 0; if (GMT_getpen (&argv[i][j], &pen[id])) { GMT_pen_syntax ('W'); error++; } if (id == 0) draw_contours = TRUE; } break; case 'Z': if (argv[i][2]) { new_z_level = atof (&argv[i][2]); set_z = TRUE; } break; /* Illegal options */ default: error = TRUE; GMT_default_error (argv[i][1]); break; } } else topofile = argv[i]; } if (!(mesh || surface || image || tiling)) mesh = TRUE; /* Default is mesh plot */ if (argc == 1 || GMT_quick) { fprintf (stderr, "grdview %s - Plot topofiles in 3-D\n\n", GMT_VERSION); fprintf (stderr, "usage: grdview -J [-B] [-C] [-E]\n"); fprintf (stderr, "\t[-G] [-I] [-K] [-L[]] [-N[/]] [-O] [-P]\n"); fprintf (stderr, "\t[-Q] [-R] [-S] [-T[s]] [-U] [-V] [-W]\n"); fprintf (stderr, "\t[-X] [-Y] [-Z] [-c]\n\n"); if (GMT_quick) exit (EXIT_FAILURE); fprintf (stderr, "\t is data set to be plotted\n"); GMT_explain_option ('j'); fprintf (stderr, "\n\tOPTIONS:\n"); GMT_explain_option ('b'); fprintf (stderr, "\t-C Color palette file\n"); fprintf (stderr, "\t-E Draw perspective from viewpoint at azimuth , elevation (degrees).\n"); fprintf (stderr, "\t Default is looking straight down [180/90].\n"); fprintf (stderr, "\t-G rather than is the data set to color-code.\n"); fprintf (stderr, "\t Use as the relief and \'drape\' the image on top.\n"); fprintf (stderr, "\t Note that -Jz and -N always refers to the \n"); fprintf (stderr, "\t-I gives name of intensity file and selects illumination\n"); fprintf (stderr, "\t-L sets boundary conditions when resampling the grid. can be either\n"); fprintf (stderr, "\t g for geographic boundary conditions\n"); fprintf (stderr, "\t or one or both of\n"); fprintf (stderr, "\t x for periodic boundary conditions on x\n"); fprintf (stderr, "\t y for periodic boundary conditions on y\n"); fprintf (stderr, "\t If no are set, use bilinear rather than bicubic [Default] resampling \n"); GMT_explain_option ('Z'); GMT_explain_option ('K'); fprintf (stderr, "\t-N will draw a plane at z = level. Append color [/r/g/b] to paint\n"); fprintf (stderr, "\t the facade between the plane and the data perimeter.\n"); GMT_explain_option ('O'); GMT_explain_option ('P'); fprintf (stderr, "\t-Q sets plot reQuest. Choose one of the following:\n"); fprintf (stderr, "\t -Qm for Mesh plot [Default]. Append color for mesh paint [%d/%d/%d]\n", fill.rgb[0], fill.rgb[1], fill.rgb[2]); fprintf (stderr, "\t -Qs[m] for colored or shaded Surface. Append m to draw meshlines on the surface.\n"); fprintf (stderr, "\t -Qi for scanline converting polygons to rasterimage. Append effective dpi [100].\n"); fprintf (stderr, "\t To force a monochrome image using the YIQ transformation, append g\n"); GMT_explain_option ('R'); fprintf (stderr, "\t-S will smooth contours first (see grdcontour for value info).\n"); fprintf(stderr, "\t-T will image the data without interpolation by painting polygonal tiles\n"); GMT_explain_option ('U'); GMT_explain_option ('V'); fprintf (stderr, "\t-W sets pen attributes. can be c for contours or m for mesh\n"); fprintf (stderr, "\t -Wc draw scontours on top of surface or mesh. [Default is no contours]\n"); fprintf (stderr, "\t Optionally append pen attributes [width = %lgp, color = (%d/%d/%d), solid line].\n", pen[0].width, pen[0].rgb[0], pen[0].rgb[1], pen[0].rgb[2]); fprintf (stderr, "\t -Wm sets attributes for mesh lines [[width = %lgp, color = (%d/%d/%d), solid line].\n", pen[1].width, pen[1].rgb[0], pen[1].rgb[1], pen[1].rgb[2]); fprintf (stderr, "\t Requires -Qm or -Qsm to take effect.\n"); GMT_explain_option ('X'); fprintf (stderr, "\t-Z For 3-D plots: Set the z-level of map [0]\n"); GMT_explain_option ('c'); GMT_explain_option ('.'); exit (EXIT_FAILURE); } if (tiling && q_set) { /* Gave both -Q and -T */ fprintf (stderr, "%s: GMT ERROR: -Q and -T are mutually exclusive options\n", GMT_program); error++; } if (!topofile) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify input file\n", GMT_program); error++; } if (drape && !drapefile) { fprintf (stderr, "%s: GMT SYNTAX ERROR option -G: Must specify drape file\n", GMT_program); error++; } if (intens && !intensfile) { fprintf (stderr, "%s: GMT SYNTAX ERROR option -I: Must specify intensity file\n", GMT_program); error++; } if ((surface || image || draw_contours) && !cpt_file) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify color palette table\n", GMT_program); error++; } if (image && dpi_i <= 0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -Qi option: Must specify positive dpi\n", GMT_program); error++; } if (smooth < 0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -S option: smooth value must be positive\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 (error) exit (EXIT_FAILURE); GMT_put_history (argc, argv); /* Update .gmtcommands */ if (cpt_file) { GMT_read_cpt (cpt_file); if (GMT_b_and_w) monochrome = TRUE; } get_contours = ( (mesh && draw_contours) || (surface && GMT_n_colors > 1) ); two = (drape && get_contours) ? 2 : 0; /* Must read topofile with 2 boundary columns/rows */ if (!strcmp (topofile, "=")) { fprintf (stderr, "%s: Piping of topofile not supported!\n", GMT_program); exit (EXIT_FAILURE); } if (GMT_read_grd_info (topofile, &t_head)) { fprintf (stderr, "%s: Error opening file %s\n", GMT_program, topofile); exit (EXIT_FAILURE); } if (intens && GMT_read_grd_info (intensfile, &i_head)) { fprintf (stderr, "%s: Error opening file %s\n", GMT_program, intensfile); exit (EXIT_FAILURE); } if (drape && GMT_read_grd_info (drapefile, &d_head)) { fprintf (stderr, "%s: Error opening file %s\n", GMT_program, drapefile); exit (EXIT_FAILURE); } off = (t_head.node_offset) ? 0 : 1; del_off = (t_head.node_offset) ? 0.5 : 0.0; GMT_boundcond_init (&edgeinfo); /* Determine what wesn to pass to map_setup */ if (!project_info.region_supplied) { west = t_head.x_min; east = t_head.x_max; south = t_head.y_min; north = t_head.y_max; } else if (!(west == t_head.x_min && east == t_head.x_max && south == t_head.y_min && north == t_head.y_max)) subset = TRUE; if (project_info.z_bottom == 0.0 && project_info.z_top == 0.0) { project_info.z_bottom = t_head.z_min; project_info.z_top = t_head.z_max; if (draw_plane && plane_level < project_info.z_bottom) project_info.z_bottom = plane_level; if (draw_plane && plane_level > project_info.z_top) project_info.z_top = plane_level; } GMT_map_setup (west, east, south, north); /* Determine the wesn to be used to read the grdfile */ GMT_grd_setregion (&t_head, &data_west, &data_east, &data_south, &data_north); /* Read data */ if (gmtdefs.verbose) fprintf (stderr, "%s: Processing shape file\n", GMT_program); nx = irint ( (data_east - data_west) / t_head.x_inc) + off; ny = irint ( (data_north - data_south) / t_head.y_inc) + off; mx = nx + 2 * two; my = ny + 2 * two; nm = nx * ny; nm2 = mx * my; topo = (float *) GMT_memory (VNULL, (size_t)nm2, sizeof (float), GMT_program); GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = two; if (GMT_read_grd (topofile, &t_head, topo, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) { fprintf (stderr, "%s: Error reading file %s\n", GMT_program, topofile); exit (EXIT_FAILURE); } GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 0; nx_f = t_head.nx; ny_f = t_head.ny; delx = (t_head.node_offset) ? 0.5 * t_head.x_inc :0.0; dely = (t_head.node_offset) ? 0.5 * t_head.y_inc :0.0; xval = (double *) GMT_memory (VNULL, (size_t)t_head.nx, sizeof (double), GMT_program); yval = (double *) GMT_memory (VNULL, (size_t)t_head.ny, sizeof (double), GMT_program); for (i = 0; i < t_head.nx; i++) xval[i] = t_head.x_min + i * t_head.x_inc + delx; for (j = 0; j < t_head.ny; j++) yval[j] = t_head.y_max - j * t_head.y_inc - dely; GMT_boundcond_param_prep (&t_head, &edgeinfo); if (drape) { /* draping wanted */ if (gmtdefs.verbose) fprintf (stderr, "%s: Processing drape file\n", GMT_program); grd = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program); if (GMT_read_grd (drapefile, &d_head, grd, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) { fprintf (stderr, "%s: Error reading file %s\n", GMT_program, drapefile); exit (EXIT_FAILURE); } if (d_head.nx != nx_f || d_head.ny != ny_f) { fprintf (stderr, "%s: Drape file has improper dimensions!\n", GMT_program); exit (EXIT_FAILURE); } header = d_head; } else { grd = topo; header = t_head; } if (!project_info.xyz_pos[2]) /* Negative z-scale, must flip */ d_swap (project_info.z_bottom, project_info.z_top); ij_inc[0] = 0; ij_inc[1] = 1; ij_inc[2] = 1 - mx; ij_inc[3] = -mx; bin_inc[0] = 0; bin_inc[1] = 1; bin_inc[2] = 1 - header.nx; bin_inc[3] = -header.nx; nw = two * mx + two; ne = nw + t_head.nx - 1; sw = (t_head.ny + two - 1) * mx + two; se = sw + t_head.nx - 1; grdview_init_setup (&t_head, topo, two, draw_plane, plane_level); /* Find projected min/max in y-direction */ i_start = (z_project.quadrant == 1 || z_project.quadrant == 2) ? 0 : header.nx - 2; i_stop = (z_project.quadrant == 1 || z_project.quadrant == 2) ? header.nx - 1 : -1; i_inc = (z_project.quadrant == 1 || z_project.quadrant == 2) ? 1 : -1; j_start = (z_project.quadrant == 1 || z_project.quadrant == 4) ? header.ny - 1 : 1; j_stop = (z_project.quadrant == 1 || z_project.quadrant == 4) ? 0 : header.ny; j_inc = (z_project.quadrant == 1 || z_project.quadrant == 4) ? -1 : 1; x_inc[0] = x_inc[3] = 0.0; x_inc[1] = x_inc[2] = header.x_inc; y_inc[0] = y_inc[1] = 0.0; y_inc[2] = y_inc[3] = header.y_inc; if (get_contours) { /* Need to find contours */ if (gmtdefs.verbose) fprintf (stderr, "%s: Find contours\n", GMT_program); n_edges = header.ny * (int )ceil (header.nx / 16.0); edge = (int *) GMT_memory (VNULL, (size_t)n_edges, sizeof (int), GMT_program); binij = (struct BIN *) GMT_memory (VNULL, (size_t)nm, sizeof (struct BIN), GMT_program); small = SMALL * (header.z_max - header.z_min); if (gmtdefs.verbose) fprintf (stderr, "%s: Trace and bin contours...\n", GMT_program); first = TRUE; for (c = 0; c < GMT_n_colors+1; c++) { /* For each color change */ /* Reset markers and set up new zero-contour*/ cval = (c == GMT_n_colors) ? GMT_lut[c-1].z_high : GMT_lut[c].z_low; if (cval < header.z_min || cval > header.z_max) continue; if (gmtdefs.verbose) fprintf (stderr, "%s: Now tracing contour interval %8lg\r", GMT_program, cval); take_out = (first) ? cval : cval - GMT_lut[c-1].z_low; first = FALSE; for (i = 0; i < nm; i++) { if (!GMT_is_fnan (grd[i])) grd[i] -= (float)take_out; if (grd[i] == 0.0) grd[i] += (float)small; } side = 0; begin = TRUE; while (side < 5) { while ((n = GMT_contours (grd, &header, smooth, gmtdefs.interpolant, &side, edge, begin, &x, &y)) > 0) { begin = FALSE; i_bin_old = j_bin_old = -1; for (i = 1; i < n; i++) { i_bin = (int)floor (((0.5 * (x[i-1] + x[i]) - header.x_min) / header.x_inc) - del_off); j_bin = (int)floor (((header.y_max - 0.5 * (y[i-1] + y[i])) / header.y_inc) - del_off) + 1; if (i_bin != i_bin_old || j_bin != j_bin_old) { /* Entering new bin */ bin = j_bin * header.nx + i_bin; this_cont = get_cont_struct (bin, cval); this_cont->value = cval; this_cont->first_point = get_point (x[i-1], y[i-1]); this_point = this_cont->first_point; i_bin_old = i_bin; j_bin_old = j_bin; } this_point->next_point = get_point (x[i], y[i]); this_point = this_point->next_point; } GMT_free ((void *)x); GMT_free ((void *)y); } begin = FALSE; } } /* Remove temporary variables */ GMT_free ((void *)edge); /* Go back to beginning and reread since grd has been destroyed */ if (drape) { if (GMT_read_grd_info (drapefile, &d_head)) { fprintf (stderr, "%s: Error opening file %s\n", GMT_program, drapefile); exit (EXIT_FAILURE); } if (GMT_read_grd (drapefile, &d_head, grd, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) { fprintf (stderr, "%s: Error reading file %s\n", GMT_program, drapefile); exit (EXIT_FAILURE); } header = d_head; } else { if (GMT_read_grd_info (topofile, &t_head)) { fprintf (stderr, "%s: Error opening file %s\n", GMT_program, topofile); exit (EXIT_FAILURE); } GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = two; if (GMT_read_grd (topofile, &t_head, topo, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) { fprintf (stderr, "%s: Error reading file %s\n", GMT_program, topofile); exit (EXIT_FAILURE); } GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 0; grd = topo; header = t_head; } if (gmtdefs.verbose) fprintf (stderr, "\n"); } if (intens) { /* Illumination wanted */ if (gmtdefs.verbose) fprintf (stderr, "%s: Processing illumination file\n", GMT_program); intensity = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program); if (GMT_read_grd (intensfile, &i_head, intensity, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) { fprintf (stderr, "%s: Error reading file %s\n", GMT_program, intensfile); exit (EXIT_FAILURE); } if (i_head.nx != nx_f || i_head.ny != ny_f) { fprintf (stderr, "%s: Intensity file has improper dimensions!\n", GMT_program); exit (EXIT_FAILURE); } } if (two) { /* Initialize bcr stuff */ GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = two; GMT_bcr_init (&t_head, GMT_pad, bilinear); /* Set boundary conditions */ GMT_boundcond_set (&t_head, &edgeinfo, GMT_pad, topo); GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 0; } dx2 = 0.5 * header.x_inc; dy2 = 0.5 * header.y_inc; no_nans = TRUE; /* Temporarily modify the outer edge coordinates */ xval[0] += dx2; xval[header.nx-1] -= dx2; yval[0] -= dy2; yval[header.ny-1] += dy2; for (j = bin = ij = 0; j < header.ny; j++) { /* Nodes outside -R is set to NaN and not used */ ij = (two) ? (j + 2) * mx + two : bin; for (i = 0; i < header.nx; i++, bin++, ij++) { if (GMT_is_fnan (topo[ij])) no_nans = FALSE; if (GMT_map_outside (xval[i], yval[j])) topo[ij] = GMT_f_NaN; } } /* Reset the outer edge coordinates */ xval[0] -= dx2; xval[header.nx-1] += dx2; yval[0] += dy2; yval[header.ny-1] -= dy2; max = 2 * (MAX(1,smooth) * (((header.nx > header.ny) ? header.nx : header.ny) + 2)) + 1; x = (double *) GMT_memory (VNULL, (size_t)max, sizeof (double), GMT_program); y = (double *) GMT_memory (VNULL, (size_t)max, sizeof (double), GMT_program); z = (double *) GMT_memory (VNULL, (size_t)max, sizeof (double), GMT_program); v = (double *) GMT_memory (VNULL, (size_t)max, sizeof (double), GMT_program); if (gmtdefs.verbose) fprintf (stderr, "%s: Start creating PostScript plot\n", GMT_program); 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])); ps_setformat (3); 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 (frame_info.plot && project_info.three_D) GMT_map_basemap (); /* Plot basemap first if 3-D */ if (project_info.z_pars[0] == 0.0) GMT_map_clip_on (GMT_no_rgb, 3); if (set_z) project_info.z_level = new_z_level; xx = (double *) GMT_memory (VNULL, (size_t)max, sizeof(double), GMT_program); yy = (double *) GMT_memory (VNULL, (size_t)max, sizeof(double), GMT_program); if (draw_plane) { ps_comment ("Plot the plane at desired level"); if (!z_project.draw[0]) { /* Southern side */ if (!project_info.region) { GMT_geoz_to_xy (z_project.corner_x[0], z_project.corner_y[0], plane_level, &xx[0], &yy[0]); GMT_geoz_to_xy (z_project.corner_x[1], z_project.corner_y[1], plane_level, &xx[1], &yy[1]); ps_line (xx, yy, 2, 3, TRUE, TRUE); } else { for (i = 0; i < header.nx; i++) GMT_geoz_to_xy (header.x_min + i * header.x_inc + delx, header.y_min + dely, plane_level, &xx[i], &yy[i]); ps_line (xx, yy, header.nx, 3, TRUE, TRUE); } } if (!z_project.draw[2]) { /* Northern side */ if (!project_info.region) { GMT_geoz_to_xy (z_project.corner_x[3], z_project.corner_y[3], plane_level, &xx[0], &yy[0]); GMT_geoz_to_xy (z_project.corner_x[2], z_project.corner_y[2], plane_level, &xx[1], &yy[1]); ps_line (xx, yy, 2, 3, TRUE, TRUE); } else { for (i = 0; i < header.nx; i++) GMT_geoz_to_xy (header.x_min + i * header.x_inc + delx, header.y_max - dely, plane_level, &xx[i], &yy[i]); ps_line (xx, yy, header.nx, 3, TRUE, TRUE); } } if (!z_project.draw[3]) { /* Western side */ if (!project_info.region) { GMT_geoz_to_xy (z_project.corner_x[0], z_project.corner_y[0], plane_level, &xx[0], &yy[0]); GMT_geoz_to_xy (z_project.corner_x[3], z_project.corner_y[3], plane_level, &xx[1], &yy[1]); ps_line (xx, yy, 2, 3, TRUE, TRUE); } else { for (j = 0; j < header.ny; j++) GMT_geoz_to_xy (header.x_min + delx, header.y_max - j * header.y_inc - dely, plane_level, &xx[j], &yy[j]); ps_line (xx, yy, header.ny, 3, TRUE, TRUE); } } if (!z_project.draw[1]) { /* Eastern side */ if (!project_info.region) { GMT_geoz_to_xy (z_project.corner_x[1], z_project.corner_y[1], plane_level, &xx[0], &yy[0]); GMT_geoz_to_xy (z_project.corner_x[2], z_project.corner_y[2], plane_level, &xx[1], &yy[1]); ps_line (xx, yy, 2, 3, TRUE, TRUE); } else { for (j = 0; j < header.ny; j++) GMT_geoz_to_xy (header.x_max - delx, header.y_max - j * header.y_inc - dely, plane_level, &xx[j], &yy[j]); ps_line (xx, yy, header.ny, 3, TRUE, TRUE); } } if (project_info.region) { GMT_geoz_to_xy (header.x_min + delx, header.y_min + dely, plane_level, &xx[0], &yy[0]); GMT_geoz_to_xy (header.x_max - delx, header.y_min + dely, plane_level, &xx[1], &yy[1]); GMT_geoz_to_xy (header.x_max - delx, header.y_max - dely, plane_level, &xx[2], &yy[2]); GMT_geoz_to_xy (header.x_min + delx, header.y_max - dely, plane_level, &xx[3], &yy[3]); if (!GMT_is_fnan (topo[nw])) { GMT_geoz_to_xy (header.x_min + delx, header.y_max - dely, (double)(topo[nw]), &x_left, &y_top); ps_plot (x_left, y_top, 3); ps_plot (xx[3], yy[3], 2); } if (!GMT_is_fnan (topo[ne])) { GMT_geoz_to_xy (header.x_max - delx, header.y_max - dely, (double)(topo[ne]), &x_right, &y_top); ps_plot (x_right, y_top, 3); ps_plot (xx[2], yy[2], 2); } if (!GMT_is_fnan (topo[se])) { GMT_geoz_to_xy (header.x_max - delx, header.y_min + dely, (double)(topo[se]), &x_right, &y_bottom); ps_plot (x_right, y_bottom, 3); ps_plot (xx[1], yy[1], 2); } if (!GMT_is_fnan (topo[sw])) { GMT_geoz_to_xy (header.x_min + delx, header.y_min + dely, (double)(topo[sw]), &x_left, &y_bottom); ps_plot (x_left, y_bottom, 3); ps_plot (xx[0], yy[0], 2); } ps_plot (0.0, 0.0, 3); /* To make sure path is stroked */ } } if (tiling) { /* Plot image as polygonal pieces */ int p; double *xx, *yy, d, xpos, ypos; d = (off) ? 0.0 : 0.5; if (gmtdefs.verbose) fprintf (stderr, "%s: Tiling without interpolation\n", GMT_program); for (j = k = 0; j < header.ny; j++) { for (i = 0; i < header.nx; i++, k++) { /* Compute rgb for each pixel */ if (GMT_is_fnan (topo[k]) && tiling == 2) continue; GMT_get_rgb24 (topo[k], rgb); if (intens) GMT_illuminate (intensity[k], rgb); n = GMT_graticule_path (&xx, &yy, 1, xval[i] - dx2, xval[i] + dx2, yval[j] - dy2, yval[j] + dy2); /* Do map projection and plot */ for (p = 0; p < n; p++) { GMT_geoz_to_xy (xx[p], yy[p], topo[k], &xpos, &ypos); xx[p] = xpos; yy[p] = ypos; } ps_polygon (xx, yy, n, rgb, FALSE); GMT_free ((void *)xx); GMT_free ((void *)yy); } } } if (image) { /* compute image */ int nx_i, ny_i, kk, ip, jp, min_i, max_i, min_j, max_j, dist, node, nm_i, layers, *ix, *iy; int *top_jp, *bottom_jp, p; BOOLEAN done; double xp, yp, sum_w, w, sum_i, x_width, y_width; double sum_r, sum_g, sum_b, intval; unsigned char *bitimage_24, *bitimage_8; if (gmtdefs.verbose) { if (GMT_cpt_pattern) fprintf (stderr, "%s: Warning: Patterns in cpt file will not work with -Qi\n", GMT_program); fprintf (stderr, "%s: get and store projected vertices\n", GMT_program); } ps_comment ("Plot 3-D surface using scanline conversion of polygons to raster image"); x_width = z_project.xmax - z_project.xmin; /* Size of image in inches */ y_width = z_project.ymax - z_project.ymin; nx_i = irint (x_width * dpi_i); /* Size of image in pixels */ ny_i = irint (y_width * dpi_i); ix = (int *) GMT_memory (VNULL, (size_t)nm, sizeof (int), GMT_program); iy = (int *) GMT_memory (VNULL, (size_t)nm, sizeof (int), GMT_program); for (j = bin = ij = 0; j < header.ny; j++) { /* Get projected coordinates converted to pixel locations */ ij = (two) ? (j + 2) * mx + two : bin; for (i = 0; i < header.nx; i++, bin++, ij++) { if (GMT_is_fnan (topo[ij])) { /* Outside -R or NaNs not used */ ix[bin] = iy[bin] = -1; } else { GMT_geoz_to_xy (xval[i], yval[j], (double)topo[ij], &xp, &yp); ix[bin] = (int)floor ((xp - z_project.xmin) * dpi_i); iy[bin] = (int)floor ((yp - z_project.ymin) * dpi_i); } } } GMT_free ((void *) xval); GMT_free ((void *) yval); /* Allocate image array and set background to PAGE_COLOR */ if (monochrome) { int gray; nm_i = nx_i * ny_i; layers = 1; bitimage_8 = (unsigned char *) GMT_memory (VNULL, (size_t)nm_i, sizeof (char), GMT_program); gray = YIQ (gmtdefs.page_rgb); memset ((void *)bitimage_8, gray, (size_t)nm_i); } else { nm_i = nx_i * ny_i * 3; layers = 3; bitimage_24 = (unsigned char *) GMT_memory (VNULL, (size_t)nm_i, sizeof (char), GMT_program); kk = 0; while (kk < nm_i) { bitimage_24[kk++] = (unsigned char)gmtdefs.page_rgb[0]; bitimage_24[kk++] = (unsigned char)gmtdefs.page_rgb[1]; bitimage_24[kk++] = (unsigned char)gmtdefs.page_rgb[2]; } } if (no_nans) { /* Set up arrays for staircase clippath and initialize them */ top_jp = (int *) GMT_memory (VNULL, (size_t)nx_i, sizeof (int), GMT_program); bottom_jp = (int *) GMT_memory (VNULL, (size_t)nx_i, sizeof (int), GMT_program); for (ip = 0; ip < nx_i; ip++) bottom_jp[ip] = ny_i; } /* Plot from back to front */ if (gmtdefs.verbose) fprintf (stderr, "%s: Start rasterization\n", GMT_program); for (j = j_start; j != j_stop; j += j_inc) { if (gmtdefs.verbose) fprintf (stderr, "%s: Scan line conversion at j-line %.4d\r", GMT_program, j); for (i = i_start; i != i_stop; i += i_inc) { bin = j * header.nx + i; ij = (two) ? (j + 2) * header.nx + i + 2 : bin; for (k = bad = 0; !bad && k < 4; k++) bad = (ix[bin+bin_inc[k]] < 0 || iy[bin+bin_inc[k]] < 0); if (bad) continue; min_i = max_i = ix[bin]; min_j = max_j = iy[bin]; for (k = 1; k < 4; k++) { p = bin+bin_inc[k]; if (ix[p] < min_i) min_i = ix[p]; if (ix[p] > max_i) max_i = ix[p]; if (iy[p] < min_j) min_j = iy[p]; if (iy[p] > max_j) max_j = iy[p]; } for (jp = min_j; jp <= max_j; jp++) { if (jp < 0 || jp >= ny_i) continue; for (ip = min_i; ip <= max_i; ip++) { if (ip < 0 || ip >= nx_i) continue; if (!pixel_inside (ip, jp, ix, iy, bin)) continue; if (no_nans) { /* Update clip mask */ if (jp > top_jp[ip]) top_jp[ip] = jp; if (jp < bottom_jp[ip]) bottom_jp[ip] = jp; } kk = layers * ((ny_i-jp-1) * nx_i + ip); sum_r = sum_g = sum_b = sum_w = sum_i = 0.0; done = FALSE; for (k = 0; !done && k < 4; k++) { node = bin + bin_inc[k]; GMT_get_rgb24 (grd[node], rgb); dist = quick_idist (ip, jp, ix[node], iy[node]); if (dist == 0) { /* Only need this corner value */ done = TRUE; if (intens) intval = intensity[node]; } else { w = 1.0 / (double)dist; sum_r += rgb[0] * w; sum_g += rgb[1] * w; sum_b += rgb[2] * w; if (intens) sum_i += intensity[node] * w; sum_w += w; } } if (!done) { /* Must get weighted value */ sum_w = 1.0 / sum_w; rgb[0] = irint (sum_r * sum_w); rgb[1] = irint (sum_g * sum_w); rgb[2] = irint (sum_b * sum_w); if (intens) intval = sum_i * sum_w; } if (intens) GMT_illuminate (intval, rgb); if (monochrome) /* YIQ transformation */ bitimage_8[kk] = (unsigned char) YIQ (rgb); else { bitimage_24[kk++] = (unsigned char) rgb[0]; bitimage_24[kk++] = (unsigned char) rgb[1]; bitimage_24[kk] = (unsigned char) rgb[2]; } } } } } if (gmtdefs.verbose) fprintf (stderr, "\n"); if (no_nans) { /* Must implement the clip path for the image */ x_pixel_size = x_width / (double)nx_i; y_pixel_size = y_width / (double)ny_i; n4 = 4 * nx_i; x_imask = (double *) GMT_memory (VNULL, (size_t)n4, sizeof (double), GMT_program); y_imask = (double *) GMT_memory (VNULL, (size_t)n4, sizeof (double), GMT_program); nk = n4 - 1; for (ip = k = 0; ip < nx_i; ip++, k+= 2) { k1 = k + 1; x_imask[k] = x_imask[nk-k] = z_project.xmin + ip * x_pixel_size; x_imask[k1] = x_imask[nk-k1] = x_imask[k] + x_pixel_size; if (top_jp[ip] < bottom_jp[ip]) { /* No pixels set in this column */ y_imask[k] = y_imask[k1] = y_imask[nk-k] = y_imask[nk-k1] = z_project.ymin; } else { /* Set top of upper pixel and bottom of lower pixel */ y_imask[k] = y_imask[k1] = z_project.ymin + (top_jp[ip] + 1) * y_pixel_size; y_imask[nk-k] = y_imask[nk-k1] = z_project.ymin + bottom_jp[ip] * y_pixel_size; } } ps_clipon (x_imask, y_imask, 4 * nx_i, GMT_no_rgb, 3); GMT_free ((void *)x_imask); GMT_free ((void *)y_imask); } if (gmtdefs.verbose) fprintf (stderr, "%s: Creating PostScript image ", GMT_program); if (monochrome) { if (gmtdefs.verbose) fprintf (stderr, "[B/W image]\n"); ps_image (z_project.xmin, z_project.ymin, x_width, y_width, bitimage_8, nx_i, ny_i, 8); GMT_free ((void *)bitimage_8); } else { if (gmtdefs.verbose) fprintf (stderr, "[%s]\n", c_method[gmtdefs.color_image]); GMT_color_image (z_project.xmin, z_project.ymin, x_width, y_width, bitimage_24, nx_i, ny_i); GMT_free ((void *)bitimage_24); } if (no_nans){ ps_clipoff (); GMT_free ((void *)top_jp); GMT_free ((void *)bottom_jp); } GMT_free ((void *)ix); GMT_free ((void *)iy); } if (mesh) { ps_comment ("Start of mesh plot"); GMT_setpen (&pen[1]); if (monochrome) fill.rgb[0] = fill.rgb[1] = fill.rgb[2] = YIQ (fill.rgb); /* Do YIQ transformation */ for (j = j_start; j != j_stop; j += j_inc) { y_bottom = yval[j]; y_top = y_bottom + abs (j_inc) * header.y_inc; for (i = i_start; i != i_stop; i += i_inc) { bin = j * header.nx + i; ij = (two) ? (j + 2) * mx + i + 2 : bin; for (k = bad = 0; !bad && k < 4; k++) bad += GMT_is_fnan (topo[ij+ij_inc[k]]); if (bad) continue; x_left = xval[i]; x_right = x_left + abs (i_inc) * header.x_inc; GMT_geoz_to_xy (x_left, y_bottom, (double)(topo[ij+ij_inc[0]]), &xx[0], &yy[0]); GMT_geoz_to_xy (x_right, y_bottom, (double)(topo[ij+ij_inc[1]]), &xx[1], &yy[1]); GMT_geoz_to_xy (x_right, y_top, (double)(topo[ij+ij_inc[2]]), &xx[2], &yy[2]); GMT_geoz_to_xy (x_left, y_top, (double)(topo[ij+ij_inc[3]]), &xx[3], &yy[3]); ps_patch (xx, yy, 4, fill.rgb, TRUE); if (draw_contours) { pen_not_set = TRUE; if (binij[bin].first_cont == NULL) continue; for (this_cont = binij[bin].first_cont->next_cont; this_cont; this_cont = this_cont->next_cont) { for (k = 0, this_point = this_cont->first_point; this_point; this_point = this_point->next_point) { z_val = (drape) ? GMT_get_bcr_z (&t_head, (double)this_point->x, (double)this_point->y, topo, &edgeinfo) : this_cont->value; if (GMT_is_dnan (z_val)) continue; GMT_geoz_to_xy ((double)this_point->x, (double)this_point->y, z_val, &xx[k], &yy[k]); k++; } if (pen_not_set) { GMT_setpen (&pen[0]); pen_not_set = FALSE; } ps_line (xx, yy, k, 3, FALSE, TRUE); } if (!pen_not_set) GMT_setpen (&pen[1]); } } } GMT_free ((void *) xval); GMT_free ((void *) yval); } else if (surface) { int start_side, entry_side, exit_side, next_side, low, ncont, nw_se_diagonal, check; int corner[2], bad_side[2][2], p, p1, p2, saddle_sign; double *xcont, *ycont, *zcont, *vcont, X_vert[4], Y_vert[4], saddle_small; xcont = (double *) GMT_memory (VNULL, (size_t)max, sizeof (double), GMT_program); ycont = (double *) GMT_memory (VNULL, (size_t)max, sizeof (double), GMT_program); zcont = (double *) GMT_memory (VNULL, (size_t)max, sizeof (double), GMT_program); vcont = (double *) GMT_memory (VNULL, (size_t)max, sizeof (double), GMT_program); ps_comment ("Start of filled surface"); if (outline) GMT_setpen (&pen[1]); for (j = j_start; j != j_stop; j += j_inc) { y_bottom = yval[j]; y_top = y_bottom + header.y_inc; for (i = i_start; i != i_stop; i += i_inc) { bin = j * header.nx + i; ij = (two) ? (j + 2) * mx + i + 2 : bin; x_left = xval[i]; x_right = x_left + header.x_inc; for (k = bad = 0; !bad && k < 4; k++) bad += GMT_is_fnan (topo[ij+ij_inc[k]]); if (bad) { if (GMT_bfn.skip[2] || project_info.three_D) continue; X_vert[0] = X_vert[3] = x_left; X_vert[1] = X_vert[2] = x_right; Y_vert[0] = Y_vert[1] = y_bottom; Y_vert[2] = Y_vert[3] = y_top; for (k = 0; k < 4; k++) GMT_geoz_to_xy (X_vert[k], Y_vert[k], 0.0, &xmesh[k], &ymesh[k]); paint_it (xmesh, ymesh, 4, GMT_d_NaN, FALSE, 0.0); if (outline) ps_patch (xmesh, ymesh, 4, GMT_no_rgb, TRUE); continue; } if (intens) { this_intensity = get_intensity (intensity, bin, header.nx); if (GMT_is_dnan (this_intensity)) continue; } /* Get mesh polygon */ X_vert[0] = X_vert[3] = x_left; X_vert[1] = X_vert[2] = x_right; Y_vert[0] = Y_vert[1] = y_bottom; Y_vert[2] = Y_vert[3] = y_top; if (get_contours && binij[bin].first_cont) { /* Contours go thru here */ /* Determine if this bin will give us saddle trouble */ start_cont = this_cont = binij[bin].first_cont->next_cont; saddle = FALSE; while (!saddle && this_cont->next_cont) { if (this_cont->next_cont->value == this_cont->value) saddle = TRUE; else this_cont = this_cont->next_cont; } if (saddle) { /* Must deal with this separately */ this_point = this_cont->first_point; entry_side = get_side (this_point->x, this_point->y, x_left, y_bottom, header.x_inc, header.y_inc, dx2, dy2); while (this_point->next_point) this_point = this_point->next_point; /* Go to end */ exit_side = get_side (this_point->x, this_point->y, x_left, y_bottom, header.x_inc, header.y_inc, dx2, dy2); if (MIN (grd[bin+bin_inc[1]], grd[bin+bin_inc[3]]) > MAX (grd[bin], grd[bin+bin_inc[2]])) { saddle_sign = +1; check = TRUE; } else if (MAX (grd[bin+bin_inc[1]], grd[bin+bin_inc[3]]) < MIN (grd[bin], grd[bin+bin_inc[2]])) { saddle_sign = -1; check = TRUE; } else if (MIN (grd[bin], grd[bin+bin_inc[2]]) > MAX (grd[bin+bin_inc[1]], grd[bin+bin_inc[3]])) { saddle_sign = +1; check = FALSE; } else { saddle_sign = -1; check = FALSE; } nw_se_diagonal = ((entry_side + exit_side) == 3); if (nw_se_diagonal != check) saddle_sign = -saddle_sign; if (nw_se_diagonal) { /* Diagonal goes NW - SE */ corner[0] = 0; bad_side[0][0] = 1; bad_side[0][1] = 2; corner[1] = 2; bad_side[1][0] = 0; bad_side[1][1] = 3; } else { /* Diagonal goes NE -SW */ corner[0] = 1; bad_side[0][0] = 2; bad_side[0][1] = 3; corner[1] = 3; bad_side[1][0] = 0; bad_side[1][1] = 1; } saddle_small = saddle_sign * small; for (p = 0; p < 2; p++) { /* For each triangular half */ /* Set this points as the start anchor */ low = corner[p]; n = 0; add_node (x, y, z, v, &n, low, X_vert, Y_vert, topo, grd, ij+ij_inc[low], bin+bin_inc[low]); start_side = next_side = low; way = 0; for (this_cont = start_cont; this_cont; this_cont = this_cont->next_cont) { /* First get all the x/y pairs for this contour */ for (k = 0, this_point = this_cont->first_point; this_point; this_point = this_point->next_point) { xcont[k] = this_point->x; ycont[k] = this_point->y; zcont[k] = (drape) ? GMT_get_bcr_z (&t_head, xcont[k], ycont[k], topo, &edgeinfo) : this_cont->value; if (GMT_is_dnan (zcont[k])) continue; vcont[k] = this_cont->value; k++; } ncont = k; entry_side = get_side (xcont[0], ycont[0], x_left, y_bottom, header.x_inc, header.y_inc, dx2, dy2); exit_side = get_side (xcont[ncont-1], ycont[ncont-1], x_left, y_bottom, header.x_inc, header.y_inc, dx2, dy2); if (entry_side == bad_side[p][0] || entry_side == bad_side[p][1]) continue; if (exit_side == bad_side[p][0] || exit_side == bad_side[p][1]) continue; /* OK, got the correct contour */ next_up = (this_cont->next_cont) ? this_cont->next_cont->value : DBL_MAX; exit_side = get_side (xcont[ncont-1], ycont[ncont-1], x_left, y_bottom, header.x_inc, header.y_inc, dx2, dy2); if (way == 0 || next_side == entry_side) { /* Just hook up */ copy_points_fw (x, y, z, v, xcont, ycont, zcont, vcont, ncont, &n); next_side = exit_side; } else if (next_side == exit_side) { /* Just hook up but reverse */ copy_points_bw (x, y, z, v, xcont, ycont, zcont, vcont, ncont, &n) ; next_side = entry_side; } /* Compute the xy from the xyz triplets */ for (k = 0; k < n; k++) GMT_geoz_to_xy (x[k], y[k], z[k], &xx[k], &yy[k]); z_ave = (GMT_continuous) ? get_z_ave (v, next_up, n) : this_cont->value; /* Now paint the polygon piece */ paint_it (xx, yy, n, z_ave-saddle_small, intens, this_intensity); /* Reset the anchor points to previous contour */ n = 0; copy_points_fw (x, y, z, v, xcont, ycont, zcont, vcont, ncont, &n); next_side = exit_side; start_side = entry_side; way = (grd[bin+bin_inc[low]] < this_cont->value) ? -1 : 1; } /* Final contour needs to add diagonal */ if (corner[p] == 0 || corner[p] == 2) { p1 = (next_side < 2) ? 1 : 3; p2 = (next_side < 2) ? 3 : 1; } else { p1 = (next_side % 3) ? 2 : 0; p2 = (next_side % 3) ? 0 : 2; } add_node (x, y, z, v, &n, p1, X_vert, Y_vert, topo, grd, ij+ij_inc[p1], bin+bin_inc[p1]); add_node (x, y, z, v, &n, p2, X_vert, Y_vert, topo, grd, ij+ij_inc[p2], bin+bin_inc[p2]); /* Compute the xy from the xyz triplets */ for (k = 0; k < n; k++) GMT_geoz_to_xy (x[k], y[k], z[k], &xx[k], &yy[k]); z_ave = (GMT_continuous) ? get_z_ave (v, next_up, n) : v[0]; /* Now paint the polygon piece */ paint_it (xx, yy, n, z_ave+saddle_small, intens, this_intensity); } /* End triangular piece */ } /* End Saddle section */ else { /* Ok, here we do not have to worry about saddles */ /* Find lowest corner (id = low) */ for (k = 1, low = 0; k < 4; k++) if (grd[bin+bin_inc[k]] < grd[bin+bin_inc[low]]) low = k; /* Set this points as the start anchor */ n = 0; add_node (x, y, z, v, &n, low, X_vert, Y_vert, topo, grd, ij+ij_inc[low], bin+bin_inc[low]); start_side = next_side = low; way = 1; this_cont = start_cont; while (this_cont) { next_up = (this_cont->next_cont) ? this_cont->next_cont->value : DBL_MAX; /* First get all the x/y pairs for this contour */ for (k = 0, this_point = this_cont->first_point; this_point; this_point = this_point->next_point) { xcont[k] = this_point->x; ycont[k] = this_point->y; zcont[k] = (drape) ? GMT_get_bcr_z (&t_head, xcont[k], ycont[k], topo, &edgeinfo) : this_cont->value; if (GMT_is_dnan (zcont[k])) continue; vcont[k] = this_cont->value; k++; } ncont = k; entry_side = get_side (xcont[0], ycont[0], x_left, y_bottom, header.x_inc, header.y_inc, dx2, dy2); exit_side = get_side (xcont[ncont-1], ycont[ncont-1], x_left, y_bottom, header.x_inc, header.y_inc, dx2, dy2); while (!(next_side == entry_side || next_side == exit_side)) { /* Must add intervening corner */ if (way == 1) next_side = (next_side + 1) % 4; add_node (x, y, z, v, &n, next_side, X_vert, Y_vert, topo, grd, ij+ij_inc[next_side], bin+bin_inc[next_side]); if (way == -1) next_side = (next_side - 1 + 4) % 4; } if (next_side == entry_side) { /* Just hook up */ copy_points_fw (x, y, z, v, xcont, ycont, zcont, vcont, ncont, &n); next_side = exit_side; } else if (next_side == exit_side) { /* Just hook up but reverse */ copy_points_bw (x, y, z, v, xcont, ycont, zcont, vcont, ncont, &n); next_side = entry_side; } /* Now we must complete the polygon if necessary */ while (!(start_side == next_side)) { /* Must add intervening corner */ if (way == 1) next_side = (next_side + 1) % 4; add_node (x, y, z, v, &n, next_side, X_vert, Y_vert, topo, grd, ij+ij_inc[next_side], bin+bin_inc[next_side]); if (way == -1) next_side = (next_side - 1 + 4) % 4; } /* Compute the xy from the xyz triplets */ for (k = 0; k < n; k++) GMT_geoz_to_xy (x[k], y[k], z[k], &xx[k], &yy[k]); z_ave = (GMT_continuous) ? get_z_ave (v, next_up, n) : this_cont->value; /* Now paint the polygon piece */ paint_it (xx, yy, n, z_ave-small, intens, this_intensity); /* Reset the anchor points to previous contour */ n = 0; copy_points_fw (x, y, z, v, xcont, ycont, zcont, vcont, ncont, &n); next_side = exit_side; start_side = entry_side; way = (grd[bin+bin_inc[start_side]] < this_cont->value) ? -1 : 1; this_cont = this_cont->next_cont; /* Goto next contour */ } /* Final contour needs to compete with corners only */ while (!(start_side == next_side)) { /* Must add intervening corner */ if (way == 1) next_side = (next_side +1) % 4; add_node (x, y, z, v, &n, next_side, X_vert, Y_vert, topo, grd, ij+ij_inc[next_side], bin+bin_inc[next_side]); if (way == -1) next_side = (next_side - 1 + 4) % 4; } /* Compute the xy from the xyz triplets */ for (k = 0; k < n; k++) GMT_geoz_to_xy (x[k], y[k], z[k], &xx[k], &yy[k]); z_ave = (GMT_continuous) ? get_z_ave (v, next_up, n) : v[0]; /* Now paint the polygon piece */ paint_it (xx, yy, n, z_ave+small, intens, this_intensity); } /* End non-saddle case */ /* Draw contour lines if desired */ pen_not_set = TRUE; for (this_cont = start_cont; draw_contours && this_cont; this_cont = this_cont->next_cont) { for (k = 0, this_point = this_cont->first_point; this_point; this_point = this_point->next_point) { z_val = (drape) ? GMT_get_bcr_z (&t_head, (double)this_point->x, (double)this_point->y, topo, &edgeinfo) : this_cont->value; if (GMT_is_dnan (z_val)) continue; GMT_geoz_to_xy ((double)this_point->x, (double)this_point->y, z_val, &xx[k], &yy[k]); k++; } if (pen_not_set) { GMT_setpen (&pen[0]); pen_not_set = FALSE; } ps_line (xx, yy, k, 3, FALSE, TRUE); } if (!pen_not_set) GMT_setpen (&pen[1]); if (outline) { for (k = 0; k < 4; k++) GMT_geoz_to_xy (X_vert[k], Y_vert[k], (double)(topo[ij+ij_inc[k]]), &xmesh[k], &ymesh[k]); ps_patch (xmesh, ymesh, 4, GMT_no_rgb, TRUE); } } else { /* No Contours, paint tile (any of the corner values will give correct color */ if (GMT_continuous) { for (n = 0, z_ave = 0.0; n < 4; n++) z_ave += grd[bin+bin_inc[n]]; z_ave *= 0.25; } else z_ave = grd[bin]; /* Now paint the polygon piece */ for (k = 0; k < 4; k++) GMT_geoz_to_xy (X_vert[k], Y_vert[k], (double)(topo[ij+ij_inc[k]]), &xmesh[k], &ymesh[k]); paint_it (xmesh, ymesh, 4, z_ave, intens, this_intensity); if (outline) ps_patch (xmesh, ymesh, 4, GMT_no_rgb, TRUE); } } } GMT_free ((void *) xval); GMT_free ((void *) yval); GMT_free ((void *) xcont); GMT_free ((void *) ycont); GMT_free ((void *) zcont); GMT_free ((void *) vcont); } if (pen[1].texture || pen[0].texture) ps_setdash (CNULL, 0); if (project_info.z_pars[0] == 0.0) GMT_map_clip_off(); if (facade) { /* Cover the two front sides */ ps_comment ("Paiting the frontal facade"); GMT_setpen (&pen[0]); if (!z_project.draw[0]) { /* Southern side */ for (i = n = 0, ij = sw; i < header.nx; i++, ij++) { if (GMT_is_fnan (topo[ij])) continue; GMT_geoz_to_xy (header.x_min+i*header.x_inc+delx, header.y_min+dely, (double)(topo[ij]), &xx[n], &yy[n]); n++; } for (i = header.nx - 1; i >= 0; i--, n++) GMT_geoz_to_xy (header.x_min+i*header.x_inc+delx, header.y_min+dely, plane_level, &xx[n], &yy[n]); ps_polygon (xx, yy, n, rgb_facade, TRUE); } if (!z_project.draw[1]) { /* Eastern side */ for (j = n = 0, ij = ne; j < header.ny; j++, ij += mx) { if (GMT_is_fnan (topo[ij])) continue; GMT_geoz_to_xy (header.x_max-delx, header.y_max-j*header.y_inc-dely, (double)(topo[ij]), &xx[n], &yy[n]); n++; } for (j = header.ny - 1; j >= 0; j--, n++) GMT_geoz_to_xy (header.x_max-delx, header.y_max-j*header.y_inc-dely, plane_level, &xx[n], &yy[n]); ps_polygon (xx, yy, n, rgb_facade, TRUE); } if (!z_project.draw[2]) { /* Northern side */ for (i = n = 0, ij = nw; i < header.nx; i++, ij++) { if (GMT_is_fnan (topo[ij])) continue; GMT_geoz_to_xy (header.x_min+i*header.x_inc+delx, header.y_max-dely, (double)(topo[i]), &xx[n], &yy[n]); n++; } for (i = header.nx - 1; i >= 0; i--, n++) GMT_geoz_to_xy (header.x_min+i*header.x_inc+delx, header.y_max-dely, plane_level, &xx[n], &yy[n]); ps_polygon (xx, yy, n, rgb_facade, TRUE); } if (!z_project.draw[3]) { /* Western side */ for (j = n = 0, ij = nw; j < header.ny; j++, ij += mx) { if (GMT_is_fnan (topo[ij])) continue; GMT_geoz_to_xy (header.x_min+delx, header.y_max-j*header.y_inc-dely, (double)(topo[ij]), &xx[n], &yy[n]); n++; } for (j = header.ny - 1; j >= 0; j--, n++) GMT_geoz_to_xy (header.x_min+delx, header.y_max-j*header.y_inc-dely, plane_level, &xx[n], &yy[n]); ps_polygon (xx, yy, n, rgb_facade, TRUE); } } if (project_info.three_D) GMT_vertical_axis (2); /* Draw background axis */ if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0); if (frame_info.plot && !project_info.three_D) GMT_map_basemap (); /* Plot basemap last if not 3-D */ ps_plotend (gmtdefs.last_page); /* Free memory */ if (get_contours) { for (ij = 0; ij < nm; ij++) { if (!binij[ij].first_cont) continue; last_cont = binij[ij].first_cont; for (this_cont = binij[ij].first_cont->next_cont; this_cont; this_cont = this_cont->next_cont) { if (this_cont->first_point) { last_point = this_cont->first_point; for (this_point = this_cont->first_point->next_point; this_point; this_point = this_point->next_point) { GMT_free ((void *)last_point); last_point = this_point; } GMT_free ((void *)last_point); } GMT_free ((void *)last_cont); last_cont = this_cont; } GMT_free ((void *)last_cont); } GMT_free ((void *)binij); } GMT_free ((void *)xx); GMT_free ((void *)yy); GMT_free ((void *)x); GMT_free ((void *)y); GMT_free ((void *)z); GMT_free ((void *)v); GMT_free ((void *)topo); if (intens) GMT_free ((void *)intensity); if (drape) GMT_free ((void *)grd); if (gmtdefs.verbose) fprintf (stderr, "%s: Done!\n", GMT_program); GMT_end (argc, argv); } struct CONT *get_cont_struct (int bin, double value) { struct CONT *cont, *new_cont; if (!binij[bin].first_cont) binij[bin].first_cont = (struct CONT *) GMT_memory (VNULL, (size_t)1, sizeof (struct CONT), GMT_program); for (cont = binij[bin].first_cont; cont->next_cont && cont->next_cont->value <= value; cont = cont->next_cont); new_cont = (struct CONT *) GMT_memory (VNULL, (size_t)1, sizeof (struct CONT), GMT_program); if (cont->next_cont) { /* Put it in the link */ new_cont->next_cont = cont->next_cont; cont->next_cont = new_cont; } else /* End of list */ cont->next_cont = new_cont; return (new_cont); } struct POINT *get_point (double x, double y) { struct POINT *point; point = (struct POINT *) GMT_memory (VNULL, (size_t)1, sizeof (struct POINT), GMT_program); point->x = x; point->y = y; return (point); } void grdview_init_setup (struct GRD_HEADER *header, float *topo, int two, BOOLEAN draw_plane, double plane_level) { int i, j, ij; double xtmp, ytmp, tmp, delx, dely; delx = (header->node_offset) ? 0.5 * header->x_inc :0.0; dely = (header->node_offset) ? 0.5 * header->y_inc :0.0; /* Find projected min/max in y-direction */ z_project.ymax = z_project.ymin; /* Reset from whatever it was */ for (j = 0; j < header->ny; j++) { ij = (j + two) * header->nx + two; for (i = 0; i < header->nx; i++, ij++) { if (GMT_is_fnan (topo[ij])) continue; GMT_geoz_to_xy (header->x_min + i * header->x_inc, header->y_max - j * header->y_inc, (double)topo[ij], &xtmp, &ytmp); z_project.ymin = MIN (z_project.ymin, ytmp); z_project.ymax = MAX (z_project.ymax, ytmp); } } if (draw_plane) { /* plane or facade may exceed the found min/max */ for (i = 0; i < header->nx; i++) { tmp = header->x_min + i * header->x_inc + delx; GMT_geoz_to_xy (tmp, header->y_min + dely, plane_level, &xtmp, &ytmp); z_project.ymin = MIN (z_project.ymin, ytmp); z_project.ymax = MAX (z_project.ymax, ytmp); GMT_geoz_to_xy (tmp, header->y_max-dely, plane_level, &xtmp, &ytmp); z_project.ymin = MIN (z_project.ymin, ytmp); z_project.ymax = MAX (z_project.ymax, ytmp); } for (j = 0; j < header->ny; j++) { tmp = header->y_max - j * header->y_inc - dely; GMT_geoz_to_xy (header->x_min+delx, tmp, plane_level, &xtmp, &ytmp); z_project.ymin = MIN (z_project.ymin, ytmp); z_project.ymax = MAX (z_project.ymax, ytmp); GMT_geoz_to_xy (header->x_max-delx, tmp, plane_level, &xtmp, &ytmp); z_project.ymin = MIN (z_project.ymin, ytmp); z_project.ymax = MAX (z_project.ymax, ytmp); } } } double get_intensity (float *intensity, int k, int nx) { /* Finds the agerage intensity for this polygon */ return (0.25 * (intensity[k] + intensity[k+1] + intensity[k-nx] + intensity[k-nx+1])); } int pixel_inside (int ip, int jp, int *ix, int *iy, int bin) { int i, what; double x[6], y[6]; for (i = 0; i < 4; i++) { x[i] = ix[bin+bin_inc[i]]; y[i] = iy[bin+bin_inc[i]]; } x[4] = x[0]; y[4] = y[0]; what = GMT_non_zero_winding ((double)ip, (double)jp, x, y, 5); return (what); } int quick_idist (int x1, int y1, int x2, int y2) { if ((x2 -= x1) < 0) x2 = -x2; if ((y2 -= y1) < 0) y2 = -y2; return (x2 + y2 - (((x2 > y2) ? y2 : x2) >> 1)); } int get_side (double x, double y, double x_left, double y_bottom, double xinc, double yinc, double dx2, double dy2) { /* Figure out on what side this point sites on */ double del_x, del_y; int side; del_x = x - x_left; if (del_x > dx2) del_x = xinc - del_x; del_y = y - y_bottom; if (del_y > dy2) del_y = yinc - del_y; if (del_x < del_y) /* Cutting N-S gridlines */ side = ((x-x_left) > dx2) ? 1 : 3; else /* Cutting E-W gridlines */ side = ((y-y_bottom) > dy2) ? 2 : 0; return (side); } void copy_points_fw (double x[], double y[], double z[], double v[], double xcont[], double ycont[], double zcont[], double vcont[], int ncont, int *n) { int k; for (k = 0; k < ncont; k++, (*n)++) { x[*n] = xcont[k]; y[*n] = ycont[k]; z[*n] = zcont[k]; v[*n] = vcont[k]; } } void copy_points_bw (double x[], double y[], double z[], double v[], double xcont[], double ycont[], double zcont[], double vcont[], int ncont, int *n) { int k; for (k = ncont - 1; k >= 0; k--, (*n)++) { x[*n] = xcont[k]; y[*n] = ycont[k]; z[*n] = zcont[k]; v[*n] = vcont[k]; } } double get_z_ave (double v[], double next_up, int n) { int k; double z_ave; for (k = 0, z_ave = 0.0; k < n; k++) z_ave += MIN (v[k], next_up); return (z_ave / n); } void add_node (double x[], double y[], double z[], double v[], int *k, int node, double X_vert[], double Y_vert[], float topo[], float grd[], int ij, int bin) { /* Adds a corner node to list of points and increments counter */ x[*k] = X_vert[node]; y[*k] = Y_vert[node]; z[*k] = topo[ij]; v[*k] = grd[bin]; (*k)++; } void paint_it (double x[], double y[], int n, double z, BOOLEAN intens, double intensity) { int index, 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 ((index >= 0 && (f = GMT_lut[index].fill)) || (index < 0 && (f = GMT_bfn.fill[index+3]))) { /* Pattern */ GMT_fill (x, y, n, f, FALSE); } else { /* Solid color/gray */ if (intens) GMT_illuminate (intensity, rgb); if (monochrome) rgb[0] = rgb[1] = rgb[2] = YIQ (rgb); /* YIQ transformation */ ps_patch (x, y, n, rgb, FALSE); /* Contours drawn separately (after this call) if desired */ } }