00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011 void TS_syntax(char * str)
00012 {
00013 if( str != NULL ){
00014 fprintf(stderr,"*** Fatal Error: %s\n",str) ; exit(1) ;
00015 }
00016
00017 printf("Usage: 3dTshift [options] dataset\n"
00018 "Shifts voxel time series from the input dataset so that the separate\n"
00019 "slices are aligned to the same temporal origin. By default, uses the\n"
00020 "slicewise shifting information in the dataset header (from the 'tpattern'\n"
00021 "input to program to3d).\n"
00022 "\n"
00023 "Method: detrend -> interpolate -> retrend (optionally)\n"
00024 "\n"
00025 "The input dataset can have a sub-brick selector attached, as documented\n"
00026 "in '3dcalc -help'.\n"
00027 "\n"
00028 "The output dataset time series will be interpolated from the input to\n"
00029 "the new temporal grid. This may not be the best way to analyze your\n"
00030 "data, but it can be convenient.\n"
00031 "\n"
00032 "Warnings:\n"
00033 "* Please recall the phenomenon of 'aliasing': frequencies above 1/(2*TR) can't\n"
00034 " be properly interpolated. For most 3D FMRI data, this means that cardiac\n"
00035 " and respiratory effects will not be treated properly by this program.\n"
00036 "\n"
00037 "* The images at the beginning of a high-speed FMRI imaging run are usually\n"
00038 " of a different quality than the later images, due to transient effects\n"
00039 " before the longitudinal magnetization settles into a steady-state value.\n"
00040 " These images should not be included in the interpolation! For example,\n"
00041 " if you wish to exclude the first 4 images, then the input dataset should\n"
00042 " be specified in the form 'prefix+orig[4..$]'. Alternatively, you can\n"
00043 " use the '-ignore ii' option.\n"
00044 "\n"
00045 "* It seems to be best to use 3dTshift before using 3dvolreg.\n"
00046 "\n"
00047 "Options:\n"
00048 " -verbose = print lots of messages while program runs\n"
00049 "\n"
00050 " -TR ddd = use 'ddd' as the TR, rather than the value\n"
00051 " stored in the dataset header using to3d.\n"
00052 " You may attach the suffix 's' for seconds,\n"
00053 " or 'ms' for milliseconds.\n"
00054 "\n"
00055 " -tzero zzz = align each slice to time offset 'zzz';\n"
00056 " the value of 'zzz' must be between the\n"
00057 " minimum and maximum slice temporal offsets.\n"
00058 " N.B.: The default alignment time is the average\n"
00059 " of the 'tpattern' values (either from the\n"
00060 " dataset header or from the -tpattern option)\n"
00061 "\n"
00062 " -slice nnn = align each slice to the time offset of slice\n"
00063 " number 'nnn' - only one of the -tzero and\n"
00064 " -slice options can be used.\n"
00065 "\n"
00066 " -prefix ppp = use 'ppp' for the prefix of the output file;\n"
00067 " the default is 'tshift'.\n"
00068 "\n"
00069 " -ignore ii = Ignore the first 'ii' points. (Default is ii=0.)\n"
00070 " The first ii values will be unchanged in the output\n"
00071 " (regardless of the -rlt option). They also will\n"
00072 " not be used in the detrending or time shifting.\n"
00073 "\n"
00074 " -rlt = Before shifting, the mean and linear trend\n"
00075 " -rlt+ = of each time series is removed. The default\n"
00076 " action is to add these back in after shifting.\n"
00077 " -rlt means to leave both of these out of the output\n"
00078 " -rlt+ means to add only the mean back into the output\n"
00079 " (cf. '3dTcat -help')\n"
00080 "\n"
00081 " -Fourier = Use a Fourier method (the default: most accurate; slowest).\n"
00082 " -linear = Use linear (1st order polynomial) interpolation (least accurate).\n"
00083 " -cubic = Use the cubic (3rd order) Lagrange polynomial interpolation.\n"
00084 " -quintic = Use the quintic (5th order) Lagrange polynomial interpolation.\n"
00085 " -heptic = Use the heptic (7th order) Lagrange polynomial interpolation.\n"
00086 "\n"
00087 " -tpattern ttt = use 'ttt' as the slice time pattern, rather\n"
00088 " than the pattern in the input dataset header;\n"
00089 " 'ttt' can have any of the values that would\n"
00090 " go in the 'tpattern' input to to3d, described below:\n"
00091 "\n"
00092 " alt+z = altplus = alternating in the plus direction\n"
00093 " alt+z2 = alternating, starting at slice #1 instead of #0\n"
00094 " alt-z = altminus = alternating in the minus direction\n"
00095 " alt-z2 = alternating, starting at slice #nz-2 instead of #nz-1\n"
00096 " seq+z = seqplus = sequential in the plus direction\n"
00097 " seq-z = seqminus = sequential in the minus direction\n"
00098 " @filename = read temporal offsets from 'filename'\n"
00099 "\n"
00100 " For example if nz = 5 and TR = 1000, then the inter-slice\n"
00101 " time is taken to be dt = TR/nz = 200. In this case, the\n"
00102 " slices are offset in time by the following amounts:\n"
00103 "\n"
00104 " S L I C E N U M B E R\n"
00105 " tpattern 0 1 2 3 4 Comment\n"
00106 " --------- --- --- --- --- --- -------------------------------\n"
00107 " altplus 0 600 200 800 400 Alternating in the +z direction\n"
00108 " alt+z2 400 0 600 200 800 Alternating, but starting at #1\n"
00109 " altminus 400 800 200 600 0 Alternating in the -z direction\n"
00110 " alt-z2 800 200 600 0 400 Alternating, starting at #nz-2 \n"
00111 " seqplus 0 200 400 600 800 Sequential in the -z direction\n"
00112 " seqplus 800 600 400 200 0 Sequential in the -z direction\n"
00113 "\n"
00114 " If @filename is used for tpattern, then nz ASCII-formatted numbers\n"
00115 " are read from the file. These indicate the time offsets for each\n"
00116 " slice. For example, if 'filename' contains\n"
00117 " 0 600 200 800 400\n"
00118 " then this is equivalent to 'altplus' in the above example.\n"
00119 " (nz = number of slices in the input dataset)\n"
00120 "\n"
00121 "N.B.: if you are using -tpattern, make sure that the units supplied\n"
00122 " match the units of TR in the dataset header, or provide a\n"
00123 " new TR using the -TR option.\n"
00124 "\n"
00125 "As a test of how well 3dTshift interpolates, you can take a dataset\n"
00126 "that was created with '-tpattern alt+z', run 3dTshift on it, and\n"
00127 "then run 3dTshift on the new dataset with '-tpattern alt-z' -- the\n"
00128 "effect will be to reshift the dataset back to the original time\n"
00129 "grid. Comparing the original dataset to the shifted-then-reshifted\n"
00130 "output will show where 3dTshift does a good job and where it does\n"
00131 "a bad job.\n"
00132 "\n"
00133 "-- RWCox - 31 October 1999\n"
00134 ) ;
00135
00136 printf("\n" MASTER_SHORTHELP_STRING ) ;
00137
00138 exit(0) ;
00139 }
00140
00141
00142
00143 float * TS_parse_tpattern( int nzz , float TR , char * tpattern )
00144 {
00145 int ii ;
00146 float tframe , tsl ;
00147 float * tpat ;
00148
00149 tpat = (float *) malloc( sizeof(float) * nzz ) ;
00150 for( ii=0 ; ii < nzz ; ii++ ) tpat[ii] = 0.0 ;
00151
00152 tframe = TR / nzz ;
00153
00154 if( nzz > 1 && tpattern[0] == '@' ){
00155 FILE * fp ;
00156
00157
00158
00159 fp = fopen( tpattern+1 , "r" ) ;
00160 if( fp == NULL ){
00161 fprintf(stderr,"*** Cannot open tpattern file %s\n",tpattern+1) ;
00162 exit(1) ;
00163 } else {
00164 for( ii=0 ; ii < nzz ; ii++ ){
00165 fscanf( fp , "%f" , tpat + ii ) ;
00166 if( tpat[ii] < 0.0 || tpat[ii] > TR ){
00167 fprintf(stderr,"*** Illegal value %g in tpattern file %s\n",
00168 tpat[ii] , tpattern+1 ) ;
00169 exit(1) ;
00170 }
00171 }
00172 fclose( fp ) ;
00173 }
00174
00175 } else if( nzz > 1 &&
00176 (strcmp(tpattern,"alt+z")==0 || strcmp(tpattern,"altplus")==0) ){
00177
00178
00179
00180 tsl = 0.0 ;
00181 for( ii=0 ; ii < nzz ; ii+=2 ){
00182 tpat[ii] = tsl ; tsl += tframe ;
00183 }
00184 for( ii=1 ; ii < nzz ; ii+=2 ){
00185 tpat[ii] = tsl ; tsl += tframe ;
00186 }
00187
00188 } else if( nzz > 1 && strcmp(tpattern,"alt+z2")==0 ){
00189
00190
00191
00192 tsl = 0.0 ;
00193 for( ii=1 ; ii < nzz ; ii+=2 ){
00194 tpat[ii] = tsl ; tsl += tframe ;
00195 }
00196 for( ii=0 ; ii < nzz ; ii+=2 ){
00197 tpat[ii] = tsl ; tsl += tframe ;
00198 }
00199
00200 } else if( nzz > 1 &&
00201 (strcmp(tpattern,"alt-z")==0 || strcmp(tpattern,"altminus")==0) ){
00202
00203
00204
00205 tsl = 0.0 ;
00206 for( ii=nzz-1 ; ii >=0 ; ii-=2 ){
00207 tpat[ii] = tsl ; tsl += tframe ;
00208 }
00209 for( ii=nzz-2 ; ii >=0 ; ii-=2 ){
00210 tpat[ii] = tsl ; tsl += tframe ;
00211 }
00212
00213 } else if( nzz > 1 && strcmp(tpattern,"alt-z2") == 0 ){
00214
00215
00216
00217 tsl = 0.0 ;
00218 for( ii=nzz-2 ; ii >=0 ; ii-=2 ){
00219 tpat[ii] = tsl ; tsl += tframe ;
00220 }
00221 for( ii=nzz-1 ; ii >=0 ; ii-=2 ){
00222 tpat[ii] = tsl ; tsl += tframe ;
00223 }
00224
00225 } else if( nzz > 1 &&
00226 (strcmp(tpattern,"seq+z")==0 || strcmp(tpattern,"seqplus")==0) ){
00227
00228
00229
00230 tsl = 0.0 ;
00231 for( ii=0 ; ii < nzz ; ii++ ){
00232 tpat[ii] = tsl ; tsl += tframe ;
00233 }
00234
00235 } else if( nzz > 1 &&
00236 (strcmp(tpattern,"seq-z")==0 || strcmp(tpattern,"seqminus")==0) ){
00237
00238
00239
00240 tsl = 0.0 ;
00241 for( ii=nzz-1 ; ii >=0 ; ii-- ){
00242 tpat[ii] = tsl ; tsl += tframe ;
00243 }
00244
00245 } else if( nzz == 1 ||
00246 (strcmp(tpattern,"zero")==0 || strcmp(tpattern,"simult")==0) ){
00247
00248
00249
00250 } else {
00251 fprintf(stderr,"*** Unknown tpattern = %s\n",tpattern) ;
00252 exit(1) ;
00253 }
00254
00255 return tpat ;
00256 }
00257
00258
00259
00260 static float TS_TR = 0.0 ;
00261 static int TS_tunits = UNITS_SEC_TYPE ;
00262 static float * TS_tpat = NULL ;
00263 static float TS_tzero = -1.0 ;
00264 static int TS_slice = -1 ;
00265 static int TS_rlt = 0 ;
00266
00267 static int TS_verbose = 0 ;
00268
00269 static int TS_ignore = 0 ;
00270
00271 static THD_3dim_dataset * TS_dset = NULL , * TS_oset = NULL ;
00272
00273 static char * TS_tpattern = NULL ;
00274
00275 static char TS_prefix[THD_MAX_PREFIX] = "tshift" ;
00276
00277
00278
00279 int main( int argc , char *argv[] )
00280 {
00281 int nopt=1 ;
00282 int nzz, ii,jj,kk , ntt,nxx,nyy,nxy , nup ;
00283 float tomax,tomin , tshift , fmin,fmax , gmin,gmax , f0,f1 , g0,g1 ;
00284 float ffmin,ffmax , ggmin,ggmax ;
00285 MRI_IMAGE * flim , * glim ;
00286 float * far , * gar ;
00287 int ignore=0 , BAD=0 ;
00288
00289
00290
00291 if( argc < 2 || strcmp(argv[1],"-help") == 0 ) TS_syntax(NULL) ;
00292
00293 mainENTRY("3dTshift main"); machdep(); AFNI_logger("3dTshift",argc,argv);
00294 PRINT_VERSION("3dTshift");
00295
00296 SHIFT_set_method( MRI_FOURIER ) ;
00297
00298 while( nopt < argc && argv[nopt][0] == '-' ){
00299
00300 if( strcmp(argv[nopt],"-BAD") == 0 ){
00301 BAD = 1 ; nopt++ ; continue ;
00302 }
00303
00304 if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00305 TS_verbose++ ;
00306 nopt++ ; continue ;
00307 }
00308
00309 if( strcmp(argv[nopt],"-ignore") == 0 ){
00310 TS_ignore = (int) strtod(argv[++nopt],NULL) ;
00311 if( TS_ignore < 0 ) TS_syntax("-ignore value is negative!") ;
00312 nopt++ ; continue ;
00313 }
00314
00315 if( strncmp(argv[nopt],"-Fourier",4) == 0 || strncmp(argv[nopt],"-fourier",4) == 0 ){
00316 SHIFT_set_method( MRI_FOURIER ) ;
00317 nopt++ ; continue ;
00318 }
00319
00320 if( strncmp(argv[nopt],"-cubic",4) == 0 || strncmp(argv[nopt],"-Cubic",4) == 0 ){
00321 SHIFT_set_method( MRI_CUBIC ) ;
00322 nopt++ ; continue ;
00323 }
00324
00325 if( strncmp(argv[nopt],"-quintic",4) == 0 || strncmp(argv[nopt],"-Quintic",4) == 0 ){
00326 SHIFT_set_method( MRI_QUINTIC ) ;
00327 nopt++ ; continue ;
00328 }
00329
00330 if( strncmp(argv[nopt],"-heptic",4) == 0 || strncmp(argv[nopt],"-Heptic",4) == 0 ){
00331 SHIFT_set_method( MRI_HEPTIC ) ;
00332 nopt++ ; continue ;
00333 }
00334
00335 if( strncmp(argv[nopt],"-linear",4) == 0 || strncmp(argv[nopt],"-Linear",4) == 0 ){
00336 SHIFT_set_method( MRI_LINEAR ) ;
00337 nopt++ ; continue ;
00338 }
00339
00340 if( strcmp(argv[nopt],"-TR") == 0 ){
00341 char * eptr ;
00342 if( ++nopt >= argc ) TS_syntax("-TR needs an argument!") ;
00343 TS_TR = strtod( argv[nopt] , &eptr ) ;
00344 if( TS_TR <= 0.0 ) TS_syntax("illegal value after -TR!") ;
00345
00346 if( strcmp(eptr,"ms")==0 || strcmp(eptr,"msec")==0 ){
00347 TS_tunits = UNITS_MSEC_TYPE ;
00348 } else if( strcmp(eptr,"s")==0 || strcmp(eptr,"sec")==0 ){
00349 TS_tunits = UNITS_SEC_TYPE ;
00350 } else if( strcmp(eptr,"Hz")==0 || strcmp(eptr,"Hertz")==0 ){
00351 TS_tunits = UNITS_HZ_TYPE ;
00352 }
00353
00354 nopt++ ; continue ;
00355 }
00356
00357 if( strcmp(argv[nopt],"-tzero") == 0 ){
00358 if( ++nopt >= argc ) TS_syntax("-tzero needs an argument!") ;
00359 TS_tzero = strtod( argv[nopt] , NULL ) ;
00360 if( TS_tzero < 0.0 ) TS_syntax("illegal value after -tzero!") ;
00361 nopt++ ; continue ;
00362 }
00363
00364 if( strcmp(argv[nopt],"-slice") == 0 ){
00365 if( ++nopt >= argc ) TS_syntax("-slice needs an argument!") ;
00366 TS_slice = strtod( argv[nopt] , NULL ) ;
00367 if( TS_slice < 0 ) TS_syntax("illegal value after -slice!") ;
00368 nopt++ ; continue ;
00369 }
00370
00371 if( strcmp(argv[nopt],"-prefix") == 0 ){
00372 if( ++nopt >= argc ) TS_syntax("-prefix needs an argument!") ;
00373 MCW_strncpy( TS_prefix , argv[nopt] , THD_MAX_PREFIX ) ;
00374 if( !THD_filename_ok(TS_prefix) ) TS_syntax("illegal value after -prefix") ;
00375 nopt++ ; continue ;
00376 }
00377
00378 if( strcmp(argv[nopt],"-rlt") == 0 ){
00379 TS_rlt = 1 ;
00380 nopt++ ; continue ;
00381 }
00382
00383 if( strcmp(argv[nopt],"-rlt+") == 0 ){
00384 TS_rlt = 2 ;
00385 nopt++ ; continue ;
00386 }
00387
00388 if( strcmp(argv[nopt],"-tpattern") == 0 ){
00389 if( ++nopt >= argc ) TS_syntax("-tpattern needs an argument!") ;
00390 TS_tpattern = argv[nopt] ;
00391 nopt++ ; continue ;
00392 }
00393
00394 fprintf(stderr,"*** Unknown option: %s\n",argv[nopt]) ; exit(1) ;
00395
00396 }
00397
00398
00399
00400 if( nopt >= argc ) TS_syntax("Need a dataset input?!") ;
00401 if( TS_verbose ) printf("++ opening input dataset header\n") ;
00402 TS_dset = THD_open_dataset( argv[nopt] ) ;
00403 if( TS_dset == NULL ) TS_syntax("Can't open input dataset!") ;
00404
00405 nxx = DSET_NX(TS_dset) ;
00406 nyy = DSET_NY(TS_dset) ; nxy = nxx * nyy ;
00407 nzz = DSET_NZ(TS_dset) ;
00408 ntt = DSET_NVALS(TS_dset) ;
00409
00410 if( DSET_NVALS(TS_dset) < 2 ) TS_syntax("Dataset has only 1 value per voxel!") ;
00411 if( TS_slice >= nzz ) TS_syntax("-slice value is too large") ;
00412
00413 if( TS_ignore > ntt-5 ) TS_syntax("-ignore value is too large") ;
00414
00415 if( TS_dset->taxis == NULL ){
00416 if( TS_TR == 0.0 || TS_tpattern == NULL )
00417 TS_syntax("dataset has no time axis => you must supply -TR and -tpattern!") ;
00418
00419 } else if( TS_tpattern == NULL && TS_dset->taxis->toff_sl == NULL ){
00420 TS_syntax("dataset is already aligned in time!") ;
00421 }
00422
00423 if( TS_TR == 0.0 ){
00424 TS_TR = DSET_TIMESTEP(TS_dset) ;
00425 TS_tunits = TS_dset->taxis->units_type ;
00426 if( TS_verbose )
00427 printf("++ using dataset TR = %g %s\n",TS_TR,UNITS_TYPE_LABEL(TS_tunits)) ;
00428 }
00429
00430 if( TS_tpattern != NULL ){
00431 TS_tpat = TS_parse_tpattern( nzz , TS_TR , TS_tpattern ) ;
00432 } else {
00433 if( TS_dset->taxis->nsl != nzz )
00434 TS_syntax("dataset temporal pattern is malformed!") ;
00435
00436 TS_tpat = (float *) malloc( sizeof(float) * nzz ) ;
00437 memcpy( TS_tpat , TS_dset->taxis->toff_sl , sizeof(float)*nzz ) ;
00438 }
00439 if( TS_verbose ){
00440 printf("++ using tpattern = ") ;
00441 for( ii=0 ; ii < nzz ; ii++ ) printf("%g ",TS_tpat[ii]) ;
00442 printf("%s\n",UNITS_TYPE_LABEL(TS_tunits)) ;
00443 }
00444
00445 tomin = WAY_BIG ; tomax = -WAY_BIG ;
00446 for( ii=0 ; ii < nzz ; ii++ ){
00447 if( TS_tpat[ii] > tomax ) tomax = TS_tpat[ii] ;
00448 if( TS_tpat[ii] < tomin ) tomin = TS_tpat[ii] ;
00449 }
00450 if( tomin < 0.0 || tomax > TS_TR )
00451 TS_syntax("some value in tpattern is outside 0..TR") ;
00452 else if( tomin >= tomax )
00453 TS_syntax("temporal pattern is already aligned in time!") ;
00454
00455 if( TS_slice >= 0 && TS_slice < nzz ){
00456 TS_tzero = TS_tpat[TS_slice] ;
00457 } else if( TS_tzero < 0.0 ){
00458 TS_tzero = 0.0 ;
00459 for( ii=0 ; ii < nzz ; ii++ ) TS_tzero += TS_tpat[ii] ;
00460 TS_tzero /= nzz ;
00461 }
00462 if( TS_verbose ) printf("++ common time point set to %g\n",TS_tzero) ;
00463
00464
00465
00466 if( TS_verbose ) printf("++ copying input dataset bricks\n") ;
00467
00468 TS_oset = EDIT_full_copy( TS_dset , TS_prefix ) ;
00469 if( TS_oset == NULL ) TS_syntax("Can't copy input dataset!") ;
00470 DSET_unload( TS_dset ) ;
00471
00472 if( THD_is_file(DSET_HEADNAME(TS_oset)) )
00473 TS_syntax("output dataset already exists!") ;
00474
00475 tross_Copy_History( TS_dset , TS_oset ) ;
00476 tross_Make_History( "3dTshift" , argc,argv , TS_oset ) ;
00477
00478
00479
00480 EDIT_dset_items( TS_oset ,
00481 ADN_ntt , ntt ,
00482 ADN_ttdel , TS_TR ,
00483 ADN_tunits , TS_tunits ,
00484 ADN_nsl , 0 ,
00485 #if 0
00486 ADN_ttorg , 0.0 ,
00487 ADN_ttdur , 0.0 ,
00488 #endif
00489 ADN_none ) ;
00490
00491
00492
00493 nup = csfft_nextup_one35( ntt+4 ) ;
00494 ignore = TS_ignore ;
00495
00496 if( TS_verbose && SHIFT_get_method() == MRI_FOURIER )
00497 printf("++ Time series length = %d; FFT length set to %d\n",ntt,nup) ;
00498
00499 for( kk=0 ; kk < nzz ; kk++ ){
00500
00501 tshift = (TS_tzero - TS_tpat[kk]) / TS_TR ;
00502 #if 1
00503 if( !BAD ) tshift = -tshift ;
00504 #endif
00505
00506 if( TS_verbose )
00507 printf("++ slice %d: fractional shift = %g\n",kk,tshift) ;
00508
00509 if( fabs(tshift) < 0.001 ) continue ;
00510
00511 for( ii=0 ; ii < nxy ; ii+=2 ){
00512
00513 flim = THD_extract_series( ii+kk*nxy , TS_oset , 0 ) ;
00514 far = MRI_FLOAT_PTR(flim) ;
00515
00516 if( TS_rlt == 0 ){
00517 for( ffmin=ffmax=far[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
00518 if( far[jj] < ffmin ) ffmin = far[jj] ;
00519 else if( far[jj] > ffmax ) ffmax = far[jj] ;
00520 }
00521 }
00522
00523 THD_linear_detrend( ntt-ignore , far+ignore , &f0,&f1 ) ;
00524
00525 for( fmin=fmax=far[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
00526 if( far[jj] < fmin ) fmin = far[jj] ;
00527 else if( far[jj] > fmax ) fmax = far[jj] ;
00528 }
00529
00530 if( ii < nxy-1 ){
00531 glim = THD_extract_series( ii+kk*nxy+1 , TS_oset , 0 ) ;
00532 gar = MRI_FLOAT_PTR(glim) ;
00533
00534 if( TS_rlt == 0 ){
00535 for( ggmin=ggmax=gar[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
00536 if( gar[jj] < ggmin ) ggmin = gar[jj] ;
00537 else if( gar[jj] > ggmax ) ggmax = gar[jj] ;
00538 }
00539 }
00540
00541 THD_linear_detrend( ntt-ignore , gar+ignore , &g0,&g1 ) ;
00542
00543 for( gmin=gmax=gar[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
00544 if( gar[jj] < gmin ) gmin = gar[jj] ;
00545 else if( gar[jj] > gmax ) gmax = gar[jj] ;
00546 }
00547 } else {
00548 gar = NULL ;
00549 }
00550
00551 if( gar != NULL )
00552 SHIFT_two_rows( ntt-ignore,nup, tshift,far+ignore , tshift, gar+ignore ) ;
00553 else
00554 SHIFT_two_rows( ntt-ignore,nup, tshift,far+ignore , tshift, NULL ) ;
00555
00556 for( jj=ignore ; jj < ntt ; jj++ ){
00557 if( far[jj] < fmin ) far[jj] = fmin ;
00558 else if( far[jj] > fmax ) far[jj] = fmax ;
00559 switch( TS_rlt ){
00560 case 0:
00561 far[jj] += (f0 + (jj-ignore)*f1) ;
00562 if( far[jj] < ffmin ) far[jj] = ffmin ;
00563 else if( far[jj] > ffmax ) far[jj] = ffmax ;
00564 break ;
00565
00566 case 2:
00567 far[jj] += f0 ;
00568 break ;
00569 }
00570 }
00571
00572 if( gar != NULL ){
00573 for( jj=ignore ; jj < ntt ; jj++ ){
00574 if( gar[jj] < gmin ) gar[jj] = gmin ;
00575 else if( gar[jj] > gmax ) gar[jj] = gmax ;
00576 switch( TS_rlt ){
00577 case 0:
00578 gar[jj] += (g0 + (jj-ignore)*g1) ;
00579 if( gar[jj] < ggmin ) gar[jj] = ggmin ;
00580 else if( gar[jj] > ggmax ) gar[jj] = ggmax ;
00581 break ;
00582
00583 case 2:
00584 gar[jj] += g0 ;
00585 break ;
00586 }
00587 }
00588 }
00589
00590
00591
00592 THD_insert_series( ii+kk*nxy , TS_oset , ntt , MRI_float , far , 0 ) ;
00593 if( gar != NULL )
00594 THD_insert_series( ii+kk*nxy+1 , TS_oset , ntt , MRI_float , gar , 0 ) ;
00595
00596
00597
00598 mri_free(flim) ; if( gar != NULL ) mri_free(glim) ;
00599 }
00600 }
00601
00602 DSET_write( TS_oset ) ;
00603 if( TS_verbose ) fprintf(stderr,"++ Wrote output: %s\n",DSET_BRIKNAME(TS_oset)) ;
00604 exit(0) ;
00605 }
