00001 #include "mrilib.h"
00002 
00003 #ifndef DONT_INCLUDE_ANALYZE_STRUCT
00004 #define DONT_INCLUDE_ANALYZE_STRUCT
00005 #endif
00006 #include "nifti1_io.h"   /** will include nifti1.h **/
00007 
00008 /*******************************************************************/
00009 /********** 26 Aug 2003: read a NIFTI-1 file as a dataset **********/
00010 /*******************************************************************/
00011 
00012 THD_3dim_dataset * THD_open_nifti( char *pathname )
00013 {
00014    THD_3dim_dataset *dset=NULL ;
00015    nifti_image *nim ;
00016    int ntt , nbuc , nvals ;
00017    int use_qform = 0 , use_sform = 0 ;
00018    int statcode = 0 , datum , iview , ibr ;
00019    THD_ivec3 orixyz , nxyz ;
00020    THD_fvec3 dxyz , orgxyz ;
00021    THD_mat33 R ;
00022    char *ppp , prefix[THD_MAX_PREFIX] ;
00023 
00024 ENTRY("THD_open_nifti") ;
00025 
00026    /*-- open input file --*/
00027 
00028    { /* set the nifti_io debug level       8 Apr 2005 [rickr] */
00029       char * ept = my_getenv("AFNI_NIFTI_DEBUG");
00030       if( ept != NULL ) nifti_set_debug_level(atoi(ept));
00031    }
00032 
00033    nim = nifti_image_read( pathname, 0 ) ;
00034 
00035    if( nim == NULL || nim->nifti_type == 0 ) RETURN(NULL) ;
00036 
00037    /*-- extract some useful AFNI-ish information from the nim struct --*/
00038 
00039    /* we must have at least 2 spatial dimensions */
00040 
00041    if( nim->nx < 2 || nim->ny < 2 ) RETURN(NULL) ;
00042 
00043    /* 4th dimension = time; 5th dimension = bucket:
00044       these are mutually exclusive in AFNI at present */
00045 
00046    ntt = nim->nt ; nbuc = nim->nu ;
00047    if( ntt > 1 && nbuc > 1 ){
00048      fprintf(stderr,
00049              "** AFNI can't deal with 5 dimensional NIfTI dataset file %s\n",
00050              pathname ) ;
00051      RETURN(NULL) ;
00052    }
00053    nvals = MAX(ntt,nbuc) ;
00054 
00055    /* determine type of dataset values:
00056       if we are scaling, or if the data type in the NIfTI file
00057       is something AFNI can't handle, then the result will be floats */
00058 
00059    switch( nim->datatype ){
00060      default:
00061        fprintf(stderr,
00062                "** AFNI can't handle NIFTI datatype=%d (%s) in file %s\n",
00063                nim->datatype, nifti_datatype_string(nim->datatype), pathname );
00064        RETURN(NULL) ;
00065      break ;
00066 
00067      case DT_UINT8:     datum = (nim->scl_slope != 0.0) ? MRI_float : MRI_byte  ; break ;
00068      case DT_INT16:     datum = (nim->scl_slope != 0.0) ? MRI_float : MRI_short ; break ;
00069 
00070      case DT_FLOAT32:   datum = MRI_float   ; break ;
00071      case DT_COMPLEX64: datum = MRI_complex ; break ;
00072      case DT_RGB24:     datum = MRI_rgb     ; break ;
00073 
00074      case DT_INT8:      /* NIfTI-1 data types that AFNI can't handle directly */
00075      case DT_UINT16:
00076      case DT_INT32:
00077      case DT_UINT32:
00078      case DT_FLOAT64:
00079        fprintf(stderr,
00080               "** AFNI converts NIFTI_datatype=%d (%s) in file %s to FLOAT32\n",
00081               nim->datatype, nifti_datatype_string(nim->datatype), pathname );
00082        datum = MRI_float ;
00083      break ;
00084 
00085 #if 0
00086      case DT_COMPLEX128:  /* this case would be too much like real work */
00087        fprintf(stderr,
00088                "** AFNI convert NIFTI_datatype=%d (%s) in file %s to COMPLEX64\n",
00089                nim->datatype, nifti_datatype_string(nim->datatype), pathname );
00090        datum = MRI_complex ;
00091      break ;
00092 #endif
00093    }
00094 
00095    /* check for statistics code */
00096 
00097    if( nim->intent_code >= NIFTI_FIRST_STATCODE &&
00098        nim->intent_code <= NIFTI_LAST_STATCODE    ){
00099 
00100      if( nim->intent_code > FUNC_PT_TYPE ){
00101        fprintf(stderr,
00102                "** AFNI doesn't understand NIFTI statistic type %d (%s) in file %s\n",
00103                nim->intent_code , nifti_intent_string(nim->intent_code) , pathname ) ;
00104      } else {
00105        statcode = nim->intent_code ;
00106        if( nbuc > 1 ){
00107          fprintf(stderr,
00108                  "** AFNI doesn't support NIFTI voxel-dependent statistic parameters"
00109                  " in file %s\n" , pathname ) ;
00110          statcode = 0 ;
00111        }
00112      }
00113    }
00114 
00115 
00116   /* KRH 07/11/05 -- adding ability to choose spatial transform
00117       from the options of qform, sform, bothform, or noform.
00118 
00119      If qform is present, it will be used.
00120 
00121      If qform is absent, but sform present, then the sform
00122        will be modified to be an orthogonal rotation and used.
00123 
00124      If both qform and sform are absent, then we will have
00125        an error.
00126 
00127      Previously assumed qform present.  */
00128 
00129    if ((nim->qform_code > 0) && (nim->sform_code > 0) ) {
00130      use_qform = 1 ;
00131      use_sform = 0 ;
00132    } else if (nim->qform_code > 0) {
00133      use_qform = 1 ;
00134      use_sform = 0 ;
00135    } else if (nim->sform_code > 0) { 
00136      use_qform = 0 ;
00137      use_sform = 1 ;
00138    } else {
00139      use_qform = 0 ;
00140      use_sform = 0 ;
00141      fprintf(stderr,
00142              "** NO spatial transform (neither preferred qform nor sform), in NIfTI dataset file:\n %s\nUsing senseless defaults.\n",
00143              pathname ) ;
00144    }
00145 
00146 
00147    if (use_qform) {
00148 
00149      float orgx, orgy, orgz ;
00150 
00151      /* determine orientation from the qto_xyz matrix,
00152       which transforms (i,j,k) voxel indexes to (x,y,z) LPI coordinates */
00153 
00154      LOAD_MAT(R, -nim->qto_xyz.m[0][0] ,  /* negate x and y   */
00155                  -nim->qto_xyz.m[0][1] ,  /* coefficients,    */
00156                  -nim->qto_xyz.m[0][2] ,  /* since AFNI works */
00157                  -nim->qto_xyz.m[1][0] ,  /* with RAI coords, */
00158                  -nim->qto_xyz.m[1][1] ,  /* but NIFTI uses   */
00159                  -nim->qto_xyz.m[1][2] ,  /* LPI coordinates. */
00160                   nim->qto_xyz.m[2][0] ,  /* [Which is my own] */
00161                   nim->qto_xyz.m[2][1] ,  /* [damn fault!!!!!] */
00162                   nim->qto_xyz.m[2][2]  ) ;
00163 
00164      orixyz = THD_matrix_to_orientation( R ) ;   /* compute orientation codes */
00165    
00166      iview = ((nim->qform_code == NIFTI_XFORM_TALAIRACH ) ||
00167               (nim->qform_code == NIFTI_XFORM_MNI_152 ))
00168              ? VIEW_TALAIRACH_TYPE : VIEW_ORIGINAL_TYPE ;
00169 
00170 
00171      /* load the offsets and the grid spacings */
00172 
00173      if (ORIENT_xyz[orixyz.ijk[0]] == 'z' )  {
00174        orgx = nim->qto_xyz.m[ORIENT_xyzint[orixyz.ijk[0]] - 1][3] ;
00175      } else {
00176        orgx = - nim->qto_xyz.m[ORIENT_xyzint[orixyz.ijk[0]] - 1][3] ;
00177      }
00178 
00179      if (ORIENT_xyz[orixyz.ijk[1]] == 'z' )  {
00180        orgy = nim->qto_xyz.m[ORIENT_xyzint[orixyz.ijk[1]] - 1][3] ;
00181      } else {
00182        orgy = - nim->qto_xyz.m[ORIENT_xyzint[orixyz.ijk[1]] - 1][3] ;
00183      }
00184 
00185      if (ORIENT_xyz[orixyz.ijk[2]] == 'z' )  {
00186        orgz = nim->qto_xyz.m[ORIENT_xyzint[orixyz.ijk[2]] - 1][3] ;
00187      } else {
00188        orgz = - nim->qto_xyz.m[ORIENT_xyzint[orixyz.ijk[2]] - 1][3] ;
00189      }
00190 
00191 
00192      LOAD_FVEC3( orgxyz ,  orgx ,
00193                            orgy ,
00194                            orgz  ) ;
00195 #if 0
00196      LOAD_FVEC3( orgxyz , -nim->qto_xyz.m[0][3] ,    /* again, negate  */
00197                         -nim->qto_xyz.m[1][3] ,    /* x and y coords */
00198                          nim->qto_xyz.m[2][3]  ) ;
00199 #endif
00200 
00201      /* AFNI space units are always mm */
00202 
00203      if( nim->xyz_units == NIFTI_UNITS_METER ){
00204        nim->dx *= 1000.0 ; nim->dy *= 1000.0 ; nim->dz *= 1000.0 ;
00205      } else if(  nim->xyz_units == NIFTI_UNITS_MICRON ){
00206        nim->dx *= 0.001  ; nim->dy *= 0.001  ; nim->dz *= 0.001  ;
00207      }
00208 
00209      LOAD_FVEC3( dxyz , (ORIENT_sign[orixyz.ijk[0]]=='+') ? nim->dx : -nim->dx ,
00210                         (ORIENT_sign[orixyz.ijk[1]]=='+') ? nim->dy : -nim->dy ,
00211                         (ORIENT_sign[orixyz.ijk[2]]=='+') ? nim->dz : -nim->dz  ) ;
00212 
00213    } else if (use_sform) {
00214 
00215      int orimap[7] = { 6 , 1 , 0 , 2 , 3 , 4 , 5 } ;
00216      int oritmp[3] ;
00217      float dxtmp, dytmp, dztmp ;
00218      float xmax, ymax, zmax ;
00219      float orgx, orgy, orgz ;
00220      float fig_merit, ang_merit ;
00221 
00222      /* convert sform to nifti orientation codes */
00223 
00224      nifti_mat44_to_orientation(nim->sto_xyz, &oritmp[0], &oritmp[1], &oritmp[2] ) ;
00225 
00226      /* convert nifti orientation codes to AFNI codes and store in vector */
00227 
00228      LOAD_IVEC3( orixyz , orimap[oritmp[0]] ,
00229                           orimap[oritmp[1]] ,
00230                           orimap[oritmp[2]] ) ;
00231 
00232      /* assume original view if there's no talairach id present */
00233 
00234      iview = ((nim->sform_code == NIFTI_XFORM_TALAIRACH ) ||
00235               (nim->sform_code == NIFTI_XFORM_MNI_152 ))
00236              ? VIEW_TALAIRACH_TYPE : VIEW_ORIGINAL_TYPE ;
00237 
00238      /* load the offsets and the grid spacings */
00239 
00240      if (ORIENT_xyz[orixyz.ijk[0]] == 'z' )  {
00241        orgx = nim->sto_xyz.m[ORIENT_xyzint[orixyz.ijk[0]] - 1][3] ;
00242      } else {
00243        orgx = - nim->sto_xyz.m[ORIENT_xyzint[orixyz.ijk[0]] - 1][3] ;
00244      }
00245 
00246      if (ORIENT_xyz[orixyz.ijk[1]] == 'z' )  {
00247        orgy = nim->sto_xyz.m[ORIENT_xyzint[orixyz.ijk[1]] - 1][3] ;
00248      } else {
00249        orgy = - nim->sto_xyz.m[ORIENT_xyzint[orixyz.ijk[1]] - 1][3] ;
00250      }
00251 
00252      if (ORIENT_xyz[orixyz.ijk[2]] == 'z' )  {
00253        orgz = nim->sto_xyz.m[ORIENT_xyzint[orixyz.ijk[2]] - 1][3] ;
00254      } else {
00255        orgz = - nim->sto_xyz.m[ORIENT_xyzint[orixyz.ijk[2]] - 1][3] ;
00256      }
00257 
00258 
00259      LOAD_FVEC3( orgxyz ,  orgx ,
00260                            orgy ,
00261                            orgz  ) ;
00262 
00263 #if 0
00264      LOAD_FVEC3( orgxyz , -nim->sto_xyz.m[0][3] ,    /* again, negate  */
00265                           -nim->sto_xyz.m[1][3] ,    /* x and y coords */
00266                            nim->sto_xyz.m[2][3] ) ;
00267 #endif
00268 
00269 #define MAXNUM(a,b) ( (a) > (b) ? (a):(b))
00270 #define MAX3(a,b,c) ( (MAXNUM(a,b)) > (MAXNUM(a,c)) ? (MAXNUM(a,b)):(MAXNUM(a,c)))
00271 #define MINNUM(a,b) ( (a) < (b) ? (a):(b))
00272 #define MIN3(a,b,c) ( (MINNUM(a,b)) < (MINNUM(a,c)) ? (MINNUM(a,b)):(MINNUM(a,c)))
00273 
00274      dxtmp = sqrt ( nim->sto_xyz.m[0][0] * nim->sto_xyz.m[0][0] +
00275                     nim->sto_xyz.m[1][0] * nim->sto_xyz.m[1][0] +
00276                     nim->sto_xyz.m[2][0] * nim->sto_xyz.m[2][0] ) ;
00277 
00278      xmax = MAX3(fabs(nim->sto_xyz.m[0][0]),fabs(nim->sto_xyz.m[1][0]),fabs(nim->sto_xyz.m[2][0])) / dxtmp ;
00279  
00280      dytmp = sqrt ( nim->sto_xyz.m[0][1] * nim->sto_xyz.m[0][1] +
00281                     nim->sto_xyz.m[1][1] * nim->sto_xyz.m[1][1] +
00282                     nim->sto_xyz.m[2][1] * nim->sto_xyz.m[2][1] ) ;
00283  
00284      ymax = MAX3(fabs(nim->sto_xyz.m[0][1]),fabs(nim->sto_xyz.m[1][1]),fabs(nim->sto_xyz.m[2][1])) / dytmp ;
00285 
00286      dztmp = sqrt ( nim->sto_xyz.m[0][2] * nim->sto_xyz.m[0][2] +
00287                     nim->sto_xyz.m[1][2] * nim->sto_xyz.m[1][2] +
00288                     nim->sto_xyz.m[2][2] * nim->sto_xyz.m[2][2] ) ;
00289 
00290      zmax = MAX3(fabs(nim->sto_xyz.m[0][2]),fabs(nim->sto_xyz.m[1][2]),fabs(nim->sto_xyz.m[2][2])) / dztmp ;
00291 
00292      fig_merit = MIN3(xmax,ymax,zmax) ;
00293      ang_merit = acos (fig_merit) * 180.0 / 3.141592653 ;
00294 
00295      if (fabs(ang_merit) > .01) {
00296        fprintf(stderr, "qform not present, sform used.\n"
00297                        "sform was not exact, and the worst axis is\n"
00298                        "%f degrees from plumb.\n",ang_merit ) ;
00299      }
00300 
00301      if( nim->xyz_units == NIFTI_UNITS_METER ){
00302        dxtmp *= 1000.0 ; dytmp *= 1000.0 ; dztmp *= 1000.0 ;
00303      } else if(  nim->xyz_units == NIFTI_UNITS_MICRON ){
00304        dxtmp *= 0.001  ; dytmp *= 0.001  ; dztmp *= 0.001  ;
00305      }
00306 
00307      LOAD_FVEC3( dxyz , (ORIENT_sign[orixyz.ijk[0]]=='+') ? dxtmp : -dxtmp ,
00308                         (ORIENT_sign[orixyz.ijk[1]]=='+') ? dytmp : -dytmp ,
00309                         (ORIENT_sign[orixyz.ijk[2]]=='+') ? dztmp : -dztmp ) ;
00310 
00311    } else { /* NO SPATIAL XFORM. BAD BAD BAD BAD BAD BAD. */
00312 
00313      float dxtmp, dytmp, dztmp ;
00314 
00315      /* if pixdim data are present, use them in order to set pixel 
00316         dimensions.  otherwise, set the dimensions to 1 unit.         */ 
00317 
00318      dxtmp = ((nim->pixdim[1] > 0) ? nim->pixdim[1] : 1) ;
00319      dytmp = ((nim->pixdim[2] > 0) ? nim->pixdim[2] : 1) ;
00320      dztmp = ((nim->pixdim[3] > 0) ? nim->pixdim[3] : 1) ;
00321 
00322      if( nim->xyz_units == NIFTI_UNITS_METER ){
00323        dxtmp *= 1000.0 ; dytmp *= 1000.0 ; dztmp *= 1000.0 ;
00324      } else if(  nim->xyz_units == NIFTI_UNITS_MICRON ){
00325        dxtmp *= 0.001  ; dytmp *= 0.001  ; dztmp *= 0.001  ;
00326      }
00327 
00328      /* set orientation to LPI by default    */
00329 
00330      LOAD_IVEC3( orixyz , 1 ,
00331                           2 ,
00332                           4 ) ;
00333 
00334      LOAD_FVEC3( dxyz , (ORIENT_sign[orixyz.ijk[0]]=='+') ? dxtmp : -dxtmp ,
00335                         (ORIENT_sign[orixyz.ijk[1]]=='+') ? dytmp : -dytmp ,
00336                         (ORIENT_sign[orixyz.ijk[2]]=='+') ? dztmp : -dztmp ) ;
00337 
00338      iview = VIEW_ORIGINAL_TYPE ;
00339 
00340      /* set origin to 0,0,0   */
00341 
00342      LOAD_FVEC3( orgxyz , 0 ,
00343                           0 ,
00344                           0 ) ;
00345    }
00346 
00347 
00348    /*-- make an AFNI dataset! --*/
00349 
00350    dset = EDIT_empty_copy(NULL) ;
00351 
00352    ppp  = THD_trailname(pathname,0) ;               /* strip directory */
00353    MCW_strncpy( prefix , ppp , THD_MAX_PREFIX ) ;   /* to make prefix */
00354 
00355    nxyz.ijk[0] = nim->nx ;                          /* grid dimensions */
00356    nxyz.ijk[1] = nim->ny ;
00357    nxyz.ijk[2] = nim->nz ;
00358 
00359    dset->idcode.str[0] = 'N' ;  /* overwrite 1st 3 bytes with something special */
00360    dset->idcode.str[1] = 'I' ;
00361    dset->idcode.str[2] = 'I' ;
00362 
00363    MCW_hash_idcode( pathname , dset ) ;  /* 06 May 2005 */
00364 
00365    EDIT_dset_items( dset ,
00366                       ADN_prefix      , prefix ,
00367                       ADN_datum_all   , datum ,
00368                       ADN_nxyz        , nxyz ,
00369                       ADN_xyzdel      , dxyz ,
00370                       ADN_xyzorg      , orgxyz ,
00371                       ADN_xyzorient   , orixyz ,
00372                       ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00373                       ADN_view_type   , iview ,
00374                       ADN_type        , (statcode != 0) ? HEAD_FUNC_TYPE
00375                                                         : HEAD_ANAT_TYPE ,
00376                     ADN_none ) ;
00377 
00378    /* not a time dependent dataset */
00379 
00380    if( ntt < 2 ){
00381      EDIT_dset_items( dset ,
00382                         ADN_nvals     , nbuc ,
00383                         ADN_datum_all , datum ,
00384                         ADN_func_type , (statcode != 0) ? FUNC_BUCK_TYPE
00385                                                         : ANAT_BUCK_TYPE ,
00386                       ADN_none ) ;
00387 
00388    } else {  /* is a time dependent dataset */
00389 
00390           if( nim->time_units == NIFTI_UNITS_MSEC ) nim->dt *= 0.001 ;
00391      else if( nim->time_units == NIFTI_UNITS_USEC ) nim->dt *= 1.e-6 ;
00392      EDIT_dset_items( dset ,
00393                         ADN_nvals     , ntt ,
00394                         ADN_ntt       , ntt ,
00395                         ADN_datum_all , datum ,
00396                         ADN_ttorg     , 0.0 ,
00397                         ADN_ttdel     , nim->dt ,
00398                         ADN_ttdur     , 0.0 ,
00399                         ADN_tunits    , UNITS_SEC_TYPE ,
00400                         ADN_func_type , (statcode != 0) ? FUNC_FIM_TYPE
00401                                                         : ANAT_EPI_TYPE ,
00402                       ADN_none ) ;
00403 
00404      /* if present, add stuff about the slice-timing offsets */
00405 
00406      if( nim->slice_dim      == 3               && /* AFNI can only deal with */
00407          nim->slice_code     >  0               && /* slice timing offsets    */
00408          nim->slice_duration >  0.0             && /* along the k-axis of     */
00409          nim->slice_start    >= 0               && /* the dataset volume      */
00410          nim->slice_start    < nim->nz          &&
00411          nim->slice_end      > nim->slice_start &&
00412          nim->slice_end      < nim->nz             ){
00413 
00414        float *toff=(float *)calloc(sizeof(float),nim->nz) , tsl ;
00415        int kk ;
00416 
00417             if( nim->time_units == NIFTI_UNITS_MSEC ) nim->slice_duration *= 0.001;
00418        else if( nim->time_units == NIFTI_UNITS_USEC ) nim->slice_duration *= 1.e-6;
00419 
00420        /* set up slice time offsets in the divers orders */
00421 
00422        switch( nim->slice_code ){
00423          case NIFTI_SLICE_SEQ_INC:
00424            tsl = 0.0 ;
00425            for( kk=nim->slice_start ; kk <= nim->slice_end ; kk++ ){
00426              toff[kk] = tsl ; tsl += nim->slice_duration ;
00427            }
00428          break ;
00429          case NIFTI_SLICE_SEQ_DEC:
00430            tsl = 0.0 ;
00431            for( kk=nim->slice_end ; kk >= nim->slice_end ; kk-- ){
00432              toff[kk] = tsl ; tsl += nim->slice_duration ;
00433            }
00434          break ;
00435          case NIFTI_SLICE_ALT_INC:
00436            tsl = 0.0 ;
00437            for( kk=nim->slice_start ; kk <= nim->slice_end ; kk+=2 ){
00438              toff[kk] = tsl ; tsl += nim->slice_duration ;
00439            }
00440            for( kk=nim->slice_start+1 ; kk <= nim->slice_end ; kk+=2 ){
00441              toff[kk] = tsl ; tsl += nim->slice_duration ;
00442            }
00443          break ;
00444          case NIFTI_SLICE_ALT_INC2:
00445            tsl = 0.0 ;
00446            for( kk=nim->slice_start+1 ; kk <= nim->slice_end ; kk+=2 ){
00447              toff[kk] = tsl ; tsl += nim->slice_duration ;
00448            }
00449            for( kk=nim->slice_start ; kk <= nim->slice_end ; kk+=2 ){
00450              toff[kk] = tsl ; tsl += nim->slice_duration ;
00451            }
00452          break ;
00453          case NIFTI_SLICE_ALT_DEC:
00454            tsl = 0.0 ;
00455            for( kk=nim->slice_end ; kk >= nim->slice_start ; kk-=2 ){
00456              toff[kk] = tsl ; tsl += nim->slice_duration ;
00457            }
00458            for( kk=nim->slice_end-1 ; kk >= nim->slice_start ; kk-=2 ){
00459              toff[kk] = tsl ; tsl += nim->slice_duration ;
00460            }
00461          break ;
00462          case NIFTI_SLICE_ALT_DEC2:
00463            tsl = 0.0 ;
00464            for( kk=nim->slice_end-1 ; kk >= nim->slice_start ; kk-=2 ){
00465              toff[kk] = tsl ; tsl += nim->slice_duration ;
00466            }
00467            for( kk=nim->slice_end ; kk >= nim->slice_start ; kk-=2 ){
00468              toff[kk] = tsl ; tsl += nim->slice_duration ;
00469            }
00470          break ;
00471        }
00472 
00473        EDIT_dset_items( dset ,
00474                           ADN_nsl     , nim->nz       ,
00475                           ADN_zorg_sl , orgxyz.xyz[2] ,
00476                           ADN_dz_sl   , dxyz.xyz[2]   ,
00477                           ADN_toff_sl , toff          ,
00478                         ADN_none ) ;
00479 
00480        free(toff) ;
00481 
00482      } /* end of slice timing stuff */
00483 
00484    } /* end of 3D+time dataset stuff */
00485 
00486    /* add statistics, if present */
00487 
00488    if( statcode != 0 ){
00489      for( ibr=0 ; ibr < nvals ; ibr++ )
00490        EDIT_STATAUX4(dset,ibr,statcode,nim->intent_p1,nim->intent_p2,nim->intent_p3,0) ;
00491    }
00492 
00493    /*-- flag to read data from disk using NIFTI functions --*/
00494 
00495    dset->dblk->diskptr->storage_mode = STORAGE_BY_NIFTI ;
00496    strcpy( dset->dblk->diskptr->brick_name , pathname ) ;
00497    dset->dblk->diskptr->byte_order = nim->byteorder ;
00498 
00499 #if 0
00500    for( ibr=0 ; ibr < nvals ; ibr++ ){     /* make sub-brick labels */
00501      sprintf(prefix,"%s[%d]",tname,ibr) ;
00502      EDIT_BRICK_LABEL( dset , ibr , prefix ) ;
00503    }
00504 #endif
00505 
00506    /** 10 May 2005: see if there is an AFNI extension;
00507                     if so, load attributes from it and
00508                     then edit the dataset appropriately **/
00509 
00510    { int ee ;  /* extension index */
00511 
00512      /* scan extension list to find the first AFNI extension */
00513 
00514      for( ee=0 ; ee < nim->num_ext ; ee++ )
00515        if( nim->ext_list[ee].ecode == NIFTI_ECODE_AFNI &&
00516            nim->ext_list[ee].esize > 32                &&
00517            nim->ext_list[ee].edata != NULL               ) break ;
00518 
00519      /* if found an AFNI extension ... */
00520 
00521      if( ee < nim->num_ext ){
00522        char *buf = nim->ext_list[ee].edata , *rhs , *cpt ;
00523        int  nbuf = nim->ext_list[ee].esize - 8 ;
00524        NI_stream ns ;
00525        void     *nini ;
00526        NI_group *ngr , *nngr ;
00527 
00528        /* if have data, it's long enough, and starts properly, then ... */
00529 
00530        if( buf != NULL && nbuf > 32 && strncmp(buf,"<?xml",5)==0 ){
00531          if( buf[nbuf-1] != '\0' ) buf[nbuf-1] = '\0' ;         /* for safety */
00532          cpt = strstr(buf,"?>") ;                    /* find XML prolog close */
00533          if( cpt != NULL ){                          /* if found it, then ... */
00534            ns = NI_stream_open( "str:" , "r" ) ;
00535            NI_stream_setbuf( ns , cpt+2 ) ;        /* start just after prolog */
00536            nini = NI_read_element(ns,1) ;                 /* get root element */
00537            NI_stream_close(ns) ;
00538            if( NI_element_type(nini) == NI_GROUP_TYPE ){   /* must be a group */
00539              ngr = (NI_group *)nini ;
00540              if( strcmp(ngr->name,"AFNI_attributes") == 0 ){    /* root is OK */
00541                nngr = ngr ;
00542              } else {                   /* search in group for proper element */
00543                int nn ; void **nnini ;
00544                nn = NI_search_group_deep( ngr , "AFNI_attributes" , &nnini ) ;
00545                if( nn <= 0 ) nngr = NULL ;
00546                else        { nngr = (NI_group *)nnini[0]; NI_free(nnini); }
00547              }
00548 
00549              if( NI_element_type(nngr) == NI_GROUP_TYPE ){ /* have  good name */
00550                rhs = NI_get_attribute( nngr , "self_idcode" ) ;
00551                if( rhs == NULL )
00552                  rhs = NI_get_attribute( nngr , "AFNI_idcode" ) ;
00553                if( rhs != NULL )    /* set dataset ID code from XML attribute */
00554                  MCW_strncpy( dset->idcode.str , rhs , MCW_IDSIZE ) ;
00555                rhs = NI_get_attribute( nngr , "NIfTI_nums" ) ;    /* check if */
00556                if( rhs != NULL ){                       /* dataset dimensions */
00557                  char buf[128] ;                              /* were altered */
00558                  sprintf(buf,"%d,%d,%d,%d,%d,%d" ,             /* 12 May 2005 */
00559                    nim->nx, nim->ny, nim->nz, nim->nt, nim->nu, nim->datatype );
00560                  if( strcmp(buf,rhs) != 0 ){
00561                    static int nnn=0 ;
00562                    if(nnn==0){fprintf(stderr,"\n"); nnn=1;}
00563                    fprintf(stderr,
00564                      "** WARNING: NIfTI file %s dimensions altered since "
00565                                  "AFNI extension was added\n",pathname ) ;
00566                  }
00567                }
00568                THD_dblkatr_from_niml( nngr , dset->dblk ); /* load attributes */
00569                THD_datablock_apply_atr( dset ) ;   /* apply to dataset struct */
00570              }
00571              NI_free_element( ngr ) ;          /* get rid of the root element */
00572 
00573            } /* end of if found a group element at the root */
00574          } /* end of if extension data array had an XML prolog close */
00575        } /* end of if had a good extension data array */
00576      } /* end of if had an AFNI extension */
00577    } /* end of processing extensions */
00578 
00579    /* return unpopulated dataset */
00580 
00581    nifti_image_free(nim) ; RETURN(dset) ;
00582 }
00583 
00584 /*-----------------------------------------------------------------
00585   Load a NIFTI dataset's data into memory
00586   (called from THD_load_datablock in thd_loaddblk.c)
00587     - RWC: modified 07 Apr 2005 to read data bricks via nifti_io.c
00588       'NBL' functions, rather than directly from disk
00589 -------------------------------------------------------------------*/
00590 
00591 void THD_load_nifti( THD_datablock *dblk )
00592 {
00593    THD_diskptr *dkptr ;
00594    int nx,ny,nz,nxy,nxyz,nxyzv , nerr=0,ibr,nv, nslice ;
00595    int datum, need_copy=0 ;
00596    void *ptr ;
00597    nifti_image *nim ;
00598    nifti_brick_list NBL ;  /* holds the data read from disk */
00599 
00600 ENTRY("THD_load_nifti") ;
00601 
00602    /*-- open and read input [these errors should never occur] --*/
00603 
00604    if( !ISVALID_DATABLOCK(dblk)                        ||
00605        dblk->diskptr->storage_mode != STORAGE_BY_NIFTI ||
00606        dblk->brick == NULL                               ) EXRETURN ;
00607 
00608    dkptr = dblk->diskptr ;
00609 
00610    STATUS("calling nifti_image_read_bricks") ;
00611    NBL.nbricks = 0 ;
00612    nim = nifti_image_read_bricks( dkptr->brick_name, 0,NULL , &NBL ) ;
00613    if( nim == NULL || NBL.nbricks <= 0 ) EXRETURN ;
00614 
00615    datum = DBLK_BRICK_TYPE(dblk,0) ;  /* destination data type */
00616 
00617    /*-- determine if we need to copy the data from the
00618         bricks as loaded above because of a type conversion --*/
00619 
00620    switch( nim->datatype ){
00621      case DT_INT16:
00622      case DT_UINT8:      need_copy = (datum == MRI_float) ; break ;
00623 
00624      case DT_FLOAT32:
00625      case DT_COMPLEX64:
00626      case DT_RGB24:      need_copy = 0 ; break ;
00627 
00628      case DT_INT8:       /* these are the cases where AFNI can't */
00629      case DT_UINT16:     /* directly handle the NIFTI datatype,  */
00630      case DT_INT32:      /* so we'll convert them to floats.     */
00631      case DT_UINT32:
00632      case DT_FLOAT64:    need_copy = 1 ; break ;
00633 #if 0
00634      case DT_COMPLEX128: need_copy = 1 ; break ;
00635 #endif
00636    }
00637 
00638    /*-- various dimensions --*/
00639 
00640    nx = dkptr->dimsizes[0] ;
00641    ny = dkptr->dimsizes[1] ;  nxy   = nx * ny   ;
00642    nz = dkptr->dimsizes[2] ;  nxyz  = nxy * nz  ;
00643    nv = dkptr->nvals       ;  if( nv > NBL.nbricks ) nv = NBL.nbricks ;
00644    nxyzv = nxyz * nv ; nslice = nz*nv ;
00645 
00646    dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
00647 
00648    /*------ don't need to copy data ==> just copy pointers from NBL ------*/
00649 
00650    if( !need_copy ){
00651 
00652      STATUS("copying brick pointers directly") ;
00653      for( ibr=0 ; ibr < nv ; ibr++ ){
00654        mri_fix_data_pointer( NBL.bricks[ibr] ,DBLK_BRICK(dblk,ibr) ) ;
00655        NBL.bricks[ibr] = NULL ;  /* so it won't be deleted later */
00656 
00657        if( AFNI_yesenv("AFNI_FLOATSCAN") ){  /*--- check float inputs? ---*/
00658          if( DBLK_BRICK_TYPE(dblk,ibr) == MRI_float )
00659            nerr += thd_floatscan( DBLK_BRICK_NVOX(dblk,ibr) ,
00660                                   DBLK_ARRAY(dblk,ibr)        ) ;
00661          else if( DBLK_BRICK_TYPE(dblk,ibr) == MRI_complex )
00662            nerr += thd_complexscan( DBLK_BRICK_NVOX(dblk,ibr) ,
00663                                     DBLK_ARRAY(dblk,ibr)        ) ;
00664        }
00665      }
00666      if( nerr > 0 ) fprintf(stderr ,
00667                     "** %s: found %d float errors -- see program float_scan\n",
00668                     dkptr->brick_name , nerr ) ;
00669 
00670    } else { /*---------- need to copy data ==> do some more work -----------*/
00671 
00672      register int ii ; void *nbuf ;
00673 
00674      STATUS("converting input bricks to floats") ;
00675      for( ibr=0 ; ibr < nv ; ibr++ ){
00676 
00677        if( DBLK_ARRAY(dblk,ibr) == NULL ){                     /* make space */
00678          ptr = AFMALL(void, DBLK_BRICK_BYTES(dblk,ibr) ) ;     /* for this   */
00679          mri_fix_data_pointer( ptr ,  DBLK_BRICK(dblk,ibr) ) ; /* sub-brick! */
00680        }
00681        ptr = DBLK_ARRAY(dblk,ibr) ; if( ptr == NULL ) break ;  /* bad news!! */
00682 
00683        nbuf = NBL.bricks[ibr] ;              /* data as read from NIfTI file */
00684 
00685        /* macro to convert data from type "ityp" in nbuf to float in dataset */
00686 
00687 #undef  CPF
00688 #define CPF(ityp) do{ ityp *sar = (ityp *)nbuf ; float *far = (float *)ptr ;   \
00689                       for( ii=0 ; ii < nxyz ; ii++ ) far[ii] = (float)sar[ii]; \
00690                   } while(0)
00691 
00692        /* load from nbuf into brick array (will be float or complex) */
00693 
00694        switch( nim->datatype ){
00695          case DT_UINT8:    CPF(unsigned char)  ; break ;
00696          case DT_INT8:     CPF(signed char)    ; break ;
00697          case DT_INT16:    CPF(signed short)   ; break ;
00698          case DT_UINT16:   CPF(unsigned short) ; break ;
00699          case DT_INT32:    CPF(signed int)     ; break ;
00700          case DT_UINT32:   CPF(unsigned int)   ; break ;
00701          case DT_FLOAT64:  CPF(double)         ; break ;
00702 #if 0
00703          case DT_COMPLEX128: break ;
00704 #endif
00705        }
00706 
00707        free(NBL.bricks[ibr]) ; NBL.bricks[ibr] = NULL ;
00708      }
00709    }
00710 
00711    nifti_free_NBL( &NBL ) ;  /* done with this */
00712 
00713    /*-- scale results? ---*/
00714 
00715    if( nim->scl_slope != 0.0 ){       /*--- scale results? ---*/
00716      STATUS("scaling sub-bricks") ;
00717      for( ibr=0 ; ibr < nv ; ibr++ ){
00718        if( DBLK_BRICK_TYPE(dblk,ibr) == MRI_float ){
00719          float *far = (float *) DBLK_ARRAY(dblk,ibr) ; int ii ;
00720          for( ii=0 ; ii < nxyz ; ii++ )
00721            far[ii] = nim->scl_slope * far[ii] + nim->scl_inter ;
00722        } else if( DBLK_BRICK_TYPE(dblk,ibr) == MRI_complex ){
00723          complex *car = (complex *) DBLK_ARRAY(dblk,ibr) ; int ii ;
00724          for( ii=0 ; ii < nxyz ; ii++ ){
00725            car[ii].r = nim->scl_slope * car[ii].r + nim->scl_inter ;
00726            car[ii].i = nim->scl_slope * car[ii].i + nim->scl_inter ;
00727          }
00728        }
00729      }
00730    }
00731 
00732    /*-- throw away the trash and return --*/
00733 
00734    nifti_image_free(nim) ; EXRETURN ;
00735 }