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 #include "mrilib.h"
8
9 /*--------------------------------------------------------------------------*/
10
11 static void TS_copy_input_to_output(void) ; /* prototype */
12
TS_syntax(char * str)13 void TS_syntax(char * str)
14 {
15 if( str != NULL ){
16 ERROR_exit("%s",str) ; exit(1) ;
17 }
18
19 printf(
20 "Usage: 3dTshift [options] dataset\n\n"
21 "* Shifts voxel time series from the input dataset so that the separate\n"
22 " slices are aligned to the same temporal origin. By default, uses the\n"
23 " slicewise shifting information in the dataset header (from the 'tpattern'\n"
24 " input to program to3d).\n"
25 "\n"
26 "Method: detrend -> interpolate -> retrend (optionally)\n"
27 "\n"
28 "* The input dataset can have a sub-brick selector attached, as documented\n"
29 " in '3dcalc -help'.\n"
30 "\n"
31 "* The output dataset time series will be interpolated from the input to\n"
32 " the new temporal grid. This may not be the best way to analyze your\n"
33 " data, but it can be convenient.\n"
34 "\n"
35 "* Slices where significant time interpolation happens will have extra\n"
36 " temporal autocorrelation introduced by the interpolation. The amount\n"
37 " of extra correlation along the time axis depends on the type of\n"
38 " interpolation used. Higher order interpolation will produce smaller\n"
39 " such 'extra' correlation; in order, from lowest (most extra correlation)\n"
40 " to highest (least extra correlation):\n"
41 " -linear -cubic -quintic -heptic\n"
42 " -wsinc5 -wsinc9 -Fourier\n"
43 "* The last two methods do not add much correlation in time. However, they\n"
44 " have the widest interpolation 'footprint' and so the output data values\n"
45 " will have contributions from data points further away in time.\n"
46 "* To properly account for these extra correlations, which vary in space,\n"
47 " we advise you to analyze the time series using 3dREMLfit, which uses\n"
48 " a voxel-dependent prewhitening (de-correlating) linear regression method,\n"
49 " unlike most other FMRI time series regression software.\n"
50 " ++ Or else use '-wsinc9' interpolation, which has a footprint of 18 time points:\n"
51 " 9 before and 9 after the intermediate time at which the value is output.\n"
52 "\n"
53 "WARNINGS:\n"
54 "--------\n"
55 "* Please recall the phenomenon of 'aliasing': frequencies above 1/(2*TR) can't\n"
56 " be properly interpolated. For most 3D FMRI data, this means that cardiac\n"
57 " and respiratory effects will not be treated properly by this program.\n"
58 "\n"
59 "* The images at the beginning of a high-speed FMRI imaging run are usually\n"
60 " of a different quality than the later images, due to transient effects\n"
61 " before the longitudinal magnetization settles into a steady-state value.\n"
62 " These images should not be included in the interpolation! For example,\n"
63 " if you wish to exclude the first 4 images, then the input dataset should\n"
64 " be specified in the form 'prefix+orig[4..$]'. Alternatively, you can\n"
65 " use the '-ignore ii' option.\n"
66 "\n"
67 "* It seems to be best to use 3dTshift before using 3dvolreg.\n"
68 " (But this statement is controversial.)\n"
69 "\n"
70 "* If the input dataset does not have any slice timing information, and\n"
71 " '-tpattern' is not given, then this program just copies the input to\n"
72 " the output. [02 Nov 2011 -- formerly, it failed]\n"
73 "\n"
74 "* Please consider the potential impact of 3dTshift on any subsequent\n"
75 " linear regression model. While the temporal resampling of 3dTshift is\n"
76 " not exact, it is attempting to interpolate the slice timing so that it\n"
77 " is as if each volume were acquired at time 'tzero' + k*TR. So with\n"
78 " -tzero 0, it becomes akin to each entire volume being acquired at the\n"
79 " very beginning of its TR. By default, the offset is the average offset\n"
80 " across the slices, which for alt+z or seq is:\n"
81 " (nslices-1)/nslices * TR/2\n"
82 " That average approaches TR/2 as the number of slices increases.\n"
83 "\n"
84 " The new slice/volume timing is intended to be the real timing from the\n"
85 " start of the run.\n"
86 "\n"
87 " How might this affect stimulus timing in 3dDeconvolve?\n"
88 " 3dDeconvolve creates regressors based on volume times of k*TR, matching\n"
89 " tzero=0. So an event at run time t=0 would start at the time of volume\n"
90 " #0. However using -tzero 1 (or the default, in the case of TR~=2s),\n"
91 " an event at run time t=0 would then be 1s *before* the first volume.\n"
92 " Note that this matches reality. An event at time t=0 happens before\n"
93 " all but the first acquired slice. In particular, a slice acquired at\n"
94 " TR offset 1s might be unaffected by 3dTshift. And an event at run time\n"
95 " t=0 seems to happen at time t=-1s from the perspective of that slice.\n"
96 "\n"
97 " To align stimulus times with the applied tzero of 3dTshift, tzero\n"
98 " should be subtracted from each stimulus event time (3dDeconvolve\n"
99 " effectively subtracts tzero from the EPI timing, so that should be\n"
100 " applied to the event times as well).\n"
101 "\n"
102 "\n"
103 "OPTIONS:\n"
104 "-------\n"
105 " -verbose = print lots of messages while program runs\n"
106 "\n"
107 " -TR ddd = use 'ddd' as the TR, rather than the value\n"
108 " stored in the dataset header using to3d.\n"
109 " You may attach the suffix 's' for seconds,\n"
110 " or 'ms' for milliseconds.\n"
111 "\n"
112 " -tzero zzz = align each slice to time offset 'zzz';\n"
113 " the value of 'zzz' must be between the\n"
114 " minimum and maximum slice temporal offsets.\n"
115 " N.B.: The default alignment time is the average\n"
116 " of the 'tpattern' values (either from the\n"
117 " dataset header or from the -tpattern option)\n"
118 "\n"
119 " -slice nnn = align each slice to the time offset of slice\n"
120 " number 'nnn' - only one of the -tzero and\n"
121 " -slice options can be used.\n"
122 "\n"
123 " -prefix ppp = use 'ppp' for the prefix of the output file;\n"
124 " the default is 'tshift'.\n"
125 "\n"
126 " -ignore ii = Ignore the first 'ii' points. (Default is ii=0.)\n"
127 " The first ii values will be unchanged in the output\n"
128 " (regardless of the -rlt option). They also will\n"
129 " not be used in the detrending or time shifting.\n"
130 "\n"
131 " -rlt = Before shifting, the mean and linear trend\n"
132 " -rlt+ = of each time series is removed. The default\n"
133 " action is to add these back in after shifting.\n"
134 " -rlt means to leave both of these out of the output\n"
135 " -rlt+ means to add only the mean back into the output\n"
136 " (cf. '3dTcat -help')\n"
137 "\n"
138 " -no_detrend = Do not remove or restore linear trend.\n"
139 " Heptic becomes the default interpolation method.\n"
140 "\n"
141 " ** Options to choose the temporal interpolation method: **\n"
142 " -Fourier = Use a Fourier method (the default: most accurate; slowest).\n"
143 " -linear = Use linear (1st order polynomial) interpolation (least accurate).\n"
144 " -cubic = Use the cubic (3rd order) Lagrange polynomial interpolation.\n"
145 " -quintic = Use the quintic (5th order) Lagrange polynomial interpolation.\n"
146 " -heptic = Use the heptic (7th order) Lagrange polynomial interpolation.\n"
147 " -wsinc5 = Use weighted sinc interpolation - plus/minus 5 [Aug 2019].\n"
148 " -wsinc9 = Use weighted sinc interpolation - plus/minus 9.\n"
149 "\n"
150 " -tpattern ttt = use 'ttt' as the slice time pattern, rather\n"
151 " than the pattern in the input dataset header;\n"
152 " 'ttt' can have any of the values that would\n"
153 " go in the 'tpattern' input to to3d, described below:\n"
154 "\n"
155 " alt+z = altplus = alternating in the plus direction\n"
156 " alt+z2 = alternating, starting at slice #1 instead of #0\n"
157 " alt-z = altminus = alternating in the minus direction\n"
158 " alt-z2 = alternating, starting at slice #nz-2 instead of #nz-1\n"
159 " seq+z = seqplus = sequential in the plus direction\n"
160 " seq-z = seqminus = sequential in the minus direction\n"
161 " @filename = read temporal offsets from 'filename'\n"
162 "\n"
163 " For example if nz = 5 and TR = 1000, then the inter-slice\n"
164 " time is taken to be dt = TR/nz = 200. In this case, the\n"
165 " slices are offset in time by the following amounts:\n"
166 "\n"
167 " S L I C E N U M B E R\n"
168 " tpattern 0 1 2 3 4 Comment\n"
169 " --------- --- --- --- --- --- -------------------------------\n"
170 " altplus 0 600 200 800 400 Alternating in the +z direction\n"
171 " alt+z2 400 0 600 200 800 Alternating, but starting at #1\n"
172 " altminus 400 800 200 600 0 Alternating in the -z direction\n"
173 " alt-z2 800 200 600 0 400 Alternating, starting at #nz-2 \n"
174 " seqplus 0 200 400 600 800 Sequential in the +z direction\n"
175 " seqminus 800 600 400 200 0 Sequential in the -z direction\n"
176 "\n"
177 " If @filename is used for tpattern, then nz ASCII-formatted numbers\n"
178 " are read from the file. These indicate the time offsets for each\n"
179 " slice. For example, if 'filename' contains\n"
180 " 0 600 200 800 400\n"
181 " then this is equivalent to 'altplus' in the above example.\n"
182 " (nz = number of slices in the input dataset)\n"
183 "\n"
184 "N.B.: if you are using -tpattern, make sure that the units supplied\n"
185 " match the units of TR in the dataset header, or provide a\n"
186 " new TR using the -TR option.\n"
187 "\n"
188 "As a test of how well 3dTshift interpolates, you can take a dataset\n"
189 "that was created with '-tpattern alt+z', run 3dTshift on it, and\n"
190 "then run 3dTshift on the new dataset with '-tpattern alt-z' -- the\n"
191 "effect will be to reshift the dataset back to the original time\n"
192 "grid. Comparing the original dataset to the shifted-then-reshifted\n"
193 "output will show where 3dTshift does a good job and where it does\n"
194 "a bad job.\n"
195 "\n"
196 "******* Voxel-Wise Shifting -- New Option [Sep 2011] *******\n"
197 "\n"
198 " -voxshift fset = Read in dataset 'fset' and use the values in there\n"
199 " to shift each input dataset's voxel's time series a\n"
200 " different amount. The values in 'fset' are NOT in\n"
201 " units of time, but rather are fractions of a TR\n"
202 " to shift -- a positive value means to shift backwards.\n"
203 " * To compute an fset-style dataset that matches the\n"
204 " time pattern of an existing dataset, try\n"
205 " set TR = 2.5\n"
206 " 3dcalc -a 'dset+orig[0..1]' -datum float -prefix Toff -expr \"t/${TR}-l\"\n"
207 " where you first set the shell variable TR to the true TR\n"
208 " of the dataset, then create a dataset Toff+orig with the\n"
209 " fractional shift of each slice stored in each voxel. Then\n"
210 " the two commands below should give identical outputs:\n"
211 " 3dTshift -ignore 2 -tzero 0 -prefix Dold -heptic dset+orig\n"
212 " 3dTshift -ignore 2 -voxshift Toff+orig -prefix Dnew -heptic dset+orig\n"
213 "\n"
214 " Use of '-voxshift' means that options such as '-tzero' and '-tpattern' are\n"
215 " ignored -- the burden is on you to encode all the shifts into the 'fset'\n"
216 " dataset somehow. (3dcalc can be your friend here.)\n"
217 "\n"
218 "-- RWCox - 31 October 1999, et cetera\n"
219 ) ;
220
221 printf("\n" MASTER_SHORTHELP_STRING ) ;
222
223 PRINT_COMPILE_DATE ; exit(0) ;
224 }
225
226 /*-----------------------------------------------------------------------------*/
227
228 static float TS_TR = 0.0 ;
229 static int TS_tunits = UNITS_SEC_TYPE ;
230 static float *TS_tpat = NULL ;
231 static float TS_tzero = -1.0 ;
232 static int TS_slice = -1 ;
233 static int TS_rlt = 0 ; /* 0=add both in; 1=add neither; 2=add mean */
234
235 static int TS_detrend = 1 ; /* do any detrend? 3 Jan 2007 [rickr] */
236 static int TS_verbose = 0 ;
237
238 static int TS_ignore = 0 ; /* 15 Feb 2001 */
239
240 static THD_3dim_dataset *TS_dset = NULL , *TS_oset = NULL ;
241
242 static THD_3dim_dataset *TS_fset = NULL ; /* 21 Sep 2011 */
243
244 static char *TS_tpattern = NULL ;
245
246 static char TS_prefix[THD_MAX_NAME] = "tshift" ;
247
248 /*-----------------------------------------------------------------------------*/
249
main(int argc,char * argv[])250 int main( int argc , char *argv[] )
251 {
252 int nopt=1 ;
253 int nzz, ii,jj,kk , ntt,nxx,nyy,nxy , nup ;
254 float tomax,tomin , fmin,fmax , gmin=0.0,gmax=0.0 , f0,f1 , g0,g1 ;
255 float ffmin=0.0,ffmax=0.0 , ggmin=0.0,ggmax=0.0 , fshift=0.0,gshift=0.0 ;
256 MRI_IMAGE *flim , *glim=NULL ;
257 float *far , *gar ;
258 int ignore=0 , BAD=0 ;
259
260 /*- scan command line -*/
261
262 if( argc < 2 || strcmp(argv[1],"-help") == 0 ) TS_syntax(NULL) ;
263
264 mainENTRY("3dTshift main"); machdep(); AFNI_logger("3dTshift",argc,argv);
265 PRINT_VERSION("3dTshift");
266
267 SHIFT_set_method( MRI_FOURIER ) ;
268
269 while( nopt < argc && argv[nopt][0] == '-' ){
270
271 if( strcmp(argv[nopt],"-voxshift") == 0 ){ /* 21 Sep 2011 */
272 if( TS_fset != NULL ) ERROR_exit("Can't use -voxshift twice!") ;
273 TS_fset = THD_open_dataset( argv[++nopt] ) ;
274 if( TS_fset == NULL )
275 ERROR_exit("Can't open -voxshift dataset '%s'",argv[nopt]) ;
276 DSET_load(TS_fset) ; CHECK_LOAD_ERROR(TS_fset) ;
277 nopt++ ; continue ;
278 }
279
280 if( strcmp(argv[nopt],"-BAD") == 0 ){
281 BAD = 1 ; nopt++ ; continue ;
282 }
283
284 if( strncmp(argv[nopt],"-verbose",5) == 0 ){
285 TS_verbose++ ;
286 nopt++ ; continue ;
287 }
288
289 if( strcmp(argv[nopt],"-ignore") == 0 ){ /* 15 Feb 2001 */
290 TS_ignore = (int) strtod(argv[++nopt],NULL) ;
291 if( TS_ignore < 0 ) ERROR_exit("-ignore value %d is negative!",TS_ignore) ;
292 nopt++ ; continue ;
293 }
294
295 if( strncmp(argv[nopt],"-Fourier",4) == 0 || strncmp(argv[nopt],"-fourier",4) == 0 ){
296 SHIFT_set_method( MRI_FOURIER ) ;
297 nopt++ ; continue ;
298 }
299
300 if( strncmp(argv[nopt],"-cubic",4) == 0 || strncmp(argv[nopt],"-Cubic",4) == 0 ){
301 SHIFT_set_method( MRI_CUBIC ) ;
302 nopt++ ; continue ;
303 }
304
305 if( strncmp(argv[nopt],"-quintic",4) == 0 || strncmp(argv[nopt],"-Quintic",4) == 0 ){
306 SHIFT_set_method( MRI_QUINTIC ) ;
307 nopt++ ; continue ;
308 }
309
310 if( strncmp(argv[nopt],"-heptic",4) == 0 || strncmp(argv[nopt],"-Heptic",4) == 0 ){
311 SHIFT_set_method( MRI_HEPTIC ) ;
312 nopt++ ; continue ;
313 }
314
315 if( strncasecmp(argv[nopt],"-wsinc5",7) == 0 ){ /* Aug 2019 */
316 SHIFT_set_method( MRI_WSINC5 ) ;
317 nopt++ ; continue ;
318 }
319
320 if( strncasecmp(argv[nopt],"-wsinc9",7) == 0 ){ /* Aug 2019 */
321 SHIFT_set_method( MRI_WSINC9 ) ;
322 nopt++ ; continue ;
323 }
324
325 if( strncmp(argv[nopt],"-linear",4) == 0 || strncmp(argv[nopt],"-Linear",4) == 0 ){
326 SHIFT_set_method( MRI_LINEAR ) ;
327 nopt++ ; continue ;
328 }
329
330 if( strcmp(argv[nopt],"-TR") == 0 ){
331 char * eptr ;
332 if( ++nopt >= argc ) ERROR_exit("-TR needs an argument!") ;
333 TS_TR = strtod( argv[nopt] , &eptr ) ;
334 if( TS_TR <= 0.0 ) ERROR_exit("illegal value '%s' after -TR!",argv[nopt]) ;
335
336 if( strcasecmp(eptr,"ms")==0 || strcasecmp(eptr,"msec")==0 ){
337 TS_tunits = UNITS_MSEC_TYPE ;
338 WARNING_message("TR expressed in milliseconds is deprecated [not wanted].") ;
339 } else if( strcasecmp(eptr,"s")==0 || strcasecmp(eptr,"sec")==0 ){
340 TS_tunits = UNITS_SEC_TYPE ;
341 } else if( strcasecmp(eptr,"Hz")==0 || strcasecmp(eptr,"Hertz")==0 ){
342 TS_tunits = UNITS_HZ_TYPE ;
343 }
344
345 nopt++ ; continue ;
346 }
347
348 if( strcmp(argv[nopt],"-tzero") == 0 ){
349 if( ++nopt >= argc ) ERROR_exit("-tzero needs an argument!") ;
350 TS_tzero = strtod( argv[nopt] , NULL ) ;
351 if( TS_tzero < 0.0 ) ERROR_exit("illegal value '%s' after -tzero!",argv[nopt]) ;
352 nopt++ ; continue ;
353 }
354
355 if( strcmp(argv[nopt],"-slice") == 0 ){
356 if( ++nopt >= argc ) ERROR_exit("-slice needs an argument!") ;
357 TS_slice = strtod( argv[nopt] , NULL ) ;
358 if( TS_slice < 0 ) ERROR_exit("illegal value '%s' after -slice!",argv[nopt]) ;
359 nopt++ ; continue ;
360 }
361
362 if( strcmp(argv[nopt],"-prefix") == 0 ){
363 if( ++nopt >= argc ) ERROR_exit("-prefix needs an argument!") ;
364 MCW_strncpy( TS_prefix , argv[nopt] , THD_MAX_NAME ) ;
365 if( !THD_filename_ok(TS_prefix) )
366 ERROR_exit("illegal value '%s' after -prefix",argv[nopt]) ;
367 nopt++ ; continue ;
368 }
369
370 if( strcmp(argv[nopt],"-rlt") == 0 ){
371 if( !TS_detrend ) ERROR_exit("cannot use both -rlt and -no_detrend");
372 TS_rlt = 1 ;
373 nopt++ ; continue ;
374 }
375
376 if( strcmp(argv[nopt],"-rlt+") == 0 ){
377 if( !TS_detrend ) ERROR_exit("cannot use both -rlt+ and -no_detrend");
378 TS_rlt = 2 ;
379 nopt++ ; continue ;
380 }
381
382 if( strcmp(argv[nopt],"-no_detrend") == 0 ){ /* 3 Jan 2007 */
383 if( TS_rlt ) ERROR_exit("cannot use both -rlt and -no_detrend");
384 if( SHIFT_get_method() == MRI_FOURIER ){
385 ERROR_exit("found -no_detrend, changing default to -heptic\n");
386 SHIFT_set_method(MRI_HEPTIC);
387 }
388 TS_detrend = 0 ;
389 TS_rlt = 2 ; /* still de-mean/re-mean, as Bob suggests */
390 nopt++ ; continue ;
391 }
392
393 if( strcmp(argv[nopt],"-tpattern") == 0 ){
394 if( ++nopt >= argc ) ERROR_exit("-tpattern needs an argument!") ;
395 TS_tpattern = argv[nopt] ;
396 nopt++ ; continue ;
397 }
398
399 ERROR_exit("Unknown option: %s",argv[nopt]) ;
400
401 } /* end of scan command line */
402
403 if( TS_detrend == 0 && SHIFT_get_method() == MRI_FOURIER )
404 WARNING_message("-no_detrend with Fourier interpolation is dangerous");
405
406 /*- open dataset; extract values, check for errors -*/
407
408 if( nopt >= argc ) ERROR_exit("Need a dataset input?!") ;
409 if( TS_verbose ) printf("++ opening input dataset header\n") ;
410 TS_dset = THD_open_dataset( argv[nopt] ) ;
411 if( TS_dset == NULL ) ERROR_exit("Can't open input dataset '%s'",argv[nopt]) ;
412
413 nxx = DSET_NX(TS_dset) ; /* get dimensions */
414 nyy = DSET_NY(TS_dset) ; nxy = nxx * nyy ;
415 nzz = DSET_NZ(TS_dset) ;
416 ntt = DSET_NVALS(TS_dset) ;
417
418 if( DSET_NVALS(TS_dset) < 2 ){
419 WARNING_message("Input dataset has only 1 value per voxel!") ;
420 TS_copy_input_to_output() ;
421 }
422
423 if( TS_slice >= nzz ) ERROR_exit("-slice value is too large (%d >= %d)",TS_slice,nzz) ;
424
425 if( TS_ignore > ntt-5 ) ERROR_exit("-ignore value %d is too large",TS_ignore) ;
426
427 if( TS_TR <= 0.0 ){ /* set TR from dataset */
428 if( TS_dset->taxis != NULL ){
429 TS_TR = DSET_TIMESTEP(TS_dset) ;
430 TS_tunits = TS_dset->taxis->units_type ;
431 }
432 if( TS_TR <= 0.0 ){
433 TS_TR = 1.0f ;
434 TS_tunits = UNITS_SEC_TYPE ;
435 }
436 if( TS_verbose )
437 printf("++ using dataset TR = %g %s\n",TS_TR,UNITS_TYPE_LABEL(TS_tunits)) ;
438 }
439
440 if( TS_fset != NULL ){
441
442 if( nxx != DSET_NX(TS_fset) ||
443 nyy != DSET_NY(TS_fset) ||
444 nzz != DSET_NZ(TS_fset) )
445 ERROR_exit("-voxshift and input datasets don't match!") ;
446
447 } else {
448
449 if( TS_dset->taxis == NULL ){
450 if( TS_TR == 0.0 || TS_tpattern == NULL ){
451 WARNING_message("dataset has no time axis!") ;
452 TS_copy_input_to_output() ;
453 }
454 } else if( TS_tpattern == NULL && TS_dset->taxis->toff_sl == NULL ){
455 WARNING_message("dataset is already aligned in time!") ;
456 TS_copy_input_to_output() ;
457 }
458
459 if( TS_tpattern != NULL ){ /* set pattern */
460 TS_tpat = TS_parse_tpattern( nzz , TS_TR , TS_tpattern ) ;
461 /* move exit() out of TS_parse_tpattern 8 Mar 2013 [rickr] */
462 if ( ! TS_tpat ) exit(1) ;
463 } else {
464 if( TS_dset->taxis->nsl != nzz ){
465 WARNING_message("dataset temporal pattern is malformed!") ; /* should not happen */
466 TS_copy_input_to_output() ;
467 }
468
469 TS_tpat = (float *) malloc( sizeof(float) * nzz ) ;
470 memcpy( TS_tpat , TS_dset->taxis->toff_sl , sizeof(float)*nzz ) ;
471 }
472 if( TS_verbose ){
473 printf("++ using tpattern = ") ;
474 for( ii=0 ; ii < nzz ; ii++ ) printf("%g ",TS_tpat[ii]) ;
475 printf("%s\n",UNITS_TYPE_LABEL(TS_tunits)) ;
476 }
477
478 tomin = WAY_BIG ; tomax = -WAY_BIG ; /* check pattern */
479 for( ii=0 ; ii < nzz ; ii++ ){
480 if( TS_tpat[ii] > tomax ) tomax = TS_tpat[ii] ;
481 if( TS_tpat[ii] < tomin ) tomin = TS_tpat[ii] ;
482 }
483 if( tomin < 0.0 || tomax > TS_TR ){
484 WARNING_message("some value in tpattern is outside range 0..TR=%g",TS_TR) ;
485 TS_copy_input_to_output() ;
486 } else if( tomin >= tomax ){
487 WARNING_message("temporal pattern is already aligned in time!") ;
488 TS_copy_input_to_output() ;
489 }
490
491
492 if( TS_slice >= 0 && TS_slice < nzz ){ /* set common time point */
493 TS_tzero = TS_tpat[TS_slice] ;
494 } else if( TS_tzero < 0.0 ){
495 TS_tzero = 0.0 ;
496 for( ii=0 ; ii < nzz ; ii++ ) TS_tzero += TS_tpat[ii] ;
497 TS_tzero /= nzz ;
498 }
499 if( TS_verbose ) printf("++ common time point set to %g\n",TS_tzero) ;
500
501 }
502
503 /*- copy input dataset, modify it to be the output -*/
504
505 if( TS_verbose ) printf("++ copying input dataset bricks\n") ;
506
507 TS_oset = EDIT_full_copy( TS_dset , TS_prefix ) ;
508 if( TS_oset == NULL )
509 ERROR_exit("Can't copy input dataset '%s'",DSET_HEADNAME(TS_dset)) ;
510 DSET_unload( TS_dset ) ;
511
512 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(TS_oset)) )
513 ERROR_exit("output dataset already exists '%s'",DSET_HEADNAME(TS_oset)) ;
514
515 tross_Copy_History( TS_dset , TS_oset ) ;
516 tross_Make_History( "3dTshift" , argc,argv , TS_oset ) ;
517
518 /*- reconfigure the time axis -*/
519
520 printf("++ updating time offset to %g\n",TS_tzero) ;
521 EDIT_dset_items( TS_oset ,
522 ADN_ntt , ntt , /* in case not already set */
523 ADN_ttdel , TS_TR , /* may have changed */
524 ADN_tunits , TS_tunits , /* may have changed */
525 ADN_nsl , 0 , /* will have no offsets when done */
526 /* set toffset to whatever is applied
527 * see: https://github.com/afni/afni/issues/297
528 * 21 Sep 2021 [rickr] */
529 ADN_ttorg , TS_tzero , /* the adjusted offset */
530 #if 0
531 ADN_ttorg , 0.0 , /* in case not already set */
532 ADN_ttdur , 0.0 , /* in case not already set */
533 #endif
534 ADN_none ) ;
535
536 /*---- do the temporal shifting! ----*/
537
538 nup = csfft_nextup_one35( ntt+4 ) ;
539 ignore = TS_ignore ;
540
541 if( TS_verbose && SHIFT_get_method() == MRI_FOURIER )
542 printf("++ Time series length = %d; FFT length set to %d\n",ntt,nup) ;
543
544 for( kk=0 ; kk < nzz ; kk++ ){ /* loop over slices */
545
546 if( TS_fset == NULL ){
547 fshift = (TS_tzero - TS_tpat[kk]) / TS_TR ; /* rightward fractional shift */
548 #if 1
549 if( !BAD ) fshift = -fshift ; /* 24 Apr 2003 -- OOG */
550 #endif
551 gshift = fshift ;
552
553 if( TS_verbose )
554 printf("++ slice %d: fractional shift = %g\n",kk,fshift) ;
555
556 /* if removing any trend, process the slice 22 May 2008 [rickr] */
557 if( fabs(fshift) < 0.001 && ! TS_rlt ) continue ; /* skip this slice */
558 }
559
560 for( ii=0 ; ii < nxy ; ii+=2 ){ /* loop over voxel pairs in slice */
561
562 flim = THD_extract_series( ii+kk*nxy , TS_oset , 0 ) ; /* get this voxel */
563 far = MRI_FLOAT_PTR(flim) ;
564
565 if( TS_fset != NULL ){ /* 21 Sep 2011 */
566 MRI_IMAGE *qim = THD_extract_series( ii+kk*nxy , TS_fset , 0 ) ;
567 float *qar = MRI_FLOAT_PTR(qim) ; fshift = qar[0] ; mri_free(qim) ;
568 }
569
570 if( TS_rlt == 0 ){
571 for( ffmin=ffmax=far[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
572 if( far[jj] < ffmin ) ffmin = far[jj] ;
573 else if( far[jj] > ffmax ) ffmax = far[jj] ;
574 }
575 }
576
577 if( TS_detrend ) /* remove trend, else mean 3 Jan 2007 [rickr] */
578 THD_linear_detrend( ntt-ignore , far+ignore , &f0,&f1 ) ;
579 else
580 THD_const_detrend( ntt-ignore , far+ignore , &f0 ) ;
581
582 for( fmin=fmax=far[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
583 if( far[jj] < fmin ) fmin = far[jj] ; /* range of data: after */
584 else if( far[jj] > fmax ) fmax = far[jj] ;
585 }
586
587 if( ii < nxy-1 ){ /* get next voxel */
588 glim = THD_extract_series( ii+kk*nxy+1 , TS_oset , 0 ) ;
589 gar = MRI_FLOAT_PTR(glim) ;
590
591 if( TS_fset != NULL ){ /* 21 Sep 2011 */
592 MRI_IMAGE *qim = THD_extract_series( ii+kk*nxy+1 , TS_fset , 0 ) ;
593 float *qar = MRI_FLOAT_PTR(qim) ; gshift = qar[0] ; mri_free(qim) ;
594 }
595
596 if( TS_rlt == 0 ){
597 for( ggmin=ggmax=gar[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
598 if( gar[jj] < ggmin ) ggmin = gar[jj] ;
599 else if( gar[jj] > ggmax ) ggmax = gar[jj] ;
600 }
601 }
602
603 if( TS_detrend )
604 THD_linear_detrend( ntt-ignore , gar+ignore , &g0,&g1 ) ;
605 else
606 THD_const_detrend( ntt-ignore , far+ignore , &g0 ) ;
607
608 for( gmin=gmax=gar[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
609 if( gar[jj] < gmin ) gmin = gar[jj] ;
610 else if( gar[jj] > gmax ) gmax = gar[jj] ;
611 }
612 } else {
613 gar = NULL ;
614 }
615
616 if( gar != NULL )
617 SHIFT_two_rows( ntt-ignore,nup, fshift,far+ignore , gshift, gar+ignore ) ;
618 else
619 SHIFT_two_rows( ntt-ignore,nup, fshift,far+ignore , gshift, NULL ) ;
620
621 for( jj=ignore ; jj < ntt ; jj++ ){
622 if( far[jj] < fmin ) far[jj] = fmin ; /* clip to input range */
623 else if( far[jj] > fmax ) far[jj] = fmax ;
624 switch( TS_rlt ){ /* restore trend? */
625 case 0:
626 far[jj] += (f0 + (jj-ignore)*f1) ;
627 if( far[jj] < ffmin ) far[jj] = ffmin ;
628 else if( far[jj] > ffmax ) far[jj] = ffmax ;
629 break ;
630
631 case 2:
632 far[jj] += f0 ;
633 break ;
634 }
635 }
636
637 if( gar != NULL ){
638 for( jj=ignore ; jj < ntt ; jj++ ){
639 if( gar[jj] < gmin ) gar[jj] = gmin ;
640 else if( gar[jj] > gmax ) gar[jj] = gmax ;
641 switch( TS_rlt ){
642 case 0:
643 gar[jj] += (g0 + (jj-ignore)*g1) ;
644 if( gar[jj] < ggmin ) gar[jj] = ggmin ;
645 else if( gar[jj] > ggmax ) gar[jj] = ggmax ;
646 break ;
647
648 case 2:
649 gar[jj] += g0 ;
650 break ;
651 }
652 }
653 }
654
655 /* put back into dataset */
656
657 THD_insert_series( ii+kk*nxy , TS_oset , ntt , MRI_float , far , 0 ) ;
658 if( gar != NULL )
659 THD_insert_series( ii+kk*nxy+1 , TS_oset , ntt , MRI_float , gar , 0 ) ;
660
661 /* throw out the trash */
662
663 mri_free(flim) ; if( gar != NULL ) mri_free(glim) ;
664 }
665 }
666
667 if( TS_fset != NULL ) DSET_delete( TS_fset ) ;
668
669 DSET_write( TS_oset ) ;
670 if( TS_verbose ) fprintf(stderr,"++ Wrote output: %s\n",DSET_BRIKNAME(TS_oset)) ;
671 exit(0) ;
672 }
673
674 /*-----------------------------------------------------------------------------*/
675
TS_copy_input_to_output(void)676 static void TS_copy_input_to_output(void)
677 {
678 WARNING_message("==>> output dataset is just a copy of input dataset") ;
679 TS_oset = EDIT_full_copy( TS_dset , TS_prefix ) ;
680 if( TS_oset == NULL ) ERROR_exit("Can't copy input dataset '%s'",DSET_HEADNAME(TS_dset)) ;
681 DSET_unload( TS_dset ) ;
682 DSET_write( TS_oset ) ;
683 if( TS_verbose ) INFO_message("Wrote output: %s",DSET_HEADNAME(TS_oset)) ;
684 exit(0) ;
685 }
686