Doxygen Source Code Documentation
thd_niftiread.c File Reference
#include "mrilib.h"#include "nifti1_io.h"Go to the source code of this file.
Defines | |
| #define | MAXNUM(a, b) ( (a) > (b) ? (a):(b)) |
| #define | MAX3(a, b, c) ( (MAXNUM(a,b)) > (MAXNUM(a,c)) ? (MAXNUM(a,b)):(MAXNUM(a,c))) |
| #define | MINNUM(a, b) ( (a) < (b) ? (a):(b)) |
| #define | MIN3(a, b, c) ( (MINNUM(a,b)) < (MINNUM(a,c)) ? (MINNUM(a,b)):(MINNUM(a,c))) |
| #define | CPF(ityp) |
Functions | |
| THD_3dim_dataset * | THD_open_nifti (char *pathname) |
| void | THD_load_nifti (THD_datablock *dblk) |
Define Documentation
|
|
Value: |
|
|
|
|
|
|
|
|
|
|
|
|
Function Documentation
|
|
10 May 2005: see if there is an AFNI extension; if so, load attributes from it and then edit the dataset appropriately * Definition at line 591 of file thd_niftiread.c. References AFMALL, AFNI_yesenv(), THD_datablock::brick, THD_diskptr::brick_name, DATABLOCK_MEM_MALLOC, DBLK_ARRAY, DBLK_BRICK, DBLK_BRICK_BYTES, DBLK_BRICK_NVOX, DBLK_BRICK_TYPE, THD_diskptr::dimsizes, THD_datablock::diskptr, ENTRY, far, free, complex::i, ISVALID_DATABLOCK, THD_datablock::malloc_type, mri_fix_data_pointer(), THD_diskptr::nvals, nz, complex::r, STATUS, STORAGE_BY_NIFTI, THD_diskptr::storage_mode, thd_complexscan(), and thd_floatscan(). Referenced by THD_load_datablock().
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 }
|
|
|
will include nifti1.h * Definition at line 12 of file thd_niftiread.c. References ADN_datum_all, ADN_dz_sl, ADN_func_type, ADN_malloc_type, ADN_none, ADN_nsl, ADN_ntt, ADN_nvals, ADN_nxyz, ADN_prefix, ADN_toff_sl, ADN_ttdel, ADN_ttdur, ADN_ttorg, ADN_tunits, ADN_type, ADN_view_type, ADN_xyzdel, ADN_xyzorg, ADN_xyzorient, ADN_zorg_sl, ANAT_BUCK_TYPE, ANAT_EPI_TYPE, THD_diskptr::brick_name, THD_diskptr::byte_order, calloc, DATABLOCK_MEM_MALLOC, THD_3dim_dataset::dblk, THD_datablock::diskptr, EDIT_BRICK_LABEL, EDIT_dset_items(), EDIT_empty_copy(), EDIT_STATAUX4, ENTRY, free, FUNC_BUCK_TYPE, FUNC_FIM_TYPE, HEAD_ANAT_TYPE, HEAD_FUNC_TYPE, THD_3dim_dataset::idcode, THD_ivec3::ijk, LOAD_FVEC3, LOAD_IVEC3, LOAD_MAT, MAX, MCW_hash_idcode(), MCW_IDSIZE, MCW_strncpy, my_getenv(), NI_group::name, NI_element_type(), NI_free, NI_free_element(), NI_get_attribute(), NI_GROUP_TYPE, NI_read_element(), NI_search_group_deep(), NI_stream_close(), NI_stream_open(), NI_stream_setbuf(), RETURN, STORAGE_BY_NIFTI, THD_diskptr::storage_mode, MCW_idcode::str, THD_datablock_apply_atr(), THD_dblkatr_from_niml(), THD_matrix_to_orientation(), THD_MAX_PREFIX, THD_trailname(), UNITS_SEC_TYPE, VIEW_ORIGINAL_TYPE, VIEW_TALAIRACH_TYPE, THD_fvec3::xyz, and zmax. Referenced by THD_init_session(), THD_open_analyze(), and THD_open_one_dataset().
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 }
|
