* { nifti_brick_list NB_orig, NB_select;
* nifti_image * nim_orig, * nim_select;
* int blist[5] = { 7, 0, 5, 5, 9 };
*
* nim_orig = nifti_image_read_bricks("myfile.nii", 0, NULL, &NB_orig);
* nim_select = nifti_image_read_bricks("myfile.nii", 5, blist, &NB_select);
* }
*
*
* Here, nim_orig gets the entire dataset, where NB_orig.nbricks = 12. But
* nim_select has NB_select.nbricks = 5.
*
* Note that the first case is not quite the same as just calling the
* nifti_image_read function, as here the data is separated into sub-bricks.
*
* Note that valid blist elements are in [0..nt*nu*nv*nw-1],
* or written [ 0 .. (dim[4]*dim[5]*dim[6]*dim[7] - 1) ].
*
* Note that, as is the case with all of the reading functions, the
* data will be allocated, read in, and properly byte-swapped, if
* necessary.
*
* \sa nifti_image_load_bricks, nifti_free_NBL, valid_nifti_brick_list,
nifti_image_read
*//*----------------------------------------------------------------------*/
nifti_image *nifti_image_read_bricks(const char * hname, int nbricks,
const int * blist, nifti_brick_list * NBL)
{
nifti_image * nim;
if( !hname || !NBL ){
fprintf(stderr,"** nifti_image_read_bricks: bad params (%p,%p)\n",
hname, (void *)NBL);
return NULL;
}
if( blist && nbricks <= 0 ){
fprintf(stderr,"** nifti_image_read_bricks: bad nbricks, %d\n", nbricks);
return NULL;
}
nim = nifti_image_read(hname, 0); /* read header, but not data */
if( !nim ) return NULL; /* errors were already printed */
/* if we fail, free image and return */
if( nifti_image_load_bricks(nim, nbricks, blist, NBL) <= 0 ){
nifti_image_free(nim);
return NULL;
}
if( blist ) update_nifti_image_for_brick_list(nim, nbricks);
return nim;
}
/*----------------------------------------------------------------------
* update_nifti_image_for_brick_list - update nifti_image
*
* When loading a specific brick list, the distinction between
* nt, nu, nv and nw is lost. So put everything in t, and set
* dim[0] = 4.
*----------------------------------------------------------------------*/
static void update_nifti_image_for_brick_list( nifti_image * nim , int nbricks )
{
int ndim;
if( g_opts.debug > 2 ){
fprintf(stderr,"+d updating image dimensions for %d bricks in list\n",
nbricks);
fprintf(stderr," ndim = %d\n",nim->ndim);
fprintf(stderr," nx,ny,nz,nt,nu,nv,nw: (%d,%d,%d,%d,%d,%d,%d)\n",
nim->nx, nim->ny, nim->nz, nim->nt, nim->nu, nim->nv, nim->nw);
}
nim->nt = nbricks;
nim->nu = nim->nv = nim->nw = 1;
nim->dim[4] = nbricks;
nim->dim[5] = nim->dim[6] = nim->dim[7] = 1;
/* compute nvox */
/* do not rely on dimensions above dim[0] 16 Nov 2005 [rickr] */
for( nim->nvox = 1, ndim = 1; ndim <= nim->dim[0]; ndim++ )
nim->nvox *= nim->dim[ndim];
/* update the dimensions to 4 or lower */
for( ndim = 4; (ndim > 1) && (nim->dim[ndim] <= 1); ndim-- )
;
if( g_opts.debug > 2 ){
fprintf(stderr,"+d ndim = %d -> %d\n",nim->ndim, ndim);
fprintf(stderr," --> (%d,%d,%d,%d,%d,%d,%d)\n",
nim->nx, nim->ny, nim->nz, nim->nt, nim->nu, nim->nv, nim->nw);
}
nim->dim[0] = nim->ndim = ndim;
}
/*----------------------------------------------------------------------*/
/*! nifti_update_dims_from_array - update nx, ny, ... from nim->dim[]
Fix all the dimension information, based on a new nim->dim[].
Note: we assume that dim[0] will not increase.
Check for updates to pixdim[], dx,..., nx,..., nvox, ndim, dim[0].
*//*--------------------------------------------------------------------*/
int nifti_update_dims_from_array( nifti_image * nim )
{
int c, ndim;
if( !nim ){
fprintf(stderr,"** update_dims: missing nim\n");
return 1;
}
if( g_opts.debug > 2 ){
fprintf(stderr,"+d updating image dimensions given nim->dim:");
for( c = 0; c < 8; c++ ) fprintf(stderr," %d", nim->dim[c]);
fputc('\n',stderr);
}
/* verify dim[0] first */
if(nim->dim[0] < 1 || nim->dim[0] > 7){
fprintf(stderr,"** invalid dim[0], dim[] = ");
for( c = 0; c < 8; c++ ) fprintf(stderr," %d", nim->dim[c]);
fputc('\n',stderr);
return 1;
}
/* set nx, ny ..., dx, dy, ..., one by one */
/* less than 1, set to 1, else copy */
if(nim->dim[1] < 1) nim->nx = nim->dim[1] = 1;
else nim->nx = nim->dim[1];
nim->dx = nim->pixdim[1];
/* if undefined, or less than 1, set to 1 */
if(nim->dim[0] < 2 || (nim->dim[0] >= 2 && nim->dim[2] < 1))
nim->ny = nim->dim[2] = 1;
else
nim->ny = nim->dim[2];
/* copy delta values, in any case */
nim->dy = nim->pixdim[2];
if(nim->dim[0] < 3 || (nim->dim[0] >= 3 && nim->dim[3] < 1))
nim->nz = nim->dim[3] = 1;
else /* just copy vals from arrays */
nim->nz = nim->dim[3];
nim->dz = nim->pixdim[3];
if(nim->dim[0] < 4 || (nim->dim[0] >= 4 && nim->dim[4] < 1))
nim->nt = nim->dim[4] = 1;
else /* just copy vals from arrays */
nim->nt = nim->dim[4];
nim->dt = nim->pixdim[4];
if(nim->dim[0] < 5 || (nim->dim[0] >= 5 && nim->dim[5] < 1))
nim->nu = nim->dim[5] = 1;
else /* just copy vals from arrays */
nim->nu = nim->dim[5];
nim->du = nim->pixdim[5];
if(nim->dim[0] < 6 || (nim->dim[0] >= 6 && nim->dim[6] < 1))
nim->nv = nim->dim[6] = 1;
else /* just copy vals from arrays */
nim->nv = nim->dim[6];
nim->dv = nim->pixdim[6];
if(nim->dim[0] < 7 || (nim->dim[0] >= 7 && nim->dim[7] < 1))
nim->nw = nim->dim[7] = 1;
else /* just copy vals from arrays */
nim->nw = nim->dim[7];
nim->dw = nim->pixdim[7];
for( c = 1, nim->nvox = 1; c <= nim->dim[0]; c++ )
nim->nvox *= nim->dim[c];
/* compute ndim, assuming it can be no larger than the old one */
for( ndim = nim->dim[0]; (ndim > 1) && (nim->dim[ndim] <= 1); ndim-- )
;
if( g_opts.debug > 2 ){
fprintf(stderr,"+d ndim = %d -> %d\n",nim->ndim, ndim);
fprintf(stderr," --> (%d,%d,%d,%d,%d,%d,%d)\n",
nim->nx, nim->ny, nim->nz, nim->nt, nim->nu, nim->nv, nim->nw);
}
nim->dim[0] = nim->ndim = ndim;
return 0;
}
/*----------------------------------------------------------------------*/
/*! Load the image data from disk into an already-prepared image struct.
*
* \param nim - initialized nifti_image, without data
* \param nbricks - the length of blist (must be 0 if blist is NULL)
* \param blist - an array of xyz volume indices to read (can be NULL)
* \param NBL - pointer to struct where resulting data will be stored
*
* If blist is NULL, read all sub-bricks.
*
* \return the number of loaded bricks (NBL->nbricks),
* 0 on failure, < 0 on error
*
* NOTE: it is likely that another function will copy the data pointers
* out of NBL, in which case the only pointer the calling function
* will want to free is NBL->bricks (not each NBL->bricks[i]).
*//*--------------------------------------------------------------------*/
int nifti_image_load_bricks( nifti_image * nim , int nbricks,
const int * blist, nifti_brick_list * NBL )
{
int * slist = NULL, * sindex = NULL, rv;
znzFile fp;
/* we can have blist == NULL */
if( !nim || !NBL ){
fprintf(stderr,"** nifti_image_load_bricks, bad params (%p,%p)\n",
(void *)nim, (void *)NBL);
return -1;
}
if( blist && nbricks <= 0 ){
if( g_opts.debug > 1 )
fprintf(stderr,"-d load_bricks: received blist with nbricks = %d,"
"ignoring blist\n", nbricks);
blist = NULL; /* pretend nothing was passed */
}
if( blist && ! valid_nifti_brick_list(nim, nbricks, blist, g_opts.debug>0) )
return -1;
/* for efficiency, let's read the file in order */
if( blist && nifti_copynsort( nbricks, blist, &slist, &sindex ) != 0 )
return -1;
/* open the file and position the FILE pointer */
fp = nifti_image_load_prep( nim );
if( !fp ){
if( g_opts.debug > 0 )
fprintf(stderr,"** nifti_image_load_bricks, failed load_prep\n");
if( blist ){ free(slist); free(sindex); }
return -1;
}
/* this will flag to allocate defaults */
if( !blist ) nbricks = 0;
if( nifti_alloc_NBL_mem( nim, nbricks, NBL ) != 0 ){
if( blist ){ free(slist); free(sindex); }
znzclose(fp);
return -1;
}
rv = nifti_load_NBL_bricks(nim, slist, sindex, NBL, fp);
if( rv != 0 ){
nifti_free_NBL( NBL ); /* failure! */
NBL->nbricks = 0; /* repetative, but clear */
}
if( slist ){ free(slist); free(sindex); }
znzclose(fp);
return NBL->nbricks;
}
/*----------------------------------------------------------------------*/
/*! nifti_free_NBL - free all pointers and clear structure
*
* note: this does not presume to free the structure pointer
*//*--------------------------------------------------------------------*/
void nifti_free_NBL( nifti_brick_list * NBL )
{
int c;
if( NBL->bricks ){
for( c = 0; c < NBL->nbricks; c++ )
if( NBL->bricks[c] ) free(NBL->bricks[c]);
free(NBL->bricks);
NBL->bricks = NULL;
}
NBL->bsize = NBL->nbricks = 0;
}
/*----------------------------------------------------------------------
* nifti_load_NBL_bricks - read the file data into the NBL struct
*
* return 0 on success, -1 on failure
*----------------------------------------------------------------------*/
static int nifti_load_NBL_bricks( nifti_image * nim , int * slist, int * sindex,
nifti_brick_list * NBL, znzFile fp )
{
size_t oposn, fposn; /* orig and current file positions */
size_t rv;
long test;
int c;
int prev, isrc, idest; /* previous and current sub-brick, and new index */
test = znztell(fp); /* store current file position */
if( test < 0 ){
fprintf(stderr,"** load bricks: ztell failed??\n");
return -1;
}
fposn = oposn = test;
/* first, handle the default case, no passed blist */
if( !slist ){
for( c = 0; c < NBL->nbricks; c++ ) {
rv = nifti_read_buffer(fp, NBL->bricks[c], NBL->bsize, nim);
if( rv != NBL->bsize ){
fprintf(stderr,"** load bricks: cannot read brick %d from '%s'\n",
c, nim->iname ? nim->iname : nim->fname);
return -1;
}
}
if( g_opts.debug > 1 )
fprintf(stderr,"+d read %d default %u-byte bricks from file %s\n",
NBL->nbricks, (unsigned int)NBL->bsize,
nim->iname ? nim->iname:nim->fname );
return 0;
}
if( !sindex ){
fprintf(stderr,"** load_NBL_bricks: missing index list\n");
return -1;
}
prev = -1; /* use prev for previous sub-brick */
for( c = 0; c < NBL->nbricks; c++ ){
isrc = slist[c]; /* this is original brick index (c is new one) */
idest = sindex[c]; /* this is the destination index for this data */
/* if this sub-brick is not the previous, we must read from disk */
if( isrc != prev ){
/* if we are not looking at the correct sub-brick, scan forward */
if( fposn != (oposn + isrc*NBL->bsize) ){
fposn = oposn + isrc*NBL->bsize;
if( znzseek(fp, (long)fposn, SEEK_SET) < 0 ){
fprintf(stderr,"** failed to locate brick %d in file '%s'\n",
isrc, nim->iname ? nim->iname : nim->fname);
return -1;
}
}
/* only 10,000 lines later and we're actually reading something! */
rv = nifti_read_buffer(fp, NBL->bricks[idest], NBL->bsize, nim);
if( rv != NBL->bsize ){
fprintf(stderr,"** failed to read brick %d from file '%s'\n",
isrc, nim->iname ? nim->iname : nim->fname);
if( g_opts.debug > 1 )
fprintf(stderr," (read %u of %u bytes)\n",
(unsigned int)rv, (unsigned int)NBL->bsize);
return -1;
}
fposn += NBL->bsize;
} else {
/* we have already read this sub-brick, just copy the previous one */
/* note that this works because they are sorted */
memcpy(NBL->bricks[idest], NBL->bricks[sindex[c-1]], NBL->bsize);
}
prev = isrc; /* in any case, note the now previous sub-brick */
}
return 0;
}
/*----------------------------------------------------------------------
* nifti_alloc_NBL_mem - allocate memory for bricks
*
* return 0 on success, -1 on failure
*----------------------------------------------------------------------*/
static int nifti_alloc_NBL_mem(nifti_image * nim, int nbricks,
nifti_brick_list * nbl)
{
int c;
/* if nbricks is not specified, use the default */
if( nbricks > 0 ) nbl->nbricks = nbricks;
else { /* I missed this one with the 1.17 change 02 Mar 2006 [rickr] */
nbl->nbricks = 1;
for( c = 4; c <= nim->ndim; c++ )
nbl->nbricks *= nim->dim[c];
}
nbl->bsize = (size_t)nim->nx * nim->ny * nim->nz * nim->nbyper;/* bytes */
nbl->bricks = (void **)malloc(nbl->nbricks * sizeof(void *));
if( ! nbl->bricks ){
fprintf(stderr,"** NANM: failed to alloc %d void ptrs\n",nbricks);
return -1;
}
for( c = 0; c < nbl->nbricks; c++ ){
nbl->bricks[c] = (void *)malloc(nbl->bsize);
if( ! nbl->bricks[c] ){
fprintf(stderr,"** NANM: failed to alloc %u bytes for brick %d\n",
(unsigned int)nbl->bsize, c);
/* so free and clear everything before returning */
while( c > 0 ){
c--;
free(nbl->bricks[c]);
}
free(nbl->bricks);
nbl->bricks = NULL;
nbl->bsize = nbl->nbricks = 0;
return -1;
}
}
if( g_opts.debug > 2 )
fprintf(stderr,"+d NANM: alloc'd %d bricks of %u bytes for NBL\n",
nbl->nbricks, (unsigned int)nbl->bsize);
return 0;
}
/*----------------------------------------------------------------------
* nifti_copynsort - copy int list, and sort with indices
*
* 1. duplicate the incoming list
* 2. create an sindex list, and init with 0..nbricks-1
* 3. do a slow insertion sort on the small slist, along with sindex list
* 4. check results, just to be positive
*
* So slist is sorted, and sindex hold original positions.
*
* return 0 on success, -1 on failure
*----------------------------------------------------------------------*/
static int nifti_copynsort(int nbricks, const int * blist, int ** slist,
int ** sindex)
{
int * stmp, * itmp; /* for ease of typing/reading */
int c1, c2, spos, tmp;
*slist = (int *)malloc(nbricks * sizeof(int));
*sindex = (int *)malloc(nbricks * sizeof(int));
if( !*slist || !*sindex ){
fprintf(stderr,"** NCS: failed to alloc %d ints for sorting\n",nbricks);
if(*slist) free(*slist); /* maybe one succeeded */
if(*sindex) free(*sindex);
return -1;
}
/* init the lists */
memcpy(*slist, blist, nbricks*sizeof(int));
for( c1 = 0; c1 < nbricks; c1++ ) (*sindex)[c1] = c1;
/* now actually sort slist */
stmp = *slist;
itmp = *sindex;
for( c1 = 0; c1 < nbricks-1; c1++ ) {
/* find smallest value, init to current */
spos = c1;
for( c2 = c1+1; c2 < nbricks; c2++ )
if( stmp[c2] < stmp[spos] ) spos = c2;
if( spos != c1 ) /* swap: fine, don't maintain sub-order, see if I care */
{
tmp = stmp[c1]; /* first swap the sorting values */
stmp[c1] = stmp[spos];
stmp[spos] = tmp;
tmp = itmp[c1]; /* then swap the index values */
itmp[c1] = itmp[spos];
itmp[spos] = tmp;
}
}
if( g_opts.debug > 2 ){
fprintf(stderr, "+d sorted indexing list:\n");
fprintf(stderr, " orig : ");
for( c1 = 0; c1 < nbricks; c1++ ) fprintf(stderr," %d",blist[c1]);
fprintf(stderr,"\n new : ");
for( c1 = 0; c1 < nbricks; c1++ ) fprintf(stderr," %d",stmp[c1]);
fprintf(stderr,"\n indices: ");
for( c1 = 0; c1 < nbricks; c1++ ) fprintf(stderr," %d",itmp[c1]);
fputc('\n', stderr);
}
/* check the sort (why not? I've got time...) */
for( c1 = 0; c1 < nbricks-1; c1++ ){
if( (stmp[c1] > stmp[c1+1]) || (blist[itmp[c1]] != stmp[c1]) ){
fprintf(stderr,"** sorting screw-up, way to go, rick!\n");
free(stmp); free(itmp); *slist = NULL; *sindex = NULL;
return -1;
}
}
if( g_opts.debug > 2 ) fprintf(stderr,"-d sorting is okay\n");
return 0;
}
/*----------------------------------------------------------------------*/
/*! valid_nifti_brick_list - check sub-brick list for image
*
* This function verifies that nbricks and blist are appropriate
* for use with this nim, based on the dimensions.
*
* \param nim nifti_image to check against
* \param nbricks number of brick indices in blist
* \param blist list of brick indices to check in nim
* \param disp_error if this flag is set, report errors to user
*
* \return 1 if valid, 0 if not
*//*--------------------------------------------------------------------*/
int valid_nifti_brick_list(nifti_image * nim , int nbricks,
const int * blist, int disp_error)
{
int c, nsubs;
if( !nim ){
if( disp_error || g_opts.debug > 0 )
fprintf(stderr,"** valid_nifti_brick_list: missing nifti image\n");
return 0;
}
if( nbricks <= 0 || !blist ){
if( disp_error || g_opts.debug > 1 )
fprintf(stderr,"** valid_nifti_brick_list: no brick list to check\n");
return 0;
}
if( nim->dim[0] < 3 ){
if( disp_error || g_opts.debug > 1 )
fprintf(stderr,"** cannot read explict brick list from %d-D dataset\n",
nim->dim[0]);
return 0;
}
/* nsubs sub-brick is nt*nu*nv*nw */
for( c = 4, nsubs = 1; c <= nim->dim[0]; c++ )
nsubs *= nim->dim[c];
if( nsubs <= 0 ){
fprintf(stderr,"** VNBL warning: bad dim list (%d,%d,%d,%d)\n",
nim->dim[4], nim->dim[5], nim->dim[6], nim->dim[7]);
return 0;
}
for( c = 0; c < nbricks; c++ )
if( (blist[c] < 0) || (blist[c] >= nsubs) ){
if( disp_error || g_opts.debug > 1 )
fprintf(stderr,
"** volume index %d (#%d) is out of range [0,%d]\n",
blist[c], c, nsubs-1);
return 0;
}
return 1; /* all is well */
}
/*----------------------------------------------------------------------*/
/* verify that NBL struct is a valid data source for the image
*
* return 1 if so, 0 otherwise
*//*--------------------------------------------------------------------*/
static int nifti_NBL_matches_nim(const nifti_image *nim,
const nifti_brick_list *NBL)
{
size_t volbytes = 0; /* bytes per volume */
int ind, errs = 0, nvols = 0;
if( !nim || !NBL ) {
if( g_opts.debug > 0 )
fprintf(stderr,"** nifti_NBL_matches_nim: NULL pointer(s)\n");
return 0;
}
/* for nim, compute volbytes and nvols */
if( nim->ndim > 0 ) {
/* first 3 indices are over a single volume */
volbytes = (size_t)nim->nbyper;
for( ind = 1; ind <= nim->ndim && ind < 4; ind++ )
volbytes *= (size_t)nim->dim[ind];
for( ind = 4, nvols = 1; ind <= nim->ndim; ind++ )
nvols *= nim->dim[ind];
}
if( volbytes != NBL->bsize ) {
if( g_opts.debug > 1 )
fprintf(stderr,"** NBL/nim mismatch, volbytes = %u, %u\n",
(unsigned)NBL->bsize, (unsigned)volbytes);
errs++;
}
if( nvols != NBL->nbricks ) {
if( g_opts.debug > 1 )
fprintf(stderr,"** NBL/nim mismatch, nvols = %d, %d\n",
NBL->nbricks, nvols);
errs++;
}
if( errs ) return 0;
else if ( g_opts.debug > 2 )
fprintf(stderr,"-- nim/NBL agree: nvols = %d, nbytes = %u\n",
nvols, (unsigned)volbytes);
return 1;
}
/* end of new nifti_image_read_bricks() functionality */
/*----------------------------------------------------------------------*/
/*! display the orientation from the quaternian fields
*
* \param mesg if non-NULL, display this message first
* \param mat the matrix to convert to "nearest" orientation
*
* \return -1 if results cannot be determined, 0 if okay
*//*--------------------------------------------------------------------*/
int nifti_disp_matrix_orient( const char * mesg, mat44 mat )
{
int i, j, k;
if ( mesg ) fputs( mesg, stderr ); /* use stdout? */
nifti_mat44_to_orientation( mat, &i,&j,&k );
if ( i <= 0 || j <= 0 || k <= 0 ) return -1;
/* so we have good codes */
fprintf(stderr, " i orientation = '%s'\n"
" j orientation = '%s'\n"
" k orientation = '%s'\n",
nifti_orientation_string(i),
nifti_orientation_string(j),
nifti_orientation_string(k) );
return 0;
}
/*----------------------------------------------------------------------*/
/*! duplicate the given string (alloc length+1)
*
* \return allocated pointer (or NULL on failure)
*//*--------------------------------------------------------------------*/
char *nifti_strdup(const char *str)
{
char *dup;
if( !str ) return NULL; /* allow calls passing NULL */
dup = (char *)malloc(strlen(str) + 1);
/* check for failure */
if( dup ) strcpy(dup, str);
else fprintf(stderr,"** nifti_strdup: failed to alloc %u bytes\n",
(unsigned int)strlen(str)+1);
return dup;
}
/*---------------------------------------------------------------------------*/
/*! Return a pointer to a string holding the name of a NIFTI datatype.
\param dt NIfTI-1 datatype
\return pointer to static string holding the datatype name
\warning Do not free() or modify this string!
It points to static storage.
\sa NIFTI1_DATATYPES group in nifti1.h
*//*-------------------------------------------------------------------------*/
char *nifti_datatype_string( int dt )
{
switch( dt ){
case DT_UNKNOWN: return "UNKNOWN" ;
case DT_BINARY: return "BINARY" ;
case DT_INT8: return "INT8" ;
case DT_UINT8: return "UINT8" ;
case DT_INT16: return "INT16" ;
case DT_UINT16: return "UINT16" ;
case DT_INT32: return "INT32" ;
case DT_UINT32: return "UINT32" ;
case DT_INT64: return "INT64" ;
case DT_UINT64: return "UINT64" ;
case DT_FLOAT32: return "FLOAT32" ;
case DT_FLOAT64: return "FLOAT64" ;
case DT_FLOAT128: return "FLOAT128" ;
case DT_COMPLEX64: return "COMPLEX64" ;
case DT_COMPLEX128: return "COMPLEX128" ;
case DT_COMPLEX256: return "COMPLEX256" ;
case DT_RGB24: return "RGB24" ;
case DT_RGBA32: return "RGBA32" ;
}
return "**ILLEGAL**" ;
}
/*----------------------------------------------------------------------*/
/*! Determine if the datatype code dt is an integer type (1=YES, 0=NO).
\return whether the given NIfTI-1 datatype code is valid
\sa NIFTI1_DATATYPES group in nifti1.h
*//*--------------------------------------------------------------------*/
int nifti_is_inttype( int dt )
{
switch( dt ){
case DT_UNKNOWN: return 0 ;
case DT_BINARY: return 0 ;
case DT_INT8: return 1 ;
case DT_UINT8: return 1 ;
case DT_INT16: return 1 ;
case DT_UINT16: return 1 ;
case DT_INT32: return 1 ;
case DT_UINT32: return 1 ;
case DT_INT64: return 1 ;
case DT_UINT64: return 1 ;
case DT_FLOAT32: return 0 ;
case DT_FLOAT64: return 0 ;
case DT_FLOAT128: return 0 ;
case DT_COMPLEX64: return 0 ;
case DT_COMPLEX128: return 0 ;
case DT_COMPLEX256: return 0 ;
case DT_RGB24: return 1 ;
case DT_RGBA32: return 1 ;
}
return 0 ;
}
/*---------------------------------------------------------------------------*/
/*! Return a pointer to a string holding the name of a NIFTI units type.
\param uu NIfTI-1 unit code
\return pointer to static string for the given unit type
\warning Do not free() or modify this string!
It points to static storage.
\sa NIFTI1_UNITS group in nifti1.h
*//*-------------------------------------------------------------------------*/
char *nifti_units_string( int uu )
{
switch( uu ){
case NIFTI_UNITS_METER: return "m" ;
case NIFTI_UNITS_MM: return "mm" ;
case NIFTI_UNITS_MICRON: return "um" ;
case NIFTI_UNITS_SEC: return "s" ;
case NIFTI_UNITS_MSEC: return "ms" ;
case NIFTI_UNITS_USEC: return "us" ;
case NIFTI_UNITS_HZ: return "Hz" ;
case NIFTI_UNITS_PPM: return "ppm" ;
case NIFTI_UNITS_RADS: return "rad/s" ;
}
return "Unknown" ;
}
/*---------------------------------------------------------------------------*/
/*! Return a pointer to a string holding the name of a NIFTI transform type.
\param xx NIfTI-1 xform code
\return pointer to static string describing xform code
\warning Do not free() or modify this string!
It points to static storage.
\sa NIFTI1_XFORM_CODES group in nifti1.h
*//*-------------------------------------------------------------------------*/
char *nifti_xform_string( int xx )
{
switch( xx ){
case NIFTI_XFORM_SCANNER_ANAT: return "Scanner Anat" ;
case NIFTI_XFORM_ALIGNED_ANAT: return "Aligned Anat" ;
case NIFTI_XFORM_TALAIRACH: return "Talairach" ;
case NIFTI_XFORM_MNI_152: return "MNI_152" ;
}
return "Unknown" ;
}
/*---------------------------------------------------------------------------*/
/*! Return a pointer to a string holding the name of a NIFTI intent type.
\param ii NIfTI-1 intent code
\return pointer to static string describing code
\warning Do not free() or modify this string!
It points to static storage.
\sa NIFTI1_INTENT_CODES group in nifti1.h
*//*-------------------------------------------------------------------------*/
char *nifti_intent_string( int ii )
{
switch( ii ){
case NIFTI_INTENT_CORREL: return "Correlation statistic" ;
case NIFTI_INTENT_TTEST: return "T-statistic" ;
case NIFTI_INTENT_FTEST: return "F-statistic" ;
case NIFTI_INTENT_ZSCORE: return "Z-score" ;
case NIFTI_INTENT_CHISQ: return "Chi-squared distribution" ;
case NIFTI_INTENT_BETA: return "Beta distribution" ;
case NIFTI_INTENT_BINOM: return "Binomial distribution" ;
case NIFTI_INTENT_GAMMA: return "Gamma distribution" ;
case NIFTI_INTENT_POISSON: return "Poisson distribution" ;
case NIFTI_INTENT_NORMAL: return "Normal distribution" ;
case NIFTI_INTENT_FTEST_NONC: return "F-statistic noncentral" ;
case NIFTI_INTENT_CHISQ_NONC: return "Chi-squared noncentral" ;
case NIFTI_INTENT_LOGISTIC: return "Logistic distribution" ;
case NIFTI_INTENT_LAPLACE: return "Laplace distribution" ;
case NIFTI_INTENT_UNIFORM: return "Uniform distribition" ;
case NIFTI_INTENT_TTEST_NONC: return "T-statistic noncentral" ;
case NIFTI_INTENT_WEIBULL: return "Weibull distribution" ;
case NIFTI_INTENT_CHI: return "Chi distribution" ;
case NIFTI_INTENT_INVGAUSS: return "Inverse Gaussian distribution" ;
case NIFTI_INTENT_EXTVAL: return "Extreme Value distribution" ;
case NIFTI_INTENT_PVAL: return "P-value" ;
case NIFTI_INTENT_LOGPVAL: return "Log P-value" ;
case NIFTI_INTENT_LOG10PVAL: return "Log10 P-value" ;
case NIFTI_INTENT_ESTIMATE: return "Estimate" ;
case NIFTI_INTENT_LABEL: return "Label index" ;
case NIFTI_INTENT_NEURONAME: return "NeuroNames index" ;
case NIFTI_INTENT_GENMATRIX: return "General matrix" ;
case NIFTI_INTENT_SYMMATRIX: return "Symmetric matrix" ;
case NIFTI_INTENT_DISPVECT: return "Displacement vector" ;
case NIFTI_INTENT_VECTOR: return "Vector" ;
case NIFTI_INTENT_POINTSET: return "Pointset" ;
case NIFTI_INTENT_TRIANGLE: return "Triangle" ;
case NIFTI_INTENT_QUATERNION: return "Quaternion" ;
case NIFTI_INTENT_DIMLESS: return "Dimensionless number" ;
}
return "Unknown" ;
}
/*---------------------------------------------------------------------------*/
/*! Return a pointer to a string holding the name of a NIFTI slice_code.
\param ss NIfTI-1 slice order code
\return pointer to static string describing code
\warning Do not free() or modify this string!
It points to static storage.
\sa NIFTI1_SLICE_ORDER group in nifti1.h
*//*-------------------------------------------------------------------------*/
char *nifti_slice_string( int ss )
{
switch( ss ){
case NIFTI_SLICE_SEQ_INC: return "sequential_increasing" ;
case NIFTI_SLICE_SEQ_DEC: return "sequential_decreasing" ;
case NIFTI_SLICE_ALT_INC: return "alternating_increasing" ;
case NIFTI_SLICE_ALT_DEC: return "alternating_decreasing" ;
case NIFTI_SLICE_ALT_INC2: return "alternating_increasing_2" ;
case NIFTI_SLICE_ALT_DEC2: return "alternating_decreasing_2" ;
}
return "Unknown" ;
}
/*---------------------------------------------------------------------------*/
/*! Return a pointer to a string holding the name of a NIFTI orientation.
\param ii orientation code
\return pointer to static string holding the orientation information
\warning Do not free() or modify the return string!
It points to static storage.
\sa NIFTI_L2R in nifti1_io.h
*//*-------------------------------------------------------------------------*/
char *nifti_orientation_string( int ii )
{
switch( ii ){
case NIFTI_L2R: return "Left-to-Right" ;
case NIFTI_R2L: return "Right-to-Left" ;
case NIFTI_P2A: return "Posterior-to-Anterior" ;
case NIFTI_A2P: return "Anterior-to-Posterior" ;
case NIFTI_I2S: return "Inferior-to-Superior" ;
case NIFTI_S2I: return "Superior-to-Inferior" ;
}
return "Unknown" ;
}
/*--------------------------------------------------------------------------*/
/*! Given a datatype code, set number of bytes per voxel and the swapsize.
\param datatype nifti1 datatype code
\param nbyper pointer to return value: number of bytes per voxel
\param swapsize pointer to return value: size of swap blocks
\return appropriate values at nbyper and swapsize
The swapsize is set to 0 if this datatype doesn't ever need swapping.
\sa NIFTI1_DATATYPES in nifti1.h
*//*------------------------------------------------------------------------*/
void nifti_datatype_sizes( int datatype , int *nbyper, int *swapsize )
{
int nb=0, ss=0 ;
switch( datatype ){
case DT_INT8:
case DT_UINT8: nb = 1 ; ss = 0 ; break ;
case DT_INT16:
case DT_UINT16: nb = 2 ; ss = 2 ; break ;
case DT_RGB24: nb = 3 ; ss = 0 ; break ;
case DT_RGBA32: nb = 4 ; ss = 0 ; break ;
case DT_INT32:
case DT_UINT32:
case DT_FLOAT32: nb = 4 ; ss = 4 ; break ;
case DT_COMPLEX64: nb = 8 ; ss = 4 ; break ;
case DT_FLOAT64:
case DT_INT64:
case DT_UINT64: nb = 8 ; ss = 8 ; break ;
case DT_FLOAT128: nb = 16 ; ss = 16 ; break ;
case DT_COMPLEX128: nb = 16 ; ss = 8 ; break ;
case DT_COMPLEX256: nb = 32 ; ss = 16 ; break ;
}
ASSIF(nbyper,nb) ; ASSIF(swapsize,ss) ; return ;
}
/*---------------------------------------------------------------------------*/
/*! Given the quaternion parameters (etc.), compute a transformation matrix.
See comments in nifti1.h for details.
- qb,qc,qd = quaternion parameters
- qx,qy,qz = offset parameters
- dx,dy,dz = grid stepsizes (non-negative inputs are set to 1.0)
- qfac = sign of dz step (< 0 is negative; >= 0 is positive)
If qx=qy=qz=0, dx=dy=dz=1, then the output is a rotation matrix. For qfac >= 0, the rotation is proper. For qfac < 0, the rotation is improper.\see "QUATERNION REPRESENTATION OF ROTATION MATRIX" in nifti1.h \see nifti_mat44_to_quatern, nifti_make_orthog_mat44, nifti_mat44_to_orientation *//*-------------------------------------------------------------------------*/ mat44 nifti_quatern_to_mat44( float qb, float qc, float qd, float qx, float qy, float qz, float dx, float dy, float dz, float qfac ) { mat44 R ; double a,b=qb,c=qc,d=qd , xd,yd,zd ; /* last row is always [ 0 0 0 1 ] */ R.m[3][0]=R.m[3][1]=R.m[3][2] = 0.0 ; R.m[3][3]= 1.0 ; /* compute a parameter from b,c,d */ a = 1.0l - (b*b + c*c + d*d) ; if( a < 1.e-7l ){ /* special case */ a = 1.0l / sqrt(b*b+c*c+d*d) ; b *= a ; c *= a ; d *= a ; /* normalize (b,c,d) vector */ a = 0.0l ; /* a = 0 ==> 180 degree rotation */ } else{ a = sqrt(a) ; /* angle = 2*arccos(a) */ } /* load rotation matrix, including scaling factors for voxel sizes */ xd = (dx > 0.0) ? dx : 1.0l ; /* make sure are positive */ yd = (dy > 0.0) ? dy : 1.0l ; zd = (dz > 0.0) ? dz : 1.0l ; if( qfac < 0.0 ) zd = -zd ; /* left handedness? */ R.m[0][0] = (a*a+b*b-c*c-d*d) * xd ; R.m[0][1] = 2.0l * (b*c-a*d ) * yd ; R.m[0][2] = 2.0l * (b*d+a*c ) * zd ; R.m[1][0] = 2.0l * (b*c+a*d ) * xd ; R.m[1][1] = (a*a+c*c-b*b-d*d) * yd ; R.m[1][2] = 2.0l * (c*d-a*b ) * zd ; R.m[2][0] = 2.0l * (b*d-a*c ) * xd ; R.m[2][1] = 2.0l * (c*d+a*b ) * yd ; R.m[2][2] = (a*a+d*d-c*c-b*b) * zd ; /* load offsets */ R.m[0][3] = qx ; R.m[1][3] = qy ; R.m[2][3] = qz ; return R ; } /*---------------------------------------------------------------------------*/ /*! Given the 3x4 upper corner of the matrix R, compute the quaternion parameters that fit it. - Any NULL pointer on input won't get assigned (e.g., if you don't want dx,dy,dz, just pass NULL in for those pointers). - If the 3 input matrix columns are NOT orthogonal, they will be orthogonalized prior to calculating the parameters, using the polar decomposition to find the orthogonal matrix closest to the column-normalized input matrix. - However, if the 3 input matrix columns are NOT orthogonal, then the matrix produced by nifti_quatern_to_mat44 WILL have orthogonal columns, so it won't be the same as the matrix input here. This "feature" is because the NIFTI 'qform' transform is deliberately not fully general -- it is intended to model a volume with perpendicular axes. - If the 3 input matrix columns are not even linearly independent, you'll just have to take your luck, won't you? \see "QUATERNION REPRESENTATION OF ROTATION MATRIX" in nifti1.h \see nifti_quatern_to_mat44, nifti_make_orthog_mat44, nifti_mat44_to_orientation *//*-------------------------------------------------------------------------*/ void nifti_mat44_to_quatern( mat44 R , float *qb, float *qc, float *qd, float *qx, float *qy, float *qz, float *dx, float *dy, float *dz, float *qfac ) { double r11,r12,r13 , r21,r22,r23 , r31,r32,r33 ; double xd,yd,zd , a,b,c,d ; mat33 P,Q ; /* offset outputs are read write out of input matrix */ ASSIF(qx,R.m[0][3]) ; ASSIF(qy,R.m[1][3]) ; ASSIF(qz,R.m[2][3]) ; /* load 3x3 matrix into local variables */ r11 = R.m[0][0] ; r12 = R.m[0][1] ; r13 = R.m[0][2] ; r21 = R.m[1][0] ; r22 = R.m[1][1] ; r23 = R.m[1][2] ; r31 = R.m[2][0] ; r32 = R.m[2][1] ; r33 = R.m[2][2] ; /* compute lengths of each column; these determine grid spacings */ xd = sqrt( r11*r11 + r21*r21 + r31*r31 ) ; yd = sqrt( r12*r12 + r22*r22 + r32*r32 ) ; zd = sqrt( r13*r13 + r23*r23 + r33*r33 ) ; /* if a column length is zero, patch the trouble */ if( xd == 0.0l ){ r11 = 1.0l ; r21 = r31 = 0.0l ; xd = 1.0l ; } if( yd == 0.0l ){ r22 = 1.0l ; r12 = r32 = 0.0l ; yd = 1.0l ; } if( zd == 0.0l ){ r33 = 1.0l ; r13 = r23 = 0.0l ; zd = 1.0l ; } /* assign the output lengths */ ASSIF(dx,xd) ; ASSIF(dy,yd) ; ASSIF(dz,zd) ; /* normalize the columns */ r11 /= xd ; r21 /= xd ; r31 /= xd ; r12 /= yd ; r22 /= yd ; r32 /= yd ; r13 /= zd ; r23 /= zd ; r33 /= zd ; /* At this point, the matrix has normal columns, but we have to allow for the fact that the hideous user may not have given us a matrix with orthogonal columns. So, now find the orthogonal matrix closest to the current matrix. One reason for using the polar decomposition to get this orthogonal matrix, rather than just directly orthogonalizing the columns, is so that inputting the inverse matrix to R will result in the inverse orthogonal matrix at this point. If we just orthogonalized the columns, this wouldn't necessarily hold. */ Q.m[0][0] = r11 ; Q.m[0][1] = r12 ; Q.m[0][2] = r13 ; /* load Q */ Q.m[1][0] = r21 ; Q.m[1][1] = r22 ; Q.m[1][2] = r23 ; Q.m[2][0] = r31 ; Q.m[2][1] = r32 ; Q.m[2][2] = r33 ; P = nifti_mat33_polar(Q) ; /* P is orthog matrix closest to Q */ r11 = P.m[0][0] ; r12 = P.m[0][1] ; r13 = P.m[0][2] ; /* unload */ r21 = P.m[1][0] ; r22 = P.m[1][1] ; r23 = P.m[1][2] ; r31 = P.m[2][0] ; r32 = P.m[2][1] ; r33 = P.m[2][2] ; /* [ r11 r12 r13 ] */ /* at this point, the matrix [ r21 r22 r23 ] is orthogonal */ /* [ r31 r32 r33 ] */ /* compute the determinant to determine if it is proper */ zd = r11*r22*r33-r11*r32*r23-r21*r12*r33 +r21*r32*r13+r31*r12*r23-r31*r22*r13 ; /* should be -1 or 1 */ if( zd > 0 ){ /* proper */ ASSIF(qfac,1.0) ; } else { /* improper ==> flip 3rd column */ ASSIF(qfac,-1.0) ; r13 = -r13 ; r23 = -r23 ; r33 = -r33 ; } /* now, compute quaternion parameters */ a = r11 + r22 + r33 + 1.0l ; if( a > 0.5l ){ /* simplest case */ a = 0.5l * sqrt(a) ; b = 0.25l * (r32-r23) / a ; c = 0.25l * (r13-r31) / a ; d = 0.25l * (r21-r12) / a ; } else { /* trickier case */ xd = 1.0 + r11 - (r22+r33) ; /* 4*b*b */ yd = 1.0 + r22 - (r11+r33) ; /* 4*c*c */ zd = 1.0 + r33 - (r11+r22) ; /* 4*d*d */ if( xd > 1.0 ){ b = 0.5l * sqrt(xd) ; c = 0.25l* (r12+r21) / b ; d = 0.25l* (r13+r31) / b ; a = 0.25l* (r32-r23) / b ; } else if( yd > 1.0 ){ c = 0.5l * sqrt(yd) ; b = 0.25l* (r12+r21) / c ; d = 0.25l* (r23+r32) / c ; a = 0.25l* (r13-r31) / c ; } else { d = 0.5l * sqrt(zd) ; b = 0.25l* (r13+r31) / d ; c = 0.25l* (r23+r32) / d ; a = 0.25l* (r21-r12) / d ; } if( a < 0.0l ){ b=-b ; c=-c ; d=-d; a=-a; } } ASSIF(qb,b) ; ASSIF(qc,c) ; ASSIF(qd,d) ; return ; } /*---------------------------------------------------------------------------*/ /*! Compute the inverse of a bordered 4x4 matrix.
- Some numerical code fragments were generated by Maple 8.
- If a singular matrix is input, the output matrix will be all zero.
- You can check for this by examining the [3][3] element, which will
be 1.0 for the normal case and 0.0 for the bad case.
The input matrix should have the form:
[ r11 r12 r13 v1 ]
[ r21 r22 r23 v2 ]
[ r31 r32 r33 v3 ]
[ 0 0 0 1 ]
*//*-------------------------------------------------------------------------*/
mat44 nifti_mat44_inverse( mat44 R )
{
double r11,r12,r13,r21,r22,r23,r31,r32,r33,v1,v2,v3 , deti ;
mat44 Q ;
/* INPUT MATRIX IS: */
r11 = R.m[0][0]; r12 = R.m[0][1]; r13 = R.m[0][2]; /* [ r11 r12 r13 v1 ] */
r21 = R.m[1][0]; r22 = R.m[1][1]; r23 = R.m[1][2]; /* [ r21 r22 r23 v2 ] */
r31 = R.m[2][0]; r32 = R.m[2][1]; r33 = R.m[2][2]; /* [ r31 r32 r33 v3 ] */
v1 = R.m[0][3]; v2 = R.m[1][3]; v3 = R.m[2][3]; /* [ 0 0 0 1 ] */
deti = r11*r22*r33-r11*r32*r23-r21*r12*r33
+r21*r32*r13+r31*r12*r23-r31*r22*r13 ;
if( deti != 0.0l ) deti = 1.0l / deti ;
Q.m[0][0] = deti*( r22*r33-r32*r23) ;
Q.m[0][1] = deti*(-r12*r33+r32*r13) ;
Q.m[0][2] = deti*( r12*r23-r22*r13) ;
Q.m[0][3] = deti*(-r12*r23*v3+r12*v2*r33+r22*r13*v3
-r22*v1*r33-r32*r13*v2+r32*v1*r23) ;
Q.m[1][0] = deti*(-r21*r33+r31*r23) ;
Q.m[1][1] = deti*( r11*r33-r31*r13) ;
Q.m[1][2] = deti*(-r11*r23+r21*r13) ;
Q.m[1][3] = deti*( r11*r23*v3-r11*v2*r33-r21*r13*v3
+r21*v1*r33+r31*r13*v2-r31*v1*r23) ;
Q.m[2][0] = deti*( r21*r32-r31*r22) ;
Q.m[2][1] = deti*(-r11*r32+r31*r12) ;
Q.m[2][2] = deti*( r11*r22-r21*r12) ;
Q.m[2][3] = deti*(-r11*r22*v3+r11*r32*v2+r21*r12*v3
-r21*r32*v1-r31*r12*v2+r31*r22*v1) ;
Q.m[3][0] = Q.m[3][1] = Q.m[3][2] = 0.0l ;
Q.m[3][3] = (deti == 0.0l) ? 0.0l : 1.0l ; /* failure flag if deti == 0 */
return Q ;
}
/*---------------------------------------------------------------------------*/
/*! Input 9 floats and make an orthgonal mat44 out of them.
Each row is normalized, then nifti_mat33_polar() is used to orthogonalize
them. If row #3 (r31,r32,r33) is input as zero, then it will be taken to
be the cross product of rows #1 and #2.
This function can be used to create a rotation matrix for transforming
an oblique volume to anatomical coordinates. For this application:
- row #1 (r11,r12,r13) is the direction vector along the image i-axis
- row #2 (r21,r22,r23) is the direction vector along the image j-axis
- row #3 (r31,r32,r33) is the direction vector along the slice direction
(if available; otherwise enter it as 0's)
The first 2 rows can be taken from the DICOM attribute (0020,0037)
"Image Orientation (Patient)".
After forming the rotation matrix, the complete affine transformation from
(i,j,k) grid indexes to (x,y,z) spatial coordinates can be computed by
multiplying each column by the appropriate grid spacing:
- column #1 (R.m[0][0],R.m[1][0],R.m[2][0]) by delta-x
- column #2 (R.m[0][1],R.m[1][1],R.m[2][1]) by delta-y
- column #3 (R.m[0][2],R.m[1][2],R.m[2][2]) by delta-z
and by then placing the center (x,y,z) coordinates of voxel (0,0,0) into
the column #4 (R.m[0][3],R.m[1][3],R.m[2][3]).
\sa nifti_quatern_to_mat44, nifti_mat44_to_quatern,
nifti_mat44_to_orientation
*//*-------------------------------------------------------------------------*/
mat44 nifti_make_orthog_mat44( float r11, float r12, float r13 ,
float r21, float r22, float r23 ,
float r31, float r32, float r33 )
{
mat44 R ;
mat33 Q , P ;
double val ;
R.m[3][0] = R.m[3][1] = R.m[3][2] = 0.0l ; R.m[3][3] = 1.0l ;
Q.m[0][0] = r11 ; Q.m[0][1] = r12 ; Q.m[0][2] = r13 ; /* load Q */
Q.m[1][0] = r21 ; Q.m[1][1] = r22 ; Q.m[1][2] = r23 ;
Q.m[2][0] = r31 ; Q.m[2][1] = r32 ; Q.m[2][2] = r33 ;
/* normalize row 1 */
val = Q.m[0][0]*Q.m[0][0] + Q.m[0][1]*Q.m[0][1] + Q.m[0][2]*Q.m[0][2] ;
if( val > 0.0l ){
val = 1.0l / sqrt(val) ;
Q.m[0][0] *= val ; Q.m[0][1] *= val ; Q.m[0][2] *= val ;
} else {
Q.m[0][0] = 1.0l ; Q.m[0][1] = 0.0l ; Q.m[0][2] = 0.0l ;
}
/* normalize row 2 */
val = Q.m[1][0]*Q.m[1][0] + Q.m[1][1]*Q.m[1][1] + Q.m[1][2]*Q.m[1][2] ;
if( val > 0.0l ){
val = 1.0l / sqrt(val) ;
Q.m[1][0] *= val ; Q.m[1][1] *= val ; Q.m[1][2] *= val ;
} else {
Q.m[1][0] = 0.0l ; Q.m[1][1] = 1.0l ; Q.m[1][2] = 0.0l ;
}
/* normalize row 3 */
val = Q.m[2][0]*Q.m[2][0] + Q.m[2][1]*Q.m[2][1] + Q.m[2][2]*Q.m[2][2] ;
if( val > 0.0l ){
val = 1.0l / sqrt(val) ;
Q.m[2][0] *= val ; Q.m[2][1] *= val ; Q.m[2][2] *= val ;
} else {
Q.m[2][0] = Q.m[0][1]*Q.m[1][2] - Q.m[0][2]*Q.m[1][1] ; /* cross */
Q.m[2][1] = Q.m[0][2]*Q.m[1][0] - Q.m[0][0]*Q.m[1][2] ; /* product */
Q.m[2][2] = Q.m[0][0]*Q.m[1][1] - Q.m[0][1]*Q.m[1][0] ;
}
P = nifti_mat33_polar(Q) ; /* P is orthog matrix closest to Q */
R.m[0][0] = P.m[0][0] ; R.m[0][1] = P.m[0][1] ; R.m[0][2] = P.m[0][2] ;
R.m[1][0] = P.m[1][0] ; R.m[1][1] = P.m[1][1] ; R.m[1][2] = P.m[1][2] ;
R.m[2][0] = P.m[2][0] ; R.m[2][1] = P.m[2][1] ; R.m[2][2] = P.m[2][2] ;
R.m[0][3] = R.m[1][3] = R.m[2][3] = 0.0 ; return R ;
}
/*----------------------------------------------------------------------*/
/*! compute the inverse of a 3x3 matrix
*//*--------------------------------------------------------------------*/
mat33 nifti_mat33_inverse( mat33 R ) /* inverse of 3x3 matrix */
{
double r11,r12,r13,r21,r22,r23,r31,r32,r33 , deti ;
mat33 Q ;
/* INPUT MATRIX: */
r11 = R.m[0][0]; r12 = R.m[0][1]; r13 = R.m[0][2]; /* [ r11 r12 r13 ] */
r21 = R.m[1][0]; r22 = R.m[1][1]; r23 = R.m[1][2]; /* [ r21 r22 r23 ] */
r31 = R.m[2][0]; r32 = R.m[2][1]; r33 = R.m[2][2]; /* [ r31 r32 r33 ] */
deti = r11*r22*r33-r11*r32*r23-r21*r12*r33
+r21*r32*r13+r31*r12*r23-r31*r22*r13 ;
if( deti != 0.0l ) deti = 1.0l / deti ;
Q.m[0][0] = deti*( r22*r33-r32*r23) ;
Q.m[0][1] = deti*(-r12*r33+r32*r13) ;
Q.m[0][2] = deti*( r12*r23-r22*r13) ;
Q.m[1][0] = deti*(-r21*r33+r31*r23) ;
Q.m[1][1] = deti*( r11*r33-r31*r13) ;
Q.m[1][2] = deti*(-r11*r23+r21*r13) ;
Q.m[2][0] = deti*( r21*r32-r31*r22) ;
Q.m[2][1] = deti*(-r11*r32+r31*r12) ;
Q.m[2][2] = deti*( r11*r22-r21*r12) ;
return Q ;
}
/*----------------------------------------------------------------------*/
/*! compute the determinant of a 3x3 matrix
*//*--------------------------------------------------------------------*/
float nifti_mat33_determ( mat33 R ) /* determinant of 3x3 matrix */
{
double r11,r12,r13,r21,r22,r23,r31,r32,r33 ;
/* INPUT MATRIX: */
r11 = R.m[0][0]; r12 = R.m[0][1]; r13 = R.m[0][2]; /* [ r11 r12 r13 ] */
r21 = R.m[1][0]; r22 = R.m[1][1]; r23 = R.m[1][2]; /* [ r21 r22 r23 ] */
r31 = R.m[2][0]; r32 = R.m[2][1]; r33 = R.m[2][2]; /* [ r31 r32 r33 ] */
return r11*r22*r33-r11*r32*r23-r21*r12*r33
+r21*r32*r13+r31*r12*r23-r31*r22*r13 ;
}
/*----------------------------------------------------------------------*/
/*! compute the max row norm of a 3x3 matrix
*//*--------------------------------------------------------------------*/
float nifti_mat33_rownorm( mat33 A ) /* max row norm of 3x3 matrix */
{
float r1,r2,r3 ;
r1 = fabs(A.m[0][0])+fabs(A.m[0][1])+fabs(A.m[0][2]) ;
r2 = fabs(A.m[1][0])+fabs(A.m[1][1])+fabs(A.m[1][2]) ;
r3 = fabs(A.m[2][0])+fabs(A.m[2][1])+fabs(A.m[2][2]) ;
if( r1 < r2 ) r1 = r2 ;
if( r1 < r3 ) r1 = r3 ;
return r1 ;
}
/*----------------------------------------------------------------------*/
/*! compute the max column norm of a 3x3 matrix
*//*--------------------------------------------------------------------*/
float nifti_mat33_colnorm( mat33 A ) /* max column norm of 3x3 matrix */
{
float r1,r2,r3 ;
r1 = fabs(A.m[0][0])+fabs(A.m[1][0])+fabs(A.m[2][0]) ;
r2 = fabs(A.m[0][1])+fabs(A.m[1][1])+fabs(A.m[2][1]) ;
r3 = fabs(A.m[0][2])+fabs(A.m[1][2])+fabs(A.m[2][2]) ;
if( r1 < r2 ) r1 = r2 ;
if( r1 < r3 ) r1 = r3 ;
return r1 ;
}
/*----------------------------------------------------------------------*/
/*! multiply 2 3x3 matrices
*//*--------------------------------------------------------------------*/
mat33 nifti_mat33_mul( mat33 A , mat33 B ) /* multiply 2 3x3 matrices */
{
mat33 C ; int i,j ;
for( i=0 ; i < 3 ; i++ )
for( j=0 ; j < 3 ; j++ )
C.m[i][j] = A.m[i][0] * B.m[0][j]
+ A.m[i][1] * B.m[1][j]
+ A.m[i][2] * B.m[2][j] ;
return C ;
}
/*---------------------------------------------------------------------------*/
/*! polar decomposition of a 3x3 matrix
This finds the closest orthogonal matrix to input A
(in both the Frobenius and L2 norms).
Algorithm is that from NJ Higham, SIAM J Sci Stat Comput, 7:1160-1174.
*//*-------------------------------------------------------------------------*/
mat33 nifti_mat33_polar( mat33 A )
{
mat33 X , Y , Z ;
float alp,bet,gam,gmi , dif=1.0 ;
int k=0 ;
X = A ;
/* force matrix to be nonsingular */
gam = nifti_mat33_determ(X) ;
while( gam == 0.0 ){ /* perturb matrix */
gam = 0.00001 * ( 0.001 + nifti_mat33_rownorm(X) ) ;
X.m[0][0] += gam ; X.m[1][1] += gam ; X.m[2][2] += gam ;
gam = nifti_mat33_determ(X) ;
}
while(1){
Y = nifti_mat33_inverse(X) ;
if( dif > 0.3 ){ /* far from convergence */
alp = sqrt( nifti_mat33_rownorm(X) * nifti_mat33_colnorm(X) ) ;
bet = sqrt( nifti_mat33_rownorm(Y) * nifti_mat33_colnorm(Y) ) ;
gam = sqrt( bet / alp ) ;
gmi = 1.0 / gam ;
} else {
gam = gmi = 1.0 ; /* close to convergence */
}
Z.m[0][0] = 0.5 * ( gam*X.m[0][0] + gmi*Y.m[0][0] ) ;
Z.m[0][1] = 0.5 * ( gam*X.m[0][1] + gmi*Y.m[1][0] ) ;
Z.m[0][2] = 0.5 * ( gam*X.m[0][2] + gmi*Y.m[2][0] ) ;
Z.m[1][0] = 0.5 * ( gam*X.m[1][0] + gmi*Y.m[0][1] ) ;
Z.m[1][1] = 0.5 * ( gam*X.m[1][1] + gmi*Y.m[1][1] ) ;
Z.m[1][2] = 0.5 * ( gam*X.m[1][2] + gmi*Y.m[2][1] ) ;
Z.m[2][0] = 0.5 * ( gam*X.m[2][0] + gmi*Y.m[0][2] ) ;
Z.m[2][1] = 0.5 * ( gam*X.m[2][1] + gmi*Y.m[1][2] ) ;
Z.m[2][2] = 0.5 * ( gam*X.m[2][2] + gmi*Y.m[2][2] ) ;
dif = fabs(Z.m[0][0]-X.m[0][0])+fabs(Z.m[0][1]-X.m[0][1])
+fabs(Z.m[0][2]-X.m[0][2])+fabs(Z.m[1][0]-X.m[1][0])
+fabs(Z.m[1][1]-X.m[1][1])+fabs(Z.m[1][2]-X.m[1][2])
+fabs(Z.m[2][0]-X.m[2][0])+fabs(Z.m[2][1]-X.m[2][1])
+fabs(Z.m[2][2]-X.m[2][2]) ;
k = k+1 ;
if( k > 100 || dif < 3.e-6 ) break ; /* convergence or exhaustion */
X = Z ;
}
return Z ;
}
/*---------------------------------------------------------------------------*/
/*! compute the (closest) orientation from a 4x4 ijk->xyz tranformation matrix
Input: 4x4 matrix that transforms (i,j,k) indexes to (x,y,z) coordinates,
where +x=Right, +y=Anterior, +z=Superior.
(Only the upper-left 3x3 corner of R is used herein.)
Output: 3 orientation codes that correspond to the closest "standard"
anatomical orientation of the (i,j,k) axes.
Method: Find which permutation of (x,y,z) has the smallest angle to the
(i,j,k) axes directions, which are the columns of the R matrix.
Errors: The codes returned will be zero.
For example, an axial volume might get return values of
*icod = NIFTI_R2L (i axis is mostly Right to Left)
*jcod = NIFTI_P2A (j axis is mostly Posterior to Anterior)
*kcod = NIFTI_I2S (k axis is mostly Inferior to Superior)
\see "QUATERNION REPRESENTATION OF ROTATION MATRIX" in nifti1.h
\see nifti_quatern_to_mat44, nifti_mat44_to_quatern,
nifti_make_orthog_mat44
*//*-------------------------------------------------------------------------*/
void nifti_mat44_to_orientation( mat44 R , int *icod, int *jcod, int *kcod )
{
float xi,xj,xk , yi,yj,yk , zi,zj,zk , val,detQ,detP ;
mat33 P , Q , M ;
int i,j,k=0,p,q,r , ibest,jbest,kbest,pbest,qbest,rbest ;
float vbest ;
if( icod == NULL || jcod == NULL || kcod == NULL ) return ; /* bad */
*icod = *jcod = *kcod = 0 ; /* error returns, if sh*t happens */
/* load column vectors for each (i,j,k) direction from matrix */
/*-- i axis --*/ /*-- j axis --*/ /*-- k axis --*/
xi = R.m[0][0] ; xj = R.m[0][1] ; xk = R.m[0][2] ;
yi = R.m[1][0] ; yj = R.m[1][1] ; yk = R.m[1][2] ;
zi = R.m[2][0] ; zj = R.m[2][1] ; zk = R.m[2][2] ;
/* normalize column vectors to get unit vectors along each ijk-axis */
/* normalize i axis */
val = sqrt( xi*xi + yi*yi + zi*zi ) ;
if( val == 0.0 ) return ; /* stupid input */
xi /= val ; yi /= val ; zi /= val ;
/* normalize j axis */
val = sqrt( xj*xj + yj*yj + zj*zj ) ;
if( val == 0.0 ) return ; /* stupid input */
xj /= val ; yj /= val ; zj /= val ;
/* orthogonalize j axis to i axis, if needed */
val = xi*xj + yi*yj + zi*zj ; /* dot product between i and j */
if( fabs(val) > 1.e-4 ){
xj -= val*xi ; yj -= val*yi ; zj -= val*zi ;
val = sqrt( xj*xj + yj*yj + zj*zj ) ; /* must renormalize */
if( val == 0.0 ) return ; /* j was parallel to i? */
xj /= val ; yj /= val ; zj /= val ;
}
/* normalize k axis; if it is zero, make it the cross product i x j */
val = sqrt( xk*xk + yk*yk + zk*zk ) ;
if( val == 0.0 ){ xk = yi*zj-zi*yj; yk = zi*xj-zj*xi ; zk=xi*yj-yi*xj ; }
else { xk /= val ; yk /= val ; zk /= val ; }
/* orthogonalize k to i */
val = xi*xk + yi*yk + zi*zk ; /* dot product between i and k */
if( fabs(val) > 1.e-4 ){
xk -= val*xi ; yk -= val*yi ; zk -= val*zi ;
val = sqrt( xk*xk + yk*yk + zk*zk ) ;
if( val == 0.0 ) return ; /* bad */
xk /= val ; yk /= val ; zk /= val ;
}
/* orthogonalize k to j */
val = xj*xk + yj*yk + zj*zk ; /* dot product between j and k */
if( fabs(val) > 1.e-4 ){
xk -= val*xj ; yk -= val*yj ; zk -= val*zj ;
val = sqrt( xk*xk + yk*yk + zk*zk ) ;
if( val == 0.0 ) return ; /* bad */
xk /= val ; yk /= val ; zk /= val ;
}
Q.m[0][0] = xi ; Q.m[0][1] = xj ; Q.m[0][2] = xk ;
Q.m[1][0] = yi ; Q.m[1][1] = yj ; Q.m[1][2] = yk ;
Q.m[2][0] = zi ; Q.m[2][1] = zj ; Q.m[2][2] = zk ;
/* at this point, Q is the rotation matrix from the (i,j,k) to (x,y,z) axes */
detQ = nifti_mat33_determ( Q ) ;
if( detQ == 0.0 ) return ; /* shouldn't happen unless user is a DUFIS */
/* Build and test all possible +1/-1 coordinate permutation matrices P;
then find the P such that the rotation matrix M=PQ is closest to the
identity, in the sense of M having the smallest total rotation angle. */
/* Despite the formidable looking 6 nested loops, there are
only 3*3*3*2*2*2 = 216 passes, which will run very quickly. */
vbest = -666.0 ; ibest=pbest=qbest=rbest=1 ; jbest=2 ; kbest=3 ;
for( i=1 ; i <= 3 ; i++ ){ /* i = column number to use for row #1 */
for( j=1 ; j <= 3 ; j++ ){ /* j = column number to use for row #2 */
if( i == j ) continue ;
for( k=1 ; k <= 3 ; k++ ){ /* k = column number to use for row #3 */
if( i == k || j == k ) continue ;
P.m[0][0] = P.m[0][1] = P.m[0][2] =
P.m[1][0] = P.m[1][1] = P.m[1][2] =
P.m[2][0] = P.m[2][1] = P.m[2][2] = 0.0 ;
for( p=-1 ; p <= 1 ; p+=2 ){ /* p,q,r are -1 or +1 */
for( q=-1 ; q <= 1 ; q+=2 ){ /* and go into rows #1,2,3 */
for( r=-1 ; r <= 1 ; r+=2 ){
P.m[0][i-1] = p ; P.m[1][j-1] = q ; P.m[2][k-1] = r ;
detP = nifti_mat33_determ(P) ; /* sign of permutation */
if( detP * detQ <= 0.0 ) continue ; /* doesn't match sign of Q */
M = nifti_mat33_mul(P,Q) ;
/* angle of M rotation = 2.0*acos(0.5*sqrt(1.0+trace(M))) */
/* we want largest trace(M) == smallest angle == M nearest to I */
val = M.m[0][0] + M.m[1][1] + M.m[2][2] ; /* trace */
if( val > vbest ){
vbest = val ;
ibest = i ; jbest = j ; kbest = k ;
pbest = p ; qbest = q ; rbest = r ;
}
}}}}}}
/* At this point ibest is 1 or 2 or 3; pbest is -1 or +1; etc.
The matrix P that corresponds is the best permutation approximation
to Q-inverse; that is, P (approximately) takes (x,y,z) coordinates
to the (i,j,k) axes.
For example, the first row of P (which contains pbest in column ibest)
determines the way the i axis points relative to the anatomical
(x,y,z) axes. If ibest is 2, then the i axis is along the y axis,
which is direction P2A (if pbest > 0) or A2P (if pbest < 0).
So, using ibest and pbest, we can assign the output code for
the i axis. Mutatis mutandis for the j and k axes, of course. */
switch( ibest*pbest ){
case 1: i = NIFTI_L2R ; break ;
case -1: i = NIFTI_R2L ; break ;
case 2: i = NIFTI_P2A ; break ;
case -2: i = NIFTI_A2P ; break ;
case 3: i = NIFTI_I2S ; break ;
case -3: i = NIFTI_S2I ; break ;
}
switch( jbest*qbest ){
case 1: j = NIFTI_L2R ; break ;
case -1: j = NIFTI_R2L ; break ;
case 2: j = NIFTI_P2A ; break ;
case -2: j = NIFTI_A2P ; break ;
case 3: j = NIFTI_I2S ; break ;
case -3: j = NIFTI_S2I ; break ;
}
switch( kbest*rbest ){
case 1: k = NIFTI_L2R ; break ;
case -1: k = NIFTI_R2L ; break ;
case 2: k = NIFTI_P2A ; break ;
case -2: k = NIFTI_A2P ; break ;
case 3: k = NIFTI_I2S ; break ;
case -3: k = NIFTI_S2I ; break ;
}
*icod = i ; *jcod = j ; *kcod = k ; return ;
}
/*---------------------------------------------------------------------------*/
/* Routines to swap byte arrays in various ways:
- 2 at a time: ab -> ba [short]
- 4 at a time: abcd -> dcba [int, float]
- 8 at a time: abcdDCBA -> ABCDdcba [long long, double]
- 16 at a time: abcdefghHGFEDCBA -> ABCDEFGHhgfedcba [long double]
-----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------*/
/*! swap each byte pair from the given list of n pairs
*
* Due to alignment of structures at some architectures (e.g. on ARM),
* stick to char varaibles.
* Fixes http://bugs.debian.org/446893 Yaroslav
\return 0 if file looks like ANALYZE 7.5 [checks sizeof_hdr field == 348]
1 if file marked as NIFTI (header+data in 1 file)
2 if file marked as NIFTI (header+data in 2 files)
-1 if it can't tell, file doesn't exist, etc.
*//*------------------------------------------------------------------------*/
int is_nifti_file( const char *hname )
{
struct nifti_1_header nhdr ;
znzFile fp ;
int ii ;
char *tmpname;
/* bad input name? */
if( !nifti_validfilename(hname) ) return -1 ;
/* open file */
tmpname = nifti_findhdrname(hname);
if( tmpname == NULL ){
if( g_opts.debug > 0 )
fprintf(stderr,"** no header file found for '%s'\n",hname);
return -1;
}
fp = znzopen( tmpname , "rb" , nifti_is_gzfile(tmpname) ) ;
free(tmpname);
if (znz_isnull(fp)) return -1 ; /* bad open? */
/* read header, close file */
ii = (int)znzread( &nhdr , 1 , sizeof(nhdr) , fp ) ;
znzclose( fp ) ;
if( ii < (int) sizeof(nhdr) ) return -1 ; /* bad read? */
/* check for NIFTI-ness */
if( NIFTI_VERSION(nhdr) != 0 ){
return ( NIFTI_ONEFILE(nhdr) ) ? 1 : 2 ;
}
/* check for ANALYZE-ness (sizeof_hdr field == 348) */
ii = nhdr.sizeof_hdr ;
if( ii == (int)sizeof(nhdr) ) return 0 ; /* matches */
/* try byte-swapping header */
swap_4(ii) ;
if( ii == (int)sizeof(nhdr) ) return 0 ; /* matches */
return -1 ; /* not good */
}
static int print_hex_vals( const char * data, int nbytes, FILE * fp )
{
int c;
if ( !data || nbytes < 1 || !fp ) return -1;
fputs("0x", fp);
for ( c = 0; c < nbytes; c++ )
fprintf(fp, " %x", data[c]);
return 0;
}
/*----------------------------------------------------------------------*/
/*! display the contents of the nifti_1_header (send to stdout)
\param info if non-NULL, print this character string
\param hp pointer to nifti_1_header
*//*--------------------------------------------------------------------*/
int disp_nifti_1_header( const char * info, const nifti_1_header * hp )
{
int c;
fputs( "-------------------------------------------------------\n", stdout );
if ( info ) fputs( info, stdout );
if ( !hp ){ fputs(" ** no nifti_1_header to display!\n",stdout); return 1; }
fprintf(stdout," nifti_1_header :\n"
" sizeof_hdr = %d\n"
" data_type[10] = ", hp->sizeof_hdr);
print_hex_vals(hp->data_type, 10, stdout);
fprintf(stdout, "\n"
" db_name[18] = ");
print_hex_vals(hp->db_name, 18, stdout);
fprintf(stdout, "\n"
" extents = %d\n"
" session_error = %d\n"
" regular = 0x%x\n"
" dim_info = 0x%x\n",
hp->extents, hp->session_error, hp->regular, hp->dim_info );
fprintf(stdout, " dim[8] =");
for ( c = 0; c < 8; c++ ) fprintf(stdout," %d", hp->dim[c]);
fprintf(stdout, "\n"
" intent_p1 = %f\n"
" intent_p2 = %f\n"
" intent_p3 = %f\n"
" intent_code = %d\n"
" datatype = %d\n"
" bitpix = %d\n"
" slice_start = %d\n"
" pixdim[8] =",
hp->intent_p1, hp->intent_p2, hp->intent_p3, hp->intent_code,
hp->datatype, hp->bitpix, hp->slice_start);
/* break pixdim over 2 lines */
for ( c = 0; c < 4; c++ ) fprintf(stdout," %f", hp->pixdim[c]);
fprintf(stdout, "\n ");
for ( c = 4; c < 8; c++ ) fprintf(stdout," %f", hp->pixdim[c]);
fprintf(stdout, "\n"
" vox_offset = %f\n"
" scl_slope = %f\n"
" scl_inter = %f\n"
" slice_end = %d\n"
" slice_code = %d\n"
" xyzt_units = 0x%x\n"
" cal_max = %f\n"
" cal_min = %f\n"
" slice_duration = %f\n"
" toffset = %f\n"
" glmax = %d\n"
" glmin = %d\n",
hp->vox_offset, hp->scl_slope, hp->scl_inter, hp->slice_end,
hp->slice_code, hp->xyzt_units, hp->cal_max, hp->cal_min,
hp->slice_duration, hp->toffset, hp->glmax, hp->glmin);
fprintf(stdout,
" descrip = '%.80s'\n"
" aux_file = '%.24s'\n"
" qform_code = %d\n"
" sform_code = %d\n"
" quatern_b = %f\n"
" quatern_c = %f\n"
" quatern_d = %f\n"
" qoffset_x = %f\n"
" qoffset_y = %f\n"
" qoffset_z = %f\n"
" srow_x[4] = %f, %f, %f, %f\n"
" srow_y[4] = %f, %f, %f, %f\n"
" srow_z[4] = %f, %f, %f, %f\n"
" intent_name = '%-.16s'\n"
" magic = '%-.4s'\n",
hp->descrip, hp->aux_file, hp->qform_code, hp->sform_code,
hp->quatern_b, hp->quatern_c, hp->quatern_d,
hp->qoffset_x, hp->qoffset_y, hp->qoffset_z,
hp->srow_x[0], hp->srow_x[1], hp->srow_x[2], hp->srow_x[3],
hp->srow_y[0], hp->srow_y[1], hp->srow_y[2], hp->srow_y[3],
hp->srow_z[0], hp->srow_z[1], hp->srow_z[2], hp->srow_z[3],
hp->intent_name, hp->magic);
fputs( "-------------------------------------------------------\n", stdout );
fflush(stdout);
return 0;
}
#undef ERREX
#define ERREX(msg) \
do{ fprintf(stderr,"** ERROR: nifti_convert_nhdr2nim: %s\n", (msg) ) ; \
return NULL ; } while(0)
/*----------------------------------------------------------------------*/
/*! convert a nifti_1_header into a nift1_image
\return an allocated nifti_image, or NULL on failure
*//*--------------------------------------------------------------------*/
nifti_image* nifti_convert_nhdr2nim(struct nifti_1_header nhdr,
const char * fname)
{
int ii , doswap , ioff ;
int is_nifti , is_onefile ;
nifti_image *nim;
nim = (nifti_image *)calloc( 1 , sizeof(nifti_image) ) ;
if( !nim ) ERREX("failed to allocate nifti image");
/* be explicit with pointers */
nim->fname = NULL;
nim->iname = NULL;
nim->data = NULL;
/**- check if we must swap bytes */
doswap = need_nhdr_swap(nhdr.dim[0], nhdr.sizeof_hdr); /* swap data flag */
if( doswap < 0 ){
if( doswap == -1 ) ERREX("bad dim[0]") ;
ERREX("bad sizeof_hdr") ; /* else */
}
/**- determine if this is a NIFTI-1 compliant header */
is_nifti = NIFTI_VERSION(nhdr) ;
/*
* before swapping header, record the Analyze75 orient code
*/
if(!is_nifti)
{
/**- in analyze75, the orient code is at the same address as
* qform_code, but it's just one byte
* the qform_code will be zero, at which point you can check
* analyze75_orient if you care to.
*/
unsigned char c = *((char *)(&nhdr.qform_code));
nim->analyze75_orient = (analyze_75_orient_code)c;
}
if( doswap ) {
if ( g_opts.debug > 3 ) disp_nifti_1_header("-d ni1 pre-swap: ", &nhdr);
swap_nifti_header( &nhdr , is_nifti ) ;
}
if ( g_opts.debug > 2 ) disp_nifti_1_header("-d nhdr2nim : ", &nhdr);
if( nhdr.datatype == DT_BINARY ||
nhdr.datatype == DT_UNKNOWN ) ERREX("bad datatype") ;
if( nhdr.dim[1] <= 0 ) ERREX("bad dim[1]") ;
/* fix bad dim[] values in the defined dimension range */
for( ii=2 ; ii <= nhdr.dim[0] ; ii++ )
if( nhdr.dim[ii] <= 0 ) nhdr.dim[ii] = 1 ;
/* fix any remaining bad dim[] values, so garbage does not propagate */
/* (only values 0 or 1 seem rational, otherwise set to arbirary 1) */
for( ii=nhdr.dim[0]+1 ; ii <= 7 ; ii++ )
if( nhdr.dim[ii] != 1 && nhdr.dim[ii] != 0) nhdr.dim[ii] = 1 ;
#if 0 /* rely on dim[0], do not attempt to modify it 16 Nov 2005 [rickr] */
/**- get number of dimensions (ignoring dim[0] now) */
for( ii=7 ; ii >= 2 ; ii-- ) /* loop backwards until we */
if( nhdr.dim[ii] > 1 ) break ; /* find a dim bigger than 1 */
ndim = ii ;
#endif
/**- set bad grid spacings to 1.0 */
for( ii=1 ; ii <= nhdr.dim[0] ; ii++ )
if( nhdr.pixdim[ii] == 0.0 )
nhdr.pixdim[ii] = 1.0;
is_onefile = is_nifti && NIFTI_ONEFILE(nhdr) ;
if( is_nifti ) nim->nifti_type = (is_onefile) ? NIFTI_FTYPE_NIFTI1_1
: NIFTI_FTYPE_NIFTI1_2 ;
else nim->nifti_type = NIFTI_FTYPE_ANALYZE ;
ii = nifti_short_order() ;
if( doswap ) nim->byteorder = REVERSE_ORDER(ii) ;
else nim->byteorder = ii ;
/**- set dimensions of data array */
nim->ndim = nim->dim[0] = nhdr.dim[0];
nim->nx = nim->dim[1] = nhdr.dim[1];
nim->ny = nim->dim[2] = nhdr.dim[2];
nim->nz = nim->dim[3] = nhdr.dim[3];
nim->nt = nim->dim[4] = nhdr.dim[4];
nim->nu = nim->dim[5] = nhdr.dim[5];
nim->nv = nim->dim[6] = nhdr.dim[6];
nim->nw = nim->dim[7] = nhdr.dim[7];
for( ii=1, nim->nvox=1; ii <= nhdr.dim[0]; ii++ )
nim->nvox *= nhdr.dim[ii];
/**- set the type of data in voxels and how many bytes per voxel */
nim->datatype = nhdr.datatype ;
nifti_datatype_sizes( nim->datatype , &(nim->nbyper) , &(nim->swapsize) ) ;
if( nim->nbyper == 0 ){ free(nim); ERREX("bad datatype"); }
/**- set the grid spacings */
nim->dx = nim->pixdim[1] = nhdr.pixdim[1] ;
nim->dy = nim->pixdim[2] = nhdr.pixdim[2] ;
nim->dz = nim->pixdim[3] = nhdr.pixdim[3] ;
nim->dt = nim->pixdim[4] = nhdr.pixdim[4] ;
nim->du = nim->pixdim[5] = nhdr.pixdim[5] ;
nim->dv = nim->pixdim[6] = nhdr.pixdim[6] ;
nim->dw = nim->pixdim[7] = nhdr.pixdim[7] ;
/**- compute qto_xyz transformation from pixel indexes (i,j,k) to (x,y,z) */
if( !is_nifti || nhdr.qform_code <= 0 ){
/**- if not nifti or qform_code <= 0, use grid spacing for qto_xyz */
nim->qto_xyz.m[0][0] = nim->dx ; /* grid spacings */
nim->qto_xyz.m[1][1] = nim->dy ; /* along diagonal */
nim->qto_xyz.m[2][2] = nim->dz ;
/* off diagonal is zero */
nim->qto_xyz.m[0][1]=nim->qto_xyz.m[0][2]=nim->qto_xyz.m[0][3] = 0.0;
nim->qto_xyz.m[1][0]=nim->qto_xyz.m[1][2]=nim->qto_xyz.m[1][3] = 0.0;
nim->qto_xyz.m[2][0]=nim->qto_xyz.m[2][1]=nim->qto_xyz.m[2][3] = 0.0;
/* last row is always [ 0 0 0 1 ] */
nim->qto_xyz.m[3][0]=nim->qto_xyz.m[3][1]=nim->qto_xyz.m[3][2] = 0.0;
nim->qto_xyz.m[3][3]= 1.0 ;
nim->qform_code = NIFTI_XFORM_UNKNOWN ;
if( g_opts.debug > 1 ) fprintf(stderr,"-d no qform provided\n");
} else {
/**- else NIFTI: use the quaternion-specified transformation */
nim->quatern_b = nhdr.quatern_b;
nim->quatern_c = nhdr.quatern_c;
nim->quatern_d = nhdr.quatern_d;
nim->qoffset_x = nhdr.qoffset_x;
nim->qoffset_y = nhdr.qoffset_y;
nim->qoffset_z = nhdr.qoffset_z;
nim->qfac = (nhdr.pixdim[0] < 0.0) ? -1.0 : 1.0 ; /* left-handedness? */
nim->qto_xyz = nifti_quatern_to_mat44(
nim->quatern_b, nim->quatern_c, nim->quatern_d,
nim->qoffset_x, nim->qoffset_y, nim->qoffset_z,
nim->dx , nim->dy , nim->dz ,
nim->qfac ) ;
nim->qform_code = nhdr.qform_code ;
if( g_opts.debug > 1 )
nifti_disp_matrix_orient("-d qform orientations:\n", nim->qto_xyz);
}
/**- load inverse transformation (x,y,z) -> (i,j,k) */
nim->qto_ijk = nifti_mat44_inverse( nim->qto_xyz ) ;
/**- load sto_xyz affine transformation, if present */
if( !is_nifti || nhdr.sform_code <= 0 ){
/**- if not nifti or sform_code <= 0, then no sto transformation */
nim->sform_code = NIFTI_XFORM_UNKNOWN ;
if( g_opts.debug > 1 ) fprintf(stderr,"-d no sform provided\n");
} else {
/**- else set the sto transformation from srow_*[] */
nim->sto_xyz.m[0][0] = nhdr.srow_x[0] ;
nim->sto_xyz.m[0][1] = nhdr.srow_x[1] ;
nim->sto_xyz.m[0][2] = nhdr.srow_x[2] ;
nim->sto_xyz.m[0][3] = nhdr.srow_x[3] ;
nim->sto_xyz.m[1][0] = nhdr.srow_y[0] ;
nim->sto_xyz.m[1][1] = nhdr.srow_y[1] ;
nim->sto_xyz.m[1][2] = nhdr.srow_y[2] ;
nim->sto_xyz.m[1][3] = nhdr.srow_y[3] ;
nim->sto_xyz.m[2][0] = nhdr.srow_z[0] ;
nim->sto_xyz.m[2][1] = nhdr.srow_z[1] ;
nim->sto_xyz.m[2][2] = nhdr.srow_z[2] ;
nim->sto_xyz.m[2][3] = nhdr.srow_z[3] ;
/* last row is always [ 0 0 0 1 ] */
nim->sto_xyz.m[3][0]=nim->sto_xyz.m[3][1]=nim->sto_xyz.m[3][2] = 0.0;
nim->sto_xyz.m[3][3]= 1.0 ;
nim->sto_ijk = nifti_mat44_inverse( nim->sto_xyz ) ;
nim->sform_code = nhdr.sform_code ;
if( g_opts.debug > 1 )
nifti_disp_matrix_orient("-d sform orientations:\n", nim->sto_xyz);
}
/**- set miscellaneous NIFTI stuff */
if( is_nifti ){
nim->scl_slope = nhdr.scl_slope;
nim->scl_inter = nhdr.scl_inter;
nim->intent_code = nhdr.intent_code ;
nim->intent_p1 = nhdr.intent_p1;
nim->intent_p2 = nhdr.intent_p2;
nim->intent_p3 = nhdr.intent_p3;
nim->toffset = nhdr.toffset;
memcpy(nim->intent_name,nhdr.intent_name,15); nim->intent_name[15] = '\0';
nim->xyz_units = XYZT_TO_SPACE(nhdr.xyzt_units) ;
nim->time_units = XYZT_TO_TIME (nhdr.xyzt_units) ;
nim->freq_dim = DIM_INFO_TO_FREQ_DIM ( nhdr.dim_info ) ;
nim->phase_dim = DIM_INFO_TO_PHASE_DIM( nhdr.dim_info ) ;
nim->slice_dim = DIM_INFO_TO_SLICE_DIM( nhdr.dim_info ) ;
nim->slice_code = nhdr.slice_code ;
nim->slice_start = nhdr.slice_start ;
nim->slice_end = nhdr.slice_end ;
nim->slice_duration = nhdr.slice_duration;
}
/**- set Miscellaneous ANALYZE stuff */
nim->cal_min = nhdr.cal_min;
nim->cal_max = nhdr.cal_max;
memcpy(nim->descrip ,nhdr.descrip ,79) ; nim->descrip [79] = '\0' ;
memcpy(nim->aux_file,nhdr.aux_file,23) ; nim->aux_file[23] = '\0' ;
/**- set ioff from vox_offset (but at least sizeof(header)) */
is_onefile = is_nifti && NIFTI_ONEFILE(nhdr) ;
if( is_onefile ){
ioff = (int)nhdr.vox_offset ;
if( ioff < (int) sizeof(nhdr) ) ioff = (int) sizeof(nhdr) ;
} else {
ioff = (int)nhdr.vox_offset ;
}
nim->iname_offset = ioff ;
/**- deal with file names if set */
if (fname!=NULL) {
nifti_set_filenames(nim,fname,0,0);
if (nim->iname==NULL) { ERREX("bad filename"); }
} else {
nim->fname = NULL;
nim->iname = NULL;
}
/* clear extension fields */
nim->num_ext = 0;
nim->ext_list = NULL;
return nim;
}
#undef ERREX
#define ERREX(msg) \
do{ fprintf(stderr,"** ERROR: nifti_image_open(%s): %s\n", \
(hname != NULL) ? hname : "(null)" , (msg) ) ; \
return fptr ; } while(0)
/***************************************************************
* nifti_image_open
***************************************************************/
/*! znzFile nifti_image_open( char *hname, char *opts , nifti_image **nim)
\brief Read in NIFTI-1 or ANALYZE-7.5 file (pair) header information into a nifti_image struct.
- The image data is not read from disk (it may be read later using
nifti_image_load(), for example).
- The image data will be stored in whatever data format the
input data is; no scaling will be applied.
- DT_BINARY data is not supported.
- nifti_image_free() can be used to delete the returned struct,
when you are done with it.
\param hname filename of dataset .hdr or .nii file
\param opts options string for opening the header file
\param nim pointer to pointer to nifti_image struct
(this routine allocates the nifti_image struct)
\return file pointer (gzippable) to the file with the image data,
ready for reading.
nifti_1_header my_header;
my_header = nifti_convert_nim2nhdr(my_nim_pointer);
*//*--------------------------------------------------------------------*/
struct nifti_1_header nifti_convert_nim2nhdr(const nifti_image * nim)
{
struct nifti_1_header nhdr;
memset(&nhdr,0,sizeof(nhdr)) ; /* zero out header, to be safe */
/**- load the ANALYZE-7.5 generic parts of the header struct */
nhdr.sizeof_hdr = sizeof(nhdr) ;
nhdr.regular = 'r' ; /* for some stupid reason */
nhdr.dim[0] = nim->ndim ;
nhdr.dim[1] = nim->nx ; nhdr.dim[2] = nim->ny ; nhdr.dim[3] = nim->nz ;
nhdr.dim[4] = nim->nt ; nhdr.dim[5] = nim->nu ; nhdr.dim[6] = nim->nv ;
nhdr.dim[7] = nim->nw ;
nhdr.pixdim[0] = 0.0 ;
nhdr.pixdim[1] = nim->dx ; nhdr.pixdim[2] = nim->dy ;
nhdr.pixdim[3] = nim->dz ; nhdr.pixdim[4] = nim->dt ;
nhdr.pixdim[5] = nim->du ; nhdr.pixdim[6] = nim->dv ;
nhdr.pixdim[7] = nim->dw ;
nhdr.datatype = nim->datatype ;
nhdr.bitpix = 8 * nim->nbyper ;
if( nim->cal_max > nim->cal_min ){
nhdr.cal_max = nim->cal_max ;
nhdr.cal_min = nim->cal_min ;
}
if( nim->scl_slope != 0.0 ){
nhdr.scl_slope = nim->scl_slope ;
nhdr.scl_inter = nim->scl_inter ;
}
if( nim->descrip[0] != '\0' ){
memcpy(nhdr.descrip ,nim->descrip ,79) ; nhdr.descrip[79] = '\0' ;
}
if( nim->aux_file[0] != '\0' ){
memcpy(nhdr.aux_file ,nim->aux_file ,23) ; nhdr.aux_file[23] = '\0' ;
}
/**- Load NIFTI specific stuff into the header */
if( nim->nifti_type > NIFTI_FTYPE_ANALYZE ){ /* then not ANALYZE */
if( nim->nifti_type == NIFTI_FTYPE_NIFTI1_1 ) strcpy(nhdr.magic,"n+1") ;
else strcpy(nhdr.magic,"ni1") ;
nhdr.pixdim[1] = fabs(nhdr.pixdim[1]) ; nhdr.pixdim[2] = fabs(nhdr.pixdim[2]) ;
nhdr.pixdim[3] = fabs(nhdr.pixdim[3]) ; nhdr.pixdim[4] = fabs(nhdr.pixdim[4]) ;
nhdr.pixdim[5] = fabs(nhdr.pixdim[5]) ; nhdr.pixdim[6] = fabs(nhdr.pixdim[6]) ;
nhdr.pixdim[7] = fabs(nhdr.pixdim[7]) ;
nhdr.intent_code = nim->intent_code ;
nhdr.intent_p1 = nim->intent_p1 ;
nhdr.intent_p2 = nim->intent_p2 ;
nhdr.intent_p3 = nim->intent_p3 ;
if( nim->intent_name[0] != '\0' ){
memcpy(nhdr.intent_name,nim->intent_name,15) ;
nhdr.intent_name[15] = '\0' ;
}
nhdr.vox_offset = (float) nim->iname_offset ;
nhdr.xyzt_units = SPACE_TIME_TO_XYZT( nim->xyz_units, nim->time_units ) ;
nhdr.toffset = nim->toffset ;
if( nim->qform_code > 0 ){
nhdr.qform_code = nim->qform_code ;
nhdr.quatern_b = nim->quatern_b ;
nhdr.quatern_c = nim->quatern_c ;
nhdr.quatern_d = nim->quatern_d ;
nhdr.qoffset_x = nim->qoffset_x ;
nhdr.qoffset_y = nim->qoffset_y ;
nhdr.qoffset_z = nim->qoffset_z ;
nhdr.pixdim[0] = (nim->qfac >= 0.0) ? 1.0 : -1.0 ;
}
if( nim->sform_code > 0 ){
nhdr.sform_code = nim->sform_code ;
nhdr.srow_x[0] = nim->sto_xyz.m[0][0] ;
nhdr.srow_x[1] = nim->sto_xyz.m[0][1] ;
nhdr.srow_x[2] = nim->sto_xyz.m[0][2] ;
nhdr.srow_x[3] = nim->sto_xyz.m[0][3] ;
nhdr.srow_y[0] = nim->sto_xyz.m[1][0] ;
nhdr.srow_y[1] = nim->sto_xyz.m[1][1] ;
nhdr.srow_y[2] = nim->sto_xyz.m[1][2] ;
nhdr.srow_y[3] = nim->sto_xyz.m[1][3] ;
nhdr.srow_z[0] = nim->sto_xyz.m[2][0] ;
nhdr.srow_z[1] = nim->sto_xyz.m[2][1] ;
nhdr.srow_z[2] = nim->sto_xyz.m[2][2] ;
nhdr.srow_z[3] = nim->sto_xyz.m[2][3] ;
}
nhdr.dim_info = FPS_INTO_DIM_INFO( nim->freq_dim ,
nim->phase_dim , nim->slice_dim ) ;
nhdr.slice_code = nim->slice_code ;
nhdr.slice_start = nim->slice_start ;
nhdr.slice_end = nim->slice_end ;
nhdr.slice_duration = nim->slice_duration ;
}
return nhdr;
}
/*----------------------------------------------------------------------*/
/*! \fn int nifti_copy_extensions(nifti_image * nim_dest, nifti_image * nim_src)
\brief copy the nifti1_extension list from src to dest
Duplicate the list of nifti1_extensions. The dest structure must
be clear of extensions.
\return 0 on success, -1 on failure
\sa nifti_add_extension, nifti_free_extensions
*/
int nifti_copy_extensions(nifti_image * nim_dest, const nifti_image * nim_src)
{
char * data;
size_t bytes;
int c, size, old_size;
if( nim_dest->num_ext > 0 || nim_dest->ext_list != NULL ){
fprintf(stderr,"** will not copy extensions over existing ones\n");
return -1;
}
if( g_opts.debug > 1 )
fprintf(stderr,"+d duplicating %d extension(s)\n", nim_src->num_ext);
if( nim_src->num_ext <= 0 ) return 0;
bytes = nim_src->num_ext * sizeof(nifti1_extension); /* I'm lazy */
nim_dest->ext_list = (nifti1_extension *)malloc(bytes);
if( !nim_dest->ext_list ){
fprintf(stderr,"** failed to allocate %d nifti1_extension structs\n",
nim_src->num_ext);
return -1;
}
/* copy the extension data */
nim_dest->num_ext = 0;
for( c = 0; c < nim_src->num_ext; c++ ){
size = old_size = nim_src->ext_list[c].esize;
if( size & 0xf ) size = (size + 0xf) & ~0xf; /* make multiple of 16 */
if( g_opts.debug > 2 )
fprintf(stderr,"+d dup'ing ext #%d of size %d (from size %d)\n",
c, size, old_size);
/* data length is size-8, as esize includes space for esize and ecode */
data = (char *)calloc(size-8,sizeof(char)); /* maybe size > old */
if( !data ){
fprintf(stderr,"** failed to alloc %d bytes for extention\n", size);
if( c == 0 ) { free(nim_dest->ext_list); nim_dest->ext_list = NULL; }
/* otherwise, keep what we have (a.o.t. deleting them all) */
return -1;
}
/* finally, fill the new structure */
nim_dest->ext_list[c].esize = size;
nim_dest->ext_list[c].ecode = nim_src->ext_list[c].ecode;
nim_dest->ext_list[c].edata = data;
memcpy(data, nim_src->ext_list[c].edata, old_size-8);
nim_dest->num_ext++;
}
return 0;
}
/*----------------------------------------------------------------------*/
/*! compute the total size of all extensions
\return the total of all esize fields
Note that each esize includes 4 bytes for ecode, 4 bytes for esize,
and the bytes used for the data. Each esize also needs to be a
multiple of 16, so it may be greater than the sum of its 3 parts.
*//*--------------------------------------------------------------------*/
int nifti_extension_size(nifti_image *nim)
{
int c, size = 0;
if( !nim || nim->num_ext <= 0 ) return 0;
if( g_opts.debug > 2 ) fprintf(stderr,"-d ext sizes:");
for ( c = 0; c < nim->num_ext; c++ ){
size += nim->ext_list[c].esize;
if( g_opts.debug > 2 ) fprintf(stderr," %d",nim->ext_list[c].esize);
}
if( g_opts.debug > 2 ) fprintf(stderr," (total = %d)\n",size);
return size;
}
/*----------------------------------------------------------------------*/
/*! set the nifti_image iname_offset field, based on nifti_type
- if writing to 2 files, set offset to 0
- if writing to a single NIFTI-1 file, set the offset to
352 + total extension size, then align to 16-byte boundary
- if writing an ASCII header, set offset to -1
*//*--------------------------------------------------------------------*/
void nifti_set_iname_offset(nifti_image *nim)
{
int offset;
switch( nim->nifti_type ){
default: /* writing into 2 files */
/* we only write files with 0 offset in the 2 file format */
nim->iname_offset = 0 ;
break ;
/* NIFTI-1 single binary file - always update */
case NIFTI_FTYPE_NIFTI1_1:
offset = nifti_extension_size(nim)+sizeof(struct nifti_1_header)+4;
/* be sure offset is aligned to a 16 byte boundary */
if ( ( offset % 16 ) != 0 ) offset = ((offset + 0xf) & ~0xf);
if( nim->iname_offset != offset ){
if( g_opts.debug > 1 )
fprintf(stderr,"+d changing offset from %d to %d\n",
nim->iname_offset, offset);
nim->iname_offset = offset;
}
break ;
/* non-standard case: NIFTI-1 ASCII header + binary data (single file) */
case NIFTI_FTYPE_ASCII:
nim->iname_offset = -1 ; /* compute offset from filesize */
break ;
}
}
/*----------------------------------------------------------------------*/
/*! write the nifti_image dataset to disk, optionally including data
This is just a front-end for nifti_image_write_hdr_img2.
\param nim nifti_image to write to disk
\param write_data write options (see nifti_image_write_hdr_img2)
\param opts file open options ("wb" from nifti_image_write)
\sa nifti_image_write, nifti_image_write_hdr_img2, nifti_image_free,
nifti_set_filenames
*//*--------------------------------------------------------------------*/
znzFile nifti_image_write_hdr_img( nifti_image *nim , int write_data ,
const char* opts )
{
return nifti_image_write_hdr_img2(nim,write_data,opts,NULL,NULL);
}
#undef ERREX
#define ERREX(msg) \
do{ fprintf(stderr,"** ERROR: nifti_image_write_hdr_img: %s\n",(msg)) ; \
return fp ; } while(0)
/* ----------------------------------------------------------------------*/
/*! This writes the header (and optionally the image data) to file
*
* If the image data file is left open it returns a valid znzFile handle.
* It also uses imgfile as the open image file is not null, and modifies
* it inside.
*
* \param nim nifti_image to write to disk
* \param write_opts flags whether to write data and/or close file (see below)
* \param opts file-open options, probably "wb" from nifti_image_write()
* \param imgfile optional open znzFile struct, for writing image data
(may be NULL)
* \param NBL optional nifti_brick_list, containing the image data
(may be NULL)
*
* Values for write_opts mode are based on two binary flags
* ( 0/1 for no-write/write data, and 0/2 for close/leave-open files ) :
* - 0 = do not write data and close (do not open data file)
* - 1 = write data and close
* - 2 = do not write data and leave data file open
* - 3 = write data and leave data file open
*
* \sa nifti_image_write, nifti_image_write_hdr_img, nifti_image_free,
* nifti_set_filenames
*//*---------------------------------------------------------------------*/
znzFile nifti_image_write_hdr_img2(nifti_image *nim, int write_opts,
const char * opts, znzFile imgfile, const nifti_brick_list * NBL)
{
struct nifti_1_header nhdr ;
znzFile fp=NULL;
size_t ss ;
int write_data, leave_open;
char func[] = { "nifti_image_write_hdr_img2" };
write_data = write_opts & 1; /* just separate the bits now */
leave_open = write_opts & 2;
if( ! nim ) ERREX("NULL input") ;
if( ! nifti_validfilename(nim->fname) ) ERREX("bad fname input") ;
if( write_data && ! nim->data && ! NBL ) ERREX("no image data") ;
if( write_data && NBL && ! nifti_NBL_matches_nim(nim, NBL) )
ERREX("NBL does not match nim");
nifti_set_iname_offset(nim);
if( g_opts.debug > 1 ){
fprintf(stderr,"-d writing nifti file '%s'...\n", nim->fname);
if( g_opts.debug > 2 )
fprintf(stderr,"-d nifti type %d, offset %d\n",
nim->nifti_type, nim->iname_offset);
}
if( nim->nifti_type == NIFTI_FTYPE_ASCII ) /* non-standard case */
return nifti_write_ascii_image(nim,NBL,opts,write_data,leave_open);
nhdr = nifti_convert_nim2nhdr(nim); /* create the nifti1_header struct */
/* if writing to 2 files, make sure iname is set and different from fname */
if( nim->nifti_type != NIFTI_FTYPE_NIFTI1_1 ){
if( nim->iname && strcmp(nim->iname,nim->fname) == 0 ){
free(nim->iname) ; nim->iname = NULL ;
}
if( nim->iname == NULL ){ /* then make a new one */
nim->iname = nifti_makeimgname(nim->fname,nim->nifti_type,0,0);
if( nim->iname == NULL ) return NULL;
}
}
/* if we have an imgfile and will write the header there, use it */
if( ! znz_isnull(imgfile) && nim->nifti_type == NIFTI_FTYPE_NIFTI1_1 ){
if( g_opts.debug > 2 ) fprintf(stderr,"+d using passed file for hdr\n");
fp = imgfile;
}
else {
if( g_opts.debug > 2 )
fprintf(stderr,"+d opening output file %s [%s]\n",nim->fname,opts);
fp = znzopen( nim->fname , opts , nifti_is_gzfile(nim->fname) ) ;
if( znz_isnull(fp) ){
LNI_FERR(func,"cannot open output file",nim->fname);
return fp;
}
}
/* write the header and extensions */
ss = znzwrite(&nhdr , 1 , sizeof(nhdr) , fp); /* write header */
if( ss < sizeof(nhdr) ){
LNI_FERR(func,"bad header write to output file",nim->fname);
znzclose(fp); return fp;
}
/* partial file exists, and errors have been printed, so ignore return */
if( nim->nifti_type != NIFTI_FTYPE_ANALYZE )
(void)nifti_write_extensions(fp,nim);
/* if the header is all we want, we are done */
if( ! write_data && ! leave_open ){
if( g_opts.debug > 2 ) fprintf(stderr,"-d header is all we want: done\n");
znzclose(fp); return(fp);
}
if( nim->nifti_type != NIFTI_FTYPE_NIFTI1_1 ){ /* get a new file pointer */
znzclose(fp); /* first, close header file */
if( ! znz_isnull(imgfile) ){
if(g_opts.debug > 2) fprintf(stderr,"+d using passed file for img\n");
fp = imgfile;
}
else {
if( g_opts.debug > 2 )
fprintf(stderr,"+d opening img file '%s'\n", nim->iname);
fp = znzopen( nim->iname , opts , nifti_is_gzfile(nim->iname) ) ;
if( znz_isnull(fp) ) ERREX("cannot open image file") ;
}
}
znzseek(fp, nim->iname_offset, SEEK_SET); /* in any case, seek to offset */
if( write_data ) nifti_write_all_data(fp,nim,NBL);
if( ! leave_open ) znzclose(fp);
return fp;
}
/*----------------------------------------------------------------------*/
/*! write a nifti_image to disk in ASCII format
*//*--------------------------------------------------------------------*/
znzFile nifti_write_ascii_image(nifti_image *nim, const nifti_brick_list * NBL,
const char *opts, int write_data, int leave_open)
{
znzFile fp;
char * hstr;
hstr = nifti_image_to_ascii( nim ) ; /* get header in ASCII form */
if( ! hstr ){ fprintf(stderr,"** failed image_to_ascii()\n"); return NULL; }
fp = znzopen( nim->fname , opts , nifti_is_gzfile(nim->fname) ) ;
if( znz_isnull(fp) ){
free(hstr);
fprintf(stderr,"** failed to open '%s' for ascii write\n",nim->fname);
return fp;
}
znzputs(hstr,fp); /* header */
nifti_write_extensions(fp,nim); /* extensions */
if ( write_data ) { nifti_write_all_data(fp,nim,NBL); } /* data */
if ( ! leave_open ) { znzclose(fp); }
free(hstr);
return fp; /* returned but may be closed */
}
/*--------------------------------------------------------------------------*/
/*! Write a nifti_image to disk.
Since data is properly byte-swapped upon reading, it is assumed
to be in the byte-order of the current CPU at write time. Thus,
nim->byte_order should match that of the current CPU. Note that
the nifti_set_filenames() function takes the flag, set_byte_order.
The following fields of nim affect how the output appears:
- nifti_type = 0 ==> ANALYZE-7.5 format file pair will be written
- nifti_type = 1 ==> NIFTI-1 format single file will be written
(data offset will be 352+extensions)
- nifti_type = 2 ==> NIFTI_1 format file pair will be written
- nifti_type = 3 ==> NIFTI_1 ASCII single file will be written
- fname is the name of the output file (header or header+data)
- if a file pair is being written, iname is the name of the data file
- existing files WILL be overwritten with extreme prejudice
- if qform_code > 0, the quatern_*, qoffset_*, and qfac fields determine
the qform output, NOT the qto_xyz matrix; if you want to compute these
fields from the qto_xyz matrix, you can use the utility function
nifti_mat44_to_quatern()
\sa nifti_image_write_bricks, nifti_image_free, nifti_set_filenames,
nifti_image_write_hdr_img
*//*------------------------------------------------------------------------*/
void nifti_image_write( nifti_image *nim )
{
znzFile fp = nifti_image_write_hdr_img(nim,1,"wb");
if( fp ){
if( g_opts.debug > 2 ) fprintf(stderr,"-d niw: done with znzFile\n");
free(fp);
}
if( g_opts.debug > 1 ) fprintf(stderr,"-d nifti_image_write: done\n");
}
/*----------------------------------------------------------------------*/
/*! similar to nifti_image_write, but data is in NBL struct, not nim->data
\sa nifti_image_write, nifti_image_free, nifti_set_filenames, nifti_free_NBL
*//*--------------------------------------------------------------------*/
void nifti_image_write_bricks( nifti_image *nim, const nifti_brick_list * NBL )
{
znzFile fp = nifti_image_write_hdr_img2(nim,1,"wb",NULL,NBL);
if( fp ){
if( g_opts.debug > 2 ) fprintf(stderr,"-d niwb: done with znzFile\n");
free(fp);
}
if( g_opts.debug > 1 ) fprintf(stderr,"-d niwb: done writing bricks\n");
}
/*----------------------------------------------------------------------*/
/*! copy the nifti_image structure, without data
Duplicate the structure, including fname, iname and extensions.
Leave the data pointer as NULL.
*//*--------------------------------------------------------------------*/
nifti_image * nifti_copy_nim_info(const nifti_image * src)
{
nifti_image *dest;
dest = (nifti_image *)calloc(1,sizeof(nifti_image));
if( !dest ){
fprintf(stderr,"** NCNI: failed to alloc nifti_image\n");
return NULL;
}
memcpy(dest, src, sizeof(nifti_image));
if( src->fname ) dest->fname = nifti_strdup(src->fname);
if( src->iname ) dest->iname = nifti_strdup(src->iname);
dest->num_ext = 0;
dest->ext_list = NULL;
/* errors will be printed in NCE(), continue in either case */
(void)nifti_copy_extensions(dest, src);
dest->data = NULL;
return dest;
}
/*------------------------------------------------------------------------*/
/* Un-escape a C string in place -- that is, convert XML escape sequences
back into their characters. (This can be done in place since the
replacement is always smaller than the input.) Escapes recognized are:
- < -> <
- > -> >
- " -> "
- ' -> '
- & -> &
Also replace CR LF pair (Microsoft), or CR alone (Macintosh) with
LF (Unix), per the XML standard.
Return value is number of replacements made (if you care).
--------------------------------------------------------------------------*/
#undef CR
#undef LF
#define CR 0x0D
#define LF 0x0A
static int unescape_string( char *str )
{
int ii,jj , nn,ll ;
if( str == NULL ) return 0 ; /* no string? */
ll = (int)strlen(str) ; if( ll == 0 ) return 0 ;
/* scan for escapes: &something; */
for( ii=jj=nn=0 ; ii
This function may be used to read parts of a nifti dataset, such as
the time series for a single voxel, or perhaps a slice. It is similar
to nifti_image_load(), though the passed 'data' parameter is used for
returning the image, not nim->data.
\param nim given nifti_image struct, corresponding to the data file
\param dims given list of dimensions (see below)
\param data pointer to data pointer (if *data is NULL, data will be
allocated, otherwise not)
Here, dims is an array of 8 ints, similar to nim->dim[8]. While dims[0]
is unused at this point, the other indices specify which dimensions to
collapse (and at which index), and which not to collapse. If dims[i] is
set to -1, then that entire dimension will be read in, from index 0 to
index (nim->dim[i] - 1). If dims[i] >= 0, then only that index will be
read in (so dims[i] must also be < nim->dim[i]).
Example: given nim->dim[8] = { 4, 64, 64, 21, 80, 1, 1, 1 } (4-D dataset)
if dims[8] = { 0, 5, 4, 17, -1, -1, -1, -1 }
-> read time series for voxel i,j,k = 5,4,17
if dims[8] = { 0, -1, -1, -1, 17, -1, -1, -1 }
-> read single volume at time point 17
Example: given nim->dim[8] = { 6, 64, 64, 21, 80, 4, 3, 1 } (6-D dataset)
if dims[8] = { 0, 5, 4, 17, -1, 2, 1, 0 }
-> read time series for the voxel i,j,k = 5,4,17, and dim 5,6 = 2,1
if dims[8] = { 0, 5, 4, -1, -1, 0, 0, 0 }
-> read time series for slice at i,j = 5,4, and dim 5,6,7 = 0,0,0
(note that dims[7] is not relevant, but must be 0 or -1)
If *data is NULL, then *data will be set as a pointer to new memory,
allocated here for the resulting collapsed image data.
e.g. { int dims[8] = { 0, 5, 4, 17, -1, -1, -1, -1 };
void * data = NULL;
ret_val = nifti_read_collapsed_image(nim, dims, &data);
if( ret_val > 0 ){
process_time_series(data);
if( data != NULL ) free(data);
}
}
NOTE: If *data is not NULL, then it will be assumed that it points to
valid memory, sufficient to hold the results. This is done for
speed and possibly repeated calls to this function.
e.g. { int dims[8] = { 0, -1, -1, -1, -1, -1, -1, -1 };
void * data = NULL;
for( zslice = 0; zslice < nzslices; zslice++ ){
dims[3] = zslice;
ret_val = nifti_read_collapsed_image(nim, dims, &data);
if( ret_val > 0 ) process_slice(zslice, data);
}
if( data != NULL ) free(data);
}
\return
- the total number of bytes read, or < 0 on failure
- the read and byte-swapped data, in 'data'
\sa nifti_image_read, nifti_image_free, nifti_image_read_bricks
nifti_image_load
*//*-------------------------------------------------------------------------*/
int nifti_read_collapsed_image( nifti_image * nim, const int dims [8],
void ** data )
{
znzFile fp;
int pivots[8], prods[8], nprods; /* sizes are bounded by dims[], so 8 */
int c, bytes;
/** - check pointers for sanity */
if( !nim || !dims || !data ){
fprintf(stderr,"** nifti_RCI: bad params %p, %p, %p\n",
(void *)nim, (void *)dims, (void *)data);
return -1;
}
if( g_opts.debug > 2 ){
fprintf(stderr,"-d read_collapsed_image:\n dims =");
for(c = 0; c < 8; c++) fprintf(stderr," %3d", dims[c]);
fprintf(stderr,"\n nim->dims =");
for(c = 0; c < 8; c++) fprintf(stderr," %3d", nim->dim[c]);
fputc('\n', stderr);
}
/** - verify that dim[] makes sense */
if( ! nifti_nim_is_valid(nim, g_opts.debug > 0) ){
fprintf(stderr,"** invalid nim (file is '%s')\n", nim->fname );
return -1;
}
/** - verify that dims[] makes sense for this dataset */
for( c = 1; c <= nim->dim[0]; c++ ){
if( dims[c] >= nim->dim[c] ){
fprintf(stderr,"** nifti_RCI: dims[%d] >= nim->dim[%d] (%d,%d)\n",
c, c, dims[c], nim->dim[c]);
return -1;
}
}
/** - prepare pivot list - pivots are fixed indices */
if( make_pivot_list(nim, dims, pivots, prods, &nprods) < 0 ) return -1;
bytes = rci_alloc_mem(data, prods, nprods, nim->nbyper);
if( bytes < 0 ) return -1;
/** - open the image file for reading at the appropriate offset */
fp = nifti_image_load_prep( nim );
if( ! fp ){ free(*data); *data = NULL; return -1; } /* failure */
/** - call the recursive reading function, passing nim, the pivot info,
location to store memory, and file pointer and position */
c = rci_read_data(nim, pivots,prods,nprods,dims,
(char *)*data, fp, znztell(fp));
znzclose(fp); /* in any case, close the file */
if( c < 0 ){ free(*data); *data = NULL; return -1; } /* failure */
if( g_opts.debug > 1 )
fprintf(stderr,"+d read %d bytes of collapsed image from %s\n",
bytes, nim->fname);
return bytes;
}
/* local function to find strides per dimension. assumes 7D size and
** stride array.
*/
static void
compute_strides(int *strides,const int *size,int nbyper)
{
int i;
strides[0] = nbyper;
for(i = 1; i < 7; i++)
{
strides[i] = size[i-1] * strides[i-1];
}
}
/*---------------------------------------------------------------------------*/
/*! read an arbitrary subregion from a nifti image
This function may be used to read a single arbitary subregion of any
rectangular size from a nifti dataset, such as a small 5x5x5 subregion
around the center of a 3D image.
\param nim given nifti_image struct, corresponding to the data file
\param start_index the index location of the first voxel that will be returned
\param region_size the size of the subregion to be returned
\param data pointer to data pointer (if *data is NULL, data will be
allocated, otherwise not)
Example: given nim->dim[8] = {3, 64, 64, 64, 1, 1, 1, 1 } (3-D dataset)
if start_index[7] = { 29, 29, 29, 0, 0, 0, 0 } and
region_size[7] = { 5, 5, 5, 1, 1, 1, 1 }
-> read 5x5x5 region starting with the first voxel location at (29,29,29)
NOTE: If *data is not NULL, then it will be assumed that it points to
valid memory, sufficient to hold the results. This is done for
speed and possibly repeated calls to this function.
\return
- the total number of bytes read, or < 0 on failure
- the read and byte-swapped data, in 'data'
\sa nifti_image_read, nifti_image_free, nifti_image_read_bricks
nifti_image_load, nifti_read_collapsed_image
*//*-------------------------------------------------------------------------*/
int nifti_read_subregion_image( nifti_image * nim,
int *start_index,
int *region_size,
void ** data )
{
znzFile fp; /* file to read */
int i,j,k,l,m,n; /* indices for dims */
long int bytes = 0; /* total # bytes read */
int total_alloc_size; /* size of buffer allocation */
char *readptr; /* where in *data to read next */
int strides[7]; /* strides between dimensions */
int collapsed_dims[8]; /* for read_collapsed_image */
int *image_size; /* pointer to dimensions in header */
long int initial_offset;
long int offset; /* seek offset for reading current row */
/* probably ignored, but set to ndim for consistency*/
collapsed_dims[0] = nim->ndim;
/* build a dims array for collapsed image read */
for(i = 0; i < nim->ndim; i++)
{
/* if you take the whole extent in this dimension */
if(start_index[i] == 0 &&
region_size[i] == nim->dim[i+1])
{
collapsed_dims[i+1] = -1;
}
/* if you specify a single element in this dimension */
else if(region_size[i] == 1)
{
collapsed_dims[i+1] = start_index[i];
}
else
{
collapsed_dims[i+1] = -2; /* sentinel value */
}
}
/* fill out end of collapsed_dims */
for(i = nim->ndim ; i < 7; i++)
{
collapsed_dims[i+1] = -1;
}
/* check to see whether collapsed read is possible */
for(i = 1; i <= nim->ndim; i++)
{
if(collapsed_dims[i] == -2)
{
break;
}
}
/* if you get through all the dimensions without hitting
** a subrange of size > 1, a collapsed read is possible
*/
if(i > nim->ndim)
{
return nifti_read_collapsed_image(nim, collapsed_dims, data);
}
/* point past first element of dim, which holds nim->ndim */
image_size = &(nim->dim[1]);
/* check region sizes for sanity */
for(i = 0; i < nim->ndim; i++)
{
if(start_index[i] + region_size[i] > image_size[i])
{
if(g_opts.debug > 1)
{
fprintf(stderr,"region doesn't fit within image size\n");
}
return -1;
}
}
/* get the file open */
fp = nifti_image_load_prep( nim );
/* the current offset is just past the nifti header, save
* location so that SEEK_SET can be used below
*/
initial_offset = znztell(fp);
/* get strides*/
compute_strides(strides,image_size,nim->nbyper);
total_alloc_size = nim->nbyper; /* size of pixel */
/* find alloc size */
for(i = 0; i < nim->ndim; i++)
{
total_alloc_size *= region_size[i];
}
/* allocate buffer, if necessary */
if(*data == 0)
{
*data = (void *)malloc(total_alloc_size);
}
if(*data == 0)
{
if(g_opts.debug > 1)
{
fprintf(stderr,"allocation of %d bytes failed\n",total_alloc_size);
return -1;
}
}
/* point to start of data buffer as char * */
readptr = *((char **)data);
{
/* can't assume that start_index and region_size have any more than
** nim->ndim elements so make local copies, filled out to seven elements
*/
int si[7], rs[7];
for(i = 0; i < nim->ndim; i++)
{
si[i] = start_index[i];
rs[i] = region_size[i];
}
for(i = nim->ndim; i < 7; i++)
{
si[i] = 0;
rs[i] = 1;
}
/* loop through subregion and read a row at a time */
for(i = si[6]; i < (si[6] + rs[6]); i++)
{
for(j = si[5]; j < (si[5] + rs[5]); j++)
{
for(k = si[4]; k < (si[4] + rs[4]); k++)
{
for(l = si[3]; l < (si[3] + rs[3]); l++)
{
for(m = si[2]; m < (si[2] + rs[2]); m++)
{
for(n = si[1]; n < (si[1] + rs[1]); n++)
{
int nread,read_amount;
offset = initial_offset +
(i * strides[6]) +
(j * strides[5]) +
(k * strides[4]) +
(l * strides[3]) +
(m * strides[2]) +
(n * strides[1]) +
(si[0] * strides[0]);
znzseek(fp, offset, SEEK_SET); /* seek to current row */
read_amount = rs[0] * nim->nbyper; /* read a row of the subregion*/
nread = (int)nifti_read_buffer(fp, readptr, read_amount, nim);
if(nread != read_amount)
{
if(g_opts.debug > 1)
{
fprintf(stderr,"read of %d bytes failed\n",read_amount);
return -1;
}
}
bytes += nread;
readptr += read_amount;
}
}
}
}
}
}
}
return bytes;
}
/* read the data from the file pointed to by fp
- this a recursive function, so start with the base case
- data is now (char *) for easy incrementing
return 0 on success, < 0 on failure
*/
static int rci_read_data(nifti_image * nim, int * pivots, int * prods,
int nprods, const int dims[], char * data, znzFile fp, size_t base_offset)
{
size_t sublen, offset, read_size;
int c;
/* bad check first - base_offset may not have been checked */
if( nprods <= 0 ){
fprintf(stderr,"** rci_read_data, bad prods, %d\n", nprods);
return -1;
}
/* base case: actually read the data */
if( nprods == 1 ){
size_t nread, bytes;
/* make sure things look good here */
if( *pivots != 0 ){
fprintf(stderr,"** rciRD: final pivot == %d!\n", *pivots);
return -1;
}
/* so just seek and read (prods[0] * nbyper) bytes from the file */
znzseek(fp, (long)base_offset, SEEK_SET);
bytes = (size_t)prods[0] * nim->nbyper;
nread = nifti_read_buffer(fp, data, bytes, nim);
if( nread != bytes ){
fprintf(stderr,"** rciRD: read only %u of %u bytes from '%s'\n",
(unsigned)nread, (unsigned)bytes, nim->fname);
return -1;
} else if( g_opts.debug > 3 )
fprintf(stderr,"+d successful read of %u bytes at offset %u\n",
(unsigned)bytes, (unsigned)base_offset);
return 0; /* done with base case - return success */
}
/* not the base case, so do a set of reduced reads */
/* compute size of sub-brick: all dimensions below pivot */
for( c = 1, sublen = 1; c < *pivots; c++ ) sublen *= nim->dim[c];
/* compute number of values to read, i.e. remaining prods */
for( c = 1, read_size = 1; c < nprods; c++ ) read_size *= prods[c];
read_size *= nim->nbyper; /* and multiply by bytes per voxel */
/* now repeatedly compute offsets, and recursively read */
for( c = 0; c < prods[0]; c++ ){
/* offset is (c * sub-block size (including pivot dim)) */
/* + (dims[] index into pivot sub-block) */
/* the unneeded multiplication is to make this more clear */
offset = (size_t)c * sublen * nim->dim[*pivots] +
(size_t)sublen * dims[*pivots];
offset *= nim->nbyper;
if( g_opts.debug > 3 )
fprintf(stderr,"-d reading %u bytes, foff %u + %u, doff %u\n",
(unsigned)read_size, (unsigned)base_offset, (unsigned)offset,
(unsigned)(c*read_size));
/* now read the next level down, adding this offset */
if( rci_read_data(nim, pivots+1, prods+1, nprods-1, dims,
data + c * read_size, fp, base_offset + offset) < 0 )
return -1;
}
return 0;
}
/* allocate memory for all collapsed image data
If *data is already set, do not allocate, but still calculate
size for debug report.
return total size on success, and < 0 on failure
*/
static int rci_alloc_mem(void ** data, int prods[8], int nprods, int nbyper )
{
int size, index;
if( nbyper < 0 || nprods < 1 || nprods > 8 ){
fprintf(stderr,"** rci_am: bad params, %d, %d\n", nbyper, nprods);
return -1;
}
for( index = 0, size = 1; index < nprods; index++ )
size *= prods[index];
size *= nbyper;
if( ! *data ){ /* then allocate what is needed */
if( g_opts.debug > 1 )
fprintf(stderr,"+d alloc %d (= %d x %d) bytes for collapsed image\n",
size, size/nbyper, nbyper);
*data = malloc(size); /* actually allocate the memory */
if( ! *data ){
fprintf(stderr,"** rci_am: failed to alloc %d bytes for data\n", size);
return -1;
}
} else if( g_opts.debug > 1 )
fprintf(stderr,"-d rci_am: *data already set, need %d (%d x %d) bytes\n",
size, size/nbyper, nbyper);
return size;
}
/* prepare a pivot list for reading
The pivot points are the indices into dims where the calling function
wants to collapse a dimension. The last pivot should always be zero
(note that we have space for that in the lists).
*/
static int make_pivot_list(nifti_image * nim, const int dims[], int pivots[],
int prods[], int * nprods )
{
int len, index;
len = 0;
index = nim->dim[0];
while( index > 0 ){
prods[len] = 1;
while( index > 0 && (nim->dim[index] == 1 || dims[index] == -1) ){
prods[len] *= nim->dim[index];
index--;
}
pivots[len] = index;
len++;
index--; /* fine, let it drop out at -1 */
}
/* make sure to include 0 as a pivot (instead of just 1, if it is) */
if( pivots[len-1] != 0 ){
pivots[len] = 0;
prods[len] = 1;
len++;
}
*nprods = len;
if( g_opts.debug > 2 ){
fprintf(stderr,"+d pivot list created, pivots :");
for(index = 0; index < len; index++) fprintf(stderr," %d", pivots[index]);
fprintf(stderr,", prods :");
for(index = 0; index < len; index++) fprintf(stderr," %d", prods[index]);
fputc('\n',stderr);
}
return 0;
}
#undef ISEND
#define ISEND(c) ( (c)==']' || (c)=='}' || (c)=='\0' )
/*---------------------------------------------------------------------*/
/*! Get an integer list in the range 0..(nvals-1), from the
character string str. If we call the output pointer fred,
then fred[0] = number of integers in the list (> 0), and
fred[i] = i-th integer in the list for i=1..fred[0].
If on return, fred == NULL or fred[0] == 0, then something is
wrong, and the caller must deal with that.
Syntax of input string:
- initial '{' or '[' is skipped, if present
- ends when '}' or ']' or end of string is found
- contains entries separated by commas
- entries have one of these forms:
- a single number
- a dollar sign '$', which means nvals-1
- a sequence of consecutive numbers in the form "a..b" or
"a-b", where "a" and "b" are single numbers (or '$')
- a sequence of evenly spaced numbers in the form
"a..b(c)" or "a-b(c)", where "c" encodes the step
- Example: "[2,7..4,3..9(2)]" decodes to the list
2 7 6 5 4 3 5 7 9
- entries should be in the range 0..nvals-1
(borrowed, with permission, from thd_intlist.c)
*//*-------------------------------------------------------------------*/
int * nifti_get_intlist( int nvals , const char * str )
{
int *subv = NULL ;
int ii , ipos , nout , slen ;
int ibot,itop,istep , nused ;
char *cpt ;
/* Meaningless input? */
if( nvals < 1 ) return NULL ;
/* No selection list? */
if( str == NULL || str[0] == '\0' ) return NULL ;
/* skip initial '[' or '{' */
subv = (int *)malloc( sizeof(int) * 2 ) ;
if( !subv ) {
fprintf(stderr,"** nifti_get_intlist: failed alloc of 2 ints\n");
return NULL;
}
subv[0] = nout = 0 ;
ipos = 0 ;
if( str[ipos] == '[' || str[ipos] == '{' ) ipos++ ;
if( g_opts.debug > 1 )
fprintf(stderr,"-d making int_list (vals = %d) from '%s'\n", nvals, str);
/**- for each sub-selector until end of input... */
slen = (int)strlen(str) ;
while( ipos < slen && !ISEND(str[ipos]) ){
while( isspace((int) str[ipos]) ) ipos++ ; /* skip blanks */
if( ISEND(str[ipos]) ) break ; /* done */
/**- get starting value */
if( str[ipos] == '$' ){ /* special case */
ibot = nvals-1 ; ipos++ ;
} else { /* decode an integer */
ibot = strtol( str+ipos , &cpt , 10 ) ;
if( ibot < 0 ){
fprintf(stderr,"** ERROR: list index %d is out of range 0..%d\n",
ibot,nvals-1) ;
free(subv) ; return NULL ;
}
if( ibot >= nvals ){
fprintf(stderr,"** ERROR: list index %d is out of range 0..%d\n",
ibot,nvals-1) ;
free(subv) ; return NULL ;
}
nused = (cpt-(str+ipos)) ;
if( ibot == 0 && nused == 0 ){
fprintf(stderr,"** ERROR: list syntax error '%s'\n",str+ipos) ;
free(subv) ; return NULL ;
}
ipos += nused ;
}
while( isspace((int) str[ipos]) ) ipos++ ; /* skip blanks */
/**- if that's it for this sub-selector, add one value to list */
if( str[ipos] == ',' || ISEND(str[ipos]) ){
nout++ ;
subv = (int *)realloc( (char *)subv , sizeof(int) * (nout+1) ) ;
if( !subv ) {
fprintf(stderr,"** nifti_get_intlist: failed realloc of %d ints\n",
nout+1);
return NULL;
}
subv[0] = nout ;
subv[nout] = ibot ;
if( ISEND(str[ipos]) ) break ; /* done */
ipos++ ; continue ; /* re-start loop at next sub-selector */
}
/**- otherwise, must have '..' or '-' as next inputs */
if( str[ipos] == '-' ){
ipos++ ;
} else if( str[ipos] == '.' && str[ipos+1] == '.' ){
ipos++ ; ipos++ ;
} else {
fprintf(stderr,"** ERROR: index list syntax is bad: '%s'\n",
str+ipos) ;
free(subv) ; return NULL ;
}
/**- get ending value for loop now */
if( str[ipos] == '$' ){ /* special case */
itop = nvals-1 ; ipos++ ;
} else { /* decode an integer */
itop = strtol( str+ipos , &cpt , 10 ) ;
if( itop < 0 ){
fprintf(stderr,"** ERROR: index %d is out of range 0..%d\n",
itop,nvals-1) ;
free(subv) ; return NULL ;
}
if( itop >= nvals ){
fprintf(stderr,"** ERROR: index %d is out of range 0..%d\n",
itop,nvals-1) ;
free(subv) ; return NULL ;
}
nused = (cpt-(str+ipos)) ;
if( itop == 0 && nused == 0 ){
fprintf(stderr,"** ERROR: index list syntax error '%s'\n",str+ipos) ;
free(subv) ; return NULL ;
}
ipos += nused ;
}
/**- set default loop step */
istep = (ibot <= itop) ? 1 : -1 ;
while( isspace((int) str[ipos]) ) ipos++ ; /* skip blanks */
/**- check if we have a non-default loop step */
if( str[ipos] == '(' ){ /* decode an integer */
ipos++ ;
istep = strtol( str+ipos , &cpt , 10 ) ;
if( istep == 0 ){
fprintf(stderr,"** ERROR: index loop step is 0!\n") ;
free(subv) ; return NULL ;
}
nused = (cpt-(str+ipos)) ;
ipos += nused ;
if( str[ipos] == ')' ) ipos++ ;
if( (ibot-itop)*istep > 0 ){
fprintf(stderr,"** WARNING: index list '%d..%d(%d)' means nothing\n",
ibot,itop,istep ) ;
}
}
/**- add values to output */
for( ii=ibot ; (ii-itop)*istep <= 0 ; ii += istep ){
nout++ ;
subv = (int *)realloc( (char *)subv , sizeof(int) * (nout+1) ) ;
if( !subv ) {
fprintf(stderr,"** nifti_get_intlist: failed realloc of %d ints\n",
nout+1);
return NULL;
}
subv[0] = nout ;
subv[nout] = ii ;
}
/**- check if we have a comma to skip over */
while( isspace((int) str[ipos]) ) ipos++ ; /* skip blanks */
if( str[ipos] == ',' ) ipos++ ; /* skip commas */
} /* end of loop through selector string */
if( g_opts.debug > 1 ) {
fprintf(stderr,"+d int_list (vals = %d): ", subv[0]);
for( ii = 1; ii <= subv[0]; ii++ ) fprintf(stderr,"%d ", subv[ii]);
fputc('\n',stderr);
}
if( subv[0] == 0 ){ free(subv); subv = NULL; }
return subv ;
}
/*---------------------------------------------------------------------*/
/*! Given a NIFTI_TYPE string, such as "NIFTI_TYPE_INT16", return the
* corresponding integral type code. The type code is the macro
* value defined in nifti1.h.
*//*-------------------------------------------------------------------*/
int nifti_datatype_from_string( const char * name )
{
int tablen = sizeof(nifti_type_list)/sizeof(nifti_type_ele);
int c;
if( !name ) return DT_UNKNOWN;
for( c = tablen-1; c > 0; c-- )
if( !strcmp(name, nifti_type_list[c].name) )
break;
return nifti_type_list[c].type;
}
/*---------------------------------------------------------------------*/
/*! Given a NIFTI_TYPE value, such as NIFTI_TYPE_INT16, return the
* corresponding macro label as a string. The dtype code is the
* macro value defined in nifti1.h.
*//*-------------------------------------------------------------------*/
char * nifti_datatype_to_string( int dtype )
{
int tablen = sizeof(nifti_type_list)/sizeof(nifti_type_ele);
int c;
for( c = tablen-1; c > 0; c-- )
if( nifti_type_list[c].type == dtype )
break;
return nifti_type_list[c].name;
}
/*---------------------------------------------------------------------*/
/*! Determine whether dtype is a valid NIFTI_TYPE.
*
* DT_UNKNOWN is considered invalid
*
* The only difference 'for_nifti' makes is that DT_BINARY
* should be invalid for a NIfTI dataset.
*//*-------------------------------------------------------------------*/
int nifti_datatype_is_valid( int dtype, int for_nifti )
{
int tablen = sizeof(nifti_type_list)/sizeof(nifti_type_ele);
int c;
/* special case */
if( for_nifti && dtype == DT_BINARY ) return 0;
for( c = tablen-1; c > 0; c-- )
if( nifti_type_list[c].type == dtype )
return 1;
return 0;
}
/*---------------------------------------------------------------------*/
/*! Only as a test, verify that the new nifti_type_list table matches
* the the usage of nifti_datatype_sizes (which could be changed to
* use the table, if there were interest).
*
* return the number of errors (so 0 is success, as usual)
*//*-------------------------------------------------------------------*/
int nifti_test_datatype_sizes(int verb)
{
int tablen = sizeof(nifti_type_list)/sizeof(nifti_type_ele);
int nbyper, ssize;
int c, errs = 0;
for( c = 0; c < tablen; c++ )
{
nbyper = ssize = -1;
nifti_datatype_sizes(nifti_type_list[c].type, &nbyper, &ssize);
if( nbyper < 0 || ssize < 0 ||
nbyper != nifti_type_list[c].nbyper ||
ssize != nifti_type_list[c].swapsize )
{
if( verb || g_opts.debug > 2 )
fprintf(stderr, "** type mismatch: %s, %d, %d, %d : %d, %d\n",
nifti_type_list[c].name, nifti_type_list[c].type,
nifti_type_list[c].nbyper, nifti_type_list[c].swapsize,
nbyper, ssize);
errs++;
}
}
if( errs )
fprintf(stderr,"** nifti_test_datatype_sizes: found %d errors\n",errs);
else if( verb || g_opts.debug > 1 )
fprintf(stderr,"-- nifti_test_datatype_sizes: all OK\n");
return errs;
}
/*---------------------------------------------------------------------*/
/*! Display the nifti_type_list table.
*
* if which == 1 : display DT_*
* if which == 2 : display NIFTI_TYPE*
* else : display all
*//*-------------------------------------------------------------------*/
int nifti_disp_type_list( int which )
{
char * style;
int tablen = sizeof(nifti_type_list)/sizeof(nifti_type_ele);
int lwhich, c;
if ( which == 1 ){ lwhich = 1; style = "DT_"; }
else if( which == 2 ){ lwhich = 2; style = "NIFTI_TYPE_"; }
else { lwhich = 3; style = "ALL"; }
printf("nifti_type_list entries (%s) :\n"
" name type nbyper swapsize\n"
" --------------------- ---- ------ --------\n", style);
for( c = 0; c < tablen; c++ )
if( (lwhich & 1 && nifti_type_list[c].name[0] == 'D') ||
(lwhich & 2 && nifti_type_list[c].name[0] == 'N') )
printf(" %-22s %5d %3d %5d\n",
nifti_type_list[c].name,
nifti_type_list[c].type,
nifti_type_list[c].nbyper,
nifti_type_list[c].swapsize);
return 0;
}