home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grd2cpt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-08  |  10.2 KB  |  330 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system: @(#)grd2cpt.c    2.12  08 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.  * grd2cpt reads a 2d binary gridded grdfile and creates a continous-color-
  9.  * pallette cpt file, with a non-linear histogram-equalized mapping between
  10.  * hue and data value.  (The linear mapping can be made with makecpt.)
  11.  *
  12.  * Creates a cumulative distribution function f(z) describing the data
  13.  * in the grdfile.  f(z) is sampled at z values supplied by the user
  14.  * [with -S option] or guessed from the sample mean and standard deviation.
  15.  * f(z) is then found by looping over the grd array for each z and counting
  16.  * data values <= z.  Once f(z) is found then
  17.  * hue(z) = hue_min + (hue_max - hue_min) * f(z), where hue_min and
  18.  * hue_max can be set by user or default to 300, 0.
  19.  *
  20.  * The goal here is to have something that gets the inexperienced user
  21.  * started with a .cpt file that will do nice image processing.
  22.  *
  23.  * Author:    Walter H. F. Smith
  24.  * Date:    12-JAN-1994
  25.  * 
  26.  */
  27.  
  28. #include "gmt.h"
  29.  
  30. struct CDF_CPT {
  31.     double    z;    /* Data value  */
  32.     double    f;    /* Cumulative distribution function f(z)  */
  33. } *cdf_cpt = NULL;
  34.  
  35. main (argc, argv)
  36. int argc;
  37. char **argv;
  38. {
  39.  
  40.     int i, j, nxy, nfound, ngood, ncdf, pad[4];
  41.     int    lastr, lastg, lastb, thisr, thisg, thisb;
  42.     
  43.     BOOLEAN error = FALSE, set_limits = FALSE, set_z_vals = FALSE, monochrome = FALSE;
  44.     
  45.     double min_limit, max_limit, min_hue, max_hue;
  46.     double    lasthue, thishue, sat, val;
  47.     double    z_start, z_stop, z_inc;
  48.     double    mean, sd, min_gray, max_gray;
  49.     
  50.     char *grdfile;
  51.     
  52.     struct GRD_HEADER grd;
  53.  
  54.     float    *zdata;
  55.     
  56.     min_hue = 300.0;
  57.     max_hue = 0.0;
  58.     min_gray = 0.0;
  59.     max_gray = 1.0;
  60.     grdfile = CNULL;
  61.     pad[3] = pad[2] = pad[1] = pad[0] = 0;
  62.     
  63.     argc = gmt_begin (argc, argv);
  64.  
  65.     for (i = 1; !error && i < argc; i++) {
  66.         if (argv[i][0] == '-') {
  67.             switch (argv[i][1]) {
  68.             
  69.                 /* Common parameters */
  70.             
  71.                 case 'V':
  72.                 case '\0':
  73.                     error += get_common_args (argv[i], 0, 0, 0, 0);
  74.                     break;
  75.  
  76.                 /* Supplemental parameters */
  77.                 
  78.                 case 'C':
  79.                     if (sscanf(&argv[i][2], "%lf/%lf", &min_hue, &max_hue) != 2) {
  80.                         fprintf(stderr,"%s: GMT SYNTAX ERROR -C option:  Cannot decode hues\n", gmt_program);
  81.                         error++;
  82.                     }
  83.                     else {
  84.                         if (min_hue < 0.0 || min_hue >= 360.0 || max_hue < 0.0 || max_hue > 360.0) {
  85.                             fprintf(stderr,"%s: GMT SYNTAX ERROR -C option:  Hues our of 0-360 range\n", gmt_program);
  86.                             error++;
  87.                         }
  88.                     }
  89.                     break;
  90.                 case 'L':
  91.                     if (sscanf(&argv[i][2], "%lf/%lf", &min_limit, &max_limit) != 2) {
  92.                         fprintf(stderr,"%s: GMT SYNTAX ERROR -L option:  Cannot decode limits\n", gmt_program);
  93.                         error++;
  94.                     }
  95.                     else {
  96.                         if (min_limit >= max_limit) {
  97.                             fprintf(stderr,"%s: GMT SYNTAX ERROR -L option:  min_limit must be less than max_limit.\n", gmt_program);
  98.                             error++;
  99.                         }
  100.                     }
  101.                     if (!error) set_limits = TRUE;
  102.                     break;
  103.                 case 'M':
  104.                     monochrome = TRUE;
  105.                     if (argv[i][2] && sscanf(&argv[i][2], "%lf/%lf", &min_gray, &max_gray) != 2) {
  106.                         fprintf(stderr,"%s: GMT SYNTAX ERROR -M option:  Cannot decode gray values\n", gmt_program);
  107.                         error++;
  108.                     }
  109.                     break;
  110.                 case 'S':
  111.                     if (sscanf(&argv[i][2], "%lf/%lf/%lf", &z_start, &z_stop, &z_inc) != 3) {
  112.                         fprintf(stderr,"%s: GMT SYNTAX ERROR -S option:  Cannot decode values\n", gmt_program);
  113.                         error++;
  114.                     }
  115.                     else {
  116.                         if (z_stop <= z_start || z_inc <= 0.0) {
  117.                             fprintf(stderr,"%s: GMT SYNTAX ERROR -S option:  Bad arguments\n", gmt_program);
  118.                             error++;
  119.                         }
  120.                     }
  121.                     if (!error) set_z_vals = TRUE;
  122.                     break;
  123.  
  124.                 default:
  125.                     error = TRUE;
  126.                     gmt_default_error (argv[i][1]);
  127.                     break;
  128.             }
  129.         }
  130.         else {
  131.             grdfile = argv[i];
  132.             if (read_grd_info (grdfile, &grd)) {
  133.                 fprintf (stderr, "grd2cpt: Error opening file %s\n", grdfile);
  134.                 error++;
  135.             }
  136.         }
  137.     }
  138.     
  139.     if (argc == 1 || gmt_quick) {
  140.         fprintf (stderr,"grd2cpt %s - Make a histogram-equalized color palette table from a grdfile\n\n", GMT_VERSION);
  141.         fprintf (stderr, "usage: grd2cpt <grdfile> [-C<min_hue>/<max_hue>] [-L<min_limit>/<max_limit>]\n");
  142.         fprintf (stderr, "    [-M[<mingray>/<maxgray>]] [-S<z_start>/<z_stop>/<z_inc>] [-V]\n");
  143.         
  144.         if (gmt_quick) exit(-1);
  145.         
  146.         fprintf (stderr, "    <grdfile> is the name of the 2-D binary data set\n");
  147.         fprintf (stderr, "\n\tOPTIONS:\n");
  148.         explain_option ('H');
  149.         fprintf (stderr, "    -C Set the hues to use at min and max data values.  [300/0].\n");
  150.         fprintf (stderr, "    -L Limit the range of the data [Default uses actual min,max of data].\n");
  151.         fprintf (stderr, "    -M Set Monochrome (gray).  Optionally set the gray levels for min,max data,\n");
  152.         fprintf (stderr, "             using a scale where 0 = black and 1 = white.  [0/1]\n");
  153.         fprintf (stderr, "    -S Sample points should Step from z_start to z_stop by z_inc [Default guesses some values].\n");
  154.         explain_option ('V');
  155.         exit(-1);
  156.     }
  157.     
  158.     if (error) exit (-1);
  159.     
  160.     nxy = grd.nx * grd.ny;
  161.     zdata = (float *) memory (CNULL, nxy, sizeof (float), "grd2cpt");
  162.  
  163.     if (read_grd (grdfile, &grd, zdata, 0.0, 0.0, 0.0, 0.0, pad, FALSE) ) {
  164.         fprintf(stderr,"grd2cpt:  Error reading grdfile.\n");
  165.         free((char *)zdata);
  166.         exit(-1);
  167.     }
  168.  
  169.     /* Loop over the file and find NaNs.  If set limits, may create more NaNs  */
  170.     nfound = 0;
  171.     mean = 0.0;
  172.     sd = 0.0;
  173.     if (set_limits) {
  174.         /* Loop over the grdfile, and set anything outside the limiting values to NaN.  */
  175.  
  176.         grd.z_min = min_limit;
  177.         grd.z_max = max_limit;
  178.         for (i = 0; i < nxy; i++) {
  179.             if (bad_float ((double)zdata[i]))
  180.                 nfound++;
  181.             else {
  182.                 if (zdata[i] < min_limit || zdata[i] > max_limit) {
  183.                     nfound++;
  184.                     zdata[i] = gmt_NaN;
  185.                 }
  186.                 else {
  187.                     mean += zdata[i];
  188.                     sd += zdata[i] * zdata[i];
  189.                 }
  190.             }
  191.         }
  192.     }
  193.     else {
  194.         min_limit = grd.z_max;    /* This is just to double check grd.z_min, grd.z_max  */
  195.         max_limit = grd.z_min;
  196.         for (i = 0; i < nxy; i++) {
  197.             if (bad_float ((double)zdata[i]))
  198.                 nfound++;
  199.             else {
  200.                 if (zdata[i] < min_limit) min_limit = zdata[i];
  201.                 if (zdata[i] > max_limit) max_limit = zdata[i];
  202.                 mean += zdata[i];
  203.                 sd += zdata[i] * zdata[i];
  204.             }
  205.         }
  206.         grd.z_min = min_limit;
  207.         grd.z_max = max_limit;
  208.     }
  209.     ngood = nxy - nfound;    /* This is the number of non-NaN points for the cdf function  */
  210.     mean /= ngood;
  211.     sd /= ngood;
  212.     sd = sqrt(sd - mean*mean);
  213.     if (gmtdefs.verbose) fprintf(stderr,"grd2cpt:  Mean and S.D. of data are %.8lg %.8lg\n", mean, sd);
  214.  
  215.     /* Now the zdata are ready.  Decide how to make steps in z.  */
  216.     if (set_z_vals) {
  217.         ncdf =  (grd.z_min < z_start) ? 1 : 0;
  218.         ncdf += floor((z_stop - z_start)/z_inc) + 1;
  219.         if (grd.z_max > z_stop) ncdf++;
  220.         cdf_cpt = (struct CDF_CPT *)memory(CNULL, ncdf, sizeof(struct CDF_CPT), "grd2cpt");
  221.         if (grd.z_min < z_start) {
  222.             cdf_cpt[0].z = grd.z_min;
  223.             cdf_cpt[1].z = z_start;
  224.             i = 2;
  225.         }
  226.         else {
  227.             cdf_cpt[0].z = z_start;
  228.             i = 1;
  229.         }
  230.         j = (grd.z_max > z_stop) ? ncdf - 1 : ncdf;
  231.         while (i < j) {
  232.             cdf_cpt[i].z = cdf_cpt[i-1].z + z_inc;
  233.             i++;
  234.         }
  235.         if (j == ncdf-1) cdf_cpt[j].z = grd.z_max;
  236.     }
  237.     else {
  238.         /* This is completely ad-hoc.  It chooses z based on steps
  239.             of 0.1 for a Gaussian CDF:  */
  240.         ncdf = 11;
  241.         cdf_cpt = (struct CDF_CPT *)memory(CNULL, ncdf, sizeof(struct CDF_CPT), "grd2cpt");
  242.         cdf_cpt[0].z = grd.z_min;
  243.         cdf_cpt[1].z = mean - 1.28155 * sd;
  244.         cdf_cpt[2].z = mean - 0.84162 * sd;
  245.         cdf_cpt[3].z = mean - 0.52440 * sd;
  246.         cdf_cpt[4].z = mean - 0.25335 * sd;
  247.         cdf_cpt[5].z = mean;
  248.         cdf_cpt[6].z = mean + 0.25335 * sd;
  249.         cdf_cpt[7].z = mean + 0.52440 * sd;
  250.         cdf_cpt[8].z = mean + 0.84162 * sd;
  251.         cdf_cpt[9].z = mean + 1.28155 * sd;
  252.         cdf_cpt[10].z = grd.z_max;
  253.     }
  254.  
  255.     /* Get here when we are ready to go.  cdf_cpt[].z contains the sample points.  */
  256.     for (j = 0; j < ncdf; j++) {
  257.         if (cdf_cpt[j].z <= grd.z_min)
  258.             cdf_cpt[j].f = 0.0;
  259.         else if (cdf_cpt[j].z >= grd.z_max)
  260.             cdf_cpt[j].f = 1.0;
  261.         else {
  262.             nfound = 0;
  263.             for (i = 0; i < nxy; i++) {
  264.                 if (!bad_float((double)zdata[i]) && zdata[i] <= cdf_cpt[j].z) nfound++;
  265.             }
  266.             cdf_cpt[j].f = (double)(nfound-1)/(double)(ngood-1);
  267.         }
  268.         if (gmtdefs.verbose) fprintf (stderr, "grd2cpt: z = %.8lg and CDF(z) = %.8lg\n",
  269.             cdf_cpt[j].z, cdf_cpt[j].f);
  270.     }
  271.  
  272.     /* Now the cdf function has been found.  We get hsv and then covert to rgb as needed  */
  273.     if (monochrome) {
  274.         if (gmtdefs.color_model) {
  275.             /* Make a monochrome HSV file.  */
  276.             for (j = 1; j < ncdf; j++) {
  277.                 printf("%.6lg\t0\t0\t%.6lg\t%.6lg\t0\t0\t%.6lg\n",
  278.                     cdf_cpt[j-1].z, min_gray + (max_gray - min_gray) * cdf_cpt[j-1].f,
  279.                         cdf_cpt[j].z, min_gray + (max_gray - min_gray) * cdf_cpt[j].f);
  280.             }
  281.         }
  282.         else {
  283.             /* Make an r,g,b gray file  */
  284.             lasthue = min_gray + (max_gray - min_gray) * cdf_cpt[0].f;
  285.             lastr = floor(255.0 * lasthue);
  286.             for (j = 1; j < ncdf; j++) {
  287.                 thishue = min_gray + (max_gray - min_gray) * cdf_cpt[j].f;
  288.                 thisr = floor(255.0 * thishue);
  289.                 printf("%.6lg\t%d\t%d\t%d\t%.6lg\t%d\t%d\t%d\n",
  290.                     cdf_cpt[j-1].z, lastr, lastr, lastr, cdf_cpt[j].z, thisr, thisr, thisr);
  291.                 lasthue = thishue;
  292.                 lastr = thisr;
  293.             }
  294.         }
  295.     }
  296.     else {
  297.         if (gmtdefs.color_model) {
  298.             sat = gmtdefs.hsv_min_saturation;
  299.             val = gmtdefs.hsv_max_value;
  300.         }
  301.         else {
  302.             sat = val = 1.0;
  303.         }
  304.         lasthue = min_hue;
  305.         if (!gmtdefs.color_model) gmt_hsv_to_rgb (&lastr, &lastg, &lastb, lasthue, sat, val);
  306.  
  307.         for (j = 1; j < ncdf; j++) {
  308.             thishue = min_hue + (max_hue - min_hue) * cdf_cpt[j].f;
  309.             if (gmtdefs.color_model) {
  310.                 printf("%.6lg\t%.6lg\t%.6lg\t%.6lg\t%.6lg\t%.6lg\t%.6lg\t%.6lg\n",
  311.                     cdf_cpt[j-1].z, lasthue, sat, val, cdf_cpt[j].z, thishue, sat, val);
  312.             }
  313.             else {
  314.                 gmt_hsv_to_rgb (&thisr, &thisg, &thisb, thishue, sat, val);
  315.                 printf("%.6lg\t%d\t%d\t%d\t%.6lg\t%d\t%d\t%d\n",
  316.                 cdf_cpt[j-1].z, lastr, lastg, lastb, cdf_cpt[j].z, thisr, thisg, thisb);
  317.                 lastr = thisr;
  318.                 lastg = thisg;
  319.                 lastb = thisb;
  320.             }
  321.             lasthue = thishue;
  322.         }
  323.     }
  324.     
  325.     free ((char *)cdf_cpt);
  326.     free ((char *)zdata);
  327.     
  328.     gmt_end (argc, argv);
  329. }
  330.