00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #define ALLOW_BUCKETS
00036 #define ALLOW_SUBV
00037
00038 #include "mrilib.h"
00039 #include "parser.h"
00040
00041 #ifndef myXtFree
00042 #define myXtFree(xp) (XtFree((char *)(xp)) , (xp)=NULL)
00043 #endif
00044
00045
00046
00047 static int CALC_datum = ILLEGAL_TYPE ;
00048 static int CALC_nvox = -1 ;
00049 static PARSER_code * CALC_code = NULL ;
00050 static int ntime[26] ;
00051 static int ntime_max = 0 ;
00052 static int CALC_fscale = 0 ;
00053 static int CALC_gscale = 0 ;
00054 static int CALC_nscale = 0 ;
00055
00056 static int CALC_histpar = -1 ;
00057
00058
00059
00060 #define DSHIFT_MODE_STOP 0
00061 #define DSHIFT_MODE_WRAP 1
00062 #define DSHIFT_MODE_ZERO 2
00063
00064 static int CALC_dshift [26] ;
00065 static int CALC_dshift_i [26] ;
00066 static int CALC_dshift_j [26] ;
00067 static int CALC_dshift_k [26] ;
00068 static int CALC_dshift_l [26] ;
00069 static int CALC_dshift_mode[26] ;
00070
00071 static int CALC_dshift_mode_current = DSHIFT_MODE_STOP ;
00072 static int CALC_has_timeshift = 0 ;
00073
00074
00075
00076 static int CALC_has_sym[26] ;
00077 static char abet[] = "abcdefghijklmnopqrstuvwxyz" ;
00078
00079 #define HAS_I CALC_has_sym[ 8]
00080 #define HAS_J CALC_has_sym[ 9]
00081 #define HAS_K CALC_has_sym[10]
00082 #define HAS_X CALC_has_sym[23]
00083 #define HAS_Y CALC_has_sym[24]
00084 #define HAS_Z CALC_has_sym[25]
00085 #define HAS_T CALC_has_sym[19]
00086 #define HAS_L CALC_has_sym[11]
00087
00088 #define PREDEFINED_MASK ((1<< 8)|(1<< 9)|(1<<10)|(1<<11)| \
00089 (1<<19)|(1<<23)|(1<<24)|(1<<25) )
00090
00091 static int CALC_has_predefined = 0 ;
00092 static int CALC_has_xyz = 0 ;
00093 static int CALC_mangle_xyz = 0 ;
00094
00095 #define MANGLE_NONE 0
00096 #define MANGLE_RAI 1
00097 #define MANGLE_LPI 2
00098
00099 static THD_3dim_dataset * CALC_dset[26] ;
00100 static int CALC_type[26] ;
00101 static byte ** CALC_byte[26] ;
00102 static short ** CALC_short[26] ;
00103 static float ** CALC_float[26] ;
00104 static float * CALC_ffac[26] ;
00105 static int CALC_noffac[26] ;
00106
00107 static int CALC_verbose = 0 ;
00108
00109 static char CALC_output_prefix[THD_MAX_PREFIX] = "calc" ;
00110
00111 static char CALC_session[THD_MAX_NAME] = "./" ;
00112
00113 static MRI_IMAGE * TS_flim[26] ;
00114 static float * TS_flar[26] ;
00115 static int TS_nmax = 0 ;
00116 static int TS_make = 0 ;
00117 static float TS_dt = 1.0 ;
00118
00119 static MRI_IMAGE * IJKAR_flim[26] ;
00120 static float * IJKAR_flar[26] ;
00121 static int IJKAR_dcod[26] ;
00122
00123
00124
00125
00126 #define VAR_DEFINED(kv) \
00127 (TS_flim[kv] != NULL || IJKAR_flim[kv] != NULL || \
00128 CALC_dset[kv] != NULL || CALC_dshift[kv] >= 0 )
00129
00130 static float Rfac = 0.299 ;
00131 static float Gfac = 0.587 ;
00132 static float Bfac = 0.114 ;
00133
00134 static int CALC_taxis_num = 0 ;
00135
00136
00137 void CALC_read_opts( int , char ** ) ;
00138 void CALC_Syntax(void) ;
00139 int TS_reader( int , char * ) ;
00140 int IJKAR_reader( int , char * ) ;
00141
00142
00143
00144
00145
00146
00147 int TS_reader( int ival , char *fname )
00148 {
00149 MRI_IMAGE *tsim ;
00150
00151 if( ival < 0 || ival >= 26 ) return -1 ;
00152
00153 tsim = mri_read_1D( fname ) ;
00154 if( tsim == NULL ) return -1 ;
00155 if( tsim->nx < 2 ){ mri_free(tsim) ; return -1 ; }
00156
00157 TS_flim[ival] = tsim ;
00158 TS_nmax = MAX( TS_nmax , TS_flim[ival]->nx ) ;
00159 TS_flar[ival] = MRI_FLOAT_PTR( TS_flim[ival] ) ;
00160 return 0 ;
00161 }
00162
00163
00164
00165
00166
00167
00168 int IJKAR_reader( int ival , char *fname )
00169 {
00170 MRI_IMAGE *tsim ;
00171
00172 if( ival < 0 || ival >= 26 ) return -1 ;
00173
00174 tsim = mri_read_1D( fname ) ;
00175 if( tsim == NULL ) return -1 ;
00176 if( tsim->nx < 2 ){ mri_free(tsim) ; return -1 ; }
00177
00178 IJKAR_flim[ival] = tsim ;
00179 IJKAR_flar[ival] = MRI_FLOAT_PTR( IJKAR_flim[ival] ) ;
00180 return 0 ;
00181 }
00182
00183
00184
00185
00186
00187 void CALC_read_opts( int argc , char * argv[] )
00188 {
00189 int nopt = 1 ;
00190 int ids ;
00191 int ii, kk;
00192
00193 for( ids=0 ; ids < 26 ; ids++ ){
00194 CALC_dset[ids] = NULL ;
00195 CALC_type[ids] = -1 ;
00196 TS_flim[ids] = NULL ;
00197 IJKAR_flim[ids] = NULL ;
00198
00199 CALC_dshift[ids] = -1 ;
00200 CALC_dshift_mode[ids] = CALC_dshift_mode_current ;
00201
00202 CALC_noffac[ids] = 1 ;
00203 }
00204
00205 while( nopt < argc && argv[nopt][0] == '-' ){
00206
00207
00208
00209 if( strcasecmp(argv[nopt],"-dicom") == 0 || strcasecmp(argv[nopt],"-rai") == 0 ){
00210 CALC_mangle_xyz = MANGLE_RAI ;
00211 nopt++ ; continue ;
00212 }
00213
00214 if( strcasecmp(argv[nopt],"-spm") == 0 || strcasecmp(argv[nopt],"-lpi") == 0 ){
00215 CALC_mangle_xyz = MANGLE_LPI ;
00216 nopt++ ; continue ;
00217 }
00218
00219
00220
00221 if( strncasecmp(argv[nopt],"-rgbfac",7) == 0 ){
00222 if( ++nopt >= argc )
00223 ERROR_exit("need an argument after -rgbfac!\n") ;
00224 Rfac = strtod( argv[nopt++] , NULL ) ;
00225 Gfac = strtod( argv[nopt++] , NULL ) ;
00226 Bfac = strtod( argv[nopt++] , NULL ) ;
00227 if( Rfac == 0.0 && Gfac == 0.0 && Bfac == 0.0 )
00228 ERROR_exit("All 3 factors after -rgbfac are zero!?\n");
00229 continue ;
00230 }
00231
00232
00233
00234 if( strncasecmp(argv[nopt],"-taxis",6) == 0 ){
00235 char *cpt ;
00236 if( ++nopt >= argc )
00237 ERROR_exit("need an argument after -taxis!\n") ;
00238 CALC_taxis_num = strtod( argv[nopt] , &cpt ) ;
00239 if( CALC_taxis_num < 2 )
00240 ERROR_exit("N value after -taxis must be bigger than 1!\n");
00241 if( *cpt == ':' ){
00242 float dt = strtod( cpt+1 , &cpt ) ;
00243 if( dt > 0.0 ){
00244 TS_dt = dt ;
00245 if( *cpt == 'm' && *(cpt+1) == 's' ) TS_dt *= 0.001 ;
00246 } else {
00247 WARNING_message("time step value in '-taxis %s' not legal!\n",argv[nopt]);
00248 }
00249 }
00250 nopt++ ; continue ;
00251 }
00252
00253
00254
00255 if( strncasecmp(argv[nopt],"-dt",3) == 0 || strncmp(argv[nopt],"-TR",3) == 0 ){
00256 char *cpt ;
00257 if( ++nopt >= argc )
00258 ERROR_exit("need an argument after -dt!\n") ;
00259 TS_dt = strtod( argv[nopt] , &cpt ) ;
00260 if( TS_dt <= 0.0 )
00261 ERROR_exit("Illegal time step value after -dt!\n");
00262 if( *cpt == 'm' && *(cpt+1) == 's' ) TS_dt *= 0.001 ;
00263 nopt++ ; continue ;
00264 }
00265
00266
00267
00268 if( strncasecmp(argv[nopt],"-histpar",5) == 0 ){
00269 if( ++nopt >= argc )
00270 ERROR_exit("need an argument after -histpar!\n") ;
00271 if( argv[nopt][0] < 'a' || argv[nopt][0] > 'z')
00272 ERROR_exit("argument after -histpar is illegal!\n");
00273 CALC_histpar = (int) (argv[nopt][0] - 'a') ;
00274
00275 nopt++ ; continue ;
00276 }
00277
00278
00279
00280 if( strncasecmp(argv[nopt],"-datum",6) == 0 ){
00281 if( ++nopt >= argc )
00282 ERROR_exit("need an argument after -datum!\n") ;
00283 if( strcasecmp(argv[nopt],"short") == 0 ){
00284 CALC_datum = MRI_short ;
00285 } else if( strcasecmp(argv[nopt],"float") == 0 ){
00286 CALC_datum = MRI_float ;
00287 } else if( strcasecmp(argv[nopt],"byte") == 0 ){
00288 CALC_datum = MRI_byte ;
00289 } else if( strcasecmp(argv[nopt],"complex") == 0 ){
00290 CALC_datum = MRI_complex ;
00291 } else {
00292 ERROR_exit("-datum of type '%s' not supported in 3dcalc!\n",argv[nopt]) ;
00293 }
00294 nopt++ ; continue ;
00295 }
00296
00297
00298
00299 if( strncasecmp(argv[nopt],"-verbose",5) == 0 ){
00300 CALC_verbose = 1 ;
00301 nopt++ ; continue ;
00302 }
00303
00304
00305
00306 if( strncasecmp(argv[nopt],"-nscale",6) == 0 ){
00307 CALC_gscale = CALC_fscale = 0 ;
00308 CALC_nscale = 1 ;
00309 nopt++ ; continue ;
00310 }
00311
00312
00313
00314 if( strncasecmp(argv[nopt],"-fscale",6) == 0 ){
00315 CALC_fscale = 1 ;
00316 CALC_nscale = 0 ;
00317 nopt++ ; continue ;
00318 }
00319
00320
00321
00322 if( strncasecmp(argv[nopt],"-gscale",6) == 0 ){
00323 CALC_gscale = CALC_fscale = 1 ;
00324 CALC_nscale = 0 ;
00325 nopt++ ; continue ;
00326 }
00327
00328
00329
00330 if( strncasecmp(argv[nopt],"-prefix",6) == 0 ){
00331 nopt++ ;
00332 if( nopt >= argc )
00333 ERROR_exit("need argument after -prefix!\n") ;
00334 MCW_strncpy( CALC_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00335 continue ;
00336 }
00337
00338
00339
00340 if( strncasecmp(argv[nopt],"-session",6) == 0 ){
00341 nopt++ ;
00342 if( nopt >= argc )
00343 ERROR_exit("need argument after -session!\n") ;
00344 MCW_strncpy( CALC_session , argv[nopt++] , THD_MAX_NAME ) ;
00345 continue ;
00346 }
00347
00348
00349
00350 if( strncasecmp(argv[nopt],"-expr",4) == 0 ){
00351 if( CALC_code != NULL )
00352 ERROR_exit("cannot have 2 -expr options!\n") ;
00353 nopt++ ;
00354 if( nopt >= argc )
00355 ERROR_exit("need argument after -expr!\n") ;
00356 PARSER_set_printout(1) ;
00357 CALC_code = PARSER_generate_code( argv[nopt++] ) ;
00358 if( CALC_code == NULL )
00359 ERROR_exit("illegal expression!\n") ;
00360 PARSER_mark_symbols( CALC_code , CALC_has_sym ) ;
00361 continue ;
00362 }
00363
00364
00365
00366 if( strncasecmp(argv[nopt],"-dsSTOP",6) == 0 ){
00367 CALC_dshift_mode_current = DSHIFT_MODE_STOP ;
00368 nopt++ ; continue ;
00369 }
00370
00371
00372
00373 if( strncasecmp(argv[nopt],"-dsWRAP",6) == 0 ){
00374 CALC_dshift_mode_current = DSHIFT_MODE_WRAP ;
00375 nopt++ ; continue ;
00376 }
00377
00378
00379
00380 if( strncasecmp(argv[nopt],"-dsZERO",6) == 0 ){
00381 CALC_dshift_mode_current = DSHIFT_MODE_ZERO ;
00382 nopt++ ; continue ;
00383 }
00384
00385
00386
00387 ids = strlen( argv[nopt] ) ;
00388
00389 if( (argv[nopt][1] >= 'a' && argv[nopt][1] <= 'z') &&
00390 (ids == 2 ||
00391 (ids > 2 && argv[nopt][2] >= '0' && argv[nopt][2] <= '9')) ){
00392
00393 int ival , nxyz , isub , ll ;
00394 THD_3dim_dataset * dset ;
00395
00396 ival = argv[nopt][1] - 'a' ;
00397 if( VAR_DEFINED(ival) )
00398 ERROR_exit("Can't define %c symbol twice\n",argv[nopt][1]);
00399
00400 isub = (ids == 2) ? 0 : strtol(argv[nopt]+2,NULL,10) ;
00401 if( isub < 0 )
00402 ERROR_exit("Illegal sub-brick value: %s\n",argv[nopt]) ;
00403
00404 nopt++ ;
00405 if( nopt >= argc )
00406 ERROR_exit("need argument after %s\n",argv[nopt-1]);
00407
00408
00409
00410 ll = strlen(argv[nopt]) ;
00411 if( ll >= 4 &&
00412 strstr(argv[nopt],"1D") != NULL &&
00413 argv[nopt][1] == ':' &&
00414 (argv[nopt][0] == 'I' || argv[nopt][0] == 'i' ||
00415 argv[nopt][0] == 'J' || argv[nopt][0] == 'j' ||
00416 argv[nopt][0] == 'K' || argv[nopt][0] == 'k' ) ){
00417
00418 ll = IJKAR_reader( ival , argv[nopt]+2 ) ;
00419 if( ll == 0 ){
00420 switch( argv[nopt][0] ){
00421 case 'I': case 'i': IJKAR_dcod[ival] = 8 ; break ;
00422 case 'J': case 'j': IJKAR_dcod[ival] = 9 ; break ;
00423 case 'K': case 'k': IJKAR_dcod[ival] = 10 ; break ;
00424 }
00425 nopt++ ; goto DSET_DONE ;
00426 }
00427 }
00428
00429
00430
00431 ll = strlen(argv[nopt]) ;
00432 if( ll >= 4 && ( strstr(argv[nopt],".1D") != NULL ||
00433 strstr(argv[nopt],"1D:") != NULL ) ){
00434
00435 ll = TS_reader( ival , argv[nopt] ) ;
00436 if( ll == 0 ){ nopt++ ; goto DSET_DONE ; }
00437
00438
00439 }
00440
00441
00442
00443
00444 if( (argv[nopt][0] >= 'a' && argv[nopt][0] <= 'z') &&
00445 ( (ll >= 3 && argv[nopt][1] == '[') ||
00446 (ll == 3 &&
00447 (argv[nopt][1] == '+' || argv[nopt][1] == '-'))
00448 ) ){
00449
00450 int jds = argv[nopt][0] - 'a' ;
00451 int * ijkl ;
00452
00453
00454
00455 if( ids > 2 )
00456 ERROR_exit("Can't combine %s with differential subscripting %s\n",
00457 argv[nopt-1],argv[nopt]) ;
00458 if( CALC_dset[jds] == NULL )
00459 ERROR_exit("Must define dataset %c before using it in %s\n",
00460 argv[nopt][0] , argv[nopt] ) ;
00461
00462
00463
00464 if( argv[nopt][1] == '[' ){
00465 MCW_intlist_allow_negative(1) ;
00466 ijkl = MCW_get_intlist( 9999 , argv[nopt]+1 ) ;
00467 MCW_intlist_allow_negative(0) ;
00468 if( ijkl == NULL || ijkl[0] != 4 )
00469 ERROR_exit("Illegal differential subscripting %s\n",
00470 argv[nopt] ) ;
00471 } else {
00472 ijkl = (int *) malloc( sizeof(int) * 5 ) ;
00473 ijkl[1] = ijkl[2] = ijkl[3] = ijkl[4] = 0 ;
00474 switch( argv[nopt][2] ){
00475 default:
00476 ERROR_exit("Bad differential subscripting %s\n",argv[nopt]);
00477
00478 case 'i': ijkl[1] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
00479 case 'j': ijkl[2] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
00480 case 'k': ijkl[3] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
00481 case 'l': ijkl[4] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
00482 }
00483 }
00484
00485
00486
00487 if( ijkl[1]==0 && ijkl[2]==0 && ijkl[3]==0 && ijkl[4]==0 )
00488 WARNING_message("differential subscript %s is all zero\n",argv[nopt]);
00489
00490 if( ntime[jds] == 1 && ijkl[4] != 0 ){
00491 WARNING_message(
00492 "differential subscript %s has nonzero time\n"
00493 " + shift on base dataset with 1 sub-brick!\n"
00494 " + Setting time shift to 0.\n" ,
00495 argv[nopt] ) ;
00496 ijkl[4] = 0 ;
00497 }
00498
00499
00500
00501 CALC_dshift [ival] = jds ;
00502 CALC_dshift_i[ival] = ijkl[1] ;
00503 CALC_dshift_j[ival] = ijkl[2] ;
00504 CALC_dshift_k[ival] = ijkl[3] ;
00505 CALC_dshift_l[ival] = ijkl[4] ;
00506
00507 CALC_dshift_mode[ival] = CALC_dshift_mode_current ;
00508
00509 CALC_has_timeshift = CALC_has_timeshift || (ijkl[4] != 0) ;
00510
00511
00512
00513 free(ijkl) ; nopt++ ; goto DSET_DONE ;
00514
00515 }
00516
00517
00518
00519 #ifndef ALLOW_SUBV
00520 dset = THD_open_one_dataset( argv[nopt++] ) ;
00521 if( dset == NULL )
00522 ERROR_exit("can't open dataset %s\n",argv[nopt-1]) ;
00523 if( isub >= DSET_NVALS(dset) )
00524 ERROR_exit("dataset %s only has %d sub-bricks\n",
00525 argv[nopt-1],DSET_NVALS(dset)) ;
00526 #else
00527 { char dname[512] ;
00528
00529 if( ids > 2 ){
00530 if( strstr(argv[nopt],"[") != NULL ){
00531 ERROR_exit(
00532 "Illegal combination of sub-brick specifiers: "
00533 "%s %s\n" ,
00534 argv[nopt-1] , argv[nopt] ) ;
00535 }
00536 sprintf(dname,"%s[%d]",argv[nopt++],isub) ;
00537 isub = 0 ;
00538 } else {
00539 strcpy(dname,argv[nopt++]) ;
00540 }
00541 dset = THD_open_dataset( dname ) ;
00542 if( dset == NULL )
00543 ERROR_exit("can't open dataset %s\n",dname) ;
00544 }
00545 #endif
00546
00547
00548
00549 #ifdef ALLOW_BUCKETS
00550 ntime[ival] = DSET_NVALS(dset) ;
00551 #else
00552 ntime[ival] = DSET_NUM_TIMES(dset);
00553 #endif
00554 if ( ids > 2 ) ntime[ival] = 1 ;
00555 ntime_max = MAX( ntime_max, ntime[ival] );
00556
00557 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz ;
00558 if( CALC_nvox < 0 ){
00559 CALC_nvox = nxyz ;
00560 } else if( nxyz != CALC_nvox ){
00561 ERROR_exit("dataset %s differs in size from others\n",argv[nopt-1]);
00562 }
00563
00564 if( !DSET_datum_constant(dset) ){
00565 float *far ;
00566
00567 WARNING_message("dataset %s has sub-bricks with different types\n"
00568 " + ==> converting all sub-bricks to floats\n",
00569 argv[nopt-1]);
00570
00571 DSET_mallocize(dset) ; DSET_load(dset) ;
00572 if( ! DSET_LOADED(dset) )
00573 ERROR_exit("can't load %s from disk!\n",argv[nopt-1]);
00574
00575 for( ii=0 ; ii < ntime[ival] ; ii++ ){
00576 if( DSET_BRICK_TYPE(dset,ii) != MRI_float ){
00577 far = calloc( sizeof(float) , nxyz ) ;
00578 if( far == NULL )
00579 ERROR_exit("can't malloc space for conversion\n");
00580 EDIT_coerce_scale_type( nxyz , DSET_BRICK_FACTOR(dset,ii) ,
00581 DSET_BRICK_TYPE(dset,ii), DSET_ARRAY(dset,ii),
00582 MRI_float , far ) ;
00583 EDIT_substitute_brick( dset , ii , MRI_float , far ) ;
00584 DSET_BRICK_FACTOR(dset,ii) = 0.0 ;
00585 }
00586 }
00587 }
00588
00589 CALC_type[ival] = DSET_BRICK_TYPE(dset,isub) ;
00590 CALC_dset[ival] = dset ;
00591
00592
00593
00594
00595
00596 CALC_ffac[ival] = (float *) malloc( sizeof(float) * ntime[ival] ) ;
00597 if ( ntime[ival] == 1 ) {
00598 CALC_ffac[ival][0] = DSET_BRICK_FACTOR( dset , isub) ;
00599 if (CALC_ffac[ival][0] == 0.0 ) CALC_ffac[ival][0] = 1.0 ;
00600 if( CALC_ffac[ival][0] != 1.0 ) CALC_noffac[ival] = 0 ;
00601 } else {
00602 for (ii = 0 ; ii < ntime[ival] ; ii ++ ) {
00603 CALC_ffac[ival][ii] = DSET_BRICK_FACTOR(dset, ii) ;
00604 if (CALC_ffac[ival][ii] == 0.0 ) CALC_ffac[ival][ii] = 1.0;
00605 if( CALC_ffac[ival][ii] != 1.0 ) CALC_noffac[ival] = 0 ;
00606 }
00607 }
00608
00609
00610
00611 if( CALC_verbose ){
00612 int iv , nb ;
00613 for( iv=nb=0 ; iv < DSET_NVALS(dset) ; iv++ )
00614 nb += DSET_BRICK_BYTES(dset,iv) ;
00615 INFO_message("Reading dataset %s (%d bytes)\n",argv[nopt-1],nb);
00616 }
00617
00618 if( ! DSET_LOADED(dset) ){
00619 THD_load_datablock( dset->dblk ) ;
00620 if( ! DSET_LOADED(dset) )
00621 ERROR_exit("Can't read data brick for dataset %s\n",argv[nopt-1]) ;
00622 }
00623
00624
00625 switch (CALC_type[ival]) {
00626 case MRI_short:
00627 CALC_short[ival] = (short **) malloc( sizeof(short *) * ntime[ival] ) ;
00628 if (ntime[ival] == 1 )
00629 CALC_short[ival][0] = (short *) DSET_ARRAY(dset, isub) ;
00630 else
00631 for (ii=0; ii < ntime[ival]; ii++)
00632 CALC_short[ival][ii] = (short *) DSET_ARRAY(dset, ii);
00633 break;
00634
00635 case MRI_float:
00636 CALC_float[ival] = (float **) malloc( sizeof(float *) * ntime[ival] ) ;
00637 if (ntime[ival] == 1 )
00638 CALC_float[ival][0] = (float *) DSET_ARRAY(dset, isub) ;
00639 else
00640 for (ii=0; ii < ntime[ival]; ii++)
00641 CALC_float[ival][ii] = (float *) DSET_ARRAY(dset, ii);
00642 break;
00643
00644 case MRI_byte:
00645 CALC_byte[ival] = (byte **) malloc( sizeof(byte *) * ntime[ival] ) ;
00646 if (ntime[ival] == 1 )
00647 CALC_byte[ival][0] = (byte *) DSET_ARRAY(dset, isub) ;
00648 else
00649 for (ii=0; ii < ntime[ival]; ii++)
00650 CALC_byte[ival][ii] = (byte *) DSET_ARRAY(dset, ii);
00651 break;
00652
00653 case MRI_rgb:
00654 CALC_byte[ival] = (byte **) malloc( sizeof(byte *) * ntime[ival] ) ;
00655 if (ntime[ival] == 1 )
00656 CALC_byte[ival][0] = (byte *) DSET_ARRAY(dset, isub) ;
00657 else
00658 for (ii=0; ii < ntime[ival]; ii++)
00659 CALC_byte[ival][ii] = (byte *) DSET_ARRAY(dset, ii);
00660 break ;
00661
00662 default:
00663 ERROR_exit("Dataset %s has illegal data type: %s\n" ,
00664 argv[nopt-1] , MRI_type_name[CALC_type[ival]] ) ;
00665
00666 }
00667 if( CALC_datum < 0 && CALC_type[ival] != MRI_rgb ) CALC_datum = CALC_type[ival] ;
00668
00669 DSET_DONE: continue;
00670
00671 }
00672
00673 ERROR_exit("Unknown option: %s\n",argv[nopt]) ;
00674
00675 }
00676
00677
00678
00679
00680 if( nopt < argc )
00681 ERROR_exit("Extra command line arguments puzzle me! argv[%d]=%s ...\n",nopt,argv[nopt]) ;
00682
00683 for( ids=0 ; ids < 26 ; ids++ ) if( CALC_dset[ids] != NULL ) break ;
00684 if( ids == 26 )
00685 ERROR_exit("No actual input datasets given!\n") ;
00686
00687
00688
00689 for( ii=0 ; ii < 26 ; ii++ ){
00690 if( IJKAR_flim[ii] != NULL ){
00691 int siz ;
00692 switch( IJKAR_dcod[ii] ){
00693 case 8: siz = DSET_NX(CALC_dset[ids]) ; break ;
00694 case 9: siz = DSET_NY(CALC_dset[ids]) ; break ;
00695 case 10: siz = DSET_NZ(CALC_dset[ids]) ; break ;
00696 }
00697 if( IJKAR_flim[ii]->nx != siz )
00698 ERROR_message("dimension mismatch between '-%c' and '%-c'\n", 'a'+ii , 'a'+ids ) ;
00699 }
00700 }
00701
00702 if( CALC_code == NULL ) ERROR_exit("No expression given!\n") ;
00703
00704 if( CALC_histpar >= 0 && CALC_dset[CALC_histpar] == NULL ){
00705 WARNING_message("-histpar dataset not defined!\n") ;
00706 CALC_histpar = -1 ;
00707 }
00708
00709 for (ids=0; ids < 26; ids ++)
00710 if (ntime[ids] > 1 && ntime[ids] != ntime_max ) {
00711 #ifdef ALLOW_BUCKETS
00712 ERROR_exit("Multi-brick datasets don't match!\n") ;
00713 #else
00714 ERROR_exit("3D+time datasets don't match!\n") ;
00715 #endif
00716 }
00717
00718
00719
00720
00721
00722 if( ntime_max == 1 && TS_nmax > 0 ){
00723 ntime_max = TS_nmax ;
00724 TS_make = 1 ;
00725 INFO_message(
00726 "Calculating 3D+time[%d]"
00727 " dataset from 3D datasets and time series, with dt=%g s\n" ,
00728 ntime_max , TS_dt ) ;
00729 }
00730
00731 if( CALC_taxis_num > 0 ){
00732 if( ntime_max > 1 ){
00733 WARNING_message("-taxis %d overriden by dataset input(s)\n",
00734 CALC_taxis_num) ;
00735 } else {
00736 ntime_max = CALC_taxis_num ;
00737 TS_make = 1 ;
00738 INFO_message("Calculating 3D+time[%d]"
00739 " dataset from 3D datasets and -taxis with dt=%g s\n" ,
00740 ntime_max , TS_dt ) ;
00741 }
00742 }
00743
00744
00745
00746
00747 for (ids=0; ids < 26; ids ++){
00748 if( VAR_DEFINED(ids) && !CALC_has_sym[ids] )
00749 WARNING_message("input '%c' is not used in the expression\n" ,
00750 abet[ids] ) ;
00751
00752 else if( !VAR_DEFINED(ids) && CALC_has_sym[ids] ){
00753
00754 if( ((1<<ids) & PREDEFINED_MASK) == 0 ){
00755 WARNING_message( "symbol %c is used but not defined\n" , abet[ids] ) ;
00756 } else {
00757 CALC_has_predefined++ ;
00758 INFO_message("Symbol %c using predefined value\n",abet[ids] ) ;
00759 if( ids >= 23 ) CALC_has_xyz = 1 ;
00760 }
00761 }
00762 }
00763
00764 return ;
00765 }
00766
00767
00768
00769 void CALC_Syntax(void)
00770 {
00771 printf(
00772 "Program: 3dcalc \n"
00773 "Author: RW Cox et al \n"
00774 " \n"
00775 "3dcalc - AFNI's calculator program \n"
00776 " \n"
00777 " This program does voxel-by-voxel arithmetic on 3D datasets \n"
00778 " (limited to inter-voxel computation). \n"
00779 " \n"
00780 " The program assumes that the voxel-by-voxel computations are being \n"
00781 " performed on datasets that occupy the same space and have the same \n"
00782 " orientations. \n"
00783 " \n"
00784 "------------------------------------------------------------------------\n"
00785 "Usage: \n"
00786 "----- \n"
00787 " 3dcalc -a dsetA [-b dsetB...] \\ \n"
00788 " -expr EXPRESSION \\ \n"
00789 " [options] \n"
00790 " \n"
00791 "Examples: \n"
00792 "-------- \n"
00793 "1. Average datasets together, on a voxel-by-voxel basis: \n"
00794 " \n"
00795 " 3dcalc -a fred+tlrc -b ethel+tlrc -c lucy+tlrc \\ \n"
00796 " -expr '(a+b+c)/3' -prefix subjects_mean \n"
00797 " \n"
00798 "2. Perform arithmetic calculations between the sub-bricks of a single \n"
00799 " dataset by noting the sub-brick number on the command line: \n"
00800 " \n"
00801 " 3dcalc -a 'func+orig[2]' -b 'func+orig[4]' -expr 'sqrt(a*b)' \n"
00802 " \n"
00803 "3. Create a simple mask that consists only of values in sub-brick #0 \n"
00804 " that are greater than 3.14159: \n"
00805 " \n"
00806 " 3dcalc -a 'func+orig[0]' -expr 'ispositive(a-3.14159)' \\ \n"
00807 " -prefix mask \n"
00808 " \n"
00809 "4. Normalize subjects' time series datasets to percent change values in \n"
00810 " preparation for group analysis: \n"
00811 " \n"
00812 " Voxel-by-voxel, the example below divides each intensity value in \n"
00813 " the time series (epi_r1+orig) with the voxel's mean value (mean+orig)\n"
00814 " to get a percent change value. The 'ispositive' command will ignore \n"
00815 " voxels with mean values less than 167 (i.e., they are labeled as \n"
00816 " 'zero' in the output file 'percent_change+orig') and are most likely \n"
00817 " background/noncortical voxels. \n"
00818 " \n"
00819 " 3dcalc -a epi_run1+orig -b mean+orig \\ \n"
00820 " -expr '100 * a/b * ispositive(b-167)' -prefix percent_chng \n"
00821 " \n"
00822 "5. Create a compound mask from a statistical dataset, where 3 stimuli \n"
00823 " show activation. \n"
00824 " NOTE: 'step' and 'ispositive' are identical expressions that can \n"
00825 " be used interchangeably: \n"
00826 " \n"
00827 " 3dcalc -a 'func+orig[12]' -b 'func+orig[15]' -c 'func+orig[18]' \\ \n"
00828 " -expr 'step(a-4.2)*step(b-2.9)*step(c-3.1)' \\ \n"
00829 " -prefix compound_mask \n"
00830 " \n"
00831 "6. Same as example #5, but this time create a mask of 8 different values\n"
00832 " showing all combinations of activations (i.e., not only where \n"
00833 " everything is active, but also each stimulus individually, and all \n"
00834 " combinations). The output mask dataset labels voxel values as such: \n"
00835 " 0 = none active 1 = A only active 2 = B only active \n"
00836 " 3 = A and B only 4 = C only active 5 = A and C only \n"
00837 " 6 = B and C only 7 = all A, B, and C active \n"
00838 " \n"
00839 " 3dcalc -a 'func+orig[12]' -b 'func+orig[15]' -c 'func+orig[18]' \\ \n"
00840 " -expr 'step(a-4.2)+2*step(b-2.9)+4*step(c-3.1)' \\ \n"
00841 " -prefix mask_8 \n"
00842 " \n"
00843 "7. Create a region-of-interest mask comprised of a 3-dimensional sphere.\n"
00844 " Values within the ROI sphere will be labeled as '1' while values \n"
00845 " outside the mask will be labeled as '0'. Statistical analyses can \n"
00846 " then be done on the voxels within the ROI sphere. \n"
00847 " \n"
00848 " The example below puts a solid ball (sphere) of radius 3=sqrt(9) \n"
00849 " about the point with coordinates (x,y,z)=(20,30,70): \n"
00850 " \n"
00851 " 3dcalc -a anat+tlrc \\\n"
00852 " -expr 'step((9-(x-20)*(x-20)-(y-30)*(y-30)-(z-70)*(z-70))'\\\n"
00853 " -prefix ball \n"
00854 " \n"
00855 " 8. Some datsets are 'short' (16 bit) integers with a scalar attached, \n"
00856 " which allow them to be smaller than float datasets and to contain \n"
00857 " fractional values. \n"
00858 " \n"
00859 " Dataset 'a' is always used as a template for the output dataset. For\n"
00860 " the examples below, assume that datasets d1+orig and d2+orig consist\n"
00861 " of small integers. \n"
00862 " \n"
00863 " a) When dividing 'a' by 'b', the result should be scaled, so that a \n"
00864 " value of 2.4 is not truncated to '2'. To avoid this truncation, \n"
00865 " force scaling with the -fscale option: \n"
00866 " \n"
00867 " 3dcalc -a d1+orig -b d2+orig -expr 'a/b' -prefix quot -fscale \n"
00868 " \n"
00869 " b) If it is preferable that the result is of type 'float', then set \n"
00870 " the output data type (datum) to float: \n"
00871 " \n"
00872 " 3dcalc -a d1+orig -b d2+orig -expr 'a/b' -prefix quot \\ \n"
00873 " -datum float \n"
00874 " \n"
00875 " c) Perhaps an integral division is desired, so that 9/4=2, not 2.24.\n"
00876 " Force the results not to be scaled (opposite of example 8b) using\n"
00877 " the -nscale option: \n"
00878 " \n"
00879 " 3dcalc -a d1+orig -b d2+orig -expr 'a/b' -prefix quot -nscale \n"
00880 " \n"
00881 "------------------------------------------------------------------------\n"
00882 " \n"
00883 "ARGUMENTS for 3dcalc (must be included on command line): \n"
00884 "-------------------- ---- \n"
00885 " \n"
00886 " -a dname = Read dataset 'dname' and call the voxel values 'a' in the\n"
00887 " expression (-expr) that is input below. Up to 24 dnames \n"
00888 " (-a, -b, -c, ... -z) can be included in a single 3dcalc \n"
00889 " calculation/expression. \n"
00890 " ** If some letter name is used in the expression, but \n"
00891 " not present in one of the dataset options here, then \n"
00892 " that variable is set to 0. \n"
00893 " ** If the letter is followed by a number, then that \n"
00894 " number is used to select the sub-brick of the dataset \n"
00895 " which will be used in the calculations. \n"
00896 " E.g., '-b3 dname' specifies that the variable 'b' \n"
00897 " refers to sub-brick '3' of that dataset \n"
00898 " (indexes in AFNI start at 0). \n"
00899 " \n"
00900 " -expr = Apply the expression - within quotes - to the input \n"
00901 " datasets (dnames), one voxel at time, to produce the \n"
00902 " output dataset. \n"
00903 "------------------------------------------------------------------------\n"
00904 ) ;
00905 printf(
00906 " OPTIONS for 3dcalc: \n"
00907 " ------- \n"
00908 " \n"
00909 " -verbose = Makes the program print out various information as it \n"
00910 " progresses. \n"
00911 " \n"
00912 " -datum type= Coerce the output data to be stored as the given type, \n"
00913 " which may be byte, short, or float. \n"
00914 " [default = datum of first input dataset] \n"
00915 " \n"
00916 " -fscale = Force scaling of the output to the maximum integer \n"
00917 " range. This only has effect if the output datum is byte \n"
00918 " or short (either forced or defaulted). This option is \n"
00919 " often necessary to eliminate unpleasant truncation \n"
00920 " artifacts. \n"
00921 " [The default is to scale only if the computed values \n"
00922 " seem to need it -- are all <= 1.0 or there is at \n"
00923 " least one value beyond the integer upper limit.] \n"
00924 " \n"
00925 " ** In earlier versions of 3dcalc, scaling (if used) was \n"
00926 " applied to all sub-bricks equally -- a common scale \n"
00927 " factor was used. This would cause trouble if the \n"
00928 " values in different sub-bricks were in vastly \n"
00929 " different scales. In this version, each sub-brick \n"
00930 " gets its own scale factor. To override this behavior,\n"
00931 " use the '-gscale' option. \n"
00932 " \n"
00933 " -gscale = Same as '-fscale', but also forces each output sub-brick \n"
00934 " to get the same scaling factor. This may be desirable \n"
00935 " for 3D+time datasets, for example. \n"
00936 " \n"
00937 " -nscale = Don't do any scaling on output to byte or short datasets.\n"
00938 " This may be especially useful when operating on mask \n"
00939 " datasets whose output values are only 0's and 1's. \n"
00940 #ifndef ALLOW_SUBV
00941 " ** The type and number of sub-bricks in a dataset can be \n"
00942 " printed out using the '3dinfo' program. \n"
00943 #else
00944 " ** Another way to achieve the effect of '-b3' is described\n"
00945 " below in the dataset 'INPUT' specification section. \n"
00946 #endif
00947 " \n"
00948 " -prefix pname = Use 'pname' for the output dataset prefix name. \n"
00949 " [default='calc'] \n"
00950 " \n"
00951 " -session dir = Use 'dir' for the output dataset session directory. \n"
00952 " [default='./'=current working directory] \n"
00953 " \n"
00954 " -dt tstep = Use 'tstep' as the TR for manufactured 3D+time datasets.\n"
00955 " \n"
00956 " -TR tstep = If not given, defaults to 1 second. \n"
00957 " \n"
00958 " -taxis N = If only 3D datasets are input (no 3D+time or .1D files),\n"
00959 " *OR* then normally only a 3D dataset is calculated. With \n"
00960 " -taxis N:tstep: this option, you can force the creation of a time axis\n"
00961 " of length 'N', optionally using time step 'tstep'. In\n"
00962 " such a case, you will probably want to use the pre- \n"
00963 " defined time variables 't' and/or 'k' in your \n"
00964 " expression, or each resulting sub-brick will be \n"
00965 " identical. For example: \n"
00966 " '-taxis 121:0.1' will produce 121 points in time, \n"
00967 " spaced with TR 0.1. \n"
00968 " \n"
00969 " N.B.: You can also specify the TR using the -dt option. \n"
00970 " N.B.: You can specify 1D input datasets using the \n"
00971 " '1D:n@val,n@val' notation to get a similar effect. \n"
00972 " For example: \n"
00973 " -dt 0.1 -w '1D:121@0' \n"
00974 " will have pretty much the same effect as \n"
00975 " -taxis 121:0.1\n"
00976 " N.B.: For both '-dt' and '-taxis', the 'tstep' value is in \n"
00977 " seconds. You can suffix it with 'ms' to specify that\n"
00978 " the value is in milliseconds instead; e.g., '-dt 2000ms'.\n"
00979 " \n"
00980 " -rgbfac A B C = For RGB input datasets, the 3 channels (r,g,b) are \n"
00981 " collapsed to one for the purposes of 3dcalc, using the\n"
00982 " formula value = A*r + B*g + C*b \n"
00983 " \n"
00984 " The default values are A=0.299 B=0.587 C=0.114, which \n"
00985 " gives the grayscale intensity. To pick out the Green \n"
00986 " channel only, use '-rgbfac 0 1 0', for example. Note \n"
00987 " that each channel in an RGB dataset is a byte in the \n"
00988 " range 0..255. Thus, '-rgbfac 0.001173 0.002302 0.000447'\n"
00989 " will compute the intensity rescaled to the range 0..1.0\n"
00990 " (i.e., 0.001173=0.299/255, etc.) \n"
00991 " \n"
00992 "------------------------------------------------------------------------\n"
00993 "DATASET TYPES: \n"
00994 "------------- \n"
00995 " \n"
00996 " The most common AFNI dataset types are 'byte', 'short', and 'float'. \n"
00997 " \n"
00998 " A byte value is an 8-bit signed integer (0..255), a short value ia a \n"
00999 " 16-bit signed integer (-32768..32767), and a float value is a 32-bit \n"
01000 " real number. A byte value has almost 3 decimals of accuracy, a short \n"
01001 " has almost 5, and a float has approximately 7 (from a 23+1 bit \n"
01002 " mantissa). \n"
01003 " \n"
01004 " Datasets can also have a scalar attached to each sub-brick. The main \n"
01005 " use of this is allowing a short type dataset to take on non-integral \n"
01006 " values, while being half the size of a float dataset. \n"
01007 " \n"
01008 " As an example, consider a short dataset with a scalar of 0.0001. This \n"
01009 " could represent values between -32.768 and +32.767, at a resolution of \n"
01010 " 0.001. One could represnt the difference between 4.916 and 4.917, for \n"
01011 " instance, but not 4.9165. Each number has 15 bits of accuracy, plus a \n"
01012 " sign bit, which gives 4-5 decimal places of accuracy. If this is not \n"
01013 " enough, then it makes sense to use the larger type, float. \n"
01014 " \n"
01015 "------------------------------------------------------------------------\n"
01016 "3D+TIME DATASETS: \n"
01017 "---------------- \n"
01018 " \n"
01019 " This version of 3dcalc can operate on 3D+time datasets. Each input \n"
01020 " dataset will be in one of these conditions: \n"
01021 " \n"
01022 " (A) Is a regular 3D (no time) dataset; or \n"
01023 " (B) Is a 3D+time dataset with a sub-brick index specified ('-b3'); or\n"
01024 " (C) Is a 3D+time dataset with no sub-brick index specified ('-b'). \n"
01025 " \n"
01026 " If there is at least one case (C) dataset, then the output dataset will\n"
01027 " also be 3D+time; otherwise it will be a 3D dataset with one sub-brick. \n"
01028 " When producing a 3D+time dataset, datasets in case (A) or (B) will be \n"
01029 " treated as if the particular brick being used has the same value at each\n"
01030 " point in time. \n"
01031 " \n"
01032 #ifdef ALLOW_BUCKETS
01033 " Multi-brick 'bucket' datasets may also be used. Note that if multi-brick\n"
01034 " (bucket or 3D+time) datasets are used, the lowest letter dataset will \n"
01035 " serve as the template for the output; that is, '-b fred+tlrc' takes \n"
01036 " precedence over '-c wilma+tlrc'. (The program 3drefit can be used to \n"
01037 " alter the .HEAD parameters of the output dataset, if desired.) \n"
01038 #endif
01039
01040 #ifdef ALLOW_SUBV
01041 " \n"
01042 "------------------------------------------------------------------------\n"
01043 MASTER_HELP_STRING
01044 " \n"
01045 "** WARNING: you cannot combine sub-brick selection of the form \n"
01046 " -b3 bambam+orig (the old method) \n"
01047 " with sub-brick selection of the form \n"
01048 " -b 'bambam+orig[3]' (the new method) \n"
01049 " If you try, the Doom of Mandos will fall upon you! \n"
01050 #endif
01051 " \n"
01052 "------------------------------------------------------------------------\n"
01053 "1D TIME SERIES: \n"
01054 "-------------- \n"
01055 " \n"
01056 " You can also input a '*.1D' time series file in place of a dataset. \n"
01057 " In this case, the value at each spatial voxel at time index n will be \n"
01058 " the same, and will be the n-th value from the time series file. \n"
01059 " At least one true dataset must be input. If all the input datasets \n"
01060 " are 3D (single sub-brick) or are single sub-bricks from multi-brick \n"
01061 " datasets, then the output will be a 'manufactured' 3D+time dataset. \n"
01062 " \n"
01063 " For example, suppose that 'a3D+orig' is a 3D dataset: \n"
01064 " \n"
01065 " 3dcalc -a a3D+orig -b b.1D -expr \"a*b\" \n"
01066 " \n"
01067 " The output dataset will 3D+time with the value at (x,y,z,t) being \n"
01068 " computed by a3D(x,y,z)*b(t). The TR for this dataset will be set \n"
01069 " to 'tstep' seconds -- this could be altered later with program 3drefit.\n"
01070 " Another method to set up the correct timing would be to input an \n"
01071 " unused 3D+time dataset -- 3dcalc will then copy that dataset's time \n"
01072 " information, but simply do not use that dataset's letter in -expr. \n"
01073 " \n"
01074 " If the *.1D file has multiple columns, only the first read will be \n"
01075 " used in this program. You can select a column to be the first by \n"
01076 " using a sub-vector selection of the form 'b.1D[3]', which will \n"
01077 " choose the 4th column (since counting starts at 0). \n"
01078 " \n"
01079 " '{...}' row selectors can also be used - see the output of '1dcat -help'\n"
01080 " for more details on these. Note that if multiple timeseries or 3D+time\n"
01081 " or 3D bucket datasets are input, they must all have the same number of \n"
01082 " points along the 'time' dimension. \n"
01083 " \n"
01084 "------------------------------------------------------------------------\n"
01085 "'1D:' INPUT: \n"
01086 "----------- \n"
01087 " \n"
01088 " You can input a 1D time series 'dataset' directly on the command line, \n"
01089 " without an external file. The 'filename for such input takes the \n"
01090 " general format \n"
01091 " \n"
01092 " '1D:n_1@val_1,n_2@val_2,n_3@val_3,...' \n"
01093 " \n"
01094 " where each 'n_i' is an integer and each 'val_i' is a float. For \n"
01095 " example \n"
01096 " \n"
01097 " -a '1D:5@0,10@1,5@0,10@1,5@0' \n"
01098 " \n"
01099 " specifies that variable 'a' be assigned to a 1D time series of 35, \n"
01100 " alternating in blocks between values 0 and value 1. \n"
01101 " \n"
01102 "------------------------------------------------------------------------\n"
01103 "'I:*.1D' and 'J:*.1D' and 'K:*.1D' INPUT: \n"
01104 "---------------------------------------- \n"
01105 " \n"
01106 " You can input a 1D time series 'dataset' to be defined as spatially \n"
01107 " dependent instead of time dependent using a syntax like: \n"
01108 " \n"
01109 " -c I:fred.1D \n"
01110 " \n"
01111 " This indicates that the n-th value from file fred.1D is to be associated\n"
01112 " with the spatial voxel index i=n (respectively j=n and k=n for 'J: and \n"
01113 " K: input dataset names). This technique can be useful if you want to \n"
01114 " scale each slice by a fixed constant; for example: \n"
01115 " \n"
01116 " -a dset+orig -b K:slicefactor.1D -expr 'a*b' \n"
01117 " \n"
01118 " In this example, the '-b' value only varies in the k-index spatial \n"
01119 " direction. \n"
01120 " \n"
01121 "------------------------------------------------------------------------\n"
01122 "COORDINATES and PREDEFINED VALUES: \n"
01123 "--------------------------------- \n"
01124 " \n"
01125 " If you don't use '-x', '-y', or '-z' for a dataset, then the voxel \n"
01126 " spatial coordinates will be loaded into those variables. For example, \n"
01127 " the expression 'a*step(x*x+y*y+z*z-100)' will zero out all the voxels \n"
01128 " inside a 10 mm radius of the origin x=y=z=0. \n"
01129 " \n"
01130 " Similarly, the '-t' value, if not otherwise used by a dataset or *.1D \n"
01131 " input, will be loaded with the voxel time coordinate, as determined \n"
01132 " from the header file created for the OUTPUT. Please note that the units\n"
01133 " of this are variable; they might be in milliseconds, seconds, or Hertz.\n"
01134 " In addition, slices of the dataset might be offset in time from one \n"
01135 " another, and this is allowed for in the computation of 't'. Use program\n"
01136 " 3dinfo to find out the structure of your datasets, if you are not sure.\n"
01137 " If no input datasets are 3D+time, then the effective value of TR is \n"
01138 " tstep in the output dataset, with t=0 at the first sub-brick. \n"
01139 " \n"
01140 " Similarly, the '-i', '-j', and '-k' values, if not otherwise used, \n"
01141 " will be loaded with the voxel spatial index coordinates. The '-l' \n"
01142 " (letter 'ell') value will be loaded with the temporal index coordinate.\n"
01143 " \n"
01144 " Otherwise undefined letters will be set to zero. In the future, \n"
01145 " new default values for other letters may be added. \n"
01146 " \n"
01147 " NOTE WELL: By default, the coordinate order of (x,y,z) is the order in \n"
01148 " ********* which the data array is stored on disk; this order is output\n"
01149 " by 3dinfo. The options below control can change this order:\n"
01150 " \n"
01151 " -dicom }= Sets the coordinates to appear in DICOM standard (RAI) order,\n"
01152 " -RAI }= (the AFNI standard), so that -x=Right, -y=Anterior , -z=Inferior,\n"
01153 " +x=Left , +y=Posterior, +z=Superior.\n"
01154 " \n"
01155 " -SPM }= Sets the coordinates to appear in SPM (LPI) order, \n"
01156 " -LPI }= so that -x=Left , -y=Posterior, -z=Inferior,\n"
01157 " +x=Right, +y=Anterior , +z=Superior.\n"
01158 " \n"
01159 "------------------------------------------------------------------------\n"
01160 "DIFFERENTIAL SUBSCRIPTS [22 Nov 1999]: \n"
01161 "----------------------- \n"
01162 " \n"
01163 " Normal calculations with 3dcalc are strictly on a per-voxel basis:\n"
01164 " there is no 'cross-talk' between spatial or temporal locations.\n"
01165 " The differential subscript feature allows you to specify variables\n"
01166 " that refer to different locations, relative to the base voxel.\n"
01167 " For example,\n"
01168 " -a fred+orig -b 'a[1,0,0,0]' -c 'a[0,-1,0,0]' -d 'a[0,0,2,0]'\n"
01169 " means: symbol 'a' refers to a voxel in dataset fred+orig,\n"
01170 " symbol 'b' refers to the following voxel in the x-direction,\n"
01171 " symbol 'c' refers to the previous voxel in the y-direction\n"
01172 " symbol 'd' refers to the 2nd following voxel in the z-direction\n"
01173 "\n"
01174 " To use this feature, you must define the base dataset (e.g., 'a')\n"
01175 " first. Then the differentially subscripted symbols are defined\n"
01176 " using the base dataset symbol followed by 4 integer subscripts,\n"
01177 " which are the shifts in the x-, y-, z-, and t- (or sub-brick index)\n"
01178 " directions. For example,\n"
01179 "\n"
01180 " -a fred+orig -b 'a[0,0,0,1]' -c 'a[0,0,0,-1]' -expr 'median(a,b,c)'\n"
01181 "\n"
01182 " will produce a temporal median smoothing of a 3D+time dataset (this\n"
01183 " can be done more efficiently with program 3dTsmooth).\n"
01184 "\n"
01185 " Note that the physical directions of the x-, y-, and z-axes depend\n"
01186 " on how the dataset was acquired or constructed. See the output of\n"
01187 " program 3dinfo to determine what direction corresponds to what axis.\n"
01188 "\n"
01189 " For convenience, the following abbreviations may be used in place of\n"
01190 " some common subscript combinations:\n"
01191 "\n"
01192 " [1,0,0,0] == +i [-1, 0, 0, 0] == -i\n"
01193 " [0,1,0,0] == +j [ 0,-1, 0, 0] == -j\n"
01194 " [0,0,1,0] == +k [ 0, 0,-1, 0] == -k\n"
01195 " [0,0,0,1] == +l [ 0, 0, 0,-1] == -l\n"
01196 "\n"
01197 " The median smoothing example can thus be abbreviated as\n"
01198 "\n"
01199 " -a fred+orig -b a+l -c a-l -expr 'median(a,b,c)'\n"
01200 "\n"
01201 " When a shift calls for a voxel that is outside of the dataset range,\n"
01202 " one of three things can happen:\n"
01203 "\n"
01204 " STOP => shifting stops at the edge of the dataset\n"
01205 " WRAP => shifting wraps back to the opposite edge of the dataset\n"
01206 " ZERO => the voxel value is returned as zero\n"
01207 "\n"
01208 " Which one applies depends on the setting of the shifting mode at the\n"
01209 " time the symbol using differential subscripting is defined. The mode\n"
01210 " is set by one of the switches '-dsSTOP', '-dsWRAP', or '-dsZERO'. The\n"
01211 " default mode is STOP. Suppose that a dataset has range 0..99 in the\n"
01212 " x-direction. Then when voxel 101 is called for, the value returned is\n"
01213 "\n"
01214 " STOP => value from voxel 99 [didn't shift past edge of dataset]\n"
01215 " WRAP => value from voxel 1 [wrapped back through opposite edge]\n"
01216 " ZERO => the number 0.0 \n"
01217 "\n"
01218 " You can set the shifting mode more than once - the most recent setting\n"
01219 " on the command line applies when a differential subscript symbol is\n"
01220 " encountered.\n"
01221 "\n"
01222 "------------------------------------------------------------------------\n"
01223 "PROBLEMS:\n"
01224 "-------- \n"
01225 "\n"
01226 " * Complex-valued datasets cannot be processed.\n"
01227 " * This program is not very efficient (but is faster than it once was).\n"
01228 " * Differential subscripts slow the program down even more.\n"
01229 "\n"
01230 "------------------------------------------------------------------------\n"
01231 "EXPRESSIONS:\n"
01232 "----------- \n"
01233 "\n"
01234 " Arithmetic expressions are allowed, using + - * / ** and parentheses.\n"
01235 " As noted above, datasets are referred to by single letter variable names.\n"
01236 " At this time, C relational, boolean, and conditional expressions are\n"
01237 " NOT implemented. Built in functions include:\n"
01238 "\n"
01239 " sin , cos , tan , asin , acos , atan , atan2, \n"
01240 " sinh , cosh , tanh , asinh , acosh , atanh , exp , \n"
01241 " log , log10, abs , int , sqrt , max , min , \n"
01242 " J0 , J1 , Y0 , Y1 , erf , erfc , qginv, qg , \n"
01243 " rect , step , astep, bool , and , or , mofn , \n"
01244 " sind , cosd , tand , median, lmode , hmode , mad , \n"
01245 " gran , uran , iran , eran , lran , orstat, \n"
01246 " mean , stdev, sem , Pleg\n"
01247 "\n"
01248 " where:\n"
01249 " * qg(x) = reversed cdf of a standard normal distribution\n"
01250 " * qginv(x) = inverse function to qg\n"
01251 " * min, max, atan2 each take 2 arguments ONLY\n"
01252 " * J0, J1, Y0, Y1 are Bessel functions (see Watson)\n"
01253 " * Pleg(m,x) is the m'th Legendre polynomial evaluated at x\n"
01254 " * erf, erfc are the error and complementary error functions\n"
01255 " * sind, cosd, tand take arguments in degrees (vs. radians)\n"
01256 " * median(a,b,c,...) computes the median of its arguments\n"
01257 " * mad(a,b,c,...) computes the MAD of its arguments\n"
01258 " * mean(a,b,c,...) computes the mean of its arguments\n"
01259 " * stdev(a,b,c,...) computes the standard deviation of its arguments\n"
01260 " * sem(a,b,c,...) computes the standard error of the mean of its arguments,\n"
01261 " where sem(n arguments) = stdev(same)/sqrt(n)\n"
01262 " * orstat(n,a,b,c,...) computes the n-th order statistic of\n"
01263 " {a,b,c,...} - that is, the n-th value in size, starting\n"
01264 " at the bottom (e.g., orstat(1,a,b,c) is the minimum)\n"
01265 " * lmode(a,b,c,...) and hmode(a,b,c,...) compute the mode\n"
01266 " of their arguments - lmode breaks ties by choosing the\n"
01267 " smallest value with the maximal count, hmode breaks ties by\n"
01268 " choosing the largest value with the maximal count\n"
01269 " [median,lmode,hmode take a variable number of arguments]\n"
01270 " * gran(m,s) returns a Gaussian deviate with mean=m, stdev=s\n"
01271 " * uran(r) returns a uniform deviate in the range [0,r]\n"
01272 " * iran(t) returns a random integer in the range [0..t]\n"
01273 " * eran(s) returns an exponentially distributed deviate\n"
01274 " * lran(t) returns a logistically distributed deviate\n"
01275 "\n"
01276 " You may use the symbol 'PI' to refer to the constant of that name.\n"
01277 " This is the only 2 letter symbol defined; all input files are\n"
01278 " referred to by 1 letter symbols. The case of the expression is\n"
01279 " ignored (in fact, it is converted to uppercase as the first step\n"
01280 " in the parsing algorithm).\n"
01281 "\n"
01282 " The following functions are designed to help implement logical\n"
01283 " functions, such as masking of 3D volumes against some criterion:\n"
01284 " step(x) = {1 if x>0 , 0 if x<=0},\n"
01285 " astep(x,y) = {1 if abs(x) > y , 0 otherwise} = step(abs(x)-y)\n"
01286 " rect(x) = {1 if abs(x)<=0.5, 0 if abs(x)>0.5},\n"
01287 " bool(x) = {1 if x != 0.0 , 0 if x == 0.0},\n"
01288 " notzero(x) = bool(x),\n"
01289 " iszero(x) = 1-bool(x) = { 0 if x != 0.0, 1 if x == 0.0 },\n"
01290 " equals(x,y) = 1-bool(x-y) = { 1 if x == y , 0 if x != y },\n"
01291 " ispositive(x) = { 1 if x > 0; 0 if x <= 0 },\n"
01292 " isnegative(x) = { 1 if x < 0; 0 if x >= 0 },\n"
01293 " and(a,b,...,c) = {1 if all arguments are nonzero, 0 if any are zero}\n"
01294 " or(a,b,...,c) = {1 if any arguments are nonzero, 0 if all are zero}\n"
01295 " mofn(m,a,...,c) = {1 if at least 'm' arguments are nonzero, 0 otherwise}\n"
01296 " argmax(a,b,...) = index of largest argument; = 0 if all args are 0\n"
01297 " argnum(a,b,...) = number of nonzero arguments\n"
01298 "\n"
01299 " [These last 5 functions take a variable number of arguments.]\n"
01300 "\n"
01301 " The following 27 new [Mar 1999] functions are used for statistical\n"
01302 " conversions, as in the program 'cdf':\n"
01303 " fico_t2p(t,a,b,c), fico_p2t(p,a,b,c), fico_t2z(t,a,b,c),\n"
01304 " fitt_t2p(t,a) , fitt_p2t(p,a) , fitt_t2z(t,a) ,\n"
01305 " fift_t2p(t,a,b) , fift_p2t(p,a,b) , fift_t2z(t,a,b) ,\n"
01306 " fizt_t2p(t) , fizt_p2t(p) , fizt_t2z(t) ,\n"
01307 " fict_t2p(t,a) , fict_p2t(p,a) , fict_t2z(t,a) ,\n"
01308 " fibt_t2p(t,a,b) , fibt_p2t(p,a,b) , fibt_t2z(t,a,b) ,\n"
01309 " fibn_t2p(t,a,b) , fibn_p2t(p,a,b) , fibn_t2z(t,a,b) ,\n"
01310 " figt_t2p(t,a,b) , figt_p2t(p,a,b) , figt_t2z(t,a,b) ,\n"
01311 " fipt_t2p(t,a) , fipt_p2t(p,a) , fipt_t2z(t,a) .\n"
01312 "\n"
01313 " See the output of 'cdf -help' for documentation on the meanings of\n"
01314 " and arguments to these functions. (After using one of these, you\n"
01315 " may wish to use program '3drefit' to modify the dataset statistical\n"
01316 " auxiliary parameters.)\n"
01317 "\n"
01318 " Computations are carried out in double precision before being\n"
01319 " truncated to the final output 'datum'.\n"
01320 "\n"
01321 " Note that the quotes around the expression are needed so the shell\n"
01322 " doesn't try to expand * characters, or interpret parentheses.\n"
01323 "\n"
01324 " (Try the 'ccalc' program to see how the expression evaluator works.\n"
01325 " The arithmetic parser and evaluator is written in Fortran-77 and\n"
01326 " is derived from a program written long ago by RW Cox to facilitate\n"
01327 " compiling on an array processor hooked up to a VAX. It's a mess,\n"
01328 " but it works - somewhat slowly.)\n"
01329 ) ;
01330 exit(0) ;
01331 }
01332
01333
01334
01335 int main( int argc , char * argv[] )
01336 {
01337 #define VSIZE 1024
01338
01339 double * atoz[26] ;
01340 int ii , ids , jj, kk, kt, ll, jbot, jtop ;
01341 THD_3dim_dataset * new_dset=NULL ;
01342 float ** buf;
01343 double temp[VSIZE];
01344 int nbad ;
01345
01346 THD_ivec3 iv ;
01347 THD_fvec3 fv ;
01348 float xxx[VSIZE], yyy[VSIZE], zzz[VSIZE] ;
01349 int iii,jjj,kkk , nx,nxy ;
01350 THD_dataxes * daxes ;
01351
01352
01353
01354 if( argc < 2 || strncasecmp(argv[1],"-help",4) == 0 ) CALC_Syntax() ;
01355
01356
01357
01358 mainENTRY("3dcalc main"); machdep() ; PRINT_VERSION("3dcalc") ;
01359
01360 { int new_argc ; char ** new_argv ;
01361 addto_args( argc , argv , &new_argc , &new_argv ) ;
01362 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
01363 }
01364
01365 AFNI_logger("3dcalc",argc,argv) ;
01366
01367 for (ii=0; ii<26; ii++) ntime[ii] = 0 ;
01368
01369 CALC_read_opts( argc , argv ) ;
01370
01371
01372
01373 if( ntime_max == 1 || TS_make == 1 ){
01374 for( ids=0 ; ids < 26 ; ids++ ) if( CALC_dset[ids] != NULL ) break ;
01375 } else {
01376 for( ids=0 ; ids < 26 ; ids++ ) if( CALC_dset[ids] != NULL &&
01377 ntime[ids] > 1 ) break ;
01378 }
01379 if( ids == 26 ) ERROR_exit("Can't find template dataset?!\n") ;
01380
01381 new_dset = EDIT_empty_copy( CALC_dset[ids] ) ;
01382
01383
01384
01385 for( iii=0 ; iii < 26 ; iii++ ){
01386 if( iii != ids &&
01387 CALC_dset[iii] != NULL &&
01388 !EQUIV_DATAXES(new_dset->daxes,CALC_dset[iii]->daxes) )
01389 WARNING_message("dataset '%c'=%s grid mismatch with %s\n",
01390 'a'+iii , DSET_BRIKNAME(CALC_dset[iii]) ,
01391 DSET_BRIKNAME(CALC_dset[ids]) ) ;
01392 }
01393
01394
01395
01396 if( CALC_histpar < 0 ){
01397 for( iii=jjj=0 ; iii < 26 ; iii++ )
01398 if( CALC_dset[iii] != NULL ) jjj++ ;
01399 } else {
01400 ids = CALC_histpar ;
01401 jjj = 1 ;
01402 }
01403
01404 if( jjj == 1 ){
01405 tross_Copy_History( CALC_dset[ids] , new_dset ) ;
01406 } else {
01407 char hbuf[64] ;
01408 tross_Append_History( new_dset ,
01409 "===================================" ) ;
01410 tross_Append_History( new_dset ,
01411 "=== History of inputs to 3dcalc ===" ) ;
01412 for( iii=0 ; iii < 26 ; iii++ ){
01413 if( CALC_dset[iii] != NULL ){
01414 sprintf(hbuf,"=== Input %c:", 'a'+iii ) ;
01415 tross_Append_History( new_dset , hbuf ) ;
01416 tross_Addto_History( CALC_dset[iii] , new_dset ) ;
01417 }
01418 }
01419 tross_Append_History( new_dset ,
01420 "===================================" ) ;
01421 }
01422 tross_Make_History( "3dcalc" , argc,argv , new_dset ) ;
01423
01424 if( CALC_datum < 0 ) CALC_datum = MRI_float ;
01425
01426 EDIT_dset_items( new_dset ,
01427 ADN_prefix , CALC_output_prefix ,
01428 ADN_directory_name , CALC_session ,
01429 ADN_datum_all , CALC_datum ,
01430 ADN_none ) ;
01431
01432 if( DSET_NVALS(new_dset) != ntime_max )
01433 EDIT_dset_items( new_dset , ADN_nvals , ntime_max , ADN_none ) ;
01434
01435
01436
01437
01438 if( TS_make ){
01439 EDIT_dset_items( new_dset ,
01440 ADN_ntt , ntime_max ,
01441 ADN_ttdel , TS_dt ,
01442 ADN_ttorg , 0.0 ,
01443 ADN_ttdur , 0.0 ,
01444 ADN_tunits , UNITS_SEC_TYPE ,
01445 ADN_none ) ;
01446 }
01447
01448 if( ISFUNC(new_dset) && ! ISFUNCBUCKET(new_dset) && new_dset->taxis != NULL )
01449 EDIT_dset_items( new_dset , ADN_func_type , FUNC_FIM_TYPE , ADN_none ) ;
01450 else if( ISANATBUCKET(new_dset) )
01451 EDIT_dset_items( new_dset , ADN_func_type , ANAT_EPI_TYPE , ADN_none ) ;
01452
01453 if( THD_is_file(new_dset->dblk->diskptr->header_name) )
01454 ERROR_exit("Output file %s already exists -- cannot continue!\n",
01455 new_dset->dblk->diskptr->header_name ) ;
01456
01457 for (ids=0; ids<26; ids++)
01458 atoz[ids] = (double *) malloc(sizeof(double) * VSIZE ) ;
01459
01460 for( ids=0 ; ids < 26 ; ids++ )
01461 for (ii=0; ii<VSIZE; ii++)
01462 atoz[ids][ii] = 0.0 ;
01463
01464
01465
01466 nx = DSET_NX(new_dset) ;
01467 nxy = nx * DSET_NY(new_dset) ; daxes = new_dset->daxes ;
01468
01469 buf = (float **) malloc(sizeof(float *) * ntime_max);
01470
01471 for ( kt = 0 ; kt < ntime_max ; kt ++ ) {
01472
01473 if( CALC_verbose )
01474 INFO_message("Computing sub-brick %d\n",kt) ;
01475
01476
01477
01478 buf[kt] = (float *) malloc(sizeof(float) * CALC_nvox);
01479 if( buf[kt] == NULL )
01480 ERROR_exit("Can't malloc output dataset sub-brick %d!\n",kt) ;
01481
01482
01483
01484 for ( ii = 0 ; ii < CALC_nvox ; ii += VSIZE ) {
01485
01486 jbot = ii ;
01487 jtop = MIN( ii + VSIZE , CALC_nvox ) ;
01488
01489
01490
01491 if( CALC_has_xyz ){
01492 for( jj=jbot ; jj < jtop ; jj++ ){
01493 LOAD_IVEC3( iv , jj%nx , (jj%nxy)/nx , jj/nxy ) ;
01494 fv = THD_3dind_to_3dmm( new_dset , iv ) ;
01495 if( CALC_mangle_xyz )
01496 fv = THD_3dmm_to_dicomm(new_dset,fv) ;
01497 UNLOAD_FVEC3(fv,xxx[jj-jbot],yyy[jj-jbot],zzz[jj-jbot]) ;
01498 if( CALC_mangle_xyz == MANGLE_LPI ){
01499 xxx[jj-jbot] = -xxx[jj-jbot] ; yyy[jj-jbot] = -yyy[jj-jbot] ;
01500 }
01501 }
01502 }
01503
01504
01505
01506 for (ids = 0 ; ids < 26 ; ids ++ ) {
01507
01508
01509
01510
01511 if( TS_flim[ids] != NULL ){
01512 if( jbot == 0 ){
01513 double tval ;
01514 if( kt < TS_flim[ids]->nx ) tval = TS_flar[ids][kt] ;
01515 else tval = 0.0 ;
01516
01517 for (jj =jbot ; jj < jtop ; jj ++ )
01518 atoz[ids][jj-ii] = tval ;
01519 }
01520 }
01521
01522
01523
01524 else if( IJKAR_flim[ids] != NULL ){
01525 int ss , ix=IJKAR_flim[ids]->nx ;
01526
01527 switch( IJKAR_dcod[ids] ){
01528 case 8:
01529 for( jj=jbot ; jj < jtop ; jj++ ){
01530 ss = (jj%nx) ;
01531 atoz[ids][jj-jbot] = (ss < ix) ? IJKAR_flar[ids][ss] : 0.0 ;
01532 }
01533 break ;
01534 case 9:
01535 for( jj=jbot ; jj < jtop ; jj++ ){
01536 ss = ((jj%nxy)/nx) ;
01537 atoz[ids][jj-jbot] = (ss < ix) ? IJKAR_flar[ids][ss] : 0.0 ;
01538 }
01539 break ;
01540 case 10:
01541 for( jj=jbot ; jj < jtop ; jj++ ){
01542 ss = (jj/nxy) ;
01543 atoz[ids][jj-jbot] = (ss < ix) ? IJKAR_flar[ids][ss] : 0.0 ;
01544 }
01545 break ;
01546 }
01547 }
01548
01549
01550
01551 else if( CALC_dshift[ids] >= 0 ){
01552 int jds = CALC_dshift[ids] ;
01553 int kts , jjs , ix,jy,kz ;
01554 int id=CALC_dshift_i[ids] , jd=CALC_dshift_j[ids] ,
01555 kd=CALC_dshift_k[ids] , ld=CALC_dshift_l[ids] ;
01556 int ijkd = ((id!=0) || (jd!=0) || (kd!=0)) ;
01557 int dsx = DSET_NX(CALC_dset[jds]) - 1 ;
01558 int dsy = DSET_NY(CALC_dset[jds]) - 1 ;
01559 int dsz = DSET_NZ(CALC_dset[jds]) - 1 ;
01560 int dst = ntime[jds] - 1 ;
01561 int mode = CALC_dshift_mode[ids] , dun=0 ;
01562
01563 kts = kt + ld ;
01564 if( kts < 0 || kts > dst ){
01565 switch( mode ){
01566 case DSHIFT_MODE_ZERO:
01567 for( jj=jbot ; jj < jtop ; jj++ ) atoz[ids][jj-ii] = 0.0 ;
01568 dun = 1 ;
01569 break ;
01570 default:
01571 case DSHIFT_MODE_STOP:
01572 if( kts < 0 ) kts = 0 ;
01573 else if( kts > dst ) kts = dst ;
01574 break ;
01575 case DSHIFT_MODE_WRAP:
01576 while( kts < 0 ) kts += (dst+1) ;
01577 while( kts > dst ) kts -= (dst+1) ;
01578 break ;
01579 }
01580 }
01581
01582 if( !dun ){
01583 for( dun=0,jj=jbot ; jj < jtop ; jj++ ){
01584 jjs = jj ;
01585 if( ijkd ){
01586 ix = DSET_index_to_ix(CALC_dset[jds],jj) ;
01587 jy = DSET_index_to_jy(CALC_dset[jds],jj) ;
01588 kz = DSET_index_to_kz(CALC_dset[jds],jj) ;
01589
01590 ix += id ;
01591 if( ix < 0 || ix > dsx ){
01592 switch( mode ){
01593 case DSHIFT_MODE_ZERO:
01594 atoz[ids][jj-ii] = 0.0 ; dun = 1 ;
01595 break ;
01596 default:
01597 case DSHIFT_MODE_STOP:
01598 if( ix < 0 ) ix = 0 ;
01599 else if( ix > dsx ) ix = dsx ;
01600 break ;
01601 case DSHIFT_MODE_WRAP:
01602 while( ix < 0 ) ix += (dsx+1) ;
01603 while( ix > dsx ) ix -= (dsx+1) ;
01604 break ;
01605 }
01606 }
01607 if( dun ){ dun=0; continue; }
01608
01609 jy += jd ;
01610 if( jy < 0 || jy > dsy ){
01611 switch( mode ){
01612 case DSHIFT_MODE_ZERO:
01613 atoz[ids][jj-ii] = 0.0 ; dun = 1 ;
01614 break ;
01615 default:
01616 case DSHIFT_MODE_STOP:
01617 if( jy < 0 ) jy = 0 ;
01618 else if( jy > dsy ) jy = dsy ;
01619 break ;
01620 case DSHIFT_MODE_WRAP:
01621 while( jy < 0 ) jy += (dsy+1) ;
01622 while( jy > dsy ) jy -= (dsy+1) ;
01623 break ;
01624 }
01625 }
01626 if( dun ){ dun=0; continue; }
01627
01628 kz += kd ;
01629 if( kz < 0 || kz > dsz ){
01630 switch( mode ){
01631 case DSHIFT_MODE_ZERO:
01632 atoz[ids][jj-ii] = 0.0 ; dun = 1 ;
01633 break ;
01634 default:
01635 case DSHIFT_MODE_STOP:
01636 if( kz < 0 ) kz = 0 ;
01637 else if( kz > dsz ) kz = dsz ;
01638 break ;
01639 case DSHIFT_MODE_WRAP:
01640 while( kz < 0 ) kz += (dsz+1) ;
01641 while( kz > dsz ) kz -= (dsz+1) ;
01642 break ;
01643 }
01644 }
01645 if( dun ){ dun=0; continue; }
01646
01647 jjs = DSET_ixyz_to_index(CALC_dset[jds],ix,jy,kz) ;
01648 }
01649 switch( CALC_type[jds] ) {
01650 case MRI_short:
01651 atoz[ids][jj-ii] = CALC_short[jds][kts][jjs]
01652 * CALC_ffac[jds][kts];
01653 break ;
01654 case MRI_float:
01655 atoz[ids][jj-ii] = CALC_float[jds][kts][jjs]
01656 * CALC_ffac[jds][kts];
01657 break ;
01658 case MRI_byte:
01659 atoz[ids][jj-ii] = CALC_byte[jds][kts][jjs]
01660 * CALC_ffac[jds][kts];
01661 break ;
01662 case MRI_rgb:
01663 atoz[ids][jj-ii] = Rfac*CALC_byte[jds][kts][3*jjs ]
01664 +Gfac*CALC_byte[jds][kts][3*jjs+1]
01665 +Bfac*CALC_byte[jds][kts][3*jjs+2] ;
01666 break ;
01667 }
01668 }
01669 }
01670 }
01671
01672
01673
01674 else if ( ntime[ids] == 1 && CALC_type[ids] >= 0 ) {
01675 switch( CALC_type[ids] ) {
01676 case MRI_short:
01677 if( CALC_noffac[ids] )
01678 for (jj =jbot ; jj < jtop ; jj ++ )
01679 atoz[ids][jj-ii] = CALC_short[ids][0][jj] ;
01680 else
01681 for (jj =jbot ; jj < jtop ; jj ++ )
01682 atoz[ids][jj-ii] = CALC_short[ids][0][jj] * CALC_ffac[ids][0] ;
01683 break;
01684
01685 case MRI_float:
01686 if( CALC_noffac[ids] )
01687 for (jj =jbot ; jj < jtop ; jj ++ )
01688 atoz[ids][jj-ii] = CALC_float[ids][0][jj] ;
01689 else
01690 for (jj =jbot ; jj < jtop ; jj ++ )
01691 atoz[ids][jj-ii] = CALC_float[ids][0][jj] * CALC_ffac[ids][0] ;
01692 break;
01693
01694 case MRI_byte:
01695 if( CALC_noffac[ids] )
01696 for (jj =jbot ; jj < jtop ; jj ++ )
01697 atoz[ids][jj-ii] = CALC_byte[ids][0][jj] ;
01698 else
01699 for (jj =jbot ; jj < jtop ; jj ++ )
01700 atoz[ids][jj-ii] = CALC_byte[ids][0][jj] * CALC_ffac[ids][0] ;
01701 break;
01702
01703 case MRI_rgb:
01704 for (jj =jbot ; jj < jtop ; jj ++ )
01705 atoz[ids][jj-ii] = Rfac*CALC_byte[ids][0][3*jj ]
01706 +Gfac*CALC_byte[ids][0][3*jj+1]
01707 +Bfac*CALC_byte[ids][0][3*jj+2] ;
01708 break;
01709 }
01710 }
01711
01712
01713
01714 else if( ntime[ids] > 1 && CALC_type[ids] >= 0 ) {
01715 switch ( CALC_type[ids] ) {
01716 case MRI_short:
01717 if( CALC_noffac[ids] )
01718 for (jj = jbot ; jj < jtop ; jj ++ )
01719 atoz[ids][jj-ii] = CALC_short[ids][kt][jj] ;
01720 else
01721 for (jj = jbot ; jj < jtop ; jj ++ )
01722 atoz[ids][jj-ii] = CALC_short[ids][kt][jj] * CALC_ffac[ids][kt];
01723 break;
01724
01725 case MRI_float:
01726 if( CALC_noffac[ids] )
01727 for (jj = jbot ; jj < jtop ; jj ++ )
01728 atoz[ids][jj-ii] = CALC_float[ids][kt][jj] ;
01729 else
01730 for (jj = jbot ; jj < jtop ; jj ++ )
01731 atoz[ids][jj-ii] = CALC_float[ids][kt][jj] * CALC_ffac[ids][kt];
01732 break;
01733
01734 case MRI_byte:
01735 if( CALC_noffac[ids] )
01736 for (jj = jbot ; jj < jtop ; jj ++ )
01737 atoz[ids][jj-ii] = CALC_byte[ids][kt][jj] ;
01738 else
01739 for (jj = jbot ; jj < jtop ; jj ++ )
01740 atoz[ids][jj-ii] = CALC_byte[ids][kt][jj] * CALC_ffac[ids][kt];
01741 break;
01742
01743 case MRI_rgb:
01744 for (jj =jbot ; jj < jtop ; jj ++ )
01745 atoz[ids][jj-ii] = Rfac*CALC_byte[ids][kt][3*jj ]
01746 +Gfac*CALC_byte[ids][kt][3*jj+1]
01747 +Bfac*CALC_byte[ids][kt][3*jj+2] ;
01748 break;
01749 }
01750 }
01751
01752
01753
01754 else if( CALC_has_predefined ) {
01755
01756 switch( ids ){
01757 case 23:
01758 if( HAS_X )
01759 for( jj=jbot ; jj < jtop ; jj++ )
01760 atoz[ids][jj-ii] = xxx[jj-ii] ;
01761 break ;
01762
01763 case 24:
01764 if( HAS_Y )
01765 for( jj=jbot ; jj < jtop ; jj++ )
01766 atoz[ids][jj-ii] = yyy[jj-ii] ;
01767 break ;
01768
01769 case 25:
01770 if( HAS_Z )
01771 for( jj=jbot ; jj < jtop ; jj++ )
01772 atoz[ids][jj-ii] = zzz[jj-ii] ;
01773 break ;
01774
01775 case 8:
01776 if( HAS_I )
01777 for( jj=jbot ; jj < jtop ; jj++ )
01778 atoz[ids][jj-ii] = (jj%nx) ;
01779 break ;
01780
01781 case 9:
01782 if( HAS_J )
01783 for( jj=jbot ; jj < jtop ; jj++ )
01784 atoz[ids][jj-ii] = ((jj%nxy)/nx) ;
01785 break ;
01786
01787 case 10:
01788 if( HAS_K )
01789 for( jj=jbot ; jj < jtop ; jj++ )
01790 atoz[ids][jj-ii] = (jj/nxy) ;
01791 break ;
01792
01793 case 19:
01794 if( HAS_T )
01795 for( jj=jbot ; jj < jtop ; jj++ )
01796 atoz[ids][jj-ii] = THD_timeof_vox(kt,jj,new_dset) ;
01797 break ;
01798
01799 case 11:
01800 if( HAS_L )
01801 for( jj=jbot ; jj < jtop ; jj++ )
01802 atoz[ids][jj-ii] = kt ;
01803 break ;
01804 }
01805
01806 }
01807 }
01808
01809
01810
01811 PARSER_evaluate_vector(CALC_code, atoz, jtop-jbot, temp);
01812 for ( jj = jbot ; jj < jtop ; jj ++ )
01813 buf[kt][jj] = temp[jj-ii];
01814
01815 }
01816
01817
01818
01819 nbad = thd_floatscan( CALC_nvox , buf[kt] ) ;
01820 if( nbad > 0 )
01821 WARNING_message("%d bad floats replaced by 0 in sub-brick %d\n\a",
01822 nbad , kt ) ;
01823
01824
01825
01826 if( ! CALC_has_timeshift ){
01827 for( ids=0 ; ids < 26 ; ids ++ ){
01828 if( CALC_dset[ids] != NULL && ntime[ids] > 1 &&
01829 CALC_dset[ids]->dblk->malloc_type == DATABLOCK_MEM_MALLOC ){
01830
01831 void * ptr = DSET_ARRAY(CALC_dset[ids],kt) ;
01832 if( ptr != NULL ) free(ptr) ;
01833 mri_clear_data_pointer( DSET_BRICK(CALC_dset[ids],kt) ) ;
01834 }
01835 }
01836 }
01837
01838 }
01839
01840 for( ids=0 ; ids < 26 ; ids++ ){
01841 if( CALC_dset[ids] != NULL ) PURGE_DSET( CALC_dset[ids] ) ;
01842 if( TS_flim[ids] != NULL ) mri_free( TS_flim[ids] ) ;
01843 if( IJKAR_flim[ids]!= NULL ) mri_free( IJKAR_flim[ids] ) ;
01844 }
01845
01846
01847
01848 switch( CALC_datum ){
01849
01850 default:
01851 ERROR_exit("Somehow ended up with CALC_datum = %d\n",CALC_datum) ;
01852 exit(1) ;
01853
01854
01855
01856 case MRI_float:{
01857 for( ii=0 ; ii < ntime_max ; ii++ ){
01858 EDIT_substitute_brick(new_dset, ii, MRI_float, buf[ii]);
01859 DSET_BRICK_FACTOR(new_dset, ii) = 0.0;
01860 }
01861 }
01862 break ;
01863
01864
01865
01866 case MRI_byte:
01867 case MRI_short:{
01868 void ** dfim ;
01869 float gtop , fimfac , gtemp ;
01870
01871 if( CALC_verbose )
01872 INFO_message("Scaling output to type %s brick(s)\n",
01873 MRI_TYPE_name[CALC_datum] ) ;
01874
01875 dfim = (void ** ) malloc( sizeof( void * ) * ntime_max ) ;
01876
01877 if( CALC_gscale ){
01878 gtop = 0.0 ;
01879 for( ii=0 ; ii < ntime_max ; ii++ ){
01880 gtemp = MCW_vol_amax( CALC_nvox , 1 , 1 , MRI_float, buf[ii] ) ;
01881 gtop = MAX( gtop , gtemp ) ;
01882 if( gtemp == 0.0 )
01883 WARNING_message("output sub-brick %d is all zeros!\n",ii) ;
01884 }
01885 }
01886
01887 for (ii = 0 ; ii < ntime_max ; ii ++ ) {
01888
01889
01890
01891 if( ! CALC_gscale ){
01892 gtop = MCW_vol_amax( CALC_nvox , 1 , 1 , MRI_float, buf[ii] ) ;
01893 if( gtop == 0.0 )
01894 WARNING_message("output sub-brick %d is all zeros!\n",ii) ;
01895 }
01896
01897
01898
01899 if( CALC_fscale ){
01900 fimfac = (gtop > 0.0) ? MRI_TYPE_maxval[CALC_datum] / gtop : 0.0 ;
01901
01902 } else if( !CALC_nscale ){
01903
01904 fimfac = (gtop > MRI_TYPE_maxval[CALC_datum] || (gtop > 0.0 && gtop <= 1.0) )
01905 ? MRI_TYPE_maxval[CALC_datum]/ gtop : 0.0 ;
01906
01907 if( fimfac == 0.0 && gtop > 0.0 ){
01908 float fv,iv ; int kk ;
01909 for( kk=0 ; kk < CALC_nvox ; kk++ ){
01910 fv = buf[ii][kk] ; iv = rint(fv) ;
01911 if( fabs(fv-iv) >= 0.01 ){
01912 fimfac = MRI_TYPE_maxval[CALC_datum]/ gtop ; break ;
01913 }
01914 }
01915 }
01916
01917 } else {
01918 fimfac = 0.0 ;
01919 }
01920
01921 if( CALC_verbose ){
01922 if( fimfac != 0.0 )
01923 INFO_message("Sub-brick %d scale factor = %f\n",ii,fimfac) ;
01924 else
01925 INFO_message("Sub-brick %d: no scale factor\n" ,ii) ;
01926 }
01927
01928
01929
01930 dfim[ii] = (void *) malloc( mri_datum_size(CALC_datum) * CALC_nvox ) ;
01931 if( dfim[ii] == NULL ) ERROR_exit("malloc fails at output[%d]\n",ii);
01932
01933 if( CALC_datum == MRI_byte ){
01934 int nneg ;
01935 for( nneg=jj=0 ; jj < CALC_nvox ; jj++ ) nneg += (buf[ii][jj] < 0.0f) ;
01936 if( nneg > 0 )
01937 WARNING_message(
01938 "sub-brick #%d has %d negative values set=0 in conversion to bytes\n",
01939 ii , nneg ) ;
01940 }
01941
01942 EDIT_coerce_scale_type( CALC_nvox , fimfac ,
01943 MRI_float, buf[ii] , CALC_datum,dfim[ii] ) ;
01944 free( buf[ii] ) ;
01945
01946
01947
01948 EDIT_substitute_brick(new_dset, ii, CALC_datum, dfim[ii] );
01949
01950 DSET_BRICK_FACTOR(new_dset,ii) = (fimfac != 0.0) ? 1.0/fimfac : 0.0 ;
01951 }
01952 }
01953 break ;
01954 }
01955
01956 if( CALC_verbose ) INFO_message("Computing output statistics\n") ;
01957 THD_load_statistics( new_dset ) ;
01958
01959 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01960 if( CALC_verbose ) WROTE_DSET(new_dset) ;
01961
01962 exit(0) ;
01963 }
