* { 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 = 11. 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 values for blist are in [0..nt*nu*nv*nw-1],
* or written [ 0 .. (dim[4]*dim[5]*dim[6]*dim[7] - 1) ].
*//*----------------------------------------------------------------------*/
nifti_image *nifti_image_read_bricks( char *hname, int nbricks, 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->nvox = nim->nx * nim->ny * nim->nz
* nim->nt * nim->nu * nim->nv * nim->nw;
nim->dim[4] = nbricks;
nim->dim[5] = nim->dim[6] = nim->dim[7] = 1;
/* 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[].
*//*--------------------------------------------------------------------*/
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);
}
/* set nx, ..., one by one (abusing promotion, just a little) */
if(nim->dim[1] <= 1){ nim->pixdim[1] = nim->dx = nim->nx = nim->dim[1] = 1; }
if(nim->dim[2] <= 1){ nim->pixdim[2] = nim->dy = nim->ny = nim->dim[2] = 1; }
if(nim->dim[3] <= 1){ nim->pixdim[3] = nim->dz = nim->nz = nim->dim[3] = 1; }
if(nim->dim[4] <= 1){ nim->pixdim[4] = nim->dt = nim->nt = nim->dim[4] = 1; }
if(nim->dim[5] <= 1){ nim->pixdim[5] = nim->du = nim->nu = nim->dim[5] = 1; }
if(nim->dim[6] <= 1){ nim->pixdim[6] = nim->dv = nim->nv = nim->dim[6] = 1; }
if(nim->dim[7] <= 1){ nim->pixdim[7] = nim->dw = nim->nw = nim->dim[7] = 1; }
nim->nvox = nim->nx * nim->ny * nim->nz
* nim->nt * nim->nu * nim->nv * nim->nw;
/* compute ndim */
for( ndim = 7; (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
* \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 the dataset normally.
*
* \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 NBL->bricks[i]).
*//*--------------------------------------------------------------------*/
int nifti_image_load_bricks( nifti_image * nim , int nbricks, 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->nbricks = NBL->bsize = 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 )
{
int c, rv;
int oposn, fposn; /* orig and current file positions */
int prev, isrc, idest; /* previous and current sub-brick, and new index */
oposn = znztell(fp); /* store current file position */
fposn = oposn;
if( fposn < 0 ){
fprintf(stderr,"** load bricks: ztell failed??\n");
return -1;
}
/* 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 bricks from file %s\n",
NBL->nbricks, 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, 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);
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 nbl->nbricks = nim->nt * nim->nu * nim->nv * nim->nw;
nbl->bsize = 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 %d bytes for brick %d\n",
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->nbricks = nbl->bsize = 0;
return -1;
}
}
if( g_opts.debug > 2 )
fprintf(stderr,"+d NANM: alloc'd %d bricks of %d bytes for NBL\n",
nbl->nbricks, 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, 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.
*
* \return 1 if valid, 0 if not
*//*--------------------------------------------------------------------*/
int valid_nifti_brick_list(nifti_image * nim , int nbricks, 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;
}
/* nsubs sub-brick is nt*nu*nv*nw */
nsubs = nim->dim[4] * nim->dim[5] * nim->dim[6] * nim->dim[7];
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]);
int_force_positive(nim->dim+4, 4);
nsubs = nim->dim[4] * nim->dim[5] * nim->dim[6] * nim->dim[7];
}
for( c = 0; c < nbricks; c++ )
if( (blist[c] < 0) || (blist[c] >= nsubs) ){
if( disp_error || g_opts.debug > 1 )
fprintf(stderr,
"-d ** bad sub-brick chooser %d (#%d), valid range is [0,%d]\n",
blist[c], c, nsubs-1);
return 0;
}
return 1; /* all is well */
}
/* set any non-positive values to 1 */
static int int_force_positive( int * list, int nel )
{
int c;
if( !list || nel < 0 ){
if( g_opts.debug > 0 )
fprintf(stderr,"** int_force_positive: bad params (%p,%d)\n",
(void *)list,nel);
return -1;
}
for( c = 0; c < nel; c++ )
if( list[c] <= 0 ) list[c] = 1;
return 0;
}
/* end of new nifti_image_read_bricks() functionality */
/*----------------------------------------------------------------------*/
/*! display the orientation from the quaternian fields
*
* \return -1 if results cannot be determined, 0 if okay
*//*--------------------------------------------------------------------*/
int nifti_disp_matrix_orient( 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= (char *)malloc( strlen(str)+1 );
if (dup) strcpy(dup,str);
return dup;
}
/*---------------------------------------------------------------------------*/
/*! Return a pointer to a string holding the name of a NIFTI datatype.
Don't free() or modify this string! It points to static storage.
*//*-------------------------------------------------------------------------*/
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" ;
}
return "**ILLEGAL**" ;
}
/*----------------------------------------------------------------------*/
/*! Determine if the datatype code dt is an integer type (1=YES, 0=NO).
*//*--------------------------------------------------------------------*/
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 ;
}
return 0 ;
}
/*---------------------------------------------------------------------------*/
/*! Return a pointer to a string holding the name of a NIFTI units type.
Don't free() or modify this string! It points to static storage.
*//*-------------------------------------------------------------------------*/
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.
Don't free() or modify this string! It points to static storage.
*//*-------------------------------------------------------------------------*/
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.
Don't free() or modify this string! It points to static storage.
*//*-------------------------------------------------------------------------*/
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.
Don't free() or modify this string! It points to static storage.
*//*-------------------------------------------------------------------------*/
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" ;
}
return "Unknown" ;
}
/*---------------------------------------------------------------------------*/
/*! Return a pointer to a string holding the name of a NIFTI orientation.
Don't free() or modify this string! It points to static storage.
*//*-------------------------------------------------------------------------*/
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.
The swapsize is set to 0 if this datatype doesn't ever need swapping.
*//*------------------------------------------------------------------------*/
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_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.*//*-------------------------------------------------------------------------*/ 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. See comments in nifti1.h for details. - 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? *//*-------------------------------------------------------------------------*/ 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. *//*-------------------------------------------------------------------------*/ 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]). *//*-------------------------------------------------------------------------*/ 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)
*//*-------------------------------------------------------------------------*/
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]
-----------------------------------------------------------------------------*/
typedef struct { unsigned char a,b ; } twobytes ;
/*----------------------------------------------------------------------*/
/*! swap each byte pair from the given list of n pairs
*//*--------------------------------------------------------------------*/
void nifti_swap_2bytes( int n , void *ar ) /* 2 bytes at a time */
{
register int ii ;
register twobytes *tb = (twobytes *)ar ;
register unsigned char tt ;
for( ii=0 ; ii < n ; ii++ ){
tt = tb[ii].a ; tb[ii].a = tb[ii].b ; tb[ii].b = tt ;
}
return ;
}
/*---------------------------------------------------------------------------*/
typedef struct { unsigned char a,b,c,d ; } fourbytes ;
/*----------------------------------------------------------------------*/
/*! swap 4 bytes at a time from the given list of n sets of 4 bytes
*//*--------------------------------------------------------------------*/
void nifti_swap_4bytes( int n , void *ar ) /* 4 bytes at a time */
{
register int ii ;
register fourbytes *tb = (fourbytes *)ar ;
register unsigned char tt ;
for( ii=0 ; ii < n ; ii++ ){
tt = tb[ii].a ; tb[ii].a = tb[ii].d ; tb[ii].d = tt ;
tt = tb[ii].b ; tb[ii].b = tb[ii].c ; tb[ii].c = tt ;
}
return ;
}
/*---------------------------------------------------------------------------*/
typedef struct { unsigned char a,b,c,d , D,C,B,A ; } eightbytes ;
/*----------------------------------------------------------------------*/
/*! swap 8 bytes at a time from the given list of n sets of 8 bytes
*//*--------------------------------------------------------------------*/
void nifti_swap_8bytes( int n , void *ar ) /* 8 bytes at a time */
{
register int ii ;
register eightbytes *tb = (eightbytes *)ar ;
register unsigned char tt ;
for( ii=0 ; ii < n ; ii++ ){
tt = tb[ii].a ; tb[ii].a = tb[ii].A ; tb[ii].A = tt ;
tt = tb[ii].b ; tb[ii].b = tb[ii].B ; tb[ii].B = tt ;
tt = tb[ii].c ; tb[ii].c = tb[ii].C ; tb[ii].C = tt ;
tt = tb[ii].d ; tb[ii].d = tb[ii].D ; tb[ii].D = tt ;
}
return ;
}
/*---------------------------------------------------------------------------*/
typedef struct { unsigned char a,b,c,d,e,f,g,h ,
H,G,F,E,D,C,B,A ; } sixteenbytes ;
/*----------------------------------------------------------------------*/
/*! swap 16 bytes at a time from the given list of n sets of 16 bytes
*//*--------------------------------------------------------------------*/
void nifti_swap_16bytes( int n , void *ar ) /* 16 bytes at a time */
{
register int ii ;
register sixteenbytes *tb = (sixteenbytes *)ar ;
register unsigned char tt ;
for( ii=0 ; ii < n ; ii++ ){
tt = tb[ii].a ; tb[ii].a = tb[ii].A ; tb[ii].A = tt ;
tt = tb[ii].b ; tb[ii].b = tb[ii].B ; tb[ii].B = tt ;
tt = tb[ii].c ; tb[ii].c = tb[ii].C ; tb[ii].C = tt ;
tt = tb[ii].d ; tb[ii].d = tb[ii].D ; tb[ii].D = tt ;
tt = tb[ii].e ; tb[ii].e = tb[ii].E ; tb[ii].E = tt ;
tt = tb[ii].f ; tb[ii].f = tb[ii].F ; tb[ii].F = tt ;
tt = tb[ii].g ; tb[ii].g = tb[ii].G ; tb[ii].G = tt ;
tt = tb[ii].h ; tb[ii].h = tb[ii].H ; tb[ii].H = tt ;
}
return ;
}
/*---------------------------------------------------------------------------*/
/*----------------------------------------------------------------------*/
/*! based on siz, call the appropriate nifti_swap_Nbytes() function
*//*--------------------------------------------------------------------*/
void nifti_swap_Nbytes( int n , int siz , void *ar ) /* subsuming case */
{
switch( siz ){
case 2: nifti_swap_2bytes ( n , ar ) ; break ;
case 4: nifti_swap_4bytes ( n , ar ) ; break ;
case 8: nifti_swap_8bytes ( n , ar ) ; break ;
case 16: nifti_swap_16bytes( n , ar ) ; break ;
}
return ;
}
/*-------------------------------------------------------------------------*/
/*! Byte swap NIFTI-1 file header in various places and ways.
If is_nifti is nonzero, will also swap the NIFTI-specific
components of the header; otherwise, only the components
common to NIFTI and ANALYZE will be swapped.
*//*---------------------------------------------------------------------- */
void swap_nifti_header( struct nifti_1_header *h , int is_nifti )
{
#if 0 /* ANALYZE fields not used by this software */
swap_4(h->sizeof_hdr) ;
swap_4(h->extents) ;
swap_2(h->session_error) ;
swap_4(h->compressed) ;
swap_4(h->glmax) ; swap_4(h->glmin) ;
#endif
/* this stuff is always present, for ANALYZE and NIFTI */
swap_4(h->sizeof_hdr) ;
nifti_swap_2bytes( 8 , h->dim ) ;
nifti_swap_4bytes( 8 , h->pixdim ) ;
swap_2(h->datatype) ;
swap_2(h->bitpix) ;
swap_4(h->vox_offset); swap_4(h->cal_max); swap_4(h->cal_min);
/* this stuff is NIFTI specific */
if( is_nifti ){
swap_4(h->intent_p1); swap_4(h->intent_p2); swap_4(h->intent_p3);
swap_2(h->intent_code);
swap_2(h->slice_start); swap_2(h->slice_end);
swap_4(h->scl_slope); swap_4(h->scl_inter);
swap_4(h->slice_duration); swap_4(h->toffset);
swap_2(h->qform_code); swap_2(h->sform_code);
swap_4(h->quatern_b); swap_4(h->quatern_c); swap_4(h->quatern_d);
swap_4(h->qoffset_x); swap_4(h->qoffset_y); swap_4(h->qoffset_z);
nifti_swap_4bytes(4,h->srow_x);
nifti_swap_4bytes(4,h->srow_y);
nifti_swap_4bytes(4,h->srow_z);
}
return ;
}
#define USE_STAT
#ifdef USE_STAT
/*---------------------------------------------------------------------------*/
/* Return the file length (0 if file not found or has no contents).
This is a Unix-specific function, since it uses stat().
-----------------------------------------------------------------------------*/
#include
nifti_1_header my_header;
my_header = nifti_convert_nim2nhdr(my_nim_pointer);
*//*--------------------------------------------------------------------*/
struct nifti_1_header nifti_convert_nim2nhdr(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
*/
int nifti_copy_extensions(nifti_image * nim_dest, 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 = (char *)calloc(size,sizeof(char)); /* calloc, 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);
nim_dest->num_ext++;
}
return 0;
}
/*----------------------------------------------------------------------*/
/*! compute the total size of all extensions
*//*--------------------------------------------------------------------*/
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
*//*--------------------------------------------------------------------*/
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
*//*--------------------------------------------------------------------*/
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.
*
* 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
*//*---------------------------------------------------------------------*/
znzFile nifti_image_write_hdr_img2( nifti_image *nim , int write_opts ,
const char * opts, znzFile imgfile, 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") ;
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,(char*)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'\n",nim->fname);
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, nifti_brick_list * NBL,
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); }
return fp; /* returned but may be closed */
}
/*--------------------------------------------------------------------------*/
/*! Write a nifti_image to disk.
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()
*//*------------------------------------------------------------------------*/
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 */
void nifti_image_write_bricks( nifti_image *nim, 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(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);
/* 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 = 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.
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);
}
}
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.
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
*//*-------------------------------------------------------------------------*/
int nifti_read_collapsed_image( nifti_image * nim, 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;
}
/* 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, int dims[], char * data, znzFile fp, int base_offset)
{
int c, sublen, offset, read_size;
/* bad check first - base_offset may not have been checked */
if( base_offset < 0 || nprods <= 0 ){
fprintf(stderr,"** rci_read_data, bad params, %d,%d\n",
nprods, base_offset);
return -1;
}
/* base case: actually read the data */
if( nprods == 1 ){
int 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, base_offset, SEEK_SET);
bytes = prods[0] * nim->nbyper;
nread = nifti_read_buffer(fp, data, bytes, nim);
if( nread != bytes ){
fprintf(stderr,"** rciRD: read only %d of %d bytes from '%s'\n",
nread, bytes, nim->fname);
return -1;
} else if( g_opts.debug > 3 )
fprintf(stderr,"+d successful read of %d bytes at offset %d\n",
bytes, 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 = c * sublen * nim->dim[*pivots] + sublen * dims[*pivots];
offset *= nim->nbyper;
if( g_opts.debug > 3 )
fprintf(stderr,"-d reading %d bytes, foff %d + %d, doff %d\n",
read_size, base_offset, offset, 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, 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 , 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 ) ;
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 = strlen(str) ;
while( ipos < slen && !ISEND(str[ipos]) ){
while( isspace(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(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) ) ;
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(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) ) ;
subv[0] = nout ;
subv[nout] = ii ;
}
/** check if we have a comma to skip over */
while( isspace(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 ;
}