#include "mrilib.h"
#include "netcdf.h"
/*******************************************************************/
/******* 29 Oct 2001: Read a 3D MINC file as an AFNI dataset *******/
/*******************************************************************/
static int first_err = 1 ;
static char *fname_err ;
#define EPR(x,s1,s2) \
do{ if( first_err ){ fprintf(stderr,"\n"); first_err=0; } \
fprintf(stderr,"** NetCDF error [%s,%s]: %s\n",s1,s2,nc_strerror(x)); \
} while(0)
/*-----------------------------------------------------------------
Read the attributes relevant to a MINC dimension variable
-------------------------------------------------------------------*/
typedef struct {
int dimid , varid , good , afni_orient ;
size_t len ;
float start , step , xcos,ycos,zcos ;
char spacetype[32] ;
} mincdim ;
static mincdim read_mincdim( int ncid , char *dname )
{
mincdim ddd ;
int code ;
float ccc[3] ;
size_t lll ;
static char *ename=NULL ;
ddd.good = 0 ; /* flag for bad result */
ddd.step = 1.0 ; /* defaults */
ddd.xcos = ddd.ycos = ddd.zcos = 0.0 ;
ddd.spacetype[0] = '\0' ;
lll = strlen(fname_err) + strlen(dname) + 4 ;
ename = AFREALL(ename, char, lll) ;
sprintf(ename,"%s:%s",fname_err,dname) ;
/* get ID of this dimension name */
code = nc_inq_dimid( ncid , dname , &ddd.dimid ) ;
if( code != NC_NOERR ){ EPR(code,ename,"dimid"); return ddd; }
/* get length of this dimension */
code = nc_inq_dimlen( ncid , ddd.dimid , &ddd.len ) ;
if( code != NC_NOERR ){ EPR(code,ename,"dimlen"); return ddd; }
/* get ID of corresponding variable */
code = nc_inq_varid( ncid , dname , &ddd.varid ) ;
if( code != NC_NOERR ){ /* this is bad: try to fake it */
if( first_err ){ fprintf(stderr,"\n"); first_err=0; }
fprintf(stderr,"** MINC warning: %s variable missing\n",ename);
ddd.start = -0.5*ddd.step*ddd.len ;
ddd.good = 1 ; return ddd ;
}
/* get step attribute of this variable */
code = nc_get_att_float( ncid , ddd.varid , "step" , &ddd.step ) ;
if( code != NC_NOERR ){
if( first_err ){ fprintf(stderr,"\n"); first_err=0; }
fprintf(stderr,"** MINC warning: %s:step missing\n",ename);
} else if( ddd.step == 0.0 ){
if( first_err ){ fprintf(stderr,"\n"); first_err=0; }
fprintf(stderr,"** MINC warning: %s:step=0\n",ename);
ddd.step = 1.0 ;
}
/* get start attribute of this variable */
code = nc_get_att_float( ncid , ddd.varid , "start" , &ddd.start ) ;
if( code != NC_NOERR ){
if( first_err ){ fprintf(stderr,"\n"); first_err=0; }
fprintf(stderr,"** MINC warning: %s:start missing\n",ename) ;
ddd.start = -0.5*ddd.step*ddd.len ;
}
/* get direction_cosines attribute of this variable [not used yet] */
code = nc_get_att_float( ncid,ddd.varid , "direction_cosines" , ccc ) ;
if( code == NC_NOERR ){
ddd.xcos = ccc[0] ; ddd.ycos = ccc[1] ; ddd.zcos = ccc[2] ;
}
/* get spacetype attribute of this variable [Talairach or not?] */
lll = 0 ;
code = nc_inq_attlen( ncid , ddd.varid , "spacetype" , &lll ) ;
if( code == NC_NOERR && lll > 0 && lll < 32 ){
code = nc_get_att_text( ncid,ddd.varid , "spacetype" , ddd.spacetype ) ;
if( code == NC_NOERR ){
ddd.spacetype[lll] = '\0' ; /* make sure ends in NUL */
} else {
ddd.spacetype[0] = '\0' ; /* make sure is empty */
}
}
ddd.good = 1 ; return ddd ;
}
/*-----------------------------------------------------------------
Open a MINC file as an AFNI dataset
-------------------------------------------------------------------*/
THD_3dim_dataset * THD_open_minc( char *pathname )
{
THD_3dim_dataset *dset=NULL ;
int ncid , code ;
mincdim xspace , yspace , zspace , *xyz[3] ;
int im_varid , im_rank , im_dimid[4] ;
nc_type im_type ;
char prefix[1024] , *ppp , tname[12] ;
THD_ivec3 nxyz , orixyz ;
THD_fvec3 dxyz , orgxyz ;
int datum , nvals=1 , ibr , iview ;
size_t len ;
float dtime=1.0 ;
ENTRY("THD_open_minc") ;
/*-- open input file --*/
if( pathname == NULL || pathname[0] == '\0' ) RETURN(NULL);
first_err = 1 ; /* put a newline before 1st error message */
fname_err = pathname ; /* filename for errors in read_mincdim() */
code = nc_open( pathname , NC_NOWRITE , &ncid ) ;
if( code != NC_NOERR ){ EPR(code,pathname,"open"); RETURN(NULL); }
/*---------- get info about spatial axes ---------*/
/*[[ MINC x and y are reversed relative to AFNI ]]*/
/*[[ MINC x=L-R y=P-A; AFNI x=R-L y=A-P; z=I-S ]]*/
xspace = read_mincdim( ncid , "xspace" ) ;
xspace.step = -xspace.step ; xspace.start = -xspace.start ;
xspace.afni_orient = (xspace.step < 0.0) ? ORI_L2R_TYPE : ORI_R2L_TYPE ;
yspace = read_mincdim( ncid , "yspace" ) ;
yspace.step = -yspace.step ; yspace.start = -yspace.start ;
yspace.afni_orient = (yspace.step < 0.0) ? ORI_P2A_TYPE : ORI_A2P_TYPE ;
zspace = read_mincdim( ncid , "zspace" ) ;
zspace.afni_orient = (zspace.step > 0.0) ? ORI_I2S_TYPE : ORI_S2I_TYPE ;
/*-- if don't have all 3 named dimensions, quit --*/
if( !xspace.good || !yspace.good || !zspace.good ){
nc_close(ncid) ; RETURN(NULL) ;
}
/*-- get info about the image data --*/
code = nc_inq_varid( ncid , "image" , &im_varid ) ;
if( code != NC_NOERR ){
EPR(code,pathname,"image varid"); nc_close(ncid); RETURN(NULL);
}
if( !AFNI_yesenv("AFNI_MINC_FLOATIZE") ){ /* determine type of dataset values */
code = nc_inq_vartype( ncid , im_varid , &im_type ) ;
if( code != NC_NOERR ){ EPR(code,pathname,"image vartype"); nc_close(ncid); RETURN(NULL); }
switch( im_type ){
default:
if( first_err ){ fprintf(stderr,"\n"); first_err=0; }
fprintf(stderr,"** MINC error: can't handle image type code %d\n",im_type);
nc_close(ncid); RETURN(NULL);
break ;
case NC_BYTE: datum = MRI_byte ; break ;
case NC_SHORT: datum = MRI_short ; break ;
case NC_FLOAT:
case NC_DOUBLE:
case NC_INT: datum = MRI_float ; break ;
}
} else { /* always read in as floats if */
datum = MRI_float ; /* AFNI_MINC_FLOATIZE is Yes */
}
/* we allow nD mages only for n=3 or n=4 */
code = nc_inq_varndims( ncid , im_varid , &im_rank ) ;
if( code != NC_NOERR ){
EPR(code,pathname,"image varndims"); nc_close(ncid); RETURN(NULL);
}
if( im_rank < 3 || im_rank > 4 ){
if( first_err ){ fprintf(stderr,"\n"); first_err=0; }
fprintf(stderr,"** MINC error: image rank=%d\n",im_rank) ;
nc_close(ncid); RETURN(NULL);
}
code = nc_inq_vardimid( ncid , im_varid , im_dimid ) ;
if( code != NC_NOERR ){
EPR(code,pathname,"image vardimid"); nc_close(ncid); RETURN(NULL);
}
/*-- find axes order --*/
/* fastest variation */
if( im_dimid[im_rank-1] == xspace.dimid ) xyz[0] = &xspace ;
else if( im_dimid[im_rank-1] == yspace.dimid ) xyz[0] = &yspace ;
else if( im_dimid[im_rank-1] == zspace.dimid ) xyz[0] = &zspace ;
else {
if( first_err ){ fprintf(stderr,"\n"); first_err=0; }
fprintf(stderr,"** MINC error: image dim[2] illegal\n") ;
nc_close(ncid) ; RETURN(NULL) ;
}
/* next fastest variation */
if( im_dimid[im_rank-2] == xspace.dimid ) xyz[1] = &xspace ;
else if( im_dimid[im_rank-2] == yspace.dimid ) xyz[1] = &yspace ;
else if( im_dimid[im_rank-2] == zspace.dimid ) xyz[1] = &zspace ;
else {
if( first_err ){ fprintf(stderr,"\n"); first_err=0; }
fprintf(stderr,"** MINC error: image dim[1] illegal\n") ;
nc_close(ncid) ; RETURN(NULL) ;
}
/* slowest variation */
if( im_dimid[im_rank-3] == xspace.dimid ) xyz[2] = &xspace ;
else if( im_dimid[im_rank-3] == yspace.dimid ) xyz[2] = &yspace ;
else if( im_dimid[im_rank-3] == zspace.dimid ) xyz[2] = &zspace ;
else {
if( first_err ){ fprintf(stderr,"\n"); first_err=0; }
fprintf(stderr,"** MINC error: image dim[0] illegal\n") ;
nc_close(ncid) ; RETURN(NULL) ;
}
/*-- determine the length of the 4th dimension, if any --*/
strcpy(tname,"MINC") ; /* default name for 4th dimension */
if( im_rank == 4 ){
char fname[NC_MAX_NAME+4] ;
code = nc_inq_dimlen( ncid , im_dimid[0] , &len ) ;
if( code != NC_NOERR ){
EPR(code,pathname,"dimid[0] dimlen"); nc_close(ncid); RETURN(NULL);
}
nvals = len ; /* number of volumes in this file */
/* get the name of this dimension */
code = nc_inq_dimname( ncid , im_dimid[0] , fname ) ;
if( code == NC_NOERR ){
MCW_strncpy(tname,fname,10) ;
}
/* get the stepsize of this dimension */
code = nc_get_att_float( ncid , im_dimid[0] , "step" , &dtime ) ;
if( code != NC_NOERR || dtime <= 0.0 ) dtime = 1.0 ;
}
/*-- make a dataset --*/
dset = EDIT_empty_copy(NULL) ;
ppp = THD_trailname(pathname,0) ; /* strip directory */
MCW_strncpy( prefix , ppp , THD_MAX_PREFIX ) ; /* to make prefix */
nxyz.ijk[0] = xyz[0]->len ; dxyz.xyz[0] = xyz[0]->step ; /* setup axes */
nxyz.ijk[1] = xyz[1]->len ; dxyz.xyz[1] = xyz[1]->step ;
nxyz.ijk[2] = xyz[2]->len ; dxyz.xyz[2] = xyz[2]->step ;
orixyz.ijk[0] = xyz[0]->afni_orient ; orgxyz.xyz[0] = xyz[0]->start ;
orixyz.ijk[1] = xyz[1]->afni_orient ; orgxyz.xyz[1] = xyz[1]->start ;
orixyz.ijk[2] = xyz[2]->afni_orient ; orgxyz.xyz[2] = xyz[2]->start ;
if( strstr(xyz[0]->spacetype,"alairach") != NULL ){ /* coord system */
iview = VIEW_TALAIRACH_TYPE ;
} else {
iview = VIEW_ORIGINAL_TYPE ;
}
dset->idcode.str[0] = 'M' ; /* overwrite 1st 3 bytes with something special */
dset->idcode.str[1] = 'N' ;
dset->idcode.str[2] = 'C' ;
MCW_hash_idcode( pathname , dset ) ; /* 06 May 2005 */
EDIT_dset_items( dset ,
ADN_prefix , prefix ,
ADN_datum_all , datum ,
ADN_nxyz , nxyz ,
ADN_xyzdel , dxyz ,
ADN_xyzorg , orgxyz ,
ADN_xyzorient , orixyz ,
ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
ADN_nvals , nvals ,
ADN_type , HEAD_ANAT_TYPE ,
ADN_view_type , iview ,
ADN_func_type , (nvals==1) ? ANAT_MRAN_TYPE
: ANAT_BUCK_TYPE ,
ADN_none ) ;
if( nvals > 9 ) /* pretend it is 3D+time */
EDIT_dset_items( dset ,
ADN_func_type, ANAT_EPI_TYPE ,
ADN_ntt , nvals ,
ADN_ttorg , 0.0 ,
ADN_ttdel , dtime ,
ADN_ttdur , 0.0 ,
ADN_tunits , UNITS_SEC_TYPE ,
ADN_none ) ;
/*-- flag to read data from disk using MINC/netCDF functions --*/
dset->dblk->diskptr->storage_mode = STORAGE_BY_MINC ;
strcpy( dset->dblk->diskptr->brick_name , pathname ) ;
for( ibr=0 ; ibr < nvals ; ibr++ ){ /* make sub-brick labels */
sprintf(prefix,"%s[%d]",tname,ibr) ;
EDIT_BRICK_LABEL( dset , ibr , prefix ) ;
}
/*-- read and save the global:history attribute --*/
#if 0
sprintf(prefix,"THD_open_minc(%s)",pathname) ;
tross_Append_History( dset , prefix ) ;
#endif
code = nc_inq_attlen( ncid , NC_GLOBAL , "history" , &len ) ;
if( code == NC_NOERR && len > 0 ){
ppp = AFMALL(char,len+4) ;
code = nc_get_att_text( ncid , NC_GLOBAL , "history" , ppp ) ;
if( code == NC_NOERR ){ /* should always happen */
ppp[len] = '\0' ;
tross_Append_History( dset , ppp ) ;
}
free(ppp) ;
}
/*-- close file and return the empty dataset */
nc_close(ncid) ; RETURN(dset) ;
}
/*-----------------------------------------------------------------
Load a MINC dataset's data into memory
(called from THD_load_datablock in thd_loaddblk.c)
-------------------------------------------------------------------*/
void THD_load_minc( THD_datablock *dblk )
{
THD_diskptr *dkptr ;
int im_varid , code , ncid ;
int nx,ny,nz,nxy,nxyz,nxyzv , nbad,ibr,nv, nslice ;
void *ptr ;
nc_type im_type=NC_SHORT ;
float im_valid_range[2], fac,denom , intop,inbot , outbot,outtop , sfac=1.0;
float *im_min=NULL , *im_max=NULL ;
int im_min_varid , im_max_varid , do_scale ;
ENTRY("THD_load_minc") ;
/*-- open input [these errors should never occur] --*/
if( !ISVALID_DATABLOCK(dblk) ||
dblk->diskptr->storage_mode != STORAGE_BY_MINC ||
dblk->brick == NULL ) EXRETURN ;
dkptr = dblk->diskptr ;
code = nc_open( dkptr->brick_name , NC_NOWRITE , &ncid ) ;
if( code != NC_NOERR ){ EPR(code,dkptr->brick_name,"open"); EXRETURN; }
code = nc_inq_varid( ncid , "image" , &im_varid ) ;
if( code != NC_NOERR ){
EPR(code,dkptr->brick_name,"image varid"); nc_close(ncid); EXRETURN;
}
/*-- allocate space for data --*/
nx = dkptr->dimsizes[0] ;
ny = dkptr->dimsizes[1] ; nxy = nx * ny ;
nz = dkptr->dimsizes[2] ; nxyz = nxy * nz ;
nv = dkptr->nvals ; nxyzv = nxyz * nv ; nslice = nz*nv ;
dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
/** malloc space for each brick separately **/
for( nbad=ibr=0 ; ibr < nv ; ibr++ ){
if( DBLK_ARRAY(dblk,ibr) == NULL ){
ptr = AFMALL(void, DBLK_BRICK_BYTES(dblk,ibr) ) ;
mri_fix_data_pointer( ptr , DBLK_BRICK(dblk,ibr) ) ;
if( ptr == NULL ) nbad++ ;
}
}
/** if couldn't get them all, take our ball and go home in a snit **/
if( nbad > 0 ){
fprintf(stderr,"\n** failed to malloc %d MINC bricks out of %d\n\a",nbad,nv) ;
for( ibr=0 ; ibr < nv ; ibr++ ){
if( DBLK_ARRAY(dblk,ibr) != NULL ){
free(DBLK_ARRAY(dblk,ibr)) ;
mri_fix_data_pointer( NULL , DBLK_BRICK(dblk,ibr) ) ;
}
}
nc_close(ncid) ; EXRETURN ;
}
/** get scaling attributes/variables **/
(void) nc_inq_vartype( ncid , im_varid , &im_type ) ;
/* N.B.: we don't scale if input data is stored as floats */
do_scale = (im_type==NC_BYTE || im_type==NC_SHORT || im_type==NC_INT) ;
if( do_scale && AFNI_noenv("AFNI_MINC_SLICESCALE") ) do_scale = 0 ;
code = nc_get_att_float( ncid,im_varid , "valid_range" , im_valid_range ) ;
/** if can't get valid_range, make something up **/
if( code != NC_NOERR || im_valid_range[0] >= im_valid_range[1] ){
im_valid_range[0] = 0.0 ;
switch( im_type ){
case NC_BYTE: im_valid_range[1] = 255.0 ; break ;
case NC_SHORT: im_valid_range[1] = 32767.0 ; break ;
case NC_INT: im_valid_range[1] = 2147483647.0 ; break ;
case NC_FLOAT:
case NC_DOUBLE: im_valid_range[1] = 1.0 ; break ;
}
}
inbot = im_valid_range[0] ;
intop = im_valid_range[1] ;
denom = intop - inbot ; /* always positive */
/** Get range of image (per 2D slice) to which valid_range will be scaled **/
/** Scaling will only be done if we get both image-min and image-max **/
if( do_scale ){
code = nc_inq_varid( ncid , "image-min" , &im_min_varid ) ;
if( code == NC_NOERR ){
im_min = (float *) calloc(sizeof(float),nslice) ;
code = nc_get_var_float( ncid , im_min_varid , im_min ) ;
if( code != NC_NOERR ){ free(im_min); im_min = NULL; }
}
/* got im_min ==> try for im_max */
if( im_min != NULL ){
code = nc_inq_varid( ncid , "image-max" , &im_max_varid ) ;
if( code == NC_NOERR ){
im_max = (float *) calloc(sizeof(float),nslice) ;
code = nc_get_var_float( ncid , im_max_varid , im_max ) ;
if( code != NC_NOERR ){
free(im_min); im_min = NULL; free(im_max); im_max = NULL;
}
}
}
/* 19 Mar 2003: make sure we don't scale out of range! */
if( im_max != NULL && MRI_IS_INT_TYPE(DBLK_BRICK_TYPE(dblk,0)) ){
float stop=0.0 , vbot,vtop ;
int ii ;
for( ii=0 ; ii < nslice ; ii++ ){
if( im_min[ii] < im_max[ii] ){
vbot = fabs(im_min[ii]) ; vtop = fabs(im_max[ii]) ;
vtop = MAX(vtop,vbot) ; stop = MAX(vtop,stop) ;
}
}
if( stop > MRI_TYPE_maxval[DBLK_BRICK_TYPE(dblk,0)] ){
sfac = MRI_TYPE_maxval[DBLK_BRICK_TYPE(dblk,0)] / stop ;
fprintf(stderr,"++ Scaling %s by %g to avoid overflow to %g of %s data!\n",
dkptr->brick_name, sfac, stop, MRI_TYPE_name[DBLK_BRICK_TYPE(dblk,0)] ) ;
} else if( stop == 0.0 ){
free(im_min); im_min = NULL; free(im_max); im_max = NULL;
}
}
} /* end of do_scale */
/** read data! **/
switch( DBLK_BRICK_TYPE(dblk,0) ){ /* output datum */
/* this case should never happen */
default:
fprintf(stderr,"\n** MINC illegal datum %d found\n",
DBLK_BRICK_TYPE(dblk,0) ) ;
for( ibr=0 ; ibr < nv ; ibr++ ){
if( DBLK_ARRAY(dblk,ibr) != NULL ){
free(DBLK_ARRAY(dblk,ibr)) ;
mri_fix_data_pointer( NULL , DBLK_BRICK(dblk,ibr) ) ;
}
}
break ;
/* data is stored as bytes in AFNI */
case MRI_byte:{
int kk,ii,qq ; byte *br ;
if( nv == 1 ){
code = nc_get_var_uchar( ncid , im_varid , DBLK_ARRAY(dblk,0) ) ;
} else {
size_t start[4] , count[4] ;
start[1] = start[2] = start[3] = 0 ;
count[0] = 1 ; count[1] = nz ; count[2] = ny ; count[3] = nx ;
for( ibr=0 ; ibr < nv ; ibr++ ){
start[0] = ibr ;
code = nc_get_vara_uchar( ncid,im_varid ,
start,count , DBLK_ARRAY(dblk,ibr) ) ;
}
}
if( code != NC_NOERR ) EPR(code,dkptr->brick_name,"image get_var_uchar") ;
if( im_max != NULL ){ /* must scale values */
for( ibr=0 ; ibr < nv ; ibr++ ){ /* loop over bricks */
br = DBLK_ARRAY(dblk,ibr) ;
for( kk=0 ; kk < nz ; kk++ ){ /* loop over slices */
outbot = im_min[kk+ibr*nz] ;
outtop = im_max[kk+ibr*nz] ;
if( outbot >= outtop) continue ; /* skip */
fac = (outtop-outbot) / denom ;
qq = kk*nxy ; outbot += 0.499 ;
for( ii=0 ; ii < nxy ; ii++ ) /* scale slice */
br[ii+qq] = (byte) ((fac*(br[ii+qq]-inbot) + outbot)*sfac) ;
}
}
}
}
break ;
/* data is stored as shorts in AFNI */
case MRI_short:{
int kk,ii,qq ; short *br ;
if( nv == 1 ){
code = nc_get_var_short( ncid , im_varid , DBLK_ARRAY(dblk,0) ) ;
} else {
size_t start[4] , count[4] ;
start[1] = start[2] = start[3] = 0 ;
count[0] = 1 ; count[1] = nz ; count[2] = ny ; count[3] = nx ;
for( ibr=0 ; ibr < nv ; ibr++ ){
start[0] = ibr ;
code = nc_get_vara_short( ncid,im_varid ,
start,count , DBLK_ARRAY(dblk,ibr) ) ;
}
}
if( code != NC_NOERR ) EPR(code,dkptr->brick_name,"image get_var_short") ;
if( im_max != NULL ){ /* must scale values */
for( ibr=0 ; ibr < nv ; ibr++ ){ /* loop over bricks */
br = DBLK_ARRAY(dblk,ibr) ;
for( kk=0 ; kk < nz ; kk++ ){ /* loop over slices */
outbot = im_min[kk+ibr*nz] ;
outtop = im_max[kk+ibr*nz] ;
if( outbot >= outtop) continue ; /* skip */
fac = (outtop-outbot) / denom ;
qq = kk*nxy ; outbot += 0.499 ;
for( ii=0 ; ii < nxy ; ii++ ) /* scale slice */
br[ii+qq] = (short) ((fac*(br[ii+qq]-inbot) + outbot)*sfac) ;
}
}
}
}
break ;
/* data is stored as floats in AFNI */
case MRI_float:{
int kk,ii,qq ; float *br ;
if( nv == 1 ){
code = nc_get_var_float( ncid , im_varid , DBLK_ARRAY(dblk,0) ) ;
} else {
size_t start[4] , count[4] ;
start[1] = start[2] = start[3] = 0 ;
count[0] = 1 ; count[1] = nz ; count[2] = ny ; count[3] = nx ;
for( ibr=0 ; ibr < nv ; ibr++ ){
start[0] = ibr ;
code = nc_get_vara_float( ncid,im_varid ,
start,count , DBLK_ARRAY(dblk,ibr) ) ;
}
}
if( code != NC_NOERR ) EPR(code,dkptr->brick_name,"image get_var_float") ;
if( im_type == NC_BYTE ){ /* fix sign of bytes */
for( ibr=0 ; ibr < nv ; ibr++ ){ /* loop over bricks */
br = DBLK_ARRAY(dblk,ibr) ;
for( ii=0 ; ii < nxyz ; ii++ )
if( br[ii] < 0.0 ) br[ii] += 256.0 ;
}
}
if( im_max != NULL ){ /* must scale values */
for( ibr=0 ; ibr < nv ; ibr++ ){ /* loop over bricks */
br = DBLK_ARRAY(dblk,ibr) ;
for( kk=0 ; kk < nz ; kk++ ){ /* loop over slices */
outbot = im_min[kk+ibr*nz] ;
outtop = im_max[kk+ibr*nz] ;
if( outbot >= outtop) continue ; /* skip */
fac = (outtop-outbot) / denom ;
qq = kk*nxy ;
for( ii=0 ; ii < nxy ; ii++ ) /* scale slice */
br[ii+qq] = (fac*(br[ii+qq]-inbot) + outbot) ;
}
}
}
}
break ;
} /* end of switch on output datum */
/*-- throw away the trash and return --*/
if( im_min != NULL ) free(im_min) ;
if( im_max != NULL ) free(im_max) ;
nc_close(ncid) ; EXRETURN ;
}
syntax highlighted by Code2HTML, v. 0.9.1