home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / psmask.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-30  |  20.0 KB  |  624 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)psmask.c    2.28  30 Jun 1995
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  * psmask tries to achieve masking using one of two different approachs:
  9.  * The default way of operation is: Instead of painting tiles where there's no
  10.  * data [or where there is data] psmask uses contouring to find the polygons
  11.  * that contain the data [or contain the regions with no data].  For many types
  12.  * of data coverage, this results in a managable path instead of thousands of
  13.  * tiles.  So, instead of painting polygons, psmask sets up actual clip paths.
  14.  * As an option, the user may specify a rgb combination to fill the clipped areas.
  15.  * To avoid having to deal with the problems that arise if the data distribution
  16.  * is such that part of the map boundary should be part of the clip paths, we
  17.  * internally enlarge the grid by one gridsize unit so no nodes along the edges
  18.  * have data.  Before using the clippaths we move points outside the real region
  19.  * onto the boundary.
  20.  * Should the paths be too long for PostScript to handle, the user may override
  21.  * the default with the -T switch.  Then, masking is achieved using tiling
  22.  *
  23.  * Author:    Paul Wessel
  24.  * Date:    21-NOV-1990
  25.  * Version:    2.0
  26.  * Limitation:    Total clip path cannot exceed PostScript max path length (~1350).
  27.  *         If it does, the -T option be used instead.  Path length can be
  28.  *        reduced by choosing a coarser grid spacing.
  29.  *
  30.  */
  31.  
  32. #include "gmt.h"
  33.  
  34. #define IJ(i,j) ((i) + (j) * nx)
  35.  
  36. /* Global variables for contour_sub routines */
  37.  
  38. int p[5], i_off[5], j_off[5], k_off[5], offset;
  39. unsigned int bit[32];
  40.  
  41. char *grd;
  42. int *edge;
  43. double *x, *y;    /* Arrays holding the contour xy values */
  44. double *x_grid, *y_grid;
  45.  
  46. main (argc, argv)
  47. int argc;
  48. char **argv; {
  49.     int i, j, ij, n, nm, n_edges, di, dj, ii, jj;
  50.     int nx, ny, section, n_alloc, ix, iy, n_fields;
  51.     
  52.     BOOLEAN error = FALSE, first = TRUE, dump = FALSE, invert = FALSE, binary = FALSE, more;
  53.     BOOLEAN end_of_clip = FALSE, find_paths = TRUE, use_poly, single_precision = FALSE;
  54.     
  55.     char dfile[80], line[100];
  56.     
  57.     double west, east, north, south, dx, dy, radius = 0.0, *in;
  58.     double x0, y0, x1, y1, xx[4], yy[4], dummy;
  59.     
  60.     FILE *fp = NULL;
  61.     
  62.     struct FILL fill;
  63.     
  64.     west = east = north = south = 0.0;
  65.  
  66.     argc = gmt_begin (argc, argv);
  67.  
  68.     gmt_init_fill (&fill, -1, -1, -1);
  69.     dfile[0] = 0;
  70.     dx = dy = 0.0;
  71.     
  72.     for (i = 1; i < argc; i++) {
  73.         if (argv[i][0] == '-') {
  74.             switch (argv[i][1]) {
  75.         
  76.                 /* Common parameters */
  77.             
  78.                 case 'B':
  79.                 case 'J':
  80.                 case 'K':
  81.                 case 'O':
  82.                 case 'P':
  83.                 case 'R':
  84.                 case 'U':
  85.                 case 'V':
  86.                 case 'X':
  87.                 case 'x':
  88.                 case 'Y':
  89.                 case 'y':
  90.                 case 'c':
  91.                 case ':':
  92.                 case '\0':
  93.                     error += get_common_args (argv[i], &west, &east, &south, &north);
  94.                     break;
  95.                 
  96.                 /* Supplemental parameters */
  97.             
  98.                 case 'b':    /* Input pairs are binary, not ascii */
  99.                     single_precision = !(argv[i][strlen(argv[i])-1] == 'd');
  100.                     binary = TRUE;
  101.                     break;
  102.                 case 'C':    /* Radius of influence */
  103.                     gmt_getinc(&argv[i][2], &radius, &radius);
  104.                     break;
  105.                 case 'D':    /* Dump the polygons to files */
  106.                     strcpy (dfile, &argv[i][2]);
  107.                     dump = TRUE;
  108.                     break;
  109.                 case 'E':
  110.                     sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
  111.                     break;
  112.                 case 'F':
  113.                     if (gmt_getrgb (&argv[i][2], &gmtdefs.basemap_frame_rgb[0], &gmtdefs.basemap_frame_rgb[1], &gmtdefs.basemap_frame_rgb[2])) {
  114.                         gmt_pen_syntax ('F');
  115.                         error++;
  116.                     }
  117.                     break;
  118.                 case 'G':
  119.                     if (gmt_getfill (&argv[i][2], &fill)) {
  120.                         gmt_fill_syntax ('G');
  121.                         error++;
  122.                     }
  123.                     break;
  124.                 case 'I':
  125.                     gmt_getinc (&argv[i][2], &dx, &dy);
  126.                     break;
  127.                 case 'N':
  128.                     invert = TRUE;
  129.                     break;
  130.                 case 'S':
  131.                     end_of_clip = TRUE;
  132.                     break;
  133.                 case 'T':
  134.                     find_paths = FALSE;
  135.                     break;
  136.                 default:
  137.                     error = TRUE;
  138.                     gmt_default_error (argv[i][1]);
  139.                     break;
  140.             }
  141.         }
  142.         else {
  143.             if ((fp = fopen(argv[i], "r")) == NULL) {
  144.                 fprintf (stderr, "psmask: Could not open file %s\n", argv[i]);
  145.                 exit(-1);
  146.             }
  147.         }
  148.     }
  149.     
  150.     if (argc == 1 || gmt_quick) {
  151.         fprintf (stderr,"psmask %s - Clipping of 2-D data sets\n\n", GMT_VERSION);
  152.         fprintf (stderr, "usage: psmask <xyz-file> -I<xinc>[m|c][/<yinc>[m|c]] -J<params> -R<west/east/south/north>\n");
  153.         fprintf (stderr, "    [-B<tickinfo>] [-C<radius>[m|c]] -D<file> [-Eaz/el] [-F<r/g/b] [-G<fill>] [-K] [-N]\n");
  154.         fprintf (stderr, "    [-O] [-P] [-S] [-T] [-U[<label>]] [-V] [-X<x_shift>] [-Y<y_shift>] [-c<ncopies>][ -:] [-b[d]]\n\n");
  155.         
  156.         if (gmt_quick) exit(-1);
  157.         
  158.         fprintf (stderr, "    <xyz-file> is the datafile.  If not given, read standard input\n");
  159.         fprintf (stderr, "    -I sets grid increments; enter xinc, optionally xinc/yinc.\n");
  160.         fprintf (stderr, "       Default is yinc = xinc.  Append an m [or s] to xinc or yinc to indicate minutes [or seconds]\n");
  161.         explain_option ('j');
  162.         explain_option ('R');
  163.         fprintf (stderr, "\n\tOPTIONS:\n");
  164.         explain_option ('b');
  165.         fprintf (stderr, "    -C Circle of influence. Nodes inside circles of <radius> centered on\n");
  166.         fprintf (stderr, "       the input data poins are considered to be reliable estimates of the surface\n");
  167.         fprintf (stderr, "       Default is -C0, i.e., only the nearest node is considered reliable\n");
  168.         fprintf (stderr, "       <radius> is in the same units as x and y\n");
  169.         fprintf (stderr, "       Append an m [or c] to indicate minutes [or seconds]\n");
  170.         fprintf (stderr, "    -D dumps the clip-paths to files using the prefix <file>_\n");
  171.         fprintf (stderr, "       Ignored if -T is specified\n");
  172.         fprintf (stderr, "    -E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n");
  173.         fprintf (stderr, "    -S means stop existing clip-path (Only with -O) \n");
  174.         fprintf (stderr, "    -F Set color used for Frame and anotation [%d/%d/%d]\n",
  175.             gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  176.         fprintf (stderr, "    -G sets fill color r/g/b [Default is no fill]\n");
  177.         explain_option ('K');
  178.         fprintf (stderr, "    -N will invert the sense of the clipping [or tiling]\n");
  179.         explain_option ('O');
  180.         explain_option ('P');
  181.         fprintf (stderr, "    -T will paint tiles.  [Default will trace data outline]\n");
  182.         fprintf (stderr, "       If set you must also specify a color/fill with -G\n");
  183.         fprintf (stderr, "       !May only be used with linear, Mercator, and basic cylindrical projections!\n");
  184.         explain_option ('U');
  185.         explain_option ('V');
  186.         explain_option ('X');
  187.         explain_option ('c');
  188.         explain_option (':');
  189.         fprintf (stderr, "\t-b means input data are binary.  [Default is ascii]\n");
  190.         fprintf (stderr, "\t   Append d for double precision [Default is single].\n");
  191.         explain_option ('.');
  192.         exit(-1);
  193.     }
  194.     
  195.     if (!end_of_clip) {
  196.         if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) {
  197.             fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  Enter azimuth in 0-360 range, elevation in 0-90 range\n", gmt_program);
  198.             error++;
  199.         }
  200.         if (!project_info.region_supplied) {
  201.             fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", gmt_program);
  202.             error++;
  203.         }
  204.         if (!find_paths && !RECT_GRATICULE) {
  205.             fprintf (stderr, "%s: GMT SYNTAX ERROR -T option:  Only available with Linear, Mercator, or basic cylindrical projections\n", gmt_program);
  206.             error++;
  207.         }
  208.         if (!find_paths && fill.r == -1) {
  209.             fprintf (stderr, "%s: GMT SYNTAX ERROR -T option:  Must also specify a tile color with -G\n", gmt_program);
  210.             error++;
  211.         }
  212.         if (find_paths && fill.use_pattern) {
  213.             fprintf (stderr, "%s: GMT SYNTAX ERROR -G option:  Pattern fill not supported\n", gmt_program);
  214.             error++;
  215.         }
  216.         if (dx <= 0.0 || dy <= 0.0) {
  217.             fprintf (stderr, "%s: GMT SYNTAX ERROR -I option:  Must specify positive increments\n", gmt_program);
  218.             error++;
  219.         }
  220.         if (radius < 0.0) {
  221.             fprintf (stderr, "%s: GMT SYNTAX ERROR -C option:  Radius must be positive\n", gmt_program);
  222.             error++;
  223.         }
  224.     }
  225.     if (binary && gmtdefs.io_header) {
  226.         fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", gmt_program);
  227.         error++;
  228.     }
  229.  
  230.     if (error) exit (-1);
  231.  
  232.     if (binary) gmt_input = (single_precision) ? bin_float_input : bin_double_input;
  233.  
  234.     if (end_of_clip) gmtdefs.overlay = TRUE;    /* Force overlay mode */
  235.     if (!project_info.x_off_supplied && gmtdefs.overlay) gmtdefs.x_origin = 0.0;
  236.     if (!project_info.y_off_supplied && gmtdefs.overlay) gmtdefs.y_origin = 0.0;
  237.  
  238.     ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
  239.         gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
  240.         gmtdefs.dpi, gmtdefs.measure_unit, gmtdefs.paper_width, gmtdefs.page_rgb, gmt_epsinfo (argv[0]));
  241.     echo_command (argc, argv);
  242.     if (gmtdefs.unix_time) timestamp (argc, argv);
  243.     
  244.     if (!end_of_clip) {    /* Start new clip_path */
  245.     
  246.         map_setup (west, east, south, north);
  247.         
  248.         if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0);
  249.  
  250.         if (gmtdefs.verbose) fprintf (stderr, "psmask: Allocate memory, read and process data file\n");
  251.     
  252.         ix = (gmtdefs.xy_toggle) ? 1 : 0;        iy = 1 - ix;              /* Set up which columns have x and y */
  253.         
  254.         /* Enlarge region by 1 row/column */
  255.  
  256.         west -= dx;    east += dx;    south -= dy;    north += dy;
  257.  
  258.         nx = rint ((east - west)/dx) + 1;
  259.         ny = rint ((north - south)/dy) + 1;
  260.     
  261.         nm = nx * ny;
  262.         grd = (char *) memory (CNULL, nm, sizeof (char), "psmask");
  263.     
  264.         di = rint (0.5 * radius / dx);
  265.         dj = rint (0.5 * radius / dy);
  266.         if (fp == NULL) fp = stdin;
  267.         if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (line, 512, fp);
  268.  
  269.         while ((n_fields = gmt_input (fp, 3, &in)) == 3) {
  270.             i = rint ((in[ix] - west)/dx);
  271.             if (i < 0 || i >= nx) continue;
  272.             j = rint ((north - in[iy])/dy);
  273.             if (j < 0 || j >= ny) continue;
  274.             x0 = west + i * dx;
  275.             y0 = north - j * dy;
  276.         
  277.             for (ii = i - di; ii <= i + di; ii++) {
  278.                 if (ii < 0 || ii >= nx) continue;
  279.                 x1 = west + ii * dx;
  280.                 for (jj = j - dj; jj <= j + dj; jj++) {
  281.                     if (jj < 0 || jj >= ny) continue;
  282.                     y1 = north - jj * dy;
  283.                     if (hypot (x1 - x0, y1 - y0) > radius) continue;
  284.                     grd[jj*nx+ii] = 1;
  285.                 }
  286.             }
  287.         }
  288.         if (fp != stdin) fclose (fp);
  289.         
  290.         if (invert) for (i = 0; i < nm; i++) grd[i] = 1 - grd[i];    /* Reverse sense of test */
  291.  
  292.         /* Force perimeter nodes to be FALSE */
  293.         
  294.         for (i = 0, ij = (ny-1) * nx; i < nx; i++) grd[i] = grd[i+ij] = FALSE;
  295.         for (j = 0; j < ny; j++) grd[j*nx] = grd[(j+1)*nx-1] = FALSE;
  296.         
  297.         if (find_paths) {    /* Must trace the outline of ON/OFF values in grd */
  298.             x = (double *) memory (CNULL, GMT_CHUNK, sizeof (double), "psmask");
  299.             y = (double *) memory (CNULL, GMT_CHUNK, sizeof (double), "psmask");
  300.             n_alloc = GMT_CHUNK;
  301.     
  302.             n_edges = ny * (int )ceil (nx / 16.0);
  303.             edge = (int *) memory (CNULL, n_edges, sizeof (int), "psmask");
  304.     
  305.             if (frame_info.plot) {
  306.                 ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  307.                 map_basemap ();
  308.                 ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  309.             }
  310.  
  311.             if (gmtdefs.verbose) fprintf (stderr, "psmask: Tracing the clip path\n");
  312.     
  313.             section = 0;
  314.             first = TRUE;
  315.             while ((n = clip_contours (grd, nx, ny, west, north, dx, dy, edge, first, &x, &y, &n_alloc)) > 0) {
  316.                 shrink_clip_contours (x, y, n);
  317.                 if (dump) dump_clip_countours (x, y, n, section, dfile);
  318.                 draw_clip_contours (x, y, n, fill.r, fill.g, fill.b, section, first);
  319.                 first = FALSE;
  320.                 section++;
  321.             }
  322.     
  323.             draw_clip_contours (x, y, 0, fill.r, fill.g, fill.b, section, 2);    /* Activate clip-path */
  324.         
  325.             free ((char *)edge);
  326.             free ((char *)x);
  327.             free ((char *)y);
  328.         }
  329.         else {    /* Just paint tiles */
  330.             if (gmtdefs.verbose) fprintf (stderr, "psmask: Tiling...\n");
  331.             x_grid = (double *) memory (CNULL, nx, sizeof (double), "psmask");
  332.             y_grid = (double *) memory (CNULL, ny, sizeof (double), "psmask");
  333.             for (i = 2; i < nx - 1; i++) geo_to_xy (west + (i - 0.5) * dx, south + dy, &x_grid[i], &dummy);
  334.             for (j = 2; j < ny - 1; j++) geo_to_xy (west + dx, north - (j - 0.5) * dy, &dummy, &y_grid[j]);
  335.             geo_to_xy (west + dx, north - dy, &x_grid[1], &y_grid[1]);
  336.             geo_to_xy (east - dx, south + dy, &x_grid[nx-1], &y_grid[ny-1]);
  337.             
  338.             use_poly = (fill.use_pattern || project_info.three_D);
  339.             
  340.             for (j = 1; j < ny; j++) {
  341.                 ij = j * nx + 1;
  342.                 for (i = 1; i < nx; i++, ij++) {
  343.                     if (!grd[ij]) continue;
  344.                     
  345.                     if (use_poly) {
  346.                         xx[0] = xx[3] = x_grid[i];    xx[1] = xx[2] = x_grid[i+1];
  347.                         yy[0] = yy[1] = y_grid[j+1];    yy[2] = yy[3] = y_grid[j];
  348.                         if (project_info.three_D) two_d_to_three_d (xx, yy, 4);
  349.                     }
  350.                     
  351.                     if (fill.use_pattern)
  352.                         ps_imagefill (xx, yy, 4, fill.pattern_no, fill.pattern, fill.inverse, fill.icon_size, FALSE);
  353.                     else if (project_info.three_D)
  354.                         ps_polygon (xx, yy, 4, fill.r, fill.g, fill.b, FALSE);
  355.                     else
  356.                         ps_rect (x_grid[i], y_grid[j+1], x_grid[i+1], y_grid[j], fill.r, fill.g, fill.b, FALSE);
  357.                 }
  358.             }
  359.             
  360.             free ((char *)x_grid);
  361.             free ((char *)y_grid);
  362.             
  363.             if (frame_info.plot) {
  364.                 ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  365.                 map_basemap ();
  366.             }
  367.         }
  368.         
  369.         if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);
  370.  
  371.         free ((char *)grd);
  372.         if (gmtdefs.verbose) fprintf (stderr, "psmask: clipping on!\n");
  373.     }
  374.     else {    /* Just undo previous clip-path */
  375.         ps_clipoff ();
  376.         if (frame_info.plot) {
  377.             ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  378.             map_basemap ();
  379.         }
  380.         if (gmtdefs.verbose) fprintf (stderr, "psmask: clipping off!\n");
  381.     }
  382.     
  383.     ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  384.     
  385.     ps_plotend (gmtdefs.last_page);
  386.     
  387.     
  388.     gmt_end (argc, argv);
  389. }
  390.  
  391. int draw_clip_contours (xx, yy, nn, r, g, b, id, flag)
  392. double *xx, *yy;
  393. int nn, r, g, b, id, flag; {
  394.     int i;
  395.     double x, y;
  396.     
  397.     if (nn < 2 && flag < 2) return;
  398.     
  399.     for (i = 0; i < nn; i++) {
  400.         x = xx[i];    y = yy[i];
  401.         geo_to_xy (x, y, &xx[i], &yy[i]);
  402.     }
  403.     nn = compact_line (xx, yy, nn, FALSE, 0);
  404.     
  405.     if (project_info.three_D) two_d_to_three_d (xx, yy, nn);
  406.     
  407.     if (nn > 0) printf ("%% Start of clip path sub-segment %d\n", id);
  408.     ps_clipon (xx, yy, nn, r, g, b, flag);
  409.     if (nn > 0) printf ("%% End of clip path sub-segment %d\n", id);
  410.         
  411. }
  412.  
  413. int dump_clip_countours (xx, yy, nn, id, file)
  414. double *xx, *yy;
  415. char *file;
  416. int nn, id; {
  417.     int i;
  418.     char fname[512], format[80];
  419.     FILE *fp;
  420.     
  421.     if (nn < 2) return;
  422.     
  423.     sprintf (format, "%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format);
  424.     sprintf (fname, "%s_%d.xy\0", file, id);
  425.     fp = fopen (fname, "w");
  426.     for (i = 0; i < nn; i++) fprintf (fp, format, xx[i], yy[i]);
  427.     fclose (fp);
  428. }
  429.  
  430. int clip_contours (grd, nx, ny, west, north, dx, dy, edge, first, x, y, max)
  431. char grd[];
  432. int nx, ny, first, *max;
  433. double west, north, dx, dy;
  434. int edge[];
  435. double *x[], *y[]; {
  436.     /* The routine finds the zero-contour in the grd dataset.  it assumes that
  437.      * no node has a value exactly == 0.0.  If more than max points are found
  438.      * trace_clip_contours will try to allocate more memory in blocks of 1000 points
  439.      */
  440.      
  441.      static int i0, j0;
  442.      int i, j, ij, n = 0, n_edges, edge_word, edge_bit;
  443.      BOOLEAN go_on = TRUE;
  444.      
  445.      
  446.      n_edges = ny * (int) ceil (nx / 16.0);
  447.      offset = n_edges / 2;
  448.      
  449.      /* Reset edge-flags to zero, if necessary */
  450.      if (first) {
  451.          i0 = j0 = 0;
  452.         p[0] = p[4] = 0;    p[1] = 1;    p[2] = 1 - nx;    p[3] = -nx;
  453.         i_off[0] = i_off[2] = i_off[3] = i_off[4] = 0;    i_off[1] =  1;
  454.         j_off[0] = j_off[1] = j_off[3] = j_off[4] = 0;    j_off[2] = -1;
  455.         k_off[0] = k_off[2] = k_off[4] = 0;    k_off[1] = k_off[3] = 1;
  456.         for (i = 1, bit[0] = 1; i < 32; i++) bit[i] = bit[i-1] << 1;
  457.      }
  458.  
  459.     /* Loop over interior boxes */
  460.         
  461.     for (j = j0; go_on && j < ny; j++) {
  462.         ij = i0 + j * nx;
  463.         for (i = i0; go_on && i < nx-1; i++, ij++) {
  464.             edge_word = ij / 32 + offset;
  465.             edge_bit = ij % 32;
  466.             if (!(edge[edge_word] & bit[edge_bit]) && ((grd[ij]+grd[ij-nx]) == 1)) { /* Start tracing countour */
  467.                 *x[0] = west + i*dx;
  468.                 *y[0] = north - (j-0.5) * dy;
  469.                 edge[edge_word] |= bit[edge_bit];
  470.                 n = trace_clip_contours (grd, nx, edge, west, north, dx, dy, x, y, i, j, 3, max);
  471.                 go_on = FALSE;
  472.                 i0 = i + 1;
  473.                 j0 = j;
  474.             }
  475.         }
  476.         if (!go_on) i0 = 1;
  477.     }
  478.     
  479.     return (n);
  480. }
  481.  
  482. int trace_clip_contours (grd, nx, edge, west, north, dx, dy, xx, yy, i, j, kk, max)
  483. char grd[];
  484. int nx, i, j, kk, *max;
  485. int edge[];
  486. double west, north, dx, dy;
  487. double **xx, **yy; {
  488.     int n = 1, k, k0, ij, ij0, n_cuts, kk_opposite, first_k, more;
  489.     int edge_word, edge_bit, m;
  490.     double xk[4], yk[4], x0, y0;
  491.     
  492.     m = *max - 2;
  493.     
  494.     more = TRUE;
  495.     do {
  496.         ij = i + j * nx;
  497.         x0 = west + i * dx;
  498.         y0 = north - j * dy;
  499.         n_cuts = 0;
  500.         k0 = kk;
  501.  
  502.         for (k = 0; k < 4; k++) {    /* Loop over box sides */
  503.  
  504.             /* Skip where we already have a cut (k == k0) */
  505.             
  506.             if (k == k0) continue;
  507.             
  508.             /* Skip edge already has been used */
  509.             
  510.             ij0 = IJ (i + i_off[k], j + j_off[k]);
  511.             edge_word = ij0 / 32 + k_off[k] * offset;
  512.             edge_bit = ij0 % 32;
  513.             if (edge[edge_word] & bit[edge_bit]) continue;
  514.             
  515.             /* Skip if no zero-crossing on this edge */
  516.             
  517.             if ((grd[ij+p[k+1]] + grd[ij+p[k]]) != 1) continue;
  518.             
  519.             /* Here we have a crossing */
  520.             
  521.             if (k%2) {    /* Cutting a S-N line */
  522.                 if (k == 1) {
  523.                     xk[1] = x0 + dx;
  524.                     yk[1] = y0 + 0.5*dy;
  525.                 }
  526.                 else {
  527.                     xk[3] = x0;
  528.                     yk[3] = y0 + 0.5*dy;
  529.                 }
  530.             }
  531.             else {    /* Cutting a E-W line */
  532.                 if (k == 0) {
  533.                     xk[0] = x0 + 0.5*dx;
  534.                     yk[0] = y0;
  535.                 }
  536.                 else {
  537.                     xk[2] = x0 + 0.5*dx;
  538.                     yk[2] = y0 + dy;
  539.                 }
  540.             }
  541.             kk = k;
  542.             n_cuts++;
  543.         }
  544.         
  545.         if (n > m) {    /* Must try to allocate more memory */
  546.             *max += 1000;
  547.             m += 1000;
  548.             *xx = (double *) memory ((char *)*xx, (*max), sizeof (double), "trace_clip_contours");
  549.             *yy = (double *) memory ((char *)*yy, (*max), sizeof (double), "trace_clip_contours");
  550.         }
  551.         if (n_cuts == 0) {    /* Close interior contour and return */
  552.             if (fmod (*xx[0], dx) == 0.0)    /* On side 1 or 3 */
  553.                 first_k = ((*xx[0] - x0) == 0.0) ? 3 : 1;
  554.             else     /* On side 0 or 2 */
  555.                 first_k = ((*yy[0] - y0) == 0.0) ? 0 : 2;
  556.             kk_opposite = (first_k + 2) % 4;
  557.             if (k0 != kk_opposite) {
  558.                 (*xx)[n] = x0 + 0.5*dx;
  559.                 (*yy)[n] = y0 + 0.5*dy;
  560.                 n++;
  561.             }
  562.             (*xx)[n] = (*xx)[0];
  563.             (*yy)[n] = (*yy)[0];
  564.             n++;
  565.             more = FALSE;
  566.         }
  567.         else if (n_cuts == 1) {    /* Draw a line to this point and keep tracing */
  568.             /* Add center of box if this and previous cut NOT on opposite edges */
  569.             kk_opposite = (k0 + 2) % 4;
  570.             if (kk != kk_opposite) {
  571.                 (*xx)[n] = x0 + 0.5*dx;
  572.                 (*yy)[n] = y0 + 0.5*dy;
  573.                 n++;
  574.             }
  575.             (*xx)[n] = xk[kk];
  576.             (*yy)[n] = yk[kk];
  577.             n++;
  578.         }
  579.         else {    /* Saddle point, we decide to connect to the point nearest previous point */
  580.             kk = (k0 + 1)%4;    /* Pick next edge since it is arbitrarely where we go */
  581.             /* First add center of box */
  582.             (*xx)[n] = x0 + 0.5*dx;
  583.             (*yy)[n] = y0 + 0.5*dy;
  584.             n++;
  585.             (*xx)[n] = xk[kk];
  586.             (*yy)[n] = yk[kk];
  587.             n++;
  588.         }
  589.         if (more) {    /* Mark this edge as used */
  590.             ij0 = IJ (i + i_off[kk], j + j_off[kk]);
  591.             edge_word = ij0 / 32 + k_off[kk] * offset;
  592.             edge_bit = ij0 % 32;
  593.             edge[edge_word] |= bit[edge_bit];    
  594.         }
  595.         
  596.         /* Get next box (i,j,kk) */
  597.         
  598.         i -= (kk-2)%2;
  599.         j -= (kk-1)%2;
  600.         kk = (kk+2)%4;
  601.         
  602.     } while (more);
  603.     return (n);
  604. }
  605.  
  606. int shrink_clip_contours (x, y, n)
  607. double x[], y[];
  608. int n; {
  609.     /* Moves outside points to boundary */
  610.     int i;
  611.  
  612.     for (i = 0; i < n; i++) {
  613.         if (x[i] < project_info.w)
  614.             x[i] = project_info.w;    
  615.         if (x[i] > project_info.e)
  616.             x[i] = project_info.e;
  617.         if (y[i] < project_info.s)
  618.             y[i] = project_info.s;
  619.         if (y[i] > project_info.n)
  620.             y[i] = project_info.n;
  621.     }
  622. }
  623.  
  624.