-
Notifications
You must be signed in to change notification settings - Fork 337
/
filter1d.c
1112 lines (999 loc) · 46.7 KB
/
filter1d.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-2020 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
*--------------------------------------------------------------------*/
/*
* API functions to support the filter1d application.
*
* Author: Walter H. F. Smith
* Date: 1-JAN-2010
* Version: 6 API
*
* Brief synopsis: filter1d_func will read N columns of data from file
* or stdin and return filtered output at user-selected positions.
* The time variable can be in any specified column of the input. Several
* filters are available, with robustness as an option.
*
* Filters:
* Convolutions: Boxcar, Cosine Arch, Gaussian, Median, or Mode.
* Geospatial: Median, Mode, Extreme values.
* Robust: Option replaces outliers with medians in filter.
* Output: At input times, or from t_start to t_stop by t_int.
* Lack: Option checks for data gaps in input series.
* Symmetry: Option checks for asymmetry in filter window.
* Quality: Option checks for low mean weight in window.
*
*/
#include "gmt_dev.h"
#define THIS_MODULE_CLASSIC_NAME "filter1d"
#define THIS_MODULE_MODERN_NAME "filter1d"
#define THIS_MODULE_LIB "core"
#define THIS_MODULE_PURPOSE "Time domain filtering of 1-D data tables"
#define THIS_MODULE_KEYS "<D{,>D},FD(1"
#define THIS_MODULE_NEEDS ""
#define THIS_MODULE_OPTIONS "-:>Vabdefghijoq" GMT_OPT("HMm")
/* Control structure for filter1d */
struct FILTER1D_CTRL {
struct FILTER1D_D { /* -D<inc> */
bool active;
double inc;
} D;
struct FILTER1D_E { /* -E */
bool active;
} E;
struct FILTER1D_F { /* -F<type><width>[<mode>] */
bool active;
bool highpass;
bool variable;
char filter; /* Character codes for the filter */
double width;
int mode; /* -1/0/+1 */
char *file; /* File with custom filter coefficients or w(t) variable widths */
} F;
struct FILTER1D_L { /* -L<lackwidth> */
bool active;
double value;
} L;
struct FILTER1D_N { /* -N<t_col> or -Nc|<unit>[+a] */
bool active;
bool add_col;
char unit;
unsigned int mode;
unsigned int spatial;
int col;
} N;
struct FILTER1D_Q { /* -Q<factor> */
bool active;
double value;
} Q;
struct FILTER1D_S { /* -S<symmetry> */
bool active;
double value;
} S;
struct FILTER1D_T { /* -T<tmin/tmax/tinc>[+a|n] */
bool active;
struct GMT_ARRAY T;
} T;
};
#define FILTER1D_BOXCAR 0
#define FILTER1D_COS_ARCH 1
#define FILTER1D_GAUSSIAN 2
#define FILTER1D_CUSTOM 3
#define FILTER1D_MEDIAN 4
#define FILTER1D_MODE 5
#define FILTER1D_LOWER_ALL 6
#define FILTER1D_LOWER_POS 7
#define FILTER1D_UPPER_ALL 8
#define FILTER1D_UPPER_NEG 9
#define FILTER1D_N_FILTERS 10
#define FILTER1D_CONVOLVE 3 /* If filter_type > FILTER1D_CONVOLVE then a FILTER1D_MEDIAN, FILTER1D_MODE, or EXTREME filter is selected */
struct FILTER1D_INFO { /* Control structure for all aspects of the filter setup */
bool use_ends; /* True to start/stop at ends of series instead of 1/2 width inside */
bool check_asym; /* true to test whether the data are asymmetric about the output time */
bool check_lack; /* true to test for lack of data (gap) in the filter window */
bool check_q; /* true to test average weight or N in median */
bool robust; /* Look for outliers in data when true */
bool equidist; /* Data is evenly sampled in t */
bool out_at_time; /* true when output is required at evenly spaced intervals */
bool f_operator; /* true if custom weights coefficients sum to zero */
bool highpass; /* true if we are doing a highpass filter */
bool variable; /* true if we are doing a variable width filter */
uint64_t *n_this_col; /* Pointer to array of counters [one per column] */
uint64_t *n_left; /* Pointer to array of counters [one per column] */
uint64_t *n_right; /* Pointer to array of counters [one per column] */
uint64_t n_cols; /* Number of columns of input */
uint64_t t_col; /* Column of time abscissae (independent variable) */
uint64_t n_f_wts; /* Number of filter weights */
uint64_t half_n_f_wts; /* Half the number of filter weights */
uint64_t n_rows; /* Number of rows of input */
size_t n_row_alloc; /* Number of rows of data to allocate */
size_t n_work_alloc; /* Number of rows of workspace to allocate */
int filter_type; /* Flag indicating desired filter type */
int kind; /* -1 skip +ve, +1 skip -ve, else use all [for the l|L|u|U filter] */
int way; /* -1 find minimum, +1 find maximum [for the l|L|u|U filter] */
unsigned int mode_selection;
unsigned int n_multiples;
unsigned int nw; /* Length of variable filter width array */
double *ft; /* Pointer for array of variable filter times */
double *fw; /* Pointer for array of variable filter widths */
double *f_wt; /* Pointer for array of filter coefficients */
double *min_loc; /* Pointer for array of values, one per [column] */
double *max_loc;
double *last_loc;
double *this_loc;
double *min_scl;
double *max_scl;
double *last_scl;
double *this_scl;
double **work; /* Pointer to array of pointers to doubles for work */
double **data; /* Pointer to array of pointers to doubles for data */
double dt; /* Delta time resolution for filter computation */
double q_factor; /* Quality level for mean weights or n in median */
double filter_width; /* Full width of filter in user's units */
double half_width;
double t_start; /* x-value of first output point */
double t_stop; /* x-value of last output point */
double t_start_t; /* user specified x-value of first output point if out_at_time == true */
double t_stop_t; /* user specified x-value of last output point if out_at_time == true */
double t_int; /* Output interval */
double sym_coeff; /* Symmetry coefficient */
double lack_width; /* Lack of data width */
double extreme; /* Extreme value [for the l|L|u|U filter] */
struct GMT_DATASET *Fin; /* Pointer to table with custom weight coefficients (optional) */
struct GMT_DATASET *W; /* Pointer to table with variable filter widths (optional) */
struct GMT_ARRAY T;
};
static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
struct FILTER1D_CTRL *C;
C = gmt_M_memory (GMT, NULL, 1, struct FILTER1D_CTRL);
/* Initialize values whose defaults are not 0/false/NULL */
return (C);
}
static void Free_Ctrl (struct GMT_CTRL *GMT, struct FILTER1D_CTRL *C) { /* Deallocate control structure */
if (!C) return;
gmt_M_str_free (C->F.file);
gmt_free_array (GMT, &(C->T.T));
gmt_M_free (GMT, C);
}
static int 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 [<table>] -F<type><width>[<modifiers>] [-D<increment>] [-E]\n", name);
GMT_Message (API, GMT_TIME_NONE, "\t[-L<lack_width>] [-N<t_col>] [-Q<q_factor>] [-S<symmetry>] [-T[<min>/<max>/]<inc>[+e|n|a]]\n");
GMT_Message (API, GMT_TIME_NONE, "\t[%s] [%s] [%s] [%s] [%s]\n\t[%s] [%s]\n\t[%s] [%s] [%s]\n\t[%s] [%s] [%s]\n\n",
GMT_V_OPT, GMT_b_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_g_OPT, GMT_h_OPT, GMT_i_OPT, GMT_j_OPT, GMT_o_OPT, GMT_q_OPT, GMT_colon_OPT, GMT_PAR_OPT);
if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
GMT_Message (API, GMT_TIME_NONE, "\t-F Set filtertype. Choose from convolution and non-convolution filters\n");
GMT_Message (API, GMT_TIME_NONE, "\t and append filter <width> in same units as time column.\n");
GMT_Message (API, GMT_TIME_NONE, "\t If <width> is a file (and <type> is not f) it must have a time-series of filter widths.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Append +h select high-pass filtering [Default is low-pass filtering].\n");
GMT_Message (API, GMT_TIME_NONE, "\t Convolution filters:\n");
GMT_Message (API, GMT_TIME_NONE, "\t b: Boxcar : Weights are equal.\n");
GMT_Message (API, GMT_TIME_NONE, "\t c: Cosine arch : Weights given by cosine arch.\n");
GMT_Message (API, GMT_TIME_NONE, "\t g: Gaussian : Weights given by Gaussian function.\n");
GMT_Message (API, GMT_TIME_NONE, "\t f<name>: Custom : Weights given in one-column file <name>.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Non-convolution filters:\n");
GMT_Message (API, GMT_TIME_NONE, "\t m: Median : Return the median value.\n");
GMT_Message (API, GMT_TIME_NONE, "\t p: Maximum likelihood probability (mode) estimator : Return the mode.\n");
GMT_Message (API, GMT_TIME_NONE, "\t By default, we return the average mode if more than one is found.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Append +l to return the lowest mode if multiple modes are found [return average].\n");
GMT_Message (API, GMT_TIME_NONE, "\t Append +u to return the uppermost mode if multiple modes are found [return average].\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 positive points.\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 negative points.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Upper case type B, C, G, M, P, F will use robust filter versions,\n");
GMT_Message (API, GMT_TIME_NONE, "\t i.e., replace outliers (2.5 L1 scale (MAD) of median) with median during filtering.\n");
GMT_Message (API, GMT_TIME_NONE, "\n\tOPTIONS:\n");
GMT_Option (API, "<");
GMT_Message (API, GMT_TIME_NONE, "\t-D Set fixed increment when series is NOT equidistantly sampled.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Then <increment> will be the abscissae resolution, i.e., all abscissae\n");
GMT_Message (API, GMT_TIME_NONE, "\t will be rounded off to a multiple of <increment>.\n");
GMT_Message (API, GMT_TIME_NONE, "\t-E Include ends of time series in output [Default loses half_width at each end].\n");
GMT_Message (API, GMT_TIME_NONE, "\t-L Check for lack of data condition. If input data has a gap exceeding\n");
GMT_Message (API, GMT_TIME_NONE, "\t <width> then no output will be given at that point [Default does not check Lack].\n");
GMT_Message (API, GMT_TIME_NONE, "\t-N Set the column that contains the independent variable (time) [0].\n");
GMT_Message (API, GMT_TIME_NONE, "\t The left-most column is # 0, the right-most is # (<n_cols> - 1).\n");
GMT_Message (API, GMT_TIME_NONE, "\t-Q Assess quality of output value by checking mean weight in convolution.\n");
GMT_Message (API, GMT_TIME_NONE, "\t Enter <q_factor> between 0 and 1. If mean weight < q_factor, output is\n");
GMT_Message (API, GMT_TIME_NONE, "\t suppressed at this point [Default does not check quality].\n");
GMT_Message (API, GMT_TIME_NONE, "\t-S Check symmetry of data about window center. Enter a factor\n");
GMT_Message (API, GMT_TIME_NONE, "\t between 0 and 1. If ( (abs(n_left - n_right)) / (n_left + n_right) ) > factor,\n");
GMT_Message (API, GMT_TIME_NONE, "\t then no output will be given at this point [Default does not check Symmetry].\n");
GMT_Message (API, GMT_TIME_NONE, "\t-T Make evenly spaced output time steps from <min> to <max> by <inc> [Default uses input times].\n");
GMT_Message (API, GMT_TIME_NONE, "\t Append +n to indicate <inc> is the number of t-values to produce instead.\n");
GMT_Message (API, GMT_TIME_NONE, "\t If only <inc> is given, optionally append +e to keep increment exact [Default will adjust to fit range].\n");
GMT_Message (API, GMT_TIME_NONE, "\t For absolute time filtering, append a valid time unit (%s) to the increment.\n", GMT_TIME_UNITS_DISPLAY);
GMT_Message (API, GMT_TIME_NONE, "\t For spatial filtering with distance computed from the first two columns, specify increment as\n");
GMT_Message (API, GMT_TIME_NONE, "\t <inc> and append a geospatial distance unit (%s) or c (for Cartesian distances).\n", GMT_LEN_UNITS_DISPLAY);
GMT_Message (API, GMT_TIME_NONE, "\t Optionally, append +a to add such internal distances as a final output column [no distances added].\n");
GMT_Message (API, GMT_TIME_NONE, "\t Alternatively, give a file with output times in the first column or a comma-separated list.\n");
GMT_Option (API, "V,bi,bo,d,e,f,g,h,i,j,o,q,:,.");
return (GMT_MODULE_USAGE);
}
GMT_LOCAL char filter1d_set_unit_and_mode (const char *arg, unsigned int *mode) {
unsigned int k = 0;
*mode = GMT_GREATCIRCLE; /* Default is great circle distances */
switch (arg[0]) {
case '-': *mode = GMT_FLATEARTH; k = 1; break;
case '+': *mode = GMT_GEODESIC; k = 1; break;
}
return (arg[k]);
}
static int parse (struct GMT_CTRL *GMT, struct FILTER1D_CTRL *Ctrl, struct GMT_OPTION *options) {
/* This parses the options provided to filter1d and sets parameters in CTRL.
* 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;
int sval = 0;
char *c = NULL, p, txt[GMT_LEN64] = {""}, *t_arg = NULL;
struct GMT_OPTION *opt = NULL;
struct GMTAPI_CTRL *API = GMT->parent;
for (opt = options; opt; opt = opt->next) {
switch (opt->option) {
case '<': /* Skip input files */
if (GMT_Get_FilePath (GMT->parent, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;;
break;
/* Processes program-specific parameters */
case 'D': /* Get fixed increment */
Ctrl->D.inc = atof (opt->arg);
Ctrl->D.active = true;
break;
case 'E': /* Include ends of series */
Ctrl->E.active = true;
break;
case 'F': /* Filter selection */
Ctrl->F.active = true;
if (opt->arg[0] && strchr ("BbCcGgMmPpLlUuFf", opt->arg[0])) { /* OK filter code */
Ctrl->F.filter = opt->arg[0];
if ((c = strstr (opt->arg, "+h"))) { /* Want high-pass filter */
Ctrl->F.highpass = true;
c[0] = '\0'; /* Chop off +h */
}
if (gmt_access (GMT, &opt->arg[1], R_OK))
Ctrl->F.width = atof (&opt->arg[1]);
else if (Ctrl->F.filter != 'f') { /* Got variable filter width series */
Ctrl->F.variable = true;
Ctrl->F.file = strdup (&opt->arg[1]);
}
switch (Ctrl->F.filter) { /* Get some further info from some filters */
case 'P':
case 'p':
p = opt->arg[strlen(opt->arg-1)]; /* Last character */
if (strstr (opt->arg, "+l")) Ctrl->F.mode = -1;
else if (strstr (opt->arg, "+u")) Ctrl->F.mode = +1;
else if (p == '-') Ctrl->F.mode = -1; /* Old style */
else if (p == '+') Ctrl->F.mode = +1; /* Old style */
break;
case 'F':
case 'f':
Ctrl->F.width = DBL_MAX; /* To avoid range test errors before reading coefficients */
if (opt->arg[1] && !gmt_access (GMT, &opt->arg[1], R_OK))
Ctrl->F.file = strdup (&opt->arg[1]);
else {
GMT_Report (API, GMT_MSG_ERROR, "Option -F[Ff]: Could not find file %s.\n", &opt->arg[1]);
++n_errors;
}
break;
}
n_errors += gmt_M_check_condition (GMT, Ctrl->F.file == NULL && Ctrl->F.width <= 0.0,
"Option -F: Filterwidth must be positive\n");
}
else {
GMT_Report (API, GMT_MSG_ERROR, "Option -F: Correct syntax: -FX<width>, X one of BbCcGgMmPpFflLuU\n");
++n_errors;
}
break;
case 'I': /* DEPRECATED: Activate the ignore option. This is now -di<value> and happens during reading */
sprintf (txt, "i%s", opt->arg);
n_errors += gmt_parse_d_option (GMT, txt);
break;
case 'L': /* Check for lack of data */
Ctrl->L.active = true;
Ctrl->L.value = atof (opt->arg);
break;
case 'N': /* Select column with independent coordinate [0] */
Ctrl->N.active = true;
if (gmt_M_compat_check (GMT, 4)) { /* GMT4 LEVEL */
if (strchr (opt->arg, '/')) { /* Gave obsolete format */
int sval0;
GMT_Report (API, GMT_MSG_COMPAT, "-N<ncol>/<tcol> option is deprecated; use -N<tcol> instead.\n");
if (sscanf (opt->arg, "%d/%d", &sval0, &sval) != 2) {
GMT_Report (API, GMT_MSG_ERROR, "Option -N: Syntax is -N<tcol>\n");
++n_errors;
}
}
else if (!strchr (opt->arg, '/'))
sval = atoi (opt->arg);
}
else
sval = atoi (opt->arg);
if (gmt_M_compat_check (GMT, 5) && (strstr (opt->arg, "+a") || (opt->arg[0] == 'g' || opt->arg[0] == 'c'))) { /* Deprecated syntax */
GMT_Report (API, GMT_MSG_COMPAT, "-Nc|g[+a] option is deprecated; use -T[<min>/<max>/]<int>[+a|n] instead.\n");
if ((c = strstr (opt->arg, "+a"))) { /* Want to output any spatial distances computed */
c[0] = '\0'; /* Chop off the modifier */
Ctrl->N.add_col = true;
}
if ((opt->arg[0] == 'g' || opt->arg[0] == 'c')) { /* Spatial filtering */
Ctrl->N.spatial = 1;
if (opt->arg[0] == 'g') { /* Geospatial filtering */
if (opt->arg[1])
Ctrl->N.unit = filter1d_set_unit_and_mode (&opt->arg[1], &Ctrl->N.mode);
else /* Default to meter and great circle distances */
Ctrl->N.unit = filter1d_set_unit_and_mode ("e", &Ctrl->N.mode);
Ctrl->N.spatial = 2;
}
else
Ctrl->N.unit = 'X'; /* For Cartesian distances */
if (c) c[0] = '+'; /* Restore */
}
}
else {
n_errors += gmt_M_check_condition (GMT, sval < 0, "Option -N: Time column cannot be negative.\n");
Ctrl->N.col = sval;
}
break;
case 'Q': /* Assess quality of output */
Ctrl->Q.value = atof (opt->arg);
Ctrl->Q.active = true;
break;
case 'S': /* Activate symmetry test */
Ctrl->S.active = true;
Ctrl->S.value = atof (opt->arg);
break;
case 'T': /* Set output knots */
Ctrl->T.active = true;
t_arg = opt->arg;
break;
default: /* Report bad options */
n_errors += gmt_default_error (GMT, opt->option);
break;
}
}
if (Ctrl->T.active) /* Do this one here since we need Ctrl->N.col to be set first, if selected */
n_errors += gmt_parse_array (GMT, 'T', t_arg, &(Ctrl->T.T), GMT_ARRAY_TIME | GMT_ARRAY_DIST | GMT_ARRAY_ROUND | GMT_ARRAY_UNIQUE, Ctrl->N.col);
if (Ctrl->N.spatial) { /* Obsolete -N settings propagated to -T */
Ctrl->T.T.spatial = Ctrl->N.spatial;
Ctrl->T.T.unit = Ctrl->N.unit;
Ctrl->T.T.distmode = Ctrl->N.mode;
}
if (Ctrl->N.add_col) Ctrl->T.T.add = true; /* Obsolete -N+a settings propagated to -T */
/* Check arguments */
n_errors += gmt_M_check_condition (GMT, !Ctrl->F.active, "Option -F is required\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && Ctrl->D.inc <= 0.0, "Option -D: must give positive increment\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->T.active && Ctrl->T.T.set == 3 && (Ctrl->T.T.max - Ctrl->T.T.min) < Ctrl->F.width, "Option -T: Output interval < filterwidth\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->L.active && (Ctrl->L.value < 0.0 || Ctrl->L.value > Ctrl->F.width) , "Option -L: Unreasonable lack-of-data interval\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->S.active && (Ctrl->S.value < 0.0 || Ctrl->S.value > 1.0) , "Option -S: Enter a factor between 0 and 1\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->Q.active && (Ctrl->Q.value < 0.0 || Ctrl->Q.value > 1.0), "Option -Q: Enter a factor between 0 and 1\n");
return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
}
/* Various functions which will be accessed via pointers depending on chosen filter */
GMT_LOCAL double filter1d_boxcar_weight (double radius, double half_width) {
return ((radius > half_width) ? 0.0 : 1.0);
}
GMT_LOCAL double filter1d_cosine_weight (double radius, double half_width) {
return ((radius > half_width) ? 0.0 : 1.0 + cos (radius * M_PI / half_width));
}
GMT_LOCAL double filter1d_gaussian_weight (double radius, double half_width) {
return ((radius > half_width) ? 0.0 : exp (-4.5 * radius * radius / (half_width * half_width)));
}
GMT_LOCAL void filter1d_allocate_data_space (struct GMT_CTRL *GMT, struct FILTER1D_INFO *F) {
uint64_t i;
for (i = 0; i < F->n_cols; ++i) F->data[i] = gmt_M_memory (GMT, F->data[i], F->n_row_alloc, double);
}
GMT_LOCAL void filter1d_allocate_more_work_space (struct GMT_CTRL *GMT, struct FILTER1D_INFO *F) {
uint64_t i;
for (i = 0; i < F->n_cols; ++i) F->work[i] = gmt_M_memory (GMT, F->work[i], F->n_work_alloc, double);
}
GMT_LOCAL int filter1d_set_up_filter (struct GMT_CTRL *GMT, struct FILTER1D_INFO *F) {
uint64_t i, i1, i2;
double t_0, t_1, time, w_sum;
double (*get_weight[3]) (double, double); /* Pointers to desired weight function. */
t_0 = F->data[F->t_col][0];
t_1 = F->data[F->t_col][F->n_rows-1];
if (F->equidist) F->dt = (t_1 - t_0) / (F->n_rows - 1);
if (F->filter_type == FILTER1D_CUSTOM) { /* Use coefficients we read from file */
F->n_f_wts = F->Fin->n_records;
F->f_wt = gmt_M_memory (GMT, F->f_wt, F->n_f_wts, double);
gmt_M_memcpy (F->f_wt, F->Fin->table[0]->segment[0]->data[GMT_X], F->n_f_wts, double);
for (i = 0, w_sum = 0.0; i < F->n_f_wts; ++i) w_sum += F->f_wt[i];
F->f_operator = (gmt_M_is_zero (w_sum)); /* If weights sum to zero it is an operator like {-1 1] or [1 -2 1] */
if (w_sum > 1.0) { /* Must normalize filter weights */
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Must normalize custom filter since weight sum > 1 [%g]\n", w_sum);
for (i = 0; i < F->n_f_wts; ++i) F->f_wt[i] /= w_sum;
}
F->half_n_f_wts = F->n_f_wts / 2;
F->half_width = F->half_n_f_wts * F->dt;
F->filter_width = 2.0 * F->half_width;
}
else if (F->filter_type <= FILTER1D_CONVOLVE) {
get_weight[FILTER1D_BOXCAR] = &filter1d_boxcar_weight;
get_weight[FILTER1D_COS_ARCH] = &filter1d_cosine_weight;
get_weight[FILTER1D_GAUSSIAN] = &filter1d_gaussian_weight;
F->half_width = 0.5 * F->filter_width;
F->half_n_f_wts = lrint (F->half_width / F->dt);
F->n_f_wts = 2 * F->half_n_f_wts + 1;
F->f_wt = gmt_M_memory (GMT, F->f_wt, F->n_f_wts, double);
for (i = 0; i <= F->half_n_f_wts; ++i) {
time = i * F->dt;
i1 = F->half_n_f_wts - i;
i2 = F->half_n_f_wts + i;
F->f_wt[i1] = F->f_wt[i2] = ( *get_weight[F->filter_type]) (time, F->half_width);
}
}
else
F->half_width = 0.5 * F->filter_width;
/* Initialize start/stop time */
if (F->out_at_time) {
/* populate F->t_start and F->t_stop */
/* Here, F->T.set == 1 means we were only given an increment, otherwise we got min/max/inc.
* Thus, for the latter we want whatever stop we pick to be an integer steps of inc from min. */
if (F->T.set == 1 || F->t_start_t < t_0) /* Not set or user defined t_start_t outside bounds */
F->t_start = (F->T.set == 1) ? floor (t_0 / F->t_int) * F->t_int : F->t_start_t + floor ((t_0 - F->t_start_t) / F->t_int) * F->t_int;
else
F->t_start = F->t_start_t;
if (F->T.set == 1 || F->t_stop_t > t_1) /* Not set or user defined t_stop_t outside bounds */
F->t_stop = (F->T.set == 1) ? ceil (t_1 / F->t_int) * F->t_int : F->t_stop_t - floor ((F->t_stop_t - t_1) / F->t_int) * F->t_int;
else
F->t_stop = F->t_stop_t;
if (!F->use_ends) {
/* remove filter half width from bounds */
F->t_start += F->half_width;
F->t_stop -= F->half_width;
}
}
else {
if (F->use_ends) {
F->t_start = t_0;
F->t_stop = t_1;
}
else {
uint64_t row;
double small = (F->data[F->t_col][1] - F->data[F->t_col][0]) * GMT_CONV8_LIMIT;
for (row = 0; (F->data[F->t_col][row] - t_0 + small) < F->half_width; ++row);
F->t_start = F->data[F->t_col][row];
for (row = F->n_rows - 1; row > 0 && (t_1 - F->data[F->t_col][row] + small) < F->half_width; --row);
F->t_stop = F->data[F->t_col][row];
}
}
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "F width: %g Resolution: %g Start: %g Stop: %g\n", F->filter_width, F->dt, F->t_start, F->t_stop);
return (0);
}
GMT_LOCAL int filter1d_lack_check (struct FILTER1D_INFO *F, uint64_t i_col, uint64_t left, uint64_t right) {
uint64_t last_row, this_row;
bool lacking = false;
double last_t;
last_row = left;
while (!(gmt_M_is_dnan (F->data[i_col][last_row])) && last_row < (right - 1)) ++last_row;
last_t = F->data[F->t_col][last_row];
this_row = last_row + 1;
while (!(lacking) && this_row < (right - 1)) {
while (!(gmt_M_is_dnan (F->data[i_col][this_row])) && this_row < (right - 1)) ++this_row;
if ( (F->data[F->t_col][this_row] - last_t) > F->lack_width)
lacking = true;
else {
last_t = F->data[F->t_col][this_row];
++this_row;
}
}
return (lacking);
}
GMT_LOCAL void filter1d_get_robust_estimates (struct GMT_CTRL *GMT, struct FILTER1D_INFO *F, uint64_t j, uint64_t n, int both) {
uint64_t i, n_smooth;
bool sort_me = true;
double low, high, last, temp;
if (F->filter_type > FILTER1D_MODE)
temp = gmt_extreme (GMT, F->work[j], n, F->extreme, F->kind, F->way);
else if (F->filter_type == FILTER1D_MODE) {
n_smooth = n / 2;
gmt_mode (GMT, F->work[j], n, n_smooth, sort_me, F->mode_selection, &F->n_multiples, &temp);
}
else {
low = F->min_loc[j];
high = F->max_loc[j];
last = F->last_loc[j];
gmt_median (GMT, F->work[j], n, low, high, last, &temp);
}
F->last_loc[j] = F->this_loc[j] = temp;
if (both) {
for (i = 0; i < n; ++i) F->work[j][i] = fabs (F->work[j][i] - F->this_loc[j]);
low = F->min_scl[j];
high = F->max_scl[j];
last = F->last_scl[j];
gmt_median (GMT, F->work[j], n, low, high, last, &temp);
F->last_scl[j] = F->this_scl[j] = temp;
}
}
GMT_LOCAL int filter1d_do_the_filter (struct GMTAPI_CTRL *C, struct FILTER1D_INFO *F) {
uint64_t i_row, left, right, n_l, n_r, k = 0, n_for_call, n_good_ones, last_k;
uint64_t iq, i_col, diff;
int64_t i_f_wt, n_in_filter;
int result;
bool *good_one = NULL; /* Pointer to array of logicals [one per column] */
double t_time, delta_time, wt, val, med, scl, small, symmetry;
double *wt_sum = NULL; /* Pointer for array of weight sums [each column] */
double *data_sum = NULL; /* Pointer for array of data * weight sums [columns] */
double *outval = NULL;
struct GMT_RECORD *Out = NULL;
struct GMT_CTRL *GMT = C->GMT;
left = right = 0; /* Left/right end of filter window */
iq = lrint (F->q_factor);
if (F->out_at_time) { /* Set up equidistant output times given the min/max/inc given */
if (gmt_create_array (GMT, 'T', &(F->T), &(F->t_start), &(F->t_stop)))
return GMT_RUNTIME_ERROR;
k = 0; /* Start at first output point */
last_k = F->T.n - 1;
}
else { /* Duplicate an output time array from the input times but exclude ends (unless -E) */
F->T.array = gmt_M_memory (GMT, NULL, F->n_rows, double);
gmt_M_memcpy (F->T.array, F->data[F->t_col], F->n_rows, double);
for (k = 0; F->data[F->t_col][k] < F->t_start; ++k); /* Bypass points outside */
for (last_k = F->n_rows-1; last_k && F->data[F->t_col][last_k] > F->t_stop; --last_k); /* Bypass points outside */
}
small = (F->T.array[1] - F->T.array[0]) * GMT_CONV8_LIMIT;
outval = gmt_M_memory (GMT, NULL, F->n_cols, double);
good_one = gmt_M_memory (GMT, NULL, F->n_cols, bool);
wt_sum = gmt_M_memory (GMT, NULL, F->n_cols, double);
data_sum = gmt_M_memory (GMT, NULL, F->n_cols, double);
Out = gmt_new_record (GMT, NULL, NULL);
while (k <= last_k) {
if (F->variable) { /* Must obtain the filter width for this time */
if (F->T.array[k] < F->ft[0] || F->T.array[k] > F->ft[F->nw-1]) {
k++;
continue; /* Outside filter width array */
}
result = gmt_intpol (GMT, F->ft, F->fw, NULL, F->nw, 1, &(F->T.array[k]), &(F->filter_width), 0.0, GMT->current.setting.interpolant);
filter1d_set_up_filter (GMT, F);
if (F->T.array[k] < F->t_start || F->T.array[k] > F->t_stop) {
k++;
continue; /* Outside full filter window */
}
}
/* Determine data window for convolution/filter */
while (left < F->n_rows && (F->T.array[k] - F->data[F->t_col][left] - small) > F->half_width) ++left;
while (right < F->n_rows && (F->data[F->t_col][right] - F->T.array[k] - small) <= F->half_width) ++right;
n_in_filter = (int64_t)(right - left);
if ( n_in_filter <= 0 || (F->check_lack && ( (F->filter_width / n_in_filter) > F->lack_width) ) ) {
k++; /* Nothing to filter, go to next output time */
continue;
}
t_time = F->T.array[k]; /* Current output time */
for (i_col = 0; i_col < F->n_cols; ++i_col) {
F->n_this_col[i_col] = 0;
wt_sum[i_col] = data_sum[i_col] = 0.0;
if (i_col == F->t_col)
good_one[i_col] = false;
else if (F->check_lack)
good_one[i_col] = !(filter1d_lack_check (F, i_col, left, right));
else
good_one[i_col] = true;
if (F->check_asym) F->n_left[i_col] = F->n_right[i_col] = 0;
}
if (F->robust || F->filter_type > FILTER1D_CONVOLVE) {
if (n_in_filter > (int64_t)F->n_work_alloc) {
F->n_work_alloc = (size_t)n_in_filter;
filter1d_allocate_more_work_space (GMT, F);
}
for (i_row = left; i_row < right; ++i_row) {
for (i_col = 0; i_col < F->n_cols; ++i_col) {
if (!(good_one[i_col])) continue;
if (!gmt_M_is_dnan (F->data[i_col][i_row])) {
F->work[i_col][F->n_this_col[i_col]] = F->data[i_col][i_row];
F->n_this_col[i_col]++;
if (F->check_asym) {
if (F->data[F->t_col][i_row] < t_time) F->n_left[i_col]++;
if (F->data[F->t_col][i_row] > t_time) F->n_right[i_col]++;
}
}
}
}
if (F->check_asym) {
for (i_col = 0; i_col < F->n_cols; ++i_col) {
if (!(good_one[i_col])) continue;
n_l = F->n_left[i_col];
n_r = F->n_right[i_col];
diff = (n_l >= n_r) ? n_l - n_r : n_r - n_l;
symmetry = ((double)diff)/(n_l + n_r);
if (symmetry > F->sym_coeff) good_one[i_col] = false;
}
}
if ((F->filter_type > FILTER1D_CONVOLVE) && F->check_q) {
for (i_col = 0; i_col < F->n_cols; ++i_col) {
if (F->n_this_col[i_col] < iq) good_one[i_col] = false;
}
}
for (i_col = 0; i_col < F->n_cols; ++i_col) {
if (good_one[i_col]) {
n_for_call = F->n_this_col[i_col];
filter1d_get_robust_estimates (GMT, F, i_col, n_for_call, F->robust);
}
}
} /* That's it for the robust work */
if (F->filter_type > FILTER1D_CONVOLVE) {
/* Need to count how many good ones; use data_sum area */
Out->data = data_sum;
n_good_ones = 0;
for (i_col = 0; i_col < F->n_cols; ++i_col) {
if (i_col == F->t_col)
data_sum[i_col] = t_time;
else if (good_one[i_col]) {
data_sum[i_col] = (F->highpass) ? F->data[i_col][k] - F->this_loc[i_col] : F->this_loc[i_col];
++n_good_ones;
}
else
data_sum[i_col] = GMT->session.d_NaN;
}
if (n_good_ones) GMT_Put_Record (C, GMT_WRITE_DATA, Out);
}
else {
if (F->robust) for (i_col = 0; i_col < F->n_cols; ++i_col) F->n_this_col[i_col] = 0;
for (i_row = left; i_row < right; ++i_row) {
delta_time = t_time - F->data[F->t_col][i_row];
i_f_wt = F->half_n_f_wts + lrint (floor (0.5 + delta_time/F->dt));
if ((i_f_wt < 0) || (i_f_wt >= (int)F->n_f_wts)) continue;
for(i_col = 0; i_col < F->n_cols; ++i_col) {
if (!good_one[i_col]) continue;
if (!gmt_M_is_dnan (F->data[i_col][i_row])) {
wt = F->f_wt[i_f_wt];
val = F->data[i_col][i_row];
if (F->robust) {
med = F->this_loc[i_col];
scl = F->this_scl[i_col];
val = ((fabs(val-med)) > (2.5 * scl)) ? med : val;
}
else if (F->check_asym) { /* This wasn't already done */
if (F->data[F->t_col][i_row] < t_time) F->n_left[i_col]++;
if (F->data[F->t_col][i_row] > t_time) F->n_right[i_col]++;
}
wt_sum[i_col] += wt;
data_sum[i_col] += (wt * val);
F->n_this_col[i_col]++;
}
}
}
n_good_ones = 0;
for (i_col = 0; i_col < F->n_cols; ++i_col) {
if (!good_one[i_col]) continue;
if (!F->n_this_col[i_col]) {
good_one[i_col] = false;
continue;
}
if (F->check_asym && !(F->robust) ) {
n_l = F->n_left[i_col];
n_r = F->n_right[i_col];
diff = (n_l >= n_r) ? n_l - n_r : n_r - n_l;
symmetry = ((double)diff)/(n_l + n_r);
if (symmetry > F->sym_coeff) {
good_one[i_col] = false;
continue;
}
}
if (F->check_q && ((wt_sum[i_col] / F->n_this_col[i_col]) < F->q_factor)) {
good_one[i_col] = false;
continue;
}
++n_good_ones;
}
if (n_good_ones) {
Out->data = outval;
for (i_col = 0; i_col < F->n_cols; ++i_col) {
if (i_col == F->t_col)
outval[i_col] = t_time;
else if (good_one[i_col]) {
outval[i_col] = (F->f_operator) ? data_sum[i_col] : data_sum[i_col] / wt_sum[i_col];
if (F->highpass) outval[i_col] = F->data[i_col][k] - outval[i_col];
}
else
outval[i_col] = GMT->session.d_NaN;
}
GMT_Put_Record (C, GMT_WRITE_DATA, Out);
}
}
k++; /* Go to next output point */
}
gmt_M_free (GMT, outval);
gmt_M_free (GMT, good_one);
gmt_M_free (GMT, wt_sum);
gmt_M_free (GMT, data_sum);
gmt_M_free (GMT, Out);
gmt_M_free (GMT, F->T.array);
return (0);
}
GMT_LOCAL int filter1d_allocate_space (struct GMT_CTRL *GMT, struct FILTER1D_INFO *F) {
F->n_this_col = gmt_M_memory (GMT, NULL, F->n_cols, uint64_t);
F->data = gmt_M_memory_aligned (GMT, NULL, F->n_cols, double *);
if (F->check_asym) F->n_left = gmt_M_memory (GMT, NULL, F->n_cols, uint64_t);
if (F->check_asym) F->n_right = gmt_M_memory (GMT, NULL, F->n_cols, uint64_t);
if (F->robust || (F->filter_type > FILTER1D_CONVOLVE) ) { /* Then we need workspace */
uint64_t i;
F->work = gmt_M_memory (GMT, NULL, F->n_cols, double *);
for (i = 0; i < F->n_cols; ++i) F->work[i] = gmt_M_memory (GMT, NULL, F->n_work_alloc, double);
F->min_loc = gmt_M_memory (GMT, NULL, F->n_cols, double);
F->max_loc = gmt_M_memory (GMT, NULL, F->n_cols, double);
F->last_loc = gmt_M_memory (GMT, NULL, F->n_cols, double);
F->this_loc = gmt_M_memory (GMT, NULL, F->n_cols, double);
F->min_scl = gmt_M_memory (GMT, NULL, F->n_cols, double);
F->max_scl = gmt_M_memory (GMT, NULL, F->n_cols, double);
F->this_scl = gmt_M_memory (GMT, NULL, F->n_cols, double);
F->last_scl = gmt_M_memory (GMT, NULL, F->n_cols, double);
}
return (0);
}
GMT_LOCAL void filter1d_free_space (struct GMT_CTRL *GMT, struct FILTER1D_INFO *F) {
uint64_t i;
if (!F) return;
if (F->robust || (F->filter_type > FILTER1D_CONVOLVE) ) {
for (i = 0; i < F->n_cols; ++i) gmt_M_free (GMT, F->work[i]);
gmt_M_free (GMT, F->work);
}
for (i = 0; i < F->n_cols; ++i) gmt_M_free (GMT, F->data[i]);
gmt_M_free (GMT, F->data);
gmt_M_free (GMT, F->n_this_col);
gmt_M_free (GMT, F->n_left);
gmt_M_free (GMT, F->n_right);
gmt_M_free (GMT, F->min_loc);
gmt_M_free (GMT, F->max_loc);
gmt_M_free (GMT, F->last_loc);
gmt_M_free (GMT, F->this_loc);
gmt_M_free (GMT, F->min_scl);
gmt_M_free (GMT, F->max_scl);
gmt_M_free (GMT, F->last_scl);
gmt_M_free (GMT, F->this_scl);
gmt_M_free (GMT, F->f_wt);
}
GMT_LOCAL void filter1d_load_parameters (struct FILTER1D_INFO *F, struct FILTER1D_CTRL *Ctrl, uint64_t n_cols) {
F->filter_width = Ctrl->F.width;
F->dt = Ctrl->D.inc;
F->equidist = !Ctrl->D.active;
F->use_ends = Ctrl->E.active;
F->check_lack = Ctrl->L.active;
F->lack_width = Ctrl->L.value;
F->n_cols = n_cols;
F->t_col = Ctrl->N.col;
F->q_factor = Ctrl->Q.value;
F->check_q = Ctrl->Q.active;
F->check_asym = Ctrl->S.active;
F->sym_coeff = Ctrl->S.value;
F->t_start_t = F->t_start = Ctrl->T.T.min;
F->t_stop_t = F->t_stop = Ctrl->T.T.max;
F->t_int = Ctrl->T.T.inc;
F->T = Ctrl->T.T;
F->out_at_time = Ctrl->T.active;
F->highpass = Ctrl->F.highpass;
F->variable = Ctrl->F.variable;
}
/* Must free allocated memory before returning */
#define bailout(code) {gmt_M_free_options (mode); return (code);}
#define Return(code,...) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); GMT_Report (API, GMT_MSG_ERROR, __VA_ARGS__); bailout (code);}
#define Return2(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
EXTERN_MSC int GMT_filter1d (void *V_API, int mode, void *args) {
uint64_t col, tbl, row, seg, tseg;
int error;
unsigned int save_col[2], n_out_cols;
double last_time, new_time, in;
struct GMT_OPTION *options = NULL;
struct FILTER1D_INFO F;
struct GMT_DATASET *D = NULL;
struct GMT_DATASEGMENT *S = NULL;
struct FILTER1D_CTRL *Ctrl = NULL;
struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
/*----------------------- Standard module initialization and parsing ----------------------*/
if (API == NULL) return (GMT_NOT_A_SESSION);
if (mode == GMT_MODULE_PURPOSE) return (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 ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
/* Parse the command-line arguments */
if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, NULL, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error, "Error parsing filter1d options\n");
Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
if ((error = parse (GMT, Ctrl, options)) != 0) Return (error, "Error parsing filter1d options\n");
/*---------------------------- This is the filter1d main code ----------------------------*/
gmt_M_memset (&F, 1, struct FILTER1D_INFO); /* Init control structure to NULL */
F.n_work_alloc = GMT_CHUNK;
F.equidist = true;
GMT_Report (API, GMT_MSG_INFORMATION, "Processing input table data\n");
/* Register all data sources */
if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_LINE, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) {
Return (API->error, "Error initializing input\n");
}
/* Read the input data into memory */
if ((D = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, 0, GMT_READ_NORMAL, NULL, NULL, NULL)) == NULL) {
Return (API->error, "Error Reading input\n");
}
if (D->n_columns < 2) {
GMT_Report (API, GMT_MSG_ERROR, "Input data have %d column(s) but at least 2 are needed\n", (int)D->n_columns);
Return (GMT_DIM_TOO_SMALL, "Not enough input columns\n");
}
n_out_cols = (unsigned int)D->n_columns;
if (Ctrl->T.T.spatial) { /* Must add extra column to store distances and then compute them */
Ctrl->N.col = (int)D->n_columns;
gmt_adjust_dataset (GMT, D, D->n_columns + 1);
if (gmt_init_distaz (GMT, Ctrl->T.T.unit, Ctrl->T.T.distmode, GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE) {
Return (GMT_NOT_A_VALID_TYPE, "Bad distance unit");
}
for (tbl = 0; tbl < D->n_tables; ++tbl) { /* For each input table */
for (seg = 0; seg < D->table[tbl]->n_segments; ++seg) { /* For each segment */
S = D->table[tbl]->segment[seg];
gmt_M_free (GMT, S->data[Ctrl->N.col]);
S->data[Ctrl->N.col] = gmt_dist_array (GMT, S->data[GMT_X], S->data[GMT_Y], S->n_rows, true);
}
}
if (Ctrl->T.T.add) n_out_cols++;
}
filter1d_load_parameters (&F, Ctrl, D->n_columns); /* Pass parameters from Control structure to Filter structure */
if (GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] < F.n_cols) Return (GMT_N_COLS_VARY,
"Binary input data must have at least %d fields\n", F.n_cols);
if (strchr ("BCGMPF", Ctrl->F.filter)) { /* First deal with robustness request */
F.robust = true;
Ctrl->F.filter = (char)tolower ((int)Ctrl->F.filter);
}
switch (Ctrl->F.filter) { /* Set filter parameters */
case 'b':
F.filter_type = FILTER1D_BOXCAR;
break;
case 'c':
F.filter_type = FILTER1D_COS_ARCH;
break;
case 'g':
F.filter_type = FILTER1D_GAUSSIAN;
break;
case 'm':
F.filter_type = FILTER1D_MEDIAN;
break;
case 'p':
F.filter_type = FILTER1D_MODE;
F.mode_selection = Ctrl->F.mode;
break;
case 'l':
F.filter_type = FILTER1D_LOWER_ALL;
F.way = -1;
F.extreme = DBL_MAX;
break;
case 'L':
F.filter_type = FILTER1D_LOWER_POS;
F.way = -1;
F.kind = +1;
break;
case 'u':
F.filter_type = FILTER1D_UPPER_ALL;
F.way = +1;
F.extreme = -DBL_MAX;
break;
case 'U':
F.filter_type = FILTER1D_UPPER_NEG;
F.way = +1;
F.kind = -1;
break;
case 'f':
F.filter_type = FILTER1D_CUSTOM;
if ((error = GMT_Set_Columns (API, GMT_IN, 1, GMT_COL_FIX_NO_TEXT)) != 0) Return (error, "Error in GMT_Set_Columns");
save_col[GMT_X] = GMT->current.io.col_type[GMT_IN][GMT_X]; /* Save col type in case it is a time column */
gmt_set_column (GMT, GMT_IN, GMT_X, GMT_IS_FLOAT); /* Always read the weights as floats */
gmt_disable_bghi_opts (GMT); /* Do not want any -b -g -h -i to affect the reading from -F files */
if ((F.Fin = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_NONE, GMT_READ_NORMAL, NULL, Ctrl->F.file, NULL)) == NULL) {
Return (API->error, "Error Reading input\n");
}
gmt_reenable_bghi_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */
gmt_set_column (GMT, GMT_IN, GMT_X, save_col[GMT_X]); /* Reset this col type to whatever it actually is */
GMT_Report (API, GMT_MSG_INFORMATION, "Read %" PRIu64 " filter weights from file %s.\n", F.Fin->n_records, Ctrl->F.file);
break;
}
if (F.filter_type > FILTER1D_CONVOLVE) F.robust = false;
if (F.variable) {
if ((error = GMT_Set_Columns (API, GMT_IN, 2, GMT_COL_FIX_NO_TEXT)) != 0) Return (error, "Error in GMT_Set_Columns");
save_col[GMT_X] = GMT->current.io.col_type[GMT_IN][GMT_X]; /* Save col type in case it is a time column */
save_col[GMT_Y] = GMT->current.io.col_type[GMT_IN][GMT_Y]; /* Save col type in case it is a time column */
gmt_set_column (GMT, GMT_IN, GMT_X, GMT->current.io.col_type[GMT_IN][F.t_col]); /* Same units as time-series "t"*/
gmt_set_column (GMT, GMT_IN, GMT_Y, GMT_IS_FLOAT); /* Always read the widths as floats */
gmt_disable_bghi_opts (GMT); /* Do not want any -b -g -h -i to affect the reading from -F files */