/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 1994-2000, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
#include "mrilib.h"
/*------------------------------------------------------------------------*/
static float min_float( int n , float *x ) /* 25 Feb 2005 */
{
float m=0.0 ; int i ;
if( n < 1 || x == NULL ) return m ;
m = x[0] ;
for( i=1 ; i < n ; i++ ) if( m > x[i] ) m = x[i] ;
return m ;
}
/*------------------------------------------------------------------------*/
static float max_float( int n , float *x ) /* 24 Feb 2005 */
{
float m=0.0 ; int i ;
if( n < 1 || x == NULL ) return m ;
m = x[0] ;
for( i=1 ; i < n ; i++ ) if( m < x[i] ) m = x[i] ;
return m ;
}
/*-----------------------------------------------------------------------
A quickie. (Not so quick any more, with all the options added.)
-------------------------------------------------------------------------*/
int main( int argc , char * argv[] )
{
int narg , nvox , ii , mcount , iv , mc ;
THD_3dim_dataset * mask_dset=NULL , * input_dset=NULL ;
float mask_bot=666.0 , mask_top=-666.0 ;
byte * mmm = NULL ;
int dumpit = 0 , sigmait = 0 ;
int miv = 0 ; /* 06 Aug 1998 */
int div = -1 , div_bot,div_top , drange=0; /* 16 Sep 1998 */
float data_bot=666.0 , data_top=-666.0 ;
int indump = 0 ; /* 19 Aug 1999 */
int pslice=-1 , qslice=-1 , nxy,nz ; /* 15 Sep 1999 */
int quiet=0 ; /* 23 Nov 1999 */
int medianit = 0 ; /* 06 Jul 2003 */
int maxit = 0 ; /* 24 Feb 2005 */
int minit = 0 ; /* 25 Feb 2005 */
float *exar ; /* 06 Jul 2003 */
char *sname = "Average" ; /* 06 Jul 2003 */
int self_mask = 0 ; /* 06 Dec 2004 */
if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
printf("Usage: 3dmaskave [options] dataset\n"
"\n"
"Computes average of all voxels in the input dataset\n"
"which satisfy the criterion in the options list.\n"
"If no options are given, then all voxels are included.\n"
"\n" /* examples: 13 Mar 2006 [rickr] */
"------------------------------------------------------------\n"
"Examples:\n"
"\n"
"1. compute the average timeseries in epi_r1+orig, over voxels\n"
" that are set (any non-zero value) in the dataset, ROI+orig:\n"
"\n"
" 3dmaskave -mask ROI+orig epi_r1+orig\n"
"\n"
"2. restrict the ROI to values of 3 or 4, and save (redirect)\n"
" the output to the text file run1_roi_34.txt:\n"
"\n"
" 3dmaskave -mask ROI+orig -quiet -mrange 3 4 \\\n"
" epi_r1+orig > run1_roi_34.txt\n"
"------------------------------------------------------------\n"
"\n"
"Options:\n"
" -mask mset Means to use the dataset 'mset' as a mask:\n"
" Only voxels with nonzero values in 'mset'\n"
" will be averaged from 'dataset'. Note\n"
" that the mask dataset and the input dataset\n"
" must have the same number of voxels.\n"
" SPECIAL CASE: If 'mset' is the string 'SELF',\n"
" then the input dataset will be\n"
" used to mask itself. That is,\n"
" only nonzero voxels from the\n"
" #miv sub-brick will be used.\n"
" -mindex miv Means to use sub-brick #'miv' from the mask\n"
" dataset. If not given, miv=0.\n"
" -mrange a b Means to further restrict the voxels from\n"
" 'mset' so that only those mask values\n"
" between 'a' and 'b' (inclusive) will\n"
" be used. If this option is not given,\n"
" all nonzero values from 'mset' are used.\n"
" Note that if a voxel is zero in 'mset', then\n"
" it won't be included, even if a < 0 < b.\n"
"\n"
" -dindex div Means to use sub-brick #'div' from the dataset.\n"
" If not given, all sub-bricks will be processed.\n"
" -drange a b Means to only include voxels from the dataset whose\n"
" values fall in the range 'a' to 'b' (inclusive).\n"
" Otherwise, all voxel values are included.\n"
"\n"
" -slices p q Means to only included voxels from the dataset\n"
" whose slice numbers are in the range 'p' to 'q'\n"
" (inclusive). Slice numbers range from 0 to\n"
" NZ-1, where NZ can be determined from the output\n"
" of program 3dinfo. The default is to include\n"
" data from all slices.\n"
" [There is no provision for geometrical voxel]\n"
" [selection except in the slice (z) direction]\n"
"\n"
" -sigma Means to compute the standard deviation as well\n"
" as the mean.\n"
" -median Means to compute the median instead of the mean.\n"
" -max Means to compute the max instead of the mean.\n"
" -min Means to compute the min instead of the mean.\n"
" (-sigma is ignored with -median, -max, or -min)\n"
" -dump Means to print out all the voxel values that\n"
" go into the average.\n"
" -udump Means to print out all the voxel values that\n"
" go into the average, UNSCALED by any internal\n"
" factors.\n"
" N.B.: the scale factors for a sub-brick\n"
" can be found using program 3dinfo.\n"
" -indump Means to print out the voxel indexes (i,j,k) for\n"
" each dumped voxel. Has no effect if -dump\n"
" or -udump is not also used.\n"
" N.B.: if nx,ny,nz are the number of voxels in\n"
" each direction, then the array offset\n"
" in the brick corresponding to (i,j,k)\n"
" is i+j*nx+k*nx*ny.\n"
" -q or\n"
" -quiet Means to print only the minimal results.\n"
" This is useful if you want to create a *.1D file.\n"
"\n"
"The output is printed to stdout (the terminal), and can be\n"
"saved to a file using the usual redirection operation '>'.\n"
) ;
printf("\n" MASTER_SHORTHELP_STRING ) ;
exit(0) ;
}
mainENTRY("3dmaskave main"); machdep(); AFNI_logger("3dmaskave",argc,argv);
PRINT_VERSION("3dmaskave") ;
/* scan argument list */
narg = 1 ;
while( narg < argc && argv[narg][0] == '-' ){
if( strcmp(argv[narg],"-q") == 0 || strcmp(argv[narg],"-quiet") == 0 ){
quiet++ ;
narg++ ; continue ;
}
if( strncmp(argv[narg],"-mask",5) == 0 ){
if( mask_dset != NULL || self_mask ){
fprintf(stderr,"*** Cannot have two -mask options!\n") ; exit(1) ;
}
if( narg+1 >= argc ){
fprintf(stderr,"*** -mask option requires a following argument!\n") ;
exit(1) ;
}
narg++ ;
if( strcmp(argv[narg],"SELF") == 0 ){
self_mask = 1 ;
} else {
mask_dset = THD_open_dataset( argv[narg] ) ;
if( mask_dset == NULL ){
fprintf(stderr,"*** Cannot open mask dataset!\n") ; exit(1) ;
}
if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex ){
fprintf(stderr,"*** Cannot deal with complex-valued mask dataset!\n") ;
exit(1) ;
} else if( DSET_BRICK_TYPE(mask_dset,0) == MRI_rgb ){
fprintf(stderr,"*** Cannot deal with rgb-valued mask dataset!\n") ;
exit(1) ;
}
}
narg++ ; continue ;
}
/* 06 Aug 1998 */
if( strncmp(argv[narg],"-mindex",5) == 0 ){
if( narg+1 >= argc ){
fprintf(stderr,"*** -mindex option needs 1 following argument!\n") ; exit(1) ;
}
miv = (int) strtod( argv[++narg] , NULL ) ;
if( miv < 0 ){
fprintf(stderr,"*** -mindex value is negative!\n") ; exit(1) ;
}
narg++ ; continue ;
}
if( strncmp(argv[narg],"-mrange",5) == 0 ){
if( narg+2 >= argc ){
fprintf(stderr,"*** -mrange option requires 2 following arguments!\n") ; exit(1) ;
}
mask_bot = strtod( argv[++narg] , NULL ) ;
mask_top = strtod( argv[++narg] , NULL ) ;
if( mask_top < mask_top ){
fprintf(stderr,"*** -mrange inputs are illegal!\n") ; exit(1) ;
}
narg++ ; continue ;
}
/* 16 Sep 1998 */
if( strncmp(argv[narg],"-dindex",5) == 0 ){
if( narg+1 >= argc ){
fprintf(stderr,"*** -dindex option needs 1 following argument!\n") ; exit(1) ;
}
div = (int) strtod( argv[++narg] , NULL ) ;
if( div < 0 ){
fprintf(stderr,"*** -dindex value is negative!\n") ; exit(1) ;
}
narg++ ; continue ;
}
if( strncmp(argv[narg],"-drange",5) == 0 ){
if( narg+2 >= argc ){
fprintf(stderr,
"*** -drange option requires 2 following arguments!\n") ;
exit(1) ;
}
data_bot = strtod( argv[++narg] , NULL ) ;
data_top = strtod( argv[++narg] , NULL ) ;
if( data_top < data_top ){
fprintf(stderr,"*** -mrange inputs are illegal!\n") ; exit(1) ;
}
drange = 1 ;
narg++ ; continue ;
}
if( strncmp(argv[narg],"-slices",5) == 0 ){ /* 15 Sep 1999 */
if( narg+2 >= argc ){
fprintf(stderr,
"*** -slices option requires 2 following arguments!\n") ;
exit(1) ;
}
pslice = (int) strtod( argv[++narg] , NULL ) ;
qslice = (int) strtod( argv[++narg] , NULL ) ;
if( pslice < 0 || qslice < 0 || qslice < pslice ){
fprintf(stderr, "*** Illegal values after -slices!\n") ;
exit(1) ;
}
narg++ ; continue ;
}
if( strncmp(argv[narg],"-dump",5) == 0 ){
dumpit = 1 ;
narg++ ; continue ;
}
if( strncmp(argv[narg],"-udump",5) == 0 ){
dumpit = 2 ;
narg++ ; continue ;
}
if( strncmp(argv[narg],"-indump",5) == 0 ){ /* 19 Aug 1999 */
indump = 1 ;
narg++ ; continue ;
}
if( strncmp(argv[narg],"-sigma",5) == 0 ){
sigmait = 2 ;
narg++ ; continue ;
}
if( strncmp(argv[narg],"-median",5) == 0 ){
medianit = 1 ; maxit = 0 ; minit = 0 ;
narg++ ; continue ;
}
if( strncmp(argv[narg],"-max",4) == 0 ){ /* 24 Feb 2005 */
maxit = 1 ; medianit = 0 ; minit = 0 ;
narg++ ; continue ;
}
if( strncmp(argv[narg],"-min",4) == 0 ){ /* 25 Feb 2005 */
maxit = 0 ; medianit = 0 ; minit = 1 ;
narg++ ; continue ;
}
fprintf(stderr,"*** Unknown option: %s\n",argv[narg]) ; exit(1) ;
}
if( medianit ){ sigmait = 0; sname = "Median"; }
if( maxit ){ sigmait = 0; sname = "Max" ; } /* 24 Feb 2005 */
if( minit ){ sigmait = 0; sname = "Min" ; } /* 25 Feb 2005 */
/* should have one more argument */
if( narg >= argc ){
fprintf(stderr,"*** No input dataset!?\n") ; exit(1) ;
}
#if 0
if( dumpit && mask_dset == NULL ){
fprintf(stderr,"*** Can't use dump option without -mask!\n") ; exit(1) ;
}
#endif
/* read input dataset */
input_dset = THD_open_dataset( argv[narg] ) ;
if( input_dset == NULL ){
fprintf(stderr,"*** Cannot open input dataset!\n") ; exit(1) ;
}
if( self_mask ) mask_dset = input_dset ; /* 06 Dec 2004 */
if( miv > 0 ){ /* 06 Aug 1998 */
if( mask_dset == NULL ){
fprintf(stderr,"*** -mindex option used without -mask!\n") ; exit(1) ;
}
if( miv >= DSET_NVALS(mask_dset) ){
fprintf(stderr,"*** -mindex value is too large!\n") ; exit(1) ;
}
}
if( DSET_BRICK_TYPE(input_dset,0) == MRI_complex ){
fprintf(stderr,"*** Cannot deal with complex-valued input dataset!\n") ; exit(1) ;
}
if( div >= DSET_NVALS(input_dset) ){
fprintf(stderr,"*** Not enough sub-bricks in dataset for -dindex %d!\n",div) ;
exit(1) ;
}
if( pslice >= 0 ){
nxy = DSET_NX(input_dset) * DSET_NY(input_dset) ;
nz = DSET_NZ(input_dset) ;
if( qslice >= nz ){
fprintf(stderr,
"*** There are only %d slices in the input dataset!\n",nz) ;
exit(1) ;
}
if( pslice == 0 && qslice == nz-1 )
fprintf(stderr,"+++ -slice option says to use all slices!?\n") ;
}
nvox = DSET_NVOX(input_dset) ;
/* make a byte mask from mask dataset */
mmm = (byte *) malloc( sizeof(byte) * nvox ) ;
if( mmm == NULL ){
fprintf(stderr,"*** Cannot malloc workspace!\n") ; exit(1) ;
}
if( mask_dset != NULL ){
if( DSET_NVOX(mask_dset) != nvox ){
fprintf(stderr,"*** Input and mask datasets are not same dimensions!\n") ;
exit(1) ;
}
DSET_load(mask_dset) ;
if( DSET_ARRAY(mask_dset,miv) == NULL ){
fprintf(stderr,"*** Cannot read in mask dataset BRIK!\n") ; exit(1) ;
}
switch( DSET_BRICK_TYPE(mask_dset,miv) ){
default:
fprintf(stderr,"*** Cannot deal with mask dataset datum!\n") ; exit(1) ;
case MRI_short:{
short mbot , mtop ;
short * mar = (short *) DSET_ARRAY(mask_dset,miv) ;
float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ;
if( mfac == 0.0 ) mfac = 1.0 ;
if( mask_bot <= mask_top ){
mbot = SHORTIZE(mask_bot/mfac) ;
mtop = SHORTIZE(mask_top/mfac) ;
} else {
mbot = (short) -MRI_TYPE_maxval[MRI_short] ;
mtop = (short) MRI_TYPE_maxval[MRI_short] ;
}
for( mcount=0,ii=0 ; ii < nvox ; ii++ )
if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 ){ mmm[ii] = 1 ; mcount++ ; }
else { mmm[ii] = 0 ; }
}
break ;
case MRI_byte:{
byte mbot , mtop ;
byte * mar = (byte *) DSET_ARRAY(mask_dset,miv) ;
float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ;
if( mfac == 0.0 ) mfac = 1.0 ;
if( mask_bot <= mask_top ){
mbot = BYTEIZE(mask_bot/mfac) ;
mtop = BYTEIZE(mask_top/mfac) ;
if( mtop == 0 ){
fprintf(stderr,"*** Illegal mask range for mask dataset of bytes.\n") ; exit(1) ;
}
} else {
mbot = 0 ;
mtop = (byte) MRI_TYPE_maxval[MRI_short] ;
}
for( mcount=0,ii=0 ; ii < nvox ; ii++ )
if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 ){ mmm[ii] = 1 ; mcount++ ; }
else { mmm[ii] = 0 ; }
}
break ;
case MRI_float:{
float mbot , mtop ;
float * mar = (float *) DSET_ARRAY(mask_dset,miv) ;
float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ;
if( mfac == 0.0 ) mfac = 1.0 ;
if( mask_bot <= mask_top ){
mbot = (float) (mask_bot/mfac) ;
mtop = (float) (mask_top/mfac) ;
} else {
mbot = -WAY_BIG ;
mtop = WAY_BIG ;
}
for( mcount=0,ii=0 ; ii < nvox ; ii++ )
if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 ){ mmm[ii] = 1 ; mcount++ ; }
else { mmm[ii] = 0 ; }
}
break ;
}
if( !self_mask ) DSET_unload(mask_dset) ;
if( mcount == 0 ){
fprintf(stderr,"*** No voxels survive the masking operation\n"); exit(1);
}
fprintf(stderr,"+++ %d voxels survive the mask\n",mcount) ;
} else {
mcount = nvox ;
memset( mmm , 1 , mcount ) ;
fprintf(stderr,"+++ %d voxels in the entire dataset (no mask)\n",mcount) ;
}
if( pslice >= 0 ){ /* 15 Sep 1999 */
int kz , ibot ;
mcount = 0 ;
for( kz=0 ; kz < nz ; kz++ ){ /* loop over all slices */
ibot = kz*nxy ; /* base index for this slice */
if( kz >= pslice && kz <= qslice ){ /* keepers => recount */
for( ii=0 ; ii < nxy ; ii++ )
if( mmm[ii+ibot] ) mcount++ ;
} else { /* throw them back */
for( ii=0 ; ii < nxy ; ii++ )
mmm[ii+ibot] = 0 ;
}
}
if( mcount == 0 ){
fprintf(stderr,"*** No voxels survive the slicing operation\n"); exit(1);
}
fprintf(stderr,"+++ %d voxels survive the slicing\n",mcount) ;
}
if( mcount < 2 && sigmait ){
fprintf(stderr,"+++ [too few voxels; cannot compute sigma]\n") ;
sigmait = 0 ;
}
DSET_load(input_dset) ;
if( DSET_ARRAY(input_dset,0) == NULL ){
fprintf(stderr,"*** Cannot read in input dataset BRIK!\n") ; exit(1) ;
}
/* loop over input sub-bricks */
if( div < 0 ){
div_bot = 0 ; div_top = DSET_NVALS(input_dset) ;
} else {
div_bot = div ; div_top = div+1 ;
}
exar = (float *) malloc( sizeof(float) * mcount ) ; /* 06 Jul 2003 */
for( iv=div_bot ; iv < div_top ; iv++ ){
switch( DSET_BRICK_TYPE(input_dset,iv) ){
default:
printf("*** Illegal sub-brick datum at %d\n",iv) ;
break ;
#define INRANGE(i) ( !drange || ( mfac*bar[i] >= data_bot && mfac*bar[i] <= data_top ) )
#define GOOD(i) ( mmm[i] && INRANGE(i) )
case MRI_short:{
short * bar = (short *) DSET_ARRAY(input_dset,iv) ;
double sum = 0.0 , sigma = 0.0 ;
float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
if( mfac == 0.0 || dumpit == 2 ) mfac = 1.0 ;
if( dumpit ){
int noscal = (dumpit==2) || (mfac==1.0) ;
printf("+++ Dump for sub-brick %d:\n",iv) ;
for( ii=0 ; ii < nvox ; ii++ ) if( GOOD(ii) ){
if( noscal ) printf(" %d",bar[ii]) ;
else printf(" %g",bar[ii]*mfac) ;
if( indump )
printf(" (%d,%d,%d)",
DSET_index_to_ix(input_dset,ii),
DSET_index_to_jy(input_dset,ii),
DSET_index_to_kz(input_dset,ii) ) ;
printf("\n") ;
}
}
for( ii=mc=0 ; ii < nvox ; ii++ )
if( GOOD(ii) ){ sum += bar[ii] ; exar[mc++] = bar[ii] ; }
if( mc > 0 ) sum = sum / mc ;
if( sigmait && mc > 1 ){
for( ii=0 ; ii < nvox ; ii++ ) if( GOOD(ii) ) sigma += SQR(bar[ii]-sum) ;
sigma = mfac * sqrt( sigma/(mc-1) ) ;
}
if( medianit ) sum = qmed_float( mc , exar ) ;
else if( maxit ) sum = max_float ( mc , exar ) ;
else if( minit ) sum = min_float ( mc , exar ) ;
sum = mfac * sum ;
if( dumpit ) printf("+++ %s = %g",sname,sum) ;
else printf("%g",sum) ;
if( sigmait ){
if( dumpit ) printf(" Sigma = %g",sigma) ;
else printf(" %g",sigma) ;
}
if( !quiet ) printf(" [%d voxels]\n",mc) ;
else printf("\n") ;
}
break ;
case MRI_byte:{
byte * bar = (byte *) DSET_ARRAY(input_dset,iv) ;
double sum = 0.0 , sigma = 0.0 ;
float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
if( mfac == 0.0 || dumpit == 2 ) mfac = 1.0 ;
if( dumpit ){
int noscal = (dumpit==2) || (mfac==1.0) ;
printf("+++ Dump for sub-brick %d:\n",iv) ;
for( ii=0 ; ii < nvox ; ii++ ) if( GOOD(ii) ){
if( noscal ) printf(" %d",bar[ii]) ;
else printf(" %g",bar[ii]*mfac) ;
if( indump )
printf(" (%d,%d,%d)",
DSET_index_to_ix(input_dset,ii),
DSET_index_to_jy(input_dset,ii),
DSET_index_to_kz(input_dset,ii) ) ;
printf("\n") ;
}
}
for( ii=mc=0 ; ii < nvox ; ii++ )
if( GOOD(ii) ){ sum += bar[ii] ; exar[mc++] = bar[ii] ; }
if( mc > 0 ) sum = sum / mc ;
if( sigmait && mc > 1 ){
for( ii=0 ; ii < nvox ; ii++ ) if( GOOD(ii) ) sigma += SQR(bar[ii]-sum) ;
sigma = mfac * sqrt( sigma/(mc-1) ) ;
}
if( medianit ) sum = qmed_float( mc , exar ) ;
else if( maxit ) sum = max_float ( mc , exar ) ;
else if( minit ) sum = min_float ( mc , exar ) ;
sum = mfac * sum ;
if( dumpit ) printf("+++ %s = %g",sname,sum) ;
else printf("%g",sum) ;
if( sigmait ){
if( dumpit ) printf(" Sigma = %g",sigma) ;
else printf(" %g",sigma) ;
}
if( !quiet ) printf(" [%d voxels]\n",mc) ;
else printf("\n") ;
}
break ;
case MRI_float:{
float * bar = (float *) DSET_ARRAY(input_dset,iv) ;
double sum = 0.0 , sigma = 0.0 ;
float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
if( mfac == 0.0 || dumpit == 2 ) mfac = 1.0 ;
if( dumpit ){
int noscal = (dumpit==2) || (mfac==1.0) ;
printf("+++ Dump for sub-brick %d:\n",iv) ;
for( ii=0 ; ii < nvox ; ii++ ) if( GOOD(ii) ){
if( noscal ) printf(" %g",bar[ii]) ;
else printf(" %g",bar[ii]*mfac) ;
if( indump )
printf(" (%d,%d,%d)",
DSET_index_to_ix(input_dset,ii),
DSET_index_to_jy(input_dset,ii),
DSET_index_to_kz(input_dset,ii) ) ;
printf("\n") ;
}
}
for( ii=mc=0 ; ii < nvox ; ii++ )
if( GOOD(ii) ){ sum += bar[ii] ; exar[mc++] = bar[ii] ; }
if( mc > 0 ) sum = sum / mc ;
if( sigmait && mc > 1 ){
for( ii=0 ; ii < nvox ; ii++ ) if( GOOD(ii) ) sigma += SQR(bar[ii]-sum) ;
sigma = mfac * sqrt( sigma/(mc-1) ) ;
}
if( medianit ) sum = qmed_float( mc , exar ) ;
else if( maxit ) sum = max_float ( mc , exar ) ;
else if( minit ) sum = min_float ( mc , exar ) ;
sum = mfac * sum ;
if( dumpit ) printf("+++ %s = %g",sname,sum) ;
else printf("%g",sum) ;
if( sigmait ){
if( dumpit ) printf(" Sigma = %g",sigma) ;
else printf(" %g",sigma) ;
}
if( !quiet ) printf(" [%d voxels]\n",mc) ;
else printf("\n") ;
}
break ;
}
}
free(exar) ; DSET_delete(input_dset) ;
exit(0) ;
}
syntax highlighted by Code2HTML, v. 0.9.1