00001 #include "mrilib.h"
00002 #include "thd.h"
00003 
00004 /*******************************************************************/
00005 /* Stuff to read CTF MRI and SAM datasets from the NIH MEG system. */
00006 /* 04 Dec 2002: first cut by RWCox.                                */
00007 /*******************************************************************/
00008 
00009 /****************************************************************/
00010 /*********** header structs for CTF MRI version 2.0 file        */
00011 /****************************************************************/
00012 
00013 enum { coronal=0, sagittal, axial };
00014 
00015 #define LEFT_ON_LEFT      0
00016 #define LEFT_ON_RIGHT     1
00017 
00018 enum { Modality_MRI=0, Modality_CT, Modality_PET, Modality_SPECT, Modality_OTHER };
00019 
00020 #define VERSION_1_STR   "CTF_MRI_FORMAT VER 1.0"
00021 #define VERSION_2_STR   "CTF_MRI_FORMAT VER 2.0"
00022 #define VERSION_21_STR  "CTF_MRI_FORMAT VER 2.1"
00023 #define VERSION_22_STR  "CTF_MRI_FORMAT VER 2.2"
00024 
00025 typedef struct HeadModel_Info {
00026      short               Nasion_Sag;          /* fiduciary point voxel locations */
00027      short               Nasion_Cor;          /* Sag = sagittal direction */
00028      short               Nasion_Axi;          /* Cor = coronal direction */
00029      short               LeftEar_Sag;         /* Axi = axial direction */
00030      short               LeftEar_Cor;
00031      short               LeftEar_Axi;
00032      short               RightEar_Sag;
00033      short               RightEar_Cor;
00034      short               RightEar_Axi;
00035      float               defaultSphereX;      /* default sphere parameters in mm */
00036      float               defaultSphereY;      /* (in head based coordinate system */
00037      float               defaultSphereZ;
00038      float               defaultSphereRadius;
00039 } HeadModel_Info;
00040 
00041 /* this struct isn't used in AFNI */
00042 typedef struct Image_Info {                   /* scan and/or sequence parameters */
00043      short              modality;             /* 0=MRI, 1=CT, 2=PET, 3=SPECT, 4=OTHER */
00044      char               manufacturerName[64];
00045      char               instituteName[64];
00046      char               patientID[32];
00047      char               dateAndTime[32];
00048      char               scanType[32];
00049      char               contrastAgent[32];
00050      char               imagedNucleus[32];
00051      float              Frequency;
00052      float              FieldStrength;
00053      float              EchoTime;
00054      float              RepetitionTime;
00055      float              InversionTime;
00056      float              FlipAngle;
00057      short              NoExcitations;
00058      short              NoAcquisitions;
00059      char               commentString[256];
00060      char               forFutureUse[64];
00061 } Image_Info;
00062 
00063 /* the header for the .mri files */
00064 typedef struct Version_2_Header {
00065      char               identifierString[32];   /* "CTF_MRI_FORMAT VER 2.x"            */
00066      short              imageSize;              /* always = 256                        */
00067      short              dataSize;               /* 1 = 8 bit data, 2 = 16 bit data     */
00068      short              clippingRange;          /* max. integer value in data          */
00069      short              imageOrientation;       /* 0 = left on left, 1 = left on right */
00070      float              mmPerPixel_sagittal;    /* voxel dimensions in mm              */
00071      float              mmPerPixel_coronal;     /* voxel dimensions in mm              */
00072      float              mmPerPixel_axial;       /* voxel dimensions in mm              */
00073      HeadModel_Info     headModel;              /* structure defined above (34 bytes)  */
00074      Image_Info         imageInfo;              /* structure defined above (638 bytes) */
00075      float              headOrigin_sagittal;    /* voxel location of head origin       */
00076      float              headOrigin_coronal;     /* voxel location of head origin       */
00077      float              headOrigin_axial;       /* voxel location of head origin       */
00078 
00079                                                 /* Euler angles to align MR to head          */
00080                                                 /* coordinate system (angles in degrees!)    */
00081      float              rotate_coronal;         /* 1. rotate in coronal plane by this angle  */
00082      float              rotate_sagittal;        /* 2. rotate in sagittal plane by this angle */
00083      float              rotate_axial;           /* 3. rotate in axial plane by this angle    */
00084 
00085      short              orthogonalFlag;         /* 1 if image is orthogonalized to head frame */
00086      short              interpolatedFlag;       /* 1 if slices interpolated during conversion */
00087      float              originalSliceThickness;
00088      float              transformMatrix[4][4];  /* 4x4 transformation matrix (head to mri) */
00089      unsigned char      unused[202];            /* pad header to 1028 bytes                */
00090 } Version_2_Header;
00091 
00092 /*---------------------------------------------------------------*/
00093 /*! Swap the 4 bytes pointed to by ppp: abcd -> dcba. */
00094 
00095 static void swap_4(void *ppp)
00096 {
00097    unsigned char *pntr = (unsigned char *) ppp ;
00098    unsigned char b0, b1, b2, b3;
00099 
00100    b0 = *pntr; b1 = *(pntr+1); b2 = *(pntr+2); b3 = *(pntr+3);
00101    *pntr = b3; *(pntr+1) = b2; *(pntr+2) = b1; *(pntr+3) = b0;
00102 }
00103 
00104 /*---------------------------------------------------------------*/
00105 /*! Swap the 8 bytes pointed to by ppp: abcdefgh -> hgfedcba. */
00106 
00107 static void swap_8(void *ppp)
00108 {
00109    unsigned char *pntr = (unsigned char *) ppp ;
00110    unsigned char b0, b1, b2, b3;
00111    unsigned char b4, b5, b6, b7;
00112 
00113    b0 = *pntr    ; b1 = *(pntr+1); b2 = *(pntr+2); b3 = *(pntr+3);
00114    b4 = *(pntr+4); b5 = *(pntr+5); b6 = *(pntr+6); b7 = *(pntr+7);
00115 
00116    *pntr     = b7; *(pntr+1) = b6; *(pntr+2) = b5; *(pntr+3) = b4;
00117    *(pntr+4) = b3; *(pntr+5) = b2; *(pntr+6) = b1; *(pntr+7) = b0;
00118 }
00119 
00120 /*---------------------------------------------------------------*/
00121 /*! Swap the 2 bytes pointed to by ppp: ab -> ba. */
00122 
00123 static void swap_2(void *ppp)
00124 {
00125    unsigned char *pntr = (unsigned char *) ppp ;
00126    unsigned char b0, b1;
00127 
00128    b0 = *pntr; b1 = *(pntr+1);
00129    *pntr = b1; *(pntr+1) = b0;
00130 }
00131 
00132 /*-------------------------*/
00133 /*! Macro for bad return.  */
00134 
00135 #undef  BADBAD
00136 #define BADBAD(s)                                                \
00137   do{ fprintf(stderr,"** THD_open_ctfmri(%s): %s\n",fname,s);    \
00138       RETURN(NULL);                                              \
00139   } while(0)
00140 
00141 /*-----------------------------------------------------------------*/
00142 /*! Function to count slices like CTF does. */
00143 
00144 int CTF_count( double start, double end , double delta )
00145 {
00146    int nn=0 ; double cc ;
00147    for( cc=start ; cc <= (end+1.0e-6) ; cc += delta ) nn++ ;
00148    return nn ;
00149 }
00150 
00151 /*-----------------------------------------------------------------*/
00152 /*! Open a CTF .mri file as an unpopulated AFNI dataset.
00153     It will be populated later, in THD_load_ctfmri().
00154 -------------------------------------------------------------------*/
00155 
00156 THD_3dim_dataset * THD_open_ctfmri( char *fname )
00157 {
00158    FILE *fp ;
00159    Version_2_Header hh ;
00160    int ii,nn , swap ;
00161    THD_3dim_dataset *dset=NULL ;
00162    char prefix[THD_MAX_PREFIX] , *ppp , tname[12] , ori[4] ;
00163    THD_ivec3 nxyz , orixyz ;
00164    THD_fvec3 dxyz , orgxyz ;
00165    int iview ;
00166    int ngood , length , datum_type , datum_len , oxx,oyy,ozz ;
00167    int   nx,ny,nz ;
00168    float dx,dy,dz , xorg,yorg,zorg ;
00169 
00170 ENTRY("THD_open_ctfmri") ;
00171 
00172    /* open input file */
00173 
00174    if( fname == NULL || *fname == '\0' ) BADBAD("bad input filename");
00175    fp = fopen( fname , "rb" ) ;
00176    if( fp == NULL )                      BADBAD("can't open input file");
00177 
00178    /* read 1028 byte header */
00179 
00180    nn = fread( &hh , sizeof(hh) , 1 , fp ) ;
00181    fclose(fp) ;
00182    if( nn <= 0 )                         BADBAD("can't read input file");
00183 
00184    /* check if header string matches what we want */
00185 
00186    hh.identifierString[31] = '\0' ;  /* make sure is terminated */
00187    if( strcmp(hh.identifierString,VERSION_22_STR) != 0 ) BADBAD("bad version string");
00188    if( hh.imageSize == 0 )                               BADBAD("bad imageSize") ;
00189 
00190    /* determine if must swap header */
00191 
00192    swap = (hh.imageSize != 256) ;
00193 
00194    if( swap ){                          /* swap bytes in various header fields */
00195      swap_2(&hh.imageSize              ) ;
00196      swap_2(&hh.dataSize               ) ;
00197      swap_2(&hh.clippingRange          ) ;
00198      swap_2(&hh.imageOrientation       ) ;
00199      swap_4(&hh.mmPerPixel_sagittal    ) ;
00200      swap_4(&hh.mmPerPixel_coronal     ) ;
00201      swap_4(&hh.mmPerPixel_axial       ) ;
00202      swap_4(&hh.headOrigin_sagittal    ) ;
00203      swap_4(&hh.headOrigin_coronal     ) ;
00204      swap_4(&hh.headOrigin_axial       ) ;
00205      swap_4(&hh.rotate_coronal         ) ;
00206      swap_4(&hh.rotate_sagittal        ) ;
00207      swap_4(&hh.rotate_axial           ) ;
00208      swap_2(&hh.orthogonalFlag         ) ;
00209      swap_2(&hh.interpolatedFlag       ) ;
00210      swap_4(&hh.originalSliceThickness ) ;
00211      swap_2(&hh.headModel.Nasion_Sag   ) ;
00212      swap_2(&hh.headModel.Nasion_Cor   ) ;
00213      swap_2(&hh.headModel.Nasion_Axi   ) ;
00214      swap_2(&hh.headModel.LeftEar_Sag  ) ;
00215      swap_2(&hh.headModel.LeftEar_Cor  ) ;
00216      swap_2(&hh.headModel.LeftEar_Axi  ) ;
00217      swap_2(&hh.headModel.RightEar_Sag ) ;
00218      swap_2(&hh.headModel.RightEar_Cor ) ;
00219      swap_2(&hh.headModel.RightEar_Axi ) ;
00220 
00221      swap_4(&hh.transformMatrix[0][0]  ) ;   /* this stuff not used yet */
00222      swap_4(&hh.transformMatrix[0][1]  ) ;
00223      swap_4(&hh.transformMatrix[0][2]  ) ;
00224      swap_4(&hh.transformMatrix[0][3]  ) ;
00225      swap_4(&hh.transformMatrix[1][0]  ) ;
00226      swap_4(&hh.transformMatrix[1][1]  ) ;
00227      swap_4(&hh.transformMatrix[1][2]  ) ;
00228      swap_4(&hh.transformMatrix[1][3]  ) ;
00229      swap_4(&hh.transformMatrix[2][0]  ) ;
00230      swap_4(&hh.transformMatrix[2][1]  ) ;
00231      swap_4(&hh.transformMatrix[2][2]  ) ;
00232      swap_4(&hh.transformMatrix[2][3]  ) ;
00233      swap_4(&hh.transformMatrix[3][0]  ) ;
00234      swap_4(&hh.transformMatrix[3][1]  ) ;
00235      swap_4(&hh.transformMatrix[3][2]  ) ;
00236      swap_4(&hh.transformMatrix[3][3]  ) ;
00237    }
00238 
00239    /* simple checks on header stuff */
00240 
00241    if( hh.imageSize != 256           ||
00242        hh.dataSize  <  1             ||
00243        hh.dataSize  >  2             ||
00244        hh.mmPerPixel_sagittal <= 0.0 ||
00245        hh.mmPerPixel_coronal  <= 0.0 ||
00246        hh.mmPerPixel_axial    <= 0.0   ) BADBAD("bad header data") ;
00247 
00248    /*- 16 Mar 2005: instead of complaining about negative Origins,
00249                     just reset them to something semi-reasonable;
00250        I just get by with a little help from my friends
00251        - Zuxiang Li in this case.                                -*/
00252 
00253    if( hh.headOrigin_sagittal <= 0.0 ) hh.headOrigin_sagittal = hh.imageSize/2.;
00254    if( hh.headOrigin_coronal  <= 0.0 ) hh.headOrigin_coronal  = hh.imageSize/2.;
00255    if( hh.headOrigin_axial    <= 0.0 ) hh.headOrigin_axial    = hh.imageSize/2.;
00256 
00257    /* debugging code to print header information */
00258 #if 0
00259    printf("\n") ;
00260    printf("*** CTF MRI filename   = %s\n",fname                 ) ;
00261    printf("identifierString       = %s\n",hh.identifierString   ) ;
00262    printf("imageSize              = %d\n",hh.imageSize          ) ;
00263    printf("dataSize               = %d\n",hh.dataSize           ) ;
00264    printf("clippingRange          = %d\n",hh.clippingRange      ) ;
00265    printf("imageOrientation       = %d\n",hh.imageOrientation   ) ;
00266    printf("mmPerPixel_sagittal    = %f\n",hh.mmPerPixel_sagittal) ;
00267    printf("mmPerPixel_coronal     = %f\n",hh.mmPerPixel_coronal ) ;
00268    printf("mmPerPixel_axial       = %f\n",hh.mmPerPixel_axial   ) ;
00269    printf("headOrigin_sagittal    = %f\n",hh.headOrigin_sagittal) ;
00270    printf("headOrigin_coronal     = %f\n",hh.headOrigin_coronal ) ;
00271    printf("headOrigin_axial       = %f\n",hh.headOrigin_axial   ) ;
00272    printf("rotate_coronal         = %f\n",hh.rotate_coronal     ) ;
00273    printf("rotate_sagittal        = %f\n",hh.rotate_sagittal    ) ;
00274    printf("rotate_axial           = %f\n",hh.rotate_axial       ) ;
00275    printf("orthogonalFlag         = %d\n",hh.orthogonalFlag     ) ;
00276    printf("interpolatedFlag       = %d\n",hh.interpolatedFlag   ) ;
00277    printf("originalSliceThickness = %f\n",hh.originalSliceThickness ) ;
00278    printf("\n") ;
00279    printf("headModel.Nasion_Sag   = %d\n",hh.headModel.Nasion_Sag  ) ;
00280    printf("headModel.Nasion_Cor   = %d\n",hh.headModel.Nasion_Cor  ) ;
00281    printf("headModel.Nasion_Axi   = %d\n",hh.headModel.Nasion_Axi  ) ;
00282    printf("headModel.LeftEar_Sag  = %d\n",hh.headModel.LeftEar_Sag ) ;
00283    printf("headModel.LeftEar_Cor  = %d\n",hh.headModel.LeftEar_Cor ) ;
00284    printf("headModel.LeftEar_Axi  = %d\n",hh.headModel.LeftEar_Axi ) ;
00285    printf("headModel.RightEar_Sag = %d\n",hh.headModel.RightEar_Sag) ;
00286    printf("headModel.RightEar_Cor = %d\n",hh.headModel.RightEar_Cor) ;
00287    printf("headModel.RightEar_Axi = %d\n",hh.headModel.RightEar_Axi) ;
00288    printf("\n") ;
00289    printf("transformMatrix:\n"
00290           "  [ %9.4f %9.4f %9.4f %9.4f ]\n"
00291           "  [ %9.4f %9.4f %9.4f %9.4f ]\n"
00292           "  [ %9.4f %9.4f %9.4f %9.4f ]\n"
00293           "  [ %9.4f %9.4f %9.4f %9.4f ]\n" ,
00294       hh.transformMatrix[0][0] , hh.transformMatrix[0][1] ,
00295         hh.transformMatrix[0][2] , hh.transformMatrix[0][3] ,
00296       hh.transformMatrix[1][0] , hh.transformMatrix[1][1] ,
00297         hh.transformMatrix[1][2] , hh.transformMatrix[1][3] ,
00298       hh.transformMatrix[2][0] , hh.transformMatrix[2][1] ,
00299         hh.transformMatrix[2][2] , hh.transformMatrix[2][3] ,
00300       hh.transformMatrix[3][0] , hh.transformMatrix[3][1] ,
00301         hh.transformMatrix[3][2] , hh.transformMatrix[3][3]  ) ;
00302 #endif
00303 
00304    /* determine if file is big enough to hold all data it claims */
00305 
00306    nn = THD_filesize(fname) ;
00307    if( nn < hh.dataSize*hh.imageSize*hh.imageSize*hh.imageSize )
00308      BADBAD("input file too small") ;
00309 
00310    /*** from here, a lot of code is adapted from thd_analyzeread.c ***/
00311 
00312    datum_len = hh.dataSize ;
00313    switch( datum_len ){                         /* the only 2 cases */
00314      case 1:  datum_type = MRI_byte ; break ;
00315      case 2:  datum_type = MRI_short; break ;
00316    }
00317    nx = ny = nz = hh.imageSize ;              /* volumes are cubes! */
00318 
00319    /* set orientation:
00320       for now, assume (based on 1 sample) that data is stored in ASL or ASR order */
00321 
00322    ori[0] = 'A' ;          /* x is A-P */
00323    ori[1] = 'S' ;          /* y is S-I */
00324 
00325    /* determine if z is L-R or R-L from position of markers */
00326 
00327    ori[2] = (hh.headModel.LeftEar_Sag <= hh.headModel.RightEar_Sag) ? 'L' : 'R' ;
00328 
00329    oxx = ORCODE(ori[0]); oyy = ORCODE(ori[1]); ozz = ORCODE(ori[2]);
00330    if( !OR3OK(oxx,oyy,ozz) ){
00331      oxx = ORI_A2P_TYPE; oyy = ORI_S2I_TYPE; ozz = ORI_L2R_TYPE;   /** ASL? **/
00332    }
00333 
00334    /* now set grid size, keeping in mind that
00335        A-P is positive and P-A is negative,
00336        R-L is positive and L-R is negative,
00337        I-S is positive and S-I is negative.   */
00338 
00339    switch( ori[0] ){
00340      case 'A':  dx =  hh.mmPerPixel_coronal ; xorg = hh.headOrigin_coronal ; break ;
00341      case 'P':  dx = -hh.mmPerPixel_coronal ; xorg = hh.headOrigin_coronal ; break ;
00342      case 'R':  dx =  hh.mmPerPixel_sagittal; xorg = hh.headOrigin_sagittal; break ;
00343      case 'L':  dx = -hh.mmPerPixel_sagittal; xorg = hh.headOrigin_sagittal; break ;
00344      case 'I':  dx =  hh.mmPerPixel_axial   ; xorg = hh.headOrigin_axial   ; break ;
00345      case 'S':  dx = -hh.mmPerPixel_axial   ; xorg = hh.headOrigin_axial   ; break ;
00346    }
00347    switch( ori[1] ){
00348      case 'A':  dy =  hh.mmPerPixel_coronal ; yorg = hh.headOrigin_coronal ; break ;
00349      case 'P':  dy = -hh.mmPerPixel_coronal ; yorg = hh.headOrigin_coronal ; break ;
00350      case 'R':  dy =  hh.mmPerPixel_sagittal; yorg = hh.headOrigin_sagittal; break ;
00351      case 'L':  dy = -hh.mmPerPixel_sagittal; yorg = hh.headOrigin_sagittal; break ;
00352      case 'I':  dy =  hh.mmPerPixel_axial   ; yorg = hh.headOrigin_axial   ; break ;
00353      case 'S':  dy = -hh.mmPerPixel_axial   ; yorg = hh.headOrigin_axial   ; break ;
00354    }
00355    switch( ori[2] ){
00356      case 'A':  dz =  hh.mmPerPixel_coronal ; zorg = hh.headOrigin_coronal ; break ;
00357      case 'P':  dz = -hh.mmPerPixel_coronal ; zorg = hh.headOrigin_coronal ; break ;
00358      case 'R':  dz =  hh.mmPerPixel_sagittal; zorg = hh.headOrigin_sagittal; break ;
00359      case 'L':  dz = -hh.mmPerPixel_sagittal; zorg = hh.headOrigin_sagittal; break ;
00360      case 'I':  dz =  hh.mmPerPixel_axial   ; zorg = hh.headOrigin_axial   ; break ;
00361      case 'S':  dz = -hh.mmPerPixel_axial   ; zorg = hh.headOrigin_axial   ; break ;
00362    }
00363 
00364    /* At this point, (xorg,yorg,zorg) are voxel indices;
00365       now, translate them into shifts such that if voxel
00366       index ii is the location of x=0, then xorg+ii*dx=0. */
00367 
00368    xorg = -dx*xorg ; yorg = -dy*yorg ; zorg = -dz*zorg ;
00369 
00370    /*-- make a dataset --*/
00371 
00372    dset = EDIT_empty_copy(NULL) ;
00373 
00374    dset->idcode.str[0] = 'C' ;  /* overwrite 1st 3 bytes */
00375    dset->idcode.str[1] = 'T' ;
00376    dset->idcode.str[2] = 'F' ;
00377 
00378    MCW_hash_idcode( fname , dset ) ;  /* 06 May 2005 */
00379 
00380    ppp = THD_trailname(fname,0) ;                   /* strip directory */
00381    MCW_strncpy( prefix , ppp , THD_MAX_PREFIX ) ;   /* to make prefix */
00382 
00383    nxyz.ijk[0] = nx ; dxyz.xyz[0] = dx ;  /* setup axes lengths and voxel sizes */
00384    nxyz.ijk[1] = ny ; dxyz.xyz[1] = dy ;
00385    nxyz.ijk[2] = nz ; dxyz.xyz[2] = dz ;
00386 
00387    orixyz.ijk[0] = oxx ; orgxyz.xyz[0] = xorg ;
00388    orixyz.ijk[1] = oyy ; orgxyz.xyz[1] = yorg ;
00389    orixyz.ijk[2] = ozz ; orgxyz.xyz[2] = zorg ;
00390 
00391    iview = VIEW_ORIGINAL_TYPE ;
00392 
00393    /*-- actually send the values above into the dataset header --*/
00394 
00395    EDIT_dset_items( dset ,
00396                       ADN_prefix      , prefix ,
00397                       ADN_datum_all   , datum_type ,
00398                       ADN_nxyz        , nxyz ,
00399                       ADN_xyzdel      , dxyz ,
00400                       ADN_xyzorg      , orgxyz ,
00401                       ADN_xyzorient   , orixyz ,
00402                       ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00403                       ADN_nvals       , 1 ,
00404                       ADN_type        , HEAD_ANAT_TYPE ,
00405                       ADN_view_type   , iview ,
00406                       ADN_func_type   , ANAT_MRAN_TYPE ,
00407                     ADN_none ) ;
00408 
00409    /*-- set byte order (for reading from disk) --*/
00410 
00411    ii = mri_short_order() ;
00412    if( swap ) ii = REVERSE_ORDER(ii) ;
00413    dset->dblk->diskptr->byte_order = ii ;
00414 
00415    /*-- flag to read data from disk using CTF MRI mode --*/
00416 
00417    dset->dblk->diskptr->storage_mode = STORAGE_BY_CTFMRI ;
00418    strcpy( dset->dblk->diskptr->brick_name , fname ) ;
00419 
00420    /*-- for fun, add a set of tags for MEG fiducial points, if present --*/
00421 
00422    if( hh.headModel.LeftEar_Sag != hh.headModel.RightEar_Sag ){
00423      THD_usertaglist *tagset = myXtNew(THD_usertaglist) ;
00424      int nas_ii,nas_jj,nas_kk , lft_ii,lft_jj,lft_kk , rgt_ii,rgt_jj,rgt_kk ;
00425      THD_fvec3 fv ; THD_ivec3 iv ;
00426 
00427      tagset->num = 3 ;
00428      TAGLIST_SETLABEL( tagset , "CTF MEG Fiducials" ) ;
00429 
00430      /* load voxel indexes into dataset of the 3 tag points;
00431         note we have to permute these into the dataset axes order */
00432 
00433      switch( ori[0] ){
00434        case 'P':
00435        case 'A':  nas_ii = hh.headModel.Nasion_Cor ;
00436                   lft_ii = hh.headModel.LeftEar_Cor ;
00437                   rgt_ii = hh.headModel.RightEar_Cor ; break ;
00438        case 'R':
00439        case 'L':  nas_ii = hh.headModel.Nasion_Sag ;
00440                   lft_ii = hh.headModel.LeftEar_Sag ;
00441                   rgt_ii = hh.headModel.RightEar_Sag ; break ;
00442        case 'I':
00443        case 'S':  nas_ii = hh.headModel.Nasion_Axi ;
00444                   lft_ii = hh.headModel.LeftEar_Axi ;
00445                   rgt_ii = hh.headModel.RightEar_Axi ; break ;
00446      }
00447      switch( ori[1] ){
00448        case 'P':
00449        case 'A':  nas_jj = hh.headModel.Nasion_Cor ;
00450                   lft_jj = hh.headModel.LeftEar_Cor ;
00451                   rgt_jj = hh.headModel.RightEar_Cor ; break ;
00452        case 'R':
00453        case 'L':  nas_jj = hh.headModel.Nasion_Sag ;
00454                   lft_jj = hh.headModel.LeftEar_Sag ;
00455                   rgt_jj = hh.headModel.RightEar_Sag ; break ;
00456        case 'I':
00457        case 'S':  nas_jj = hh.headModel.Nasion_Axi ;
00458                   lft_jj = hh.headModel.LeftEar_Axi ;
00459                   rgt_jj = hh.headModel.RightEar_Axi ; break ;
00460      }
00461      switch( ori[2] ){
00462        case 'P':
00463        case 'A':  nas_kk = hh.headModel.Nasion_Cor ;
00464                   lft_kk = hh.headModel.LeftEar_Cor ;
00465                   rgt_kk = hh.headModel.RightEar_Cor ; break ;
00466        case 'R':
00467        case 'L':  nas_kk = hh.headModel.Nasion_Sag ;
00468                   lft_kk = hh.headModel.LeftEar_Sag ;
00469                   rgt_kk = hh.headModel.RightEar_Sag ; break ;
00470        case 'I':
00471        case 'S':  nas_kk = hh.headModel.Nasion_Axi ;
00472                   lft_kk = hh.headModel.LeftEar_Axi ;
00473                   rgt_kk = hh.headModel.RightEar_Axi ; break ;
00474      }
00475 
00476      TAG_SETLABEL( tagset->tag[0] , "Nasion" ) ;
00477      LOAD_IVEC3( iv , nas_ii,nas_jj,nas_kk ) ;  /* compute DICOM  */
00478      fv = THD_3dind_to_3dmm( dset , iv ) ;      /* coordinates of */
00479      fv = THD_3dmm_to_dicomm( dset , fv ) ;     /* this point     */
00480      UNLOAD_FVEC3( fv , tagset->tag[0].x , tagset->tag[0].y , tagset->tag[0].z ) ;
00481      tagset->tag[0].val = 0.0 ;
00482      tagset->tag[0].ti  = 0 ;
00483      tagset->tag[0].set = 1 ;
00484 
00485      TAG_SETLABEL( tagset->tag[1] , "Left Ear" ) ;
00486      LOAD_IVEC3( iv , lft_ii,lft_jj,lft_kk ) ;
00487      fv = THD_3dind_to_3dmm( dset , iv ) ;
00488      fv = THD_3dmm_to_dicomm( dset , fv ) ;
00489      UNLOAD_FVEC3( fv , tagset->tag[1].x , tagset->tag[1].y , tagset->tag[1].z ) ;
00490      tagset->tag[1].val = 0.0 ;
00491      tagset->tag[1].ti  = 0 ;
00492      tagset->tag[1].set = 1 ;
00493 
00494      TAG_SETLABEL( tagset->tag[2] , "Right Ear" ) ;
00495      LOAD_IVEC3( iv , rgt_ii,rgt_jj,rgt_kk ) ;
00496      fv = THD_3dind_to_3dmm( dset , iv ) ;
00497      fv = THD_3dmm_to_dicomm( dset , fv ) ;
00498      UNLOAD_FVEC3( fv , tagset->tag[2].x , tagset->tag[2].y , tagset->tag[2].z ) ;
00499      tagset->tag[2].val = 0.0 ;
00500      tagset->tag[2].ti  = 0 ;
00501      tagset->tag[2].set = 1 ;
00502 
00503      dset->tagset = tagset ;
00504    }
00505 
00506    RETURN(dset) ;
00507 }
00508 
00509 /*------------------------------------------------------------------*/
00510 /*! Actually load data from a CTF MRI file into a dataset.
00511     Adapted from THD_load_analyze().
00512 --------------------------------------------------------------------*/
00513 
00514 void THD_load_ctfmri( THD_datablock *dblk )
00515 {
00516    THD_diskptr *dkptr ;
00517    int nx,ny,nz,nv , nxy,nxyz,nxyzv , ibr,nbad ;
00518    FILE *fp ;
00519    void *ptr ;
00520 
00521 ENTRY("THD_load_ctfmri") ;
00522 
00523    /*-- check inputs --*/
00524 
00525    if( !ISVALID_DATABLOCK(dblk)                         ||
00526        dblk->diskptr->storage_mode != STORAGE_BY_CTFMRI ||
00527        dblk->brick == NULL                                ) EXRETURN ;
00528 
00529    dkptr = dblk->diskptr ;
00530 
00531    /* open and position file at start of data (after header) */
00532 
00533    fp = fopen( dkptr->brick_name , "rb" ) ;  /* .mri file */
00534    if( fp == NULL ) EXRETURN ;
00535 
00536    /*-- allocate space for data --*/
00537 
00538    nx = dkptr->dimsizes[0] ;
00539    ny = dkptr->dimsizes[1] ; nxy   = nx * ny   ;
00540    nz = dkptr->dimsizes[2] ; nxyz  = nxy * nz  ;
00541    nv = dkptr->nvals       ; nxyzv = nxyz * nv ;
00542 
00543    /* 26 Feb 2005: seek backwards from end,
00544                    instead of forwards from start */
00545 
00546 #if 0
00547    fseek( fp , sizeof(Version_2_Header) , SEEK_SET ) ;  /* old */
00548 #else
00549    switch( DBLK_BRICK_TYPE(dblk,0) ){
00550      case MRI_float: ibr = sizeof(float) ; break ;   /* illegal */
00551      case MRI_short: ibr = sizeof(short) ; break ;
00552      case MRI_byte:  ibr = sizeof(byte)  ; break ;
00553    }
00554    fseek( fp , -ibr*nxyzv , SEEK_END ) ;                /* new */
00555 #endif
00556 
00557    dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
00558 
00559    /*-- malloc space for each brick separately --*/
00560 
00561    for( nbad=ibr=0 ; ibr < nv ; ibr++ ){
00562      if( DBLK_ARRAY(dblk,ibr) == NULL ){
00563        ptr = AFMALL(void, DBLK_BRICK_BYTES(dblk,ibr) ) ;
00564        mri_fix_data_pointer( ptr ,  DBLK_BRICK(dblk,ibr) ) ;
00565        if( ptr == NULL ) nbad++ ;
00566      }
00567    }
00568 
00569    /*-- if couldn't get them all, take our ball and go home in a snit --*/
00570 
00571    if( nbad > 0 ){
00572      fprintf(stderr,
00573              "\n** failed to malloc %d CTR MRI bricks out of %d\n\a",nbad,nv);
00574      for( ibr=0 ; ibr < nv ; ibr++ ){
00575        if( DBLK_ARRAY(dblk,ibr) != NULL ){
00576          free(DBLK_ARRAY(dblk,ibr)) ;
00577          mri_fix_data_pointer( NULL , DBLK_BRICK(dblk,ibr) ) ;
00578        }
00579      }
00580      fclose(fp) ; EXRETURN ;
00581    }
00582 
00583 
00584    /*-- read data from .img file into sub-brick arrays! --*/
00585 
00586    for( ibr=0 ; ibr < nv ; ibr++ )
00587      fread( DBLK_ARRAY(dblk,ibr), 1, DBLK_BRICK_BYTES(dblk,ibr), fp ) ;
00588 
00589    fclose(fp) ;
00590 
00591    /*-- swap bytes? --*/
00592 
00593    if( dkptr->byte_order != mri_short_order() ){
00594      for( ibr=0 ; ibr < nv ; ibr++ ){
00595        switch( DBLK_BRICK_TYPE(dblk,ibr) ){
00596          case MRI_short:
00597            mri_swap2( DBLK_BRICK_NVOX(dblk,ibr) , DBLK_ARRAY(dblk,ibr) ) ;
00598          break ;
00599        }
00600      }
00601    }
00602 
00603    EXRETURN ;
00604 }
00605 
00606 /****************************************************************************
00607  SAM static image files are structured as follows:
00608 
00609   char     Identity[8] = "SAMIMAGE"; // uniquely identifies image file
00610   SAM_HDR  SAMHeader;                // SAM header
00611   double   Voxel[0];                 // 1st SAM voxel (units=A-m, (A-m)^2,
00612   double   Voxel[1];                 // 2nd SAM voxel        Z, T, F, or P)
00613               "
00614               "
00615   double   Voxel[V];                 // last SAM voxel
00616 
00617   Coefficients & image voxels are ordered in X,Y,Z sequence, with Z the least
00618   significant index (most rapidly changing), Y is next, and then X.
00619   Coordinate indices always advance in the positive direction. This implies
00620   that Voxel[0] is in the right, posterior, inferior position relative to
00621   the region of interest (bounding box of image).
00622 
00623   RWCox: the data storage order seems to be IRP, based on the above
00624          comments, and on the CTF head coordinates system being PRI.
00625 *****************************************************************************/
00626 
00627 #define COV_VERSION      1                  /* this is version 1 -- got it? */
00628 #define SAM_VERSION      1                  /* this, too! */
00629 
00630 /** SAM file types **/
00631 
00632 #define SAM_TYPE_IMAGE         0            /* flags file as a SAM static image file */
00633 #define SAM_TYPE_WT_ARRAY      1            /* flags file as SAM coefficients for regular target array */
00634 #define SAM_TYPE_WT_LIST       2            /* flags file as SAM coefficients for target list */
00635 
00636 /** define SAM unit types **/
00637 
00638 #define      SAM_UNIT_COEFF    0            /* SAM coefficients A-m/T */
00639 #define      SAM_UNIT_MOMENT   1            /* SAM source (or noise) strength A-m */
00640 #define      SAM_UNIT_POWER    2            /* SAM source (or noise) power (A-m)^2 */
00641 #define      SAM_UNIT_SPMZ     3            /* SAM z-deviate */
00642 #define      SAM_UNIT_SPMF     4            /* SAM F-statistic */
00643 #define      SAM_UNIT_SPMT     5            /* SAM T-statistic */
00644 #define      SAM_UNIT_SPMP     6            /* SAM probability */
00645 #define      SAM_UNIT_MUSIC    7            /* MUSIC metric */
00646 
00647 /* 'SAM_HDR' is to be used for both SAM coefficients & SAM static images */
00648 typedef struct {
00649    int    Version;       /* file version number (should be 1) */
00650    char   SetName[256];  /* name of parent dataset */
00651    int    NumChans;      /* # of channels used by SAM */
00652    int    NumWeights;    /* # of SAM virtual channels (0=static image) */
00653    int    pad_bytes1;    /* ** align next double on 8 byte boundary */
00654    double XStart;        /* x-start coordinate (m) */
00655    double XEnd;          /* x-end coordinate (m) */
00656    double YStart;        /* y-start coordinate (m) */
00657    double YEnd;          /* y-end coordinate (m) */
00658    double ZStart;        /* z-start coordinate (m) */
00659    double ZEnd;          /* z-end coordinate (m) */
00660    double StepSize;      /* voxel step size (m) */
00661    double HPFreq;        /* highpass frequency (Hz) */
00662    double LPFreq;        /* lowpass frequency (Hz) */
00663    double BWFreq;        /* bandwidth of filters (Hz) */
00664    double MeanNoise;     /* mean primary sensor noise (T) */
00665    char   MriName[256];  /* MRI image file name */
00666    int    Nasion[3];     /* MRI voxel index for nasion */
00667    int    RightPA[3];    /* MRI voxel index for right pre-auricular */
00668    int    LeftPA[3];     /* MRI voxel index for left pre-auricular */
00669    int    SAMType;       /* SAM file type */
00670    int    SAMUnit;       /* SAM units (a bit redundant, but may be useful) */
00671    int    pad_bytes2;    /* ** align end of structure on 8 byte boundary */
00672 } SAM_HDR;
00673 
00674 /*** 26 Feb 2005: version 2 of the SAM header ***/
00675 
00676 #if 0
00677 typedef struct {
00678    int     Version;         /* file version number (should be 2) */
00679    char    SetName[256];    /* name of parent dataset */
00680    int     NumChans;        /* # of channels used by SAM */
00681    int     NumWeights;      /* # of SAM virtual channels (0=static image) */
00682    int     pad_bytes1;      /* ** align next double on 8 byte boundary */
00683    double  XStart;          /* x-start coordinate (m) */
00684    double  XEnd;            /* x-end coordinate (m) */
00685    double  YStart;          /* y-start coordinate (m) */
00686    double  YEnd;            /* y-end coordinate (m) */
00687    double  ZStart;          /* z-start coordinate (m) */
00688    double  ZEnd;            /* z-end coordinate (m) */
00689    double  StepSize;        /* voxel step size (m) */
00690    double  HPFreq;          /* highpass frequency (Hz) */
00691    double  LPFreq;          /* lowpass frequency (Hz) */
00692    double  BWFreq;          /* bandwidth of filters (Hz) */
00693    double  MeanNoise;       /* mean primary sensor noise (T) */
00694    char    MriName[256];    /* MRI image file name */
00695    int     Nasion[3];       /* MRI voxel index for nasion */
00696    int     RightPA[3];      /* MRI voxel index for right pre-auricular */
00697    int     LeftPA[3];       /* MRI voxel index for left pre-auricular */
00698    int     SAMType;         /* SAM file type */
00699    int     SAMUnit;         /* SAM units (a bit redundant, but may be useful) */
00700    int     pad_bytes2;      /* ** align end of structure on 8 byte boundary */
00701    double  MegNasion[3];    /* MEG dewar coordinates for nasion (m) */
00702    double  MegRightPA[3];   /* MEG dewar coordinates for R pre-auricular (m) */
00703    double  MegLeftPA[3];    /* MEG dewar coordinates for L pre-auricular (m) */
00704    char    SAMUnitName[32]; /* SAM units (redundant, but useful too!) */
00705 } SAM_HDR_v2;
00706 #endif
00707 
00708 
00709 /*-------------------------*/
00710 /*! Macro for bad return. */
00711 
00712 #undef  BADBAD
00713 #define BADBAD(s)                                                \
00714   do{ fprintf(stderr,"** THD_open_ctfsam(%s): %s\n",fname,s);    \
00715       RETURN(NULL);                                              \
00716   } while(0)
00717 
00718 /*-----------------------------------------------------------------*/
00719 /*! Open a CTF .svl (SAM) file as an unpopulated AFNI dataset.
00720     It will be populated later, in THD_load_ctfsam().
00721 -------------------------------------------------------------------*/
00722 
00723 THD_3dim_dataset * THD_open_ctfsam( char *fname )
00724 {
00725    FILE *fp ;
00726    SAM_HDR hh ;
00727    char Identity[9] ;
00728    int ii,nn , swap ;
00729    THD_3dim_dataset *dset=NULL ;
00730    char prefix[THD_MAX_PREFIX] , *ppp , tname[12] , ori[4] ;
00731    THD_ivec3 nxyz , orixyz ;
00732    THD_fvec3 dxyz , orgxyz ;
00733    int iview ;
00734    int ngood , length , datum_type , datum_len , oxx,oyy,ozz ;
00735    int   nx,ny,nz ;
00736    float dx,dy,dz , xorg,yorg,zorg ;
00737 
00738 ENTRY("THD_open_ctfsam") ;
00739 
00740    /* open input file */
00741 
00742    if( fname == NULL || *fname == '\0' ) BADBAD("bad input filename");
00743    fp = fopen( fname , "rb" ) ;
00744    if( fp == NULL )                      BADBAD("can't open input file");
00745 
00746    /* read header [1st 8 bytes are "SAMIMAGE"] */
00747 
00748    fread( Identity , 1,8 , fp ) ; Identity[8] = '\0' ;
00749    fread( &hh , sizeof(hh) , 1 , fp ) ;
00750    fclose(fp) ;
00751 
00752    if( strcmp(Identity,"SAMIMAGE") != 0 ) BADBAD("Identity != SAMIMAGE") ;
00753    if( hh.Version                  == 0 ) BADBAD("bad header Version") ;
00754 
00755    swap = (hh.Version < 0) || (hh.Version > 3) ;    /* byte swap? */
00756 
00757    if( swap ){                   /* swap various header fields */
00758      swap_4( &hh.Version    ) ;
00759      swap_4( &hh.NumChans   ) ;
00760      swap_4( &hh.NumWeights ) ;
00761      swap_8( &hh.XStart     ) ;
00762      swap_8( &hh.XEnd       ) ;
00763      swap_8( &hh.YStart     ) ;
00764      swap_8( &hh.YEnd       ) ;
00765      swap_8( &hh.ZStart     ) ;
00766      swap_8( &hh.ZEnd       ) ;
00767      swap_8( &hh.StepSize   ) ;
00768      swap_8( &hh.HPFreq     ) ;
00769      swap_8( &hh.LPFreq     ) ;
00770      swap_8( &hh.BWFreq     ) ;
00771      swap_8( &hh.MeanNoise  ) ;
00772      swap_4( &hh.Nasion[0]  ) ;
00773      swap_4( &hh.RightPA[0] ) ;
00774      swap_4( &hh.LeftPA[0]  ) ;
00775      swap_4( &hh.Nasion[1]  ) ;
00776      swap_4( &hh.RightPA[1] ) ;
00777      swap_4( &hh.LeftPA[1]  ) ;
00778      swap_4( &hh.Nasion[2]  ) ;
00779      swap_4( &hh.RightPA[2] ) ;
00780      swap_4( &hh.LeftPA[2]  ) ;
00781      swap_4( &hh.SAMType    ) ;
00782      swap_4( &hh.SAMUnit    ) ;
00783    }
00784 
00785    /* simple checks on header values */
00786 
00787    if( hh.Version  < 0        ||
00788        hh.Version  > 3        ||  /* 26 Feb 2005 */
00789        hh.XStart   >= hh.XEnd ||
00790        hh.YStart   >= hh.YEnd ||
00791        hh.ZStart   >= hh.ZEnd ||
00792        hh.StepSize <= 0.0       ) BADBAD("bad header data") ;
00793 
00794 #if 0
00795    printf("\n") ;
00796    printf("**CTF SAM : %s\n",fname) ;
00797    printf("Version   = %d\n",hh.Version) ;
00798    printf("NumChans  = %d\n",hh.NumChans) ;
00799    printf("NumWeights= %d\n",hh.NumWeights) ;
00800    printf("XStart    = %g\n",hh.XStart) ;
00801    printf("Xend      = %g\n",hh.XEnd) ;
00802    printf("YStart    = %g\n",hh.YStart) ;
00803    printf("YEnd      = %g\n",hh.YEnd) ;
00804    printf("ZStart    = %g\n",hh.ZStart) ;
00805    printf("Zend      = %g\n",hh.ZEnd) ;
00806    printf("StepSize  = %g\n",hh.StepSize) ;
00807    printf("HPFreq    = %g\n",hh.HPFreq) ;
00808    printf("LPFreq    = %g\n",hh.LPFreq) ;
00809    printf("BWFreq    = %g\n",hh.BWFreq) ;
00810    printf("MeanNoise = %g\n",hh.MeanNoise) ;
00811    printf("Nasion[0] = %d\n",hh.Nasion[0]) ;
00812    printf("Nasion[1] = %d\n",hh.Nasion[1]) ;
00813    printf("Nasion[2] = %d\n",hh.Nasion[2]) ;
00814    printf("RightPA[0]= %d\n",hh.RightPA[0]) ;
00815    printf("RightPA[1]= %d\n",hh.RightPA[1]) ;
00816    printf("RightPA[2]= %d\n",hh.RightPA[2]) ;
00817    printf("LeftPA[0] = %d\n",hh.LeftPA[0]) ;
00818    printf("LeftPA[1] = %d\n",hh.LeftPA[1]) ;
00819    printf("LeftPA[2] = %d\n",hh.LeftPA[2]) ;
00820    printf("SAMtype   = %d\n",hh.SAMType) ;
00821    printf("SAMunit   = %d\n",hh.SAMUnit) ;
00822    printf("SetName   = %s\n",hh.SetName) ;
00823    printf("MriName   = %s\n",hh.MriName) ;
00824    printf("headersize= %d\n",sizeof(hh)+8) ;
00825 #endif
00826 
00827    hh.StepSize *= 1000.0 ;   /* convert distances from m to mm */
00828    hh.XStart   *= 1000.0 ;   /* (who the hell uses meters for brain imaging?) */
00829    hh.YStart   *= 1000.0 ;   /* (blue whales?  elephants?) */
00830    hh.ZStart   *= 1000.0 ;
00831    hh.XEnd     *= 1000.0 ;
00832    hh.YEnd     *= 1000.0 ;
00833    hh.ZEnd     *= 1000.0 ;
00834 
00835    dx = dy = dz = hh.StepSize ;  /* will be altered below */
00836 
00837 #if 0
00838    nx = (int)((hh.ZEnd - hh.ZStart)/dz + 0.99999); /* dataset is stored in Z,Y,X order */
00839    ny = (int)((hh.YEnd - hh.YStart)/dy + 0.99999); /* but AFNI calls these x,y,z       */
00840    nz = (int)((hh.XEnd - hh.XStart)/dx + 0.99999);
00841 #else
00842    nx = CTF_count( hh.ZStart , hh.ZEnd , hh.StepSize ) ;
00843    ny = CTF_count( hh.YStart , hh.YEnd , hh.StepSize ) ;
00844    nz = CTF_count( hh.XStart , hh.XEnd , hh.StepSize ) ;
00845 #endif
00846 
00847    /* determine if file is big enough to hold all data it claims */
00848 
00849    nn = THD_filesize(fname) ;
00850    if( nn < sizeof(double)*nx*ny*nz ) BADBAD("input file too small") ;
00851 
00852    datum_type = MRI_float ;  /* actually is double, but AFNI doesn't grok that */
00853                              /* will be converted to floats when reading data */
00854 
00855    /* set orientation = IRP = xyz ordering */
00856 
00857    ori[0] = 'I'; ori[1] = 'R'; ori[2] = 'P';
00858 
00859    oxx = ORCODE(ori[0]); oyy = ORCODE(ori[1]); ozz = ORCODE(ori[2]);
00860    if( !OR3OK(oxx,oyy,ozz) ){
00861      oxx = ORI_I2S_TYPE; oyy = ORI_R2L_TYPE; ozz = ORI_P2A_TYPE;   /** IRP? **/
00862    }
00863 
00864    orixyz.ijk[0] = oxx ; orixyz.ijk[1] = oyy ; orixyz.ijk[2] = ozz ;
00865 
00866    /* now set grid size, keeping in mind that
00867       A-P is positive and P-A is negative,
00868       R-L is positive and L-R is negative,
00869       I-S is positive and S-I is negative.       */
00870 
00871    switch( ori[0] ){
00872      case 'A':  dx =  hh.StepSize ; xorg = -hh.XStart ; break ;
00873      case 'P':  dx = -hh.StepSize ; xorg = -hh.XStart ; break ;
00874      case 'R':  dx =  hh.StepSize ; xorg =  hh.YStart ; break ;
00875      case 'L':  dx = -hh.StepSize ; xorg =  hh.YStart ; break ;
00876      case 'I':  dx =  hh.StepSize ; xorg =  hh.ZStart ; break ;
00877      case 'S':  dx = -hh.StepSize ; xorg =  hh.ZStart ; break ;
00878    }
00879    switch( ori[1] ){
00880      case 'A':  dy =  hh.StepSize ; yorg = -hh.XStart ; break ;
00881      case 'P':  dy = -hh.StepSize ; yorg = -hh.XStart ; break ;
00882      case 'R':  dy =  hh.StepSize ; yorg =  hh.YStart ; break ;
00883      case 'L':  dy = -hh.StepSize ; yorg =  hh.YStart ; break ;
00884      case 'I':  dy =  hh.StepSize ; yorg =  hh.ZStart ; break ;
00885      case 'S':  dy = -hh.StepSize ; yorg =  hh.ZStart ; break ;
00886    }
00887    switch( ori[2] ){
00888      case 'A':  dz =  hh.StepSize ; zorg = -hh.XStart ; break ;
00889      case 'P':  dz = -hh.StepSize ; zorg = -hh.XStart ; break ;
00890      case 'R':  dz =  hh.StepSize ; zorg =  hh.YStart ; break ;
00891      case 'L':  dz = -hh.StepSize ; zorg =  hh.YStart ; break ;
00892      case 'I':  dz =  hh.StepSize ; zorg =  hh.ZStart ; break ;
00893      case 'S':  dz = -hh.StepSize ; zorg =  hh.ZStart ; break ;
00894    }
00895 
00896    /*-- make a dataset --*/
00897 
00898    dset = EDIT_empty_copy(NULL) ;
00899 
00900    dset->idcode.str[0] = 'C' ;  /* overwrite 1st 3 bytes */
00901    dset->idcode.str[1] = 'T' ;
00902    dset->idcode.str[2] = 'F' ;
00903 
00904    MCW_hash_idcode( fname , dset ) ;  /* 06 May 2005 */
00905 
00906    ppp = THD_trailname(fname,0) ;                   /* strip directory */
00907    MCW_strncpy( prefix , ppp , THD_MAX_PREFIX ) ;   /* to make prefix */
00908 
00909    nxyz.ijk[0] = nx ; dxyz.xyz[0] = dx ;  /* setup axes lengths and voxel sizes */
00910    nxyz.ijk[1] = ny ; dxyz.xyz[1] = dy ;
00911    nxyz.ijk[2] = nz ; dxyz.xyz[2] = dz ;
00912 
00913    orixyz.ijk[0] = oxx ; orgxyz.xyz[0] = xorg ;
00914    orixyz.ijk[1] = oyy ; orgxyz.xyz[1] = yorg ;
00915    orixyz.ijk[2] = ozz ; orgxyz.xyz[2] = zorg ;
00916 
00917    iview = VIEW_ORIGINAL_TYPE ;
00918 
00919    /*-- actually send the values above into the dataset header --*/
00920 
00921    EDIT_dset_items( dset ,
00922                       ADN_prefix      , prefix ,
00923                       ADN_datum_all   , datum_type ,
00924                       ADN_nxyz        , nxyz ,
00925                       ADN_xyzdel      , dxyz ,
00926                       ADN_xyzorg      , orgxyz ,
00927                       ADN_xyzorient   , orixyz ,
00928                       ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00929                       ADN_nvals       , 1 ,
00930                       ADN_type        , HEAD_FUNC_TYPE ,
00931                       ADN_view_type   , iview ,
00932                       ADN_func_type   , FUNC_FIM_TYPE ,
00933                     ADN_none ) ;
00934 
00935    /*-- set byte order (for reading from disk) --*/
00936 
00937    ii = mri_short_order() ;
00938    if( swap ) ii = REVERSE_ORDER(ii) ;
00939    dset->dblk->diskptr->byte_order = ii ;
00940 
00941    /*-- flag to read data from disk using CTF SAM mode --*/
00942 
00943    dset->dblk->diskptr->storage_mode = STORAGE_BY_CTFSAM ;
00944    strcpy( dset->dblk->diskptr->brick_name , fname ) ;
00945 
00946    RETURN(dset) ;
00947 }
00948 
00949 /*------------------------------------------------------------------*/
00950 /*! Actually load data from a CTF SAM file into a dataset.
00951 --------------------------------------------------------------------*/
00952 
00953 void THD_load_ctfsam( THD_datablock *dblk )
00954 {
00955    THD_diskptr *dkptr ;
00956    int nx,ny,nz,nv , nxy,nxyz,nxyzv , ibr,nbad , ii,swap ;
00957    FILE *fp ;
00958    void *ptr ;
00959    double *dbar ;
00960    float  *ftar ;
00961 
00962 ENTRY("THD_load_ctfsam") ;
00963 
00964    /*-- check inputs --*/
00965 
00966    if( !ISVALID_DATABLOCK(dblk)                         ||
00967        dblk->diskptr->storage_mode != STORAGE_BY_CTFSAM ||
00968        dblk->brick == NULL                                ) EXRETURN ;
00969 
00970    dkptr = dblk->diskptr ;
00971 
00972    /*-- allocate space for data --*/
00973 
00974    nx = dkptr->dimsizes[0] ;
00975    ny = dkptr->dimsizes[1] ; nxy   = nx * ny   ;
00976    nz = dkptr->dimsizes[2] ; nxyz  = nxy * nz  ;
00977    nv = dkptr->nvals       ; nxyzv = nxyz * nv ;
00978 
00979    /* position file 8*nxyzv bytes before end of file */
00980 
00981    fp = fopen( dkptr->brick_name , "rb" ) ;  /* .svl file */
00982    if( fp == NULL ) EXRETURN ;
00983 
00984    /* 26 Feb 2005: instead of skipping the header,
00985                    whose size varies with the SAM version number,
00986                    just seek backwards from the end to the correct size */
00987 
00988    fseek( fp , -sizeof(double)*nxyzv , SEEK_END ) ;
00989 
00990    dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
00991 
00992    /*-- malloc space for each brick separately --*/
00993 
00994    for( nbad=ibr=0 ; ibr < nv ; ibr++ ){
00995      if( DBLK_ARRAY(dblk,ibr) == NULL ){
00996        ptr = AFMALL(void, DBLK_BRICK_BYTES(dblk,ibr) ) ;
00997        mri_fix_data_pointer( ptr ,  DBLK_BRICK(dblk,ibr) ) ;
00998        if( ptr == NULL ) nbad++ ;
00999      }
01000    }
01001 
01002    /*-- if couldn't get them all, take our ball and go home in a snit --*/
01003 
01004    if( nbad > 0 ){
01005      fprintf(stderr,
01006              "\n** failed to malloc %d CTR MRI bricks out of %d\n\a",nbad,nv);
01007      for( ibr=0 ; ibr < nv ; ibr++ ){
01008        if( DBLK_ARRAY(dblk,ibr) != NULL ){
01009          free(DBLK_ARRAY(dblk,ibr)) ;
01010          mri_fix_data_pointer( NULL , DBLK_BRICK(dblk,ibr) ) ;
01011        }
01012      }
01013      fclose(fp) ; EXRETURN ;
01014    }
01015 
01016    /*-- SAM data is stored as doubles,
01017         but we have to store it in AFNI as floats --*/
01018 
01019    dbar = (double *) calloc(sizeof(double),nxyz) ;     /* workspace */
01020    swap = ( dkptr->byte_order != mri_short_order() ) ;
01021 
01022    for( ibr=0 ; ibr < nv ; ibr++ ){            /* loop over sub-bricks */
01023      fread( dbar, 1, sizeof(double)*nxyz, fp ) ; /* read data to workspace */
01024      ftar = DBLK_ARRAY(dblk,ibr) ;               /* float array in dataset */
01025      for( ii=0 ; ii < nxyz ; ii++ ){             /* loop over voxels */
01026        if( swap ) swap_8(dbar+ii) ;                /* swap it */
01027        ftar[ii] = dbar[ii] ;                       /* save it as a float */
01028      }
01029    }
01030 
01031    fclose(fp) ; free(dbar) ;  /* toss out the trash */
01032    EXRETURN ;
01033 }