00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"
00008 
00009 /* prototypes */
00010 
00011 static int THD_setup_mastery( THD_3dim_dataset * , int * ) ;
00012 static THD_3dim_dataset * THD_open_3dcalc( char * ) ;
00013 
00014 /*-----------------------------------------------------------------
00015    11 Jan 1999: Open a dataset, allowing for possible mastering.
00016    21 Feb 2001: Allow for <a..b> sub-ranging as well.
00017    26 Jul 2004: Change THD_setup_mastery to return int.    [rickr]
00018                 Add THD_copy_dset_subs() function.
00019 -------------------------------------------------------------------*/
00020 
00021 THD_3dim_dataset * THD_open_dataset( char *pathname )
00022 {
00023    THD_3dim_dataset *dset ;
00024    char dname[THD_MAX_NAME] , subv[THD_MAX_NAME] ;
00025    char *cpt , *bpt ;
00026    int  *ivlist=NULL ;
00027    int    ii , jj , kk ;
00028    float  bot=1.0 , top=0.0 ;
00029 
00030 ENTRY("THD_open_dataset") ;
00031 
00032    /*-- sanity check --*/
00033 
00034    if( pathname == NULL            ||
00035        (ii=strlen(pathname)) == 0  ||
00036        pathname[ii-1]        == '/'  ) RETURN(NULL) ;
00037 
00038    /*-- 23 Mar 2001: perhaps get from across the Web --*/
00039 
00040    if( strncmp(pathname,"http://",7) == 0 ||
00041        strncmp(pathname,"ftp://" ,6) == 0   ){
00042 
00043       dset = THD_fetch_dataset( pathname ) ;
00044       RETURN(dset) ;
00045    }
00046 
00047    /*-- 17 Mar 2000: check if this is a 3dcalc() run --*/
00048 
00049    if( strncmp(pathname,"3dcalc(",7)  == 0 ){
00050      dset = THD_open_3dcalc( pathname ) ;
00051      RETURN(dset) ;
00052    }
00053 
00054    /*-- 04 Mar 2003: allow input of .1D files   --*/
00055    /*--              which deals with [] itself --*/
00056 
00057    if( strstr(pathname,".1D") != NULL ){
00058      dset = THD_open_1D( pathname ) ;
00059      if( dset != NULL ) RETURN(dset) ;
00060    }
00061 
00062    /*-- 04 Aug 2004: allow input of a list of dataset, separated by spaces --*/
00063 
00064    if( strchr(pathname,' ') != NULL ){
00065      dset = THD_open_tcat( pathname ) ;
00066      RETURN(dset) ;
00067    }
00068 
00069    /*-- find the opening "[" and/or "<" --*/
00070 
00071    cpt = strstr(pathname,"[") ;
00072    bpt = strstr(pathname,"<") ;  /* 21 Feb 2001 */
00073 
00074    if( cpt == NULL && bpt == NULL ){            /* no "[" or "<"  */
00075      dset = THD_open_one_dataset( pathname ) ;  /* ==> open      */
00076      RETURN(dset) ;                             /*     normally */
00077    }
00078 
00079    if( cpt == pathname || bpt == pathname ) RETURN(NULL);  /* error */
00080 
00081    /* copy dataset filename to dname and selector string to subv */
00082 
00083    ii = (cpt != NULL ) ? cpt - pathname : 999999 ;
00084    jj = (bpt != NULL ) ? bpt - pathname : 999999 ;
00085    kk = MIN(ii,jj) ;
00086    memcpy(dname,pathname,kk) ; dname[kk] = '\0' ;
00087 
00088    if( STRING_HAS_SUFFIX(dname,".mnc")    ||
00089        STRING_HAS_SUFFIX(dname,".hdr")    ||
00090        STRING_HAS_SUFFIX(dname,".nia")    ||
00091        STRING_HAS_SUFFIX(dname,".nii")    ||
00092        STRING_HAS_SUFFIX(dname,".nii.gz") ||
00093        STRING_HAS_SUFFIX(dname,".mri")    ||
00094        STRING_HAS_SUFFIX(dname,".svl")      ){
00095 
00096      fprintf(stderr,"** Can't use selectors on dataset: %s\n",pathname) ;
00097      RETURN(NULL) ;
00098    }
00099 
00100    /* open the dataset */
00101 
00102    dset = THD_open_one_dataset( dname ) ;
00103    if( dset == NULL ) RETURN(NULL) ;
00104 
00105    /* parse the sub-brick selector string (if any) */
00106 
00107    if( cpt != NULL ){
00108      char *qpt ;
00109      strcpy(subv,cpt) ;
00110      qpt = strstr(subv,"<") ; if( qpt != NULL ) *qpt = '\0' ;
00111      ivlist = MCW_get_intlist( DSET_NVALS(dset) , subv ) ;
00112    }
00113    if( ivlist == NULL ){
00114      if( cpt != NULL )
00115        fprintf(stderr,"** WARNING: bad sub-brick selector => using [0..%d]\n",
00116                DSET_NVALS(dset)-1) ;
00117      ivlist = (int *) malloc(sizeof(int)*(DSET_NVALS(dset)+1)) ;
00118      ivlist[0] = DSET_NVALS(dset) ;
00119      for( kk=0 ; kk < ivlist[0] ; kk++ ) ivlist[kk+1] = kk ;
00120    }
00121 
00122    /* 21 Feb 2001: if present, load the sub-range data */
00123 
00124    if( bpt != NULL ){
00125       char *dpt = strstr(bpt,"..") ;
00126 #if 0
00127 fprintf(stderr,"bpt=%s\n",bpt) ;
00128 #endif
00129       if( dpt != NULL ){
00130 #if 0
00131 fprintf(stderr,"dpt=%s\n",dpt) ;
00132 #endif
00133          kk  = sscanf( bpt+1 , "%f" , &bot ) ;
00134          kk += sscanf( dpt+2 , "%f" , &top ) ;
00135          if( kk == 2 && bot <= top ){
00136             dset->dblk->master_bot = bot ;
00137             dset->dblk->master_top = top ;
00138          } else {
00139             dset->dblk->master_bot = 1.0 ;
00140             dset->dblk->master_top = 0.0 ;
00141          }
00142       }
00143    }
00144 
00145    /* modify the dataset according to the selector string */
00146 
00147    THD_setup_mastery( dset , ivlist ) ;
00148    free(ivlist) ;
00149 
00150    RETURN(dset) ;
00151 }
00152 
00153 /*-----------------------------------------------------------------
00154    Set up a dataset for being mastered; that is, reading only
00155    a subset of sub-bricks from the master .BRIK file.
00156 -------------------------------------------------------------------*/
00157 
00158 static int THD_setup_mastery( THD_3dim_dataset *dset , int *ivlist )
00159 {
00160    int ibr , old_nvals , new_nvals ;
00161    THD_datablock *dblk ;
00162    int *btype , *ivl ;
00163 
00164    float * old_brick_fac  ;
00165    int *   old_brick_bytes ;
00166    char ** old_brick_lab  ;
00167    char ** old_brick_keywords ;
00168    int *   old_brick_statcode ;
00169    float **old_brick_stataux ;
00170 
00171 ENTRY("THD_setup_mastery") ;
00172 
00173    /** sanity checks **/
00174 
00175    if( ! ISVALID_DSET(dset) || ivlist == NULL || ivlist[0] <= 0 ) RETURN(1) ;
00176 
00177    new_nvals = ivlist[0] ;
00178    ivl       = ivlist + 1 ;
00179    dblk      = dset->dblk ;
00180    old_nvals = dblk->nvals ;
00181 
00182    ibr = THD_count_databricks(dblk) ; if( ibr > 0 ) RETURN(2) ;
00183 
00184    for( ibr=0 ; ibr < new_nvals ; ibr++ )
00185       if( ivl[ibr] < 0 || ivl[ibr] >= old_nvals ) RETURN(3) ;
00186 
00187    /** save pointers to old datablock stuff **/
00188 
00189    old_brick_fac      = dblk->brick_fac      ; dblk->brick_fac      = NULL ;
00190    old_brick_bytes    = dblk->brick_bytes    ; dblk->brick_bytes    = NULL ;
00191    old_brick_lab      = dblk->brick_lab      ; dblk->brick_lab      = NULL ;
00192    old_brick_keywords = dblk->brick_keywords ; dblk->brick_keywords = NULL ;
00193    old_brick_statcode = dblk->brick_statcode ; dblk->brick_statcode = NULL ;
00194    old_brick_stataux  = dblk->brick_stataux  ; dblk->brick_stataux  = NULL ;
00195 
00196    /** setup new dataset brick structure **/
00197 
00198    dblk->diskptr->nvals = dblk->nvals = new_nvals ;
00199    dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
00200 
00201    if( dset->taxis != NULL ){                /* must fix time axis */
00202       if( new_nvals == 1 ){                  /* no time dependence */
00203          myXtFree( dset->taxis->toff_sl ) ;
00204          myXtFree( dset->taxis ) ;
00205       } else {                               /* different number of times */
00206          dset->taxis->ntt = new_nvals ;
00207       }
00208    } else {                                  /* 21 Feb 2001: change to bucket type */
00209 
00210       if( ISANAT(dset) && !ISANATBUCKET(dset) )
00211          EDIT_dset_items( dset , ADN_func_type,ANAT_BUCK_TYPE , ADN_none ) ;
00212       else if( ISFUNC(dset) && !ISFUNCBUCKET(dset) )
00213          EDIT_dset_items( dset , ADN_func_type,FUNC_BUCK_TYPE , ADN_none ) ;
00214 
00215    }
00216 
00217    /* redo brick_fac */
00218 
00219    dblk->brick_fac = (float *) XtMalloc( sizeof(float) * new_nvals ) ;
00220    for( ibr=0 ; ibr < new_nvals ; ibr++ )
00221       dblk->brick_fac[ibr] = old_brick_fac[ivl[ibr]] ;
00222 
00223    /* redo brick and brick_bytes */
00224 
00225    btype = (int *) malloc( sizeof(int) * new_nvals ) ;
00226    for( ibr=0 ; ibr < new_nvals ; ibr++ )
00227       btype[ibr] = DBLK_BRICK_TYPE(dblk,ivl[ibr]) ;
00228    THD_init_datablock_brick( dblk , new_nvals , btype ) ;
00229    free(btype) ;
00230 
00231    /* redo brick_lab */
00232 
00233    if( old_brick_lab != NULL ){
00234      for( ibr=0 ; ibr < new_nvals ; ibr++ )
00235        THD_store_datablock_label( dblk , ibr , old_brick_lab[ivl[ibr]] ) ;
00236    }
00237 
00238    /* redo brick_keywords */
00239 
00240    if( old_brick_keywords != NULL ){
00241      for( ibr=0 ; ibr < new_nvals ; ibr++ )
00242        THD_store_datablock_keywords( dblk , ibr , old_brick_keywords[ivl[ibr]] ) ;
00243    }
00244 
00245    /* redo brick_statcode and brick_stataux */
00246 
00247    if( old_brick_statcode != NULL ){
00248       for( ibr=0 ; ibr < new_nvals ; ibr++ )
00249          THD_store_datablock_stataux( dblk, ibr, old_brick_statcode[ivl[ibr]] ,
00250                                            999 , old_brick_stataux [ivl[ibr]]  ) ;
00251    }
00252 
00253    /** setup master stuff now **/
00254 
00255    dblk->master_nvals = old_nvals ;
00256    dblk->master_bytes = old_brick_bytes ;
00257    dblk->master_ival  = (int *) XtMalloc( sizeof(int) * new_nvals ) ;
00258    for( ibr=0 ; ibr < new_nvals ; ibr++ ) dblk->master_ival[ibr] = ivl[ibr] ;
00259 
00260    /** destroy old datablock stuff now **/
00261 
00262    myXtFree( old_brick_fac ) ;
00263 
00264    if( old_brick_lab != NULL ){
00265      for( ibr=0 ; ibr < old_nvals ; ibr++ ) myXtFree( old_brick_lab[ibr] ) ;
00266      myXtFree( old_brick_lab ) ;
00267    }
00268 
00269    if( old_brick_keywords != NULL ){
00270      for( ibr=0 ; ibr < old_nvals ; ibr++ ) myXtFree( old_brick_keywords[ibr] ) ;
00271      myXtFree( old_brick_keywords ) ;
00272    }
00273 
00274    if( old_brick_statcode != NULL ) myXtFree( old_brick_statcode ) ;
00275    if( old_brick_stataux  != NULL ){
00276       for( ibr=0 ; ibr < old_nvals ; ibr++ ) myXtFree( old_brick_stataux[ibr] ) ;
00277       myXtFree( old_brick_stataux ) ;
00278    }
00279 
00280    /** if dataset has statistics, rearrange them **/
00281 
00282    if( ISVALID_STATISTIC(dset->stats) ){
00283       THD_statistics * new_stats , * old_stats ;
00284       THD_brick_stats * bsold , * bsnew ;
00285       float bot,top ;
00286 
00287       old_stats = dset->stats ;
00288       new_stats = myXtNew( THD_statistics ) ;
00289       new_stats->type   = STATISTICS_TYPE ;
00290       new_stats->parent = (XtPointer) dset ;
00291       new_stats->bstat  = NULL ;
00292 
00293       bsold = old_stats->bstat ;
00294       bsnew = new_stats->bstat =
00295          (THD_brick_stats *) XtCalloc( new_nvals , sizeof(THD_brick_stats) ) ;
00296 
00297       new_stats->nbstat = new_nvals ;
00298 
00299       for( ibr=0 ; ibr < new_nvals ; ibr++ ){
00300          if( ibr < old_stats->nbstat ) bsnew[ibr] = bsold[ivl[ibr]] ;
00301          else                          INVALIDATE_BSTAT( bsnew[ibr] ) ;
00302       }
00303 
00304       REPLACE_KILL( dset->kl , bsold     , bsnew     ) ;
00305       REPLACE_KILL( dset->kl , old_stats , new_stats ) ;
00306       dset->stats = new_stats ;
00307 
00308       myXtFree(bsold) ; myXtFree(old_stats) ;
00309 
00310       /* 21 Feb 2001: mangle statistics if sub-ranging is used */
00311 
00312       bot = dset->dblk->master_bot ; top = dset->dblk->master_top ;
00313       if( bot <= top ){
00314               if( bot > 0.0 ) bot = 0.0 ;
00315          else if( top < 0.0 ) top = 0.0 ;
00316 
00317          for( ibr=0 ; ibr < new_nvals ; ibr++ ){
00318             if( ISVALID_BSTAT(bsnew[ibr]) ){
00319                if( bsnew[ibr].min < bot ) bsnew[ibr].min = bot ;
00320                if( bsnew[ibr].max > top ) bsnew[ibr].max = top ;
00321             }
00322          }
00323       }
00324    }
00325 
00326    RETURN(0) ;
00327 }
00328 
00329 /*----------------------------------------------------------------------
00330    Run 3dcalc to create a dataset and read it in.
00331    -- RWCox - 17 Mar 2000
00332 ------------------------------------------------------------------------*/
00333 
00334 #include <sys/types.h>
00335 #include <sys/wait.h>
00336 
00337 static THD_3dim_dataset * THD_open_3dcalc( char * pname )
00338 {
00339    int    Argc=1               ,    newArgc=0 , ii,ll  ;
00340    char * Argv[1]={ "3dcalc" } , ** newArgv=NULL ;
00341    char * qname , * tdir , prefix[16] ;
00342    pid_t  child_pid ;
00343    THD_3dim_dataset * dset ;
00344    static int ibase=1 ;
00345 
00346 ENTRY("THD_open_3dcalc") ;
00347 
00348    /*-- remove the "3dcalc(" and the ")" from the input string --*/
00349 
00350    qname = (char *) malloc(sizeof(char)*(strlen(pname)+1024)) ;
00351    strcpy(qname,pname+7) ;
00352    ll = strlen(qname) ;
00353    for( ii=ll-1 ; ii > 0 && qname[ii] != ')' ; ii++ ) ; /* nada */
00354    if( ii == 0 ){ free(qname) ; RETURN(NULL) ; }
00355    qname[ii] = '\0' ;
00356 
00357    /*-- add -session to command string --*/
00358 
00359    tdir = my_getenv("TMPDIR") ;
00360    if( tdir == NULL || strlen(tdir) > 512 ) tdir = "/tmp" ;
00361    strcat(qname," -session ") ; strcat(qname,tdir) ; ll = strlen(tdir) ;
00362 
00363    /*-- add -prefix to command string --*/
00364 
00365    for( ii=ibase ; ii < 9999 ; ii++ ){                    /* dataset name   */
00366       sprintf(prefix,"3dcalc#%04d",ii) ;
00367       if( THD_is_dataset(tdir,prefix,-1) == -1 ) break ;
00368    }
00369    if( ii > 9999 ){
00370      fprintf(stderr,"*** Can't find unused 3dcalc# dataset name in %s!\n",tdir) ;
00371      free(qname) ; RETURN(NULL) ;
00372    }
00373    ibase = ii+1 ;
00374 
00375    strcat(qname," -prefix ") ; strcat(qname,prefix) ;
00376 
00377    strcat(qname," -verbose") ;
00378 
00379    /*-- add a placeholder to be the last argument --*/
00380 
00381    strcat(qname," Zork") ;
00382 
00383    /*-- create the arg list for 3dcalc, starting with program name --*/
00384 
00385    append_string_to_args( qname , Argc , Argv , &newArgc , &newArgv ) ;
00386 
00387    free(qname) ; /* not needed no more */
00388 
00389    /*-- check if arg list was created OK --*/
00390 
00391    if( newArgv == NULL ) RETURN(NULL) ;  /* something bad? */
00392 
00393    if( newArgc < 3 ){                   /* too few args to 3dcalc */
00394      for( ii=0 ; ii < newArgc ; ii++ ) free(newArgv[ii]) ;
00395      free(newArgv) ; RETURN(NULL) ;
00396    }
00397 
00398    /*-- replace placeholder in arg list with NULL pointer --*/
00399 
00400    free( newArgv[newArgc-1] ) ; newArgv[newArgc-1] = NULL ;
00401 
00402    /*-- fork and exec --*/
00403 
00404    fprintf(stderr,"+++ Executing 3dcalc()\n") ;
00405 #if 0
00406 for(ii=0; ii< newArgc-1; ii++) fprintf(stderr," argv[%d]=%s\n",ii,newArgv[ii]);
00407 #endif
00408    child_pid = fork() ;
00409 
00410    if( child_pid == (pid_t)(-1) ){
00411      perror("*** Can't fork 3dcalc()") ;
00412      for( ii=0 ; ii < newArgc-1 ; ii++ ) free(newArgv[ii]) ;
00413      free(newArgv) ; RETURN(NULL) ;
00414    }
00415 
00416    if( child_pid == 0 ){  /*-- I'm the child --*/
00417 
00418      execvp( "3dcalc" , newArgv ) ;        /* should not return */
00419      perror("*** Can't execvp 3dcalc()") ;
00420      _exit(1) ;
00421 
00422    }
00423 
00424    /*-- I'm the parent --*/
00425 
00426    (void) waitpid( child_pid , NULL , 0 ) ; /* wait for child to exit */
00427 
00428    ii = THD_is_dataset( tdir , prefix , -1 ) ;
00429    if( ii == -1 ){
00430      fprintf(stderr,"*** 3dcalc() failed - no dataset created\n") ;
00431      RETURN(NULL) ;
00432    }
00433    qname = THD_dataset_headname( tdir , prefix , ii ) ;
00434    dset = THD_open_one_dataset( qname ) ;  /* try to read result */
00435 
00436    for( ii=0 ; ii < newArgc-1 ; ii++ ) free(newArgv[ii]) ;  /* toss trash */
00437    free(newArgv) ; free(qname) ;
00438 
00439    if( dset == NULL ){                          /* read failed */
00440      fprintf(stderr,"*** 3dcalc() failed - can't read dataset\n") ;
00441      RETURN(NULL) ;
00442    }
00443 
00444    /* read dataset into memory */
00445 
00446    DSET_mallocize(dset) ; DSET_load(dset) ;
00447    if( !DSET_LOADED(dset) ){                   /* can't read it? */
00448      THD_delete_3dim_dataset( dset , True ) ; /* kill it dead */
00449      fprintf(stderr,"*** 3dcalc() failed - can't load dataset\n") ;
00450      RETURN(NULL) ;
00451    }
00452 
00453    /* lock dataset into memory, delete its files */
00454 
00455    DSET_lock(dset) ;
00456    unlink( dset->dblk->diskptr->header_name ) ;
00457    COMPRESS_unlink( dset->dblk->diskptr->brick_name ) ;
00458 
00459    /* 30 Jul 2003: changes its directory to cwd */
00460 
00461    EDIT_dset_items( dset , ADN_directory_name , "./" , ADN_none ) ;
00462 
00463    RETURN(dset) ;
00464 }
00465 
00466 /*-----------------------------------------------------------------
00467    Copy a list of sub-bricks from a dataset.    26 Jul 2004 [rickr]
00468    The first element of dlist is the number of sub-bricks to copy.
00469 -------------------------------------------------------------------*/
00470 
00471 THD_3dim_dataset * THD_copy_dset_subs( THD_3dim_dataset * din, int * dlist )
00472 {
00473     THD_3dim_dataset * dout;
00474     MRI_TYPE           kind;
00475     char             * newdata;
00476     int                sub, subs;
00477     int                dsize, nxyz, rv;
00478 
00479 ENTRY("THD_copy_dset_subs");
00480 
00481     /* validate inputs */
00482     if ( !din || !dlist )
00483     {
00484         fprintf(stderr, "** THD_copy_dset_subs: bad input (%p,%p)\n",
00485                 din,dlist);
00486         RETURN(NULL);
00487     }
00488 
00489     if ( dlist[0] <= 0 )
00490     {
00491         fprintf(stderr,"** THD_copy_dset_subs: invalid dlist length %d\n",
00492                 dlist[0]);
00493         RETURN(NULL);
00494     }
00495 
00496     dout = EDIT_empty_copy(din);
00497     rv = THD_setup_mastery(dout, dlist);
00498     if ( rv != 0 )
00499     {
00500         fprintf(stderr, "** failure: THD_setup_mastery() returned %d\n", rv);
00501         RETURN(NULL);
00502     }
00503 
00504     /* be sure that we have some data to copy */
00505     DSET_load(din);
00506     if ( ! DSET_LOADED(din) )
00507     {
00508         fprintf(stderr,"** THD_copy_dset_subs: cannot load input dataset\n");
00509         RETURN(NULL);
00510     }
00511 
00512     /* a basic warp is needed if header is written out - PLUTO_add_dset() */
00513     dout->warp  = myXtNew( THD_warp );
00514     *dout->warp = IDENTITY_WARP;
00515     ADDTO_KILL( dout->kl, dout->warp );
00516 
00517     dout->dblk->diskptr->byte_order   = mri_short_order();
00518     dout->dblk->diskptr->storage_mode = STORAGE_BY_BRICK;
00519 
00520     /* now copy all of the sub-bricks */
00521     nxyz = dout->daxes->nxx * dout->daxes->nyy * dout->daxes->nzz;
00522     subs = dlist[0];
00523     for ( sub = 0; sub < subs; sub++ )
00524     {
00525         kind = DSET_BRICK_TYPE(dout, sub);
00526         dsize = mri_datum_size( kind );
00527         if ( (newdata = (char *)malloc( nxyz * dsize )) == NULL )
00528         {
00529             fprintf( stderr, "r frdb: alloc failure: %d bytes!\n",
00530                      nxyz * dsize );
00531             DSET_delete(dout);
00532             RETURN(NULL);
00533         }
00534 
00535         memcpy(newdata,DSET_ARRAY(din,dlist[sub+1]), nxyz*dsize);
00536         EDIT_substitute_brick(dout, sub, kind, (void *)newdata);
00537     }
00538 
00539     dout->dblk->malloc_type = DATABLOCK_MEM_MALLOC;
00540     dout->wod_flag = False;             /* since data is now in memory */
00541 
00542     RETURN(dout);
00543 }