#include "mrilib.h"
#include "mayo_analyze.h"
#include "thd.h"
/*******************************************************************/
/******* 27 Aug 2002: Read an ANALYZE file as an AFNI dataset ******/
/*******************************************************************/
/*---------------------------------------------------------------*/
/*! Swap the 2 bytes pointed to by ppp: ab -> ba. */
static void swap_2(void *ppp)
{
unsigned char *pntr = (unsigned char *) ppp ;
unsigned char b0, b1;
b0 = *pntr; b1 = *(pntr+1);
*pntr = b1; *(pntr+1) = b0;
}
/*---------------------------------------------------------------*/
/*! Swap the 4 bytes pointed to by ppp: abcd -> dcba. */
static void swap_4(void *ppp)
{
unsigned char *pntr = (unsigned char *) ppp ;
unsigned char b0, b1, b2, b3;
b0 = *pntr; b1 = *(pntr+1); b2 = *(pntr+2); b3 = *(pntr+3);
*pntr = b3; *(pntr+1) = b2; *(pntr+2) = b1; *(pntr+3) = b0;
}
#if 0
/*---------------------------------------------------------------*/
/*! Swap the 8 bytes pointed to by ppp: abcdefgh -> hgfedcba. */
static void swap_8(void *ppp)
{
unsigned char *pntr = (unsigned char *) ppp ;
unsigned char b0, b1, b2, b3;
unsigned char b4, b5, b6, b7;
b0 = *pntr ; b1 = *(pntr+1); b2 = *(pntr+2); b3 = *(pntr+3);
b4 = *(pntr+4); b5 = *(pntr+5); b6 = *(pntr+6); b7 = *(pntr+7);
*pntr = b7; *(pntr+1) = b6; *(pntr+2) = b5; *(pntr+3) = b4;
*(pntr+4) = b3; *(pntr+5) = b2; *(pntr+6) = b1; *(pntr+7) = b0;
}
#endif
/*---------------------------------------------------------------*/
/*! Byte swap ANALYZE file header in various places */
static void swap_analyze_hdr( struct dsr *pntr )
{
swap_4(&pntr->hk.sizeof_hdr) ;
swap_4(&pntr->hk.extents) ;
swap_2(&pntr->hk.session_error) ;
swap_2(&pntr->dime.dim[0]) ;
swap_2(&pntr->dime.dim[1]) ;
swap_2(&pntr->dime.dim[2]) ;
swap_2(&pntr->dime.dim[3]) ;
swap_2(&pntr->dime.dim[4]) ;
swap_2(&pntr->dime.dim[5]) ;
swap_2(&pntr->dime.dim[6]) ;
swap_2(&pntr->dime.dim[7]) ;
#if 0
swap_2(&pntr->dime.unused1) ;
#endif
swap_2(&pntr->dime.datatype) ;
swap_2(&pntr->dime.bitpix) ;
swap_4(&pntr->dime.pixdim[0]) ;
swap_4(&pntr->dime.pixdim[1]) ;
swap_4(&pntr->dime.pixdim[2]) ;
swap_4(&pntr->dime.pixdim[3]) ;
swap_4(&pntr->dime.pixdim[4]) ;
swap_4(&pntr->dime.pixdim[5]) ;
swap_4(&pntr->dime.pixdim[6]) ;
swap_4(&pntr->dime.pixdim[7]) ;
swap_4(&pntr->dime.vox_offset) ;
swap_4(&pntr->dime.funused1) ;
swap_4(&pntr->dime.funused2) ;
swap_4(&pntr->dime.cal_max) ;
swap_4(&pntr->dime.cal_min) ;
swap_4(&pntr->dime.compressed) ;
swap_4(&pntr->dime.verified) ;
swap_2(&pntr->dime.dim_un0) ;
swap_4(&pntr->dime.glmax) ;
swap_4(&pntr->dime.glmin) ;
}
/*-----------------------------------------------------------------*/
#ifndef DONT_INCLUDE_ANALYZE_STRUCT
#define DONT_INCLUDE_ANALYZE_STRUCT
#endif
#include "nifti1_io.h" /*** NIFTI spec and funcs ***/
/*-----------------------------------------------------------------*/
#undef swap_2
#undef swap_4
/*-----------------------------------------------------------------*/
/*! Open an ANALYZE .hdr file as an unpopulated AFNI dataset.
It will be populated from the .img file, in THD_load_analyze().
-------------------------------------------------------------------*/
THD_3dim_dataset * THD_open_analyze( char *hname )
{
FILE * fp ;
char iname[THD_MAX_NAME] ;
int ii , doswap ;
struct dsr hdr ; /* ANALYZE .hdr format */
int ngood , length , datum_type , datum_len ;
int nx,ny,nz,nt ;
float dx,dy,dz,dt ;
float fac=0.0 ; /* brick scaling factor */
THD_3dim_dataset *dset=NULL ;
char prefix[THD_MAX_PREFIX] , *ppp ;
THD_ivec3 nxyz , orixyz ;
THD_fvec3 dxyz , orgxyz ;
int iview ;
short spmorg , spmxx,spmyy,spmzz ; /* 03 Nov 2003 */
ENTRY("THD_open_analyze") ;
/* 28 Aug 2003: check if this is a NIFTI file instead */
{ nifti_image *nim = nifti_image_read(hname,0) ;
if( nim != NULL && nim->nifti_type > 0 ){
nifti_image_free(nim); dset = THD_open_nifti(hname); RETURN(dset);
}
}
/*-- check & prepare filenames --*/
if( hname == NULL ) RETURN( NULL );
ii = strlen(hname) ; if( ii < 5 ) RETURN( NULL );
if( strcmp(hname+ii-4,".hdr") != 0 ) RETURN( NULL );
strcpy(iname,hname) ; strcpy(iname+ii-3,"img") ; /* .img filename */
/*-- read .hdr file into struct --*/
fp = fopen( hname , "rb" ) ;
if( fp == NULL ) RETURN( NULL );
hdr.dime.dim[0] = 0 ;
fread( &hdr , 1 , sizeof(struct dsr) , fp ) ;
fclose(fp) ;
if( hdr.dime.dim[0] == 0 ){ /* bad input */
ERROR_message("ANALYZE file %s has dim[0]=0!\n",hname) ;
RETURN( NULL ) ;
}
/*-- check .img file for existence and size --*/
length = THD_filesize(iname) ; /* will use this later */
if( length <= 0 ){
ERROR_message("Can't find ANALYZE file %s\n",iname) ;
RETURN( NULL );
}
/*-- check for swap-age of header --*/
doswap = (hdr.dime.dim[0] < 0 || hdr.dime.dim[0] > 15) ;
if( doswap ) swap_analyze_hdr( &hdr ) ;
/*-- get a scale factor for images --*/
if( !AFNI_noenv("AFNI_ANALYZE_SCALE") ){
fac = hdr.dime.funused1 ;
(void) thd_floatscan( 1 , &fac ) ;
if( fac < 0.0 || fac == 1.0 ) fac = 0.0 ;
}
/*-- get data type for each voxel --*/
switch( hdr.dime.datatype ){
default:
ERROR_message("File %s: Unsupported ANALYZE datatype=%d (%s)\n",
hname,hdr.dime.datatype,ANDT_string(hdr.dime.datatype) ) ;
RETURN( NULL );
case ANDT_UNSIGNED_CHAR: datum_type = MRI_byte ; break;
case ANDT_SIGNED_SHORT: datum_type = MRI_short ; break;
case ANDT_FLOAT: datum_type = MRI_float ; break;
case ANDT_COMPLEX: datum_type = MRI_complex; break;
case ANDT_RGB: datum_type = MRI_rgb ; break;
#if 0
case ANDT_SIGNED_INT: datum_type = MRI_int ; break; /* not in AFNI */
#endif
}
datum_len = mri_datum_size(datum_type) ; /* number bytes per voxel */
/*-- compute dimensions of images, and number of images --*/
nx = hdr.dime.dim[1] ;
ny = hdr.dime.dim[2] ;
if( nx < 2 || ny < 2 ) RETURN( NULL );
switch( hdr.dime.dim[0] ){
case 2: nz = 1 ; nt = 1 ; ; break ;
case 3: nz = hdr.dime.dim[3] ; nt = 1 ; ; break ;
case 4: nz = hdr.dime.dim[3] ; nt = hdr.dime.dim[4] ; break ;
default:
ERROR_message("ANALYZE file %s has %d dimensions!\n",
hname,hdr.dime.dim[0]) ;
RETURN( NULL ) ;
}
if( nz < 1 ) nz = 1 ;
if( nt < 1 ) nt = 1 ;
/*-- voxel sizes --*/
dx = fabs(hdr.dime.pixdim[1]) ; if( dx == 0.0 ) dx = 1.0 ;
dy = fabs(hdr.dime.pixdim[2]) ; if( dy == 0.0 ) dy = 1.0 ;
dz = fabs(hdr.dime.pixdim[3]) ; if( dz == 0.0 ) dz = 1.0 ;
dt = fabs(hdr.dime.pixdim[4]) ; if( dt == 0.0 || nt == 1 ) dt = 1.0 ;
ngood = datum_len*nx*ny*nz*nt ; /* number bytes needed in .img file */
if( length < ngood ){
ERROR_message(
"ANALYZE file %s is %d bytes long but must be %d bytes long\n"
"** for nx=%d ny=%d nz=%d nt=%d and %d bytes/voxel\n",
iname,length,ngood,nx,ny,nz,nt,datum_len ) ;
fclose(fp) ; RETURN( NULL );
}
/* check SPM originator field - 03 Nov 2003 */
/* 11 Mar 2004 - oops: SPM indexes start at 1 - RWCox */
spmorg = spmxx = spmyy = spmzz = 0 ;
{ short xyzuv[5] , xx,yy,zz,aa,bb ;
memcpy( xyzuv , hdr.hist.originator , 10 ) ;
if( doswap ){ swap_2(&(xyzuv[0])); swap_2(&(xyzuv[1])); swap_2(&(xyzuv[2])); swap_2(&(xyzuv[3])); swap_2(&(xyzuv[4]));}
/*fprintf(stderr, "\n"
"hdr.hist.originator: %d %d %d %d %d\n",
xyzuv[0], xyzuv[1], xyzuv[2], xyzuv[3], xyzuv[4]
);*/
if( 1 || (xyzuv[3] == 0 && xyzuv[4] == 0) ){ /* ZSS Nov 10 05. Looks like xyzuv[3] and xyzuv[4] can be complete rubbish. Not reliable. */
xx = xyzuv[0] ; yy = xyzuv[1] ; zz = xyzuv[2] ;
if( xx > 0 && xx < nx-1 &&
yy > 0 && yy < ny-1 &&
zz > 0 && zz < nz-1 ){ spmorg=1; spmxx=xx-1; spmyy=yy-1; spmzz=zz-1; }
/*fprintf(stderr, "\n"
" xx yy zz : %d %d %d \n"
"spmorg spmxx spmyy spmzz : %d %d %d %d\n",
xx, yy, zz,
spmorg, spmxx, spmyy, spmzz);*/
}
}
/*-- make a dataset --*/
dset = EDIT_empty_copy(NULL) ;
dset->idcode.str[0] = 'A' ; /* overwrite 1st 3 bytes */
dset->idcode.str[1] = 'N' ;
dset->idcode.str[2] = 'Z' ;
MCW_hash_idcode( hname , dset ) ; /* 06 May 2005 */
ppp = THD_trailname(hname,0) ; /* strip directory */
MCW_strncpy( prefix , ppp , THD_MAX_PREFIX ) ; /* to make prefix */
nxyz.ijk[0] = nx ; dxyz.xyz[0] = dx ; /* setup axes lengths and voxel sizes */
nxyz.ijk[1] = ny ; dxyz.xyz[1] = dy ;
nxyz.ijk[2] = nz ; dxyz.xyz[2] = dz ;
/*-- default orientation is LPI, but user can alter via environment --*/
{ char *ori = getenv("AFNI_ANALYZE_ORIENT") ;
int oxx,oyy,ozz ;
if( ori == NULL || strlen(ori) < 3 ) {
static int nwarn=0 ;
ori = "LPI"; /* set default LPI */
if( nwarn < 2 ){
WARNING_message("Assuming ANALYZE orientaion is LPI.\n"
"++ To change orientation or silence this message,\n"
"++ Set AFNI_ANALYZE_ORIENT to the proper orientation\n"
"++ in your .afnirc file.\n"
"++ e.g.: AFNI_ANALYZE_ORIENT = LPI");
nwarn++ ;
}
}
oxx = ORCODE(ori[0]); oyy = ORCODE(ori[1]); ozz = ORCODE(ori[2]);
if( !OR3OK(oxx,oyy,ozz) ){
oxx = ORI_L2R_TYPE; oyy = ORI_P2A_TYPE; ozz = ORI_I2S_TYPE; /* LPI? */
}
orixyz.ijk[0] = oxx ;
orixyz.ijk[1] = oyy ;
orixyz.ijk[2] = ozz ;
}
/*-- origin of coordinates --*/
orgxyz.xyz[0] = 0.5*dx ; /* FSL: (0,0,0) is at outer corner of 1st voxel */
orgxyz.xyz[1] = 0.5*dy ; /* AFNI: these are coords of center of 1st voxel */
orgxyz.xyz[2] = 0.5*dz ;
/*-- 04 Oct 2002: allow auto-centering of ANALYZE datasets --*/
if( AFNI_yesenv("AFNI_ANALYZE_AUTOCENTER") ){
static int nwarn=0 ;
if( nwarn < 2 )
WARNING_message("Autocentering datasets because"
" AFNI_ANALYZE_AUTOCENTER is set");
orgxyz.xyz[0] = -0.5 * (nx-1) * dx ;
orgxyz.xyz[1] = -0.5 * (ny-1) * dy ;
orgxyz.xyz[2] = -0.5 * (nz-1) * dz ; nwarn++ ;
}
iview = VIEW_ORIGINAL_TYPE ; /* can't tell if it is Talairach-ed (default)*/
{/* ZSS Dec 15 03 */
char *vie = getenv("AFNI_ANALYZE_VIEW") ;
if (!vie) {
static int nwarn = 0 ;
iview = VIEW_ORIGINAL_TYPE;
if( nwarn < 2 ){
WARNING_message("Assuming view is +orig.\n"
"++ To change view or silence this message,\n"
"++ Set AFNI_ANALYZE_VIEW to the proper view\n"
"++ in your .afnirc file.\n"
"++ e.g.: AFNI_ANALYZE_VIEW = orig");
nwarn++ ;
}
} else {
static int nwarn = 0 ;
if (strcmp(vie, "tlrc") == 0) iview = VIEW_TALAIRACH_TYPE;
else if (strcmp(vie, "orig") == 0) iview = VIEW_ORIGINAL_TYPE;
else if (strcmp(vie, "acpc") == 0) iview = VIEW_ACPCALIGNED_TYPE;
else if( nwarn < 2 ) {
WARNING_message("Bad value (%s) for environment \n"
"++ variable AFNI_ANALYZE_VIEW. Choose from:\n"
"++ orig or acpc or tlrc.\n"
"++ Assuming orig view.\n", vie);
iview = VIEW_ORIGINAL_TYPE; nwarn++ ;
}
}
}
if( AFNI_yesenv("AFNI_ANALYZE_ORIGINATOR") && spmorg ){ /* 03 Nov 2003 */
if ( !getenv ("AFNI_ANALYZE_VIEW") ) { /* ZSS 16 Dec 2003 */
iview = VIEW_TALAIRACH_TYPE ; /* for backward compatibility */
}
orgxyz.xyz[0] = -spmxx * dx ; /* (0,0,0) is at (spmxx,spmyy,spmzz) */
orgxyz.xyz[1] = -spmyy * dy ;
orgxyz.xyz[2] = -spmzz * dz ;
} else {
if (!spmorg){
static int nwarn=0 ;
if( nwarn < 2 )
WARNING_message("No ANALYZE origin found in file %s\n",hname);
nwarn++ ;
} else {
static int nwarn=0 ;
if( nwarn < 2 )
WARNING_message("ANALYZE origin ignored in file %s\n"
"++ If datasets are out of alignment,\n"
"++ Set AFNI_ANALYZE_ORIGINATOR = YES\n"
"++ in your .afnirc file.\n" , hname );
nwarn++ ;
}
}
#if 0
fprintf (stderr, "\n"
"orgxyz.xyz: %f %f %f\n",
orgxyz.xyz[0], orgxyz.xyz[1], orgxyz.xyz[2]);
#endif
/* 10 Oct 2002: change voxel size signs, if axis orientation is negative */
/* [above, we assumed that axes were oriented in - to + way] */
if( ORIENT_sign[orixyz.ijk[0]] == '-' ){
dxyz.xyz[0] = -dxyz.xyz[0] ;
orgxyz.xyz[0] = -orgxyz.xyz[0] ;
}
if( ORIENT_sign[orixyz.ijk[1]] == '-' ){
dxyz.xyz[1] = -dxyz.xyz[1] ;
orgxyz.xyz[1] = -orgxyz.xyz[1] ;
}
if( ORIENT_sign[orixyz.ijk[2]] == '-' ){
dxyz.xyz[2] = -dxyz.xyz[2] ;
orgxyz.xyz[2] = -orgxyz.xyz[2] ;
}
/*-- actually send the values above into the dataset header --*/
EDIT_dset_items( dset ,
ADN_prefix , prefix ,
ADN_datum_all , datum_type ,
ADN_nxyz , nxyz ,
ADN_xyzdel , dxyz ,
ADN_xyzorg , orgxyz ,
ADN_xyzorient , orixyz ,
ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
ADN_nvals , nt ,
ADN_type , HEAD_ANAT_TYPE ,
ADN_view_type , iview ,
ADN_func_type , ANAT_MRAN_TYPE ,
ADN_none ) ;
if( nt > 1 ) /** pretend it is 3D+time **/
EDIT_dset_items( dset ,
ADN_func_type, ANAT_EPI_TYPE ,
ADN_ntt , nt ,
ADN_ttorg , 0.0 ,
ADN_ttdel , dt ,
ADN_ttdur , 0.0 ,
ADN_tunits , UNITS_SEC_TYPE ,
ADN_none ) ;
/* make it a func if is is a FSL special name */
#ifdef ALLOW_FSL_FEAT
if( strstr(prefix,"stat") != NULL )
EDIT_dset_items( dset ,
ADN_type , HEAD_FUNC_TYPE ,
ADN_func_type , FUNC_FIM_TYPE ,
ADN_none ) ;
#endif
/*-- set brick factors (all the same) --*/
if( fac > 0.0 ){
float *bf = AFMALL(float, sizeof(float)*nt) ;
for( ii=0 ; ii < nt ; ii++ ) bf[ii] = fac ;
EDIT_dset_items( dset , ADN_brick_fac,bf , ADN_none ) ;
free(bf) ;
}
/*-- set byte order (for reading from disk) --*/
ii = mri_short_order() ;
if( doswap ) ii = REVERSE_ORDER(ii) ;
dset->dblk->diskptr->byte_order = ii ;
/*-- flag to read data from disk using ANALYZE mode --*/
dset->dblk->diskptr->storage_mode = STORAGE_BY_ANALYZE ;
strcpy( dset->dblk->diskptr->brick_name , iname ) ;
RETURN(dset) ;
}
/*---------------------------------------------------------------------*/
/*! Load an ANALYZE dataset from disk.
-----------------------------------------------------------------------*/
void THD_load_analyze( THD_datablock *dblk )
{
THD_diskptr *dkptr ;
int nx,ny,nz,nv , nxy,nxyz,nxyzv , ibr,nbad ;
FILE *fp ;
void *ptr ;
ENTRY("THD_load_analyze") ;
/*-- check inputs --*/
if( !ISVALID_DATABLOCK(dblk) ||
dblk->diskptr->storage_mode != STORAGE_BY_ANALYZE ||
dblk->brick == NULL ) EXRETURN ;
dkptr = dblk->diskptr ;
fp = fopen( dkptr->brick_name , "rb" ) ; /* .img file */
if( fp == NULL ) EXRETURN ;
/*-- allocate space for data --*/
nx = dkptr->dimsizes[0] ;
ny = dkptr->dimsizes[1] ; nxy = nx * ny ;
nz = dkptr->dimsizes[2] ; nxyz = nxy * nz ;
nv = dkptr->nvals ; nxyzv = nxyz * nv ;
dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
/*-- malloc space for each brick separately --*/
for( nbad=ibr=0 ; ibr < nv ; ibr++ ){
if( DBLK_ARRAY(dblk,ibr) == NULL ){
ptr = AFMALL(void, DBLK_BRICK_BYTES(dblk,ibr) ) ;
mri_fix_data_pointer( ptr , DBLK_BRICK(dblk,ibr) ) ;
if( ptr == NULL ) nbad++ ;
}
}
/*-- if couldn't get them all, take our ball and go home in a snit --*/
if( nbad > 0 ){
ERROR_message("failed to malloc %d ANALYZE bricks out of %d\n\a",nbad,nv);
for( ibr=0 ; ibr < nv ; ibr++ ){
if( DBLK_ARRAY(dblk,ibr) != NULL ){
free(DBLK_ARRAY(dblk,ibr)) ;
mri_fix_data_pointer( NULL , DBLK_BRICK(dblk,ibr) ) ;
}
}
fclose(fp) ; EXRETURN ;
}
/*-- read data from .img file into sub-brick arrays! --*/
for( ibr=0 ; ibr < nv ; ibr++ )
fread( DBLK_ARRAY(dblk,ibr), 1, DBLK_BRICK_BYTES(dblk,ibr), fp ) ;
fclose(fp) ;
/*-- swap bytes? --*/
if( dkptr->byte_order != mri_short_order() ){
for( ibr=0 ; ibr < nv ; ibr++ ){
switch( DBLK_BRICK_TYPE(dblk,ibr) ){
case MRI_short:
mri_swap2( DBLK_BRICK_NVOX(dblk,ibr) , DBLK_ARRAY(dblk,ibr) ) ;
break ;
case MRI_complex:
mri_swap4( 2*DBLK_BRICK_NVOX(dblk,ibr), DBLK_ARRAY(dblk,ibr)) ;
break ;
case MRI_float:
case MRI_int:
mri_swap4( DBLK_BRICK_NVOX(dblk,ibr) , DBLK_ARRAY(dblk,ibr) ) ;
break ;
}
}
}
/*-- check floats --*/
for( ibr=0 ; ibr < nv ; ibr++ ){
if( DBLK_BRICK_TYPE(dblk,ibr) == MRI_float ){
nbad += thd_floatscan( DBLK_BRICK_NVOX(dblk,ibr) ,
DBLK_ARRAY(dblk,ibr) ) ;
} else if( DBLK_BRICK_TYPE(dblk,ibr) == MRI_complex ) {
nbad += thd_complexscan( DBLK_BRICK_NVOX(dblk,ibr) ,
DBLK_ARRAY(dblk,ibr) ) ;
}
}
if( nbad > 0 )
ERROR_message("File %s: found %d float errors -- see program float_scan\n" ,
dkptr->brick_name , nbad ) ;
EXRETURN ;
}
syntax highlighted by Code2HTML, v. 0.9.1