home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdfilter.c < prev    next >
C/C++ Source or Header  |  1995-04-28  |  16KB  |  481 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdfilter.c    2.22  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.  * grdfilter.c  reads a grdfile and creates filtered grd file
  9.  *
  10.  * user selects boxcar, gaussian, cosine arch, median, or mode filter
  11.  * user selects distance scaling as appropriate for map work, etc.
  12.  *
  13.  * Author:  W.H.F. Smith
  14.  * Date:     16 Aug 88
  15.  *
  16.  * Modified:    27 Sep 88 by W.H.F. Smith, to use new median routine.
  17.  *
  18.  * Updated:    PW: 13-Jun-1991 to v2.0
  19.  *        PW: 13-Jul-1992 to actually do v2 i/o.
  20.  *        PW: 15-Jul-1992 to handle arbitrary new -R -I
  21.  *        PW: 3-Jan-1995 to offer mode-filter (from filter1d)
  22.  *
  23.  */
  24.  
  25. #include "gmt.h"
  26.  
  27. double    get_wt();
  28. int    median();
  29. int    mode();
  30.  
  31. int    *i_origin;
  32. float    *input, *output;
  33. double    *weight, *work_array, *x_shift;
  34. double    deg2km;
  35.  
  36. main (argc, argv)
  37. int argc;
  38. char **argv; {
  39.  
  40.     int    nx_out, ny_out, nx_fil, ny_fil, n_in_median, n_nan = 0;
  41.     int    x_half_width, y_half_width, j_origin, i_out, j_out;
  42.     int    i_in, j_in, ii, jj, i, ij_in, ij_out, ij_wt, effort_level;
  43.     int    distance_flag, filter_type, one_or_zero = 1, dummy[4];
  44.     
  45.     BOOLEAN    error, new_range, new_increment, fast_way, shift, slow;
  46.     
  47.     double    west_new, east_new, south_new, north_new, dx_new, dy_new, offset;
  48.     double    filter_width, x_scale, y_scale, x_width, y_width;
  49.     double    x_out, y_out, wt_sum, value, last_median, this_median, xincnew2, yincnew2;
  50.     double    xincold2, yincold2, y_shift, x_fix = 0.0, y_fix = 0.0;
  51.     
  52.     char    *fin = NULL, *fout = NULL;
  53.     
  54.     struct    GRD_HEADER h;
  55.  
  56.     argc = gmt_begin (argc, argv);
  57.     
  58.     deg2km = 0.001 * 2.0 * M_PI * gmtdefs.ellipse[gmtdefs.ellipsoid].eq_radius / 360.0;
  59.     error = new_range = new_increment = FALSE;
  60.     fin = fout = NULL;
  61.     filter_width = dx_new = dy_new = west_new = east_new = 0.0;
  62.     filter_type = distance_flag = -1;
  63.     dummy[3] = dummy[2] = dummy[1] = dummy[0] = 0;
  64.     
  65.     for (i = 1; i < argc; i++) {
  66.         if (argv[i][0] == '-') {
  67.             switch (argv[i][1]) {
  68.                 /* Common parameters */
  69.             
  70.                 case 'R':
  71.                 case 'V':
  72.                 case '\0':
  73.                     error += get_common_args (argv[i], &west_new, &east_new, &south_new, &north_new);
  74.                     break;
  75.                 
  76.                 case 'F':
  77.                     switch (argv[i][2]) {
  78.                         case 'b':
  79.                             filter_type = 0;
  80.                             break;
  81.                         case 'c':
  82.                             filter_type = 1;
  83.                             break;
  84.                         case 'g':
  85.                             filter_type = 2;
  86.                             break;
  87.                         case 'm':
  88.                             filter_type = 3;
  89.                             break;
  90.                         case 'p':
  91.                             filter_type = 4;
  92.                             break;
  93.                         default:
  94.                             error = TRUE;
  95.                             break;
  96.                     }
  97.                     filter_width = atof(&argv[i][3]);
  98.                     break;
  99.                 case 'G':
  100.                     fout = &argv[i][2];
  101.                     break;
  102.                 case 'D':
  103.                     distance_flag = atoi(&argv[i][2]);
  104.                     break;
  105.                 case 'I':
  106.                     gmt_getinc (&argv[i][2], &dx_new, &dy_new);
  107.                     new_increment = TRUE;
  108.                     break;
  109.                 case 'N':    /* Pixel registration */
  110.                     one_or_zero = 0;
  111.                     break;
  112.                 default:
  113.                     error = TRUE;
  114.                     gmt_default_error (argv[i][1]);
  115.                     break;
  116.             }
  117.         }
  118.         else
  119.             fin = argv[i];
  120.     }
  121.     
  122.     if (argc == 1 || gmt_quick) {
  123.         fprintf (stderr, "grdfilter %s - Filter a 2-D netCDF grdfile in the Time domain\n\n", GMT_VERSION);
  124.         fprintf(stderr,"usage: grdfilter input_file -D<distance_flag> -F<type><filter_width>\n");
  125.         fprintf(stderr,"\t-G<output_file> [-I<xinc>[m|c][/<yinc>[m|c]] ] [-R<west/east/south/north>] [-V]\n");
  126.         
  127.         if (gmt_quick) exit (-1);
  128.         
  129.         fprintf(stderr,"\tDistance flag determines how grid (x,y) maps into distance units of filter width as follows:\n");
  130.         fprintf(stderr,"\t   -D0 grid x,y same units as <filter_width>, cartesian Distances.\n");
  131.         fprintf(stderr,"\t   -D1 grid x,y in degrees, <filter_width> in km, cartesian Distances.\n");
  132.         fprintf(stderr,"\t   -D2 grid x,y in degrees, <filter_width> in km, x_scaled by cos(middle y), cartesian Distances.\n");
  133.         fprintf(stderr,"\t   These first three options are faster; they allow weight matrix to be computed only once.\n");
  134.         fprintf(stderr,"\t   Next two options are slower; weights must be recomputed for each scan line.\n");
  135.         fprintf(stderr,"\t   -D3 grid x,y in degrees, <filter_width> in km, x_scale varies as cos(y), cartesian Distances.\n");
  136.         fprintf(stderr,"\t   -D4 grid x,y in degrees, <filter_width> in km, spherical Distances.\n");
  137.         fprintf(stderr,"\t-F sets the filter type and full (6 sigma) filter-width.  Choose between\n");
  138.         fprintf(stderr,"\t   (b)oxcar, (c)osine arch, (g)aussian, (m)edian filters\n");
  139.         fprintf(stderr,"\t   or p(maximum liklihood Probability estimator -- a mode estimator)\n");
  140.         fprintf(stderr,"\t-G sets output name for filtered grdfile\n");
  141.         fprintf(stderr, "\n\tOPTIONS:\n");
  142.         fprintf(stderr,"\t-I for new Increment of output grid; enter xinc, optionally xinc/yinc.\n");
  143.         fprintf(stderr,"\t   Default is yinc = xinc.  Append an m [or c] to xinc or yinc to indicate minutes [or seconds];\n");
  144.         fprintf(stderr,"\t   The new xinc and yinc should be divisible by the old ones (new lattice is subset of old).\n");
  145.         fprintf(stderr, "\t-N Force pixel registration for output grid [Default is gridline registration]\n");
  146.         fprintf(stderr, "\t-R for new Range of output grid; enter <WESN> (xmin, xmax, ymin, ymax) separated by slashes.\n");
  147.         explain_option ('V');
  148.         exit(-1);
  149.     }
  150.  
  151.     if (!fout) {
  152.         fprintf (stderr, "%s: GMT SYNTAX ERROR -G option:  Must specify output file\n", gmt_program);
  153.         error++;
  154.     }
  155.     if (!fin) {
  156.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", gmt_program);
  157.         error++;
  158.     }
  159.     if (distance_flag < 0 || distance_flag > 4) {
  160.         fprintf (stderr, "%s: GMT SYNTAX ERROR -D option:  Choose among 1, 2, 3, or 4\n", gmt_program);
  161.         error++;
  162.     }
  163.     if (filter_type < 0 || filter_width <= 0.0) {
  164.         fprintf (stderr, "%s: GMT SYNTAX ERROR -F option:  Correct syntax:\n", gmt_program);
  165.         fprintf (stderr, "\t-FX<width>, with X one of bcgmp, width is filter fullwidth\n");
  166.         error++;
  167.     }
  168.     if (error) exit (-1);
  169.     
  170.     if (west_new != east_new) new_range = TRUE;
  171.  
  172.     if (read_grd_info (fin, &h)) {
  173.         fprintf (stderr, "grdfilter: Error opening file %s\n", fin);
  174.         exit (-1);
  175.     }
  176.     grd_init (&h, argc, argv, TRUE);    /* Update command history only */
  177.  
  178.     /* Read the input grd file and close it  */
  179.     
  180.     input = (float *) memory (CNULL, (int)(h.nx * h.ny), sizeof(float), "grdfilter");
  181.  
  182.     if (read_grd (fin, &h, input, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  183.         fprintf (stderr, "grdfilter: Error reading file %s\n", fin);
  184.         exit (-1);
  185.     }
  186.     
  187.     last_median = 0.5 * (h.z_min + h.z_max);
  188.  
  189.     /* Check range of output area and set i,j offsets, etc.  */
  190.  
  191.     if (!new_range) {
  192.         west_new = h.x_min;
  193.         east_new = h.x_max;
  194.         south_new = h.y_min;
  195.         north_new = h.y_max;
  196.     }
  197.     if (!new_increment) {
  198.         dx_new = h.x_inc;
  199.         dy_new = h.y_inc;
  200.     }
  201.  
  202.     if (west_new < h.x_min) error = TRUE;
  203.     if (east_new > h.x_max) error = TRUE;
  204.     if (south_new < h.y_min) error = TRUE;
  205.     if (north_new > h.y_max) error = TRUE;
  206.     if (dx_new <= 0.0) error = TRUE;
  207.     if (dy_new <= 0.0) error = TRUE;
  208.     
  209.     if (error) {
  210.         fprintf(stderr,"grdfilter: New WESN incompatible with old.\n");
  211.         exit(-1);
  212.     }
  213.  
  214.     /* We can save time by computing a weight matrix once [or once pr scanline] only
  215.        if new grid spacing is multiple of old spacing */
  216.        
  217.     fast_way = ((rint (dx_new / h.x_inc) == (dx_new / h.x_inc)) && (rint (dy_new / h.y_inc) == (dy_new / h.y_inc)));
  218.     
  219.     if (!fast_way) {
  220.         fprintf (stderr, "grdfilter: Warning - Your output grid spacing is such that filter-weights must\n");
  221.         fprintf (stderr, "be recomputed for every output node, so expect this run to be slow.  Calculations\n");
  222.         fprintf (stderr, "can be speeded up significantly if output grid spacing is chosen to be a multiple\n");
  223.         fprintf (stderr, "of the input grid spacing.  If the odd output grid is necessary, consider using\n");
  224.         fprintf (stderr, "a \'fast\' grid for filtering and then resample onto your desired grid with grdsample.\n");
  225.     }
  226.  
  227.     nx_out = one_or_zero + rint ( (east_new - west_new) / dx_new);
  228.     ny_out = one_or_zero + rint ( (north_new - south_new) / dy_new);
  229.     
  230.     output = (float *) memory (CNULL, nx_out*ny_out, sizeof(float), "grdfilter");
  231.     i_origin = (int *) memory (CNULL, nx_out, sizeof(int), "grdfilter");
  232.     if (!fast_way) x_shift = (double *) memory (CNULL, nx_out, sizeof(double), "grdfilter");
  233.  
  234.     xincnew2 = (one_or_zero) ? 0.0 : 0.5 * dx_new;
  235.     yincnew2 = (one_or_zero) ? 0.0 : 0.5 * dy_new;
  236.     xincold2 = (h.node_offset) ? 0.0 : 0.5 * h.x_inc;
  237.     yincold2 = (h.node_offset) ? 0.0 : 0.5 * h.y_inc;
  238.     offset = (h.node_offset) ? 0.0 : 0.5;
  239.     
  240.     if (fast_way && h.node_offset == one_or_zero) {    /* multiple grid but one is pix, other is grid */
  241.         x_fix = 0.5 * h.x_inc;
  242.         y_fix = 0.5 * h.y_inc;
  243.         shift = (x_fix != 0.0 || y_fix != 0.0);
  244.     }
  245.     
  246.     /* Set up weight matrix and i,j range to test  */
  247.  
  248.     x_scale = y_scale = 1.0;
  249.     if (distance_flag > 0) {
  250.         x_scale *= deg2km;
  251.         y_scale *= deg2km;
  252.     }
  253.     if (distance_flag == 2)
  254.         x_scale *= cos(D2R * 0.5 * ( north_new + south_new) );
  255.     else if (distance_flag > 2) {
  256.         if ( fabs(south_new) > north_new )
  257.             x_scale *= cos(D2R * south_new);
  258.         else
  259.             x_scale *= cos(D2R * north_new);
  260.     }
  261.     x_width = filter_width / (h.x_inc * x_scale);
  262.     y_width = filter_width / (h.y_inc * y_scale);
  263.     y_half_width = ceil(y_width) / 2;
  264.     x_half_width = ceil(x_width) / 2;
  265.  
  266.     nx_fil = 2 * x_half_width + 1;
  267.     ny_fil = 2 * y_half_width + 1;
  268.     weight = (double *) memory (CNULL, nx_fil*ny_fil, sizeof(double), "grdfilter");
  269.  
  270.     slow = (filter_type == 3 || filter_type == 4);    /* Will require sorting etc */
  271.     
  272.     if (slow) work_array = (double *) memory (CNULL, nx_fil*ny_fil, sizeof(double), "grdfilter");
  273.  
  274.     if (gmtdefs.verbose) {
  275.         fprintf(stderr,"grdfilter: Input nx,ny = (%d %d), output nx,ny = (%d %d), filter nx,ny = (%d %d)\n",
  276.             h.nx, h.ny, nx_out, ny_out, nx_fil, ny_fil);
  277.         if (filter_type == 0) fprintf(stderr,"grdfilter: Filter type is Boxcar.\n");
  278.         if (filter_type == 1) fprintf(stderr,"grdfilter: Filter type is Cosine Arch.\n");
  279.         if (filter_type == 2) fprintf(stderr,"grdfilter: Filter type is Gaussian.\n");
  280.         if (filter_type == 3) fprintf(stderr,"grdfilter: Filter type is Median.\n");
  281.         if (filter_type == 4) fprintf(stderr,"grdfilter: Filter type is Mode.\n");
  282.     }
  283.     
  284.     /* Compute nearest xoutput i-indices and shifts once */
  285.     
  286.     for (i_out = 0; i_out < nx_out; i_out++) {
  287.         x_out = west_new + i_out * dx_new + xincnew2;
  288.         i_origin[i_out] = floor(((x_out - h.x_min) / h.x_inc) + offset);
  289.         if (!fast_way) x_shift[i_out] = x_out - (h.x_min + i_origin[i_out] * h.x_inc + xincold2);
  290.     }
  291.     
  292.     /* Now we can do the filtering  */
  293.  
  294.     /* Determine how much effort to compute weights:
  295.         1 = Compute weights once for entire grid
  296.         2 = Compute weights once per scanline
  297.         3 = Compute weights for every output point [slow]
  298.     */
  299.     
  300.     if (fast_way && distance_flag <= 2)
  301.         effort_level = 1;
  302.     else if (fast_way && distance_flag > 2)
  303.         effort_level = 2;
  304.     else
  305.         effort_level = 3;
  306.         
  307.     if (effort_level == 1) {    /* Only need this once */
  308.         y_out = north_new - yincnew2;
  309.         set_weight_matrix (nx_fil, ny_fil, y_out, north_new, south_new, h.x_inc, h.y_inc, filter_width, filter_type, distance_flag, x_fix, y_fix, shift);
  310.     }
  311.     
  312.     for (j_out = 0; j_out < ny_out; j_out++) {
  313.     
  314.         if (gmtdefs.verbose) fprintf (stderr, "grdfilter: Processing output line %d\r", j_out);
  315.         y_out = north_new - j_out * dy_new - yincnew2;
  316.         j_origin = floor(((h.y_max - y_out) / h.y_inc) + offset);
  317.         if (effort_level == 2)
  318.             set_weight_matrix (nx_fil, ny_fil, y_out, north_new, south_new, h.x_inc, h.y_inc, filter_width, filter_type, distance_flag, x_fix, y_fix, shift);
  319.         if (!fast_way) y_shift = y_out - (h.y_max - j_origin * h.y_inc - yincold2);
  320.         
  321.         for (i_out = 0; i_out < nx_out; i_out++) {
  322.         
  323.             if (effort_level == 3)
  324.                 set_weight_matrix (nx_fil, ny_fil, y_out, north_new, south_new, h.x_inc, h.y_inc, filter_width, filter_type, distance_flag, x_shift[i_out], y_shift, fast_way);
  325.             wt_sum = value = 0.0;
  326.             n_in_median = 0;
  327.             ij_out = j_out * nx_out + i_out;
  328.             
  329.             for (ii = -x_half_width; ii <= x_half_width; ii++) {
  330.                 i_in = i_origin[i_out] + ii;
  331.                 if ( (i_in < 0) || (i_in >= h.nx) ) continue;
  332.  
  333.                 for (jj = -y_half_width; jj <= y_half_width; jj++) {
  334.                     j_in = j_origin + jj;
  335.                     if ( (j_in < 0) || (j_in >= h.ny) ) continue;
  336.                                         
  337.                     ij_wt = (jj + y_half_width) * nx_fil + ii + x_half_width;
  338.                     if (weight[ij_wt] < 0.0) continue;
  339.                     
  340.                     ij_in = j_in*h.nx + i_in;
  341.                     if (bad_float ((double)input[ij_in])) continue;
  342.  
  343.                     /* Get here when point is usable  */
  344.                     if (slow) {
  345.                         work_array[n_in_median] = input[ij_in];
  346.                         n_in_median++;
  347.                     }
  348.                     else {
  349.                         value += input[ij_in] * weight[ij_wt];
  350.                         wt_sum += weight[ij_wt];
  351.                     }
  352.                 }
  353.             }
  354.             
  355.             /* Now we have done the convolution and we can get the value  */
  356.             
  357.             if (slow) {
  358.                 if (n_in_median) {
  359.                     if (filter_type == 3)
  360.                         median (work_array, n_in_median, h.z_min, h.z_max, last_median, &this_median);
  361.                     else
  362.                         mode (work_array, n_in_median, n_in_median/2, TRUE, &this_median);
  363.                     output[ij_out] = this_median;
  364.                     last_median = this_median;
  365.                 }
  366.                 else {
  367.                     output[ij_out] = gmt_NaN;
  368.                     n_nan++;
  369.                 }
  370.             }
  371.             else {
  372.                 if (wt_sum == 0.0) {    /* Assign value = gmt_NaN */
  373.                     n_nan++;
  374.                     output[ij_out] = gmt_NaN;
  375.                 }
  376.                 else
  377.                     output[ij_out] = value / wt_sum;
  378.             }
  379.         }
  380.     }
  381.     if (gmtdefs.verbose) fprintf (stderr, "\n");
  382.     
  383.     /* At last, that's it!  Output: */
  384.  
  385.     if (n_nan) fprintf (stderr, "grdfilter: Unable to estimate value at %d nodes, set to NaN\n", n_nan);
  386.     
  387.     h.nx = nx_out;
  388.     h.ny = ny_out;
  389.     h.x_min = west_new;
  390.     h.x_max = east_new;
  391.     h.y_min = south_new;
  392.     h.y_max = north_new;
  393.     h.x_inc = dx_new;
  394.     h.y_inc = dy_new;
  395.     h.node_offset = !one_or_zero;
  396.     
  397.     if (write_grd (fout, &h, output, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  398.         fprintf (stderr, "grdfilter: Error writing file %s\n", fout);
  399.         exit (-1);
  400.     }
  401.     
  402.     free ((char *) input);
  403.     free ((char *) output);
  404.     free ((char *) weight);
  405.     free ((char *) i_origin);
  406.     if (slow) free ((char *) work_array);
  407.     if (!fast_way) free ((char *) x_shift);
  408.     
  409.     gmt_end (argc, argv);
  410. }
  411.  
  412. int    set_weight_matrix (nx_f, ny_f, y_0, north, south, dx, dy, f_wid, f_flag, d_flag, x_off, y_off, fast)
  413.  
  414. double    y_0, north, south, dx, dy, f_wid, x_off, y_off;    /* Last two gives offset between output node and 'origin' input node for this window (0,0 for integral grids) */
  415. int    nx_f, ny_f;
  416. BOOLEAN fast;    /* TRUE when input/output grids are offset by integer values in dx/dy */
  417. int    f_flag, d_flag;
  418. {
  419.  
  420.     int    i, j, ij, i_half, j_half;
  421.     double    x_scl, y_scl, f_half, r_f_half, sigma, sig_2;
  422.     double    y1, y2, theta, x, y, r;
  423.  
  424.     /* Set Scales  */
  425.  
  426.     y_scl = (d_flag) ? deg2km : 1.0;
  427.     if (d_flag < 2)
  428.         x_scl = y_scl;
  429.     else if (d_flag == 2)
  430.         x_scl = deg2km * cos (0.5 * D2R * (north + south));
  431.     else
  432.         x_scl = deg2km * cos (D2R * y_0);
  433.  
  434.     /* Get radius, weight, etc.  */
  435.  
  436.     i_half = nx_f / 2;
  437.     j_half = ny_f / 2;
  438.     f_half = 0.5 * f_wid;
  439.     r_f_half = 1.0 / f_half;
  440.     sigma = f_wid / 6.0;
  441.     sig_2 = -0.5 / (sigma * sigma);
  442.     for (i = -i_half; i <= i_half; i++) {
  443.         for (j = -j_half; j <= j_half; j++) {
  444.             ij = (j + j_half) * nx_f + i + i_half;
  445.             if (fast && i == 0)
  446.                 r = (j == 0) ? 0.0 : j * y_scl * dy;
  447.             else if (fast && j == 0)
  448.                 r = i * x_scl * dx;
  449.             else if (d_flag < 4) {
  450.                 x = x_scl * (i * dx - x_off);
  451.                 y = y_scl * (j * dy - y_off);
  452.                 r = hypot (x, y);
  453.             }
  454.             else {
  455.                 theta = D2R * (i * dx - x_off);
  456.                 y1 = D2R * (90.0 - y_0);
  457.                 y2 = D2R * (90.0 - (y_0 + (j * dy - y_off)) );
  458.                 r = d_acos(cos(y1)*cos(y2) + sin(y1)*sin(y2)*cos(theta)) * deg2km * R2D;
  459.             }
  460.             /* Now we know r in f_wid units  */
  461.             
  462.             if (r > f_half) {
  463.                 weight[ij] = -1.0;
  464.                 continue;
  465.             }
  466.             else if (f_flag == 3) {
  467.                 weight[ij] = 1.0;
  468.                 continue;
  469.             }
  470.             else {
  471.                 if (f_flag == 0)
  472.                     weight[ij] = 1.0;
  473.                 else if (f_flag == 1)
  474.                     weight[ij] = 1.0 + cos(M_PI * r * r_f_half);
  475.                 else
  476.                     weight[ij] = exp(r * r * sig_2);
  477.             }
  478.         }
  479.     }
  480. }
  481.