home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdimage.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-22  |  12.8 KB  |  414 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:   @(#)grdimage.c    2.42  22 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.  * grdimage will read a grdfile and image the area using the PostScript
  9.  * image command. If non-linear scaling is chosen, we first resample the data
  10.  * onto the new grid before calling the image command.  THe output image
  11.  * will be 1-, 8-, or 24-bit depending on colors used.
  12.  *
  13.  * Author:    Paul Wessel
  14.  * Date:    21-MAY-1991-1995
  15.  * Ver:        3.0    based on old 2.x
  16.  */
  17.  
  18. #include "gmt.h"
  19.  
  20. #define YIQ(r, g, b) rint (0.299 * (r) + 0.587 * (g) + 0.114 * (b))    /* How B/W TV's convert RGB to Gray */
  21.  
  22. float *tmp1, *tmp2, *map, *intensity;
  23. unsigned char *bitimage_8, *bitimage_24;
  24.  
  25. char *c_method[3] = {
  26.     "colorimage",
  27.     "colortiles",
  28.     "RGB-separation"
  29. };
  30.  
  31. main (argc, argv)
  32. int argc;
  33. char **argv; {
  34.     int i, j, k, kk, r, g, b, nm, nm2, byte, off, nx_f, ny_f, grid_type;
  35.     int nx, ny, pad[4], dpi = 0, nx_proj = 0, ny_proj = 0, node;
  36.     
  37.     BOOLEAN error = FALSE, intens = FALSE, monochrome = FALSE;
  38.     BOOLEAN mapped = FALSE, set_dpi = FALSE;
  39.     
  40.     double  dx, dy, x_side, y_side, x0 = 0.0, y0 = 0.0;
  41.     double west, east, south, north, data_west, data_east, data_south, data_north;
  42.     
  43.     char *grdfile, *intensfile, *cpt_file;
  44.     
  45.     struct GRD_HEADER g_head, r_head, i_head, j_head;
  46.     
  47.     argc = gmt_begin (argc, argv);
  48.  
  49.     grd_init (&g_head, argc, argv, FALSE);
  50.     grd_init (&r_head, argc, argv, FALSE);
  51.     
  52.     grdfile = intensfile = cpt_file = CNULL;
  53.     west = east = south = north = 0.0;
  54.     
  55.     for (i = 1; i < argc; i++) {
  56.         if (argv[i][0] == '-') {
  57.             switch (argv[i][1]) {
  58.                 /* Common parameters */
  59.             
  60.                 case 'B':
  61.                 case 'J':
  62.                 case 'K':
  63.                 case 'O':
  64.                 case 'P':
  65.                 case 'R':
  66.                 case 'U':
  67.                 case 'V':
  68.                 case 'X':
  69.                 case 'x':
  70.                 case 'Y':
  71.                 case 'y':
  72.                 case 'c':
  73.                 case '\0':
  74.                     error += get_common_args (argv[i], &west, &east, &south, &north);
  75.                     break;
  76.                 
  77.                 /* Supplemental parameters */
  78.             
  79.                 case 'C':
  80.                     cpt_file = &argv[i][2];
  81.                     break;
  82.                 case 'E':
  83.                     dpi = atoi (&argv[i][2]);
  84.                     set_dpi = TRUE;
  85.                     break;
  86.                 case 'F':
  87.                     if (gmt_getrgb (&argv[i][2], &gmtdefs.basemap_frame_rgb[0], &gmtdefs.basemap_frame_rgb[1], &gmtdefs.basemap_frame_rgb[2])) {
  88.                         gmt_pen_syntax ('F');
  89.                         error++;
  90.                     }
  91.                     break;
  92.                 case 'I':
  93.                     intens = TRUE;
  94.                     intensfile = &argv[i][2];
  95.                     break;
  96.                 case 'M':
  97.                     monochrome = TRUE;
  98.                     break;
  99.                 case '0':
  100.                     gmtdefs.color_image = 0;
  101.                     break;
  102.                 case '1':
  103.                     gmtdefs.color_image = 1;
  104.                     break;
  105.                 case '2':
  106.                     gmtdefs.color_image = 2;
  107.                     break;
  108.                 default:
  109.                     error = TRUE;
  110.                     gmt_default_error (argv[i][1]);
  111.                     break;
  112.             }
  113.         }
  114.         else
  115.             grdfile = argv[i];
  116.     }
  117.     
  118.     if (argc == 1 || gmt_quick) {
  119.         fprintf (stderr,"grdimage %s - Plot grdfiles in 2-D\n\n", GMT_VERSION);
  120.         fprintf (stderr, "usage: grdimage <grdfile> -C<cpt_file> -J<params> [-B<tickinfo>] [-E<dpi>] [-F<r/g/b>]\n");
  121.         fprintf (stderr, "   [-I<intensity_file>] [-K] [-M] [-O] [-P] [-R<w/e/s/n>] [-U] [-V] [-X<x_shift>]\n");
  122.         fprintf (stderr, "   [-Y<y_shift>] [-0 -1 -2] [-c<ncopies>]\n\n");
  123.         
  124.         if (gmt_quick) exit (-1);
  125.         
  126.         fprintf (stderr,"    <grdfile> is data set to be plotted\n");
  127.         fprintf (stderr,"    -C color palette file\n");
  128.         explain_option ('j');
  129.         fprintf (stderr,"\n\tOPTIONS:\n");
  130.         explain_option ('b');
  131.         fprintf(stderr, "    -E sets dpi for the projected grid which must be constructed\n");
  132.         fprintf(stderr, "       if -Jx or -Jm is not selected [Default gives same size as input grid]\n");
  133.         fprintf (stderr, "    -F Set color used for Frame and anotation [%d/%d/%d]\n",
  134.             gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  135.         fprintf (stderr,"    -I use illumination.  Append name of intensity grd file\n");
  136.         explain_option ('K');
  137.         fprintf (stderr, "    -M force monochrome image\n");
  138.         explain_option ('O');
  139.         explain_option ('P');
  140.         explain_option ('R');
  141.         explain_option ('U');
  142.         explain_option ('V');
  143.         explain_option ('X');
  144.         fprintf (stderr,"    -0 plot color image using colorimage.\n");
  145.         fprintf (stderr,"    -1 plot color image using tiles.\n");
  146.         fprintf (stderr,"    -2 plot color image using psto24.\n");
  147.         fprintf (stderr,"       [Default = %d]\n", gmtdefs.color_image);
  148.         explain_option ('c');
  149.         explain_option ('.');
  150.         exit(-1);
  151.     }
  152.     
  153.     if (!grdfile) {
  154.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", gmt_program);
  155.         error++;
  156.     }
  157.     if (!cpt_file) {
  158.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify color palette table\n", gmt_program);
  159.         error++;
  160.     }
  161.     if (intens && !intensfile) {
  162.         fprintf (stderr, "%s: GMT SYNTAX ERROR -I option:  Must specify intensity file\n", gmt_program);
  163.         error++;
  164.     }
  165.     if (set_dpi && dpi <= 0) {
  166.         fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  dpi must be positive\n", gmt_program);
  167.         error++;
  168.     }
  169.     if (error) exit (-1);
  170.     
  171.     /* Get color palette file */
  172.     
  173.     read_cpt (cpt_file);
  174.  
  175.     /* Check limits and get data file */
  176.     
  177.     if (gmtdefs.verbose) fprintf (stderr, "grdimage: Allocates memory and read data file\n");
  178.     
  179.     if (read_grd_info (grdfile, &g_head)) {
  180.         fprintf (stderr, "grdimage: Error opening file %s\n", grdfile);
  181.         exit (-1);
  182.     }
  183.     off = (g_head.node_offset) ? 0 : 1;
  184.  
  185.     /* Determine what wesn to pass to map_setup */
  186.  
  187.     if (west == east && south == north) {
  188.         west = g_head.x_min;
  189.         east = g_head.x_max;
  190.         south = g_head.y_min;
  191.         north = g_head.y_max;
  192.     }
  193.  
  194.     map_setup (west, east, south, north);
  195.  
  196.     /* Determine the wesn to be used to read the grdfile */
  197.  
  198.     grd_setregion (&g_head, &data_west, &data_east, &data_south, &data_north);
  199.  
  200.     nx_f = g_head.nx;
  201.     ny_f = g_head.ny;
  202.     pad[0] = pad[1] = pad[2] = pad[3] = 0;
  203.     
  204.     /* Read data */
  205.  
  206.     nx = rint ( (data_east - data_west) / g_head.x_inc) + off;
  207.     ny = rint ( (data_north - data_south) / g_head.y_inc) + off;
  208.     nm = nx * ny;
  209.     tmp1 = (float *) memory (CNULL, nm, sizeof (float), "grdimage");
  210.     if (read_grd (grdfile, &g_head, tmp1, data_west, data_east, data_south, data_north, pad, FALSE)) {
  211.         fprintf (stderr, "grdcontour: Error reading file %s\n", grdfile);
  212.         exit (-1);
  213.     }
  214.  
  215.     /* If given, get intensity file or compute intensities */
  216.     
  217.     if (intens) {    /* Illumination wanted */
  218.     
  219.         if (gmtdefs.verbose) fprintf (stderr, "grdimage: Allocates memory and read intensity file\n");
  220.         
  221.         grd_init (&i_head, argc, argv, FALSE);
  222.         grd_init (&j_head, argc, argv, FALSE);
  223.         
  224.         if (read_grd_info (intensfile, &i_head)) {
  225.             fprintf (stderr, "grdimage: Error opening file %s\n", intensfile);
  226.             exit (-1);
  227.         }
  228.  
  229.         if (i_head.nx != nx_f || i_head.ny != ny_f) {
  230.             fprintf (stderr, "grdimage: Intensity file has improper dimensions!\n");
  231.             exit (-1);
  232.         }
  233.         tmp2 = (float *) memory (CNULL, nm, sizeof (float), "grdimage");
  234.         
  235.         if (read_grd (intensfile, &i_head, tmp2, data_west, data_east, data_south, data_north, pad, FALSE)) {
  236.             fprintf (stderr, "grdimage: Error reading file %s\n", intensfile);
  237.             exit (-1);
  238.         }
  239.     }
  240.     
  241.     r_head.x_min = project_info.xmin;    r_head.x_max = project_info.xmax;
  242.     r_head.y_min = project_info.ymin;    r_head.y_max = project_info.ymax;
  243.         
  244.     if (dpi > 0 || (project_info.projection > 5)) {    /* Need to resample the grd file */
  245.         
  246.         mapped = TRUE;
  247.         
  248.         if (gmtdefs.verbose) fprintf (stderr, "grdimage: project grdfiles\n");
  249.     
  250.         if (dpi == 0) {    /* Use input # of nodes as # of projected nodes */
  251.             nx_proj = g_head.nx;
  252.             ny_proj = g_head.ny;
  253.         }
  254.         grid_type = (dpi > 0) ? 1 : g_head.node_offset;    /* Force pixel if dpi is set */
  255.         grdproject_init (&r_head, 0.0, 0.0, nx_proj, ny_proj, dpi, grid_type);
  256.         nm2 = r_head.nx * r_head.ny;
  257.         map = (float *) memory (CNULL, nm2, sizeof (float), "grdproject");
  258.  
  259.         grd_forward (tmp1, &g_head, map, &r_head, 0.0, FALSE);
  260.         
  261.         free ((char *)tmp1);
  262.         if (intens) {
  263.             j_head.x_min = r_head.x_min;    j_head.x_max = r_head.x_max;
  264.             j_head.y_min = r_head.y_min;    j_head.y_max = r_head.y_max;
  265.         
  266.             if (dpi == 0) {    /* Use input # of nodes as # of projected nodes */
  267.                 nx_proj = i_head.nx;
  268.                 ny_proj = i_head.ny;
  269.             }
  270.             grdproject_init (&j_head, 0.0, 0.0, nx_proj, ny_proj, dpi, grid_type);
  271.             intensity = (float *) memory (CNULL, nm2, sizeof (float), "grdproject");
  272.             grd_forward (tmp2, &i_head, intensity, &j_head, 0.0, FALSE);
  273.             free ((char *)tmp2);
  274.         }
  275.         nm = nm2;
  276.     }
  277.     else {    /* Simply copy g_head info to r_head */
  278.         map = tmp1;
  279.         r_head.nx = g_head.nx;        r_head.ny = g_head.ny;
  280.         r_head.x_inc = g_head.x_inc;    r_head.y_inc = g_head.y_inc;
  281.         if (intens) {
  282.             j_head.nx = i_head.nx;        j_head.ny = i_head.ny;
  283.             j_head.x_inc = i_head.x_inc;    j_head.y_inc = i_head.y_inc;
  284.             intensity = tmp2;
  285.         }
  286.         grid_type = g_head.node_offset;
  287.     }
  288.     
  289.     ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
  290.         gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
  291.         gmtdefs.dpi, gmtdefs.measure_unit , gmtdefs.paper_width, gmtdefs.page_rgb, gmt_epsinfo (argv[0]));
  292.     echo_command (argc, argv);
  293.     if (gmtdefs.unix_time) timestamp (argc, argv);
  294.     
  295.     if (gmtdefs.verbose) fprintf (stderr, "grdimage: Evaluate pixel colors\n");
  296.     
  297.     if (monochrome || gmt_gray)
  298.         bitimage_8 = (unsigned char *) memory (CNULL, nm, sizeof (char), "grdimage");
  299.     else
  300.         bitimage_24 = (unsigned char *) memory (CNULL, 3 * nm, sizeof (char), "grdimage");
  301.     
  302.     nm2 = 2 * nm;
  303.     
  304.     for (j = byte = 0; j < r_head.ny; j++) {
  305.         kk = r_head.nx * (project_info.xyz_pos[1] ? j : r_head.ny - j - 1);
  306.         for (i = 0; i < r_head.nx; i++) {    /* Compute rgb for each pixel */
  307.             node = kk + (project_info.xyz_pos[0] ? i : r_head.nx - i - 1);
  308.             get_rgb24 (map[node], &r, &g, &b);
  309.             if (intens) illuminate (intensity[node], &r, &g, &b);
  310.  
  311.             if (gmt_gray)    /* Color table only has grays, pick r */
  312.                 bitimage_8[byte++] = (unsigned char) r;
  313.             else if (monochrome)    /* Convert rgb to gray using the YIQ transformation */
  314.                 bitimage_8[byte++] = (unsigned char) YIQ (r, g, b);
  315.             else if (gmtdefs.color_image == 2) {    /* RGB separation */
  316.                 bitimage_24[byte] = r;
  317.                 bitimage_24[byte+nm] = g;
  318.                 bitimage_24[byte+nm2] = b;
  319.                 byte++;
  320.             }
  321.             else {
  322.                 bitimage_24[byte++] = r;
  323.                 bitimage_24[byte++] = g;
  324.                 bitimage_24[byte++] = b;
  325.             }
  326.         }
  327.     }
  328.     
  329.     free ((char *)map);
  330.     if (intens) free ((char *)intensity);
  331.     
  332.     dx = (r_head.x_max - r_head.x_min) / (r_head.nx - off);
  333.     dy = (r_head.y_max - r_head.y_min) / (r_head.ny - off);
  334.  
  335.     x0 = r_head.x_min;    y0 = r_head.y_min;
  336.     if (mapped || grid_type == 0) {    /* Grid registration */
  337.         x0 -= 0.5 * dx;
  338.         y0 -= 0.5 * dy;
  339.         map_clip_on (-1, -1, -1, 3);
  340.     }
  341.     x_side = dx * r_head.nx;
  342.     y_side = dy * r_head.ny;
  343.     
  344.     if (gmtdefs.verbose) fprintf (stderr, "grdimage: Creating PostScript image ");
  345.  
  346.     if (gmt_gray) for (k = 0; gmt_b_and_w && k < nm; k++) if (!(bitimage_8[k] == 0 || bitimage_8[k] == 255)) gmt_b_and_w = FALSE;
  347.     
  348.     if (gmt_b_and_w) {    /* Can get away with 1 bit image */
  349.         int nx8, shift, byte, b_or_w, k8;
  350.         unsigned char *bit;
  351.         
  352.         if (gmtdefs.verbose) fprintf (stderr, "[1-bit B/W image]\n");
  353.         nx8 = ceil (r_head.nx / 8.0);
  354.         bit = (unsigned char *) memory (CNULL, (int)(nx8 * r_head.ny), sizeof (char), "grdimage");
  355.         
  356.         for (j = k = k8 = 0; j < r_head.ny; j++) {
  357.             shift = byte = 0;
  358.             for (i = 0; i < r_head.nx; i++, k++) {
  359.                 b_or_w = (bitimage_8[k] == 255);
  360.                 byte |= b_or_w;
  361.                 shift++;
  362.                 if (shift == 8) {    /* Time to dump out byte */
  363.                     bit[k8++] = byte;
  364.                     byte = shift = 0;
  365.                 }
  366.                 else
  367.                     byte <<= 1;
  368.             }
  369.             if (shift) {
  370.                 byte |= 1;
  371.                 shift++;
  372.                 while (shift < 8) {
  373.                     byte <<= 1;
  374.                     byte |= 1;
  375.                     shift++;
  376.                 }
  377.                 bit[k8++] = byte;
  378.             }
  379.         }
  380.         free ((char *)bitimage_8);
  381.         
  382.         x_side = nx8 * 8.0 * dx;
  383.         if ((nx8 * 8) > r_head.nx && grid_type == 1) {    /* must clip since 1-bit image is wider */
  384.             map_clip_on (-1, -1, -1, 3);
  385.             mapped = TRUE;
  386.         }
  387.         ps_image (x0, y0, x_side, y_side, bit, nx8, r_head.ny, 1);
  388.         free ((char *)bit);
  389.     }
  390.     else if (gmt_gray || monochrome) {
  391.         if (gmtdefs.verbose) fprintf (stderr, "[8-bit grayshade image]\n");
  392.         ps_image (x0, y0, x_side, y_side, bitimage_8, r_head.nx, r_head.ny, 8);
  393.         free ((char *)bitimage_8);
  394.     }
  395.     else {
  396.         if (gmtdefs.verbose) fprintf (stderr, "24-bit [%s]\n", c_method[gmtdefs.color_image]);
  397.         if (gmtdefs.color_image == 2)
  398.             fprintf (stderr, "%s: GMT Warning: PostScript output may only be processed with psto24!\n", gmt_program);
  399.         color_image (x0, y0, x_side, y_side, bitimage_24, r_head.nx, r_head.ny);
  400.         free ((char *)bitimage_24);
  401.     }
  402.         
  403.     if (mapped || grid_type == 0) map_clip_off();
  404.     
  405.     if (frame_info.plot) {
  406.         ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  407.         map_basemap ();
  408.         ps_setpaint (0, 0, 0);
  409.     }
  410.     ps_plotend (gmtdefs.last_page);
  411.     
  412.     gmt_end (argc, argv);
  413. }
  414.