00001 #include "mrilib.h"
00002 #include "mayo_analyze.h"
00003 #include "thd.h"
00004 
00005 /*******************************************************************/
00006 /******* 27 Aug 2002: Read an ANALYZE file as an AFNI dataset ******/
00007 /*******************************************************************/
00008 
00009 /*---------------------------------------------------------------*/
00010 /*! Swap the 2 bytes pointed to by ppp: ab -> ba. */
00011 
00012 static void swap_2(void *ppp)
00013 {
00014    unsigned char *pntr = (unsigned char *) ppp ;
00015    unsigned char b0, b1;
00016 
00017    b0 = *pntr; b1 = *(pntr+1);
00018    *pntr = b1; *(pntr+1) = b0;
00019 }
00020 
00021 /*---------------------------------------------------------------*/
00022 /*! Swap the 4 bytes pointed to by ppp: abcd -> dcba. */
00023 
00024 static void swap_4(void *ppp)
00025 {
00026    unsigned char *pntr = (unsigned char *) ppp ;
00027    unsigned char b0, b1, b2, b3;
00028 
00029    b0 = *pntr; b1 = *(pntr+1); b2 = *(pntr+2); b3 = *(pntr+3);
00030    *pntr = b3; *(pntr+1) = b2; *(pntr+2) = b1; *(pntr+3) = b0;
00031 }
00032 
00033 #if 0
00034 /*---------------------------------------------------------------*/
00035 /*! Swap the 8 bytes pointed to by ppp: abcdefgh -> hgfedcba. */
00036 
00037 static void swap_8(void *ppp)
00038 {
00039    unsigned char *pntr = (unsigned char *) ppp ;
00040    unsigned char b0, b1, b2, b3;
00041    unsigned char b4, b5, b6, b7;
00042 
00043    b0 = *pntr    ; b1 = *(pntr+1); b2 = *(pntr+2); b3 = *(pntr+3);
00044    b4 = *(pntr+4); b5 = *(pntr+5); b6 = *(pntr+6); b7 = *(pntr+7);
00045 
00046    *pntr     = b7; *(pntr+1) = b6; *(pntr+2) = b5; *(pntr+3) = b4;
00047    *(pntr+4) = b3; *(pntr+5) = b2; *(pntr+6) = b1; *(pntr+7) = b0;
00048 }
00049 #endif
00050 
00051 /*---------------------------------------------------------------*/
00052 /*! Byte swap ANALYZE file header in various places */
00053 
00054 static void swap_analyze_hdr( struct dsr *pntr )
00055 {
00056    swap_4(&pntr->hk.sizeof_hdr) ;
00057    swap_4(&pntr->hk.extents) ;
00058    swap_2(&pntr->hk.session_error) ;
00059    swap_2(&pntr->dime.dim[0]) ;
00060    swap_2(&pntr->dime.dim[1]) ;
00061    swap_2(&pntr->dime.dim[2]) ;
00062    swap_2(&pntr->dime.dim[3]) ;
00063    swap_2(&pntr->dime.dim[4]) ;
00064    swap_2(&pntr->dime.dim[5]) ;
00065    swap_2(&pntr->dime.dim[6]) ;
00066    swap_2(&pntr->dime.dim[7]) ;
00067 #if 0
00068    swap_2(&pntr->dime.unused1) ;
00069 #endif
00070    swap_2(&pntr->dime.datatype) ;
00071    swap_2(&pntr->dime.bitpix) ;
00072    swap_4(&pntr->dime.pixdim[0]) ;
00073    swap_4(&pntr->dime.pixdim[1]) ;
00074    swap_4(&pntr->dime.pixdim[2]) ;
00075    swap_4(&pntr->dime.pixdim[3]) ;
00076    swap_4(&pntr->dime.pixdim[4]) ;
00077    swap_4(&pntr->dime.pixdim[5]) ;
00078    swap_4(&pntr->dime.pixdim[6]) ;
00079    swap_4(&pntr->dime.pixdim[7]) ;
00080    swap_4(&pntr->dime.vox_offset) ;
00081    swap_4(&pntr->dime.funused1) ;
00082    swap_4(&pntr->dime.funused2) ;
00083    swap_4(&pntr->dime.cal_max) ;
00084    swap_4(&pntr->dime.cal_min) ;
00085    swap_4(&pntr->dime.compressed) ;
00086    swap_4(&pntr->dime.verified) ;
00087    swap_2(&pntr->dime.dim_un0) ;
00088    swap_4(&pntr->dime.glmax) ;
00089    swap_4(&pntr->dime.glmin) ;
00090 }
00091 
00092 /*-----------------------------------------------------------------*/
00093 #ifndef DONT_INCLUDE_ANALYZE_STRUCT
00094 #define DONT_INCLUDE_ANALYZE_STRUCT
00095 #endif
00096 #include "nifti1_io.h"                 /*** NIFTI spec and funcs ***/
00097 /*-----------------------------------------------------------------*/
00098 
00099 #undef swap_2
00100 #undef swap_4
00101 
00102 /*-----------------------------------------------------------------*/
00103 /*! Open an ANALYZE .hdr file as an unpopulated AFNI dataset.
00104     It will be populated from the .img file, in THD_load_analyze().
00105 -------------------------------------------------------------------*/
00106 
00107 THD_3dim_dataset * THD_open_analyze( char *hname )
00108 {
00109    FILE * fp ;
00110    char iname[THD_MAX_NAME] ;
00111    int ii , doswap ;
00112    struct dsr hdr ;    /* ANALYZE .hdr format */
00113    int ngood , length , datum_type , datum_len ;
00114    int   nx,ny,nz,nt ;
00115    float dx,dy,dz,dt ;
00116    float fac=0.0 ;     /* brick scaling factor */
00117 
00118    THD_3dim_dataset *dset=NULL ;
00119    char prefix[THD_MAX_PREFIX] , *ppp ;
00120    THD_ivec3 nxyz , orixyz ;
00121    THD_fvec3 dxyz , orgxyz ;
00122    int iview ;
00123    short spmorg , spmxx,spmyy,spmzz ;  /* 03 Nov 2003 */
00124 
00125 ENTRY("THD_open_analyze") ;
00126 
00127    /* 28 Aug 2003: check if this is a NIFTI file instead */
00128 
00129    { nifti_image *nim = nifti_image_read(hname,0) ;
00130      if( nim != NULL && nim->nifti_type > 0 ){
00131        nifti_image_free(nim); dset = THD_open_nifti(hname); RETURN(dset);
00132      }
00133    }
00134 
00135    /*-- check & prepare filenames --*/
00136 
00137    if( hname == NULL ) RETURN( NULL );
00138    ii = strlen(hname) ; if( ii < 5 ) RETURN( NULL );
00139    if( strcmp(hname+ii-4,".hdr") != 0 ) RETURN( NULL );
00140    strcpy(iname,hname) ; strcpy(iname+ii-3,"img") ;  /* .img filename */
00141 
00142    /*-- read .hdr file into struct --*/
00143 
00144    fp = fopen( hname , "rb" ) ;
00145    if( fp == NULL ) RETURN( NULL );
00146    hdr.dime.dim[0] = 0 ;
00147    fread( &hdr , 1 , sizeof(struct dsr) , fp ) ;
00148    fclose(fp) ;
00149    if( hdr.dime.dim[0] == 0 ){    /* bad input */
00150      fprintf(stderr,"*** ANALYZE file %s has dim[0]=0!\n",hname) ;
00151      RETURN( NULL ) ;
00152    }
00153 
00154    /*-- check .img file for existence and size --*/
00155 
00156    length = THD_filesize(iname) ;  /* will use this later */
00157    if( length <= 0 ){
00158      fprintf(stderr,"*** Can't find ANALYZE file %s\n",iname) ;
00159      RETURN( NULL );
00160    }
00161 
00162    /*-- check for swap-age of header --*/
00163 
00164    doswap = (hdr.dime.dim[0] < 0 || hdr.dime.dim[0] > 15) ;
00165    if( doswap ) swap_analyze_hdr( &hdr ) ;
00166 
00167    /*-- get a scale factor for images --*/
00168 
00169    if( !AFNI_noenv("AFNI_ANALYZE_SCALE") ){
00170      fac = hdr.dime.funused1 ;
00171      (void) thd_floatscan( 1 , &fac ) ;
00172      if( fac < 0.0 || fac == 1.0 ) fac = 0.0 ;
00173    }
00174 
00175    /*-- get data type for each voxel --*/
00176 
00177    switch( hdr.dime.datatype ){
00178      default:
00179         fprintf(stderr,"*** %s: Unsupported ANALYZE datatype=%d (%s)\n",
00180                 hname,hdr.dime.datatype,ANDT_string(hdr.dime.datatype) ) ;
00181      RETURN( NULL );
00182 
00183      case ANDT_UNSIGNED_CHAR: datum_type = MRI_byte   ; break;
00184      case ANDT_SIGNED_SHORT:  datum_type = MRI_short  ; break;
00185      case ANDT_FLOAT:         datum_type = MRI_float  ; break;
00186      case ANDT_COMPLEX:       datum_type = MRI_complex; break;
00187      case ANDT_RGB:           datum_type = MRI_rgb    ; break;
00188 #if 0
00189      case ANDT_SIGNED_INT:    datum_type = MRI_int    ; break; /* not in AFNI */
00190 #endif
00191    }
00192 
00193    datum_len = mri_datum_size(datum_type) ;  /* number bytes per voxel */
00194 
00195    /*-- compute dimensions of images, and number of images --*/
00196 
00197    nx = hdr.dime.dim[1] ;
00198    ny = hdr.dime.dim[2] ;
00199    if( nx < 2 || ny < 2 ) RETURN( NULL );
00200 
00201    switch( hdr.dime.dim[0] ){
00202      case 2:  nz = 1 ; nt = 1 ;                           ; break ;
00203      case 3:  nz = hdr.dime.dim[3] ; nt = 1 ;             ; break ;
00204      case 4:  nz = hdr.dime.dim[3] ; nt = hdr.dime.dim[4] ; break ;
00205 
00206      default:
00207        fprintf(stderr,"*** ANALYZE file %s has %d dimensions!\n",
00208                hname,hdr.dime.dim[0]) ;
00209        RETURN( NULL ) ;
00210    }
00211    if( nz < 1 ) nz = 1 ;
00212    if( nt < 1 ) nt = 1 ;
00213 
00214    /*-- voxel sizes --*/
00215 
00216    dx = fabs(hdr.dime.pixdim[1]) ; if( dx == 0.0 )            dx = 1.0 ;
00217    dy = fabs(hdr.dime.pixdim[2]) ; if( dy == 0.0 )            dy = 1.0 ;
00218    dz = fabs(hdr.dime.pixdim[3]) ; if( dz == 0.0 )            dz = 1.0 ;
00219    dt = fabs(hdr.dime.pixdim[4]) ; if( dt == 0.0 || nt == 1 ) dt = 1.0 ;
00220 
00221    ngood = datum_len*nx*ny*nz*nt ;  /* number bytes needed in .img file */
00222    if( length < ngood ){
00223      fprintf( stderr,
00224        "*** ANALYZE file %s is %d bytes long but must be %d bytes long\n"
00225        "*** for nx=%d ny=%d nz=%d nt=%d and %d bytes/voxel\n",
00226        iname,length,ngood,nx,ny,nz,nt,datum_len ) ;
00227      fclose(fp) ; RETURN( NULL );
00228    }
00229 
00230    /* check SPM originator field - 03 Nov 2003 */
00231    /* 11 Mar 2004 - oops: SPM indexes start at 1 - RWCox */
00232 
00233    spmorg = spmxx = spmyy = spmzz = 0 ;
00234    { short xyzuv[5] , xx,yy,zz ;
00235      memcpy( xyzuv , hdr.hist.originator , 10 ) ;
00236      if( xyzuv[3] == 0 && xyzuv[4] == 0 ){
00237        xx = xyzuv[0] ; yy = xyzuv[1] ; zz = xyzuv[2] ;
00238        if( doswap ){ swap_2(&xx); swap_2(&yy); swap_2(&zz); }
00239        if( xx > 0 && xx < nx-1 &&
00240            yy > 0 && yy < ny-1 &&
00241            zz > 0 && zz < nz-1   ){ spmorg=1; spmxx=xx-1; spmyy=yy-1; spmzz=zz-1; }
00242      }
00243    }
00244 
00245    /*-- make a dataset --*/
00246 
00247    dset = EDIT_empty_copy(NULL) ;
00248 
00249    dset->idcode.str[0] = 'A' ;  /* overwrite 1st 3 bytes */
00250    dset->idcode.str[1] = 'N' ;
00251    dset->idcode.str[2] = 'Z' ;
00252 
00253    MCW_hash_idcode( hname , dset ) ;  /* 06 May 2005 */
00254 
00255    ppp = THD_trailname(hname,0) ;                   /* strip directory */
00256    MCW_strncpy( prefix , ppp , THD_MAX_PREFIX ) ;   /* to make prefix */
00257 
00258    nxyz.ijk[0] = nx ; dxyz.xyz[0] = dx ;  /* setup axes lengths and voxel sizes */
00259    nxyz.ijk[1] = ny ; dxyz.xyz[1] = dy ;
00260    nxyz.ijk[2] = nz ; dxyz.xyz[2] = dz ;
00261 
00262    /*-- default orientation is LPI, but user can alter via environment --*/
00263 
00264    { char *ori = getenv("AFNI_ANALYZE_ORIENT") ;
00265      int oxx,oyy,ozz ;
00266      if( ori == NULL || strlen(ori) < 3 ) ori = "LPI"; /* set default LPI */
00267 
00268      oxx = ORCODE(ori[0]); oyy = ORCODE(ori[1]); ozz = ORCODE(ori[2]);
00269      if( !OR3OK(oxx,oyy,ozz) ){
00270        oxx = ORI_L2R_TYPE; oyy = ORI_P2A_TYPE; ozz = ORI_I2S_TYPE; /* LPI? */
00271      }
00272 
00273      orixyz.ijk[0] = oxx ;
00274      orixyz.ijk[1] = oyy ;
00275      orixyz.ijk[2] = ozz ;
00276    }
00277 
00278    /*-- origin of coordinates --*/
00279 
00280    orgxyz.xyz[0] = 0.5*dx ;  /* FSL:  (0,0,0) is at outer corner of 1st voxel */
00281    orgxyz.xyz[1] = 0.5*dy ;  /* AFNI: these are coords of center of 1st voxel */
00282    orgxyz.xyz[2] = 0.5*dz ;
00283 
00284    /*-- 04 Oct 2002: allow auto-centering of ANALYZE datasets --*/
00285 
00286    if( AFNI_yesenv("AFNI_ANALYZE_AUTOCENTER") ){
00287      orgxyz.xyz[0] = -0.5 * (nx-1) * dx ;
00288      orgxyz.xyz[1] = -0.5 * (ny-1) * dy ;
00289      orgxyz.xyz[2] = -0.5 * (nz-1) * dz ;
00290    }
00291 
00292    
00293    iview = VIEW_ORIGINAL_TYPE ;   /* can't tell if it is Talairach-ed (default)*/
00294    {/* ZSS Dec 15 03 */
00295       char *vie = getenv("AFNI_ANALYZE_VIEW") ;
00296       if (!vie) iview = VIEW_ORIGINAL_TYPE; 
00297       else {
00298          if (strcmp(vie, "tlrc") == 0) iview = VIEW_TALAIRACH_TYPE; 
00299          else if (strcmp(vie, "orig") == 0) iview = VIEW_ORIGINAL_TYPE;
00300          else {
00301             fprintf (stderr,  "\nWarning: Bad value (%s) for environment \n"
00302                               "variable AFNI_ANALYZE_VIEW. Choose from:\n"
00303                               "orig or tlrc.\n"
00304                               "Assuming orig view.\n", vie);
00305             iview = VIEW_ORIGINAL_TYPE;  
00306          }
00307       }
00308    }
00309    
00310    if( AFNI_yesenv("AFNI_ANALYZE_ORIGINATOR") && spmorg ){  /* 03 Nov 2003 */
00311      if ( !getenv ("AFNI_ANALYZE_VIEW") ) { /* ZSS Dec. 16 03 */
00312        iview         = VIEW_TALAIRACH_TYPE ;     /* for backward compatibility */
00313      } 
00314      orgxyz.xyz[0] = -spmxx * dx ; /* (0,0,0) is at (spmxx,spmyy,spmzz) */
00315      orgxyz.xyz[1] = -spmyy * dy ;
00316      orgxyz.xyz[2] = -spmzz * dz ;
00317    }
00318 
00319    /* 10 Oct 2002: change voxel size signs, if axis orientation is negative */
00320    /*              [above, we assumed that axes were oriented in - to + way] */
00321 
00322    if( ORIENT_sign[orixyz.ijk[0]] == '-' ){
00323      dxyz.xyz[0]   = -dxyz.xyz[0]   ;
00324      orgxyz.xyz[0] = -orgxyz.xyz[0] ;
00325    }
00326 
00327    if( ORIENT_sign[orixyz.ijk[1]] == '-' ){
00328      dxyz.xyz[1]   = -dxyz.xyz[1]   ;
00329      orgxyz.xyz[1] = -orgxyz.xyz[1] ;
00330    }
00331 
00332    if( ORIENT_sign[orixyz.ijk[2]] == '-' ){
00333      dxyz.xyz[2]   = -dxyz.xyz[2]   ;
00334      orgxyz.xyz[2] = -orgxyz.xyz[2] ;
00335    }
00336 
00337    /*-- actually send the values above into the dataset header --*/
00338 
00339    EDIT_dset_items( dset ,
00340                       ADN_prefix      , prefix ,
00341                       ADN_datum_all   , datum_type ,
00342                       ADN_nxyz        , nxyz ,
00343                       ADN_xyzdel      , dxyz ,
00344                       ADN_xyzorg      , orgxyz ,
00345                       ADN_xyzorient   , orixyz ,
00346                       ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00347                       ADN_nvals       , nt ,
00348                       ADN_type        , HEAD_ANAT_TYPE ,
00349                       ADN_view_type   , iview ,
00350                       ADN_func_type   , ANAT_MRAN_TYPE ,
00351                     ADN_none ) ;
00352 
00353    if( nt > 1 )              /** pretend it is 3D+time **/
00354       EDIT_dset_items( dset ,
00355                          ADN_func_type, ANAT_EPI_TYPE ,
00356                          ADN_ntt      , nt ,
00357                          ADN_ttorg    , 0.0 ,
00358                          ADN_ttdel    , dt ,
00359                          ADN_ttdur    , 0.0 ,
00360                          ADN_tunits   , UNITS_SEC_TYPE ,
00361                        ADN_none ) ;
00362 
00363    /* make it a func if is is a FSL special name */
00364 
00365 #ifdef ALLOW_FSL_FEAT
00366    if( strstr(prefix,"stat") != NULL )
00367      EDIT_dset_items( dset ,
00368                         ADN_type      , HEAD_FUNC_TYPE ,
00369                         ADN_func_type , FUNC_FIM_TYPE ,
00370                       ADN_none ) ;
00371 #endif
00372 
00373    /*-- set brick factors (all the same) --*/
00374 
00375    if( fac > 0.0 ){
00376      float *bf = AFMALL(float, sizeof(float)*nt) ;
00377      for( ii=0 ; ii < nt ; ii++ ) bf[ii] = fac ;
00378      EDIT_dset_items( dset , ADN_brick_fac,bf , ADN_none ) ;
00379      free(bf) ;
00380    }
00381 
00382    /*-- set byte order (for reading from disk) --*/
00383 
00384    ii = mri_short_order() ;
00385    if( doswap ) ii = REVERSE_ORDER(ii) ;
00386    dset->dblk->diskptr->byte_order = ii ;
00387 
00388    /*-- flag to read data from disk using ANALYZE mode --*/
00389 
00390    dset->dblk->diskptr->storage_mode = STORAGE_BY_ANALYZE ;
00391    strcpy( dset->dblk->diskptr->brick_name , iname ) ;
00392 
00393    RETURN(dset) ;
00394 }
00395 
00396 /*---------------------------------------------------------------------*/
00397 /*! Load an ANALYZE dataset from disk.
00398 -----------------------------------------------------------------------*/
00399 
00400 void THD_load_analyze( THD_datablock *dblk )
00401 {
00402    THD_diskptr *dkptr ;
00403    int nx,ny,nz,nv , nxy,nxyz,nxyzv , ibr,nbad ;
00404    FILE *fp ;
00405    void *ptr ;
00406 
00407 ENTRY("THD_load_analyze") ;
00408 
00409    /*-- check inputs --*/
00410 
00411    if( !ISVALID_DATABLOCK(dblk)                          ||
00412        dblk->diskptr->storage_mode != STORAGE_BY_ANALYZE ||
00413        dblk->brick == NULL                                 ) EXRETURN ;
00414 
00415    dkptr = dblk->diskptr ;
00416 
00417    fp = fopen( dkptr->brick_name , "rb" ) ;  /* .img file */
00418    if( fp == NULL ) EXRETURN ;
00419 
00420    /*-- allocate space for data --*/
00421 
00422    nx = dkptr->dimsizes[0] ;
00423    ny = dkptr->dimsizes[1] ;  nxy   = nx * ny   ;
00424    nz = dkptr->dimsizes[2] ;  nxyz  = nxy * nz  ;
00425    nv = dkptr->nvals       ;  nxyzv = nxyz * nv ;
00426 
00427    dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
00428 
00429    /*-- malloc space for each brick separately --*/
00430 
00431    for( nbad=ibr=0 ; ibr < nv ; ibr++ ){
00432       if( DBLK_ARRAY(dblk,ibr) == NULL ){
00433          ptr = AFMALL(void, DBLK_BRICK_BYTES(dblk,ibr) ) ;
00434          mri_fix_data_pointer( ptr ,  DBLK_BRICK(dblk,ibr) ) ;
00435          if( ptr == NULL ) nbad++ ;
00436       }
00437    }
00438 
00439    /*-- if couldn't get them all, take our ball and go home in a snit --*/
00440 
00441    if( nbad > 0 ){
00442       fprintf(stderr,
00443               "\n** failed to malloc %d ANALYZE bricks out of %d\n\a",nbad,nv);
00444       for( ibr=0 ; ibr < nv ; ibr++ ){
00445         if( DBLK_ARRAY(dblk,ibr) != NULL ){
00446            free(DBLK_ARRAY(dblk,ibr)) ;
00447            mri_fix_data_pointer( NULL , DBLK_BRICK(dblk,ibr) ) ;
00448         }
00449       }
00450       fclose(fp) ; EXRETURN ;
00451    }
00452 
00453    /*-- read data from .img file into sub-brick arrays! --*/
00454 
00455    for( ibr=0 ; ibr < nv ; ibr++ )
00456      fread( DBLK_ARRAY(dblk,ibr), 1, DBLK_BRICK_BYTES(dblk,ibr), fp ) ;
00457 
00458    fclose(fp) ;
00459 
00460    /*-- swap bytes? --*/
00461 
00462    if( dkptr->byte_order != mri_short_order() ){
00463      for( ibr=0 ; ibr < nv ; ibr++ ){
00464        switch( DBLK_BRICK_TYPE(dblk,ibr) ){
00465          case MRI_short:
00466            mri_swap2( DBLK_BRICK_NVOX(dblk,ibr) , DBLK_ARRAY(dblk,ibr) ) ;
00467          break ;
00468 
00469          case MRI_complex:
00470            mri_swap4( 2*DBLK_BRICK_NVOX(dblk,ibr), DBLK_ARRAY(dblk,ibr)) ;
00471          break ;
00472 
00473          case MRI_float:
00474          case MRI_int:
00475            mri_swap4( DBLK_BRICK_NVOX(dblk,ibr) , DBLK_ARRAY(dblk,ibr) ) ;
00476          break ;
00477        }
00478      }
00479    }
00480 
00481    /*-- check floats --*/
00482 
00483    for( ibr=0 ; ibr < nv ; ibr++ ){
00484      if( DBLK_BRICK_TYPE(dblk,ibr) == MRI_float ){
00485        nbad += thd_floatscan( DBLK_BRICK_NVOX(dblk,ibr) ,
00486                               DBLK_ARRAY(dblk,ibr)        ) ;
00487 
00488      } else if( DBLK_BRICK_TYPE(dblk,ibr) == MRI_complex ) {
00489        nbad += thd_complexscan( DBLK_BRICK_NVOX(dblk,ibr) ,
00490                                 DBLK_ARRAY(dblk,ibr)        ) ;
00491      }
00492    }
00493    if( nbad > 0 )
00494      fprintf(stderr ,
00495              "*** %s: found %d float errors -- see program float_scan\n" ,
00496              dkptr->brick_name , nbad ) ;
00497 
00498    EXRETURN ;
00499 }