00001 #include "mrilib.h"
00002 
00003 /*******************************************************************/
00004 /********** 04 Mar 2003: Read a 1D file as an AFNI dataset *********/
00005 /*******************************************************************/
00006 
00007 /*-----------------------------------------------------------------*/
00008 /*! Open a 1D file as an AFNI dataset.  Each column is a separate
00009     sub-brick.  Ignore comments at this time.  Data is not loaded
00010     into bricks.
00011 -------------------------------------------------------------------*/
00012 
00013 THD_3dim_dataset * THD_open_1D( char *pathname )
00014 {
00015    THD_3dim_dataset *dset=NULL ;
00016    MRI_IMAGE *flim ;
00017    char prefix[1024] , *ppp , tname[12] , *pn ;
00018    THD_ivec3 nxyz ;
00019    THD_fvec3 dxyz , orgxyz ;
00020    int nvals , lpn , flip=0 ;
00021    FILE *fp ;
00022 
00023 ENTRY("THD_open_1D") ;
00024 
00025    /*-- open input file --*/
00026 
00027    if( pathname == NULL || pathname[0] == '\0' ) RETURN(NULL);
00028 
00029    /*-- check if it is a NIML-ish AFNI dataset;
00030         if so, read it in that way instead of the 1D way [21 Mar 2003] --*/
00031 
00032    if( strchr(pathname,'[') == NULL &&
00033        strchr(pathname,'{') == NULL && strchr(pathname,'\'') == NULL ){
00034 
00035      pn = strdup(pathname) ; fp = fopen(pn,"r") ;
00036 
00037      if( fp == NULL ){
00038        char *p1 = strstr(pn,"1D") ;   /* if can't open .1D, try .3D */
00039        if( p1 != NULL ){
00040          *p1 = '3' ; fp = fopen(pn,"r") ;
00041        }
00042        if( fp == NULL ){
00043          fprintf(stderr,"** THD_open_1D(%s): can't open file\n",pathname);
00044          free(pn); RETURN(NULL);
00045        }
00046      }
00047      memset(prefix,0,133) ; fread(prefix,1,128,fp) ; fclose(fp) ;
00048      if( strstr(prefix,"<AFNI_") != NULL && strstr(prefix,"dataset") != NULL ){
00049        dset = THD_open_3D(pn) ;
00050        if( dset != NULL && strcmp(pathname,pn) != 0 )
00051          fprintf(stderr,"** THD_open_1D(%s): substituted %s\n",pathname,pn) ;
00052        free(pn) ; return dset ;
00053      }
00054    }
00055 
00056    /*-- otherwise, read it into an MRI_IMAGE, then mangle image into dataset --*/
00057 
00058    lpn = strlen(pathname) ; pn = strdup(pathname) ;
00059 
00060    flip = (pn[lpn-1] == '\'') ;     /* 12 Jul 2005: allow for tranposing input */
00061    if( flip ) pn[lpn-1] = '\0' ;
00062 
00063    flim = mri_read_1D(pn) ;
00064    if( flim == NULL ){
00065      fprintf(stderr,"** Can't read 1D dataset file %s\n",pn); free(pn); RETURN(NULL);
00066    }
00067    if( flip ){
00068      MRI_IMAGE *qim = mri_transpose(flim); mri_free(flim); flim = qim;
00069    }
00070 
00071    /*-- make a dataset --*/
00072 
00073    dset = EDIT_empty_copy(NULL) ;        /* default (useless) dataset */
00074 
00075    ppp  = THD_trailname(pn,0) ;                   /* strip directory */
00076    MCW_strncpy( prefix , ppp , THD_MAX_PREFIX ) ;  /* to make prefix */
00077    ppp  = strchr( prefix , '[' ) ;
00078    if( ppp != NULL ) *ppp = '\0' ;
00079 
00080    nxyz.ijk[0] = flim->nx ; dxyz.xyz[0] = 1.0 ;    /* setup axes */
00081    nxyz.ijk[1] = 1        ; dxyz.xyz[1] = 1.0 ;
00082    nxyz.ijk[2] = 1        ; dxyz.xyz[2] = 1.0 ;
00083 
00084    orgxyz.xyz[0] = 0.0 ;                           /* arbitrary origin */
00085    orgxyz.xyz[1] = 0.0 ;
00086    orgxyz.xyz[2] = 0.0 ;
00087 
00088    nvals = flim->ny ;                              /* number of sub-bricks */
00089 
00090    MCW_hash_idcode( pathname , dset ) ; /* 06 May 2005 */
00091    dset->idcode.str[0] = 'A' ;          /* overwrite 1st 3 bytes of IDcode */
00092    dset->idcode.str[1] = '1' ;
00093    dset->idcode.str[2] = 'D' ;
00094 
00095    /* note we accept default orientation (RAI) here
00096       (no use of ADN_xyzorient), since we actually
00097       don't have any spatial structure to speak of */
00098 
00099    EDIT_dset_items( dset ,
00100                       ADN_prefix      , prefix ,
00101                       ADN_datum_all   , MRI_float ,
00102                       ADN_nxyz        , nxyz ,
00103                       ADN_xyzdel      , dxyz ,
00104                       ADN_xyzorg      , orgxyz ,
00105                       ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00106                       ADN_nvals       , nvals ,
00107                       ADN_type        , HEAD_ANAT_TYPE ,
00108                       ADN_func_type   , ANAT_BUCK_TYPE ,
00109                     ADN_none ) ;
00110 
00111    if( nvals > 1 && AFNI_yesenv("AFNI_1D_TIME") ){ /* pretend it is 3D+time */
00112      char *eee = getenv("AFNI_1D_TIME_TR") ;
00113      float dt = 0.0 ;
00114      if( eee != NULL ) dt = strtod(eee,NULL) ;
00115      if( dt <= 0.0   ) dt = 1.0 ;
00116      EDIT_dset_items( dset ,
00117                         ADN_func_type, ANAT_EPI_TYPE ,
00118                         ADN_ntt      , nvals ,
00119                         ADN_ttorg    , 0.0 ,
00120                         ADN_ttdel    , dt ,    /* TR */
00121                         ADN_ttdur    , 0.0 ,
00122                         ADN_tunits   , UNITS_SEC_TYPE ,
00123                       ADN_none ) ;
00124    }
00125 
00126    dset->dblk->diskptr->storage_mode = STORAGE_BY_1D ;
00127    strcpy( dset->dblk->diskptr->brick_name , pathname ) ;
00128 
00129    /*-- purge image data and return the empty dataset */
00130 
00131    mri_free(flim) ; free(pn) ; RETURN(dset) ;
00132 }
00133 
00134 /*-----------------------------------------------------------------*/
00135 /*!  Load a 1D dataset's data into memory.
00136      Called from THD_load_datablock() in thd_loaddblk.c.
00137 -------------------------------------------------------------------*/
00138 
00139 void THD_load_1D( THD_datablock *dblk )
00140 {
00141    THD_diskptr *dkptr ;
00142    MRI_IMAGE *flim ;
00143    int nxyz , nbad,iv,nv ;
00144    float *bar , *flar ;
00145    char *pn ; int lpn , flip=0 ;
00146 
00147 ENTRY("THD_load_1D") ;
00148 
00149    /*-- open input [these errors should never occur] --*/
00150 
00151    if( !ISVALID_DATABLOCK(dblk)                     ||
00152        dblk->diskptr->storage_mode != STORAGE_BY_1D ||
00153        dblk->brick == NULL                            ) EXRETURN ;
00154 
00155    dkptr = dblk->diskptr      ;
00156    nxyz  = dkptr->dimsizes[0] ;
00157    nv    = dkptr->nvals       ;
00158 
00159    if( nxyz*nv > 1000000 ) fprintf(stderr,"++ Reading %s\n",dkptr->brick_name) ;
00160 
00161    pn = strdup(dkptr->brick_name) ; lpn = strlen(pn) ;
00162    flip = (pn[lpn-1] == '\'') ;     /* 12 Jul 2005: allow for tranposing input */
00163    if( flip ) pn[lpn-1] = '\0' ;
00164 
00165    flim = mri_read_1D(pn) ; free(pn) ;
00166 
00167    if( flim == NULL ){
00168      fprintf(stderr,"** THD_load_1D(%s): can't read file!\n",dkptr->brick_name) ;
00169      EXRETURN ;
00170    }
00171    if( flip ){
00172      MRI_IMAGE *qim = mri_transpose(flim); mri_free(flim); flim = qim;
00173    }
00174 
00175    /*-- allocate space for data --*/
00176 
00177    if( nxyz != flim->nx || nv != flim->ny ){
00178      fprintf(stderr,"** THD_load_1D(%s): nx or ny mismatch!\n",dkptr->brick_name) ;
00179      fprintf(stderr,"**  expect nx=%d; got nx=%d\n",nxyz,flim->nx) ;
00180      fprintf(stderr,"**  expect ny=%d; got ny=%d\n",nv  ,flim->ny) ;
00181      mri_free(flim) ; EXRETURN ;
00182    }
00183 
00184    dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
00185 
00186    /** malloc space for each brick separately **/
00187 
00188    for( nbad=iv=0 ; iv < nv ; iv++ ){
00189      if( DBLK_ARRAY(dblk,iv) == NULL ){
00190        bar = AFMALL( float,  DBLK_BRICK_BYTES(dblk,iv) ) ;
00191        mri_fix_data_pointer( bar ,  DBLK_BRICK(dblk,iv) ) ;
00192        if( bar == NULL ) nbad++ ;
00193      }
00194    }
00195 
00196    /** if couldn't get them all, take our ball and go home in a snit **/
00197 
00198    if( nbad > 0 ){
00199      fprintf(stderr,"\n** failed to malloc %d 1D bricks out of %d\n\a",nbad,nv) ;
00200      for( iv=0 ; iv < nv ; iv++ ){
00201        if( DBLK_ARRAY(dblk,iv) != NULL ){
00202          free(DBLK_ARRAY(dblk,iv)) ;
00203          mri_fix_data_pointer( NULL , DBLK_BRICK(dblk,iv) ) ;
00204        }
00205      }
00206      mri_free(flim) ; EXRETURN ;
00207    }
00208 
00209    /** copy data from image to bricks **/
00210 
00211    flar = MRI_FLOAT_PTR(flim) ;
00212 
00213    for( iv=0 ; iv < nv ; iv++ ){
00214      bar = DBLK_ARRAY(dblk,iv) ;
00215      memcpy( bar , flar + iv*nxyz , sizeof(float)*nxyz ) ;
00216    }
00217 
00218    mri_free(flim) ; EXRETURN ;
00219 }
00220 
00221 /*------------------------------------------------------------------*/
00222 /*! Write a dataset to disk as a 1D file.
00223     Called from THD_write_3dim_dataset().
00224 --------------------------------------------------------------------*/
00225 
00226 void THD_write_1D( char *sname, char *pname , THD_3dim_dataset *dset )
00227 {
00228    char fname[THD_MAX_NAME] , *cpt ;
00229    int iv,nv , nx,ny,nz,nxyz,ii,jj,kk ;
00230    FILE *fp=NULL ;
00231    int binflag ; char shp ; float val[1] ;
00232 
00233 ENTRY("THD_write_1D") ;
00234 
00235    if( !ISVALID_DSET(dset) || !DSET_LOADED(dset) ) EXRETURN ;
00236 
00237    nv = DSET_NVALS(dset) ;  /* number of columns */
00238    nx = DSET_NX(dset)    ;
00239    ny = DSET_NY(dset)    ;
00240    nz = DSET_NZ(dset)    ; nxyz = nx*ny*nz ;  /* number of rows */
00241 
00242    /* make up a filename for output (into fname) */
00243 
00244    cpt = DSET_PREFIX(dset) ;
00245    if( (pname != NULL && *pname == '-') ||
00246        (pname == NULL && *cpt   == '-')   ){  /* 12 Jul 2005: to stdout */
00247      fp = stdout ; strcpy(fname,"stdout") ; binflag = 0 ;
00248    }
00249 
00250    if( fp == NULL ){
00251      if( pname != NULL ){        /* have input prefix */
00252        if( sname != NULL ){      /* and input session (directory) */
00253          strcpy(fname,sname) ;
00254          ii = strlen(fname) ; if( fname[ii-1] != '/' ) strcat(fname,"/") ;
00255        } else {
00256          strcpy(fname,"./") ;    /* don't have input session */
00257        }
00258        strcat(fname,pname) ;
00259      } else {                    /* don't have input prefix */
00260        cpt = DSET_PREFIX(dset) ;
00261        if( STRING_HAS_SUFFIX(cpt,".3D") || STRING_HAS_SUFFIX(cpt,".1D") )
00262          strcpy(fname,cpt) ;
00263        else
00264          strcpy(fname,DSET_BRIKNAME(dset)) ;
00265 
00266        cpt = strchr(fname,'[') ;
00267        if( cpt != NULL ) *cpt = '\0' ;                  /* without subscripts! */
00268      }
00269      ii = strlen(fname) ;
00270      if( ii > 10 && strstr(fname,".BRIK") != NULL ){    /* delete .BRIK! */
00271        fname[ii-10] = '\0' ;
00272        if( DSET_IS_1D(dset) || (ny==1 && nz==1) )
00273          strcat(fname,".1D");
00274        else
00275          strcat(fname,".3D");                           /* 21 Mar 2003 */
00276      }
00277 
00278      fp = fopen( fname , "w" ) ; if( fp == NULL ) EXRETURN ;
00279      binflag = STRING_HAS_SUFFIX(fname,".3D") && AFNI_yesenv("AFNI_3D_BINARY") ;
00280    }
00281 
00282    strcpy( dset->dblk->diskptr->brick_name , fname ); /* 12 Jul 2005 */
00283 
00284    /* are we going to write in binary? [03 Jun 2005] */
00285 
00286    shp = (binflag) ? ' ' : '#' ;
00287 
00288    /* write some dataset info as NIML-style header/comments */
00289 
00290    if( fp != stdout )
00291      fprintf(fp,
00292                 "%c <AFNI_3D_dataset\n"
00293                 "%c  self_idcode = \"%s\"\n"
00294                 "%c  ni_type     = \"%d*float\"\n"    /* all columns are floats! */
00295                 "%c  ni_dimen    = \"%d,%d,%d\"\n"
00296                 "%c  ni_delta    = \"%g,%g,%g\"\n"
00297                 "%c  ni_origin   = \"%g,%g,%g\"\n"
00298                 "%c  ni_axes     = \"%s,%s,%s\"\n"
00299              ,
00300                 shp ,
00301                 shp , dset->idcode.str ,
00302                 shp , nv ,
00303                 shp , nx,ny,nz ,
00304                 shp , DSET_DX(dset)  , DSET_DY(dset)  , DSET_DZ(dset)  ,
00305                 shp , DSET_XORG(dset), DSET_YORG(dset), DSET_ZORG(dset),
00306                 shp , ORIENT_shortstr[dset->daxes->xxorient] ,
00307                        ORIENT_shortstr[dset->daxes->yyorient] ,
00308                          ORIENT_shortstr[dset->daxes->zzorient]
00309             ) ;
00310 
00311    if( fp != stdout && HAS_TIMEAXIS(dset) ){
00312      float dt = DSET_TR(dset) ;
00313      if( DSET_TIMEUNITS(dset) == UNITS_MSEC_TYPE ) dt *= 0.001 ;
00314      fprintf(fp , "%c  ni_timestep = \"%g\"\n" , shp,dt ) ;
00315    }
00316 
00317    if( fp != stdout && binflag )
00318      fprintf(fp , "   ni_form     = \"binary.%s\"\n" ,
00319                   (NI_byteorder()==NI_LSB_FIRST) ? "lsbfirst" : "msbfirst" ) ;
00320 
00321    /* do stataux for bricks, if any are present */
00322 
00323    for( ii=iv=0 ; iv < nv ; iv++ )
00324      if( DSET_BRICK_STATCODE(dset,iv) > 0 ) ii++ ;
00325 
00326    if( fp != stdout && ii > 0 ){
00327       fprintf(fp, "%c  ni_stat     = \"",shp) ;
00328       for( iv=0 ; iv < nv ; iv++ ){
00329         ii = DSET_BRICK_STATCODE(dset,iv) ;
00330         if( ii <=0 ){
00331           fprintf(fp,"none") ;
00332         } else {
00333           fprintf(fp,"%s(",NI_stat_distname(ii)) ;
00334           kk = NI_stat_numparam(ii) ;
00335           for( jj=0 ; jj < kk ; jj++ ){
00336             fprintf(fp,"%g",DSET_BRICK_STATPAR(dset,iv,jj)) ;
00337             if( jj < kk-1 ) fprintf(fp,",") ;
00338           }
00339           fprintf(fp,")") ;
00340         }
00341         if( ii < nv-1 ) fprintf(fp,";") ;
00342       }
00343       fprintf(fp,"\"\n") ;
00344    }
00345 
00346    /* close NIML-style header */
00347 
00348    if( fp != stdout ){
00349      if( binflag ) fprintf(fp," >") ;
00350      else          fprintf(fp,"# >\n") ;
00351      fflush(fp) ;
00352    }
00353 
00354    /* now write data */
00355 
00356    for( ii=0 ; ii < nxyz ; ii++ ){  /* loop over voxels */
00357 
00358      for( iv=0 ; iv < nv ; iv++ ){            /* loop over sub-bricks = columns */
00359        switch( DSET_BRICK_TYPE(dset,iv) ){
00360           default:
00361             val[0] = 0.0f ;
00362           break ;
00363 
00364           case MRI_float:{
00365             float *bar = DSET_ARRAY(dset,iv) ; val[0] = bar[ii] ;
00366           }
00367           break ;
00368 
00369           case MRI_short:{
00370             short *bar = DSET_ARRAY(dset,iv) ;
00371             float fac = DSET_BRICK_FACTOR(dset,iv) ; if( fac == 0.0 ) fac = 1.0 ;
00372             val[0] = fac*bar[ii] ;
00373           }
00374           break ;
00375 
00376           case MRI_byte:{
00377             byte *bar = DSET_ARRAY(dset,iv) ;
00378             float fac = DSET_BRICK_FACTOR(dset,iv) ; if( fac == 0.0 ) fac = 1.0 ;
00379             val[0] = fac*bar[ii] ;
00380           }
00381           break ;
00382 
00383           /* below here, we convert complicated types to floats, losing data! */
00384 
00385           case MRI_complex:{
00386             complex *bar = DSET_ARRAY(dset,iv) ;
00387             val[0] = CABS(bar[ii]) ;
00388           }
00389           break ;
00390 
00391           case MRI_rgb:{
00392             rgbyte *bar = DSET_ARRAY(dset,iv) ;
00393             val[0] = (0.299*bar[ii].r+0.587*bar[ii].g+0.114*bar[ii].g) ;
00394           }
00395           break ;
00396        } /* end of switch on sub-brick data type */
00397 
00398        if( binflag ) fwrite( val , sizeof(float) , 1 , fp ) ;
00399        else          fprintf( fp , " %g" , val[0] ) ;
00400 
00401      } /* end of loop over sub-bricks */
00402 
00403      if( !binflag ) fprintf(fp,"\n") ;
00404 
00405    } /* end of loop over voxels */
00406 
00407    /* NIML-style trailer */
00408 
00409    fflush(fp) ;
00410 
00411    if( fp != stdout ){
00412      fprintf(fp,"%c </AFNI_3D_dataset>\n",shp) ;
00413      fclose(fp) ;
00414    }
00415 
00416    EXRETURN ;
00417 }