home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdview.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-07-05  |  45.1 KB  |  1,231 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdview.c    2.43  05 Jul 1995
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  * grdview will read a topofile and produce a 3-D perspective plot
  9.  * of the surface z = f(x,y) using PostScript. The surface can
  10.  * be represented as:
  11.  *    1) A Mesh plot
  12.  *    2) A shaded (or colored) surface w/wo contourlines and w/wo
  13.  *       illumination by artificial sun(s).
  14.  *
  15.  * grdview calls contours to find the line segments that make up the
  16.  * contour lines. This allows the user to specify that the contours
  17.  * should be smoothed before plotting. This will make the resulting
  18.  * image smoother, especially if nx and ny are relatively small.
  19.  *
  20.  * Author:    Paul Wessel
  21.  * Date:    8-DEC-1988
  22.  * Modified:    5-JUL-1994    New version 3.0
  23.  *
  24.  */
  25.  
  26. #include "gmt.h"
  27. #include "gmt_bcr.h"
  28.  
  29. #define YIQ(r, g, b) rint (0.299 * (r) + 0.587 * (g) + 0.114 * (b))    /* How B/W TV's convert RGB to Gray */
  30.  
  31. /* Declarations needed for binning of smooth contours */
  32.  
  33. struct BIN {
  34.     struct CONT *first_cont;
  35. } *binij;
  36.  
  37. struct CONT {
  38.     struct POINT *first_point;
  39.     struct CONT *next_cont;
  40.     double value;
  41. };
  42.  
  43. struct POINT {
  44.     double x, y;
  45.     struct POINT *next_point;
  46. };
  47.  
  48. /* Global variables */
  49.  
  50. float *grd, *intensity, *topo;
  51.  
  52. int *edge;
  53. int bin_inc[4], ij_inc[4];
  54. double x_inc[4], y_inc[4], *x, *y, *z, *v, *xx, *yy;
  55.  
  56. char *c_method[3] = {
  57.         "colorimage",
  58.         "colortiles",
  59.         "RGB-separation"
  60. };
  61.  
  62. main (argc, argv)
  63. int argc;
  64. char **argv; {
  65.     int     i, j, ij, n_edges, k, n, max, i_bin, j_bin, i_bin_old, j_bin_old;
  66.     int    smooth = 1, side, way, nm, dummy[4], nx_f, ny_f, off, nx, ny;
  67.     int    red_facade = 0, green_facade = 0, blue_facade = 0, bin, two, mx, my;
  68.     int    c, r, g, b, i_start, i_stop, j_start, j_stop, i_inc, j_inc;
  69.     int    n_saddle, prev_side, entry, check, n_added, n_decimals = 3, dpi_i = 100;
  70.     int    sw, se, nw, ne, pad[4], index;
  71.     
  72.     BOOLEAN mesh = FALSE, draw_plane = FALSE, facade = FALSE, get_contours, bad;
  73.     BOOLEAN error = FALSE, outline = FALSE, surface = FALSE, pen_not_set;
  74.     BOOLEAN first, begin, draw_contours = FALSE, moved[4], intens = FALSE, drape = FALSE;
  75.     BOOLEAN saddle, cross_ok, go, image = FALSE, monochrome = FALSE, simple_clip = TRUE;
  76.     
  77.     double cval, x_left, x_right, y_top, y_bottom, small, z_ave, this_intensity;
  78.     double plane_level = 0.0, dx2, dy2, del_x, del_y, take_out, get_intensity();
  79.     double west = 0.0, east = 0.0, south = 0.0, north = 0.0, new_z_level = 0.0, del_off;
  80.     double data_west, data_east, data_south, data_north, delx, dely, z_val;
  81.     
  82.     float was[4];
  83.     
  84.     char *topofile, *intensfile, *drapefile, *cpt_file;
  85.     
  86.     struct CONT *get_cont_struct(), *start_cont, *this_cont, *last_cont;
  87.     struct POINT *get_point(), *this_point, *last_point;
  88.     
  89.     struct GRD_HEADER header, t_head, i_head, d_head;
  90.     
  91.     struct PEN cpen, mpen;
  92.     struct FILL fill;
  93.     
  94.     cpt_file = topofile = intensfile = drapefile = CNULL;
  95.     
  96.     argc = gmt_begin (argc, argv);
  97.     
  98.     gmt_init_pen (&mpen, 1);
  99.     gmt_init_pen (&cpen, 3);
  100.     gmt_init_fill (&fill, 255 - gmtdefs.basemap_frame_rgb[0], 255 - gmtdefs.basemap_frame_rgb[0], 255 - gmtdefs.basemap_frame_rgb[0]);
  101.     
  102.     for (i = 1; i < argc; i++) {
  103.         if (argv[i][0] == '-') {
  104.             switch (argv[i][1]) {
  105.         
  106.                 /* Common parameters */
  107.             
  108.                 case 'B':
  109.                 case 'J':
  110.                 case 'K':
  111.                 case 'O':
  112.                 case 'P':
  113.                 case 'R':
  114.                 case 'U':
  115.                 case 'V':
  116.                 case 'X':
  117.                 case 'x':
  118.                 case 'Y':
  119.                 case 'y':
  120.                 case 'c':
  121.                 case '\0':
  122.                     error += get_common_args (argv[i], &west, &east, &south, &north);
  123.                     break;
  124.                 
  125.                 /* Supplemental parameters */
  126.             
  127.                 case 'C':
  128.                     cpt_file = &argv[i][2];
  129.                     break;
  130.                 case 'D':
  131.                     n_decimals = atoi (&argv[i][2]);
  132.                     break;
  133.                 case 'E':
  134.                     sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
  135.                     break;
  136.                 case 'F':
  137.                     if (gmt_getrgb (&argv[i][2], &gmtdefs.basemap_frame_rgb[0], &gmtdefs.basemap_frame_rgb[1], &gmtdefs.basemap_frame_rgb[2])) {
  138.                         gmt_pen_syntax ('F');
  139.                         error++;
  140.                     }
  141.                     break;
  142.                 case 'G':
  143.                     drapefile = &argv[i][2];
  144.                     drape = TRUE;
  145.                     break;
  146.                 case 'I':
  147.                     intensfile = &argv[i][2];
  148.                     intens = TRUE;
  149.                     break;
  150.                 case 'M':
  151.                     if (gmt_getpen (&argv[i][2], &mpen)) {
  152.                         gmt_pen_syntax ('M');
  153.                         error++;
  154.                     }
  155.                     break;
  156.                 case 'N':
  157.                     sscanf (&argv[i][2], "%lf/%d/%d/%d", &plane_level, &red_facade, &green_facade, &blue_facade);
  158.                     draw_plane = TRUE;
  159.                     if (strchr (&argv[i][2], '/')) facade = TRUE;
  160.                     if (check_rgb (red_facade, green_facade, blue_facade)) {
  161.                         fprintf (stderr, "%s: GMT SYNTAX ERROR option -N:  RGB components must all be in 0-255 range\n", gmt_program);
  162.                         error = TRUE;
  163.                     }
  164.                     break;
  165.                 case 'Q':
  166.                     switch (argv[i][2]) {
  167.                         case 'm':    /* Mesh plot */
  168.                             mesh = TRUE;
  169.                             if (argv[i][3] == '/' && gmt_getfill (&argv[i][4], &fill)) {
  170.                                 fprintf (stderr, "%s: GMT SYNTAX ERROR -Qm option: RGB components must all be in 0-255 range\n", gmt_program);
  171.                                 error = TRUE;
  172.                             }
  173.                             break;
  174.                         case 's':    /* Color wo/ contours */
  175.                             surface = TRUE;
  176.                             if (argv[i][3] == 'm') outline = TRUE;
  177.                             break;
  178.                         case 'I':    /* image w/o simple clip*/
  179.                             simple_clip = FALSE;
  180.                         case 'i':    /* image w/ simple clip*/
  181.                             image = TRUE;
  182.                             if (argv[i][3]) dpi_i = atoi (&argv[i][3]);
  183.                             break;
  184.                         default:
  185.                             fprintf (stderr, "%s: GMT SYNTAX ERROR:  Unrecognized qualifier (%c) for option -%c\n", gmt_program, argv[i][2], argv[i][1]);
  186.                             error = TRUE;
  187.                             break;
  188.                     }
  189.                     monochrome = (argv[i][strlen(argv[i])-1] == 'g');
  190.                     break;
  191.                 case 'S':
  192.                     smooth = atoi (&argv[i][2]);
  193.                     break;
  194.                 case 'W':
  195.                     if (gmt_getpen (&argv[i][2], &cpen)) {
  196.                         gmt_pen_syntax ('W');
  197.                         error++;
  198.                     }
  199.                     draw_contours = TRUE;
  200.                     break;
  201.                 case 'Z':
  202.                     new_z_level = atof (&argv[i][2]);
  203.                     break;
  204.                     
  205.                 /* Illegal options */
  206.             
  207.                 default:
  208.                     error = TRUE;
  209.                     gmt_default_error (argv[i][1]);
  210.                     break;
  211.             }
  212.         }
  213.         else
  214.             topofile = argv[i];
  215.     }
  216.  
  217.     if (!(mesh || surface || image)) mesh = TRUE;    /* Default is mesh plot */
  218.     
  219.     if (argc == 1 || gmt_quick) {
  220.         fprintf (stderr,"grdview %s - Plot topofiles in 3-D\n\n", GMT_VERSION);
  221.         fprintf (stderr,"usage: grdview <topofile> -J<params> [-B<tickinfo>] [-C<cpt_file>] [-D<ndec>] [-E<v_az/v_el>]\n");
  222.         fprintf (stderr,"    [-F<r/g/b>] [-G<drapefile>] [-I<intensfile>] [-K] [-M<pen>] [-N<level>] [-O] [-P]\n");
  223.         fprintf (stderr,"    [-Q<type>] [-R<region>] [-S<smooth>] [-U] [-V] [-W<pen>]\n");
  224.         fprintf (stderr,"    [-X<x_shift>] [-Y<y_shift>] [-Z<z_level>] [-c<ncopies>]\n\n");
  225.         
  226.         if (gmt_quick) exit (-1);
  227.         
  228.         fprintf (stderr,"    <topofile> is data set to be plotted\n");
  229.         explain_option ('j');
  230.         fprintf (stderr,"\n\tOPTIONS:\n");
  231.         explain_option ('b');
  232.         fprintf (stderr,"    -C Color palette file\n");
  233.         fprintf (stderr,"    -D Number of decimals for color/gray specification in PostScript output [3]\n");
  234.         fprintf (stderr,"    -E Draw perspective from viewpoint at azimuth <v_az>, elevation <v_el> (degrees).\n");
  235.         fprintf (stderr,"       Default is looking straight down [180/90].\n");
  236.         fprintf (stderr, "    -F Set color used for Frame and anotation [%d/%d/%d]\n",
  237.             gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  238.         fprintf (stderr,"    -G <drapefile> rather than <topofile> is the data set to color-code.\n");
  239.         fprintf (stderr,"       Use <topofile> as the relief and \'drape\' the image on top.\n");
  240.         fprintf (stderr,"       Note that -Jz and -N always refers to the <topofile>\n");
  241.         fprintf (stderr,"    -I gives name of intensity file and selects illumination\n");
  242.         explain_option ('Z');
  243.         explain_option ('K');
  244.         fprintf (stderr,"    -M sets mesh pen attributes. [width = %d, color = (%d/%d/%d), solid line].\n", 
  245.             mpen.width, mpen.r, mpen.g, mpen.b);
  246.         fprintf (stderr,"    -N<level> will draw a plane at z = level.  Append color [/r/g/b] to paint\n");
  247.         fprintf (stderr,"       the facade between the plane and the data perimeter.\n");
  248.         explain_option ('O');
  249.         explain_option ('P');
  250.         fprintf (stderr,"    -Q sets plot reQuest. Choose one of the following:\n");
  251.         fprintf (stderr,"       -Qm for Mesh plot [Default].  Append color for mesh paint [%d/%d/%d]\n", fill.r, fill.g, fill.b);
  252.         fprintf (stderr,"       -Qs[m] for colored or shaded Surface.  Append m to draw meshlines on the surface.\n");
  253.         fprintf (stderr,"       -Qi for scanline converting polygons to rasterimage.  Append effective dpi [100].\n");
  254.         fprintf (stderr,"          -QI will Ignore clipping [Default will try to minimize damage to background by using naive clippath]\n");
  255.         fprintf (stderr,"      To force a monochrome image using the YIQ transformation, append g\n");
  256.         explain_option ('R');
  257.         fprintf (stderr,"    -S will smooth contours first (see grdcontour for <smooth> value info).\n");
  258.         explain_option ('U');
  259.         explain_option ('V');
  260.         fprintf (stderr,"    -W draw contours on top of surface or mesh.  [Default is no contours]\n");
  261.         fprintf (stderr, "       Also sets contour line pen attributes [width = %d, color = (%d/%d/%d), solid line].\n", 
  262.             cpen.width, cpen.r, cpen.g, cpen.b);
  263.         explain_option ('X');
  264.         fprintf (stderr, "      -Z For 3-D plots: Set the z-level of map [0]\n");
  265.         explain_option ('c');
  266.         explain_option ('.');
  267.         exit (-1);
  268.     }
  269.         
  270.     if (!topofile) {
  271.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", gmt_program);
  272.         error++;
  273.     }
  274.     if (drape && !drapefile) {
  275.         fprintf (stderr, "%s: GMT SYNTAX ERROR option -G:  Must specify drape file\n", gmt_program);
  276.         error++;
  277.     }
  278.     if (intens && !intensfile) {
  279.         fprintf (stderr, "%s: GMT SYNTAX ERROR option -I:  Must specify intensity file\n", gmt_program);
  280.         error++;
  281.     }
  282.     if ((surface || image || draw_contours) && !cpt_file) {
  283.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify color palette table\n", gmt_program);
  284.         error++;
  285.     }
  286.     if (image && dpi_i <= 0) {
  287.         fprintf (stderr, "%s: GMT SYNTAX ERROR -Qi option:  Must specify positive dpi\n", gmt_program);
  288.         error++;
  289.     }
  290.     if (n_decimals < 1 || n_decimals > 3) {
  291.         fprintf (stderr, "%s: GMT SYNTAX ERROR -D option:  Enter number of decimals as 1, 2, or 3\n", gmt_program);
  292.         error++;
  293.     }
  294.     if (smooth <= 0) {
  295.         fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  smooth value must be positive\n", gmt_program);
  296.         error++;
  297.     }
  298.     if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) {
  299.         fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  Enter azimuth in 0-360 range, elevation in 0-90 range\n", gmt_program);
  300.         error++;
  301.     }
  302.     if (error) exit (-1);
  303.     
  304.     if (cpt_file) {
  305.         read_cpt (cpt_file);
  306.         if (gmt_b_and_w) monochrome = TRUE;
  307.     }
  308.     
  309.     get_contours = ( (mesh && draw_contours) || (surface && gmt_n_colors > 1) );
  310.     
  311.     two = (drape && get_contours) ? 2 : 0;    /* Must read topofile with 2 boundary columns/rows */
  312.     
  313.     if (!strcmp (topofile, "=")) {
  314.         fprintf (stderr, "grdview: Piping of topofile not supported!\n");
  315.         exit (-1);
  316.     }
  317.  
  318.     if (read_grd_info (topofile, &t_head)) {
  319.         fprintf (stderr, "grdview: Error opening file %s\n", topofile);
  320.         exit (-1);
  321.     }
  322.     
  323.     if (intens && read_grd_info (intensfile, &i_head)) {
  324.         fprintf (stderr, "grdview: Error opening file %s\n", intensfile);
  325.         exit (-1);
  326.     }
  327.     
  328.     if (drape && read_grd_info (drapefile, &d_head)) {
  329.         fprintf (stderr, "grdview: Error opening file %s\n", drapefile);
  330.         exit (-1);
  331.     }
  332.     
  333.     off = (t_head.node_offset) ? 0 : 1;
  334.     del_off = (t_head.node_offset) ? 0.5 : 0.0;
  335.     
  336.     
  337.     /* Determine what wesn to pass to map_setup */
  338.  
  339.     if (west == east && south == north) {
  340.         west = header.x_min;
  341.         east = header.x_max;
  342.         south = header.y_min;
  343.         north = header.y_max;
  344.     }
  345.  
  346.     map_setup (west, east, south, north);
  347.  
  348.     /* Determine the wesn to be used to read the grdfile */
  349.  
  350.     grd_setregion (&t_head, &data_west, &data_east, &data_south, &data_north);
  351.  
  352.     /* Read data */
  353.  
  354.     pad[0] = pad[1] = pad[2] = pad[3] = two;
  355.     dummy[0] = dummy[1] = dummy[2] = dummy[3] = 0;
  356.     nx = rint ( (data_east - data_west) / t_head.x_inc) + off;
  357.     ny = rint ( (data_north - data_south) / t_head.y_inc) + off;
  358.     mx = nx + 2 * two;
  359.     my = ny + 2 * two;
  360.     nm = nx * ny;
  361.     topo = (float *) memory (CNULL, mx * my, sizeof (float), "grdview");
  362.     if (read_grd (topofile, &t_head, topo, data_west, data_east, data_south, data_north, pad, FALSE)) {
  363.         fprintf (stderr, "grdview: Error reading file %s\n", topofile);
  364.         exit (-1);
  365.     }
  366.  
  367.     nx_f = t_head.nx;
  368.     ny_f = t_head.ny;
  369.  
  370.     if (gmtdefs.verbose) fprintf (stderr, "grdview: Processing topofile\n");
  371.     
  372.     delx = (t_head.node_offset) ? 0.5 * t_head.x_inc :0.0;
  373.     dely = (t_head.node_offset) ? 0.5 * t_head.y_inc :0.0;
  374.     
  375.     if (drape) {    /* draping wanted */
  376.     
  377.         if (d_head.nx != nx_f || d_head.ny != ny_f) {
  378.             fprintf (stderr, "grdview: Drape file has improper dimensions!\n");
  379.             exit (-1);
  380.         }
  381.         grd = (float *) memory (CNULL, nm, sizeof (float), "grdview");
  382.         if (read_grd (drapefile, &d_head, grd, data_west, data_east, data_south, data_north, dummy, FALSE)) {
  383.             fprintf (stderr, "grdview: Error reading file %s\n", drapefile);
  384.             exit (-1);
  385.         }
  386.         header = d_head;
  387.     }
  388.     else {
  389.         grd = topo;
  390.         header = t_head;
  391.     }
  392.     
  393.     if (project_info.z_bottom == 0.0 && project_info.z_top == 0.0) {
  394.         project_info.z_bottom = t_head.z_min;
  395.         project_info.z_top = t_head.z_max;
  396.         if (draw_plane && plane_level < project_info.z_bottom) project_info.z_bottom = plane_level;
  397.         if (draw_plane && plane_level > project_info.z_top) project_info.z_top = plane_level;
  398.     }
  399.     
  400.     if (!project_info.xyz_pos[2])    /* Negative z-scale, must flip */
  401.         d_swap (project_info.z_bottom, project_info.z_top);
  402.  
  403.     ij_inc[0] = 0;        ij_inc[1] = 1;    ij_inc[2] = 1 - mx;    ij_inc[3] = -mx;
  404.     bin_inc[0] = 0;        bin_inc[1] = 1;    bin_inc[2] = 1 - header.nx;    bin_inc[3] = -header.nx;
  405.     nw = two * mx + two;
  406.     ne = nw + t_head.nx - 1;
  407.     sw = (t_head.ny + two - 1) * mx + two;
  408.     se = sw + t_head.nx - 1;
  409.     
  410.     init_setup (&t_head, topo, two, draw_plane, plane_level);    /* Find projected min/max in y-direction */
  411.             
  412.     i_start = (z_project.quadrant == 1 || z_project.quadrant == 2) ? 0 : header.nx - 2;
  413.     i_stop  = (z_project.quadrant == 1 || z_project.quadrant == 2) ? header.nx - 1 : -1;
  414.     i_inc   = (z_project.quadrant == 1 || z_project.quadrant == 2) ? 1 : -1;
  415.     j_start = (z_project.quadrant == 1 || z_project.quadrant == 4) ? header.ny - 1 : 1;
  416.     j_stop  = (z_project.quadrant == 1 || z_project.quadrant == 4) ? 0 : header.ny;
  417.     j_inc   = (z_project.quadrant == 1 || z_project.quadrant == 4) ? -1 : 1;
  418.     x_inc[0] = x_inc[3] = 0.0;    x_inc[1] = x_inc[2] = header.x_inc;
  419.     y_inc[0] = y_inc[1] = 0.0;    y_inc[2] = y_inc[3] = header.y_inc;
  420.  
  421.     if (get_contours) {    /* Need to find contours */
  422.         if (gmtdefs.verbose) fprintf (stderr, "grdview: Find contours\n");
  423.         n_edges = header.ny * (int )ceil (header.nx / 16.0);
  424.         edge = (int *) memory (CNULL, n_edges, sizeof (int), "grdview");
  425.         binij = (struct BIN *) memory (CNULL, nm, sizeof (struct BIN), "grdview");
  426.         small = SMALL * (header.z_max - header.z_min);
  427.         if (gmtdefs.verbose) fprintf (stderr, "grdview: Trace and bin contours...\n");
  428.         first = TRUE;
  429.         for (c = 0; c < gmt_n_colors+1; c++) {    /* For each color change */
  430.         
  431.             /* Reset markers and set up new zero-contour*/
  432.         
  433.             cval = (c == gmt_n_colors) ? gmt_lut[c-1].z_high : gmt_lut[c].z_low;
  434.             if (gmtdefs.verbose) fprintf (stderr, "grdview: Now tracing contour interval %8lg\r", cval);
  435.             take_out = (first) ? cval : cval - gmt_lut[c-1].z_low;
  436.             first = FALSE;
  437.             for (i = 0; i < nm; i++) {
  438.                 if (!bad_float ((double)grd[i])) grd[i] -= take_out;
  439.                 if (grd[i] == 0.0) grd[i] += small;
  440.             }
  441.         
  442.             side = 0;
  443.             begin = TRUE;
  444.             while (side < 5) {
  445.                 while ((n = contours (grd, &header, smooth, gmtdefs.interpolant, &side, edge, begin, &x, &y)) > 0) {
  446.         
  447.                     begin = FALSE;
  448.         
  449.                     i_bin_old = j_bin_old = -1;
  450.                     for (i = 1; i < n; i++) {
  451.                         i_bin = floor (((0.5 * (x[i-1] + x[i]) - header.x_min) / header.x_inc) - del_off);
  452.                         j_bin = floor (((header.y_max - 0.5 * (y[i-1] + y[i])) / header.y_inc) - del_off) + 1;
  453.                         if (i_bin != i_bin_old || j_bin != j_bin_old) {    /* Entering new bin */
  454.                             bin = j_bin * header.nx + i_bin;
  455.                             this_cont = get_cont_struct (bin, cval);
  456.                             this_cont->value = cval;
  457.                             this_cont->first_point = get_point (x[i-1], y[i-1]);
  458.                             this_point = this_cont->first_point;
  459.                             i_bin_old = i_bin;
  460.                             j_bin_old = j_bin;
  461.                         }
  462.                         this_point->next_point = get_point (x[i], y[i]);
  463.                         this_point = this_point->next_point;
  464.                     }
  465.                     free ((char *)x);
  466.                     free ((char *)y);
  467.                 }
  468.                 begin = FALSE;
  469.             }
  470.         }
  471.     
  472.         /* Remove temporary variables */
  473.         
  474.         free ((char*)edge);
  475.         
  476.         /* Go back to beginning and reread since grd has been destroyed */
  477.         
  478.         if (drape) {
  479.             if (read_grd_info (drapefile, &d_head)) {
  480.                 fprintf (stderr, "grdview: Error opening file %s\n", drapefile);
  481.                 exit (-1);
  482.             }
  483.             if (read_grd (drapefile, &d_head, grd, data_west, data_east, data_south, data_north, dummy, FALSE)) {
  484.                 fprintf (stderr, "grdview: Error reading file %s\n", drapefile);
  485.                 exit (-1);
  486.             }
  487.             header = d_head;
  488.         }
  489.         else {
  490.             if (read_grd_info (topofile, &t_head)) {
  491.                 fprintf (stderr, "grdview: Error opening file %s\n", topofile);
  492.                 exit (-1);
  493.             }
  494.             if (read_grd (topofile, &t_head, topo, data_west, data_east, data_south, data_north, pad, FALSE)) {
  495.                 fprintf (stderr, "grdview: Error reading file %s\n", topofile);
  496.                 exit (-1);
  497.             }
  498.             grd = topo;
  499.             header = t_head;
  500.         }
  501.         if (gmtdefs.verbose) fprintf (stderr, "\n");
  502.     }
  503.     
  504.     if (intens) {    /* Illumination wanted */
  505.     
  506.         if (i_head.nx != nx_f || i_head.ny != ny_f) {
  507.             fprintf (stderr, "grdview: Intensity file has improper dimensions!\n");
  508.             exit (-1);
  509.         }
  510.         intensity = (float *) memory (CNULL, nm, sizeof (float), "grdview");
  511.         if (read_grd (intensfile, &i_head, intensity, data_west, data_east, data_south, data_north, dummy, FALSE)) {
  512.             fprintf (stderr, "grdview: Error reading file %s\n", intensfile);
  513.             exit (-1);
  514.         }
  515.     }
  516.     
  517.     if (two) {    /* Initialize bcr stuff */
  518.         bcr_init (&t_head, pad, FALSE);
  519.  
  520.         /* Set boundary conditions  */
  521.     
  522.         set_bcr_boundaries (&t_head, topo);
  523.     }
  524.  
  525.     max = 2 * (smooth * (((header.nx > header.ny) ? header.nx : header.ny) + 2)) + 1;
  526.     x = (double *) memory (CNULL, max, sizeof (double), "grdview");
  527.     y = (double *) memory (CNULL, max, sizeof (double), "grdview");
  528.     z = (double *) memory (CNULL, max, sizeof (double), "grdview");
  529.     v = (double *) memory (CNULL, max, sizeof (double), "grdview");
  530.     
  531.     dx2 = 0.5 * header.x_inc;        dy2 = 0.5 * header.y_inc;
  532.             
  533.     if (gmtdefs.verbose) fprintf (stderr, "grdview: Start creating PostScript plot\n");
  534.  
  535.     ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
  536.         gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
  537.         gmtdefs.dpi, gmtdefs.measure_unit , gmtdefs.paper_width, gmtdefs.page_rgb, gmt_epsinfo (argv[0]));
  538.  
  539.     ps_setformat (n_decimals);
  540.     echo_command (argc, argv);
  541.     if (gmtdefs.unix_time) timestamp (argc, argv);
  542.  
  543.     if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0);
  544.     if (new_z_level != 0.0) project_info.z_level = new_z_level;
  545.     
  546.     if (frame_info.plot && !image) {    /* Plot basemap */
  547.         frame_info.header[0] = 0;    /* No header for grdview and psxyz */
  548.         ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  549.         map_basemap ();
  550.         ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  551.     }
  552.  
  553.     max = 2 * (smooth * (((header.nx > header.ny) ? header.nx : header.ny) + 2)) + 1;
  554.     xx = (double *) memory (CNULL, max, sizeof(double), "grdview");
  555.     yy = (double *) memory (CNULL, max, sizeof(double), "grdview");
  556.     if (draw_plane) {
  557.         ps_comment ("Plot the plane at desired level");
  558.         if (!z_project.draw[0])    {    /* Southern side */
  559.             for (i = 0; i < header.nx; i++) geoz_to_xy (header.x_min + i * header.x_inc + delx, header.y_min + dely, plane_level, &xx[i], &yy[i]);
  560.             ps_line (xx, yy, header.nx, 3, TRUE, TRUE);
  561.         }
  562.         if (!z_project.draw[2])    {    /* Northern side */
  563.             for (i = 0; i < header.nx; i++) geoz_to_xy (header.x_min + i * header.x_inc + delx, header.y_max - dely, plane_level, &xx[i], &yy[i]);
  564.             ps_line (xx, yy, header.nx, 3, TRUE, TRUE);
  565.         }
  566.         if (!z_project.draw[3])    {    /* Western side */
  567.             for (j = 0; j < header.ny; j++) geoz_to_xy (header.x_min + delx, header.y_max - j * header.y_inc - dely, plane_level, &xx[j], &yy[j]);
  568.             ps_line (xx, yy, header.ny, 3, TRUE, TRUE);
  569.         }
  570.         if (!z_project.draw[1])    {    /* Eastern side */
  571.             for (j = 0; j < header.ny; j++) geoz_to_xy (header.x_max - delx, header.y_max - j * header.y_inc - dely, plane_level, &xx[j], &yy[j]);
  572.             ps_line (xx, yy, header.ny, 3, TRUE, TRUE);
  573.         }
  574.         
  575.         geoz_to_xy (header.x_min + delx, header.y_min + dely, plane_level, &xx[0], &yy[0]);
  576.         geoz_to_xy (header.x_max - delx, header.y_min + dely, plane_level, &xx[1], &yy[1]);
  577.         geoz_to_xy (header.x_max - delx, header.y_max - dely, plane_level, &xx[2], &yy[2]);
  578.         geoz_to_xy (header.x_min + delx, header.y_max - dely, plane_level, &xx[3], &yy[3]);
  579.         if (!bad_float ((double)topo[nw])) {
  580.             geoz_to_xy (header.x_min + delx, header.y_max - dely, (double)(topo[nw]), &x_left, &y_top);
  581.             ps_plot (x_left, y_top, 3);    ps_plot (xx[3], yy[3], 2);
  582.         }
  583.         if (!bad_float ((double)topo[ne])) {
  584.             geoz_to_xy (header.x_max - delx, header.y_max - dely, (double)(topo[ne]), &x_right, &y_top);
  585.             ps_plot (x_right, y_top, 3);    ps_plot (xx[2], yy[2], 2);
  586.         }
  587.         if (!bad_float ((double)topo[se])) {
  588.             geoz_to_xy (header.x_max - delx, header.y_min + dely, (double)(topo[se]), &x_right, &y_bottom);
  589.             ps_plot (x_right, y_bottom, 3);    ps_plot (xx[1], yy[1], 2);
  590.         }
  591.         if (!bad_float ((double)topo[sw])) {
  592.             geoz_to_xy (header.x_min + delx, header.y_min + dely, (double)(topo[sw]), &x_left, &y_bottom);
  593.             ps_plot (x_left, y_bottom, 3);    ps_plot (xx[0], yy[0], 2);
  594.         }
  595.         ps_plot (0.0, 0.0, 3);    /* To make sure path is stroked */
  596.     }
  597.     
  598.     if (image) {    /* compute image */
  599.         int nx_i, ny_i, kk, ip, jp, min_i, max_i, min_j, max_j, dist, node, nm_i, layers, *ix, *iy;
  600.         BOOLEAN done;
  601.         double *xval, *yval, xp, yp, sum_w, w, sum_i, x_width, y_width;
  602.         double sum_r, sum_g, sum_b, intval;
  603.         unsigned char *bitimage_24, *bitimage_8;
  604.         
  605.         if (gmtdefs.verbose) fprintf (stderr, "grdview: get and store projected vertices\n");
  606.         ps_comment ("Plot 3-D surface using scanline conversion of polygons to raster image");
  607.  
  608.         x_width = z_project.xmax - z_project.xmin;    /* Size of image in inches */
  609.         y_width = z_project.ymax - z_project.ymin;
  610.         nx_i = rint (x_width * dpi_i);    /* Size of image in pixels */
  611.         ny_i = rint (y_width * dpi_i);
  612.         
  613.         ix = (int *) memory (CNULL, nm, sizeof (int), "grdview");
  614.         iy = (int *) memory (CNULL, nm, sizeof (int), "grdview");
  615.         
  616.         xval = (double *) memory (CNULL, header.nx, sizeof (double), "grdview");
  617.         yval = (double *) memory (CNULL, header.ny, sizeof (double), "grdview");
  618.         for (i = 0; i < header.nx; i++, bin++) xval[i] = header.x_min + i * header.x_inc;
  619.         for (j = 0; j < header.ny; j++) yval[j] = header.y_max - j * header.y_inc;
  620.         
  621.         for (j = bin = ij = 0; j < header.ny; j++) {    /* Get projected coordinates converted to pixel locations */
  622.             ij = (two) ? (j + 2) * mx + two : bin;
  623.             for (i = 0; i < header.nx; i++, bin++, ij++) {
  624.                 geoz_to_xy (xval[i], yval[j], (double)topo[ij], &xp, &yp);
  625.                 ix[bin] = floor ((xp - z_project.xmin) * dpi_i);
  626.                 iy[bin] = floor ((yp - z_project.ymin) * dpi_i);
  627.             }
  628.         }
  629.         free ((char *) xval);
  630.         free ((char *) yval);
  631.         
  632.         /* Allocate image array and set background to PAGE_COLOR */
  633.         
  634.         if (monochrome) {
  635.             int gray;
  636.             
  637.             nm_i = nx_i * ny_i;
  638.             layers = 1;
  639.             bitimage_8 = (unsigned char *) memory (CNULL, nm_i, sizeof (char), "grdview");
  640.             gray = YIQ (gmtdefs.page_rgb[0], gmtdefs.page_rgb[1], gmtdefs.page_rgb[2]);
  641.             memset ((char *)bitimage_8, gray, nm_i);
  642.         }
  643.         else {
  644.             nm_i = nx_i * ny_i * 3;
  645.             layers = 3;
  646.             bitimage_24 = (unsigned char *) memory (CNULL, nm_i, sizeof (char), "grdview");
  647.             kk = 0;
  648.             while (kk < nm_i) {
  649.                 bitimage_24[kk++] = gmtdefs.page_rgb[0];
  650.                 bitimage_24[kk++] = gmtdefs.page_rgb[1];
  651.                 bitimage_24[kk++] = gmtdefs.page_rgb[2];
  652.             }
  653.         }
  654.                 
  655.         /* Plot from back to front */
  656.         
  657.         if (gmtdefs.verbose) fprintf (stderr, "grdview: Start rasterization\n");
  658.         for (j = j_start; j != j_stop; j += j_inc) {
  659. #ifdef DEBUG
  660.             fprintf (stderr, "grdview: Scan line conversion at j-line %.4d\r", j);
  661. #endif
  662.             for (i = i_start; i != i_stop; i += i_inc) {
  663.                 bin = j * header.nx + i;
  664.                 ij = (two) ? (j + 2) * header.nx + i + 2 : bin;
  665.                 for (k = bad = 0; !bad && k < 4; k++) bad += bad_float ((double)topo[ij+ij_inc[k]]);
  666.                 if (bad) continue;
  667.                 min_i = MIN (ix[bin+bin_inc[0]], MIN (ix[bin+bin_inc[1]], MIN (ix[bin+bin_inc[2]], ix[bin+bin_inc[3]])));
  668.                 max_i = MAX (ix[bin+bin_inc[0]], MAX (ix[bin+bin_inc[1]], MAX (ix[bin+bin_inc[2]], ix[bin+bin_inc[3]])));
  669.                 min_j = MIN (iy[bin+bin_inc[0]], MIN (iy[bin+bin_inc[1]], MIN (iy[bin+bin_inc[2]], iy[bin+bin_inc[3]])));
  670.                 max_j = MAX (iy[bin+bin_inc[0]], MAX (iy[bin+bin_inc[1]], MAX (iy[bin+bin_inc[2]], iy[bin+bin_inc[3]])));
  671.                 
  672.                 for (jp = min_j; jp < max_j; jp++) {
  673.                     for (ip = min_i; ip < max_i; ip++) {
  674.                         if (!pixel_inside (ip, jp, ix, iy, bin)) continue;
  675.                         kk = layers * ((ny_i-jp-1) * nx_i + ip);
  676.                         sum_r = sum_g = sum_b = sum_w = sum_i = 0.0;
  677.                         done = FALSE;
  678.                         for (k = 0; !done && k < 4; k++) {
  679.                             node = bin + bin_inc[k];
  680.                             get_rgb24 (grd[node], &r, &g, &b);
  681.                             dist = idist (ip, jp, ix[node], iy[node]);
  682.                             if (dist == 0) {    /* Only need this corner value */
  683.                                 done = TRUE;
  684.                                 if (intens) intval = intensity[node];
  685.                             }
  686.                             else {
  687.                                 w = 1.0 / (double)dist;
  688.                                 sum_r += r * w;
  689.                                 sum_g += g * w;
  690.                                 sum_b += b * w;
  691.                                 if (intens) sum_i += intensity[node] * w;
  692.                                 sum_w += w;
  693.                             }
  694.                         }
  695.                         if (!done) {    /* Must get weighted value */
  696.                             sum_w = 1.0 / sum_w;
  697.                             r = rint (sum_r * sum_w);
  698.                             g = rint (sum_g * sum_w);
  699.                             b = rint (sum_b * sum_w);
  700.                             if (intens) intval = sum_i * sum_w;
  701.                         }
  702.                         if (intens) illuminate (intval, &r, &g, &b);
  703.                         if (monochrome) /* YIQ transformation */
  704.                             bitimage_8[kk] = (unsigned char) YIQ (r, g, b);
  705.                         else {
  706.                             bitimage_24[kk++] = r;
  707.                             bitimage_24[kk++] = g;
  708.                             bitimage_24[kk] = b;
  709.                         }
  710.                     }
  711.                 }
  712.             }
  713.         }
  714.         
  715.         if (!project_info.three_D)
  716.             map_clip_on (-1, -1, -1, 3);
  717.         else if (simple_clip) {    /* Will attempt a simple clip path to minimize damage to background plots */
  718.             int jm1, jp1;
  719.             
  720.             geoz_to_xy (header.x_min, header.y_min, project_info.z_bottom, &xx[0], &yy[0]);
  721.             geoz_to_xy (header.x_max, header.y_min, project_info.z_bottom, &xx[1], &yy[1]);
  722.             geoz_to_xy (header.x_max, header.y_max, project_info.z_bottom, &xx[2], &yy[2]);
  723.             geoz_to_xy (header.x_min, header.y_max, project_info.z_bottom, &xx[3], &yy[3]);
  724.         
  725.             for (i = 1, j = 0; i < 4; i++) if (yy[i] > yy[j]) j = i;
  726.             for (i = 3; i >= j+1; i--) xx[i+1] = xx[i], yy[i+1] = yy[i];
  727.             jm1 = j - 1;
  728.             if (jm1 < 0) jm1 += 5;
  729.             jp1 = j + 1;
  730.             if (jp1 > 3) jp1 -= 4;
  731.             xx[j] = xx[jm1];
  732.             yy[j] = z_project.ymax;        
  733.             xx[j+1] = xx[jp1];
  734.             yy[j+1] = z_project.ymax;
  735.             ps_clipon (xx, yy, 5, -1, -1, -1, 3);
  736.         }
  737.         
  738.         if (gmtdefs.verbose) fprintf (stderr, "grdview: Creating PostScript image ");
  739.         if (monochrome) {
  740.             if (gmtdefs.verbose) fprintf (stderr, "[B/W image]\n");
  741.             ps_image (z_project.xmin, z_project.ymin, x_width, y_width, bitimage_8, nx_i, ny_i, 8);
  742.             free ((char *)bitimage_8);
  743.         }
  744.         else {
  745.             if (gmtdefs.verbose) fprintf (stderr, "[%s]\n", c_method[gmtdefs.color_image]);
  746.             if (gmtdefs.color_image == 2)
  747.                 fprintf (stderr, "%s: GMT Warning: PostScript output may only be processed with psto24!\n", gmt_program);
  748.             color_image (z_project.xmin, z_project.ymin, x_width, y_width, bitimage_24, nx_i, ny_i);
  749.             free ((char *)bitimage_24);
  750.         }
  751.         if (!project_info.three_D || simple_clip) ps_clipoff ();
  752.         
  753.         free ((char *)ix);
  754.         free ((char *)iy);
  755.     }
  756.     if (mesh) {
  757.         ps_comment ("Start of mesh plot");
  758.         ps_setline (mpen.width);
  759.         ps_setpaint (mpen.r, mpen.g, mpen.b);
  760.         if (mpen.texture) ps_setdash (mpen.texture, mpen.offset);
  761.         for (j = j_start; j != j_stop; j += j_inc) {
  762.             y_bottom = header.y_max - j * header.y_inc - dely;
  763.             y_top = y_bottom + abs (j_inc) * header.y_inc;
  764.             for (i = i_start; i != i_stop; i += i_inc) {
  765.                 bin = j * header.nx + i;
  766.                 ij = (two) ? (j + 2) * mx + i + 2 : bin;
  767.                 for (k = bad = 0; !bad && k < 4; k++) bad += bad_float ((double)topo[ij+ij_inc[k]]);
  768.                 if (bad) continue;
  769.                 x_left = header.x_min + i * header.x_inc + delx;
  770.                 x_right = x_left + abs (i_inc) * header.x_inc;
  771.                 geoz_to_xy (x_left, y_bottom, (double)(topo[ij+ij_inc[0]]), &xx[0], &yy[0]);
  772.                 geoz_to_xy (x_right, y_bottom, (double)(topo[ij+ij_inc[1]]), &xx[1], &yy[1]);
  773.                 geoz_to_xy (x_right, y_top, (double)(topo[ij+ij_inc[2]]), &xx[2], &yy[2]);
  774.                 geoz_to_xy (x_left, y_top, (double)(topo[ij+ij_inc[3]]), &xx[3], &yy[3]);
  775.                 ps_patch (xx, yy, 4, fill.r, fill.g, fill.b, TRUE);
  776.                 if (draw_contours) {
  777.                     pen_not_set = TRUE;
  778.                     if (binij[bin].first_cont == NULL) continue;
  779.                     for (this_cont = binij[bin].first_cont->next_cont; this_cont; this_cont = this_cont->next_cont) {
  780.                         for (k = 0, this_point = this_cont->first_point; this_point; this_point = this_point->next_point, k++) {
  781.                             z_val = (drape) ? get_bcr_z (&t_head, (double)this_point->x, (double)this_point->y, topo) : this_cont->value;
  782.                             if (z_val < -5000.0)
  783.                                 error++;
  784.                             geoz_to_xy ((double)this_point->x, (double)this_point->y, z_val, &xx[k], &yy[k]);
  785.                         }
  786.                         if (pen_not_set) {
  787.                             if (mpen.texture) ps_setdash (CNULL, 0);
  788.                             ps_setline (cpen.width);
  789.                             ps_setpaint (cpen.r, cpen.g, cpen.b);
  790.                             if (cpen.texture) ps_setdash (cpen.texture, cpen.offset);
  791.                             pen_not_set = FALSE;
  792.                         }
  793.                         ps_line (xx, yy, k, 3, FALSE, TRUE);
  794.                     }
  795.                     if (!pen_not_set) {
  796.                         if (cpen.texture) ps_setdash (CNULL, 0);
  797.                         ps_setline (mpen.width);
  798.                         ps_setpaint (mpen.r, mpen.g, mpen.b);
  799.                         if (mpen.texture) ps_setdash (mpen.texture, mpen.offset);
  800.                     }
  801.                 }                    
  802.             }
  803.         }
  804.         if (mpen.texture) ps_setdash (CNULL, 0);
  805.     }
  806.     else if (surface) {
  807.         ps_comment ("Start of filled surface");
  808.         if (outline) {
  809.             ps_setline (mpen.width);
  810.             ps_setpaint (mpen.r, mpen.g, mpen.b);
  811.             if (mpen.texture) ps_setdash (mpen.texture, mpen.offset);
  812.         }
  813.         
  814.         for (j = j_start; j != j_stop; j += j_inc) {
  815.             y_bottom = header.y_max - j * header.y_inc - dely;
  816.             y_top = y_bottom + header.y_inc;
  817.             for (i = i_start; i != i_stop; i += i_inc) {
  818. #ifdef DEBUG
  819.                 printf ("%% Doing bin (%d,%d)\n", i, j);
  820. #endif
  821.                 bin = j * header.nx + i;
  822.                 ij = (two) ? (j + 2) * mx + i + 2 : bin;
  823.                 for (k = bad = 0; !bad && k < 4; k++) bad += bad_float ((double)topo[ij+ij_inc[k]]);
  824.                 if (bad) continue;
  825.                 x_left = header.x_min + i * header.x_inc + delx;
  826.                 x_right = x_left + header.x_inc;
  827.                 if (intens) {
  828.                     this_intensity = get_intensity (intensity, bin, header.nx);
  829.                     if (bad_float (this_intensity)) continue;
  830.                 }
  831.                 if (get_contours && binij[bin].first_cont) {    /* Contours go thru here */
  832.                     start_cont = binij[bin].first_cont->next_cont;
  833.                     
  834.                     /* Paint the part of the tile that has a value below the contour value */
  835.                     
  836.                     geoz_to_xy (x_left, y_bottom, (double)(topo[ij+ij_inc[0]]), &xx[0], &yy[0]);
  837.                     geoz_to_xy (x_right, y_bottom, (double)(topo[ij+ij_inc[1]]), &xx[1], &yy[1]);
  838.                     geoz_to_xy (x_right, y_top, (double)(topo[ij+ij_inc[2]]), &xx[2], &yy[2]);
  839.                     geoz_to_xy (x_left, y_top, (double)(topo[ij+ij_inc[3]]), &xx[3], &yy[3]);
  840.                     
  841.                     /* Set the right color */
  842.                     
  843.                     if (gmt_continuous) {
  844.                         z_ave = 0.0;
  845.                         for (n = k = 0; n < 4; n++) {
  846.                             if (grd[bin+bin_inc[n]] >= start_cont->value) continue;
  847.                             k++;
  848.                             z_ave += grd[bin+bin_inc[n]];
  849.                         }
  850.                         z_ave /= k;
  851.                     }
  852.                     else
  853.                         z_ave = start_cont->value - small;
  854.         
  855.                     index = get_rgb24 (z_ave, &r, &g, &b);
  856.                     if (index >= 0 && !gmt_lut[index].skip) {
  857.                         if (intens ) illuminate (this_intensity, &r, &g, &b);
  858.                         ps_patch (xx, yy, 4, r, g, b, outline);
  859.                     }
  860.                     
  861.                     /* Now loop over all the contours that crossed this tile */
  862.                     
  863.                     saddle = FALSE;
  864.                     cross_ok = TRUE;
  865.                     n_saddle = 0;
  866.                     
  867.                     for (this_cont = start_cont; this_cont; this_cont = this_cont->next_cont) {
  868.                     
  869.                         /* First get all the x/y pairs for this contour */
  870.                         
  871.                         for (k = 0, this_point = this_cont->first_point; this_point; this_point = this_point->next_point, k++) {
  872.                             x[k] = this_point->x;
  873.                             y[k] = this_point->y;
  874.                             z[k] = (drape) ? get_bcr_z (&t_head, x[k], y[k], topo) : this_cont->value;
  875.                             v[k] = this_cont->value;
  876.                         }
  877.                         
  878.                         /* Check for saddle situation */
  879.                         
  880.                         if (!saddle && this_cont->next_cont && this_cont->next_cont->value == this_cont->value) saddle = TRUE;
  881.                         if (saddle) n_saddle++;
  882.                         
  883.                         /* First make sure no grd-value equals the contour value */
  884.                         
  885.                         
  886.                         for (side = 0; side < 4; side++) {
  887.                             if (grd[bin+bin_inc[side]] == this_cont->value) {
  888.                                 was[side] = grd[bin+bin_inc[side]];
  889.                                 grd[bin+bin_inc[side]] += small;
  890.                                 moved[side] = TRUE;
  891.                             }
  892.                             else
  893.                                 moved[side] = FALSE;
  894.                         }
  895.                         
  896.                         /* Then figure out on what side the contour exits */
  897.                         
  898.                         del_x = fabs (fmod (x[k-1] - header.x_min, header.x_inc));
  899.                         if (del_x > dx2) del_x = header.x_inc - del_x;
  900.                         del_y = fabs (fmod (y[k-1] - header.y_min, header.y_inc));
  901.                         if (del_y > dy2) del_y = header.y_inc - del_y;
  902.                         if (del_x < del_y) /* Cutting N-S gridlines */
  903.                             side = ((x[k-1]-x_left) > dx2) ? 1 : 3;
  904.                         else /* Cutting E-W gridlines */
  905.                             side = ((y[k-1]-y_bottom) > dy2) ? 2 : 0;
  906.                             
  907.                         /* If saddle, figure out on what side the contour enters */
  908.                         
  909.                         if (saddle && n_saddle == 1) {
  910.                             del_x = fabs (fmod (x[0] - header.x_min, header.x_inc));
  911.                             if (del_x > dx2) del_x = header.x_inc - del_x;
  912.                             del_y = fabs (fmod (y[0] - header.y_min, header.y_inc));
  913.                             if (del_y > dy2) del_y = header.y_inc - del_y;
  914.                             if (del_x < del_y) /* Cutting N-S gridlines */
  915.                                 entry = ((x[0]-x_left) > dx2) ? 1 : 3;
  916.                             else /* Cutting E-W gridlines */
  917.                                 entry = ((y[0]-y_bottom) > dy2) ? 2 : 0;
  918.                             
  919.                             if ((entry == 0 && side == 3) || (entry == 3 && side == 0))
  920.                                 check = 0;
  921.                             else if ((entry == 1 && side == 2) || (entry == 2 && side == 1))
  922.                                 check = 2;
  923.                             else if ((entry == 0 && side == 1) || (entry == 1 && side == 0))
  924.                                 check = 1;
  925.                             else
  926.                                 check = 3;
  927.                             
  928.                             cross_ok = (grd[bin+bin_inc[check]] < this_cont->value);
  929.                                 
  930.                         }
  931.                             
  932.                         /* Now determine which way we have to go around the tile
  933.                            outline to close the polygon */
  934.                            
  935.                         way = (grd[bin+bin_inc[side]] >= this_cont->value) ? -1 : 1;
  936.                         
  937.                         /* Add the corner points that are needed to close the polygon */
  938.                         
  939.                         prev_side = side;
  940.                         if (way == 1) side++;
  941.                         for (n = n_added = 0; n < 4; n++) {
  942.                             side %= 4;
  943.                             if (cross_ok)
  944.                                 go = TRUE;
  945.                             else if (n_added == 0)
  946.                                 go = TRUE;
  947.                             else
  948.                                 go = (side != (prev_side + 2)%4);
  949.                             if (go && grd[bin+bin_inc[side]] > this_cont->value) {
  950.                                 x[k] = x_left + x_inc[side];
  951.                                 y[k] = y_bottom + y_inc[side];
  952.                                 z[k] = topo[ij+ij_inc[side]];
  953.                                 v[k] = grd[bin+bin_inc[side]];
  954.                                 k++;
  955.                                 n_added++;
  956.                                 prev_side = side;
  957.                             }
  958.                             side += way;
  959.                             if (side < 0) side += 4;
  960.                         }
  961.                         
  962.                         /* Compute the xy from the xyz triplets */
  963.                         
  964.                         for (n = 0; n < k; n++) geoz_to_xy (x[n], y[n], z[n], &xx[n], &yy[n]);
  965.                         
  966.                         if (gmt_continuous) {
  967.                             z_ave = 0.0;
  968.                             for (n = 0; n < k; n++) z_ave += v[n];
  969.                             z_ave /= k;
  970.                         }
  971.                         else
  972.                             z_ave = this_cont->value;
  973.                         
  974.                         /* Set the right color */
  975.                         
  976.                         get_rgb24 (z_ave+small, &r, &g, &b);
  977.                         if (intens) illuminate (this_intensity, &r, &g, &b);
  978.  
  979.                         ps_patch (xx, yy, k, r, g, b, outline);
  980.                         
  981.                         /* Reset nodes that was moved off contour-value */
  982.                         
  983.                         for (side = 0; side < 4; side++) if (moved[side]) grd[bin+bin_inc[side]] = was[side];
  984.                         
  985.                         if (n_saddle == 2) {
  986.                             n_saddle = 0;
  987.                             saddle = FALSE;
  988.                             cross_ok = TRUE;
  989.                         }
  990.                     }
  991.                     
  992.                     /* Draw contour lines if desired */
  993.                     
  994.                     pen_not_set = TRUE;
  995.                     for (this_cont = start_cont; draw_contours && this_cont; this_cont = this_cont->next_cont) {
  996.                         for (k = 0, this_point = this_cont->first_point; this_point; this_point = this_point->next_point, k++) {
  997.                             z_val = (drape) ? get_bcr_z (&t_head, (double)this_point->x, (double)this_point->y, topo) : this_cont->value;
  998.                             geoz_to_xy ((double)this_point->x, (double)this_point->y, z_val, &xx[k], &yy[k]);
  999.                         }
  1000.                         if (pen_not_set) {
  1001.                             if (mpen.texture) ps_setdash (CNULL, 0);
  1002.                             ps_setline (cpen.width);
  1003.                             ps_setpaint (cpen.r, cpen.g, cpen.b);
  1004.                             if (cpen.texture) ps_setdash (cpen.texture, cpen.offset);
  1005.                             pen_not_set = FALSE;
  1006.                         }
  1007.                         ps_line (xx, yy, k, 3, FALSE, TRUE);
  1008.                     }
  1009.                     if (!pen_not_set) {                        
  1010.                         if (cpen.texture) ps_setdash (CNULL, 0);
  1011.                         ps_setline (mpen.width);
  1012.                         ps_setpaint (mpen.r, mpen.g, mpen.b);
  1013.                         if (mpen.texture) ps_setdash (mpen.texture, mpen.offset);
  1014.                     }
  1015.                 }
  1016.                 else {    /* No Contours, paint tile (any of the corner values will give correct color */
  1017.                     geoz_to_xy (x_left, y_bottom, (double)(topo[ij]), &xx[0], &yy[0]);
  1018.                     geoz_to_xy (x_right, y_bottom, (double)(topo[ij+ij_inc[1]]), &xx[1], &yy[1]);
  1019.                     geoz_to_xy (x_right, y_top, (double)(topo[ij+ij_inc[2]]), &xx[2], &yy[2]);
  1020.                     geoz_to_xy (x_left, y_top, (double)(topo[ij+ij_inc[3]]), &xx[3], &yy[3]);
  1021.                     
  1022.                     if (gmt_continuous) {
  1023.                         for (n = 0, z_ave = 0.0; n < 4; n++) z_ave += grd[bin+bin_inc[n]];
  1024.                         z_ave *= 0.25;
  1025.                     }
  1026.                     else
  1027.                         z_ave = grd[bin];
  1028.                     index = get_rgb24 (z_ave, &r, &g, &b);
  1029.                     if (index >= 0 && !gmt_lut[index].skip) {
  1030.                         if (intens ) illuminate (this_intensity, &r, &g, &b);
  1031.                         ps_patch (xx, yy, 4, r, g, b, outline);
  1032.                     }
  1033.                 }
  1034.             }
  1035.         }
  1036.     }            
  1037.     if (mpen.texture || cpen.texture) ps_setdash (CNULL, 0);
  1038.     
  1039.     if (project_info.three_D) {    /* Redraw front axes */
  1040.         ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  1041.         vertical_axis (2);
  1042.         ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  1043.     }
  1044.     
  1045.     if (facade) {    /* Cover the two front sides */
  1046.         ps_comment ("Paiting the frontal facade");
  1047.         ps_setline (1);
  1048.         if (!z_project.draw[0])    {    /* Southern side */
  1049.             for (i = n = 0, ij = sw; i < header.nx; i++, ij++) {
  1050.                 if (bad_float ((double)topo[ij])) continue;
  1051.                 geoz_to_xy (header.x_min+i*header.x_inc+delx, header.y_min+dely, (double)(topo[ij]), &xx[n], &yy[n]);
  1052.                 n++;
  1053.             }
  1054.             for (i = header.nx - 1; i >= 0; i--, n++) geoz_to_xy (header.x_min+i*header.x_inc+delx, header.y_min+dely, plane_level, &xx[n], &yy[n]);
  1055.             ps_polygon (xx, yy, n, red_facade, green_facade, blue_facade, TRUE);
  1056.         }
  1057.         if (!z_project.draw[1]) {    /*    Eastern side */
  1058.             for (j = n = 0, ij = ne; j < header.ny; j++, ij += mx) {
  1059.                 if (bad_float ((double)topo[ij])) continue;
  1060.                 geoz_to_xy (header.x_max-delx, header.y_max-j*header.y_inc-dely, (double)(topo[ij]), &xx[n], &yy[n]);
  1061.                 n++;
  1062.             }
  1063.             for (j = header.ny - 1; j >= 0; j--, n++) geoz_to_xy (header.x_max-delx, header.y_max-j*header.y_inc-dely, plane_level, &xx[n], &yy[n]);
  1064.             ps_polygon (xx, yy, n, red_facade, green_facade, blue_facade, TRUE);
  1065.         }
  1066.         if (!z_project.draw[2])    {    /* Northern side */
  1067.             for (i = n = 0, ij = nw; i < header.nx; i++, ij++) {
  1068.                 if (bad_float ((double)topo[ij])) continue;
  1069.                 geoz_to_xy (header.x_min+i*header.x_inc+delx, header.y_max-dely, (double)(topo[i]), &xx[n], &yy[n]);
  1070.                 n++;
  1071.             }
  1072.             for (i = header.nx - 1; i >= 0; i--, n++) geoz_to_xy (header.x_min+i*header.x_inc+delx, header.y_max-dely, plane_level, &xx[n], &yy[n]);
  1073.             ps_polygon (xx, yy, n, red_facade, green_facade, blue_facade, TRUE);
  1074.         }
  1075.         if (!z_project.draw[3]) {    /*    Western side */
  1076.             for (j = n = 0, ij = nw; j < header.ny; j++, ij += mx) {
  1077.                 if (bad_float ((double)topo[ij])) continue;
  1078.                 geoz_to_xy (header.x_min+delx, header.y_max-j*header.y_inc-dely, (double)(topo[ij]), &xx[n], &yy[n]);
  1079.                 n++;
  1080.             }
  1081.             for (j = header.ny - 1; j >= 0; j--, n++) geoz_to_xy (header.x_min+delx, header.y_max-j*header.y_inc-dely, plane_level, &xx[n], &yy[n]);
  1082.             ps_polygon (xx, yy, n, red_facade, green_facade, blue_facade, TRUE);
  1083.         }
  1084.     }
  1085.     
  1086.     if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);
  1087.     
  1088.     ps_plotend (gmtdefs.last_page);
  1089.     
  1090.     /* Free memory */
  1091.     
  1092.     if (get_contours) {
  1093.         for (ij = 0; ij < nm; ij++) {
  1094.             if (!binij[ij].first_cont) continue;
  1095.             last_cont = binij[ij].first_cont;
  1096.             for (this_cont = binij[ij].first_cont->next_cont; this_cont; this_cont = this_cont->next_cont) {
  1097.                 if (this_cont->first_point) {
  1098.                     last_point = this_cont->first_point;
  1099.                     for (this_point = this_cont->first_point->next_point; this_point; this_point = this_point->next_point) {
  1100.                         free ((char *)last_point);
  1101.                         last_point = this_point;
  1102.                     }
  1103.                     free ((char *)last_point);
  1104.                 }
  1105.                 free ((char *)last_cont);
  1106.                 last_cont = this_cont;
  1107.             }
  1108.             free ((char *)last_cont);
  1109.         }
  1110.         free ((char *)binij);
  1111.     }
  1112.     
  1113.     free ((char *)xx);
  1114.     free ((char *)yy);
  1115.     free ((char *)x);
  1116.     free ((char *)y);
  1117.     free ((char *)z);
  1118.     free ((char *)v);
  1119.     free ((char *)topo);
  1120.     if (intens) free ((char *)intensity);
  1121.     if (drape) free ((char *)grd);
  1122.     
  1123.     if (gmtdefs.verbose) fprintf (stderr, "grdview: Done!\n");
  1124.     
  1125.     gmt_end (argc, argv);
  1126. }
  1127.  
  1128. struct CONT *get_cont_struct (bin, value)
  1129. int bin;
  1130. double value; {
  1131.     struct CONT *cont, *new_cont;
  1132.     
  1133.     if (!binij[bin].first_cont) binij[bin].first_cont = (struct CONT *) memory (CNULL, 1, sizeof (struct CONT), "grdview");
  1134.  
  1135.     for (cont = binij[bin].first_cont; cont->next_cont && cont->next_cont->value <= value; cont = cont->next_cont);
  1136.     
  1137.     new_cont = (struct CONT *) memory (CNULL, 1, sizeof (struct CONT), "grdview");
  1138.     if (cont->next_cont) {    /* Put it in the link */
  1139.         new_cont->next_cont = cont->next_cont;
  1140.         cont->next_cont = new_cont;
  1141.     }
  1142.     else    /* End of list */
  1143.         cont->next_cont = new_cont;
  1144.     return (new_cont);
  1145. }
  1146.  
  1147. struct POINT *get_point (x, y)
  1148. double x, y; {
  1149.     struct POINT *point;
  1150.     
  1151.     point = (struct POINT *) memory (CNULL, 1, sizeof (struct POINT), "grdview");
  1152.     point->x = x;
  1153.     point->y = y;
  1154.     return (point);
  1155. }
  1156.  
  1157. int init_setup (header, topo, two, draw_plane, plane_level)
  1158. struct GRD_HEADER *header;
  1159. float topo[];
  1160. int two;
  1161. BOOLEAN draw_plane; 
  1162. double plane_level; {
  1163.     int i, j, ij;
  1164.     double xtmp, ytmp, tmp, delx, dely;
  1165.     
  1166.     delx = (header->node_offset) ? 0.5 * header->x_inc :0.0;
  1167.     dely = (header->node_offset) ? 0.5 * header->y_inc :0.0;
  1168.     /* Find projected min/max in y-direction */
  1169.     
  1170.     z_project.ymax = z_project.ymin;    /* Reset from whatever it was */
  1171.     
  1172.     for (j = 0; j < header->ny; j++) {
  1173.         ij = (j + two) * header->nx + two;
  1174.         for (i = 0; i < header->nx; i++, ij++) {
  1175.             if (bad_float ((double)topo[ij])) continue;
  1176.             geoz_to_xy (header->x_min + i * header->x_inc, header->y_max - j * header->y_inc, (double)topo[ij], &xtmp, &ytmp);
  1177.             z_project.ymin = MIN (z_project.ymin, ytmp);
  1178.             z_project.ymax = MAX (z_project.ymax, ytmp);
  1179.         }
  1180.     }
  1181.     if (draw_plane) {    /* plane or facade may exceed the found min/max */
  1182.         for (i = 0; i < header->nx; i++) {
  1183.             tmp = header->x_min + i * header->x_inc + delx;
  1184.             geoz_to_xy (tmp, header->y_min + dely, plane_level, &xtmp, &ytmp);
  1185.             z_project.ymin = MIN (z_project.ymin, ytmp);
  1186.             z_project.ymax = MAX (z_project.ymax, ytmp);
  1187.             geoz_to_xy (tmp, header->y_max-dely, plane_level, &xtmp, &ytmp);
  1188.             z_project.ymin = MIN (z_project.ymin, ytmp);
  1189.             z_project.ymax = MAX (z_project.ymax, ytmp);
  1190.         }
  1191.         for (j = 0; j < header->ny; j++) {
  1192.             tmp = header->y_max - j * header->y_inc - dely;
  1193.             geoz_to_xy (header->x_min+delx, tmp, plane_level, &xtmp, &ytmp);
  1194.             z_project.ymin = MIN (z_project.ymin, ytmp);
  1195.             z_project.ymax = MAX (z_project.ymax, ytmp);
  1196.             geoz_to_xy (header->x_max-delx, tmp, plane_level, &xtmp, &ytmp);
  1197.             z_project.ymin = MIN (z_project.ymin, ytmp);
  1198.             z_project.ymax = MAX (z_project.ymax, ytmp);
  1199.         }
  1200.     }
  1201. }
  1202.  
  1203. double get_intensity (intensity, k, nx)
  1204. float intensity[];
  1205. int k, nx; {
  1206.     /* Finds the agerage intensity for this polygon */
  1207.     return (0.25 * (intensity[k] + intensity[k+1] + intensity[k-nx] + intensity[k-nx+1]));
  1208. }
  1209.  
  1210. int pixel_inside (ip, jp, ix, iy, bin)
  1211. int ip, jp, ix[], iy[], bin; {
  1212.     int i, what;
  1213.     double x[6], y[6];
  1214.     
  1215.     for (i = 0; i < 4; i++) {
  1216.         x[i] = ix[bin+bin_inc[i]];
  1217.         y[i] = iy[bin+bin_inc[i]];
  1218.     }
  1219.     x[4] = x[0];    y[4] = y[0];
  1220.     what = non_zero_winding ((double)ip, (double)jp, x, y, 5);
  1221.     return (what);
  1222.     /* return (inside ((double)ip, (double)jp, x, y, 4)); */
  1223. }
  1224.  
  1225. int idist (x1, y1, x2, y2)
  1226. int x1, y1, x2, y2; {
  1227.     if ((x2 -= x1) < 0) x2 = -x2;
  1228.     if ((y2 -= y1) < 0) y2 = -y2;
  1229.     return (x2 + y2 - (((x2 > y2) ? y2 : x2) >> 1));
  1230. }
  1231.