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 #define PROGRAM_NAME "3dmerge"
00032 #define LAST_MOD_DATE "02 Nov 2001"
00033
00034 #include "mrilib.h"
00035 #include "parser.h"
00036
00037 #define MAIN
00038 #define MEGA 1048576
00039
00040 #ifndef myXtFree
00041 #define myXtFree(xp) (XtFree((char *)(xp)) , (xp)=NULL)
00042 #endif
00043
00044 #define ALLOW_SUBV
00045
00046 #ifdef ALLOW_SUBV
00047 # define DSET_OPEN THD_open_dataset
00048 #else
00049 # define DSET_OPEN THD_open_one_dataset
00050 #endif
00051
00052
00053
00054
00055
00056 #define CFLAG_MEAN 1
00057 #define CFLAG_NZMEAN 2
00058 #define CFLAG_MMAX 3
00059 #define CFLAG_COUNT 4
00060 #define CFLAG_AMAX 5
00061 #define CFLAG_SMAX 6
00062 #define CFLAG_ORDER 7
00063 #define CFLAG_FISHER 8
00064
00065 #define THFLAG_NONE 0
00066 #define THFLAG_FICO 71
00067
00068 #define TANH(z) tanh(z)
00069 #define ATANH(z) atanh(z)
00070
00071
00072
00073 static EDIT_options MRG_edopt ;
00074 static int MRG_have_edopt = 0 ;
00075
00076 static int MRG_hits_g = 0 ;
00077 static int MRG_cflag_g = CFLAG_MEAN ;
00078 static int MRG_keepthr = 0 ;
00079 static float MRG_clust_rmm_g = 0.0 ;
00080 static float MRG_clust_vmul_g = 0.0 ;
00081 static int MRG_datum = ILLEGAL_TYPE ;
00082 static int MRG_thdatum = ILLEGAL_TYPE ;
00083 static int MRG_be_quiet = 0 ;
00084 static int MRG_cflag_gthr = THFLAG_NONE ;
00085 static int MRG_doall = 0;
00086 static int MRG_verbose = 0;
00087
00088 static char MRG_output_session[THD_MAX_NAME] = "./" ;
00089 static char MRG_output_prefix [THD_MAX_PREFIX] = "mrg" ;
00090 #if 0
00091 static char MRG_output_label [THD_MAX_LABEL] = "\0" ;
00092 #endif
00093
00094 static int MRG_ivfim = -1 ;
00095 static int MRG_ivthr = -1 ;
00096
00097 static int MRG_nscale = 0 ;
00098
00099
00100 int MRG_read_opts( int , char ** ) ;
00101 void MRG_Syntax(void) ;
00102
00103
00104
00105
00106
00107
00108
00109 #ifdef MERGE_DEBUG
00110 # define DUMP1 fprintf(stderr,"ARG: %s\n",argv[nopt])
00111 # define DUMP2 fprintf(stderr,"ARG: %s %s\n",argv[nopt],argv[nopt+1])
00112 # define DUMP3 fprintf(stderr,"ARG: %s %s %s\n",argv[nopt],argv[nopt+1],argv[nopt+2])
00113 #else
00114 # define DUMP1
00115 # define DUMP2
00116 # define DUMP3
00117 #endif
00118
00119 int MRG_read_opts( int argc , char * argv[] )
00120 {
00121 int nopt = 1 ;
00122 float val ;
00123 int ival ;
00124
00125 INIT_EDOPT( &MRG_edopt ) ;
00126
00127 while( nopt < argc && argv[nopt][0] == '-' ){
00128
00129
00130
00131 ival = EDIT_check_argv( argc , argv , nopt , &MRG_edopt ) ;
00132 if( ival > 0 ){
00133 nopt += ival ; MRG_have_edopt = 1 ;
00134 continue ;
00135 }
00136
00137
00138
00139 if( strncmp(argv[nopt],"-verb",5) == 0 ){
00140 MRG_verbose = MRG_edopt.verbose = 1 ;
00141 nopt++ ; continue ;
00142 }
00143
00144
00145
00146 if( strncmp(argv[nopt],"-1dindex",5) == 0 ){
00147 if( ++nopt >= argc ){
00148 fprintf(stderr,"need an argument after -1dindex!\n") ; exit(1) ;
00149 }
00150 MRG_ivfim = MRG_edopt.iv_fim = (int) strtod( argv[nopt++] , NULL ) ;
00151 continue ;
00152 }
00153
00154 if( strncmp(argv[nopt],"-1tindex",5) == 0 ){
00155 if( ++nopt >= argc ){
00156 fprintf(stderr,"need an argument after -1tindex!\n") ; exit(1) ;
00157 }
00158 MRG_ivthr = MRG_edopt.iv_thr = (int) strtod( argv[nopt++] , NULL ) ;
00159 continue ;
00160 }
00161
00162
00163
00164 if( strncmp(argv[nopt],"-1fmask",6) == 0 ){
00165 THD_3dim_dataset * mset ; int nn ;
00166
00167 if( MRG_edopt.nfmask > 0 ){
00168 fprintf(stderr,"Can't have 2 -fmask options!\n") ;
00169 exit(1) ;
00170 }
00171 if( ++nopt >= argc ){
00172 fprintf(stderr,"No argument after %s?\n",argv[nopt-1]) ;
00173 exit(1);
00174 }
00175
00176 mset = THD_open_dataset( argv[nopt++] ) ;
00177 if( mset == NULL ){
00178 fprintf(stderr,"*** Can't open -fmask dataset\n") ;
00179 exit(1) ;
00180 }
00181 if( DSET_BRICK_TYPE(mset,0) == MRI_complex ){
00182 fprintf(stderr,"*** Can't deal with complex-valued -fmask dataset\n");
00183 exit(1) ;
00184 }
00185
00186 MRG_edopt.fmask = THD_makemask( mset , 0 , 666.0,-666.0 ) ;
00187 MRG_edopt.nfmask = DSET_NVOX(mset) ;
00188
00189 nn = THD_countmask(MRG_edopt.nfmask,MRG_edopt.fmask) ;
00190 if( nn < 2 ){
00191 fprintf(stderr,"*** Too few (%d) nonzero voxels in -fmask!\n",nn) ;
00192 exit(1) ;
00193 }
00194
00195 DSET_delete(mset) ;
00196
00197 if( MRG_edopt.filter_opt == FCFLAG_AVER )
00198 MRG_edopt.filter_opt = FCFLAG_MEAN ;
00199
00200 if( MRG_edopt.thrfilter_opt == FCFLAG_AVER )
00201 MRG_edopt.thrfilter_opt = FCFLAG_MEAN ;
00202
00203 continue ;
00204 }
00205
00206
00207
00208 if( strcmp(argv[nopt],"-1filter_winsor") == 0 ){
00209 int nwin ;
00210
00211 nopt++ ;
00212 if( nopt+1 >= argc ){
00213 fprintf(stderr,"*** Need 2 arguments after -1filter_winsor\n") ;
00214 exit(1) ;
00215 }
00216 MRG_edopt.filter_rmm = strtod( argv[nopt++] , NULL ) ;
00217 if( MRG_edopt.filter_rmm <= 0.0 ){
00218 fprintf(stderr,"*** Illegal rmm value after -1filter_winsor\n");
00219 exit(1) ;
00220 }
00221 nwin = (int) strtod( argv[nopt++] , NULL ) ;
00222 if( nwin <= 0 || nwin >= FCFLAG_ONE_STEP ){
00223 fprintf(stderr,"*** Illegal nw value after -1filter_winsor\n");
00224 exit(1) ;
00225 }
00226 MRG_edopt.filter_opt = FCFLAG_WINSOR + nwin ;
00227 MRG_have_edopt = 1 ; continue ;
00228 }
00229
00230
00231
00232 if( strncmp(argv[nopt],"-1filter_expr",13) == 0 ){
00233 PARSER_code * pcode ; int hsym[26] , aa , naa , nee ;
00234 double atoz[26] , val ;
00235
00236 nopt++ ;
00237 if( nopt+1 >= argc ){
00238 fprintf(stderr,"*** Need 2 arguments after -1filter_expr\n") ;
00239 exit(1) ;
00240 }
00241 MRG_edopt.filter_opt = FCFLAG_EXPR;
00242 MRG_edopt.filter_rmm = strtod( argv[nopt++] , NULL ) ;
00243 if( MRG_edopt.filter_rmm <= 0.0 ){
00244 fprintf(stderr,"*** Illegal rmm value after -1filter_expr\n");
00245 exit(1) ;
00246 }
00247
00248 MRG_edopt.fexpr = argv[nopt++] ;
00249 pcode = PARSER_generate_code( MRG_edopt.fexpr ) ;
00250
00251 if( pcode == NULL ){
00252 fprintf(stderr,"*** Illegal expr after -1filter_expr\n");
00253 exit(1);
00254 }
00255
00256 #undef MASK
00257 #undef PREDEFINED_MASK
00258 #define MASK(x) (1 << (x-'a'))
00259 #define PREDEFINED_MASK ( MASK('r') | MASK('x') | MASK('y') | MASK('z') \
00260 | MASK('i') | MASK('j') | MASK('k') )
00261
00262
00263
00264 PARSER_mark_symbols( pcode , hsym ) ;
00265
00266 for( nee=naa=aa=0 ; aa < 26 ; aa++ ){
00267 if( hsym[aa] ){
00268 if( ((1<<aa) & PREDEFINED_MASK) == 0 ){
00269 nee++ ;
00270 fprintf(stderr,"*** Symbol %c in -1filter_expr is illegal\n",'a'+aa) ;
00271 } else {
00272 naa++ ;
00273 }
00274 }
00275 }
00276 if( nee > 0 ) exit(1) ;
00277 if( naa == 0 ){
00278 double atoz[26] , val ;
00279 val = PARSER_evaluate_one( pcode , atoz ) ;
00280 if( val != 0.0 ){
00281 fprintf(stderr,"+++ Warning: -1filter_expr is constant = %g\n",val) ;
00282 } else {
00283 fprintf(stderr,"*** -1filter_expr is identically zero\n") ;
00284 exit(1) ;
00285 }
00286 }
00287
00288 free(pcode) ;
00289
00290 MRG_have_edopt = 1 ; continue ;
00291 }
00292
00293
00294
00295 if( strncmp(argv[nopt],"-quiet",6) == 0 ){
00296 MRG_be_quiet = 1 ;
00297 nopt++ ; continue ;
00298 }
00299
00300
00301
00302
00303 if( strncmp(argv[nopt],"-keepthr",6) == 0 ){
00304 MRG_keepthr = 1 ;
00305 nopt++ ; continue ;
00306 }
00307
00308
00309
00310 if( strncmp(argv[nopt],"-doall",6) == 0 ){
00311 MRG_doall = 1 ;
00312 nopt++ ; continue ;
00313 }
00314
00315
00316
00317 if( strncmp(argv[nopt],"-datum",6) == 0 ){
00318 DUMP2 ;
00319 if( ++nopt >= argc ){
00320 fprintf(stderr,"need an argument after -datum!\n") ; exit(1) ;
00321 }
00322
00323 if( strcmp(argv[nopt],"short") == 0 ){
00324 MRG_datum = MRI_short ;
00325 } else if( strcmp(argv[nopt],"float") == 0 ){
00326 MRG_datum = MRI_float ;
00327 } else if( strcmp(argv[nopt],"byte") == 0 ){
00328 MRG_datum = MRI_byte ;
00329 } else if( strcmp(argv[nopt],"complex") == 0 ){
00330 MRG_datum = MRI_complex ;
00331 } else {
00332 fprintf(stderr,"-datum of type '%s' is not supported in 3dmerge!\n",
00333 argv[nopt] ) ;
00334 exit(1) ;
00335 }
00336 nopt++ ; continue ;
00337 }
00338
00339
00340
00341 if( strncmp(argv[nopt],"-thdatum",6) == 0 ){
00342 DUMP2 ;
00343 if( ++nopt >= argc ){
00344 fprintf(stderr,"need an argument after -thdatum!\n") ; exit(1) ;
00345 }
00346
00347 if( strcmp(argv[nopt],"short") == 0 ){
00348 MRG_thdatum = MRI_short ;
00349 } else if( strcmp(argv[nopt],"float") == 0 ){
00350 MRG_thdatum = MRI_float ;
00351 } else if( strcmp(argv[nopt],"byte") == 0 ){
00352 MRG_thdatum = MRI_byte ;
00353 } else {
00354 fprintf(stderr,"-thdatum of type '%s' is not supported in 3dmerge!\n",
00355 argv[nopt] ) ;
00356 exit(1) ;
00357 }
00358 MRG_keepthr = 1 ;
00359 nopt++ ; continue ;
00360 }
00361
00362
00363
00364 if( strncmp(argv[nopt],"-ghits",6) == 0 ){
00365 DUMP2 ;
00366 nopt++ ;
00367 if( nopt >= argc ){
00368 fprintf(stderr,"need argument after -ghits!\n") ; exit(1) ;
00369 }
00370 MRG_hits_g = strtod( argv[nopt++] , NULL ) ;
00371 if( MRG_hits_g <= 0 ){
00372 fprintf(stderr,"illegal value after -ghits\n") ;
00373 exit(1) ;
00374 }
00375 continue ;
00376 }
00377
00378
00379
00380 if( strncmp(argv[nopt],"-gclust",6) == 0 ){
00381 DUMP3 ;
00382 nopt++ ;
00383 if( nopt+1 >= argc ){
00384 fprintf(stderr,"need 2 arguments after -gclust!\n") ;
00385 exit(1) ;
00386 }
00387 MRG_clust_rmm_g = strtod( argv[nopt++] , NULL ) ;
00388 MRG_clust_vmul_g = strtod( argv[nopt++] , NULL ) ;
00389 if( MRG_clust_rmm_g <= 0 || MRG_clust_vmul_g <= 0 ){
00390 fprintf(stderr,"illegal value after -gclust\n") ;
00391 exit(1) ;
00392 }
00393 continue ;
00394 }
00395
00396
00397
00398 if( strncmp(argv[nopt],"-session",6) == 0 ){
00399 DUMP2 ;
00400 nopt++ ;
00401 if( nopt >= argc ){
00402 fprintf(stderr,"need argument after -session!\n") ;
00403 exit(1) ;
00404 }
00405 MCW_strncpy( MRG_output_session , argv[nopt++] , THD_MAX_NAME ) ;
00406 continue ;
00407 }
00408
00409
00410
00411 if( strncmp(argv[nopt],"-prefix",6) == 0 ){
00412 DUMP2 ;
00413 nopt++ ;
00414 if( nopt >= argc ){
00415 fprintf(stderr,"need argument after -prefix!\n") ;
00416 exit(1) ;
00417 }
00418 MCW_strncpy( MRG_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00419 continue ;
00420 }
00421
00422 #if 0
00423
00424
00425 if( strncmp(argv[nopt],"-label",6) == 0 ){
00426 DUMP2 ;
00427 nopt++ ;
00428 if( nopt >= argc ){
00429 fprintf(stderr,"need argument after -label!\n") ;
00430 exit(1) ;
00431 }
00432 MCW_strncpy( MRG_output_label , argv[nopt++] , THD_MAX_LABEL ) ;
00433 continue ;
00434 }
00435 #endif
00436
00437
00438
00439 if( strcmp(argv[nopt],"-nscale") == 0 ){
00440 DUMP1 ;
00441 MRG_nscale = 1 ;
00442 nopt++ ; continue ;
00443 }
00444
00445
00446
00447
00448 if( strncmp(argv[nopt],"-gmean",6) == 0 ){
00449 DUMP1 ;
00450 MRG_cflag_g = CFLAG_MEAN ;
00451 nopt++ ; continue ;
00452 }
00453
00454
00455
00456 if( strncmp(argv[nopt],"-gfisher",6) == 0 ){
00457 DUMP1 ;
00458 MRG_cflag_g = CFLAG_FISHER ;
00459 nopt++ ; continue ;
00460 }
00461
00462
00463
00464 if( strncmp(argv[nopt],"-tgfisher",6) == 0 ){
00465 DUMP1 ;
00466 MRG_cflag_gthr = THFLAG_FICO ;
00467 nopt++ ; continue ;
00468 }
00469
00470
00471
00472 if( strncmp(argv[nopt],"-gnzmean",6) == 0 ){
00473 DUMP1 ;
00474 MRG_cflag_g = CFLAG_NZMEAN ;
00475 nopt++ ; continue ;
00476 }
00477
00478
00479
00480 if( strncmp(argv[nopt],"-gmax",6) == 0 ){
00481 DUMP1 ;
00482 MRG_cflag_g = CFLAG_MMAX ;
00483 nopt++ ; continue ;
00484 }
00485
00486
00487
00488 if( strncmp(argv[nopt],"-gamax",6) == 0 ){
00489 DUMP1 ;
00490 MRG_cflag_g = CFLAG_AMAX ;
00491 nopt++ ; continue ;
00492 }
00493
00494
00495
00496 if( strncmp(argv[nopt],"-gsmax",6) == 0 ){
00497 DUMP1 ;
00498 MRG_cflag_g = CFLAG_SMAX ;
00499 nopt++ ; continue ;
00500 }
00501
00502
00503
00504 if( strncmp(argv[nopt],"-gcount",6) == 0 ){
00505 DUMP1 ;
00506 MRG_cflag_g = CFLAG_COUNT ;
00507 nopt++ ; continue ;
00508 }
00509
00510
00511
00512 if( strncmp(argv[nopt],"-gorder",6) == 0 ){
00513 DUMP1 ;
00514 MRG_cflag_g = CFLAG_ORDER ;
00515 nopt++ ; continue ;
00516 }
00517
00518
00519
00520 fprintf(stderr,"unrecognized option %s\n",argv[nopt]) ;
00521 exit(1) ;
00522
00523 }
00524
00525
00526
00527 #if 0
00528 if( strlen(MRG_output_label) == 0 ){
00529 MCW_strncpy( MRG_output_label , MRG_output_prefix , THD_MAX_LABEL ) ;
00530 }
00531 #endif
00532
00533 return( nopt );
00534 }
00535
00536
00537
00538 void MRG_Syntax(void)
00539 {
00540 printf(
00541 "Edit and/or merge 3D datasets\n"
00542 "Usage: 3dmerge [options] datasets ...\n"
00543 "where the options are:\n"
00544 ) ;
00545
00546 printf( "%s\n" , EDIT_options_help() ) ;
00547
00548 printf(
00549 "OTHER OPTIONS:\n"
00550 " -datum type = Coerce the output data to be stored as the given type,\n"
00551 " which may be byte, short, or float.\n"
00552 " N.B.: Byte data cannot be negative. If this datum type is chosen,\n"
00553 " any negative values in the edited and/or merged dataset\n"
00554 " will be set to zero.\n"
00555 " -keepthr = When using 3dmerge to edit exactly one dataset of a\n"
00556 " functional type with a threshold statistic attached,\n"
00557 " normally the resulting dataset is of the 'fim'\n"
00558 " (intensity only) type. This option tells 3dmerge to\n"
00559 " copy the threshold data (unedited in any way) into\n"
00560 " the output dataset.\n"
00561 " N.B.: This option is ignored if 3dmerge is being used to\n"
00562 " combine 2 or more datasets.\n"
00563 " N.B.: The -datum option has no effect on the storage of the\n"
00564 " threshold data. Instead use '-thdatum type'.\n"
00565 "\n"
00566 " -doall = Apply editing and merging options to ALL sub-bricks \n"
00567 " uniformly in a dataset.\n"
00568 " N.B.: All input datasets must have the same number of sub-bricks\n"
00569 " when using the -doall option. \n"
00570 " N.B.: The threshold specific options (such as -1thresh, \n"
00571 " -keepthr, -tgfisher, etc.) are not compatible with \n"
00572 " the -doall command. Neither are the -1dindex or\n"
00573 " the -1tindex options.\n"
00574 " N.B.: All labels and statistical parameters for individual \n"
00575 " sub-bricks are copied from the first dataset. It is \n"
00576 " the responsibility of the user to verify that these \n"
00577 " are appropriate. Note that sub-brick auxiliary data \n"
00578 " can be modified using program 3drefit. \n"
00579 "\n"
00580 " -1dindex j = Uses sub-brick #j as the data source , and uses sub-brick\n"
00581 " -1tindex k = #k as the threshold source. With these, you can operate\n"
00582 " on any given sub-brick of the inputs dataset(s) to produce\n"
00583 " as output a 1 brick dataset. If desired, a collection\n"
00584 " of 1 brick datasets can later be assembled into a\n"
00585 " multi-brick bucket dataset using program '3dbucket'\n"
00586 " or into a 3D+time dataset using program '3dTcat'.\n"
00587 " N.B.: If these options aren't used, j=0 and k=1 are the defaults\n"
00588 "\n"
00589 " The following option allows you to specify a mask dataset that\n"
00590 " limits the action of the 'filter' options to voxels that are\n"
00591 " nonzero in the mask:\n"
00592 "\n"
00593 " -1fmask mset = Read dataset 'mset' (which can include a\n"
00594 " sub-brick specifier) and use the nonzero\n"
00595 " voxels as a mask for the filter options.\n"
00596 " Filtering calculations will not use voxels\n"
00597 " that are outside the mask. If an output\n"
00598 " voxel does not have ANY masked voxels inside\n"
00599 " the rmm radius, then that output voxel will\n"
00600 " be set to 0.\n"
00601 " N.B.: * Only the -1filter_* and -t1filter_* options are\n"
00602 " affected by -1fmask.\n"
00603 " * In the linear averaging filters (_mean, _nzmean,\n"
00604 " and _expr), voxels not in the mask will not be used\n"
00605 " or counted in either the numerator or denominator.\n"
00606 " This can give unexpected results. If the mask is\n"
00607 " designed to exclude the volume outside the brain,\n"
00608 " then voxels exterior to the brain, but within 'rmm',\n"
00609 " will have a few voxels inside the brain included\n"
00610 " in the filtering. Since the sum of weights (the\n"
00611 " denominator) is only over those few intra-brain\n"
00612 " voxels, the effect will be to extend the significant\n"
00613 " part of the result outward by rmm from the surface\n"
00614 " of the brain. In contrast, without the mask, the\n"
00615 " many small-valued voxels outside the brain would\n"
00616 " be included in the numerator and denominator sums,\n"
00617 " which would barely change the numerator (since the\n"
00618 " voxel values are small outside the brain), but would\n"
00619 " increase the denominator greatly (by including many\n"
00620 " more weights). The effect in this case (no -1fmask)\n"
00621 " is to make the filtering taper off gradually in the\n"
00622 " rmm-thickness shell around the brain.\n"
00623 " * Thus, if the -1fmask is intended to clip off non-brain\n"
00624 " data from the filtering, its use should be followed by\n"
00625 " masking operation using 3dcalc:\n"
00626 " 3dmerge -1filter_aver 12 -1fmask mask+orig -prefix x input+orig\n"
00627 " 3dcalc -a x -b mask+orig -prefix y -expr 'a*step(b)'\n"
00628 " rm -f x+orig.*\n"
00629 " The desired result is y+orig - filtered using only\n"
00630 " brain voxels (as defined by mask+orig), and with\n"
00631 " the output confined to the brain voxels as well.\n"
00632 "\n"
00633 " The following option allows you to specify an almost arbitrary\n"
00634 " weighting function for 3D linear filtering:\n"
00635 "\n"
00636 " -1filter_expr rmm expr\n"
00637 " Defines a linear filter about each voxel of radius 'rmm' mm.\n"
00638 " The filter weights are proportional to the expression evaluated\n"
00639 " at each voxel offset in the rmm neighborhood. You can use only\n"
00640 " these symbols in the expression:\n"
00641 " r = radius from center\n"
00642 " x = dataset x-axis offset from center\n"
00643 " y = dataset y-axis offset from center\n"
00644 " z = dataset z-axis offset from center\n"
00645 " i = x-axis index offset from center\n"
00646 " j = y-axis index offset from center\n"
00647 " k = z-axis index offset from center\n"
00648 " Example:\n"
00649 " -1filter_expr 12.0 'exp(-r*r/36.067)'\n"
00650 " This does a Gaussian filter over a radius of 12 mm. In this\n"
00651 " example, the FWHM of the filter is 10 mm. [in general, the\n"
00652 " denominator in the exponent would be 0.36067 * FWHM * FWHM.\n"
00653 " This is the only way to get a Gaussian blur combined with the\n"
00654 " -1fmask option. The radius rmm=12 is chosen where the weights\n"
00655 " get smallish.] Another example:\n"
00656 " -1filter_expr 20.0 'exp(-(x*x+16*y*y+z*z)/36.067)'\n"
00657 " which is a non-spherical Gaussian filter.\n"
00658 "\n"
00659 " The following option lets you apply a 'Winsor' filter to the data:\n"
00660 "\n"
00661 " -1filter_winsor rmm nw\n"
00662 " The data values within the radius rmm of each voxel are sorted.\n"
00663 " Suppose there are 'N' voxels in this group. We index the\n"
00664 " sorted voxels as s[0] <= s[1] <= ... <= s[N-1], and we call the\n"
00665 " value of the central voxel 'v' (which is also in array s[]).\n"
00666 " If v < s[nw] , then v is replaced by s[nw]\n"
00667 " otherwise If v > s[N-1-nw], then v is replace by s[N-1-nw]\n"
00668 " otherwise v is unchanged\n"
00669 " The effect is to increase 'too small' values up to some\n"
00670 " middling range, and to decrease 'too large' values.\n"
00671 " If N is odd, and nw=(N-1)/2, this would be a median filter.\n"
00672 " In practice, I recommend that nw be about N/4; for example,\n"
00673 " -dxyz=1 -1filter_winsor 2.5 19\n"
00674 " is a filter with N=81 that gives nice results.\n"
00675 " N.B.: This option is NOT affected by -1fmask\n"
00676 " N.B.: This option is slow!\n"
00677 "\n"
00678 "MERGING OPTIONS APPLIED TO FORM THE OUTPUT DATASET:\n"
00679 " [That is, different ways to combine results. The]\n"
00680 " [following '-g' options are mutually exclusive! ]\n"
00681 " -gmean = Combine datasets by averaging intensities\n"
00682 " (including zeros) -- this is the default\n"
00683 " -gnzmean = Combine datasets by averaging intensities\n"
00684 " (not counting zeros)\n"
00685 " -gmax = Combine datasets by taking max intensity\n"
00686 " (e.g., -7 and 2 combine to 2)\n"
00687 " -gamax = Combine datasets by taking max absolute intensity\n"
00688 " (e.g., -7 and 2 combine to 7)\n"
00689 " -gsmax = Combine datasets by taking max signed intensity\n"
00690 " (e.g., -7 and 2 combine to -7)\n"
00691 " -gcount = Combine datasets by counting number of 'hits' in\n"
00692 " each voxel (see below for defintion of 'hit')\n"
00693 " -gorder = Combine datasets in order of input:\n"
00694 " * If a voxel is nonzero in dataset #1, then\n"
00695 " that value goes into the voxel.\n"
00696 " * If a voxel is zero in dataset #1 but nonzero\n"
00697 " in dataset #2, then the value from #2 is used.\n"
00698 " * And so forth: the first dataset with a nonzero\n"
00699 " entry in a given voxel 'wins'\n"
00700 " -gfisher = Takes the arctanh of each input, averages these,\n"
00701 " and outputs the tanh of the average. If the input\n"
00702 " datum is 'short', then input values are scaled by\n"
00703 " 0.0001 and output values by 10000. This option\n"
00704 " is for merging bricks of correlation coefficients.\n"
00705 "\n"
00706 " -nscale = If the output datum is shorts, don't do the scaling\n"
00707 " to the max range [similar to 3dcalc's -nscale option]\n"
00708 "\n"
00709 "MERGING OPERATIONS APPLIED TO THE THRESHOLD DATA:\n"
00710 " [That is, different ways to combine the thresholds. If none of these ]\n"
00711 " [are given, the thresholds will not be merged and the output dataset ]\n"
00712 " [will not have threshold data attached. Note that the following '-tg']\n"
00713 " [command line options are mutually exclusive, but are independent of ]\n"
00714 " [the '-g' options given above for merging the intensity data values. ]\n"
00715 " -tgfisher = This option is only applicable if each input dataset\n"
00716 " is of the 'fico' or 'fith' types -- functional\n"
00717 " intensity plus correlation or plus threshold.\n"
00718 " (In the latter case, the threshold values are\n"
00719 " interpreted as correlation coefficients.)\n"
00720 " The correlation coefficients are averaged as\n"
00721 " described by -gfisher above, and the output\n"
00722 " dataset will be of the fico type if all inputs\n"
00723 " are fico type; otherwise, the output datasets\n"
00724 " will be of the fith type.\n"
00725 " N.B.: The difference between the -tgfisher and -gfisher\n"
00726 " methods is that -tgfisher applies to the threshold\n"
00727 " data stored with a dataset, while -gfisher\n"
00728 " applies to the intensity data. Thus, -gfisher\n"
00729 " would normally be applied to a dataset created\n"
00730 " from correlation coefficients directly, or from\n"
00731 " the application of the -1thtoin option to a fico\n"
00732 " or fith dataset.\n"
00733 "\n"
00734 "OPTIONAL WAYS TO POSTPROCESS THE COMBINED RESULTS:\n"
00735 " [May be combined with the above methods.]\n"
00736 " [Any combination of these options may be used.]\n"
00737 " -ghits count = Delete voxels that aren't !=0 in at least\n"
00738 " count datasets (!=0 is a 'hit')\n"
00739 " -gclust rmm vmul = Form clusters with connection distance rmm\n"
00740 " and clip off data not in clusters of\n"
00741 " volume at least vmul microliters\n"
00742 "\n"
00743 "The '-g' and '-tg' options apply to the entire group of input datasets.\n"
00744 "\n"
00745
00746 "OPTIONS THAT CONTROL THE NAMES OF THE OUTPUT DATASET:\n"
00747 " -session dirname = write output into given directory (default=./)\n"
00748 " -prefix pname = use 'pname' for the output directory prefix\n"
00749 " (default=mrg)\n"
00750 #if 0
00751 " -label string = use 'string' for the label in the output\n"
00752 " dataset (the label is used for switching\n"
00753 " between datasets in AFNI)\n"
00754 #endif
00755 "\n"
00756
00757 "NOTES:\n"
00758 " ** If only one dataset is read into this program, then the '-g'\n"
00759 " options do not apply, and the output dataset is simply the\n"
00760 " '-1' options applied to the input dataset (i.e., edited).\n"
00761 " ** A merged output dataset is ALWAYS of the intensity-only variety.\n"
00762 " ** You can combine the outputs of 3dmerge with other sub-bricks\n"
00763 " using the program 3dbucket.\n"
00764 " ** Complex-valued datasets cannot be merged.\n"
00765 " ** This program cannot handle time-dependent datasets without -doall.\n"
00766 " ** Note that the input datasets are specified by their .HEAD files,\n"
00767 " but that their .BRIK files must exist also!\n"
00768
00769 #ifdef ALLOW_SUBV
00770 "\n"
00771 MASTER_SHORTHELP_STRING
00772 "\n"
00773 " ** Input datasets using sub-brick selectors are treated as follows:\n"
00774 " - 3D+time if the dataset is 3D+time and more than 1 brick is chosen\n"
00775 " - otherwise, as bucket datasets (-abuc or -fbuc)\n"
00776 " (in particular, fico, fitt, etc. datasets are converted to fbuc)\n"
00777 " ** If you are NOT using -doall, and choose more than one sub-brick\n"
00778 " with the selector, then you may need to use -1dindex to further\n"
00779 " pick out the sub-brick on which to operate (why you would do this\n"
00780 " I cannot fathom). If you are also using a thresholding operation\n"
00781 " (e.g., -1thresh), then you also MUST use -1tindex to choose which\n"
00782 " sub-brick counts as the 'threshold' value. When used with sub-brick\n"
00783 " selection, 'index' refers the dataset AFTER it has been read in:\n"
00784 " -1dindex 1 -1tindex 3 'dset+orig[4..7]'\n"
00785 " means to use the #5 sub-brick of dset+orig as the data for merging\n"
00786 " and the #7 sub-brick of dset+orig as the threshold values.\n"
00787 " ** The above example would better be done with\n"
00788 " -1tindex 1 'dset+orig[5,7]'\n"
00789 " since the default data index is 0. (You would only use -1tindex if\n"
00790 " you are actually using a thresholding operation.)\n"
00791 " ** -1dindex and -1tindex apply to all input datasets.\n"
00792 #endif
00793 ) ;
00794 exit(0) ;
00795 }
00796
00797
00798
00799 int main( int argc , char * argv[] )
00800 {
00801 int file_num , first_file , nx,ny,nz , nxyz , ii , num_dset ,
00802 file_count , ptmin , iclu,nclu , edit_type , ival,ivout , tval ;
00803 float dx,dy,dz , fac , dxyz , rmm,vmul ;
00804 THD_3dim_dataset * dset=NULL , * new_dset=NULL ;
00805 THD_3dim_dataset ** dsetar=NULL ;
00806 short * gnum=NULL ;
00807 float * gfim=NULL , * tfim=NULL , * ggfim=NULL , * ttfim=NULL ;
00808 int datum ;
00809 MCW_cluster_array * clar ;
00810 float fimfac , fimfacinv , first_fimfac , thrfac ;
00811 int output_datum , output_thdatum ;
00812 int input_datum , input_thdatum , first_datum ;
00813
00814 float thr_stataux[MAX_STAT_AUX] ;
00815 int num_fico ;
00816 int is_int=1 ;
00817
00818 int iv, iv_bot, iv_top;
00819
00820
00821 printf ("\n\nProgram %s \n", PROGRAM_NAME);
00822 printf ("Last revision: %s \n\n", LAST_MOD_DATE);
00823
00824
00825
00826 if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) MRG_Syntax() ;
00827
00828
00829
00830 mainENTRY("3dmerge main") ; machdep() ; PRINT_VERSION("3dmerge") ;
00831
00832 { int new_argc ; char ** new_argv ;
00833 addto_args( argc , argv , &new_argc , &new_argv ) ;
00834 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00835 }
00836
00837 AFNI_logger("3dmerge",argc,argv) ;
00838
00839 first_file = MRG_read_opts( argc , argv ) ;
00840 file_count = argc - first_file ;
00841
00842 if( ! MRG_be_quiet )
00843 printf("3dmerge: edit and combine 3D datasets, by RW Cox\n") ;
00844
00845 if( first_file < 1 || first_file >= argc ){
00846 fprintf(stderr,"*** ILLEGAL COMMAND LINE ***\n") ; exit(1) ;
00847 }
00848
00849
00850 if (MRG_doall)
00851 { int nerr = 0 ;
00852
00853 if( MRG_ivfim >= 0 ){
00854 fprintf(stderr,"-1dindex is not compatible with -doall option \n");
00855 nerr++ ; }
00856 if( MRG_ivthr >= 0 ){
00857 fprintf(stderr,"-1dindex is not compatible with -doall option \n");
00858 nerr++ ; }
00859
00860 if (MRG_edopt.thtoin > 0) {
00861 fprintf (stderr, "-1thtoin is not compatible with -doall option \n");
00862 nerr++ ; }
00863 if (MRG_edopt.thresh > 0.0) {
00864 fprintf (stderr, "-1thresh is not compatible with -doall option \n");
00865 nerr++ ; }
00866 if (MRG_edopt.thrfilter_opt > 0) {
00867 fprintf (stderr, "-t1filter is not compatible with -doall option \n");
00868 nerr++ ; }
00869 if (MRG_edopt.thrblur > 0.0) {
00870 fprintf (stderr, "-t1blur is not compatible with -doall option \n");
00871 nerr++ ; }
00872 if (MRG_thdatum >= 0) {
00873 fprintf (stderr, "-thdatum is not compatible with -doall option \n");
00874 nerr++ ; }
00875 if (MRG_keepthr) {
00876 fprintf (stderr, "-keepthr is not compatible with -doall option \n");
00877 nerr++ ; }
00878 if (MRG_cflag_gthr > 0) {
00879 fprintf (stderr, "-tgfisher is not compatible with -doall option \n");
00880 nerr++ ; }
00881
00882 if( nerr > 0 ) exit(1) ;
00883 }
00884
00885
00886
00887 if( (MRG_ivfim >= 0 || MRG_ivthr >=0) && MRG_keepthr ){
00888 fprintf(stderr,"-keepthr is not compatible with -1dindex or -1tindex\n") ;
00889 exit(1) ;
00890 }
00891
00892 if( (MRG_ivfim >= 0 || MRG_ivthr >=0) && MRG_cflag_gthr ){
00893 fprintf(stderr,"-tgfisher is not compatible with -1dindex or -1tindex\n") ;
00894 exit(1) ;
00895 }
00896
00897
00898 for (file_num = first_file; file_num < argc; file_num++)
00899 {
00900 dset = DSET_OPEN( argv[file_num] ) ;
00901 if( ! ISVALID_3DIM_DATASET(dset) )
00902 {
00903 fprintf(stderr,"*** cannot open dataset %s\n",argv[file_num]) ;
00904 exit(1) ;
00905 }
00906
00907
00908
00909 if( MRG_ivfim >= DSET_NVALS(dset) ){
00910 fprintf(stderr,
00911 "*** Dataset %s does not have enough bricks for -1dindex %d\n" ,
00912 argv[file_num],MRG_ivfim) ;
00913 exit(1) ;
00914 }
00915 if( MRG_ivthr >= DSET_NVALS(dset) ){
00916 fprintf(stderr,
00917 "*** Dataset %s does not have enough bricks for -1tindex %d\n" ,
00918 argv[file_num],MRG_ivthr) ;
00919 exit(1) ;
00920 }
00921
00922 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00923 }
00924
00925
00926
00927 if( MRG_ivfim < 0 ) fprintf(stderr,"++ default -1dindex = 0\n") ;
00928 if( MRG_ivthr < 0 ) fprintf(stderr,"++ default -1tindex = 1\n") ;
00929
00930
00931
00932 dset = DSET_OPEN( argv[first_file] ) ;
00933 if( ! ISVALID_3DIM_DATASET(dset) ){
00934 fprintf(stderr,"*** Unable to open first dataset %s\n",argv[first_file]) ;
00935 exit(1) ;
00936 }
00937
00938 if( DSET_NUM_TIMES(dset) > 1 && (!MRG_doall && MRG_ivfim < 0) ){
00939 fprintf(stderr, "*** Unable to merge time-dependent datasets"
00940 " without -doall or -1dindex\n") ;
00941 exit(1) ;
00942 }
00943
00944
00945
00946 nx = dset->daxes->nxx ;
00947 ny = dset->daxes->nyy ;
00948 nz = dset->daxes->nzz ; nxyz = nx*ny*nz ;
00949
00950 dx = fabs(dset->daxes->xxdel) ;
00951 dy = fabs(dset->daxes->yydel) ;
00952 dz = fabs(dset->daxes->zzdel) ;
00953
00954 if( MRG_edopt.fake_dxyz ) dx = dy = dz = 1.0 ;
00955
00956 nice(1) ;
00957
00958 if( MRG_edopt.nfmask > 0 && MRG_edopt.nfmask != nxyz ){
00959 fprintf(stderr,
00960 "*** -1fmask and 1st dataset don't have same number of voxels\n\a");
00961 exit(1) ;
00962 }
00963
00964
00965
00966
00967
00968 if( file_count == 1 ){
00969
00970 ival = DSET_PRINCIPAL_VALUE(dset) ;
00971 input_datum = DSET_BRICK_TYPE(dset,ival) ;
00972 if( MRG_datum >= 0 ) output_datum = MRG_datum ;
00973 else output_datum = input_datum ;
00974
00975
00976
00977
00978
00979
00980
00981
00982 new_dset = EDIT_empty_copy( dset ) ;
00983
00984 EDIT_dset_items( new_dset ,
00985 ADN_prefix , MRG_output_prefix ,
00986 ADN_directory_name , MRG_output_session ,
00987 ADN_none ) ;
00988
00989 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
00990 fprintf(stderr,
00991 "*** Output file %s already exists -- cannot continue!\n",
00992 new_dset->dblk->diskptr->header_name ) ;
00993 exit(1) ;
00994 }
00995
00996 THD_delete_3dim_dataset( new_dset , False ) ;
00997
00998
00999
01000 if( MRG_keepthr && !ISFUNC(dset) ){
01001 fprintf(stderr,"*** -keepthr can't be used on non-functional dataset!\n");
01002 exit(1) ;
01003 }
01004 if( MRG_keepthr && dset->func_type == FUNC_FIM_TYPE ){
01005 fprintf(stderr,"*** -keepthr can't be used on 'fim' type dataset!\n") ;
01006 exit(1) ;
01007 }
01008 if( MRG_keepthr && dset->func_type == FUNC_BUCK_TYPE ){
01009 fprintf(stderr,"*** -keepthr can't be used on 'fbuc' type dataset!\n"
01010 " You can use '3dbuc2fim' first, if needed.\n" );
01011 exit(1) ;
01012 }
01013
01014
01015
01016 if( ! MRG_be_quiet ){
01017 printf("-- editing input dataset in memory (%.1f MB)\n",
01018 ((double)dset->dblk->total_bytes) / MEGA ) ;
01019 fflush(stdout) ;
01020 }
01021
01022
01023 if (MRG_doall)
01024 { iv_bot = 0; iv_top = DSET_NVALS(dset); }
01025 else if( MRG_ivfim >= 0 )
01026 { iv_bot = MRG_ivfim ; iv_top = iv_bot+1 ; }
01027 else
01028 { iv_bot = DSET_PRINCIPAL_VALUE(dset); iv_top = iv_bot + 1; }
01029
01030
01031 for (iv = iv_bot; iv < iv_top; iv++)
01032 {
01033 if ((!MRG_be_quiet) && MRG_doall) printf ("Editing sub-brick %d\n", iv);
01034
01035 MRG_edopt.iv_fim = iv;
01036
01037 EDIT_one_dataset( dset , &MRG_edopt ) ;
01038
01039 if( !MRG_be_quiet && !MRG_doall ){ printf(".") ; fflush(stdout) ; }
01040 }
01041
01042 if( MRG_edopt.nfmask > 0 ){
01043 free(MRG_edopt.fmask) ; MRG_edopt.fmask = NULL ; MRG_edopt.nfmask = 0 ;
01044 }
01045
01046 if( !MRG_be_quiet && !MRG_doall ) printf("\n") ;
01047
01048
01049
01050 new_dset = EDIT_empty_copy( dset ) ;
01051
01052 tross_Copy_History( dset , new_dset ) ;
01053 tross_Make_History( "3dmerge" , argc , argv , new_dset ) ;
01054
01055 EDIT_dset_items( new_dset ,
01056 ADN_prefix , MRG_output_prefix ,
01057 ADN_label1 , MRG_output_prefix ,
01058 ADN_directory_name , MRG_output_session ,
01059 ADN_none ) ;
01060 strcat( new_dset->self_name , "(ED)" ) ;
01061
01062 if( (! MRG_keepthr) && (new_dset->dblk->nvals > 1) && (! MRG_doall) )
01063 EDIT_dset_items( new_dset ,
01064 ADN_nvals , 1 ,
01065 ADN_ntt , 0 ,
01066 ADN_func_type , FUNC_FIM_TYPE ,
01067 ADN_none ) ;
01068
01069 if( MRG_keepthr && ISFUNC(new_dset) && FUNC_HAVE_THR(new_dset->func_type) ){
01070 ii = FUNC_ival_thr[dset->func_type] ;
01071 input_thdatum = DSET_BRICK_TYPE(dset,ii) ;
01072 if( MRG_thdatum >= 0 ) output_thdatum = MRG_thdatum ;
01073 else output_thdatum = input_thdatum ;
01074 } else {
01075 output_thdatum = input_thdatum = ILLEGAL_TYPE ;
01076 }
01077
01078
01079
01080
01081 for (iv = iv_bot; iv < iv_top; iv++)
01082 {
01083
01084 if (MRG_doall)
01085 {
01086 ival = iv;
01087 ivout = iv;
01088 }
01089 else if( MRG_ivfim >= 0 )
01090 {
01091 ival = MRG_ivfim ;
01092 ivout = DSET_PRINCIPAL_VALUE(new_dset) ;
01093 }
01094 else
01095 {
01096 ival = DSET_PRINCIPAL_VALUE(dset) ;
01097 ivout = DSET_PRINCIPAL_VALUE(new_dset) ;
01098 }
01099
01100 if( input_datum == output_datum ){
01101
01102
01103
01104
01105 mri_fix_data_pointer( DSET_ARRAY(dset,ival) , DSET_BRICK(new_dset,ivout) ) ;
01106
01107 #if 1
01108 # if 0
01109 if( ivout != ival )
01110 # endif
01111 DSET_BRICK_FACTOR(new_dset,ivout) = DSET_BRICK_FACTOR(dset,ival) ;
01112 #endif
01113
01114 } else {
01115
01116
01117
01118 void * dfim , * efim ;
01119 float efac = DSET_BRICK_FACTOR(dset,ival) ;
01120
01121 if( ! MRG_be_quiet ){
01122 printf("-- coercing output datum to be %s\n",
01123 MRI_TYPE_name[output_datum]);
01124 }
01125
01126 efim = DSET_ARRAY(dset,ival) ;
01127 dfim = (void *) malloc( mri_datum_size(output_datum) * nxyz ) ;
01128 if( dfim == NULL ){
01129 fprintf(stderr,"*** Can't malloc output brick #%d\n",ivout); exit(1);
01130 }
01131
01132
01133
01134 if( MRI_IS_INT_TYPE(output_datum) ){
01135 fimfac = EDIT_coerce_autoscale( nxyz , input_datum , efim ,
01136 output_datum , dfim ) ;
01137 if( fimfac == 0.0 ) fimfac = 1.0 ;
01138 if( efac != 0.0 ) fimfac /= efac ;
01139
01140 DSET_BRICK_FACTOR(new_dset,ivout) = (fimfac != 0.0 && fimfac != 1.0)
01141 ? 1.0/fimfac : 0.0 ;
01142 } else {
01143
01144 EDIT_coerce_scale_type( nxyz , efac , input_datum , efim ,
01145 output_datum , dfim ) ;
01146 DSET_BRICK_FACTOR(new_dset,ivout) = 0.0 ;
01147 }
01148
01149 mri_free( DSET_BRICK(dset,ival) ) ;
01150 EDIT_substitute_brick( new_dset , ivout , output_datum , dfim ) ;
01151 }
01152
01153
01154
01155 if( output_thdatum >= 0 ){
01156
01157 ival = FUNC_ival_thr[ dset->func_type] ;
01158 ivout = FUNC_ival_thr[new_dset->func_type] ;
01159
01160 if( input_thdatum == output_thdatum ){
01161
01162 mri_fix_data_pointer( DSET_ARRAY(dset,ival),DSET_BRICK(new_dset,ivout) ) ;
01163
01164 #if 0
01165 DSET_BRICK_FACTOR(new_dset,ivout) = DSET_BRICK_FACTOR(dset,ival) ;
01166 #endif
01167
01168 } else {
01169 void * dfim , * efim ;
01170
01171 if( ! MRG_be_quiet ){
01172 printf("-- coercing threshold datum to be %s\n",
01173 MRI_TYPE_name[output_thdatum]);
01174 }
01175
01176 efim = DSET_ARRAY(dset,ival) ;
01177 dfim = (void *) XtMalloc( mri_datum_size(output_thdatum) * nxyz ) ;
01178
01179 switch( output_thdatum ){
01180 default: fprintf(stderr,"** illegal output_thdatum = %d\n",
01181 output_thdatum);
01182 exit(1) ;
01183
01184 case MRI_float:
01185 fimfacinv = 0.0 ;
01186 fimfac = DSET_BRICK_FACTOR(dset,ival) ;
01187 if( fimfac == 0.0 ){
01188 fimfac = (input_thdatum == MRI_short)
01189 ? 1.0/FUNC_scale_short[dset->func_type]
01190 : (input_thdatum == MRI_byte)
01191 ? 1.0/FUNC_scale_byte[dset->func_type] : 0.0 ;
01192 }
01193 break ;
01194
01195 case MRI_short:
01196 if( input_datum == MRI_float ){
01197 fimfac = FUNC_scale_short[new_dset->func_type] ;
01198 fimfacinv = 1.0 / fimfac ;
01199 } else if( input_datum == MRI_byte ){
01200 fimfac = ((float)FUNC_scale_short[new_dset->func_type])
01201 / FUNC_scale_byte[new_dset->func_type] ;
01202 fimfacinv = 1.0 / FUNC_scale_short[new_dset->func_type] ;
01203 } else {
01204 fprintf(stderr,"** illegal input_thdatum = %d\n",input_thdatum);
01205 exit(1) ;
01206 }
01207 break ;
01208
01209 case MRI_byte:
01210 if( input_datum == MRI_float ){
01211 fimfac = FUNC_scale_byte[new_dset->func_type] ;
01212 fimfacinv = 1.0 / fimfac ;
01213 } else if( input_datum == MRI_short ){
01214 fimfac = ((float)FUNC_scale_byte[new_dset->func_type])
01215 / FUNC_scale_short[new_dset->func_type] ;
01216 fimfacinv = 1.0 / FUNC_scale_byte[new_dset->func_type] ;
01217 } else {
01218 fprintf(stderr,"** illegal input_thdatum = %d\n",input_thdatum);
01219 exit(1) ;
01220 }
01221 break ;
01222 }
01223
01224 EDIT_coerce_scale_type( nxyz , fimfac ,
01225 DSET_BRICK_TYPE(dset,ival),efim ,
01226 output_thdatum,dfim ) ;
01227
01228 DSET_BRICK_FACTOR(new_dset,ivout) = fimfacinv ;
01229 EDIT_substitute_brick( new_dset , ivout , output_thdatum , dfim ) ;
01230 mri_free( DSET_BRICK(dset,ival) ) ;
01231 }
01232 }
01233
01234 }
01235
01236
01237 THD_load_statistics( new_dset ) ;
01238 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01239 if( ! MRG_be_quiet )
01240 fprintf(stderr,"-- Wrote edited dataset: %s\n" , DSET_BRIKNAME(new_dset) ) ;
01241 exit(0) ;
01242 }
01243
01244
01245
01246
01247
01248
01249 dsetar = (THD_3dim_dataset **) malloc(sizeof(THD_3dim_dataset *)*file_count);
01250 for( ii=0 ; ii < file_count ; ii++ ) dsetar[ii] = NULL ;
01251 dsetar[0] = dset ;
01252
01253
01254
01255 new_dset = EDIT_empty_copy( dset ) ;
01256 ival = DSET_PRINCIPAL_VALUE(dset) ;
01257
01258 if( MRG_datum >= 0 ) output_datum = MRG_datum ;
01259 else output_datum = DSET_BRICK_TYPE(dset,ival) ;
01260
01261 EDIT_dset_items( new_dset ,
01262 ADN_prefix , MRG_output_prefix ,
01263 ADN_label1 , MRG_output_prefix ,
01264 ADN_directory_name , MRG_output_session ,
01265 ADN_none ) ;
01266 strcat( new_dset->self_name , "(MG)" ) ;
01267
01268
01269
01270 if (! MRG_doall)
01271 switch( MRG_cflag_gthr ){
01272 default:
01273 EDIT_dset_items( new_dset , ADN_nvals,1 , ADN_ntt,0 , ADN_none ) ;
01274 if( ISFUNC(dset) )
01275 EDIT_dset_items( new_dset , ADN_func_type,FUNC_FIM_TYPE , ADN_none ) ;
01276 break ;
01277
01278 case THFLAG_FICO:
01279 num_fico = 0 ;
01280 break ;
01281 }
01282
01283 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
01284 fprintf(stderr,
01285 "*** Output file %s already exists -- cannot continue!\n",
01286 new_dset->dblk->diskptr->header_name ) ;
01287 exit(1) ;
01288 }
01289
01290 if( ! MRG_be_quiet && MRG_keepthr )
01291 printf("-- ignoring -keepthr option\n") ;
01292
01293
01294
01295 if (MRG_doall)
01296 { iv_bot = 0; iv_top = DSET_NVALS(dset); }
01297 else if( MRG_ivfim >= 0 )
01298 { iv_bot = MRG_ivfim ; iv_top = iv_bot+1 ; }
01299 else
01300 { iv_bot = DSET_PRINCIPAL_VALUE(dset); iv_top = iv_bot + 1; }
01301
01302
01303
01304 ivout = 0 ;
01305
01306 for (iv = iv_bot; iv < iv_top; iv++)
01307 {
01308 if ((!MRG_be_quiet) ) printf ("-- Editing sub-brick %d \n", iv);
01309
01310 MRG_edopt.iv_fim = iv;
01311
01312
01313
01314 tfim = (float *) XtMalloc( sizeof(float) * nxyz ) ;
01315 gfim = (float *) XtMalloc( sizeof(float) * nxyz ) ;
01316 gnum = (short *) XtMalloc( sizeof(short) * nxyz ) ;
01317 for( ii=0 ; ii < nxyz ; ii++ ) gfim[ii] = 0.0 ;
01318 for( ii=0 ; ii < nxyz ; ii++ ) gnum[ii] = 0 ;
01319
01320
01321
01322 if( MRG_cflag_gthr != THFLAG_NONE ){
01323 ttfim = (float *) XtMalloc( sizeof(float) * nxyz ) ;
01324 ggfim = (float *) XtMalloc( sizeof(float) * nxyz ) ;
01325 for( ii=0 ; ii < nxyz ; ii++ ) ggfim[ii] = 0.0 ;
01326
01327 for( ii=0 ; ii < MAX_STAT_AUX ; ii++ ) thr_stataux[ii] = 0 ;
01328 }
01329
01330 if( ! MRG_be_quiet ){
01331 float nbytes = (2.0*sizeof(float)+sizeof(short))*nxyz ;
01332 if( MRG_cflag_gthr != THFLAG_NONE ) nbytes += 2.0*sizeof(float)*nxyz ;
01333 printf("-- allocated %.1f MB scratch memory\n", nbytes/MEGA ) ;
01334 }
01335
01336
01337
01338 num_dset = 0 ;
01339 for( file_num=first_file; file_num < argc ; file_num++ ){
01340
01341
01342
01343 dset = dsetar[ file_num - first_file ] ;
01344 if( dset == NULL ){
01345 dset = dsetar[ file_num - first_file ] = DSET_OPEN( argv[file_num] ) ;
01346 if( ! ISVALID_3DIM_DATASET(dset) ){
01347 fprintf(stderr,"*** cannot open dataset %s\n",argv[file_num]); exit(1);
01348 }
01349 }
01350
01351
01352
01353 if( dset->daxes->nxx != nx ||
01354 dset->daxes->nyy != ny || dset->daxes->nzz != nz ){
01355
01356 fprintf(stderr,"*** dataset brick size mismatch at file %s\n",
01357 argv[file_num] ) ;
01358 exit(1) ;
01359 }
01360
01361
01362 if ( (MRG_doall) && (DSET_NVALS(dset) != iv_top) )
01363 fprintf (stderr, "*** -doall dataset nvals mismatch at file %s\n",
01364 argv[file_num]);
01365
01366
01367 if( DSET_NUM_TIMES(dset) > 1 && (!MRG_doall && MRG_ivfim < 0) ){
01368 fprintf(stderr,
01369 "*** cannot use time-dependent dataset %s"
01370 " without -doall or -1dindex\n",argv[file_num]) ;
01371 exit(1) ;
01372 }
01373
01374
01375
01376 if( MRG_cflag_gthr == THFLAG_FICO ){
01377
01378 if( !ISFUNC(dset) ){
01379 fprintf(stderr,
01380 "*** dataset from file %s is anatomical using '-tgfisher'!\n",
01381 argv[file_num] ) ;
01382 exit(1) ;
01383 }
01384
01385 switch( dset->func_type ){
01386 default:
01387 fprintf(stderr,
01388 "*** dataset from file %s is illegal type using '-tgfisher'!\n",
01389 argv[file_num] ) ;
01390 exit(1) ;
01391
01392 case FUNC_COR_TYPE:
01393 num_fico ++ ;
01394 for( ii=0 ; ii < FUNC_need_stat_aux[FUNC_COR_TYPE] ; ii++ )
01395 thr_stataux[ii] += dset->stat_aux[ii] ;
01396 break ;
01397
01398 case FUNC_THR_TYPE:
01399 break ;
01400 }
01401 }
01402
01403
01404
01405 #if 1
01406 ival = iv ;
01407 #else
01408 if (MRG_doall) ival = iv;
01409 else ival = DSET_PRINCIPAL_VALUE(dset) ;
01410 #endif
01411 datum = DSET_BRICK_TYPE(dset,ival) ;
01412
01413 if( ! AFNI_GOOD_FUNC_DTYPE(datum) ){
01414 fprintf(stderr,"*** Illegal datum for 3dmerge: %s in file %s ***\n" ,
01415 MRI_TYPE_name[datum] , argv[file_num] ) ;
01416 exit(1) ;
01417 }
01418
01419 if( ! MRG_be_quiet ){
01420 printf("-- processing file %s" , argv[file_num] ) ;
01421 fflush(stdout) ;
01422 }
01423
01424
01425
01426 if( MRG_have_edopt )
01427 EDIT_one_dataset( dset , &MRG_edopt ) ;
01428 else
01429 THD_load_datablock( dset->dblk ) ;
01430
01431
01432
01433 if( !DSET_LOADED(dset) ){
01434 fprintf(stderr,"** Can't get data from %s\n",DSET_BRIKNAME(dset)) ;
01435 exit(1) ;
01436 }
01437
01438 if( ! MRG_be_quiet ){ printf(".") ; fflush(stdout) ; }
01439
01440
01441
01442 fimfac = DSET_BRICK_FACTOR(dset,ival) ;
01443
01444 if( MRG_cflag_g == CFLAG_FISHER &&
01445 DSET_BRICK_TYPE(dset,ival) == MRI_short &&
01446 (fimfac==0.0 || fimfac==1.0) ){
01447
01448 fimfac = 1.0 / FUNC_COR_SCALE_SHORT ;
01449 }
01450
01451 if( num_dset == 0 ){
01452 first_fimfac = fimfac ;
01453 first_datum = datum ;
01454 }
01455
01456
01457
01458 is_int = is_int && (MRI_IS_INT_TYPE(datum) && fimfac == 0.0) ;
01459
01460
01461
01462 EDIT_coerce_scale_type( nxyz , fimfac ,
01463 DSET_BRICK_TYPE(dset,ival) , DSET_ARRAY(dset,ival) ,
01464 MRI_float , tfim ) ;
01465
01466
01467
01468 if( MRG_cflag_gthr != THFLAG_NONE && (tval=DSET_THRESH_VALUE(dset)) >= 0 ){
01469
01470 int thdatum = DSET_BRICK_TYPE(dset,tval) ;
01471
01472 if( ! AFNI_GOOD_FUNC_DTYPE(thdatum) ){
01473 fprintf(stderr,"*** Illegal threshold for 3dmerge: %s in file %s ***\n" ,
01474 MRI_TYPE_name[thdatum] , argv[file_num] ) ;
01475 exit(1) ;
01476 }
01477
01478 thrfac = DSET_BRICK_FACTOR(dset,tval) ;
01479
01480 if( MRG_cflag_gthr == THFLAG_FICO &&
01481 DSET_BRICK_TYPE(dset,tval) == MRI_short &&
01482 (thrfac==0.0 || thrfac==1.0) ){
01483
01484 thrfac = 1.0 / FUNC_COR_SCALE_SHORT ;
01485 }
01486
01487 EDIT_coerce_scale_type( nxyz , thrfac ,
01488 thdatum , DSET_ARRAY(dset,tval) ,
01489 MRI_float , ttfim ) ;
01490 }
01491
01492 DSET_unload(dset) ;
01493
01494
01495
01496 if( MRG_cflag_g == CFLAG_MMAX ){
01497 for( ii=0 ; ii < nxyz ; ii++ ){
01498 if( tfim[ii] != 0 ){
01499 gnum[ii]++ ;
01500 if( tfim[ii] > gfim[ii] ) gfim[ii] = tfim[ii] ;
01501 }
01502 }
01503 } else if( MRG_cflag_g == CFLAG_AMAX ){
01504 float dab ;
01505 for( ii=0 ; ii < nxyz ; ii++ ){
01506 if( tfim[ii] != 0 ){
01507 gnum[ii]++ ;
01508 dab = fabs(tfim[ii]) ;
01509 if( dab > gfim[ii] ) gfim[ii] = dab ;
01510 }
01511 }
01512 } else if( MRG_cflag_g == CFLAG_SMAX ){
01513 float dab ;
01514 for( ii=0 ; ii < nxyz ; ii++ ){
01515 if( tfim[ii] != 0 ){
01516 gnum[ii]++ ;
01517 dab = fabs(tfim[ii]) ;
01518 if( dab > fabs(gfim[ii]) ) gfim[ii] = tfim[ii] ;
01519 }
01520 }
01521 } else if( MRG_cflag_g == CFLAG_ORDER ){
01522 for( ii=0 ; ii < nxyz ; ii++ ){
01523 if( tfim[ii] != 0 ){
01524 gnum[ii]++ ;
01525 if( gfim[ii] == 0 ) gfim[ii] = tfim[ii] ;
01526 }
01527 }
01528 } else if( MRG_cflag_g == CFLAG_FISHER ){
01529 for( ii=0 ; ii < nxyz ; ii++ ){
01530 if( tfim[ii] != 0 ){
01531 gnum[ii]++ ; gfim[ii] += ATANH(tfim[ii]) ;
01532 }
01533 }
01534 } else {
01535 for( ii=0 ; ii < nxyz ; ii++ ){
01536 if( tfim[ii] != 0 ){
01537 gnum[ii]++ ; gfim[ii] += tfim[ii] ;
01538 }
01539 }
01540 }
01541
01542
01543
01544 if( MRG_cflag_gthr == THFLAG_FICO ){
01545 for( ii=0 ; ii < nxyz ; ii++ ) ggfim[ii] += ATANH(ttfim[ii]) ;
01546 }
01547
01548 if( ! MRG_be_quiet ){ printf(".\n") ; fflush(stdout) ; }
01549
01550 num_dset++ ;
01551 }
01552
01553 myXtFree(tfim) ;
01554 if( ttfim != NULL ) myXtFree(ttfim) ;
01555
01556 if( MRG_edopt.nfmask > 0 ){
01557 free(MRG_edopt.fmask) ; MRG_edopt.fmask = NULL ; MRG_edopt.nfmask = 0 ;
01558 }
01559
01560
01561
01562 if( num_dset <= 1 ){
01563 fprintf(stderr,"*** Only found 1 dataset -- computations aborted!\n") ;
01564 exit(1) ;
01565 }
01566
01567 if( ! MRG_be_quiet ) printf("-- merging results into sub-brick %d\n",ivout) ;
01568
01569
01570
01571
01572 if( MRG_hits_g > 0 ){
01573 for( ii=0 ; ii < nxyz ; ii++ )
01574 if( gnum[ii] < MRG_hits_g ) { gfim[ii] = 0 ; gnum[ii] = 0 ; }
01575 }
01576
01577
01578
01579
01580
01581 switch( MRG_cflag_g ){
01582 default: is_int = 0 ; break ;
01583
01584 case CFLAG_COUNT: is_int = 1 ; break ;
01585
01586 case CFLAG_MMAX:
01587 case CFLAG_SMAX:
01588 case CFLAG_AMAX:
01589 case CFLAG_ORDER: break ;
01590 }
01591
01592
01593
01594 switch( MRG_cflag_g ){
01595 default: break ;
01596
01597 case CFLAG_COUNT:
01598 first_fimfac = 0.0 ;
01599 for( ii=0 ; ii < nxyz ; ii++ ) gfim[ii] = gnum[ii] ;
01600 break ;
01601
01602 case CFLAG_MEAN:
01603 fac = 1.0 / num_dset ;
01604 for( ii=0 ; ii < nxyz ; ii++ ) gfim[ii] *= fac ;
01605 break ;
01606
01607 case CFLAG_NZMEAN:
01608 for( ii=0 ; ii < nxyz ; ii++ )
01609 if( gnum[ii] > 0 ) gfim[ii] /= gnum[ii] ;
01610 break ;
01611
01612 case CFLAG_FISHER:
01613 fac = 1.0 / num_dset ;
01614 for( ii=0 ; ii < nxyz ; ii++ ) gfim[ii] = TANH( fac * gfim[ii] ) ;
01615 break ;
01616 }
01617
01618
01619
01620 switch( MRG_cflag_gthr ){
01621 default: break ;
01622
01623 case THFLAG_FICO:
01624 fac = 1.0 / num_dset ;
01625 for( ii=0 ; ii < nxyz ; ii++ ) ggfim[ii] = TANH( fac * ggfim[ii] ) ;
01626 break ;
01627 }
01628
01629
01630
01631
01632 myXtFree( gnum ) ;
01633
01634
01635
01636 rmm = MRG_clust_rmm_g ;
01637 vmul = MRG_clust_vmul_g ;
01638 dxyz = dx*dy*dz ;
01639 ptmin = vmul / dxyz + 0.99 ;
01640
01641 if( (rmm >= dx || rmm >= dy || rmm >= dz) && ptmin > 1 ){
01642 if( ! MRG_be_quiet ) printf("-- editing merger for cluster size\n") ;
01643
01644 clar = MCW_find_clusters( nx,ny,nz , dx,dy,dz , MRI_float,gfim , rmm ) ;
01645 nclu = 0 ;
01646 if( clar != NULL ){
01647 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
01648 if( clar->clar[iclu] != NULL && clar->clar[iclu]->num_pt < ptmin ){
01649 KILL_CLUSTER(clar->clar[iclu]) ;
01650 } else if( clar->clar[iclu] != NULL ){
01651 nclu++ ;
01652 }
01653 }
01654 }
01655
01656 if( nclu > 0 ){
01657 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
01658 if( clar->clar[iclu] != NULL && clar->clar[iclu]->num_pt > 0 )
01659 MCW_cluster_to_vol( nx,ny,nz , MRI_float,gfim , clar->clar[iclu] ) ;
01660 }
01661 }
01662 }
01663
01664
01665
01666 if( !MRG_doall ){
01667 for( ii=0 ; ii < nxyz ; ii++ ) if( gfim[ii] != 0 ) break ;
01668 if( ii == nxyz ){
01669 fprintf(stderr,
01670 "*** Merged dataset has no nonzero entries -- will not write\n" ) ;
01671 exit(0) ;
01672 }
01673 }
01674
01675
01676
01677
01678 switch( output_datum ){
01679
01680 default:
01681 fprintf(stderr,
01682 "*** Fatal Error ***\n"
01683 "*** Somehow ended up with output_datum = %d\n",output_datum) ;
01684 exit(1) ;
01685
01686 case MRI_complex:{
01687 void * dfim ;
01688 dfim = (void *) XtMalloc( sizeof(complex) * nxyz ) ;
01689 EDIT_coerce_type( nxyz , MRI_float,gfim , MRI_complex,dfim ) ;
01690 myXtFree( gfim ) ;
01691 EDIT_substitute_brick( new_dset , ivout , MRI_complex , dfim ) ;
01692 DSET_BRICK_FACTOR(new_dset,ivout) = 0.0 ;
01693 }
01694 break ;
01695
01696 case MRI_float:
01697 EDIT_substitute_brick( new_dset , ivout , MRI_float , gfim ) ;
01698 DSET_BRICK_FACTOR(new_dset,ivout) = 0.0 ;
01699 break ;
01700
01701 case MRI_byte:
01702 case MRI_short:{
01703 void * dfim ;
01704 float gtop ;
01705
01706 gtop = MCW_vol_amax( nx,ny,nz , MRI_float,gfim ) ;
01707
01708 if( MRG_cflag_g == CFLAG_FISHER ){
01709 fimfac = FUNC_COR_SCALE_SHORT ;
01710 } else if( gtop == 0.0 ||
01711 MRG_nscale ||
01712 (is_int && gtop <= MRI_TYPE_maxval[output_datum]) ){
01713 fimfac = 0.0 ;
01714 } else {
01715 fimfac = MRI_TYPE_maxval[output_datum] / gtop ;
01716 }
01717
01718 dfim = (void *) XtMalloc( mri_datum_size(output_datum) * nxyz ) ;
01719 EDIT_coerce_scale_type( nxyz,fimfac , MRI_float,gfim , output_datum,dfim ) ;
01720 myXtFree( gfim ) ;
01721 EDIT_substitute_brick( new_dset , ivout , output_datum , dfim ) ;
01722 DSET_BRICK_FACTOR(new_dset,ivout) = (fimfac != 0.0) ? 1.0/fimfac : 0.0 ;
01723 }
01724 break ;
01725 }
01726
01727
01728
01729 if( MRG_cflag_gthr != THFLAG_NONE ){
01730 short * dfim ;
01731
01732 dfim = (short *) XtMalloc( sizeof(short) * nxyz ) ;
01733 thrfac = FUNC_COR_SCALE_SHORT ;
01734 EDIT_coerce_scale_type( nxyz,thrfac , MRI_float,ggfim , MRI_short,dfim ) ;
01735 myXtFree( ggfim ) ;
01736 EDIT_substitute_brick( new_dset , DSET_THRESH_VALUE(new_dset) ,
01737 MRI_short , dfim ) ;
01738 DSET_BRICK_FACTOR(new_dset,1) = 1.0/thrfac ;
01739
01740
01741
01742
01743 if( num_fico == num_dset ){
01744 (void) EDIT_dset_items( new_dset , ADN_stat_aux,thr_stataux , ADN_none ) ;
01745
01746
01747
01748 } else {
01749 EDIT_dset_items( new_dset , ADN_func_type,FUNC_THR_TYPE , ADN_none ) ;
01750 }
01751 }
01752
01753 ivout++ ;
01754
01755 }
01756
01757
01758
01759
01760 tross_Make_History( "3dmerge" , argc , argv , new_dset ) ;
01761 THD_load_statistics( new_dset ) ;
01762 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01763 if( ! MRG_be_quiet )
01764 fprintf(stderr,"-- Wrote merged dataset: %s\n" , DSET_BRIKNAME(new_dset) ) ;
01765 exit(0) ;
01766 }
