00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008 #include "parser.h"
00009
00010
00011
00012
00013
00014
00015
00016 #define ALLOW_BYSLICE
00017
00018 static THD_3dim_dataset * DT_dset = NULL ;
00019 static MRI_IMARR * DT_imar = NULL ;
00020 static char ** DT_expr = NULL ;
00021 static PARSER_code ** DT_excode = NULL ;
00022 static float * DT_exdel = NULL ;
00023 static int * DT_exvar = NULL ;
00024 static int DT_exnum = 0 ;
00025 static int DT_verb = 0 ;
00026 static int DT_replace = 0 ;
00027 static int DT_norm = 0 ;
00028 static int DT_byslice = 0 ;
00029 static int DT_nvector = 0 ;
00030
00031 static float DT_current_del = -1.0 ;
00032
00033 static char DT_output_prefix[THD_MAX_PREFIX] = "detrend" ;
00034 static char DT_session[THD_MAX_NAME] = "./" ;
00035
00036
00037
00038 void DT_read_opts( int , char ** ) ;
00039 void DT_Syntax(void) ;
00040
00041
00042
00043
00044
00045 void DT_read_opts( int argc , char * argv[] )
00046 {
00047 int nopt = 1 , nvals , ii , nvcheck ;
00048
00049 INIT_IMARR(DT_imar) ;
00050
00051 while( nopt < argc && argv[nopt][0] == '-' ){
00052
00053
00054
00055 if( strncmp(argv[nopt],"-prefix",6) == 0 ){
00056 nopt++ ;
00057 if( nopt >= argc ){
00058 fprintf(stderr,"*** need argument after -prefix!\n") ; exit(1) ;
00059 }
00060 MCW_strncpy( DT_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00061 continue ;
00062 }
00063
00064
00065
00066 if( strncmp(argv[nopt],"-session",6) == 0 ){
00067 nopt++ ;
00068 if( nopt >= argc ){
00069 fprintf(stderr,"*** need argument after -session!\n") ; exit(1) ;
00070 }
00071 MCW_strncpy( DT_session , argv[nopt++] , THD_MAX_NAME ) ;
00072 continue ;
00073 }
00074
00075
00076
00077 if( strncmp(argv[nopt],"-verb",5) == 0 ){
00078 DT_verb++ ;
00079 nopt++ ; continue ;
00080 }
00081
00082
00083
00084 if( strncmp(argv[nopt],"-replace",5) == 0 ){
00085 DT_replace++ ;
00086 nopt++ ; continue ;
00087 }
00088
00089 #ifdef ALLOW_BYSLICE
00090
00091
00092 if( strncmp(argv[nopt],"-byslice",5) == 0 ){
00093 DT_byslice++ ;
00094 nopt++ ; continue ;
00095 }
00096 #endif
00097
00098
00099
00100 if( strncmp(argv[nopt],"-normalize",5) == 0 ){
00101 DT_norm++ ;
00102 nopt++ ; continue ;
00103 }
00104
00105
00106
00107 if( strncmp(argv[nopt],"-vector",4) == 0 ){
00108 MRI_IMAGE * flim ;
00109 nopt++ ;
00110 if( nopt >= argc ){
00111 fprintf(stderr,"*** need argument after -vector!\n"); exit(1);
00112 }
00113 flim = mri_read_1D( argv[nopt++] ) ;
00114 if( flim == NULL ){
00115 fprintf(stderr,"*** can't read -vector %s\n",argv[nopt-1]); exit(1);
00116 }
00117 ADDTO_IMARR(DT_imar,flim) ;
00118 if( DT_verb )
00119 fprintf(stderr,"+++ Read in %s: rows=%d cols=%d\n",
00120 argv[nopt-1],flim->ny,flim->nx ) ;
00121 continue ;
00122 }
00123
00124
00125
00126 if( strncmp(argv[nopt],"-del",4) == 0 ){
00127 nopt++ ;
00128 if( nopt >= argc ){
00129 fprintf(stderr,"*** need argument after -del!\n"); exit(1);
00130 }
00131 DT_current_del = strtod( argv[nopt++] , NULL ) ;
00132 if( DT_verb )
00133 fprintf(stderr,"+++ Set expression stepsize = %g\n",DT_current_del) ;
00134 continue ;
00135 }
00136
00137
00138
00139 if( strncmp(argv[nopt],"-expr",4) == 0 ){
00140 int nexp , qvar , kvar ;
00141 char sym[4] ;
00142
00143 nopt++ ;
00144 if( nopt >= argc ){
00145 fprintf(stderr,"*** need argument after -expr!\n"); exit(1);
00146 }
00147
00148 nexp = DT_exnum + 1 ;
00149 if( DT_exnum == 0 ){
00150 DT_expr = (char **) malloc( sizeof(char *) ) ;
00151 DT_excode = (PARSER_code **) malloc( sizeof(PARSER_code *) ) ;
00152 DT_exdel = (float *) malloc( sizeof(float) ) ;
00153 DT_exvar = (int *) malloc( sizeof(int) ) ;
00154 } else {
00155 DT_expr = (char **) realloc( DT_expr ,
00156 sizeof(char *)*nexp ) ;
00157 DT_excode = (PARSER_code **) realloc( DT_excode ,
00158 sizeof(PARSER_code *)*nexp ) ;
00159 DT_exdel = (float *) realloc( DT_exdel ,
00160 sizeof(float)*nexp) ;
00161 DT_exvar = (int *) realloc( DT_exvar ,
00162 sizeof(int)*nexp) ;
00163 }
00164 DT_expr[DT_exnum] = argv[nopt] ;
00165 DT_exdel[DT_exnum] = DT_current_del ;
00166 DT_excode[DT_exnum] = PARSER_generate_code( argv[nopt] ) ;
00167 if( DT_excode[DT_exnum] == NULL ){
00168 fprintf(stderr,"*** Illegal expression: %s\n",argv[nopt]); exit(1);
00169 }
00170
00171 qvar = 0 ; kvar = -1 ;
00172 for( ii=0 ; ii < 26 ; ii++ ){
00173 sym[0] = 'A' + ii ; sym[1] = '\0' ;
00174 if( PARSER_has_symbol(sym,DT_excode[DT_exnum]) ){
00175 qvar++ ; if( kvar < 0 ) kvar = ii ;
00176 if( DT_verb )
00177 fprintf(stderr,"+++ Found expression symbol %s\n",sym) ;
00178 }
00179 }
00180 if( qvar > 1 ){
00181 fprintf(stderr,"*** Expression %s should have just one symbol!\n",
00182 DT_expr[DT_exnum] ) ;
00183 exit(1) ;
00184 }
00185 DT_exvar[DT_exnum] = kvar ;
00186
00187 DT_exnum = nexp ; nopt++ ; continue ;
00188 }
00189
00190
00191
00192 fprintf(stderr,"*** Unknown option: %s\n",argv[nopt]) ; exit(1) ;
00193
00194 }
00195
00196
00197
00198 if( nopt >= argc ){
00199 fprintf(stderr,"*** No input dataset!?\n") ; exit(1) ;
00200 }
00201
00202 DT_nvector = IMARR_COUNT(DT_imar) ;
00203 if( DT_nvector + DT_exnum == 0 ){
00204 fprintf(stderr,"*** No -vector or -expr options!?\n") ; exit(1) ;
00205 }
00206 #ifdef ALLOW_BYSLICE
00207 if( DT_nvector == 0 && DT_byslice ){
00208 fprintf(stderr,"*** No -vector option supplied with -byslice!?\n"); exit(1);
00209 }
00210 #endif
00211
00212
00213
00214 DT_dset = THD_open_dataset( argv[nopt] ) ;
00215 if( DT_dset == NULL ){
00216 fprintf(stderr,"*** Can't open dataset %s\n",argv[nopt]) ; exit(1) ;
00217 }
00218
00219 DT_current_del = DSET_TR(DT_dset) ;
00220 if( DT_current_del <= 0.0 ){
00221 DT_current_del = 1.0 ;
00222 if( DT_verb )
00223 fprintf(stderr,"+++ Input has no TR value; setting TR=1.0\n") ;
00224 } else if( DT_verb ){
00225 fprintf(stderr,"+++ Input has TR=%g\n",DT_current_del) ;
00226 }
00227
00228
00229
00230 nvcheck = nvals = DSET_NVALS(DT_dset) ;
00231 #ifdef ALLOW_BYSLICE
00232 if( DT_byslice ) nvcheck *= DSET_NZ(DT_dset) ;
00233 #endif
00234 for( ii=0 ; ii < DT_nvector ; ii++ ){
00235 if( IMARR_SUBIMAGE(DT_imar,ii)->nx < nvcheck ){
00236 fprintf(stderr,"*** %d-th -vector is shorter than dataset!\n",ii+1) ;
00237 exit(1) ;
00238 }
00239 }
00240
00241
00242
00243 if( DT_exnum > 0 ){
00244 double atoz[26] , del ;
00245 int kvar , jj ;
00246 MRI_IMAGE * flim ;
00247 float * flar ;
00248
00249 for( jj=0 ; jj < DT_exnum ; jj++ ){
00250 if( DT_verb ) fprintf(stderr,"+++ Evaluating %d-th -expr\n",jj+1) ;
00251 kvar = DT_exvar[jj] ;
00252 del = DT_exdel[jj] ;
00253 if( del <= 0.0 ) del = DT_current_del ;
00254 flim = mri_new( nvals , 1 , MRI_float ) ;
00255 flar = MRI_FLOAT_PTR(flim) ;
00256 for( ii=0 ; ii < 26 ; ii++ ) atoz[ii] = 0.0 ;
00257 for( ii=0 ; ii < nvals ; ii++ ){
00258 if( kvar >= 0 ) atoz[kvar] = ii * del ;
00259 flar[ii] = PARSER_evaluate_one( DT_excode[jj] , atoz ) ;
00260 }
00261 ADDTO_IMARR( DT_imar , flim ) ;
00262 }
00263 }
00264
00265 return ;
00266 }
00267
00268
00269
00270 void DT_Syntax(void)
00271 {
00272 printf(
00273 "Usage: 3dDetrend [options] dataset\n"
00274 "This program removes components from voxel time series using\n"
00275 "linear least squares. Each voxel is treated independently.\n"
00276 "The input dataset may have a sub-brick selector string; otherwise,\n"
00277 "all sub-bricks will be used.\n\n"
00278 ) ;
00279
00280 printf(
00281 "General Options:\n"
00282 " -prefix pname = Use 'pname' for the output dataset prefix name.\n"
00283 " [default='detrend']\n"
00284 " -session dir = Use 'dir' for the output dataset session directory.\n"
00285 " [default='./'=current working directory]\n"
00286 " -verb = Print out some verbose output as the program runs.\n"
00287 " -replace = Instead of subtracting the fit from each voxel,\n"
00288 " replace the voxel data with the time series fit.\n"
00289 " -normalize = Normalize each output voxel time series; that is,\n"
00290 " make the sum-of-squares equal to 1.\n"
00291 " N.B.: This option is only valid if the input dataset is\n"
00292 " stored as floats!\n"
00293 #ifdef ALLOW_BYSLICE
00294 " -byslice = Treat each input vector (infra) as describing a set of\n"
00295 " time series interlaced across slices. If NZ is the\n"
00296 " number of slices and NT is the number of time points,\n"
00297 " then each input vector should have NZ*NT values when\n"
00298 " this option is used (usually, they only need NT values).\n"
00299 " The values must be arranged in slice order, then time\n"
00300 " order, in each vector column, as shown here:\n"
00301 " f(z=0,t=0) // first slice, first time\n"
00302 " f(z=1,t=0) // second slice, first time\n"
00303 " ...\n"
00304 " f(z=NZ-1,t=0) // last slice, first time\n"
00305 " f(z=0,t=1) // first slice, second time\n"
00306 " f(z=1,t=1) // second slice, second time\n"
00307 " ...\n"
00308 " f(z=NZ-1,t=NT-1) // last slice, last time\n"
00309 #endif
00310 "\n"
00311 "Component Options:\n"
00312 "These options determine the components that will be removed from\n"
00313 "each dataset voxel time series. They may be repeated to specify\n"
00314 "multiple regression. At least one component must be specified.\n"
00315 "\n"
00316 " -vector vvv = Remove components proportional to the columns vectors\n"
00317 " of the ASCII *.1D file 'vvv'. You may use a\n"
00318 " sub-vector selector string to specify which columns\n"
00319 " to use; otherwise, all columns will be used.\n"
00320 " For example:\n"
00321 " -vector 'xyzzy.1D[3,5]'\n"
00322 " will remove the 4th and 6th columns of file xyzzy.1D\n"
00323 " from the dataset (sub-vector indexes start at 0).\n"
00324 "\n"
00325 " -expr eee = Remove components proportional to the function\n"
00326 " specified in the expression string 'eee'.\n"
00327 " Any single letter from a-z may be used as the\n"
00328 " independent variable in 'eee'. For example:\n"
00329 " -expr 'cos(2*PI*t/40)' -expr 'sin(2*PI*t/40)'\n"
00330 " will remove sine and cosine waves of period 40\n"
00331 " from the dataset. Another example:\n"
00332 " -expr '1' -expr 't' -expr 't*t'\n"
00333 " will remove a quadratic trend from the data.\n"
00334 "\n"
00335 " -del ddd = Use the numerical value 'ddd' for the stepsize\n"
00336 " in subsequent -expr options. If no -del option\n"
00337 " is ever given, then the TR given in the dataset\n"
00338 " header is used for 'ddd'; if that isn't available,\n"
00339 " then 'ddd'=1.0 is assumed. The j-th time point\n"
00340 " will have independent variable = j * ddd, starting\n"
00341 " at j=0. For example:\n"
00342 " -expr 'sin(x)' -del 2.0 -expr 'z**3'\n"
00343 " means that the stepsize in 'sin(x)' is delta-x=TR,\n"
00344 " but the stepsize in 'z**3' is delta-z = 2.\n"
00345 #ifdef ALLOW_BYSLICE
00346 "\n"
00347 " N.B.: expressions are NOT calculated on a per-slice basis when the\n"
00348 " -byslice option is used. If you want to do this, you could\n"
00349 " compute vectors with the required time series using 1devel.\n"
00350 #endif
00351 ) ;
00352
00353 printf("\n" MASTER_SHORTHELP_STRING ) ;
00354
00355 exit(0) ;
00356 }
00357
00358
00359
00360 int main( int argc , char * argv[] )
00361 {
00362 int iv,nvals , nvec , ii,jj,kk , nvox ;
00363 THD_3dim_dataset * new_dset=NULL ;
00364 double * choleski ;
00365 float ** refvec , * fv , * fc , * fit ;
00366 MRI_IMAGE * flim ;
00367
00368
00369
00370 if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) DT_Syntax() ;
00371
00372
00373
00374 mainENTRY("3dDetrend main"); machdep() ; PRINT_VERSION("3dDetrend");
00375 { int new_argc ; char ** new_argv ;
00376 addto_args( argc , argv , &new_argc , &new_argv ) ;
00377 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00378 }
00379
00380 DT_read_opts( argc , argv ) ;
00381
00382
00383
00384 new_dset = EDIT_empty_copy( DT_dset ) ;
00385
00386
00387
00388 tross_Copy_History( DT_dset , new_dset ) ;
00389 tross_Make_History( "3dDetrend" , argc,argv , new_dset ) ;
00390
00391 EDIT_dset_items( new_dset ,
00392 ADN_prefix , DT_output_prefix ,
00393 ADN_directory_name, DT_session ,
00394 ADN_none ) ;
00395
00396
00397
00398 if( THD_is_file(DSET_HEADNAME(new_dset)) ){
00399 fprintf(stderr,"*** Fatal error: file %s already exists!\n",
00400 DSET_HEADNAME(new_dset) ) ;
00401 exit(1) ;
00402 }
00403
00404
00405
00406
00407 if( DT_verb ) fprintf(stderr,"+++ Loading input dataset .BRIK\n") ;
00408
00409 DSET_mallocize( new_dset ) ;
00410 DSET_mallocize( DT_dset ) ;
00411 DSET_load( DT_dset ) ;
00412 if( !DSET_LOADED(DT_dset) ){
00413 fprintf(stderr,"*** Can't read input dataset .BRIK!\n") ;
00414 exit(1) ;
00415 }
00416
00417 nvals = DSET_NVALS(new_dset) ;
00418 for( iv=0 ; iv < nvals ; iv++ )
00419 EDIT_substitute_brick( new_dset , iv ,
00420 DSET_BRICK_TYPE(DT_dset,iv) ,
00421 DSET_ARRAY(DT_dset,iv) ) ;
00422
00423 if( DT_norm && DSET_BRICK_TYPE(new_dset,0) != MRI_float ){
00424 fprintf(stderr,"+++ Warning: turning -normalize option off since\n"
00425 " input dataset is not in float format!\n");
00426 DT_norm = 0 ;
00427 }
00428
00429
00430
00431
00432 nvec = 0 ;
00433 for( ii=0 ; ii < IMARR_COUNT(DT_imar) ; ii++ )
00434 nvec += IMARR_SUBIMAGE(DT_imar,ii)->ny ;
00435
00436 refvec = (float **) malloc( sizeof(float *)*nvec ) ;
00437 for( kk=ii=0 ; ii < IMARR_COUNT(DT_imar) ; ii++ ){
00438 fv = MRI_FLOAT_PTR( IMARR_SUBIMAGE(DT_imar,ii) ) ;
00439 for( jj=0 ; jj < IMARR_SUBIMAGE(DT_imar,ii)->ny ; jj++ )
00440 refvec[kk++] = fv + ( jj * IMARR_SUBIMAGE(DT_imar,ii)->nx ) ;
00441 }
00442
00443 fit = (float *) malloc( sizeof(float) * nvals ) ;
00444
00445
00446
00447 if( !DT_byslice ){
00448 choleski = startup_lsqfit( nvals , NULL , nvec , refvec ) ;
00449 if( choleski == NULL ){
00450 fprintf(stderr,"*** Linearly dependent vectors can't be used!\n") ;
00451 exit(1) ;
00452 }
00453
00454
00455
00456 nvox = DSET_NVOX(new_dset) ;
00457
00458 if( DT_verb ) fprintf(stderr,"+++ Computing voxel fits\n") ;
00459
00460 for( kk=0 ; kk < nvox ; kk++ ){
00461
00462 flim = THD_extract_series( kk , new_dset , 0 ) ;
00463 fv = MRI_FLOAT_PTR(flim) ;
00464 fc = delayed_lsqfit( nvals, fv, nvec, refvec, choleski ) ;
00465
00466 for( ii=0 ; ii < nvals ; ii++ ) fit[ii] = 0.0 ;
00467
00468 for( jj=0 ; jj < nvec ; jj++ )
00469 for( ii=0 ; ii < nvals ; ii++ )
00470 fit[ii] += fc[jj] * refvec[jj][ii] ;
00471
00472 if( !DT_replace )
00473 for( ii=0 ; ii < nvals ; ii++ ) fit[ii] = fv[ii] - fit[ii] ;
00474
00475 if( DT_norm ) THD_normalize( nvals , fit ) ;
00476
00477 THD_insert_series( kk, new_dset, nvals, MRI_float, fit, 0 ) ;
00478
00479 free(fc) ; mri_free(flim) ;
00480 }
00481
00482 free(choleski) ;
00483
00484
00485
00486 }
00487 #ifdef ALLOW_BYSLICE
00488 else {
00489 int ksl , nslice , tt , nx,ny , nxy , kxy ;
00490 MRI_IMAGE * vim ;
00491
00492
00493
00494 for( kk=ii=0 ; ii < DT_nvector ; ii++ ){
00495 for( jj=0 ; jj < IMARR_SUBIMAGE(DT_imar,ii)->ny ; jj++ )
00496 refvec[kk++] = (float *) malloc( sizeof(float) * nvals ) ;
00497 }
00498
00499 nslice = DSET_NZ(new_dset) ;
00500 nxy = DSET_NX(new_dset) * DSET_NY(new_dset) ;
00501
00502
00503
00504 for( ksl=0 ; ksl < nslice ; ksl++ ){
00505
00506 if( DT_verb ) fprintf(stderr,"+++ Computing voxel fits for slice %d\n",ksl) ;
00507
00508
00509
00510 for( kk=ii=0 ; ii < DT_nvector ; ii++ ){
00511 vim = IMARR_SUBIMAGE(DT_imar,ii) ;
00512 nx = vim->nx ; ny = vim->ny ;
00513 for( jj=0 ; jj < ny ; jj++ ){
00514 fv = MRI_FLOAT_PTR(vim) + (jj*nx) ;
00515 for( tt=0 ; tt < nvals ; tt++ )
00516 refvec[kk][tt] = fv[ksl+tt*nslice] ;
00517 }
00518 }
00519
00520
00521
00522 choleski = startup_lsqfit( nvals , NULL , nvec , refvec ) ;
00523 if( choleski == NULL ){
00524 fprintf(stderr,"*** Linearly dependent vectors found at slice %d\n",ksl) ;
00525 exit(1) ;
00526 }
00527
00528
00529
00530 for( kxy=0 ; kxy < nxy ; kxy++ ){
00531
00532 kk = kxy + ksl*nxy ;
00533 flim = THD_extract_series( kk , new_dset , 0 ) ;
00534 fv = MRI_FLOAT_PTR(flim) ;
00535 fc = delayed_lsqfit( nvals, fv, nvec, refvec, choleski ) ;
00536
00537 for( ii=0 ; ii < nvals ; ii++ ) fit[ii] = 0.0 ;
00538
00539 for( jj=0 ; jj < nvec ; jj++ )
00540 for( ii=0 ; ii < nvals ; ii++ )
00541 fit[ii] += fc[jj] * refvec[jj][ii] ;
00542
00543 if( !DT_replace )
00544 for( ii=0 ; ii < nvals ; ii++ ) fit[ii] = fv[ii] - fit[ii] ;
00545
00546 if( DT_norm ) THD_normalize( nvals , fit ) ;
00547
00548 THD_insert_series( kk, new_dset, nvals, MRI_float, fit, 0 ) ;
00549
00550 free(fc) ; mri_free(flim) ;
00551 }
00552
00553 free(choleski) ;
00554
00555 }
00556
00557 }
00558 #endif
00559
00560
00561
00562 DSET_write(new_dset) ;
00563 if( DT_verb ) fprintf(stderr,"+++ Wrote output dataset: %s\n",DSET_BRIKNAME(new_dset)) ;
00564 exit(0) ;
00565 }
