1 #include "mrilib.h"
2 #include "mri_bport.c"
3
4 #undef BMAX
5 #define BMAX 666
6
main(int argc,char * argv[])7 int main( int argc , char *argv[] )
8 {
9 MRI_IMAGE *bpim ;
10 int nopt=1 , nbad=0 ;
11 float dt=1.0f ;
12 int dtforce=0 , ntime=0 , tt , nblock=0 , *blocklist=NULL , *blocklen=NULL ;
13 int bbot, btop ;
14 int nband=0 , nozero=0 , quad=0 ; float fbot[BMAX] , ftop[BMAX] ;
15
16 /*-----*/
17
18 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
19 printf(
20 "Usage: 1dBport [options]\n"
21 "\n"
22 "Creates a set of columns of sines and cosines for the purpose of\n"
23 "bandpassing via regression (e.g., in 3dDeconvolve). Various option\n"
24 "are given to specify the duration and structure of the time series\n"
25 "to be created. Results are written to stdout, and usually should be\n"
26 "redirected appropriately (cf. EXAMPLES, infra). The file produced\n"
27 "could be used with the '-ortvec' option to 3dDeconvolve, for example.\n"
28 "\n"
29 "OPTIONS\n"
30 "-------\n"
31 " -band fbot ftop = Specify lowest and highest frequencies in the passband.\n"
32 " fbot can be 0 if you want to do a highpass filter only;\n"
33 " on the other hand, if ftop > Nyquist frequency, then\n"
34 " it's a lowpass filter only.\n"
35 " ** This 'option' is actually mandatory! (At least once.)\n"
36 " * For the un-enlightened, the Nyquist frequency is the\n"
37 " highest frequency supported on the given grid, and\n"
38 " is equal to 0.5/TR (units are Hz if TR is in s).\n"
39 " * The lowest nonzero frequency supported on the grid\n"
40 " is equaly to 1/(N*TR), where N=number of time points.\n"
41 " ** Multiple -band options can be used, if needed.\n"
42 " If the bands overlap, regressors will NOT be duplicated.\n"
43 " * That is, '-band 0.01 0.05 -band 0.03 0.08' is the same\n"
44 " as using '-band 0.01 0.08'.\n"
45 " ** Note that if fbot==0 and ftop>=Nyquist frequency, you\n"
46 " get a 'complete' set of trig functions, meaning that\n"
47 " using these in regression is effectively a 'no-pass'\n"
48 " filter -- probably not what you want!\n"
49 " ** It is legitimate to set fbot = ftop.\n"
50 " ** The 0 frequency (fbot = 0) component is all 1, of course.\n"
51 " But unless you use the '-quad' option, nothing generated\n"
52 " herein will deal well with linear-ish or quadratic-ish\n"
53 " trends, which fall below the lowest nonzero frequency\n"
54 " representable in a full cycle on the grid:\n"
55 " f_low = 1 / ( NT * TR )\n"
56 " where NT = number of time points.\n"
57 " ** See the fourth EXAMPLE to learn how to use 3dDeconvolve\n"
58 " to generate a file of polynomials for regression fun.\n"
59 "\n"
60 " -invert = After computing which frequency indexes correspond to the\n"
61 " input band(s), invert the selection -- that is, output\n"
62 " all those frequencies NOT selected by the -band option(s).\n"
63 " See the fifth EXAMPLE.\n"
64 "\n"
65 " -nozero } Do NOT generate the 0 frequency (constant) component\n"
66 " *OR } when fbot = 0; this has the effect of setting fbot to\n"
67 " -noconst } 1/(N*TR), and is essentially a convenient way to say\n"
68 " 'eliminate all oscillations below the ftop frequency'.\n"
69 "\n"
70 " -quad = Add regressors for linear and quadratic trends.\n"
71 " (These will be the last columns in the output.)\n"
72 "\n"
73 " -input dataset } One of these options is used to specify the number of\n"
74 " *OR* } time points to be created, as in 3dDeconvolve.\n"
75 " -input1D 1Dfile } ** '-input' allow catenated datasets, as in 3dDeconvolve.\n"
76 " *OR* } ** '-input1D' assumes TR=1 unless you use the '-TR' option.\n"
77 " -nodata NT [TR] } ** One of these options is mandatory, to specify the length\n"
78 " of the time series file to generate.\n"
79 "\n"
80 " -TR del = Set the time step to 'del' rather than use the one\n"
81 " given in the input dataset (if any).\n"
82 " ** If TR is not specified by the -input dataset or by\n"
83 " -nodata or by -TR, the program will assume it is 1.0 s.\n"
84 "\n"
85 " -concat rname = As in 3dDeconvolve, used to specify the list of start\n"
86 " indexes for concatenated runs.\n"
87 " ** Also as in 3dDeconvolve, if the -input dataset is auto-\n"
88 " catenated (by providing a list of more than one dataset),\n"
89 " the run start list is automatically generated. Otherwise,\n"
90 " this option is needed if more than one run is involved.\n"
91 "\n"
92 "EXAMPLES\n"
93 "--------\n"
94 " The first example provides basis functions to filter out all frequency\n"
95 " components from 0 to 0.25 Hz:\n"
96 " 1dBport -nodata 100 1 -band 0 0.25 > highpass.1D\n"
97 "\n"
98 " The second example provides basis functions to filter out all frequency\n"
99 " components from 0.25 Hz up to the Nyquist freqency:\n"
100 " 1dBport -nodata 100 1 -band 0.25 666 > lowpass.1D\n"
101 "\n"
102 " The third example shows how to examine the results visually, for fun:\n"
103 " 1dBport -nodata 100 1 -band 0.41 0.43 | 1dplot -stdin -thick\n"
104 "\n"
105 " The fourth example shows how to use 3dDeconvolve to generate a file of\n"
106 " polynomial 'orts', in case you find yourself needing this ability someday\n"
107 " (e.g., when stranded on a desert isle, with Gilligan, the Skipper, et al.):\n"
108 " 3dDeconvolve -nodata 100 1 -polort 2 -x1D_stop -x1D stdout: | 1dcat stdin: > pol3.1D\n"
109 "\n"
110 " The fifth example shows how to use 1dBport to generate a set of regressors to\n"
111 " eliminate all frequencies EXCEPT those in the selected range:\n"
112 " 1dBport -nodata 100 1 -band 0.03 0.13 -nozero -invert | 1dplot -stdin\n"
113 " In this example, the '-nozero' flag is used because the next step will be to\n"
114 " 3dDeconvolve with '-polort 2' and '-ortvec' to get rid of the undesirable stuff.\n"
115 "\n"
116 "ETYMOLOGICAL NOTES\n"
117 "------------------\n"
118 " * The word 'ort' was coined by Andrzej Jesmanowicz, as a shorthand name for\n"
119 " a timeseries to which you want to 'orthogonalize' your data.\n"
120 " * 'Ort' actually IS an English word, and means 'a scrap of food left from a meal'.\n"
121 " As far as I know, its only usage in modern English is in crossword puzzles,\n"
122 " and in Scrabble.\n"
123 " * For other meanings of 'ort', see http://en.wikipedia.org/wiki/Ort\n"
124 " * Do not confuse 'ort' with 'Oort': http://en.wikipedia.org/wiki/Oort_cloud\n"
125 "\n"
126 "AUTHOR -- RWCox -- Jan 2012\n"
127 ) ;
128 PRINT_COMPILE_DATE ; exit(0) ;
129 }
130
131 mainENTRY("1dBport") ; machdep() ; AFNI_logger("1dBport",argc,argv) ;
132
133 /*---- scan args ----*/
134
135 while( nopt < argc ){
136
137 /*-----*/
138
139 if( strcasecmp(argv[nopt],"-nozero") == 0 ||
140 strcasecmp(argv[nopt],"-noconst") == 0 ||
141 strcasecmp(argv[nopt],"-nonzero") == 0 ){
142
143 mri_bport_set_ffbot(1) ; nopt++ ; continue ;
144 }
145
146 if( strcasecmp(argv[nopt],"-quad") == 0 ){
147 mri_bport_set_quad(1) ; quad = 1 ; nopt++ ; continue ;
148 }
149
150 /*-----*/
151
152 if( strcasecmp(argv[nopt],"-invert") == 0 ){
153 mri_bport_set_invert(1) ; nopt++ ; continue ;
154 }
155
156 /*-----*/
157
158 if( strcasecmp(argv[nopt],"-band") == 0 ){
159 float fb,ft ;
160 if( nopt+2 >= argc ) ERROR_exit("Need 2 arguments after %s",argv[nopt]) ;
161 if( nband >= BMAX ) ERROR_exit("Too many -band options given :-(") ;
162 fb = (float)strtod(argv[++nopt],NULL) ;
163 ft = (float)strtod(argv[++nopt],NULL) ;
164 if( fb < 0.0f || ft < fb )
165 ERROR_exit("Illegal values after -band: fbot=%g ftop=%g",fbot,ftop) ;
166 fbot[nband] = fb ; ftop[nband] = ft ; nband++ ;
167 nopt++ ; continue ;
168 }
169
170 /*-----*/
171
172 if( strcasecmp(argv[nopt],"-TR") == 0 || strcasecmp(argv[nopt],"-del") == 0 ){
173 if( nopt+1 >= argc ) ERROR_exit("Need 1 argument after %s",argv[nopt]) ;
174 dt = (float)strtod(argv[++nopt],NULL) ; dtforce = 1 ;
175 if( dt <= 0.0f ){
176 WARNING_message("Illegal TR=%g replaced with TR=1",dt); dt = 1.0f; dtforce = 0;
177 }
178 nopt++ ; continue ;
179 }
180
181 /*-----*/
182
183 if( strcasecmp(argv[nopt],"-concat") == 0 ){
184 MRI_IMAGE *bim ; float *bar ;
185 if( nopt+1 >= argc ) ERROR_exit("Need 2 arguments after %s",argv[nopt]) ;
186 if( nblock > 0 )
187 WARNING_message("-concat over-riding existing block list") ;
188 bim = mri_read_1D(argv[++nopt]) ;
189 if( bim == NULL ) ERROR_exit("Can't read -concat file '%s",argv[nopt]) ;
190 nblock = bim->nx ; bar = MRI_FLOAT_PTR(bim) ;
191 blocklist = (int *)malloc(sizeof(int)*nblock) ;
192 for( tt=0 ; tt < nblock ; tt++ ) blocklist[tt] = (int)floorf(bar[tt]+0.5f) ;
193 nopt++ ; continue ;
194 }
195
196 /*-----*/
197
198 if( strcasecmp(argv[nopt],"-nodata") == 0 ){
199 if( nopt+1 >= argc ) ERROR_exit("Need at least 1 argument after %s",argv[nopt]) ;
200 ntime = (int)strtod(argv[++nopt],NULL) ;
201 if( ntime < 9 ) ERROR_exit("Value after -nodata is not allowable: %d",ntime) ;
202 nopt++ ;
203 if( nopt < argc && isdigit(argv[nopt][0]) ){
204 dt = (float)strtod(argv[nopt++],NULL) ;
205 if( dt <= 0.0f ){
206 WARNING_message("Illegal TR=%g replaced with TR=1",dt) ; dt = 1.0f ;
207 }
208 }
209 continue ;
210 }
211
212 /*-----*/
213
214 if( strcasecmp(argv[nopt],"-input1D") == 0 ){
215 MRI_IMAGE *tim ;
216 if( nopt+1 >= argc ) ERROR_exit("Need 1 argument after %s",argv[nopt]) ;
217 tim = mri_read_1D(argv[++nopt]) ;
218 if( tim == NULL ) ERROR_exit("Can't read -input1D file '%s'",argv[nopt]) ;
219 ntime = tim->nx ; mri_free(tim) ;
220 if( ntime < 9 ) ERROR_exit("Length of -input1D file is not allowable: %d",ntime) ;
221 nopt++ ; continue ;
222 }
223
224 if( strcasecmp(argv[nopt],"-input") == 0 ){
225 int iopt , slen=0 ; char *fname ; THD_3dim_dataset *dset ;
226 nopt++;
227 if( nopt >= argc ) ERROR_exit("Need 1 or more arguments after -input") ;
228 for( iopt=nopt ; iopt < argc && argv[iopt][0] != '-' ; iopt++ )
229 slen += strlen(argv[iopt])+8 ;
230 fname = calloc(sizeof(char),slen) ;
231 for( iopt=nopt ; iopt < argc && argv[iopt][0] != '-' ; iopt++ ){
232 strcat( fname , argv[iopt] ) ;
233 strcat( fname , " " ) ;
234 }
235 slen = strlen(fname) ; fname[slen-1] = '\0' ; /* trim last blank */
236 nopt = iopt ;
237 dset = THD_open_dataset(fname) ;
238 CHECK_OPEN_ERROR(dset,fname) ;
239 if( !dtforce ){ DSET_UNMSEC(dset) ; dt = DSET_TR(dset) ; }
240 ntime = DSET_NVALS(dset) ;
241 if( DSET_IS_TCAT(dset) ){
242 nblock = dset->tcat_num ;
243 blocklist = (int *)malloc(sizeof(int)*nblock) ;
244 blocklist[0] = 0 ;
245 for( tt=0 ; tt < nblock-1 ; tt++ )
246 blocklist[tt+1] = blocklist[tt] + dset->tcat_len[tt] ;
247 }
248 DSET_delete(dset) ; continue ;
249 }
250
251 /*-----*/
252
253 ERROR_message("Unknown option: %s",argv[nopt]) ;
254 suggest_best_prog_option(argv[0], argv[nopt]) ;
255 exit(1) ;
256
257 } /* end of scanning options */
258
259 /*----- check inputs -----*/
260
261 if( ntime <= 0 ){
262 ERROR_message("No way to determine time series length!") ; nbad++ ;
263 }
264 if( nband <= 0 ){
265 ERROR_message("No -band option given! What bandpass do you want?") ; nbad++ ;
266 }
267
268 if( nblock == 0 ){
269 nblock = 1 ;
270 blocklist = (int *)malloc(sizeof(int)) ;
271 blocklist[0] = 0 ;
272 }
273
274 blocklen = (int *)malloc(sizeof(int)*nblock) ;
275 for( tt=0 ; tt < nblock ; tt++ ){
276 bbot = blocklist[tt] ;
277 btop = (tt+1 < nblock) ? blocklist[tt+1] : ntime ;
278 blocklen[tt] = btop - bbot ;
279 if( blocklen[tt] < 9 ){
280 ERROR_message("Block length #%d = %d ==> too small!",tt+1,blocklen[tt]) ; nbad++ ;
281 }
282 }
283
284 if( nbad > 0 )
285 ERROR_exit("Cannot continue after above error%s :-(", (nbad==1) ? "\0" : "s" ) ;
286
287 fprintf(stderr,"++ Block length%s:",(nblock==1)?"\0":"s") ;
288 for( tt=0 ; tt < nblock ; tt++ ) fprintf(stderr," %d",blocklen[tt]) ;
289 fprintf(stderr,"\n") ;
290
291 /*----- work -----*/
292
293 bpim = mri_bport_multi( dt,nband,fbot,ftop , nblock,blocklen ) ;
294 if( bpim == NULL )
295 ERROR_exit("Computation fails! :-(((") ;
296 mri_write_1D( "-" , bpim ) ;
297 exit(0) ;
298 }
299