00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011
00012
00013
00014 #ifndef myXtFree
00015 # define myXtFree(xp) (XtFree((char *)(xp)) , (xp)=NULL)
00016 #endif
00017
00018
00019
00020 static THD_3dim_dataset_array *TCAT_dsar = NULL ;
00021 static XtPointer_array *TCAT_subv = NULL ;
00022 static int TCAT_nvox = -1 ;
00023 static int TCAT_dry = 0 ;
00024 static int TCAT_verb = 0 ;
00025 static int TCAT_type = -1 ;
00026 static int TCAT_glue = 0 ;
00027 static int TCAT_rlt = 0 ;
00028
00029 static int TCAT_rqt = 0 ;
00030 static int TCAT_rct = 0 ;
00031
00032 static char TCAT_output_prefix[THD_MAX_PREFIX] = "tcat" ;
00033 static char TCAT_session[THD_MAX_NAME] = "./" ;
00034
00035
00036
00037
00038
00039
00040 #define DSUB(id) DSET_IN_3DARR(TCAT_dsar,(id))
00041 #define NSUBV(id) ( ((int *)TCAT_subv->ar[(id)])[0] )
00042 #define SUBV(id,jj) ( ((int *)TCAT_subv->ar[(id)])[(jj)+1] )
00043
00044
00045
00046 void TCAT_read_opts( int , char ** ) ;
00047 void TCAT_Syntax(void) ;
00048 int * TCAT_get_subv( int , char * ) ;
00049
00050
00051
00052
00053
00054 void TCAT_read_opts( int argc , char *argv[] )
00055 {
00056 int nopt = 1 , ii ;
00057 char dname[THD_MAX_NAME] ;
00058 char subv[THD_MAX_NAME] ;
00059 char * cpt ;
00060 THD_3dim_dataset * dset ;
00061 int * svar ;
00062 char * str;
00063 int ok, ilen, nlen , max_nsub=0 ;
00064
00065 INIT_3DARR(TCAT_dsar) ;
00066 INIT_XTARR(TCAT_subv) ;
00067
00068 while( nopt < argc ){
00069
00070
00071
00072 if( strncmp(argv[nopt],"-prefix",6) == 0 ||
00073 strncmp(argv[nopt],"-output",6) == 0 ){
00074 if(TCAT_glue){
00075 fprintf(stderr,"*** -prefix and -glueto options are not compatible\n");
00076 exit(1) ;
00077 }
00078 nopt++ ;
00079 if( nopt >= argc ){
00080 fprintf(stderr,"*** need argument after -prefix!\n") ; exit(1) ;
00081 }
00082 MCW_strncpy( TCAT_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00083 continue ;
00084 }
00085
00086
00087
00088 if( strncmp(argv[nopt],"-session",6) == 0 ){
00089 if(TCAT_glue){
00090 fprintf(stderr,
00091 "*** -session and -glueto options are not compatible\n");
00092 exit(1) ;
00093 }
00094 nopt++ ;
00095 if( nopt >= argc ){
00096 fprintf(stderr,"*** need argument after -session!\n") ; exit(1) ;
00097 }
00098 MCW_strncpy( TCAT_session , argv[nopt++] , THD_MAX_NAME ) ;
00099 continue ;
00100 }
00101
00102
00103
00104 if( strncmp(argv[nopt],"-dry",3) == 0 ){
00105 TCAT_dry = 1 ; TCAT_verb++ ;
00106 nopt++ ; continue ;
00107 }
00108
00109
00110
00111 if( strncmp(argv[nopt],"-verb",5) == 0 ){
00112 TCAT_verb++ ;
00113 nopt++ ; continue ;
00114 }
00115
00116
00117
00118 if( strcmp(argv[nopt],"-rlt") == 0 ){
00119 TCAT_rlt = 1 ;
00120 nopt++ ; continue ;
00121 }
00122
00123
00124
00125 if( strcmp(argv[nopt],"-rlt+") == 0 ){
00126 TCAT_rlt = 2 ;
00127 nopt++ ; continue ;
00128 }
00129
00130
00131
00132 if( strcmp(argv[nopt],"-rlt++") == 0 ){
00133 TCAT_rlt = 3 ;
00134 nopt++ ; continue ;
00135 }
00136
00137
00138
00139 if( strcmp(argv[nopt],"-rqt") == 0 ){
00140 TCAT_rqt = 1 ;
00141 nopt++ ; continue ;
00142 }
00143
00144
00145
00146 if( strcmp(argv[nopt],"-rct") == 0 ){
00147 TCAT_rct = 1 ;
00148 nopt++ ; continue ;
00149 }
00150
00151
00152
00153 if( strncmp(argv[nopt],"-glueto",5) == 0 ){
00154 if( strncmp(TCAT_output_prefix, "tcat", 5) != 0 ){
00155 fprintf(stderr,"*** -prefix and -glueto options are not compatible\n");
00156 exit(1) ;
00157 }
00158 if( strncmp(TCAT_session, "./", 5) != 0 ){
00159 fprintf(stderr,
00160 "*** -session and -glueto options are not compatible\n");
00161 exit(1) ;
00162 }
00163 TCAT_glue = 1 ;
00164 nopt++ ;
00165 if( nopt >= argc ){
00166 fprintf(stderr,"*** need argument after -glueto!\n") ; exit(1) ;
00167 }
00168
00169
00170 ok = 1;
00171 nlen = strlen(argv[nopt]);
00172 if (nlen <= 5) ok = 0;
00173
00174 #define BACKASS
00175
00176 if (ok)
00177 {
00178 #ifndef BACKASS
00179 for (ilen = 0; ilen < nlen; ilen++)
00180 #else
00181 for( ilen=nlen-3 ; ilen >= 0 ; ilen-- )
00182 #endif
00183 {
00184 str = argv[nopt] + ilen;
00185 if (str[0] == '+') break;
00186 }
00187 #ifndef BACKASS
00188 if (ilen == nlen) ok = 0;
00189 #else
00190 if (ilen <= 0 ) ok = 0;
00191 #endif
00192 }
00193
00194 if (ok)
00195 {
00196 str = argv[nopt] + ilen + 1;
00197
00198 for (ii=FIRST_VIEW_TYPE ; ii <= LAST_VIEW_TYPE ; ii++)
00199 if (! strncmp(str,VIEW_codestr[ii],4)) break ;
00200
00201 if( ii > LAST_VIEW_TYPE ) ok = 0;
00202 }
00203
00204 if (! ok)
00205 {
00206 fprintf(stderr,
00207 "*** File name must end in +orig, +acpc, or +tlrc after -glueto\n");
00208 exit(1);
00209 }
00210
00211
00212 MCW_strncpy( TCAT_output_prefix , argv[nopt] , ilen+1) ;
00213
00214
00215
00216 }
00217
00218 if( argv[nopt][0] == '-' ){
00219 fprintf(stderr,"*** Unknown option: %s\n",argv[nopt]) ; exit(1) ;
00220 }
00221
00222
00223
00224 cpt = strstr(argv[nopt],"[") ;
00225
00226 if( cpt == NULL ){
00227 strcpy(dname,argv[nopt]) ;
00228 subv[0] = '\0' ;
00229 } else if( cpt == argv[nopt] ){
00230 fprintf(stderr,"*** Illegal dataset specifier: %s\n",argv[nopt]) ;
00231 exit(1) ;
00232 } else {
00233 ii = cpt - argv[nopt] ;
00234 memcpy(dname,argv[nopt],ii) ; dname[ii] = '\0' ;
00235 strcpy(subv,cpt) ;
00236 }
00237 nopt++ ;
00238
00239 dset = THD_open_one_dataset( dname ) ;
00240 if( dset == NULL ){
00241 fprintf(stderr,"*** Can't open dataset %s\n",dname) ; exit(1) ;
00242 }
00243 THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ;
00244
00245 if( TCAT_type < 0 ) TCAT_type = dset->type ;
00246
00247
00248
00249 ii = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz ;
00250 if( TCAT_nvox < 0 ){
00251 TCAT_nvox = ii ;
00252 } else if( ii != TCAT_nvox ){
00253 fprintf(stderr,"*** Dataset %s differs in size from others\n",dname);
00254 exit(1) ;
00255 }
00256 ADDTO_3DARR(TCAT_dsar,dset) ;
00257
00258
00259
00260
00261
00262
00263 svar = TCAT_get_subv( DSET_NVALS(dset) , subv ) ;
00264 if( svar == NULL || svar[0] <= 0 ){
00265 fprintf(stderr,"*** Can't decipher index codes from %s%s\n",dname,subv) ;
00266 exit(1) ;
00267 }
00268 ADDTO_XTARR(TCAT_subv,svar) ;
00269
00270 max_nsub = MAX( max_nsub , svar[0] ) ;
00271
00272 if( TCAT_rlt == 3 && svar[0] < 3 )
00273 fprintf(stderr,
00274 "*** Warning: -rlt++ option won't work properly with\n"
00275 " less than 3 sub-bricks per input dataset!\n") ;
00276
00277 }
00278
00279
00280
00281 if( max_nsub < 3 && TCAT_rlt ){
00282 fprintf(stderr,"*** Warning: can't apply -rlt option -- "
00283 "Not enough points per input dataset.\n" ) ;
00284 TCAT_rlt = 0 ;
00285 }
00286
00287 if( TCAT_rlt && TCAT_dry ){
00288 fprintf(stderr,"*** Warning: -rlt option does nothing with -dry!\n") ;
00289 TCAT_rlt = 0 ;
00290 }
00291
00292 return ;
00293 }
00294
00295
00296
00297
00298
00299 int * TCAT_get_subv( int nvals , char *str )
00300 {
00301 int * subv = NULL ;
00302 int ii , ipos , nout , slen ;
00303 int ibot,itop,istep , nused ;
00304 char * cpt ;
00305
00306
00307
00308 if( nvals < 1 ) return NULL ;
00309
00310
00311
00312 if( str == NULL || str[0] == '\0' ){
00313 subv = (int *) XtMalloc( sizeof(int) * (nvals+1) ) ;
00314 subv[0] = nvals ;
00315 for( ii=0 ; ii < nvals ; ii++ ) subv[ii+1] = ii ;
00316 return subv ;
00317 }
00318
00319
00320
00321 subv = (int *) XtMalloc( sizeof(int) * 2 ) ;
00322 subv[0] = nout = 0 ;
00323
00324 ipos = 0 ;
00325 if( str[ipos] == '[' ) ipos++ ;
00326
00327
00328
00329 slen = strlen(str) ;
00330 while( ipos < slen && str[ipos] != ']' ){
00331
00332
00333
00334 if( str[ipos] == '$' ){
00335 ibot = nvals-1 ; ipos++ ;
00336 } else {
00337 ibot = strtol( str+ipos , &cpt , 10 ) ;
00338 if( ibot < 0 ){ myXtFree(subv) ; return NULL ; }
00339 if( ibot >= nvals ) ibot = nvals-1 ;
00340 nused = (cpt-(str+ipos)) ;
00341 if( ibot == 0 && nused == 0 ){ myXtFree(subv) ; return NULL ; }
00342 ipos += nused ;
00343 }
00344
00345
00346
00347 if( str[ipos] == ',' || str[ipos] == '\0' || str[ipos] == ']' ){
00348 nout++ ;
00349 subv = (int *) XtRealloc( (char *)subv , sizeof(int) * (nout+1) ) ;
00350 subv[0] = nout ;
00351 subv[nout] = ibot ;
00352 ipos++ ; continue ;
00353 }
00354
00355
00356
00357 if( str[ipos] == '-' ){
00358 ipos++ ;
00359 } else if( str[ipos] == '.' && str[ipos+1] == '.' ){
00360 ipos++ ; ipos++ ;
00361 } else {
00362 myXtFree(subv) ; return NULL ;
00363 }
00364
00365
00366
00367 if( str[ipos] == '$' ){
00368 itop = nvals-1 ; ipos++ ;
00369 } else {
00370 itop = strtol( str+ipos , &cpt , 10 ) ;
00371 if( itop < 0 ){ myXtFree(subv) ; return NULL ; }
00372 if( itop >= nvals ) itop = nvals-1 ;
00373 nused = (cpt-(str+ipos)) ;
00374 if( itop == 0 && nused == 0 ){ myXtFree(subv) ; return NULL ; }
00375 ipos += nused ;
00376 }
00377
00378
00379
00380 istep = (ibot <= itop) ? 1 : -1 ;
00381
00382
00383
00384 if( str[ipos] == '(' ){
00385 ipos++ ;
00386 istep = strtol( str+ipos , &cpt , 10 ) ;
00387 if( istep == 0 ){ myXtFree(subv) ; return NULL ; }
00388 nused = (cpt-(str+ipos)) ;
00389 ipos += nused ;
00390 if( str[ipos] == ')' ) ipos++ ;
00391 }
00392
00393
00394
00395 for( ii=ibot ; (ii-itop)*istep <= 0 ; ii += istep ){
00396 nout++ ;
00397 subv = (int *) XtRealloc( (char *)subv , sizeof(int) * (nout+1) ) ;
00398 subv[0] = nout ;
00399 subv[nout] = ii ;
00400 }
00401
00402
00403
00404 if( str[ipos] == ',' ) ipos++ ;
00405
00406 }
00407
00408 return subv ;
00409 }
00410
00411
00412
00413 void TCAT_Syntax(void)
00414 {
00415 printf(
00416 "Concatenate sub-bricks from input datasets into one big 3D+time dataset.\n"
00417 "Usage: 3dTcat options\n"
00418 "where the options are:\n"
00419 ) ;
00420
00421 printf(
00422 " -prefix pname = Use 'pname' for the output dataset prefix name.\n"
00423 " OR -output pname [default='tcat']\n"
00424 "\n"
00425 " -session dir = Use 'dir' for the output dataset session directory.\n"
00426 " [default='./'=current working directory]\n"
00427 " -glueto fname = Append bricks to the end of the 'fname' dataset.\n"
00428 " This command is an alternative to the -prefix \n"
00429 " and -session commands. \n"
00430 " -dry = Execute a 'dry run'; that is, only print out\n"
00431 " what would be done. This is useful when\n"
00432 " combining sub-bricks from multiple inputs.\n"
00433 " -verb = Print out some verbose output as the program\n"
00434 " proceeds (-dry implies -verb).\n"
00435 " Using -verb twice results in quite lengthy output.\n"
00436 " -rlt = Remove linear trends in each voxel time series loaded\n"
00437 " from each input dataset, SEPARATELY. That is, the\n"
00438 " data from each dataset is detrended separately.\n"
00439 " At least 3 sub-bricks from a dataset must be input\n"
00440 " for this option to apply.\n"
00441 " Notes: (1) -rlt removes the least squares fit of 'a+b*t'\n"
00442 " to each voxel time series; this means that\n"
00443 " the mean is removed as well as the trend.\n"
00444 " This effect makes it impractical to compute\n"
00445 " the %% Change using AFNI's internal FIM.\n"
00446 " (2) To have the mean of each dataset time series added\n"
00447 " back in, use this option in the form '-rlt+'.\n"
00448 " In this case, only the slope 'b*t' is removed.\n"
00449 " (3) To have the overall mean of all dataset time\n"
00450 " series added back in, use this option in the\n"
00451 " form '-rlt++'. In this case, 'a+b*t' is removed\n"
00452 " from each input dataset separately, and the\n"
00453 " mean of all input datasets is added back in at\n"
00454 " the end. (This option will work properly only\n"
00455 " if all input datasets use at least 3 sub-bricks!)\n"
00456 " (4) -rlt can be used on datasets that contain shorts\n"
00457 " or floats, but not on complex- or byte-valued\n"
00458 " datasets.\n"
00459 "\n"
00460 "Command line arguments after the above are taken as input datasets.\n"
00461 "A dataset is specified using one of these forms:\n"
00462 " 'prefix+view', 'prefix+view.HEAD', or 'prefix+view.BRIK'.\n"
00463 "\n"
00464 "SUB-BRICK SELECTION:\n"
00465 "You can also add a sub-brick selection list after the end of the\n"
00466 "dataset name. This allows only a subset of the sub-bricks to be\n"
00467 "included into the output (by default, all of the input dataset\n"
00468 "is copied into the output). A sub-brick selection list looks like\n"
00469 "one of the following forms:\n"
00470 " fred+orig[5] ==> use only sub-brick #5\n"
00471 " fred+orig[5,9,17] ==> use #5, #9, and #12\n"
00472 " fred+orig[5..8] or [5-8] ==> use #5, #6, #7, and #8\n"
00473 " fred+orig[5..13(2)] or [5-13(2)] ==> use #5, #7, #9, #11, and #13\n"
00474 "Sub-brick indexes start at 0. You can use the character '$'\n"
00475 "to indicate the last sub-brick in a dataset; for example, you\n"
00476 "can select every third sub-brick by using the selection list\n"
00477 " fred+orig[0..$(3)]\n"
00478 "\n"
00479 "NOTES:\n"
00480 "* The TR and other time-axis properties are taken from the\n"
00481 " first input dataset that is itself 3D+time. If no input\n"
00482 " datasets contain such information, then TR is set to 1.0.\n"
00483 " This can be altered using the 3drefit program.\n"
00484 "\n"
00485 "* The sub-bricks are output in the order specified, which may\n"
00486 " not be the order in the original datasets. For example, using\n"
00487 " fred+orig[0..$(2),1..$(2)]\n"
00488 " will cause the sub-bricks in fred+orig to be output into the\n"
00489 " new dataset in an interleaved fashion. Using\n"
00490 " fred+orig[$..0]\n"
00491 " will reverse the order of the sub-bricks in the output.\n"
00492 " If the -rlt option is used, the sub-bricks selected from each\n"
00493 " input dataset will be re-ordered into the output dataset, and\n"
00494 " then this sequence will be detrended.\n"
00495 "\n"
00496 "* You can use the '3dinfo' program to see how many sub-bricks\n"
00497 " a 3D+time or a bucket dataset contains.\n"
00498 "\n"
00499 "* The '$', '(', ')', '[', and ']' characters are special to\n"
00500 " the shell, so you will have to escape them. This is most easily\n"
00501 " done by putting the entire dataset plus selection list inside\n"
00502 " single quotes, as in 'fred+orig[5..7,9]'.\n"
00503 "\n"
00504 "* You may wish to use the 3drefit program on the output dataset\n"
00505 " to modify some of the .HEAD file parameters.\n"
00506 ) ;
00507
00508 exit(0) ;
00509 }
00510
00511
00512
00513 int main( int argc , char * argv[] )
00514 {
00515 int ninp , ids , nv , iv,jv,kv , ivout , new_nvals , ivbot,ivtop ;
00516 THD_3dim_dataset * new_dset=NULL , * dset ;
00517 char buf[256] ;
00518 float * rlt0=NULL , *rlt1=NULL ;
00519 float *rltsum=NULL ;
00520 int nrltsum ;
00521
00522
00523
00524 if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) TCAT_Syntax() ;
00525
00526
00527
00528 mainENTRY("3dTcat main"); machdep() ; PRINT_VERSION("3dTcat") ;
00529
00530 { int new_argc ; char ** new_argv ;
00531 addto_args( argc , argv , &new_argc , &new_argv ) ;
00532 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00533 }
00534
00535 AFNI_logger("3dTcat",argc,argv) ;
00536
00537 TCAT_read_opts( argc , argv ) ;
00538
00539
00540
00541 ninp = TCAT_dsar->num ;
00542 if( ninp < 1 ){
00543 fprintf(stderr,"*** No input datasets?\n") ; exit(1) ;
00544 }
00545
00546 new_nvals = 0 ;
00547 for( ids=0 ; ids < ninp ; ids++ ) new_nvals += NSUBV(ids) ;
00548
00549 if( new_nvals < 2 ){
00550 fprintf(stderr,
00551 "*** Can't create 3D+time dataset with only %d sub-bricks!\n",
00552 new_nvals) ;
00553 exit(1) ;
00554 }
00555
00556 if( TCAT_verb ) printf("-verb: output will have %d sub-bricks\n",new_nvals) ;
00557
00558
00559
00560 for( ids=0 ; ids < ninp ; ids++ ){
00561 dset = DSUB(ids) ;
00562 if( DSET_TIMESTEP(dset) > 0.0 ) break ;
00563 }
00564 if( ids == ninp ){ ids = 0 ; dset = DSUB(0) ; }
00565
00566 new_dset = EDIT_empty_copy( dset ) ;
00567
00568
00569
00570 for( iv=0 ; iv < ninp ; iv++ ){
00571 if( iv != ids && !EQUIV_DATAXES(new_dset->daxes,DSUB(iv)->daxes) )
00572 fprintf(stderr,"++ WARNING: %s grid mismatch with %s\n",
00573 DSET_BRIKNAME(dset) , DSET_BRIKNAME(DSUB(iv)) ) ;
00574 }
00575
00576 tross_Make_History( "3dTcat" , argc,argv , new_dset ) ;
00577
00578
00579
00580 EDIT_dset_items( new_dset ,
00581 ADN_prefix , TCAT_output_prefix ,
00582 ADN_directory_name, TCAT_session ,
00583 ADN_type , TCAT_type ,
00584 ADN_func_type , ISANATTYPE(TCAT_type) ? ANAT_EPI_TYPE
00585 : FUNC_FIM_TYPE ,
00586 ADN_ntt , new_nvals ,
00587 ADN_nvals , new_nvals ,
00588 ADN_none ) ;
00589
00590
00591
00592 if( DSET_TIMESTEP(new_dset) <= 0.0 ){
00593 float TR = 1.0 ;
00594 float torg = 0.0 , tdur = 0.0 ;
00595 int tunits = UNITS_SEC_TYPE ;
00596
00597 #if 0
00598 for( ids=0 ; ids < ninp ; ids++ ){
00599 dset = DSUB(ids) ;
00600 if( DSET_TIMESTEP(dset) > 0.0 ){
00601 TR = DSET_TIMESTEP(dset) ; tunits = DSET_TIMEUNITS(dset) ;
00602 torg = DSET_TIMEORIGIN(dset) ; tdur = DSET_TIMEDURATION(dset) ;
00603 break ;
00604 }
00605 }
00606 #endif
00607
00608 EDIT_dset_items( new_dset ,
00609 ADN_tunits , tunits ,
00610 ADN_ttdel , TR ,
00611 ADN_ttorg , torg ,
00612 ADN_ttdur , tdur ,
00613 ADN_none ) ;
00614 }
00615
00616
00617
00618 if (! TCAT_glue){
00619 if( THD_is_file(DSET_HEADNAME(new_dset)) ){
00620 fprintf(stderr,"*** Fatal error: file %s already exists!\n",
00621 DSET_HEADNAME(new_dset) ) ;
00622 exit(1) ;
00623 }
00624 } else {
00625
00626
00627 new_dset->idcode = DSUB(0) -> idcode ;
00628 }
00629
00630 THD_force_malloc_type( new_dset->dblk , DATABLOCK_MEM_MALLOC ) ;
00631
00632
00633
00634 if( TCAT_rlt ){
00635 rlt0 = (float *) malloc( sizeof(float) * TCAT_nvox ) ;
00636 rlt1 = (float *) malloc( sizeof(float) * TCAT_nvox ) ;
00637 if( rlt0 == NULL || rlt1 == NULL ){
00638 fprintf(stderr,"*** Error: can't malloc memory for detrending!\n") ;
00639 exit(1) ;
00640 }
00641
00642 if( TCAT_rlt == 3 ){
00643 rltsum = (float *) malloc( sizeof(float) * TCAT_nvox ) ;
00644 if( rltsum == NULL ){
00645 fprintf(stderr,"*** Error: can't malloc memory for detrending!\n") ;
00646 exit(1) ;
00647 }
00648 for( iv=0 ; iv < TCAT_nvox ; iv++ ) rltsum[iv] = 0.0 ;
00649 nrltsum = 0 ;
00650 }
00651 }
00652
00653
00654
00655 if( ninp > 1 ) myXtFree( new_dset->keywords ) ;
00656
00657 ivout = 0 ;
00658 for( ids=0 ; ids < ninp ; ids++ ){
00659 dset = DSUB(ids) ;
00660 nv = NSUBV(ids) ;
00661
00662 if( ! TCAT_dry ){
00663 if( TCAT_verb ) printf("-verb: loading %s\n",DSET_FILECODE(dset)) ;
00664 DSET_load(dset) ;
00665 if( ! DSET_LOADED(dset) ){
00666 fprintf(stderr,"*** Fatal error: can't load data from %s\n",
00667 DSET_FILECODE(dset)) ;
00668 exit(1) ;
00669 }
00670 }
00671
00672
00673
00674 ivbot = ivout ;
00675 for( iv=0 ; iv < nv ; iv++ ){
00676 jv = SUBV(ids,iv) ;
00677
00678 if( ! TCAT_dry ){
00679 EDIT_substitute_brick( new_dset , ivout ,
00680 DSET_BRICK_TYPE(dset,jv) , DSET_ARRAY(dset,jv) ) ;
00681
00682
00683
00684
00685 if( ISBUCKET(dset) )
00686 sprintf (buf, "%s", DSET_BRICK_LABEL(dset,jv));
00687 else
00688 sprintf(buf,"%.12s[%d]",DSET_PREFIX(dset),jv) ;
00689 EDIT_dset_items( new_dset, ADN_brick_label_one+ivout, buf, ADN_none );
00690
00691 sprintf(buf,"%s[%d]",DSET_FILECODE(dset),jv) ;
00692 EDIT_dset_items(
00693 new_dset, ADN_brick_keywords_replace_one+ivout, buf, ADN_none ) ;
00694
00695 EDIT_dset_items(
00696 new_dset ,
00697 ADN_brick_fac_one +ivout, DSET_BRICK_FACTOR(dset,jv),
00698 ADN_brick_keywords_append_one+ivout, DSET_BRICK_KEYWORDS(dset,jv),
00699 ADN_none ) ;
00700
00701
00702
00703 kv = DSET_BRICK_STATCODE(dset,jv) ;
00704
00705 if( FUNC_IS_STAT(kv) ){
00706
00707 int npar = FUNC_need_stat_aux[kv] , lv ;
00708 float * par = (float *) malloc( sizeof(float) * (npar+2) ) ;
00709 float * sax = DSET_BRICK_STATAUX(dset,jv) ;
00710 par[0] = kv ;
00711 par[1] = npar ;
00712 for( lv=0 ; lv < npar ; lv++ )
00713 par[lv+2] = (sax != NULL) ? sax[lv] : 0.0 ;
00714
00715 EDIT_dset_items(new_dset ,
00716 ADN_brick_stataux_one+ivout , par ,
00717 ADN_none ) ;
00718 free(par) ;
00719 #if 0
00720
00721
00722 } else if( ISFUNC(dset) &&
00723 FUNC_IS_STAT(dset->func_type) &&
00724 jv == FUNC_ival_thr[dset->func_type] ){
00725
00726 int npar , lv ;
00727 float * par , * sax ;
00728 kv = dset->func_type ;
00729 npar = FUNC_need_stat_aux[kv] ;
00730 par = (float *) malloc( sizeof(float) * (npar+2) ) ;
00731 sax = dset->stat_aux ;
00732 par[0] = kv ;
00733 par[1] = npar ;
00734 for( lv=0 ; lv < npar ; lv++ )
00735 par[lv+2] = (sax != NULL) ? sax[lv] : 0.0 ;
00736
00737 EDIT_dset_items(new_dset ,
00738 ADN_brick_stataux_one+ivout , par ,
00739 ADN_none ) ;
00740 free(par) ;
00741 #endif
00742 }
00743
00744
00745
00746 if( TCAT_verb > 1 ) printf("-verb: copied %s[%d] into %s[%d]\n" ,
00747 DSET_FILECODE(dset) , jv ,
00748 DSET_FILECODE(new_dset) , ivout ) ;
00749 } else {
00750 printf("-dry: would copy %s[%d] into %s[%d]\n" ,
00751 DSET_FILECODE(dset) , jv ,
00752 DSET_FILECODE(new_dset) , ivout ) ;
00753 }
00754
00755 ivout++ ;
00756 }
00757 ivtop = ivout ;
00758
00759
00760
00761
00762
00763 if( ! TCAT_dry && nv < DSET_NVALS(dset) ){
00764
00765 for( kv=0 ; kv < DSET_NVALS(dset) ; kv++ ){
00766 for( iv=0 ; iv < nv ; iv++ ){
00767 jv = SUBV(ids,iv) ;
00768 if( jv == kv ) break ;
00769 }
00770 if( iv == nv ){
00771 mri_free( DSET_BRICK(dset,kv) ) ;
00772 #if 0
00773 if( TCAT_verb ) printf("-verb: unloaded unused %s[%d]\n" ,
00774 DSET_FILECODE(dset) , kv ) ;
00775 #endif
00776 }
00777 }
00778 }
00779
00780
00781
00782 if( TCAT_rlt ){
00783
00784
00785
00786 if( ivtop-ivbot < 3 ){
00787 if( TCAT_verb )
00788 printf("-verb: skipping -rlt for %s\n",DSET_FILECODE(dset)) ;
00789
00790 } else {
00791 float c0,c1,c2 , det , a0,a1,a2 , qq ;
00792 int iv , ns , kk , err=0 ;
00793
00794 if( TCAT_verb )
00795 printf("-verb: applying -rlt to data from %s\n",DSET_FILECODE(dset)) ;
00796
00797
00798
00799 ns = ivtop - ivbot ;
00800 c0 = ns ;
00801 c1 = 0.5 * ns * (ns-1) ;
00802 c2 = 0.16666667 * ns * (ns-1) * (2*ns-1) ;
00803 det = c0*c2 - c1*c1 ;
00804 a0 = c2 / det ;
00805 a1 = -c1 / det ;
00806 a2 = c0 / det ;
00807
00808
00809
00810 for( iv=0 ; iv < TCAT_nvox ; iv++ ) rlt0[iv] = rlt1[iv] = 0.0 ;
00811
00812
00813
00814 for( kk=ivbot ; kk < ivtop ; kk++ ){
00815 qq = kk - ivbot ;
00816 switch( DSET_BRICK_TYPE(new_dset,kk) ){
00817 default:
00818 err = 1 ;
00819 fprintf(stderr,
00820 "*** Warning: -rlt can't use datum type %s from %s\n",
00821 MRI_TYPE_name[DSET_BRICK_TYPE(new_dset,kk)] ,
00822 DSET_FILECODE(dset) ) ;
00823 break ;
00824
00825 case MRI_short:{
00826 short * bar = (short *) DSET_ARRAY(new_dset,kk) ;
00827 float fac = DSET_BRICK_FACTOR(new_dset,kk) ;
00828
00829 if( fac == 0.0 ) fac = 1.0 ;
00830 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
00831 rlt0[iv] += (fac * bar[iv]) ;
00832 rlt1[iv] += (fac * bar[iv]) * qq ;
00833 }
00834 }
00835 break ;
00836
00837 case MRI_float:{
00838 float * bar = (float *) DSET_ARRAY(new_dset,kk) ;
00839 float fac = DSET_BRICK_FACTOR(new_dset,kk) ;
00840
00841 if( fac == 0.0 ) fac = 1.0 ;
00842 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
00843 rlt0[iv] += (fac * bar[iv]) ;
00844 rlt1[iv] += (fac * bar[iv]) * qq ;
00845 }
00846 }
00847 break ;
00848 }
00849 if( err ) break ;
00850 }
00851
00852
00853
00854 if( !err ){
00855 float qmid = 0.0 ;
00856
00857 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
00858 c0 = a0 * rlt0[iv] + a1 * rlt1[iv] ;
00859 c1 = a1 * rlt0[iv] + a2 * rlt1[iv] ;
00860 rlt0[iv] = c0 ; rlt1[iv] = c1 ;
00861 }
00862
00863 if( TCAT_rlt == 2 ){
00864 qmid = 0.5 * (ns-1) ;
00865 for( iv=0 ; iv < TCAT_nvox ; iv++ ) rlt0[iv] = 0.0 ;
00866 } else if( TCAT_rlt == 3 ){
00867 nrltsum += ns ;
00868 for( iv=0 ; iv < TCAT_nvox ; iv++ )
00869 rltsum[iv] += (rlt0[iv] + (0.5*ns)*rlt1[iv])*ns ;
00870 }
00871
00872 for( kk=ivbot ; kk < ivtop ; kk++ ){
00873 qq = kk - ivbot ;
00874 switch( DSET_BRICK_TYPE(new_dset,kk) ){
00875
00876 #undef ROUND
00877 #define ROUND(qq) ((short)rint((qq)))
00878
00879 case MRI_short:{
00880 short * bar = (short *) DSET_ARRAY(new_dset,kk) ;
00881 float fac = DSET_BRICK_FACTOR(new_dset,kk) , val,finv ;
00882
00883 if( fac == 0.0 ) fac = 1.0 ;
00884 finv = 1.0 / fac ;
00885 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
00886 val = fac*bar[iv] - rlt0[iv] - rlt1[iv]*(qq-qmid) ;
00887 bar[iv] = ROUND(finv*val) ;
00888 }
00889 }
00890 break ;
00891
00892 case MRI_float:{
00893 float * bar = (float *) DSET_ARRAY(new_dset,kk) ;
00894 float fac = DSET_BRICK_FACTOR(new_dset,kk) , val,finv ;
00895
00896 if( fac == 0.0 ) fac = 1.0 ;
00897 finv = 1.0 / fac ;
00898 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
00899 val = fac*bar[iv] - rlt0[iv] - rlt1[iv]*(qq-qmid) ;
00900 bar[iv] = (finv*val) ;
00901 #if 0
00902 fprintf(stderr,"kk=%d iv=%d bar=%g rlt0=%g rlt1=%g qq=%g qmid=%g val=%g\n",
00903 kk,iv,bar[iv],rlt0[iv],rlt1[iv],qq,qmid,val) ;
00904 #endif
00905 }
00906 }
00907 break ;
00908 }
00909 }
00910 }
00911 }
00912 }
00913
00914 }
00915
00916
00917
00918 if( TCAT_rlt == 3 && rltsum != NULL && nrltsum > 0 ){
00919 float scl = 1.0/nrltsum ; int kk ;
00920
00921 for( iv=0 ; iv < TCAT_nvox ; iv++ ) rltsum[iv] *= scl ;
00922
00923 for( kk=0 ; kk < new_nvals ; kk++ ){
00924 switch( DSET_BRICK_TYPE(new_dset,kk) ){
00925 case MRI_short:{
00926 short * bar = (short *) DSET_ARRAY(new_dset,kk) ;
00927 float fac = DSET_BRICK_FACTOR(new_dset,kk) , val,finv ;
00928
00929 if( fac == 0.0 ) fac = 1.0 ;
00930 finv = 1.0 / fac ;
00931 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
00932 val = fac*bar[iv] + rltsum[iv] ; bar[iv] = ROUND(finv*val) ;
00933 }
00934 }
00935 break ;
00936
00937 case MRI_float:{
00938 float * bar = (float *) DSET_ARRAY(new_dset,kk) ;
00939 float fac = DSET_BRICK_FACTOR(new_dset,kk) , val,finv ;
00940
00941 if( fac == 0.0 ) fac = 1.0 ;
00942 finv = 1.0 / fac ;
00943 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
00944 val = fac*bar[iv] + rltsum[iv] ; bar[iv] = (finv*val) ;
00945 }
00946 }
00947 break ;
00948 }
00949 }
00950 }
00951
00952 if( TCAT_rlt ){ free(rlt0); free(rlt1); if(rltsum!=NULL)free(rltsum); }
00953
00954 if( ! TCAT_dry ){
00955 if( TCAT_verb ) fprintf(stderr,"-verb: computing sub-brick statistics\n") ;
00956 THD_load_statistics( new_dset ) ;
00957
00958 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
00959 if( TCAT_verb ) fprintf(stderr,"-verb: Wrote output to %s\n",
00960 DSET_BRIKNAME(new_dset) ) ;
00961 }
00962
00963 exit(0) ;
00964 }
