00046 {
00047    int nopt = 1 , nvals , ii , nvcheck ;
00048 
00049    INIT_IMARR(DT_imar) ;
00050 
00051    while( nopt < argc && argv[nopt][0] == '-' ){
00052 
00053       /**** -prefix prefix ****/
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       /**** -session directory ****/
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       /**** -verb ****/
00076 
00077       if( strncmp(argv[nopt],"-verb",5) == 0 ){
00078          DT_verb++ ;
00079          nopt++ ; continue ;
00080       }
00081 
00082       /**** -replace ****/
00083 
00084       if( strncmp(argv[nopt],"-replace",5) == 0 ){
00085          DT_replace++ ;
00086          nopt++ ; continue ;
00087       }
00088 
00089 #ifdef ALLOW_BYSLICE
00090       /**** -byslice [08 Dec 1999] ****/
00091 
00092       if( strncmp(argv[nopt],"-byslice",5) == 0 ){
00093          DT_byslice++ ;
00094          nopt++ ; continue ;
00095       }
00096 #endif
00097 
00098       /**** -normalize [23 Nov 1999] ****/
00099 
00100       if( strncmp(argv[nopt],"-normalize",5) == 0 ){
00101          DT_norm++ ;
00102          nopt++ ; continue ;
00103       }
00104 
00105       /**** -vector ****/
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       /**** -del ****/
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       /**** -expr ****/
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 ){   /* initialize storage */
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] ;                         /* string */
00165          DT_exdel[DT_exnum]  = DT_current_del ;                     /* delta */
00166          DT_excode[DT_exnum] = PARSER_generate_code( argv[nopt] ) ; /* compile */
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 ;                       /* find symbol */
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       /**** ERROR ****/
00191 
00192       fprintf(stderr,"*** Unknown option: %s\n",argv[nopt]) ; exit(1) ;
00193 
00194    }  /* end of scan over options */
00195 
00196    /*-- check for errors --*/
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    /*--- read input dataset ---*/
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    /*-- check vectors for good size --*/
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    /*--- create time series from expressions */
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 }

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 }

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    /*** read input options ***/
00369 
00370    if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) DT_Syntax() ;
00371 
00372    /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
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    /*** create new dataset (empty) ***/
00383 
00384    new_dset = EDIT_empty_copy( DT_dset ) ; /* make a copy of its header */
00385 
00386    /* modify its header */
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    /* can't re-write existing dataset */
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    /* read input in, and attach its bricks to the output dataset */
00405    /* (not good in a plugin, but OK in a standalone program!)    */
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    /* load reference (detrending) vectors;
00430       setup to do least squares fitting of each voxel */
00431 
00432    nvec = 0 ;
00433    for( ii=0 ; ii < IMARR_COUNT(DT_imar) ; ii++ )  /* number of detrending vectors */
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++ )          /* compute ptr */
00440          refvec[kk++] = fv + ( jj * IMARR_SUBIMAGE(DT_imar,ii)->nx ) ;  /* to vectors  */
00441    }
00442 
00443    fit = (float *) malloc( sizeof(float) * nvals ) ;  /* will get fit to voxel data */
00444 
00445    /*--- do the all-voxels-together case ---*/
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       /* loop over voxels, fitting and detrending (or replacing) */
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 ) ;              /* data */
00463          fv   = MRI_FLOAT_PTR(flim) ;
00464          fc   = delayed_lsqfit( nvals, fv, nvec, refvec, choleski ) ;  /* coef */
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] ;                    /* fit */
00471 
00472          if( !DT_replace )                                             /* remove */
00473             for( ii=0 ; ii < nvals ; ii++ ) fit[ii] = fv[ii] - fit[ii] ;
00474 
00475          if( DT_norm ) THD_normalize( nvals , fit ) ;  /* 23 Nov 1999 */
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       /*- end of all-voxels-together case -*/
00485 
00486    }
00487 #ifdef ALLOW_BYSLICE
00488      else {                                 /*- start of slice case [08 Dec 1999] -*/
00489       int ksl , nslice , tt , nx,ny , nxy , kxy ;
00490       MRI_IMAGE * vim ;
00491 
00492       /* make separate space for the slice-wise detrending vectors */
00493 
00494       for( kk=ii=0 ; ii < DT_nvector ; ii++ ){
00495          for( jj=0 ; jj < IMARR_SUBIMAGE(DT_imar,ii)->ny ; jj++ )       /* replace ptrs */
00496             refvec[kk++] = (float *) malloc( sizeof(float) * nvals ) ;  /* to vectors   */
00497       }
00498 
00499       nslice = DSET_NZ(new_dset) ;
00500       nxy    = DSET_NX(new_dset) * DSET_NY(new_dset) ;
00501 
00502       /* loop over slices */
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          /* extract slice vectors from input interlaced vectors */
00509 
00510          for( kk=ii=0 ; ii < DT_nvector ; ii++ ){        /* loop over vectors */
00511             vim = IMARR_SUBIMAGE(DT_imar,ii) ;           /* ii-th vector image */
00512             nx = vim->nx ; ny = vim->ny ;                /* dimensions */
00513             for( jj=0 ; jj < ny ; jj++ ){                /* loop over columns */
00514                fv = MRI_FLOAT_PTR(vim) + (jj*nx) ;       /* ptr to column */
00515                for( tt=0 ; tt < nvals ; tt++ )           /* loop over time */
00516                   refvec[kk][tt] = fv[ksl+tt*nslice] ;   /* data point */
00517             }
00518          }
00519 
00520          /* initialize fitting for this slice */
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          /* loop over voxels in this slice */
00529 
00530          for( kxy=0 ; kxy < nxy ; kxy++ ){
00531 
00532             kk   = kxy + ksl*nxy ;                                        /* 3D index */
00533             flim = THD_extract_series( kk , new_dset , 0 ) ;              /* data */
00534             fv   = MRI_FLOAT_PTR(flim) ;
00535             fc   = delayed_lsqfit( nvals, fv, nvec, refvec, choleski ) ;  /* coef */
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] ;                    /* fit */
00542 
00543             if( !DT_replace )                                             /* remove */
00544                for( ii=0 ; ii < nvals ; ii++ ) fit[ii] = fv[ii] - fit[ii] ;
00545 
00546             if( DT_norm ) THD_normalize( nvals , fit ) ;  /* 23 Nov 1999 */
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       } /* end of loop over slices */
00556 
00557    } /*- end of -byslice case -*/
00558 #endif
00559 
00560    /*-- done done done done --*/
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 }