home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / blockmedian.c < prev    next >
C/C++ Source or Header  |  1995-04-28  |  12KB  |  432 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)blockmedian.c    2.21  4/28/95
  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. /*
  9.    blockmedian.c
  10.    Takes lon, lat, data, [weight] on stdin or file and writes out one value
  11.    per cell, where cellular region is bounded by West East South North and
  12.    cell dimensions are delta_x, delta_y.
  13.       
  14.    Author:     Walter H. F. Smith
  15.    Date:    28 June, 1988
  16.    Modified    26 April, 1991-1995 for gmt v2.0 by whfs smith;
  17.         added dynamic memory allocation.
  18.    Modified:    3 Jan 1995 by PW for gmt 3.0
  19. */
  20.  
  21. #include "gmt.h"
  22.  
  23. int    compare_x();
  24. int    compare_y();
  25. int    compare_index_z();
  26. int    median_output();
  27.  
  28. struct DATA {
  29.     int    i;
  30.     float    w;
  31.     float    x;
  32.     float    y;
  33.     float    z;
  34. }    *data;
  35.  
  36.  
  37. main (argc, argv)
  38. int argc;
  39. char **argv;
  40. {    
  41.     
  42.     BOOLEAN    error, weighted, offset, report, nofile = TRUE, done = FALSE, first = TRUE;
  43.     BOOLEAN b_in = FALSE, b_out = FALSE, single_precision = FALSE, more, skip;
  44.     
  45.     FILE *fp = NULL;
  46.     
  47.     double    west, east, south, north, delta_x, delta_y, del_off;
  48.     double    *in, out[4], idx, idy;
  49.     
  50.     int    i, n_x, n_y, ix, iy, index, first_in_cell, first_in_new_cell, fno, n_files = 0, n_args;
  51.     int    n_lost, n_read, n_pitched, n_cells_filled, n_alloc, x, y, n_expected_fields, n_fields, n_out;
  52.     
  53.     char    modifier, buffer[512];
  54.  
  55.     argc = gmt_begin (argc, argv);
  56.     
  57.     west = east = south = north = delta_x = delta_y = 0.0;
  58.     del_off = 0.5;
  59.     error = weighted = offset = report = FALSE;
  60.  
  61.     for (i = 1; i < argc; i++) {
  62.         if (argv[i][0] == '-') {
  63.             switch (argv[i][1]) {
  64.               
  65.                 /* Common parameters */
  66.                       
  67.                 case 'H':
  68.                 case 'R':
  69.                 case 'V':
  70.                 case ':':
  71.                 case '\0':
  72.                                       error += get_common_args (argv[i], &west, &east, &south, &north);
  73.                                       break;
  74.                               
  75.                 /* Supplemental parameters */
  76.                               
  77.                 case 'b':
  78.                     single_precision = !(argv[i][strlen(argv[i])-1] == 'd');
  79.                     if (argv[i][2] == 'i')
  80.                         b_in = TRUE;
  81.                     else if (argv[i][2] == 'o')
  82.                         b_out = TRUE;
  83.                     else
  84.                         b_in = b_out = TRUE;
  85.                     break;
  86.                     
  87.                 case 'I':
  88.                     gmt_getinc (&argv[i][2], &delta_x, &delta_y);
  89.                     break;
  90.                 case 'N':
  91.                     offset = TRUE;
  92.                     break;
  93.                 case 'W':
  94.                     if ( (modifier = argv[i][2]) == 'i' || modifier == 'I') {
  95.                         weighted = TRUE;
  96.                         report = FALSE;
  97.                     }
  98.                     else if (modifier == 'O' || modifier == 'o') {
  99.                         report = TRUE;
  100.                         weighted = FALSE;
  101.                     }
  102.                     else
  103.                         weighted = report = TRUE;
  104.                     break;
  105.                     
  106.                 default:
  107.                     error = TRUE;
  108.                                         gmt_default_error (argv[i][1]);
  109.                     break;
  110.             }
  111.         }
  112.         else
  113.             n_files++;
  114.     }
  115.  
  116.     if (argc == 1 || gmt_quick) {
  117.         fprintf (stderr, "blockmedian %s - Block averaging by L1 norm\n\n", GMT_VERSION);
  118.         fprintf (stderr, "usage: blockmedian [infile(s)] -I<xinc[m]>[/<yinc>[m]] -R<west/east/south/north>\n");
  119.         fprintf (stderr, "\t[-H] [-N] [-V] [-W[i][o] ] [-:] [-b[i|o]i[d]]\n\n");
  120.         
  121.         if (gmt_quick) exit (-1);
  122.         
  123.         fprintf (stderr, "\t-I sets Increment of the grid; enter xinc, optionally xinc/yinc.\n");
  124.         fprintf (stderr, "\t   Default is yinc = xinc.  Append an m [or s] to xinc or yinc to indicate minutes [or seconds],\n");
  125.         fprintf (stderr, "\t   e.g., -I10m/5m grids longitude every 10 minutes, latitude every 5 minutes.\n");
  126.         explain_option ('R');
  127.         fprintf (stderr, "\n\tOPTIONS:\n");
  128.         explain_option ('H');
  129.         fprintf (stderr, "\t-N Offsets registration so block edges are on gridlines.\n");
  130.         explain_option ('V');
  131.         fprintf (stderr, "\t-W sets Weight options.  -WI reads Weighted Input (4 cols: x,y,z,w) but writes only (x,y,z) Output.\n");
  132.         fprintf (stderr, "\t   -WO reads unWeighted Input (3 cols: x,y,z) but reports sum (x,y,z,w) Output.\n");
  133.         fprintf (stderr, "\t   -W with no modifier has both weighted Input and Output; Default is no weights used.\n");
  134.         explain_option (':');
  135.         fprintf (stderr, "\t-b means binary (double) i/o.  Append i or o if only input OR output is binary\n");
  136.         fprintf (stderr, "\t   Append d for double precision [Default is single precision].\n");        explain_option ('.');
  137.         exit (-1);
  138.     }
  139.     
  140.     if (!project_info.region_supplied) {
  141.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", gmt_program);
  142.         error++;
  143.     }
  144.     if (delta_x <= 0.0 || delta_y <= 0.0) {
  145.         fprintf (stderr, "%s: GMT SYNTAX ERROR -I option.  Must specify positive increment(s)\n", gmt_program);
  146.         error = TRUE;
  147.     }
  148.     if (b_in && gmtdefs.io_header) {
  149.         fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", gmt_program);
  150.         error++;
  151.     }
  152.     
  153.     if (error) exit (-1);
  154.  
  155.     if (b_in) gmt_input = (single_precision) ? bin_float_input : bin_double_input;
  156.     if (b_out) gmt_output = (single_precision) ? bin_float_output : bin_double_output;
  157.  
  158.     idx = 1.0 / delta_x;
  159.     idy = 1.0 / delta_y;
  160.     n_x = rint((east - west) * idx) + 1;
  161.     n_y = rint((north - south) * idy) + 1;
  162.     
  163.     if (offset) {
  164.         n_x--;
  165.         n_y--;
  166.         del_off = 0.0;
  167.     }
  168.     
  169.     if (gmtdefs.verbose) fprintf (stderr, "W: %.5lf E: %.5lf S: %.5lf N: %.5lf nx: %d ny: %d\n",
  170.         west, east, south, north, n_x, n_y);
  171.  
  172.     n_read = 0;
  173.     n_pitched = 0;
  174.         n_alloc = GMT_CHUNK;
  175.         data = (struct DATA *) memory (CNULL, n_alloc, sizeof(struct DATA), "blockmedian");
  176.     
  177.     gmt_data[3] = 1.0;    /* Default weight if weights not supplied  */
  178.  
  179.     /* Read the input data  */
  180.     
  181.     x = (gmtdefs.xy_toggle) ? 1 : 0;        y = 1 - x;              /* Set up which columns have x and y */
  182.     n_expected_fields = (weighted) ? 4 : 3;
  183.     
  184.     if (n_files > 0)
  185.         nofile = FALSE;
  186.     else
  187.         n_files = 1;
  188.     n_args = (argc > 1) ? argc : 2;
  189.     
  190.     for (fno = 1; !done && fno < n_args; fno++) {    /* Loop over input files, if any */
  191.         if (!nofile && argv[fno][0] == '-') continue;
  192.         
  193.         if (nofile) {    /* Just read standard input */
  194.             fp = stdin;
  195.             done = TRUE;
  196.         }
  197.         else if ((fp = fopen (argv[fno], "r")) == NULL) {
  198.             fprintf (stderr, "blockmedian: Cannot open file %s\n", argv[fno]);
  199.             continue;
  200.         }
  201.  
  202.         if (!nofile && gmtdefs.verbose) fprintf (stderr, "blockmedian: Working on file %s\n", argv[fno]);
  203.         
  204.         if (gmtdefs.io_header) {
  205.             for (i = 0; i < gmtdefs.n_header_recs; i++) {
  206.                 fgets (buffer, 512, fp);
  207.                 buffer[strlen(buffer)-1] = 0;
  208.                 if (first && !b_out) (report && !(weighted)) ? printf ("%s weights\n", buffer) : printf ("%s\n", buffer);
  209.             }
  210.             first = FALSE;
  211.         }
  212.  
  213.         more = ((n_fields = gmt_input (fp,  n_expected_fields, &in)) == n_expected_fields);
  214.                 
  215.         while (more) {
  216.     
  217.             skip = FALSE;
  218.             
  219.             if (bad_float ((float)in[2])) skip = TRUE;    /* Skip when z = NaN */
  220.  
  221.             if (!skip) {    /* Check if point is inside area */
  222.                 n_read++;
  223.         
  224.                 ix = floor(((in[x] - west) * idx) + del_off);
  225.                 if ( ix < 0 || ix >= n_x ) skip = TRUE;
  226.                 iy = floor(((in[y] - south) * idy) + del_off);
  227.                 if ( iy < 0 || iy >= n_y ) skip = TRUE;
  228.             }
  229.         
  230.             if (!skip) {
  231.                 index = iy * n_x + ix;
  232.         
  233.                 data[n_pitched].i = index;
  234.                 data[n_pitched].w = in[3];
  235.                 data[n_pitched].x = in[x];
  236.                 data[n_pitched].y = in[y];
  237.                 data[n_pitched].z = in[2];
  238.         
  239.                 n_pitched++;
  240.                 if (n_pitched == n_alloc) {
  241.                     n_alloc += GMT_CHUNK;
  242.                     data = (struct DATA *) memory ((char *)data, n_alloc, sizeof(struct DATA), "blockmedian");
  243.                 }
  244.             }
  245.             
  246.             more = ((n_fields = gmt_input (fp,  n_expected_fields, &in)) == n_expected_fields);
  247.         }
  248.         if (fp != stdin) fclose(fp);
  249.         
  250.         if (n_fields == -1) n_fields = 0;    /* Ok to get EOF */
  251.         if (n_fields%n_expected_fields) {    /* Got garbage or multiple segment subheader */
  252.             fprintf (stderr, "%s:  Cannot read X Y Z [W] at line %d, aborts\n", gmt_program, n_read);
  253.             exit (-1);
  254.         }
  255.     }
  256.     
  257.     data = (struct DATA *) memory ((char *)data, n_pitched, sizeof(struct DATA), "blockmedian");
  258.     n_lost = n_read - n_pitched;
  259.     if (gmtdefs.verbose) fprintf(stderr,"N_read: %d N_used: %d N_outside_area: %d\n",
  260.         n_read,n_pitched,n_lost);
  261.  
  262.  
  263.     /* Ready to go. */
  264.     
  265.     n_out = (report) ? 4 : 3;
  266.     
  267.     /* Sort on index and Z value */
  268.     
  269.     qsort((char *)data, n_pitched, sizeof (struct DATA), compare_index_z);
  270.     
  271.     /* Find n_in_cell and write appropriate output  */
  272.  
  273.     first_in_cell = 0;
  274.     n_cells_filled = 0;
  275.     while (first_in_cell < n_pitched) {
  276.         out[3] = data[first_in_cell].w;
  277.         first_in_new_cell = first_in_cell + 1;
  278.         while ( (first_in_new_cell < n_pitched) && (data[first_in_new_cell].i == data[first_in_cell].i) ) {
  279.             out[3] += data[first_in_new_cell].w;
  280.             first_in_new_cell++;
  281.         }
  282.         median_output(first_in_cell, first_in_new_cell, out[3], &out[x], &out[y], &out[2]);
  283.         
  284.         gmt_output (stdout, n_out, out);
  285.         
  286.         n_cells_filled++;
  287.         first_in_cell = first_in_new_cell;
  288.     }
  289.  
  290.     if (gmtdefs.verbose) fprintf(stderr,"N_cells_filled: %d\n", n_cells_filled);
  291.  
  292.     free((char *)data);
  293.     
  294.         gmt_end (argc, argv);
  295. }
  296.  
  297. int median_output(first_in_cell, first_in_new_cell, weight_sum, xx, yy, zz)
  298.     int    first_in_cell, first_in_new_cell;
  299.     double    *xx, *yy, *zz, weight_sum;
  300. {
  301.     double    weight_half, weight_count;
  302.     int    index, n_in_cell;
  303.     
  304.     weight_half  = 0.5 * weight_sum;
  305.     n_in_cell = first_in_new_cell - first_in_cell;
  306.     
  307.     /* Data are already sorted on z  */
  308.  
  309.     index = first_in_cell;
  310.     weight_count = data[first_in_cell].w;
  311.     while (weight_count < weight_half) {
  312.         index++;
  313.         weight_count += data[index].w;
  314.     }
  315.     if ( weight_count == weight_half )
  316.         *zz = 0.5 * (data[index].z + data[index + 1].z);
  317.     else
  318.         *zz = data[index].z;
  319.     
  320.     /* Now do x and y */
  321.     
  322.     if (n_in_cell > 2)
  323.         qsort((char *)&data[first_in_cell], n_in_cell, sizeof (struct DATA), compare_x);
  324.     index = first_in_cell;
  325.     weight_count = data[first_in_cell].w;
  326.     while (weight_count < weight_half) {
  327.         index++;
  328.         weight_count += data[index].w;
  329.     }
  330.     if ( weight_count == weight_half )
  331.         *xx = 0.5 * (data[index].x + data[index + 1].x);
  332.     else
  333.         *xx = data[index].x;
  334.  
  335.     if (n_in_cell > 2)
  336.         qsort((char *)&data[first_in_cell], n_in_cell, sizeof (struct DATA), compare_y);
  337.     index = first_in_cell;
  338.     weight_count = data[first_in_cell].w;
  339.     while (weight_count < weight_half) {
  340.         index++;
  341.         weight_count += data[index].w;
  342.     }
  343.     if ( weight_count == weight_half )
  344.         *yy = 0.5 * (data[index].y + data[index + 1].y);
  345.     else
  346.         *yy = data[index].y;
  347.  
  348.     return;
  349. }
  350.  
  351. int    compare_index_z(point_1, point_2)
  352.  
  353.     struct DATA    *point_1, *point_2;
  354. {
  355.     int    index_1, index_2;
  356.     float    data_1, data_2;
  357.     
  358.     
  359.     index_1 = point_1->i;
  360.     index_2 = point_2->i;
  361.  
  362.     if (index_1 < index_2)
  363.         return (-1);
  364.     else if (index_1 > index_2)
  365.         return (1);
  366.     else {
  367.         data_1 = point_1->z;
  368.         data_2 = point_2->z;
  369.         if (data_1 < data_2)
  370.             return (-1);
  371.         else if (data_1 > data_2)
  372.             return (1);
  373.         else
  374.             return (0);
  375.     }
  376. }
  377. int    compare_x(point_1, point_2)
  378.  
  379.     struct DATA    *point_1, *point_2;
  380. {
  381.     int    index_1, index_2;
  382.     float    x_1, x_2;
  383.     
  384.     
  385.     index_1 = point_1->i;
  386.     index_2 = point_2->i;
  387.  
  388.     if (index_1 < index_2)
  389.         return (-1);
  390.     else if (index_1 > index_2)
  391.         return (1);
  392.     else {
  393.         x_1 = point_1->x;
  394.         x_2 = point_2->x;
  395.         if (x_1 < x_2)
  396.             return (-1);
  397.         else if (x_1 > x_2)
  398.             return (1);
  399.         else
  400.             return (0);
  401.     }
  402. }
  403.  
  404. int    compare_y(point_1, point_2)
  405.  
  406.     struct DATA    *point_1, *point_2;
  407. {
  408.     int    index_1, index_2;
  409.     float    y_1, y_2;
  410.     
  411.     
  412.     index_1 = point_1->i;
  413.     index_2 = point_2->i;
  414.  
  415.     if (index_1 < index_2)
  416.         return (-1);
  417.     else if (index_1 > index_2)
  418.         return (1);
  419.     else {
  420.         y_1 = point_1->y;
  421.         y_2 = point_2->y;
  422.         if (y_1 < y_2)
  423.             return (-1);
  424.         else if (y_1 > y_2)
  425.             return (1);
  426.         else
  427.             return (0);
  428.     }
  429. }
  430.  
  431.  
  432.