1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*******************************************************
8 * 3dFourier *
9 * T. Ross and K. Heimerl 7/99 *
10 *******************************************************/
11
12 #include "mrilib.h"
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <string.h>
16 #include <math.h>
17
18 #if !defined(TRUE)
19 #define TRUE (1==1)
20 #define FALSE !TRUE
21 #endif
22
23 static void *My_Malloc(size_t);
24
Error_Exit(char * message)25 void Error_Exit(char *message)
26 {
27 fprintf(stderr, "\n\nError in 3dFourier:\n%s\n\nTry 3dFourier -help\n",message);
28 exit(1);
29 }
30
31 static int g_argc ; /* Jul 2010 */
32 static char **g_argv ;
33 #define ARGC g_argc
34 #define ARGV g_argv
35
36 /* include the filter and filter driver */
37 #include "fourier_filter.c"
38
39 /***********************************************************************************/
40
help_message(void)41 void help_message(void)
42 {
43 printf(
44 "3dFourier \n"
45 "(c) 1999 Medical College of Wisconsin\n"
46 "by T. Ross and K. Heimerl\n"
47 "Version 0.8 last modified 8-17-99\n\n"
48 "Usage: 3dFourier [options] dataset\n\n"
49 "* Program to lowpass and/or highpass each voxel time series in a\n"
50 " dataset, via the FFT.\n"
51 "* Also see program 3dBandpass for something similar.\n\n"
52 "The parameters and options are:\n"
53 " dataset an afni compatible 3d+time dataset to be operated upon\n"
54 " -prefix name output name for new 3d+time dataset [default = fourier]\n"
55 " -lowpass f low pass filter with a cutoff of f Hz\n"
56 " -highpass f high pass filter with a cutoff of f Hz\n"
57 " -ignore n ignore the first n images [default = 1]\n"
58 /* " -autocorr compute the autocorrelation of the data\n" */
59 " -retrend Any mean and linear trend are removed before filtering.\n"
60 " This will restore the trend after filtering.\n"
61 "\nNote that by combining the lowpass and highpass options, one can construct\n"
62 "bandpass and notch filters\n"
63 );
64
65 printf("\n" MASTER_SHORTHELP_STRING ) ;
66 PRINT_COMPILE_DATE ;
67 }
68
69 /*--------------------------------------------------------------------------*/
70
main(int argc,char * argv[])71 int main(int argc, char *argv[])
72 {
73 int j, narg=1, autocorr=FALSE, retrend=FALSE, ignore=1;
74 char *prefix=NULL, *err;
75 float low_fc = 0.0, high_fc = 0.0;
76 THD_3dim_dataset *input=NULL;
77
78 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){ help_message(); exit(0); }
79
80 g_argc = argc ; g_argv = argv ; /* Jul 2010: save global copies */
81
82 machdep(); PRINT_VERSION("3dFourier"); AFNI_logger("3dFourier",argc,argv);
83
84 /* Loop over arguements and pull out what we need */
85 while( narg < argc && argv[narg][0] == '-' ){
86
87 if( strncmp(argv[narg],"-prefix",5) == 0 ) {
88 narg++;
89 if (narg==argc)
90 Error_Exit("-prefix must be followed by a file name");
91 j = strlen(argv[narg]);
92 prefix =(char *)My_Malloc((j+1)*sizeof(char));
93 strcpy(prefix, argv[narg++]);
94 continue;
95 }
96
97 if (strncmp(argv[narg], "-lowpass", 5) == 0) {
98 narg++;
99 if (narg==argc)
100 Error_Exit("-lowpass must be followed by a frequency");
101 low_fc = (float)atof(argv[narg++]);
102 continue;
103 }
104
105 if (strncmp(argv[narg], "-highpass", 5) == 0) {
106 narg++;
107 if (narg==argc)
108 Error_Exit("-highpass must be followed by a frequency");
109 high_fc = (float)atof(argv[narg++]);
110 continue;
111 }
112
113 if (strncmp(argv[narg], "-ignore", 5) == 0) {
114 narg++;
115 if (narg==argc)
116 Error_Exit("-ignore must be followed by an integer");
117 ignore = (int)atol(argv[narg++]);
118 if ((ignore<0) || (ignore>10))
119 Error_Exit("-ignore must be between 0 and 10, inclusive");
120 continue;
121 }
122 /*
123 if (strncmp(argv[narg], "-autocorr", 5) == 0) {
124 narg++;
125 autocorr = TRUE;
126 continue;
127 }
128 */
129 if (strncmp(argv[narg], "-retrend", 5) == 0) {
130 narg++;
131 retrend = TRUE;
132 continue;
133 }
134
135 if (strncmp(argv[narg], "-help", 5) == 0) {
136 help_message();
137 }
138
139
140 Error_Exit("Illegal or unrecoginized option");
141
142 } /* end of while over arguements */
143
144 #if 1
145 if( low_fc > 0.0f && high_fc > 0.0f && low_fc <= high_fc )
146 fprintf(stderr,"** WARNING: lowpass=%f is <= than highpass=%f\n"
147 "** -------: results from 3dFourier are suspect!\n" ,
148 low_fc,high_fc) ;
149 #endif
150
151 if( narg >= argc )
152 Error_Exit("No input datasets!?\n");
153
154 if(prefix==NULL) {
155 prefix=(char *)My_Malloc(8*sizeof(char));
156 strcpy(prefix,"fourier");
157 }
158
159 /* try to open input dataset */
160
161 input = THD_open_dataset( argv[narg] ) ; /* 16 Sep 1999 */
162 CHECK_OPEN_ERROR(input,argv[narg]) ;
163
164 if (!DSET_GRAPHABLE(input))
165 Error_Exit("Input dataset is not a 3d+time dataset");
166
167 if (DSET_BRICK_TYPE(input, 0)==MRI_complex)
168 Error_Exit("Sorry, I can't deal with complex 3d+time datasets");
169
170 err=Fourier_Filter_Driver(input, low_fc, high_fc, ignore, autocorr, retrend, prefix);
171 if (err!=NULL)
172 Error_Exit(err);
173
174 return 0; /* All went well exit code */
175 }
176