00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2001, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   This program takes bucket sub-bricks and creates a fim (fico, fitt, 
00010   fift, ...) dataset.  This program was adapted from 3dbucket.c.
00011 
00012 
00013   File:    3dbuc2fim.c
00014   Author:  B. Douglas Ward
00015   Date:    18 March 1998
00016 
00017   Mod:     Added call to AFNI_logger.
00018   Date:    15 August 2001
00019 
00020 */
00021 
00022 /*---------------------------------------------------------------------------*/
00023 
00024 #define PROGRAM_NAME "3dbuc2fim"                     /* name of this program */
00025 #define PROGRAM_AUTHOR "B. D. Ward"                        /* program author */
00026 #define PROGRAM_INITIAL "18 March 1998"   /* date of initial program release */
00027 #define PROGRAM_LATEST  "15 August 2001"  /* date of latest program revision */
00028 
00029 /*---------------------------------------------------------------------------*/
00030 
00031 
00032 #include "mrilib.h"
00033 
00034 #ifndef myXtFree
00035 #   define myXtFree(xp) (XtFree((char *)(xp)) , (xp)=NULL)
00036 #endif
00037 
00038 /*-------------------------- global data --------------------------*/
00039 
00040 static THD_3dim_dataset_array * B2F_dsar  = NULL ;
00041 static XtPointer_array        * B2F_subv  = NULL ;
00042 static int                      B2F_nvox  = -1 ;
00043 static int                      B2F_verb  = 0 ;
00044 static int                      B2F_func_type  = -1;
00045 
00046 static char B2F_output_prefix[THD_MAX_PREFIX] = "buc2fim" ;
00047 static char B2F_session[THD_MAX_NAME]         = "./"   ;
00048 
00049 #define NSUBV(id)   ( ((int *)B2F_subv->ar[(id)])[0]      )
00050 #define SUBV(id,jj) ( ((int *)B2F_subv->ar[(id)])[(jj)+1] )
00051 #define DSUB(id)    DSET_IN_3DARR(B2F_dsar,(id))
00052 
00053 /*--------------------------- prototypes ---------------------------*/
00054 
00055 void B2F_read_opts( int , char ** ) ;
00056 void B2F_Syntax(void) ;
00057 int * B2F_get_subv( int , char * ) ;
00058 
00059 /*--------------------------------------------------------------------
00060    read the arguments, load the global variables
00061 ----------------------------------------------------------------------*/
00062 
00063 void B2F_read_opts( int argc , char * argv[] )
00064 {
00065    int nopt = 1 , ii ;
00066    char dname[THD_MAX_NAME] ;
00067    char subv[THD_MAX_NAME] ;
00068    char * cpt ;
00069    THD_3dim_dataset * dset ;
00070    int * svar ;
00071    char * str;
00072    int ok, ilen, nlen;
00073 
00074    INIT_3DARR(B2F_dsar) ;
00075    INIT_XTARR(B2F_subv) ;
00076 
00077    while( nopt < argc ){
00078 
00079       /**** -prefix prefix ****/
00080 
00081       if( strncmp(argv[nopt],"-prefix",6) == 0 ||
00082           strncmp(argv[nopt],"-output",6) == 0   ){
00083          nopt++ ;
00084          if( nopt >= argc ){
00085             fprintf(stderr,"need argument after -prefix!\n") ; exit(1) ;
00086          }
00087          MCW_strncpy( B2F_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00088          continue ;
00089       }
00090 
00091       /**** -session directory ****/
00092 
00093       if( strncmp(argv[nopt],"-session",6) == 0 ){
00094          nopt++ ;
00095          if( nopt >= argc ){
00096             fprintf(stderr,"need argument after -session!\n") ; exit(1) ;
00097          }
00098          MCW_strncpy( B2F_session , argv[nopt++] , THD_MAX_NAME ) ;
00099          continue ;
00100       }
00101 
00102       if( strncmp(argv[nopt],"-verb",5) == 0 ){
00103          B2F_verb = 1 ;
00104          nopt++ ; continue ;
00105       }
00106 
00107       
00108       if( argv[nopt][0] == '-' ){
00109          fprintf(stderr,"Unknown option: %s\n",argv[nopt]) ; exit(1) ;
00110       }
00111 
00112       /**** read dataset ****/
00113 
00114       cpt = strstr(argv[nopt],"[") ;
00115       if( cpt == NULL ){
00116          strcpy(dname,argv[nopt]) ;
00117          subv[0] = '\0' ;
00118       } else if( cpt == argv[nopt] ){
00119          fprintf(stderr,"illegal dataset specifier: %s\n",argv[nopt]) ;
00120          exit(1) ;
00121       } else {
00122          ii = cpt - argv[nopt] ;
00123          memcpy(dname,argv[nopt],ii) ; dname[ii] = '\0' ;
00124          strcpy(subv,cpt) ;
00125       }
00126       nopt++ ;
00127 
00128       dset = THD_open_one_dataset( dname ) ;
00129       if( dset == NULL ){
00130          fprintf(stderr,"can't open dataset %s\n",dname) ; exit(1) ;
00131       }
00132       THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ;
00133 
00134       ii = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz ;
00135       if( B2F_nvox < 0 ){
00136          B2F_nvox = ii ;
00137       } else if( ii != B2F_nvox ){
00138          fprintf(stderr,"dataset %s differs in size from others\n",dname);
00139          exit(1) ;
00140       }
00141       ADDTO_3DARR(B2F_dsar,dset) ;
00142 
00143       svar = B2F_get_subv( DSET_NVALS(dset) , subv ) ;
00144       if( svar == NULL || svar[0] <= 0 ){
00145          fprintf(stderr,"can't decipher index codes from %s%s\n",dname,subv) ;
00146          exit(1) ;
00147       }
00148       ADDTO_XTARR(B2F_subv,svar) ;
00149 
00150    }  /* end of loop over command line arguments */
00151 
00152    return ;
00153 }
00154 
00155 /*------------------------------------------------------------------*/
00156 
00157 int * B2F_get_subv( int nvals , char * str )
00158 {
00159    int * subv = NULL ;
00160    int ii , ipos , nout , slen ;
00161    int ibot,itop,istep , nused ;
00162    char * cpt ;
00163 
00164    /* Meaningless input? */
00165 
00166    if( nvals < 1 ) return NULL ;
00167 
00168    /* No selection list ==> select it all */
00169 
00170    if( str == NULL || str[0] == '\0' ){
00171       subv = (int *) XtMalloc( sizeof(int) * (nvals+1) ) ;
00172       subv[0] = nvals ;
00173       for( ii=0 ; ii < nvals ; ii++ ) subv[ii+1] = ii ;
00174       return subv ;
00175    }
00176 
00177    /* skip initial '[' */
00178 
00179    subv    = (int *) XtMalloc( sizeof(int) * 2 ) ;
00180    subv[0] = nout = 0 ;
00181 
00182    ipos = 0 ;
00183    if( str[ipos] == '[' ) ipos++ ;
00184 
00185    /*** loop through each sub-selector until end of input ***/
00186 
00187    slen = strlen(str) ;
00188    while( ipos < slen && str[ipos] != ']' ){
00189 
00190       /** get starting value **/
00191 
00192       if( str[ipos] == '$' ){  /* special case */
00193          ibot = nvals-1 ; ipos++ ;
00194       } else {                 /* decode an integer */
00195          ibot = strtol( str+ipos , &cpt , 10 ) ;
00196          if( ibot < 0 ){ myXtFree(subv) ; return NULL ; }
00197          if( ibot >= nvals ) ibot = nvals-1 ;
00198          nused = (cpt-(str+ipos)) ;
00199          if( ibot == 0 && nused == 0 ){ myXtFree(subv) ; return NULL ; }
00200          ipos += nused ;
00201       }
00202 
00203       /** if that's it for this sub-selector, add one value to list **/
00204 
00205       if( str[ipos] == ',' || str[ipos] == '\0' || str[ipos] == ']' ){
00206          nout++ ;
00207          subv = (int *) XtRealloc( (char *)subv , sizeof(int) * (nout+1) ) ;
00208          subv[0]    = nout ;
00209          subv[nout] = ibot ;
00210          ipos++ ; continue ;  /* re-start loop at next sub-selector */
00211       }
00212 
00213       /** otherwise, must have '..' or '-' as next inputs **/
00214 
00215       if( str[ipos] == '-' ){
00216          ipos++ ;
00217       } else if( str[ipos] == '.' && str[ipos+1] == '.' ){
00218          ipos++ ; ipos++ ;
00219       } else {
00220          myXtFree(subv) ; return NULL ;
00221       }
00222 
00223       /** get ending value for loop now **/
00224 
00225       if( str[ipos] == '$' ){  /* special case */
00226          itop = nvals-1 ; ipos++ ;
00227       } else {                 /* decode an integer */
00228          itop = strtol( str+ipos , &cpt , 10 ) ;
00229          if( itop < 0 ){ myXtFree(subv) ; return NULL ; }
00230          if( itop >= nvals ) itop = nvals-1 ;
00231          nused = (cpt-(str+ipos)) ;
00232          if( itop == 0 && nused == 0 ){ myXtFree(subv) ; return NULL ; }
00233          ipos += nused ;
00234       }
00235 
00236       /** set default loop step **/
00237 
00238       istep = (ibot <= itop) ? 1 : -1 ;
00239 
00240       /** check if we have a non-default loop step **/
00241 
00242       if( str[ipos] == '(' ){  /* decode an integer */
00243          ipos++ ;
00244          istep = strtol( str+ipos , &cpt , 10 ) ;
00245          if( istep == 0 ){ myXtFree(subv) ; return NULL ; }
00246          nused = (cpt-(str+ipos)) ;
00247          ipos += nused ;
00248          if( str[ipos] == ')' ) ipos++ ;
00249       }
00250 
00251       /** add values to output **/
00252 
00253       for( ii=ibot ; (ii-itop)*istep <= 0 ; ii += istep ){
00254          nout++ ;
00255          subv = (int *) XtRealloc( (char *)subv , sizeof(int) * (nout+1) ) ;
00256          subv[0]    = nout ;
00257          subv[nout] = ii ;
00258       }
00259 
00260       /** check if we have a comma to skip over **/
00261 
00262       if( str[ipos] == ',' ) ipos++ ;
00263 
00264    }  /* end of loop through selector string */
00265 
00266    return subv ;
00267 }
00268 
00269 /*------------------------------------------------------------------*/
00270 /*
00271    Routine to display 3dbuc2fim help menu.
00272 */
00273 
00274 
00275 void B2F_Syntax(void)
00276 {
00277   printf 
00278     (
00279      "This program converts bucket sub-bricks to fim (fico, fitt, fift, ...)\n"
00280      "type dataset.                                                       \n\n"
00281      "Usage:                                                              \n\n"
00282      "3dbuc2fim  -prefix pname  d1+orig[index]                              \n"
00283      "     This produces a fim dataset.                                   \n\n"
00284      " -or-                                                               \n\n"
00285      "3dbuc2fim  -prefix pname  d1+orig[index1]  d2+orig[index2]            \n"
00286      "     This produces a fico (fitt, fift, ...) dataset,                  \n"
00287      "     depending on the statistic type of the 2nd subbrick,             \n"
00288      "     with   d1+orig[index1] -> intensity sub-brick of pname           \n"
00289      "            d2+orig[index2] -> threshold sub-brick of pname         \n\n"
00290      " -or-                                                               \n\n"
00291      "3dbuc2fim  -prefix pname  d1+orig[index1,index2]                      \n"
00292      "     This produces a fico (fitt, fift, ...) dataset,                  \n"
00293      "     depending on the statistic type of the 2nd subbrick,             \n"
00294      "     with   d1+orig[index1] -> intensity sub-brick of pname           \n"
00295      "            d1+orig[index2] -> threshold sub-brick of pname         \n\n"
00296      "where the options are:\n"
00297    ) ;
00298 
00299    printf(
00300     "     -prefix pname = Use 'pname' for the output dataset prefix name.\n"
00301     " OR  -output pname     [default='buc2fim']\n"
00302     "\n"
00303     "     -session dir  = Use 'dir' for the output dataset session directory.\n"
00304     "                       [default='./'=current working directory]\n"
00305     "     -verb         = Print out some verbose output as the program\n"
00306     "                       proceeds \n"
00307     "\n"
00308     "Command line arguments after the above are taken as input datasets.  \n"
00309     "A dataset is specified using one of these forms:\n"
00310     "   'prefix+view', 'prefix+view.HEAD', or 'prefix+view.BRIK'.\n"
00311     "Sub-brick indexes start at 0. \n"
00312     "\n"
00313     "N.B.: The sub-bricks are output in the order specified, which may\n"
00314     " not be the order in the original datasets.  For example, using\n"
00315     "           fred+orig[5,3]\n"
00316     " will cause the sub-brick #5 in fred+orig to be output as the intensity\n"
00317     " sub-brick, and sub-brick #3 to be output as the threshold sub-brick \n"
00318     " in the new dataset.\n"
00319     "\n"
00320     "N.B.: The '$', '(', ')', '[', and ']' characters are special to\n"
00321     " the shell, so you will have to escape them.  This is most easily\n"
00322     " done by putting the entire dataset plus selection list inside\n"
00323     " single quotes, as in 'fred+orig[5,9]'.\n"
00324     "\n"
00325     ) ;
00326 
00327    exit(0) ;
00328 }
00329 
00330 /*------------------------------------------------------------------*/
00331 
00332 int main( int argc , char * argv[] )
00333 {
00334    int ninp , ids , nv , iv,jv,kv , ivout , new_nvals ;
00335    THD_3dim_dataset * new_dset=NULL , * dset ;
00336    char buf[256] ;
00337 
00338   /*----- Identify software -----*/
00339   printf ("\n\n");
00340   printf ("Program: %s \n", PROGRAM_NAME);
00341   printf ("Author:  %s \n", PROGRAM_AUTHOR);
00342   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
00343   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
00344   printf ("\n");
00345 
00346 
00347    /*** read input options ***/
00348 
00349    if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) B2F_Syntax() ;
00350 
00351    mainENTRY("3dbuc2fim main"); machdep(); AFNI_logger(PROGRAM_NAME,argc,argv);
00352 
00353    B2F_read_opts( argc , argv ) ;
00354 
00355    /*** create new dataset (empty) ***/
00356    ninp = B2F_dsar->num ;
00357    if( ninp < 1 ){
00358       fprintf(stderr,"*** No input datasets?\n") ; exit(1) ;
00359    }
00360 
00361    new_nvals = 0 ;
00362    for( ids=0 ; ids < ninp ; ids++ ) new_nvals += NSUBV(ids) ;
00363 
00364    /*----- Check for acceptable number of sub-bricks -----*/
00365    if (new_nvals < 1)
00366      { fprintf(stderr,"*** Less than 1 sub-brick specified\n") ; exit(1) ; }
00367    if (new_nvals > 2)
00368      { fprintf(stderr,"*** More than 2 sub-bricks specified\n") ; exit(1) ; }
00369 
00370 
00371    if( B2F_verb ) printf("-verb: output will have %d sub-bricks\n",new_nvals) ;
00372 
00373    new_dset = EDIT_empty_copy( DSUB(0) ) ;
00374 
00375    if( ninp == 1 ) tross_Copy_History( DSUB(0) , new_dset ) ;
00376    tross_Make_History( "3dbuc2fim" , argc,argv , new_dset ) ;
00377 
00378    /*-----  Set default value for function type. This will be changed later,
00379             if the second sub-brick has a statistic type.  -----*/
00380    if (new_nvals == 1)  
00381      B2F_func_type = FUNC_FIM_TYPE;
00382    else
00383      B2F_func_type = FUNC_THR_TYPE;
00384 
00385 
00386    EDIT_dset_items (new_dset ,
00387                     ADN_prefix        , B2F_output_prefix ,
00388                     ADN_directory_name, B2F_session ,
00389                     ADN_type          , HEAD_FUNC_TYPE,
00390                     ADN_func_type     , B2F_func_type,
00391                     ADN_ntt           , 0 , 
00392                     ADN_nvals         , new_nvals ,
00393                     ADN_none ) ; 
00394    
00395 
00396    if( THD_is_file(DSET_HEADNAME(new_dset)) ){
00397      fprintf(stderr,"*** Fatal error: file %s already exists!\n",
00398              DSET_HEADNAME(new_dset) ) ;
00399      exit(1) ;
00400    }
00401 
00402    THD_force_malloc_type( new_dset->dblk , DATABLOCK_MEM_MALLOC ) ;
00403 
00404    /*** loop over input datasets ***/
00405 
00406    if( ninp > 1 ) myXtFree( new_dset->keywords ) ;
00407 
00408    ivout = 0 ;
00409    for( ids=0 ; ids < ninp ; ids++ ){
00410       dset = DSUB(ids) ;
00411       nv   = NSUBV(ids) ;
00412 
00413       DSET_load(dset) ;
00414       if( ! DSET_LOADED(dset) ){
00415         fprintf(stderr,"*** Fatal error: can't load data from %s\n",
00416                 DSET_FILECODE(dset)) ;
00417         exit(1) ;
00418       }
00419   
00420 
00421       /** loop over sub-bricks to output **/
00422 
00423       for( iv=0 ; iv < nv ; iv++ ){
00424          jv = SUBV(ids,iv) ;                /* which sub-brick to use */
00425 
00426          EDIT_substitute_brick( new_dset , ivout ,
00427                                 DSET_BRICK_TYPE(dset,jv) , DSET_ARRAY(dset,jv) ) ;
00428          
00429          /*----- If this sub-brick is from a bucket dataset, 
00430            preserve the label for this sub-brick -----*/
00431          if (dset->func_type == FUNC_BUCK_TYPE)
00432            sprintf (buf, "%s", DSET_BRICK_LABEL(dset,jv));
00433          else
00434            sprintf(buf,"%.12s[%d]",DSET_PREFIX(dset),jv) ;
00435          EDIT_dset_items( new_dset , ADN_brick_label_one+ivout, buf , ADN_none ) ;
00436          
00437          sprintf(buf,"%s[%d]",DSET_FILECODE(dset),jv) ;
00438          EDIT_dset_items(
00439                          new_dset, ADN_brick_keywords_replace_one+ivout, buf, ADN_none ) ;
00440          
00441          EDIT_dset_items(
00442                          new_dset ,
00443                          ADN_brick_fac_one            +ivout, DSET_BRICK_FACTOR(dset,jv),
00444                          ADN_brick_keywords_append_one+ivout, DSET_BRICK_KEYWORDS(dset,jv) ,
00445                          ADN_none ) ;
00446          
00447          /** possibly write statistical parameters for this sub-brick **/
00448          
00449          kv = DSET_BRICK_STATCODE(dset,jv) ;
00450          
00451          if( FUNC_IS_STAT(kv) ){ /* input sub-brick has stat params */
00452            
00453            int npar = MAX_STAT_AUX , lv ;
00454            float * par = (float *) malloc( sizeof(float) * (npar) ) ;
00455            float * sax = DSET_BRICK_STATAUX(dset,jv) ;
00456            for( lv=0 ; lv < npar ; lv++ )
00457              par[lv] = (sax != NULL) ? sax[lv] : 0.0 ;
00458              
00459            if (ivout == 1)
00460              {
00461                EDIT_dset_items(new_dset ,
00462                                ADN_func_type     , kv,               
00463                                ADN_stat_aux, par ,
00464                                ADN_none ) ;
00465              }
00466            
00467            free(par) ;
00468              
00469              /* 2: if the input dataset has statistical parameters */
00470 
00471          } else if( ISFUNC(dset)                        &&   /* dset has stat */
00472                     FUNC_IS_STAT(dset->func_type)       &&   /* params        */
00473                     jv == FUNC_ival_thr[dset->func_type]  ){ /* thr sub-brick */
00474            
00475            int npar , lv ;
00476            float * par , * sax ;
00477            kv  = dset->func_type ;
00478            npar = MAX_STAT_AUX ;
00479            par  = (float *) malloc( sizeof(float) * (npar+2) ) ;
00480            sax  = dset->stat_aux ;
00481            for( lv=0 ; lv < npar ; lv++ )
00482              par[lv] = (sax != NULL) ? sax[lv] : 0.0 ;
00483              
00484 
00485            if (ivout == 1)
00486              {
00487                for( lv=0 ; lv < npar+2 ; lv++ )
00488                  printf ("par[%d] = %f \n", lv, par[lv]);
00489                EDIT_dset_items(new_dset ,
00490                                ADN_func_type     , kv,               
00491                                ADN_stat_aux, par ,
00492                                ADN_none ) ;
00493              }
00494            
00495            free(par) ;
00496          }
00497          
00498          /** print a message? **/
00499          
00500          if( B2F_verb ) printf("-verb: copied %s[%d] into %s[%d]\n" ,
00501                                DSET_FILECODE(dset) , jv ,
00502                                DSET_FILECODE(new_dset) , ivout ) ;
00503          ivout++ ;
00504       }
00505 
00506       /** loop over all bricks in input dataset and
00507         unload them if they aren't going into the output
00508         (not required, but is done to economize on memory) **/
00509       
00510       if( nv < DSET_NVALS(dset) ){
00511         
00512         for( kv=0 ; kv < DSET_NVALS(dset) ; kv++ ){  /* all input sub-bricks */
00513           for( iv=0 ; iv < nv ; iv++ ){             /* all output sub-bricks */
00514             jv = SUBV(ids,iv) ;
00515             if( jv == kv ) break ;                 /* input matches output */
00516           }
00517           if( iv == nv ){
00518             mri_free( DSET_BRICK(dset,kv) ) ;
00519 #if 0
00520             if( B2F_verb ) printf("-verb: unloaded unused %s[%d]\n" ,
00521                                   DSET_FILECODE(dset) , kv ) ;
00522 #endif
00523           }
00524         }
00525       }
00526       
00527    } /* end of loop over input datasets */
00528 
00529    
00530    if( B2F_verb ) fprintf(stderr,"-verb: loading statistics\n") ;
00531    THD_load_statistics( new_dset ) ;
00532    THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
00533    if( B2F_verb ) fprintf(stderr,"-verb: wrote output: %s\n",DSET_BRIKNAME(new_dset)) ;
00534    
00535 
00536    exit(0) ;
00537 }
00538