-
Notifications
You must be signed in to change notification settings - Fork 377
/
Copy pathgrdfilter_mt.c
1289 lines (1154 loc) · 59.9 KB
/
grdfilter_mt.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*--------------------------------------------------------------------
*
* Copyright (c) 1991-2025 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
* See LICENSE.TXT file for copying and redistribution conditions.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation; version 3 or any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* Contact info: www.generic-mapping-tools.org
*--------------------------------------------------------------------*/
/*
* Brief synopsis: grdfilter.c reads a grid file and creates filtered grid file.
* User selects boxcar, gaussian, cosine arch, median, mode or spatial filters,
* and sets distance scaling as appropriate for map work, etc.
*
* Author: W.H.F. Smith
* Date: 1-JAN-2010
* Version: 6 API
*/
/*
This is an experimental multi-threaded version that uses gthreads from GLIB an hence depends on that lib
that on its turn depends on linintl (gettext)
To compile I patched src/CMakeList.txt by adding these two lines
set (HAVE_GLIB TRUE CACHE INTERNAL "System has GLIB")
include_directories (${GLIB_INCLUDE_DIR})
list (APPEND GMT_OPTIONAL_LIBRARIES ${GLIB_LIBRARY})
and added this to ConfigUserAdvanced.cmake
# Set location of GLIB ...:
set (GLIB_INCLUDE_DIR "C:/programs/compa_libs/glib-2.38.2/compileds/${VC}_${BITAGE}/include/glib-2.0")
set (GLIB_LIBRARY "C:/programs/compa_libs/glib-2.38.2/compileds/${VC}_${BITAGE}/lib/glib-2.0.lib")
maybe you don't need the above if cmake is able to find glib in your system.
Use undocumented (and temporary) option -z to set the number of threads. e.g. -z2, -z4, ... or -za to use all available
*/
#include "gmt_dev.h"
#define THIS_MODULE_CLASSIC_NAME "grdfilter"
#define THIS_MODULE_MODERN_NAME "grdfilter"
#define THIS_MODULE_LIB "core"
#define THIS_MODULE_PURPOSE "Filter a grid in the space (or time) domain"
#define THIS_MODULE_KEYS "<G{,FG(=1,GG}"
#define THIS_MODULE_NEEDS ""
#define THIS_MODULE_OPTIONS "-RVf"
struct GRDFILTER_CTRL {
struct GRDFILT_In {
bool active;
char *file;
} In;
struct GRDFILT_A { /* -A<a|r|w|c>row/col */
bool active;
char mode;
unsigned int ROW, COL;
double x, y;
char *file;
} A;
struct GRDFILT_D { /* -D<distflag> */
bool active;
int mode; /* -1 to 5 */
} D;
struct GRDFILT_F { /* <type>[-]<filter_width>[/<width2>][<mode>] */
bool active;
bool highpass;
bool custom;
bool operator;
char filter; /* Character codes for the filter */
char *file;
double width, width2, quantile;
bool rect;
int mode; /*-1 0 +1 */
} F;
struct GRDFILT_G { /* -G<file> */
bool active;
char *file;
} G;
struct GRDFILT_N { /* -Np|i|r */
bool active;
unsigned int mode; /* 0 is default (i), 1 is replace (r), 2 is preserve (p) */
} N;
struct GRDFILT_T { /* -T */
bool active;
} T;
struct GRDFILT_z { /* -z */
bool active;
int n_threads;
} z;
};
/* Forward/Inverse Spherical Mercator coordinate transforms */
#define IMG2LAT(img) (2.0*atand(exp((img)*D2R))-90.0)
#define LAT2IMG(lat) (R2D*log(tand(0.5*((lat)+90.0))))
/* Local ij calculation for weight matrix */
#define WT_IJ(F,jj,ii) ((jj) + F.y_half_width) * F.n_columns + (ii) + F.x_half_width
#define GRDFILTER_N_PARS 7
#define GRDFILTER_WIDTH 0
#define GRDFILTER_HALF_WIDTH 1
#define GRDFILTER_X_SCALE 2
#define GRDFILTER_Y_SCALE 3
#define GRDFILTER_INV_R_SCALE 4
#define GRDFILTER_X_DIST 5
#define GRDFILTER_Y_DIST 6
#define GRDFILTER_FILTERS "bcgfompLlUu"
#define GRDFILTER_MODES "012345p"
/* Distance modes: */
#define GRDFILTER_NOTSET -2
#define GRDFILTER_XY_PIXEL -1
#define GRDFILTER_XY_CARTESIAN 0
#define GRDFILTER_GEO_CARTESIAN 1
#define GRDFILTER_GEO_FLATEARTH1 2
#define GRDFILTER_GEO_FLATEARTH2 3
#define GRDFILTER_GEO_SPHERICAL 4
#define GRDFILTER_GEO_MERCATOR 5
/* Convolution filters: */
#define GRDFILTER_BOXCAR 0
#define GRDFILTER_COSINE 1
#define GRDFILTER_GAUSSIAN 2
#define GRDFILTER_CUSTOM 3
#define GRDFILTER_OPERATOR 4
/* Spatial non-convolution filters: */
#define GRDFILTER_MEDIAN 5
#define GRDFILTER_MODE 6
#define GRDFILTER_MIN 7
#define GRDFILTER_MINPOS 8
#define GRDFILTER_MAX 9
#define GRDFILTER_MAXNEG 10
#define GRDFILTER_MEDIAN_SPH 11
#define GRDFILTER_MODE_SPH 12
#define GRDFILTER_N_FILTERS 11 /* Not counting the two _SPH variants */
#define NAN_IGNORE 0
#define NAN_REPLACE 1
#define NAN_PRESERVE 2
struct FILTER_INFO {
unsigned int n_columns; /* The max number of filter weights in x-direction */
unsigned int n_rows; /* The max number of filter weights in y-direction */
int x_half_width; /* Number of filter nodes to either side needed at this latitude */
int y_half_width; /* Number of filter nodes above/below this point (ny_f/2) */
unsigned int d_flag;
bool rect; /* For 2-D rectangular filtering */
bool debug; /* Normally unused except under DEBUG */
double dx, dy; /* Grid spacing in original units */
double y_min, y_max; /* Grid limits in y(lat) */
double *x, *y; /* Distances in original units along x and y to distance nodes */
double par[5]; /* [0] is filter width, [1] is 0.5*filter_width, [2] is xscale, [3] is yscale, [4] is 1/r_half for filter */
double x_off, y_off; /* Offsets relative to original grid */
char *visit; /* true/false array to keep track of which longitude nodes we have already visited once */
double (*weight_func) (double, double *);
double (*radius_func) (struct GMT_CTRL *, double, double, double, double, double *);
};
struct THREAD_STRUCT {
bool fast_way, spherical, slow, slower;
int *col_origin, nx_wrap;
unsigned int row, r_start, r_stop, effort_level, filter_type, thread_num;
double x_fix, y_fix, last_median;
double *weight, *x_shift, *par;
struct GMT_CTRL *GMT;
struct GRDFILTER_CTRL *Ctrl;
struct GMT_GRID *Gin;
struct GMT_GRID *Gout;
struct GMT_GRID *A;
struct FILTER_INFO F;
};
static void *threading_function (void *args);
void grdfilter_threaded_function (struct THREAD_STRUCT *t);
static void *New_grdfilter_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
struct GRDFILTER_CTRL *C;
C = gmt_M_memory (GMT, NULL, 1, struct GRDFILTER_CTRL);
/* Initialize values whose defaults are not 0/false/NULL */
C->D.mode = GRDFILTER_NOTSET;
C->F.quantile = 0.5; /* Default is median */
C->z.n_threads = 1; /* Default number of threads (for the case of no multi-threading) */
return (C);
}
static void Free_grdfilter_Ctrl (struct GMT_CTRL *GMT, struct GRDFILTER_CTRL *C) { /* Deallocate control structure */
if (!C) return;
gmt_M_str_free (C->In.file);
gmt_M_str_free (C->F.file);
gmt_M_str_free (C->G.file);
gmt_M_str_free (C->A.file);
gmt_M_free (GMT, C);
}
/* -----------------------------------------------------------------------------------*/
static void *grdfilter_thread_function (void *args) {
grdfilter_threaded_function ((struct THREAD_STRUCT *)args);
return NULL;
}
void set_weight_matrix (struct GMT_CTRL *GMT, struct FILTER_INFO *F, double *weight, double output_lat, double par[], double x_off, double y_off)
{
/* x_off and y_off give offset between output node and 'origin' input node for this window (0,0 for integral grids).
* fast is true when input/output grids are offset by integer values in dx/dy.
* Here, par[0] = filter_width, par[1] = filter_width / 2, par[2] = x_scale, part[3] = y_scale, and
* par[4] is the normalization distance needed for the Cosine-bell or Gaussian weight function.
*/
int i, j;
int64_t ij;
double x, y, yc, y0, r, ry = 0.0, inv_x_half_width = 0.0, inv_y_half_width = 0.0;
yc = y0 = output_lat - y_off; /* Input latitude of central point (i,j) = (0,0) */
if (F->d_flag == GRDFILTER_GEO_MERCATOR) yc = IMG2LAT (yc); /* Recover actual latitude in IMG grid at this center point */
if (F->rect) {
inv_x_half_width = 1.0 / F->x_half_width;
inv_y_half_width = 1.0 / F->y_half_width;
}
for (j = -F->y_half_width; j <= F->y_half_width; j++) {
y = y0 + ((j < 0) ? F->y[-j] : -F->y[j]); /* y or latitude at this row */
if (F->d_flag > GRDFILTER_GEO_FLATEARTH1 && (y < F->y_min || y > F->y_max)) { /* This filter row is outside input grid domain */
for (i = -F->x_half_width, ij = (j + F->y_half_width) * F->n_columns; i <= F->x_half_width; i++, ij++) weight[ij] = -1.0;
continue; /* Done with this row */
}
if (F->d_flag == GRDFILTER_GEO_MERCATOR) y = IMG2LAT (y); /* Recover actual latitudes */
if (F->rect) ry = inv_y_half_width * j; /* -1 to +1 */
for (i = -F->x_half_width; i <= F->x_half_width; i++) {
x = (i < 0) ? -F->x[-i] : F->x[i];
ij = (j + F->y_half_width) * F->n_columns + i + F->x_half_width;
assert (ij >= 0);
if (F->rect) { /* 2-D rectangular filtering; radius not used as we use x/x_half_width and ry instead */
weight[ij] = F->weight_func (inv_x_half_width * i, par) * F->weight_func (ry, par);
}
else {
r = F->radius_func (GMT, x_off, yc, x, y, par);
weight[ij] = (r > par[GRDFILTER_HALF_WIDTH]) ? -1.0 : F->weight_func (r, par);
#ifdef DEBUG
if (F->debug) weight[ij] = (r > par[GRDFILTER_HALF_WIDTH]) ? -1.0 : r;
#endif
}
}
}
}
/* Various functions that will be accessed via pointers */
double CartRadius (struct GMT_CTRL *GMT, double x0, double y0, double x1, double y1, double par[])
{ /* Plain Cartesian distance (par is not used) */
gmt_M_unused (GMT);
gmt_M_unused (par);
return (hypot (x0 - x1, y0 - y1));
}
double CartScaledRadius (struct GMT_CTRL *GMT, double x0, double y0, double x1, double y1, double par[])
{ /* Plain scaled Cartesian distance (xscale = yscale) */
gmt_M_unused (GMT);
return (par[GRDFILTER_X_SCALE] * hypot (x0 - x1, y0 - y1));
}
double CartScaledRect (struct GMT_CTRL *GMT, double x0, double y0, double x1, double y1, double par[])
{ /* Pass dx,dy via par[GRDFILTER_X|Y_DIST] and return a r that is either in or out */
double r;
gmt_M_unused (GMT);
par[GRDFILTER_X_DIST] = par[GRDFILTER_X_SCALE] * (x0 - x1);
par[GRDFILTER_Y_DIST] = par[GRDFILTER_Y_SCALE] * (y0 - y1);
r = (fabs (par[GRDFILTER_X_DIST]) > par[GRDFILTER_HALF_WIDTH] || fabs (par[GRDFILTER_Y_DIST]) > par[GRDFILTER_HALF_WIDTH]) ? 2.0 : 0.0;
return (r);
}
double FlatEarthRadius (struct GMT_CTRL *GMT, double x0, double y0, double x1, double y1, double par[])
{ /* Cartesian radius with different scales */
gmt_M_unused (GMT);
return (hypot (par[GRDFILTER_X_SCALE] * (x0 - x1), par[GRDFILTER_Y_SCALE] * (y0 - y1)));
}
double SphericalRadius (struct GMT_CTRL *GMT, double x0, double y0, double x1, double y1, double par[])
{ /* Great circle distance in km with polar wrap-around test on 2nd point */
gmt_M_unused (par);
if (fabs (y1) > 90.0) { /* Must find the point across the pole */
y1 = copysign (180.0 - fabs (y1), y1);
x1 += 180.0;
}
return (0.001 * gmt_great_circle_dist_meter (GMT, x0, y0, x1, y1));
}
double UnitWeight (double r, double par[])
{ /* Return unit weight since we know r is inside radius (par is not used) */
gmt_M_unused (r);
gmt_M_unused (par);
return (1.0);
}
double CosBellWeight (double r, double par[])
{ /* Return the cosine-bell filter weight for given r.
* The parameter r_f_half is the 5th parameter passed via par.
*/
return (1.0 + cos (M_PI * r * par[GRDFILTER_INV_R_SCALE]));
}
double GaussianWeight (double r, double par[])
{ /* Return the Gaussian filter weight for given r.
* The parameter sig_2 is the 5th parameter passed via par.
*/
return (exp (r * r * par[GRDFILTER_INV_R_SCALE]));
}
struct GMT_GRID * init_area_weights (struct GMT_CTRL *GMT, struct GMT_GRID *G, int mode, char *file)
{
/* Precalculate the area weight of each node. There are several considerations:
* 1. Mercator img grids (d_flag == GRDFILTER_GEO_MERCATOR) has irregular latitude spacing so we
* must compute the north and south latitude bounds for each cell to get area.
* 2. Geographic poles for grid-registered grids require a separate weight formula
* 3. Grid-registered grids have boundary nodes that only apply to 1/2 the area
* (and the four corners (unless poles) only 1/4 the area of other cells).
*/
unsigned int row, col;
uint64_t ij;
double row_weight, col_weight, dy_half = 0.0, dx, y, lat, lat_s, lat_n, s2 = 0.0;
struct GMT_GRID *A = NULL;
/* Base the area weight grid on the input grid domain and increments. */
if ((A = GMT_Create_Data (GMT->parent, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, G->header->wesn, G->header->inc, \
G->header->registration, GMT_NOTSET, NULL)) == NULL) return (NULL);
if (mode > GRDFILTER_XY_CARTESIAN) { /* Geographic data */
if (mode == GRDFILTER_GEO_MERCATOR) dy_half = 0.5 * A->header->inc[GMT_Y]; /* Half img y-spacing */
dx = A->header->inc[GMT_X] * R2D; /* Longitude increment in radians */
s2 = sind (0.5 * A->header->inc[GMT_Y]); /* Holds sin (del_y/2) */
}
else /* Cartesian */
dx = A->header->inc[GMT_X];
gmt_M_row_loop (GMT, A, row) { /* Loop over the rows */
if (mode == GRDFILTER_GEO_MERCATOR) { /* Adjust lat if IMG grid. Note: these grids can never reach a pole. */
y = gmt_M_grd_row_to_y (GMT, row, A->header); /* Current input Merc y */
lat = IMG2LAT (y); /* Get actual latitude */
lat_s = IMG2LAT (y - dy_half); /* Bottom grid cell latitude */
lat_n = IMG2LAT (y + dy_half); /* Top grid cell latitude */
row_weight = sind (lat_n) - sind (lat_s);
}
else if (mode > GRDFILTER_XY_CARTESIAN) { /* Geographic data, and watch for poles */
lat = gmt_M_grd_row_to_y (GMT, row, A->header); /* Current input latitude */
if (gmt_M_is_pole (lat)) /* Poles are different */
row_weight = 1.0 - cosd (0.5 * A->header->inc[GMT_Y]);
else { /* All other points */
row_weight = 2.0 * cosd (lat) * s2;
/* Note: No need for special weight-sharing between w/e gridline-reg grids since we explicitly only use the west node */
}
}
else { /* If Cartesian then row_weight is a constant 1 except for gridline-registered grids at top or bottom row */
row_weight = (A->header->registration == GMT_GRID_NODE_REG && (row == 0 || row == (A->header->n_rows-1))) ? 0.5 : 1.0; /* Share weight with repeat point */
row_weight *= A->header->inc[GMT_Y];
}
gmt_M_col_loop (GMT, A, row, col, ij) { /* Now loop over the columns */
col_weight = dx * ((A->header->registration == GMT_GRID_NODE_REG && (col == 0 || col == (A->header->n_columns-1))) ? 0.5 : 1.0);
A->data[ij] = (gmt_grdfloat)(row_weight * col_weight);
}
}
if (file) { /* For debug purposes: Save the area weight grid */
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Write area weight grid to file %s\n", file);
if (GMT_Set_Comment (GMT->parent, GMT_IS_GRID, GMT_COMMENT_IS_REMARK, "Area weight grid for debugging purposes", A)) return (NULL);
if (GMT_Write_Data (GMT->parent, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, file, A) != GMT_NOERROR) return (NULL);
}
return (A);
}
int GMT_grdfilter_usage (struct GMTAPI_CTRL *API, int level)
{
const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
GMT_Message (API, GMT_TIME_NONE, "usage: %s <ingrid> -D<distance_flag> -F<type>[-]<filter_width>[/<width2>][<mode>] -G<outgrid>\n", name);
GMT_Message (API, GMT_TIME_NONE, "\t[%s] [-Ni|p|r] [%s]\n\t[-T] [%s] [%s] [%s]\n", GMT_I_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_f_OPT, GMT_PAR_OPT);
if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
GMT_Message (API, GMT_TIME_NONE, "\t<ingrid> is the input grid file to be filtered.\n");
GMT_Message (API, GMT_TIME_NONE, "\t-D Distance flag determines how grid (x,y) maps into distance units of filter width as follows:\n");
GMT_Message (API, GMT_TIME_NONE, "\t Cartesian (x, y) Data:\n");
GMT_Message (API, GMT_TIME_NONE, "\t -Dp grid x,y with <filter_width> in pixels (must be an odd number), Cartesian distances.\n");
GMT_Message (API, GMT_TIME_NONE, "\t -D0 grid x,y same units as <filter_width>, Cartesian distances.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Geographic (lon,lat) Data:\n");
GMT_Message (API, GMT_TIME_NONE, "\t -D1 grid x,y in degrees, <filter_width> in km, Cartesian distances.\n");
GMT_Message (API, GMT_TIME_NONE, "\t -D2 grid x,y in degrees, <filter_width> in km, x_scaled by cos(middle y), Cartesian distances.\n");
GMT_Message (API, GMT_TIME_NONE, "\t The options above are faster; they allow weight matrix to be computed only once.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Next three options are slower; weights must be recomputed for each scan line.\n");
GMT_Message (API, GMT_TIME_NONE, "\t -D3 grid x,y in degrees, <filter_width> in km, x_scale varies as cos(y), Cartesian distances.\n");
GMT_Message (API, GMT_TIME_NONE, "\t -D4 grid x,y in degrees, <filter_width> in km, spherical distances.\n");
GMT_Message (API, GMT_TIME_NONE, "\t -D5 grid x,y in Mercator units (-Jm1), <filter_width> in km, spherical distances.\n");
GMT_Message (API, GMT_TIME_NONE, "\t-F Set the low-pass filter type and full diameter (6 sigma) filter-width.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Choose between convolution-type filters which differ in how weights are assigned\n");
GMT_Message (API, GMT_TIME_NONE, "\t and geospatial filters that seek to return a representative value.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Give a negative filter width to select highpass filtering [lowpass].\n");
GMT_Message (API, GMT_TIME_NONE, "\t Filters are isotropic. For rectangular filtering append /<width2> (requires -Dp|0).\n");
GMT_Message (API, GMT_TIME_NONE, "\t Convolution filters:\n");
GMT_Message (API, GMT_TIME_NONE, "\t b: Boxcar : a simple averaging of all points inside filter domain.\n");
GMT_Message (API, GMT_TIME_NONE, "\t c: Cosine arch : a weighted averaging with cosine arc weights.\n");
GMT_Message (API, GMT_TIME_NONE, "\t g: Gaussian : weighted averaging with Gaussian weights.\n");
GMT_Message (API, GMT_TIME_NONE, "\t f: Custom : Give grid file with custom filter weights (requires -D0).\n");
GMT_Message (API, GMT_TIME_NONE, "\t o: Operator : Give grid file with custom operator weights (requires -D0).\n");
GMT_Message (API, GMT_TIME_NONE, "\t Geospatial filters:\n");
GMT_Message (API, GMT_TIME_NONE, "\t l: Lower : return minimum of all points.\n");
GMT_Message (API, GMT_TIME_NONE, "\t L: Lower+ : return minimum of all +ve points.\n");
GMT_Message (API, GMT_TIME_NONE, "\t m: Median : return the median (50%% quantile) value of all points.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Append +q<quantile> to select another quantile (in 0-1 range) [0.5].\n");
GMT_Message (API, GMT_TIME_NONE, "\t p: Maximum likelihood probability estimator : return mode of all points.\n");
GMT_Message (API, GMT_TIME_NONE, "\t By default, we return the average if more than one mode is found.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Append - or + to the width to instead return the smallest or largest mode.\n");
GMT_Message (API, GMT_TIME_NONE, "\t u: Upper : return maximum of all points.\n");
GMT_Message (API, GMT_TIME_NONE, "\t U: Upper- : return maximum of all -ve points.\n");
GMT_Message (API, GMT_TIME_NONE, "\t-G Set output filename for filtered grid.\n");
GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
#ifdef DEBUG
GMT_Message (API, GMT_TIME_NONE, "\t-A DEBUG: Use -A<mode><lon/<lat> to instead save filter specifics at that point. Choose:\n");
GMT_Message (API, GMT_TIME_NONE, "\t mode as a (area weights), c (composite weight), r (radii), or w (filter weight).\n");
GMT_Message (API, GMT_TIME_NONE, "\t-A DEBUG: Save area weights to <file> with -As<file> [no save]\n");
#endif
GMT_Option (API, "I");
GMT_Message (API, GMT_TIME_NONE, "\t The new xinc and yinc should be divisible by the old ones (new lattice is subset of old).\n");
GMT_Message (API, GMT_TIME_NONE, "\t-N Specify how NaNs in the input grid should be treated. There are three options:\n");
GMT_Message (API, GMT_TIME_NONE, "\t -Ni skips all NaN values and returns a filtered value unless all are NaN [Default].\n");
GMT_Message (API, GMT_TIME_NONE, "\t -Np sets filtered output to NaN is any NaNs are found inside filter circle.\n");
GMT_Message (API, GMT_TIME_NONE, "\t -Nr sets filtered output to NaN if the corresponding input node was NaN.\n");
GMT_Message (API, GMT_TIME_NONE, "\t (only possible if the input and output grids are coregistered).\n");
GMT_Message (API, GMT_TIME_NONE, "\t-T Toggle between grid and pixel registration for output grid [Default is same as input registration].\n");
GMT_Message (API, GMT_TIME_NONE, "\t-R For new Range of output grid; enter <WESN> (xmin, xmax, ymin, ymax) separated by slashes.\n");
GMT_Message (API, GMT_TIME_NONE, "\t [Default uses the same region as the input grid].\n");
GMT_Option (API, "V,f,.");
return (GMT_MODULE_USAGE);
}
int GMT_grdfilter_parse (struct GMT_CTRL *GMT, struct GRDFILTER_CTRL *Ctrl, struct GMT_OPTION *options)
{
/* This parses the options provided to grdfilter and sets parameters in Ctrl.
* Note Ctrl has already been initialized and non-zero default values set.
* Any GMT common options will override values set previously by other commands.
* It also replaces any file names specified as input or output with the data ID
* returned when registering these sources/destinations with the API.
*/
unsigned int n_errors = 0, n_files = 0;
char c, a[GMT_LEN64] = {""}, b[GMT_LEN64] = {""}, txt[GMT_LEN256] = {""}, *p = NULL;
struct GMT_OPTION *opt = NULL;
struct GMTAPI_CTRL *API = GMT->parent;
for (opt = options; opt; opt = opt->next) { /* Process all the options given */
switch (opt->option) {
case '<': /* Input file (only one is accepted) */
if (n_files++ > 0) {n_errors++; continue; }
Ctrl->In.active = true;
if (opt->arg[0]) Ctrl->In.file = strdup (opt->arg);
if (GMT_Get_FilePath (GMT->parent, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file))) n_errors++;
break;
/* Processes program-specific parameters */
#ifdef DEBUG
case 'A': /* Debug mode options */
Ctrl->A.active = true;
if (opt->arg[0] == 's') {
Ctrl->A.file = strdup (&opt->arg[1]);
}
else {
Ctrl->A.mode = opt->arg[0];
sscanf (&opt->arg[1], "%[^/]/%s", a, b);
Ctrl->A.x = atof (a);
Ctrl->A.y = atof (b);
}
break;
#endif
case 'D': /* Distance mode */
Ctrl->D.active = true;
if (strchr (GRDFILTER_MODES, opt->arg[0])) { /* OK filter code */
Ctrl->D.mode = (opt->arg[0] == 'p') ? GRDFILTER_XY_PIXEL : atoi (opt->arg);
}
else {
GMT_Report (API, GMT_MSG_ERROR, "Expected -D<mode>, where mode is one of %s\n", GRDFILTER_MODES);
n_errors++;
}
break;
case 'F': /* Filter */
if (strchr (GRDFILTER_FILTERS, opt->arg[0])) { /* OK filter code */
Ctrl->F.active = true;
Ctrl->F.filter = opt->arg[0];
strncpy (txt, opt->arg, GMT_LEN256); /* Work on a copy */
if (Ctrl->F.filter == 'm') {
if ((p = strchr (txt, 'q'))) { /* Requested another quantile */
*(--p) = 0; /* Chop off the +q modifier */
Ctrl->F.quantile = atof (p+2);
}
}
if (Ctrl->F.filter == 'f' || Ctrl->F.filter == 'o') {
if (gmt_check_filearg (GMT, 'F', &opt->arg[1], GMT_IN, GMT_IS_GRID))
Ctrl->F.file = strdup (&opt->arg[1]);
else {
GMT_Report (API, GMT_MSG_ERROR, "Option -F%c: Cannot access filter weight grid %s\n", Ctrl->F.filter, &opt->arg[1]);
n_errors++;
}
Ctrl->F.width = 1.0; /* To avoid error checking below */
Ctrl->F.custom = true;
Ctrl->F.operator = (Ctrl->F.filter == 'o'); /* Means weightsum is zero so no normalization, please */
}
else if (strchr (txt, '/')) { /* Gave xwidth/ywidth for rectangular Cartesian filtering */
sscanf (&txt[1], "%[^/]/%s", a, b);
Ctrl->F.width = atof (a);
Ctrl->F.width2 = atof (b);
Ctrl->F.rect = true;
}
else
Ctrl->F.width = atof (&txt[1]);
if (Ctrl->F.width < 0.0) Ctrl->F.highpass = true;
Ctrl->F.width = fabs (Ctrl->F.width);
if (Ctrl->F.filter == 'p') { /* Check for some further info in case of mode filtering */
c = opt->arg[strlen(txt)-1];
if (c == '-') Ctrl->F.mode = -1;
if (c == '+') Ctrl->F.mode = +1;
}
}
else {
GMT_Report (API, GMT_MSG_ERROR, "Expected -Fx<width>, where x is one of %s\n", GRDFILTER_FILTERS);
n_errors++;
}
break;
case 'G': /* Output file */
Ctrl->G.active = true;
if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
if (GMT_Get_FilePath (GMT->parent, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
break;
case 'I': /* New grid spacings */
n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg);
break;
case 'N': /* Treatment of NaNs */
Ctrl->N.active = true;
switch (opt->arg[0]) {
case 'i':
Ctrl->N.mode = NAN_IGNORE; /* Default */
break;
case 'r':
Ctrl->N.mode = NAN_REPLACE; /* Replace */
break;
case 'p':
Ctrl->N.mode = NAN_PRESERVE; /* Preserve */
break;
default:
GMT_Report (API, GMT_MSG_ERROR, "Expected -Ni|p|r\n");
n_errors++;
break;
}
break;
case 'T': /* Toggle registration */
Ctrl->T.active = true;
break;
#ifdef HAVE_GLIB_GTHREAD
case 'z':
Ctrl->z.active = true;
if (opt->arg && opt->arg[0] == 'a') /* Use all processors available */
Ctrl->z.n_threads = g_get_num_processors();
else if (opt->arg)
Ctrl->z.n_threads = atoi(opt->arg);
if (Ctrl->z.n_threads == 0)
Ctrl->z.n_threads = 1;
else if (Ctrl->z.n_threads < 0)
Ctrl->z.n_threads = MAX(g_get_num_processors() - 1, 1); /* Max -1 but at least one */
else
Ctrl->z.n_threads = MIN((int)g_get_num_processors(), Ctrl->z.n_threads); /* No more than maximum available */
break;
#endif
default: /* Report bad options */
n_errors += gmt_default_option_error (GMT, opt);
break;
}
}
n_errors += gmt_M_check_condition (GMT, !Ctrl->G.active, "Option -G: Must specify output file\n");
n_errors += gmt_M_check_condition (GMT, !Ctrl->In.file, "Must specify input file\n");
n_errors += gmt_M_check_condition (GMT, !Ctrl->D.active, "Option -D: Choose from p or 0-5\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && (Ctrl->D.mode < GRDFILTER_XY_PIXEL || Ctrl->D.mode > GRDFILTER_GEO_MERCATOR),
"Option -D: Choose from p or 0-5\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->D.mode > GRDFILTER_XY_CARTESIAN && Ctrl->F.rect,
"Option -F: Rectangular Cartesian filtering requires -Dp|0\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->D.mode > GRDFILTER_XY_CARTESIAN && Ctrl->F.custom,
"Option -Ff|o: Custom Cartesian convolution requires -D0\n");
n_errors += gmt_M_check_condition (GMT, !Ctrl->F.active, "Option -F is required:\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->F.quantile < 0.0 || Ctrl->F.quantile > 1.0 ,
"Option -F: The quantile must be in the 0-1 range.\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->F.active && Ctrl->F.width == 0.0,
"Option -F: filter fullwidth must be nonzero.\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->F.rect && Ctrl->F.width2 <= 0.0,
"Option -F: Rectangular y-width filter must be nonzero.\n");
n_errors += gmt_M_check_condition (GMT, GMT->common.R.active[ISET] && (GMT->common.R.inc[GMT_X] <= 0.0 || GMT->common.R.inc[GMT_Y] <= 0.0),
"Option -I: Must specify positive increment(s)\n");
n_errors += gmt_M_check_condition (GMT, GMT->common.R.active[RSET] && GMT->common.R.active[ISET] && Ctrl->F.highpass,
"Option -F: Highpass filtering requires original -R -I\n");
return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
}
#define bailout(code) {gmt_M_free_options (mode); return (code);}
#define Return(code) {Free_grdfilter_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
int GMT_grdfilter_mt (void *V_API, int mode, void *args)
{
bool fast_way, slow = false, slower = false, same_grid = false;
bool spherical = false, full_360, visit_check = false, get_weight_sum = true;
unsigned int n_nan = 0, col_out, row_out, effort_level;
unsigned int filter_type, one_or_zero = 1, GMT_n_multiples = 0;
int i, tid = 0, col_in, row_in, ii, jj, *col_origin = NULL, nx_wrap = 0, error = 0;
#ifdef DEBUG
unsigned int n_conv = 0;
#endif
uint64_t ij_in, ij_out, ij_wt;
double x_scale = 1.0, y_scale = 1.0, x_width, y_width, par[GRDFILTER_N_PARS];
double x_out, wt_sum, last_median = 0.0;
double x_fix = 0.0, y_fix = 0.0, max_lat;
double merc_range, *weight = NULL, *x_shift = NULL;
double wesn[4], inc[2];
char filter_code[GRDFILTER_N_FILTERS] = {'b', 'c', 'g', 'f', 'o', 'm', 'p', 'l', 'L', 'u', 'U'};
char *filter_name[GRDFILTER_N_FILTERS+2] = {"Boxcar", "Cosine Arch", "Gaussian", "Custom", "Operator", "Median", "Mode", "Lower", \
"Lower+", "Upper", "Upper-", "Spherical Median", "Spherical Mode"};
struct GMT_GRID *Gin = NULL, *Gout = NULL, *Fin = NULL, *A = NULL;
struct FILTER_INFO F;
struct GRDFILTER_CTRL *Ctrl = NULL;
struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
struct GMT_OPTION *options = NULL;
struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
struct THREAD_STRUCT *threadArg;
#ifdef HAVE_GLIB_GTHREAD
GThread **threads;
#endif
/*----------------------- Standard module initialization and parsing ----------------------*/
if (API == NULL) return (GMT_NOT_A_SESSION);
if (mode == GMT_MODULE_PURPOSE) return (GMT_grdfilter_usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
if (!options || options->option == GMT_OPT_USAGE) bailout (GMT_grdfilter_usage (API, GMT_USAGE));/* Return the usage message */
if (options->option == GMT_OPT_SYNOPSIS) bailout (GMT_grdfilter_usage (API, GMT_SYNOPSIS)); /* Return the synopsis */
/* Parse the command-line arguments */
if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, module_kw, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
Ctrl = New_grdfilter_Ctrl (GMT); /* Allocate and initialize a new control structure */
if ((error = GMT_grdfilter_parse (GMT, Ctrl, options)) != 0) Return (error);
/*---------------------------- This is the grdfilter main code ----------------------------*/
GMT_Report (API, GMT_MSG_INFORMATION, "Processing input grid\n");
if ((Gin = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->In.file, NULL)) == NULL) { /* Get entire grid */
Return (API->error);
}
if (Ctrl->T.active) /* Make output grid of the opposite registration */
one_or_zero = (Gin->header->registration == GMT_GRID_PIXEL_REG) ? 1 : 0;
else
one_or_zero = (Gin->header->registration == GMT_GRID_PIXEL_REG) ? 0 : 1;
full_360 = (Ctrl->D.mode > GRDFILTER_XY_CARTESIAN && gmt_grd_is_global (GMT, Gin->header)); /* Periodic geographic grid */
if (Ctrl->D.mode == GRDFILTER_XY_PIXEL) { /* Special case where widths are given in pixels */
if (!doubleAlmostEqual (fmod (Ctrl->F.width, 2.0), 1.0)) {
GMT_Report (API, GMT_MSG_ERROR, "Option -Dp requires filter width given as an odd number of pixels\n");
Return (GMT_RUNTIME_ERROR);
}
Ctrl->F.width *= Gin->header->inc[GMT_X]; /* Scale up to give width */
if (Ctrl->F.rect) {
if (!doubleAlmostEqual (fmod (Ctrl->F.width2, 2.0), 1.0)) {
GMT_Report (API, GMT_MSG_ERROR, "Option -Dp requires filter y-width given as an odd number of pixels\n");
Return (GMT_RUNTIME_ERROR);
}
Ctrl->F.width2 *= Gin->header->inc[GMT_X]; /* Rectangular rather than isotropic Cartesian filtering */
}
Ctrl->D.mode = GRDFILTER_XY_CARTESIAN; /* From now on they are the same distances */
}
/* Check range of output area and set i,j offsets, etc. */
/* Use the -R region for output (if set); otherwise match input grid domain */
gmt_M_memcpy (wesn, (GMT->common.R.active[RSET] ? GMT->common.R.wesn : Gin->header->wesn), 4, double);
/* Use the -I increments for output (if set); otherwise match input grid increments */
gmt_M_memcpy (inc, (GMT->common.R.active[ISET] ? GMT->common.R.inc : Gin->header->inc), 2, double);
if (!full_360) { /* Sanity checks on x-domain if not geographic */
if (wesn[XLO] < Gin->header->wesn[XLO]) error = true;
if (wesn[XHI] > Gin->header->wesn[XHI]) error = true;
}
/* Sanity checks on y-domain */
if (wesn[YLO] < Gin->header->wesn[YLO]) error = true;
if (wesn[YHI] > Gin->header->wesn[YHI]) error = true;
if (error) {
GMT_Report (API, GMT_MSG_ERROR, "Output grid domain incompatible with input grid domain.\n");
Return (GMT_RUNTIME_ERROR);
}
/* Allocate space and determine the header for the new grid; croak if there are issues. */
if ((Gout = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, wesn, inc, \
!one_or_zero, GMT_NOTSET, NULL)) == NULL) Return (API->error);
/* We can save time by computing a weight matrix once [or once pr scanline] only
if output grid spacing is a multiple of input grid spacing */
fast_way = (fabs (fmod (Gout->header->inc[GMT_X] / Gin->header->inc[GMT_X], 1.0)) < GMT_CONV4_LIMIT && fabs (fmod (Gout->header->inc[GMT_Y] / Gin->header->inc[GMT_Y], 1.0)) < GMT_CONV4_LIMIT);
same_grid = !(GMT->common.R.active[RSET] || GMT->common.R.active[ISET] || Gin->header->registration == one_or_zero);
if (!fast_way) { /* Not optimal... */
if (Ctrl->F.custom) {
GMT_Report (API, GMT_MSG_ERROR, "Option -F: For -Ff or -Fo the input and output grids must be coregistered.\n");
Return (GMT_RUNTIME_ERROR);
}
GMT_Report (API, GMT_MSG_WARNING, "Your output grid spacing is such that filter-weights must\n");
GMT_Report (API, GMT_MSG_WARNING, "be recomputed for every output node, so expect this run to be slow. Calculations\n");
GMT_Report (API, GMT_MSG_WARNING, "can be speeded up significantly if output grid spacing is chosen to be a multiple\n");
GMT_Report (API, GMT_MSG_WARNING, "of the input grid spacing. If the odd output grid is necessary, consider using\n");
GMT_Report (API, GMT_MSG_WARNING, "a \'fast\' grid for filtering and then resample onto your desired grid with grdsample.\n");
}
if (Ctrl->N.mode == NAN_REPLACE && !same_grid) {
GMT_Report (API, GMT_MSG_ERROR, "-Nr requires co-registered input/output grids, option is ignored\n");
Ctrl->N.mode = NAN_IGNORE;
}
col_origin = gmt_M_memory (GMT, NULL, Gout->header->n_columns, int);
if (!fast_way) x_shift = gmt_M_memory (GMT, NULL, Gout->header->n_columns, double);
if (fast_way && Gin->header->registration == one_or_zero) { /* multiple of grid spacings but one is pix, other is grid so adjust for 1/2 cell */
x_fix = 0.5 * Gin->header->inc[GMT_X];
y_fix = 0.5 * Gin->header->inc[GMT_Y];
}
if (Ctrl->D.mode > GRDFILTER_XY_CARTESIAN) { /* Data on a sphere so must check for both periodic and polar wrap-arounds */
spherical = true;
/* Compute the wrap-around delta_nx to use [may differ from n_columns unless a 360 grid] */
nx_wrap = (int)gmt_M_get_n (GMT, 0.0, 360.0, Gin->header->inc[GMT_X], GMT_GRID_PIXEL_REG); /* So we basically bypass the duplicate point at east */
}
if ((A = init_area_weights (GMT, Gin, Ctrl->D.mode, Ctrl->A.file)) == NULL) Return (API->error); /* Precalculate area weights, optionally save debug grid */
gmt_M_memset (&F, 1, struct FILTER_INFO);
/* Set up the distance scalings for lon and lat, and assign pointer to distance function */
if (Ctrl->F.custom) { /* Get filter-weight grid rather than compute one */
if ((Fin = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->F.file, NULL)) == NULL) { /* Get filter-weight grid */
Return (API->error);
}
F.n_columns = Fin->header->n_columns; F.n_rows = Fin->header->n_rows;
if ((F.n_columns % 2) == 0 || (F.n_rows % 2) == 0) {
GMT_Report (API, GMT_MSG_ERROR, "-Ff|o requires an odd number of rows and columns for filter|operator weight grid,\n");
Return (API->error);
}
F.x_half_width = (F.n_columns - 1) / 2; F.y_half_width = (F.n_rows - 1) / 2;
}
switch (Ctrl->D.mode) {
case GRDFILTER_XY_CARTESIAN: /* Plain, unscaled isotropic Cartesian distances */
x_scale = y_scale = 1.0;
F.radius_func = &CartRadius;
break;
case GRDFILTER_GEO_CARTESIAN: /* Plain, scaled (degree to km) isotropic Cartesian distances */
x_scale = y_scale = GMT->current.proj.DIST_KM_PR_DEG;
F.radius_func = &CartScaledRadius;
break;
case GRDFILTER_GEO_FLATEARTH1: /* Flat Earth Cartesian distances, xscale fixed at mid latitude */
x_scale = GMT->current.proj.DIST_KM_PR_DEG * cosd (0.5 * (Gout->header->wesn[YHI] + Gout->header->wesn[YLO]));
y_scale = GMT->current.proj.DIST_KM_PR_DEG;
F.radius_func = &FlatEarthRadius;
break;
case GRDFILTER_GEO_FLATEARTH2: /* Flat Earth Cartesian distances, xscale reset for each latitude (xscale here is max scale for max |lat|) */
x_scale = GMT->current.proj.DIST_KM_PR_DEG * ((fabs (Gout->header->wesn[YLO]) > Gout->header->wesn[YHI]) ? cosd (Gout->header->wesn[YLO]) : cosd (Gout->header->wesn[YHI]));
y_scale = GMT->current.proj.DIST_KM_PR_DEG;
F.radius_func = &FlatEarthRadius;
break;
case GRDFILTER_GEO_SPHERICAL: /* Great circle distances; set scale using input grid south/north */
x_scale = GMT->current.proj.DIST_KM_PR_DEG * ((fabs (Gin->header->wesn[YLO]) > Gin->header->wesn[YHI]) ? cosd (Gin->header->wesn[YLO]) : cosd (Gin->header->wesn[YHI]));
y_scale = GMT->current.proj.DIST_KM_PR_DEG;
F.radius_func = &SphericalRadius;
break;
case GRDFILTER_GEO_MERCATOR: /* Great circle distances with Mercator coordinates */
/* Get the max |lat| extent of the grid */
max_lat = IMG2LAT (MAX (fabs (Gin->header->wesn[YLO]), fabs (Gin->header->wesn[YHI])));
merc_range = LAT2IMG (max_lat + (0.5 * Ctrl->F.width / GMT->current.proj.DIST_KM_PR_DEG)) - LAT2IMG (max_lat);
x_scale = y_scale = 0.5 * Ctrl->F.width / merc_range;
F.radius_func = &SphericalRadius;
break;
}
/* Assign filter_type number */
for (filter_type = 0; filter_type < GRDFILTER_N_FILTERS && filter_code[filter_type] != Ctrl->F.filter; filter_type++);
switch (filter_type) { /* Set special weight functions for some filters */
case GRDFILTER_COSINE: /* Cosine-bell filter weights */
par[GRDFILTER_INV_R_SCALE] = (Ctrl->F.rect) ? 1.0 : 2.0 / Ctrl->F.width;
F.weight_func = &CosBellWeight;
break;
case GRDFILTER_GAUSSIAN: /* Gaussian filter weights */
par[GRDFILTER_INV_R_SCALE] = (Ctrl->F.rect) ? -4.5 : -18.0 / (Ctrl->F.width * Ctrl->F.width);
F.weight_func = &GaussianWeight;
break;
default: /* Everything else uses unit weights or nothing */
F.weight_func = &UnitWeight;
break;
}
/* Set up miscellaneous filter parameters needed when computing the weights */
par[GRDFILTER_WIDTH] = Ctrl->F.width;
par[GRDFILTER_HALF_WIDTH] = 0.5 * Ctrl->F.width;
par[GRDFILTER_X_SCALE] = x_scale;
par[GRDFILTER_Y_SCALE] = (Ctrl->D.mode == GRDFILTER_GEO_MERCATOR) ? GMT->current.proj.DIST_KM_PR_DEG : y_scale;
F.d_flag = Ctrl->D.mode;
F.rect = Ctrl->F.rect;
F.dx = Gin->header->inc[GMT_X];
F.dy = Gin->header->inc[GMT_Y];
F.y_min = Gin->header->wesn[YLO];
F.y_max = Gin->header->wesn[YHI];
x_width = Ctrl->F.width / (Gin->header->inc[GMT_X] * x_scale);
y_width = ((F.rect) ? Ctrl->F.width2 : Ctrl->F.width) / (Gin->header->inc[GMT_Y] * y_scale);
if (!Ctrl->F.custom) { /* Parameters computed from width and other settings */
F.x_half_width = irint (ceil (x_width / 2.0));
F.y_half_width = irint (ceil (y_width / 2.0));
F.n_columns = 2 * F.x_half_width + 1;
F.n_rows = 2 * F.y_half_width + 1;
if (gmt_M_is_zero (x_scale) || F.x_half_width < 0 || F.n_columns > Gin->header->n_columns) { /* Safety valve when x_scale -> 0.0 */
F.n_columns = Gin->header->n_columns;
F.x_half_width = (F.n_columns - 1) / 2;
if ((F.n_columns - 2 * F.x_half_width - 1) > 0) F.x_half_width++; /* When n_columns is even we may come up short by 1 */
F.n_columns = 2 * F.x_half_width + 1;
}
if (gmt_M_is_zero (y_scale) || F.y_half_width < 0 || F.n_rows > Gin->header->n_rows) { /* Safety valve when y_scale -> 0.0 */
F.n_rows = Gin->header->n_rows;
F.y_half_width = (F.n_rows - 1) / 2;
}
}
visit_check = ((2 * F.x_half_width + 1) >= (int)Gin->header->n_columns); /* Must make sure we only visit each node once along a row */
F.x = gmt_M_memory (GMT, NULL, F.x_half_width+1, double);
F.y = gmt_M_memory (GMT, NULL, F.y_half_width+1, double);
if (visit_check) F.visit = gmt_M_memory (GMT, NULL, Gin->header->n_columns, char);
for (ii = 0; ii <= F.x_half_width; ii++) F.x[ii] = ii * F.dx;
for (jj = 0; jj <= F.y_half_width; jj++) F.y[jj] = jj * F.dy;
weight = gmt_M_memory (GMT, NULL, F.n_columns*F.n_rows, double); /* Allocate space for convolution grid */
if (Ctrl->F.custom) { /* Read convolution grid from file */
ij_wt = 0; wt_sum = 0.0;
gmt_M_grd_loop (GMT, Fin, row_in, col_in, ij_in) { /* Just copy over to weight array while skipping the padding */
weight[ij_wt++] = Fin->data[ij_in];
wt_sum += Fin->data[ij_in];
}
if (GMT_Destroy_Data (API, &Fin) != GMT_NOERROR) { /* Done with this grid */
Return (API->error);
}
if (gmt_M_is_zero (wt_sum)) { /* The custom filter is an operator; should have used -Fo */
GMT_Report (API, GMT_MSG_WARNING, "Your custom filter weights sum to zero; switching to -Fo operator mode.\n");
Ctrl->F.operator = true;
}
if (Ctrl->F.operator) get_weight_sum = false; /* As weights sum to zero we don't want to add them up and divide */
}
if (filter_type > GRDFILTER_OPERATOR) { /* These filters are not convolutions; they require sorting or comparisons */
slow = true;
last_median = 0.5 * (Gin->header->z_min + Gin->header->z_max); /* Initial guess */
if (Ctrl->D.mode > GRDFILTER_XY_CARTESIAN && (filter_type == GRDFILTER_MEDIAN || filter_type == GRDFILTER_MODE)) { /* Spherical (weighted) median/modes requires even more work */
slower = true;
filter_type = (filter_type == GRDFILTER_MEDIAN) ? GRDFILTER_MEDIAN_SPH : GRDFILTER_MODE_SPH; /* Set to the weighted versions */
}
}
if (tid == 0) { /* First or only thread */
GMT_Report (API, GMT_MSG_INFORMATION, "Input n_columns,n_rows = (%d %d), output n_columns,n_rows = (%d %d), filter (max)n_columns,n_rows = (%d %d)\n",
Gin->header->n_columns, Gin->header->n_rows, Gout->header->n_columns, Gout->header->n_rows, F.n_columns, F.n_rows);
if (Ctrl->F.quantile != 0.5)
GMT_Report (API, GMT_MSG_INFORMATION, "Filter type is %s [using %g%% quantile].\n",
filter_name[filter_type], 100.0 * Ctrl->F.quantile);
else
GMT_Report (API, GMT_MSG_INFORMATION, "Filter type is %s.\n", filter_name[filter_type]);
#ifdef HAVE_GLIB_GTHREAD
GMT_Report (API, GMT_MSG_INFORMATION, "Calculations will be distributed over %d threads.\n", Ctrl->z.n_threads);
#endif
}
#ifdef DEBUG
if (Ctrl->A.active) { /* Picked a point to examine filter weights etc */
Ctrl->A.COL = gmt_M_grd_x_to_col (GMT, Ctrl->A.x, Gout->header);
Ctrl->A.ROW = gmt_M_grd_y_to_row (GMT, Ctrl->A.y, Gout->header);
if (Ctrl->A.mode == 'r') F.debug = true; /* In order to return radii instead of weights */
if (tid == 0) GMT_Report (API, GMT_MSG_INFORMATION, "ROW = %d COL = %d\n", Ctrl->A.ROW, Ctrl->A.COL);
}
#endif
/* Compute nearest xoutput i-indices and shifts once */
for (col_out = 0; col_out < Gout->header->n_columns; col_out++) {
x_out = gmt_M_grd_col_to_x (GMT, col_out, Gout->header); /* Current longitude */
col_origin[col_out] = (int)gmt_M_grd_x_to_col (GMT, x_out, Gin->header);
if (!fast_way) x_shift[col_out] = x_out - gmt_M_grd_col_to_x (GMT, col_origin[col_out], Gin->header);
}
/* Determine how much effort to compute weights:
0 = Weights read from custom weight grid
1 = Compute weights once for entire grid
2 = Compute weights once per scanline [slow]
3 = Compute weights for every output point [slower]
*/
if (Ctrl->F.custom)
effort_level = 0;
else if (fast_way && Ctrl->D.mode <= GRDFILTER_GEO_FLATEARTH1)
effort_level = 1;
else if (fast_way && Ctrl->D.mode >= GRDFILTER_GEO_FLATEARTH2)
effort_level = 2;
else
effort_level = 3;
if (effort_level == 1) set_weight_matrix (GMT, &F, weight, 0.0, par, x_fix, y_fix);
#ifdef DEBUG
if (Ctrl->A.active) for (ij_in = 0; ij_in < Gin->header->size; ij_in++) Gin->data[ij_in] = 0.0f; /* We are using Gin to store filter weights etc instead */
#endif
gmt_M_tic(GMT);
#ifdef HAVE_GLIB_GTHREAD
if (Ctrl->z.n_threads > 1)
threads = gmt_M_memory (GMT, NULL, Ctrl->z.n_threads, GThread *);
#endif
threadArg = gmt_M_memory (GMT, NULL, Ctrl->z.n_threads, struct THREAD_STRUCT);
for (i = 0; i < Ctrl->z.n_threads; i++) {
threadArg[i].GMT = GMT;
threadArg[i].Ctrl = Ctrl;
threadArg[i].Gin = Gin;
threadArg[i].Gout = Gout;
threadArg[i].A = A;
threadArg[i].F = F;
threadArg[i].weight = weight;
threadArg[i].x_shift = x_shift;
threadArg[i].col_origin = col_origin;
threadArg[i].par = par;
threadArg[i].x_fix = x_fix;
threadArg[i].y_fix = y_fix;
threadArg[i].last_median= last_median;
threadArg[i].fast_way = fast_way;
threadArg[i].spherical = spherical;
threadArg[i].slow = slow;
threadArg[i].slower = slower;
threadArg[i].nx_wrap = nx_wrap;
threadArg[i].effort_level = effort_level;
threadArg[i].filter_type = filter_type;
threadArg[i].r_start = i * irint((Gout->header->n_rows) / Ctrl->z.n_threads);
threadArg[i].thread_num = i;
if (Ctrl->z.n_threads == 1) { /* Independently of WITH_THREADS, if only one don't call the threading machine */
threadArg[i].r_stop = Gout->header->n_rows;
grdfilter_threaded_function (&threadArg[0]);
break; /* Make sure we don't go through the threads lines below */
}
#ifndef HAVE_GLIB_GTHREAD
}
#else
threadArg[i].r_stop = (i + 1) * irint((Gout->header->n_rows) / Ctrl->z.n_threads);
if (i == Ctrl->z.n_threads - 1) threadArg[i].r_stop = Gout->header->n_rows; /* Make sure last row is not left behind */
threads[i] = g_thread_new(NULL, grdfilter_thread_function, (void*)&(threadArg[i]));
}
if (Ctrl->z.n_threads > 1) { /* Otherwise g_thread_new was never called aand so no need to "join" */
for (i = 0; i < Ctrl->z.n_threads; i++)
g_thread_join(threads[i]);
}
if (Ctrl->z.n_threads > 1)
gmt_M_free (GMT, threads);
#endif
gmt_M_free (GMT, threadArg);
gmt_M_toc(GMT,""); /* Print total run time, but only if -Vt was set */
gmt_M_free (GMT, weight);
gmt_M_free (GMT, F.x);
gmt_M_free (GMT, F.y);
if (visit_check) gmt_M_free (GMT, F.visit);
gmt_M_free (GMT, col_origin);
if (!fast_way) gmt_M_free (GMT, x_shift);