00001 #include "mrilib.h"
00002 #include "thd.h"
00003 
00004 /*******************************************************************/
00005 /********** 21 Mar 2003: Read a 3D file as an AFNI dataset *********/
00006 /*******************************************************************/
00007 
00008 /*-----------------------------------------------------------------*/
00009 /*! Open a NIML 3D file as an AFNI dataset.  Each column is a
00010     separate sub-brick.  Dataset returned is empty (no data).
00011 -------------------------------------------------------------------*/
00012 
00013 THD_3dim_dataset * THD_open_3D( char *pathname )
00014 {
00015    THD_3dim_dataset *dset=NULL ;
00016    char prefix[1024] , *ppp ;
00017    THD_ivec3 nxyz , orixyz ;
00018    THD_fvec3 dxyz , orgxyz ;
00019    int nvals , ii ;
00020    NI_element *nel ;
00021    NI_stream ns ;
00022 
00023 ENTRY("THD_open_3D") ;
00024 
00025    /*-- read input file into NI_element --*/
00026 
00027    if( pathname == NULL || *pathname == '\0' ) RETURN(NULL) ;
00028 
00029    ppp = (char*)calloc( sizeof(char) , strlen(pathname)+16 ) ;
00030 
00031    strcpy(ppp,"file:") ; strcat(ppp,pathname) ;
00032    ns = NI_stream_open( ppp , "r" ) ; free(ppp) ;
00033    if( ns == NULL ) RETURN(NULL) ;
00034 
00035 STATUS("reading header") ;
00036 
00037    NI_skip_procins(1) ; NI_read_header_only(1) ;
00038    nel = NI_read_element(ns,333); NI_stream_close(ns);
00039    NI_skip_procins(0) ; NI_read_header_only(0) ;
00040 
00041    /*-- check data element for reasonability --*/
00042 
00043    if( nel == NULL                               ||   /* bad read */
00044        NI_element_type(nel) != NI_ELEMENT_TYPE   ||   /* bad element */
00045        nel->vec_num <= 0                         ||   /* no data */
00046        nel->vec_len <= 0                         ||   /* no data */
00047        strcmp(nel->name,"AFNI_3D_dataset") != 0  ||   /* incorrect data */
00048        nel->vec_rank > 3                           ){ /* weird header */
00049 
00050      fprintf(stderr,"** Can't read 3D head from %s\n",pathname) ;
00051      NI_free_element(nel) ; RETURN(NULL) ;
00052    }
00053 
00054 STATUS("checking header") ;
00055 
00056    /*-- check column types to make sure they are all numeric --*/
00057    /*   [AFNI doesn't like String or compound type datasets]   */
00058 
00059    for( ii=0 ; ii < nel->vec_num ; ii++ ){
00060 
00061      if( !NI_IS_NUMERIC_TYPE(nel->vec_typ[ii]) ){
00062        fprintf(stderr,"** 3D file %s isn't numeric!\n",pathname) ;
00063        NI_free_element(nel) ; RETURN(NULL) ;
00064      }
00065    }
00066 
00067    /*-- now have good data element ==> make a dataset --*/
00068 
00069 STATUS("making dataset") ;
00070 
00071    dset = EDIT_empty_copy(NULL) ;  /* default dataset */
00072 
00073    /* set prefix from input filename */
00074 
00075 STATUS("setting prefix") ;
00076 
00077    ppp = THD_trailname(pathname,0) ;              /* strip directory */
00078    NI_strncpy( prefix , ppp , THD_MAX_PREFIX ) ;  /* to make prefix */
00079 
00080    /* set grid sizes from element header */
00081 
00082 STATUS("setting grid sizes") ;
00083 
00084    nxyz.ijk[0] = nel->vec_len ; nxyz.ijk[1] = nxyz.ijk[2] = 1 ;
00085    if( nel->vec_axis_len != NULL ){
00086      if( nel->vec_rank >= 1 ) nxyz.ijk[0] = nel->vec_axis_len[0] ;
00087      if( nel->vec_rank >= 2 ) nxyz.ijk[1] = nel->vec_axis_len[1] ;
00088      if( nel->vec_rank >= 3 ) nxyz.ijk[2] = nel->vec_axis_len[2] ;
00089    }
00090 
00091    /* set grid spacings */
00092 
00093 STATUS("setting grid spacings") ;
00094 
00095    dxyz.xyz[0] = dxyz.xyz[1] = dxyz.xyz[2] = 1.0 ;
00096    if( nel->vec_axis_delta != NULL ){
00097      if( nel->vec_rank >= 1) dxyz.xyz[0] = nel->vec_axis_delta[0] ;
00098      if( nel->vec_rank >= 2) dxyz.xyz[1] = nel->vec_axis_delta[1] ;
00099      if( nel->vec_rank >= 3) dxyz.xyz[2] = nel->vec_axis_delta[2] ;
00100    }
00101 
00102    /* set grid origins */
00103 
00104 STATUS("setting grid origins") ;
00105 
00106    orgxyz.xyz[0] = orgxyz.xyz[1] = orgxyz.xyz[2] = 0.0 ;
00107    if( nel->vec_axis_origin != NULL ){
00108      if( nel->vec_rank >= 1) orgxyz.xyz[0] = nel->vec_axis_origin[0] ;
00109      if( nel->vec_rank >= 2) orgxyz.xyz[1] = nel->vec_axis_origin[1] ;
00110      if( nel->vec_rank >= 3) orgxyz.xyz[2] = nel->vec_axis_origin[2] ;
00111    }
00112 
00113    /* set grid orientations (default is RAI) */
00114 
00115 STATUS("setting grid orientation") ;
00116 
00117    { char orcx='R', orcy='A', orcz='I' ;
00118      int oxx,oyy,ozz ;
00119      if( nel->vec_rank == 3 && nel->vec_axis_label != NULL ){
00120        orcx = toupper(nel->vec_axis_label[0][0]) ;
00121        orcy = toupper(nel->vec_axis_label[1][0]) ;
00122        orcz = toupper(nel->vec_axis_label[2][0]) ;
00123      }
00124      oxx = ORCODE(orcx); oyy = ORCODE(orcy); ozz = ORCODE(orcz);
00125      if( !OR3OK(oxx,oyy,ozz) ){
00126        oxx = ORI_R2L_TYPE; oyy = ORI_A2P_TYPE; ozz = ORI_I2S_TYPE;
00127      }
00128 
00129      orixyz.ijk[0] = oxx ; orixyz.ijk[1] = oyy ; orixyz.ijk[2] = ozz ;
00130    }
00131 
00132    /* number of sub-bricks (one per vector/column) */
00133 
00134    nvals = nel->vec_num ;
00135 
00136    /* set idcode from element, or take random default one */
00137 
00138 STATUS("setting idcode") ;
00139 
00140    ppp = NI_get_attribute( nel , "self_idcode" ) ;
00141    if( ppp == NULL )
00142      ppp = NI_get_attribute( nel , "ni_idcode" ) ;
00143    if( ppp != NULL && *ppp != '\0' ){
00144      NI_strncpy( dset->idcode.str , ppp , MCW_IDSIZE ) ;
00145    } else {
00146      MCW_hash_idcode( pathname , dset ) ; /* 06 May 2005 */
00147      dset->idcode.str[0] = 'A' ;          /* overwrite 1st 3 bytes of idcode */
00148      dset->idcode.str[1] = '3' ;
00149      dset->idcode.str[2] = 'D' ;
00150    }
00151 
00152    /*-- now modify the default dataset --*/
00153 
00154 STATUS("Editing dataset") ;
00155 
00156    EDIT_dset_items( dset ,
00157                       ADN_prefix      , prefix ,
00158                       ADN_datum_array , nel->vec_typ ,
00159                       ADN_nxyz        , nxyz ,
00160                       ADN_xyzdel      , dxyz ,
00161                       ADN_xyzorg      , orgxyz ,
00162                       ADN_xyzorient   , orixyz ,
00163                       ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00164                       ADN_nvals       , nvals ,
00165                       ADN_type        , HEAD_ANAT_TYPE ,
00166                       ADN_func_type   , ANAT_BUCK_TYPE ,
00167                     ADN_none ) ;
00168 
00169    dset->dblk->diskptr->storage_mode = STORAGE_BY_3D ;
00170    NI_strncpy( dset->dblk->diskptr->brick_name , pathname , THD_MAX_NAME ) ;
00171 
00172    /*-- time axis? [03 Jun 2005] --*/
00173 
00174    ppp = NI_get_attribute( nel , "ni_timestep" ) ;
00175    if( ppp != NULL && nvals > 1 ){
00176      float dt = strtod(ppp,NULL) ; if( dt <= 0.0 ) dt = 1.0 ;
00177      EDIT_dset_items( dset ,
00178                         ADN_func_type , ANAT_EPI_TYPE ,
00179                         ADN_ntt       , nvals ,
00180                         ADN_ttdel     , dt ,
00181                         ADN_tunits    , UNITS_SEC_TYPE ,
00182                       ADN_none ) ;
00183    }
00184 
00185 STATUS("checking for statistics") ;
00186 
00187    /*-- see if we have any statistics bricks --*/
00188 
00189    ppp = NI_get_attribute( nel , "ni_stat" ) ;
00190    if( ppp != NULL ){
00191      NI_str_array *sar = NI_decode_string_list( ppp , ";" ) ;
00192      if( sar != NULL ){
00193        int itop=MIN(sar->num,nvals) , jj,ll ;
00194        char *dnam , qnam[64] ;
00195        for( ii=0 ; ii < itop ; ii++ ){
00196          if( strcmp(sar->str[ii],"none") == 0 ) continue ;
00197          for( jj=NI_STAT_FIRSTCODE ; jj <= NI_STAT_LASTCODE ; jj++ ){
00198            dnam = NI_stat_distname(jj) ;
00199            strcpy(qnam,dnam); strcat(qnam,"("); ll = strlen(qnam);
00200            if( strncmp(sar->str[ii],qnam,ll) == 0 ) break ;
00201          }
00202          if( jj >= AFNI_FIRST_STATCODE && jj <= AFNI_LAST_STATCODE ){
00203            float parm[4]={1,1,1,1} ; int np,kk,mm , sp ;
00204            np = NI_stat_numparam(jj) ; sp = ll ;
00205            for( kk=0 ; kk < np ; kk++ ){
00206              mm = 0 ; sscanf(sar->str[ii]+sp,"%f%n",parm+kk,&mm) ; sp += mm+1 ;
00207            }
00208            EDIT_STATAUX4( dset , ii , jj , parm[0],parm[1],parm[2],parm[3] ) ;
00209          }
00210        }
00211        NI_delete_str_array(sar) ;
00212      }
00213    }
00214 
00215    /*-- purge the NIML data element and return the new dataset --*/
00216 
00217 STATUS("freeing element") ;
00218 
00219    NI_free_element( nel ) ;
00220    RETURN(dset) ;
00221 }
00222 
00223 /*-----------------------------------------------------------------*/
00224 /*! Mark a dataset's storage mode.
00225 -------------------------------------------------------------------*/
00226 
00227 void THD_set_storage_mode( THD_3dim_dataset *dset , int mm )
00228 {
00229    if( !ISVALID_DSET(dset)         ||
00230        mm < 0                      ||
00231        mm > LAST_STORAGE_MODE      ||
00232        dset->dblk == NULL          ||
00233        dset->dblk->diskptr == NULL   ) return ;
00234 
00235    dset->dblk->diskptr->storage_mode = mm ;
00236 }
00237 
00238 /*-----------------------------------------------------------------*/
00239 /*!  Load a 3D dataset's data into memory.
00240      Called from THD_load_datablock() in thd_loaddblk.c.
00241 -------------------------------------------------------------------*/
00242 
00243 void THD_load_3D( THD_datablock *dblk )
00244 {
00245    THD_diskptr *dkptr ;
00246    NI_element *nel ;
00247    NI_stream ns ;
00248    int nxyz , iv,nv ;
00249    char *ppp ;
00250 
00251 ENTRY("THD_load_3D") ;
00252 
00253    /*-- open input [these errors should never occur] --*/
00254 
00255    if( !ISVALID_DATABLOCK(dblk)                     ||
00256        dblk->diskptr->storage_mode != STORAGE_BY_3D ||
00257        dblk->brick == NULL                            ) EXRETURN ;
00258 
00259    dkptr = dblk->diskptr ;
00260    nxyz  = dkptr->dimsizes[0] * dkptr->dimsizes[1] * dkptr->dimsizes[2] ;
00261    nv    = dkptr->nvals ;
00262 
00263    if( nxyz*nv > 1000000 ) fprintf(stderr,"++ Reading %s\n",dkptr->brick_name) ;
00264 
00265    ppp = (char*)calloc( sizeof(char) , strlen(dkptr->brick_name)+16 ) ;
00266 
00267    strcpy(ppp,"file:") ; strcat(ppp,dkptr->brick_name) ;
00268    ns = NI_stream_open( ppp , "r" ) ; free(ppp) ;
00269    if( ns == NULL ) EXRETURN ;
00270 
00271    NI_skip_procins(1) ;
00272    nel = NI_read_element(ns,333); NI_stream_close(ns);
00273    NI_skip_procins(0) ;
00274    if( nel == NULL ) EXRETURN ;
00275 
00276    /*-- allocate space for data --*/
00277 
00278    if( nxyz != nel->vec_len || nv != nel->vec_num ){
00279      fprintf(stderr,"THD_load_3D(%s): nxyz or nv mismatch!\n",dkptr->brick_name) ;
00280      fprintf(stderr,"                 expect nxyz=%d; got %d\n",nxyz, nel->vec_len);
00281      fprintf(stderr,"                 expect nv  =%d; got %d\n",nv  , nel->vec_num);
00282      NI_free_element(nel) ; EXRETURN ;
00283    }
00284 
00285    dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
00286 
00287    /** malloc space for each brick separately,
00288        and copy data from element into place  **/
00289 
00290    for( iv=0 ; iv < nv ; iv++ ){
00291      if( DBLK_ARRAY(dblk,iv) == NULL ){                    /* needs data */
00292        ppp = AFMALL(char, DBLK_BRICK_BYTES(dblk,iv) );     /* make space */
00293        if( ppp == NULL ) break ;                           /* bad bad bad */
00294        mri_fix_data_pointer( ppp, DBLK_BRICK(dblk,iv) ) ;
00295        memcpy( ppp, nel->vec[iv], DBLK_BRICK_BYTES(dblk,iv) ) ;
00296        NI_free(nel->vec[iv]) ; nel->vec[iv] = NULL ;
00297      }
00298    }
00299 
00300    NI_free_element(nel) ;
00301 
00302    /* if malloc failed, then delete all bricks */
00303 
00304    if( iv < nv ){
00305      fprintf(stderr,"\n** malloc failed for 3D dataset input!\n");
00306      for( iv=0 ; iv < nv ; iv++ ){
00307        if( DBLK_ARRAY(dblk,iv) != NULL ){
00308          free( DBLK_ARRAY(dblk,iv) ) ;
00309          mri_fix_data_pointer( NULL, DBLK_BRICK(dblk,iv) ) ;
00310        }
00311      }
00312    }
00313 
00314    EXRETURN ;
00315 }
00316 
00317 /*------------------------------------------------------------------*/
00318 /*! Write a dataset to disk as a 3D file.
00319     Called from THD_write_3dim_dataset().
00320     This is kind of cheating, since we just call THD_write_1D()
00321     instead.
00322 --------------------------------------------------------------------*/
00323 
00324 void THD_write_3D( char *sname, char *pname , THD_3dim_dataset *dset )
00325 {
00326    THD_write_1D( sname,pname,dset ) ;
00327 }