/*****************************************************************************
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 "afni.h"
#ifndef ALLOW_PLUGINS
# error "Plugins not properly set up -- see machdep.h"
#endif
/***********************************************************************
Plugin to compute power spectrum of a 3D+time dataset.
This is a moderately complex example, showing how to deal
with different data types on input and output.
************************************************************************/
/*------------- string to 'help' the user -------------*/
static char helpstring[] =
" Purpose: Compute 'Power Spectrum' of a 3D+time dataset.\n"
" Input items are:\n"
" Input = 3D+time dataset to analyze\n"
"\n"
" Output: Prefix = Filename prefix for new dataset\n"
" Datum = How to store results\n"
"\n"
" Ignore Count = How many points to ignore at start\n"
" Taper Percent = Amount of data to taper (Hamming)\n"
" FFT Length = Fourier transform size to use [N]\n"
" (If N > size of data, data will be zero)\n"
" (padded. 'shortest' means to use N just)\n"
" (above the length of the time series. )\n"
"\n"
" The output dataset will be stored in the 3D+time format, with\n"
" the 'time' index actually being frequency. The frequency grid\n"
" spacing will be 1/(N*dt), where N=FFT length and dt = input\n"
" dataset time spacing.\n"
"\n"
" The method used is the simplest known: squared periodogram.\n"
" A single FFT is done (i.e., each point has DOF=2.)\n"
;
/*------------- strings for output format -------------*/
static char * type_strings[]
= { "as Input" , "Byte" , "Short" , "Float" } ;
#define NUM_TYPE_STRINGS (sizeof(type_strings)/sizeof(char *))
/*------------- strings for FFT length -------------*/
static char * fft_strings[] =
#if 0
{ "shortest", "32", "64", "128", "256", "512", "1024", "2048", "4096" } ;
#else
/* 3* 15* 2** 5* */
{ "shortest", "32" , "48" , "60" , "64" , "80" ,
"96" , "120" , "128" , "160" ,
"192" , "240" , "256" , "320" ,
"384" , "480" , "512" , "640" ,
"768" , "960" , "1024" , "1280" ,
"1536" , "1920" , "2048" } ;
#endif
#define NUM_FFT_STRINGS (sizeof(fft_strings)/sizeof(char *))
/*--------------- prototypes for internal routines ---------------*/
char * POWER_main( PLUGIN_interface * ) ; /* the entry point */
#undef ALLOW_TESTING
#ifdef ALLOW_TESTING
PLUGIN_interface * TEST_init(void) ;
char * TEST_main( PLUGIN_interface * ) ; /* the entry point */
#endif
/***********************************************************************
Set up the interface to the user:
1) Create a new interface using "PLUTO_new_interface";
2) For each line of inputs, create the line with "PLUTO_add_option"
(this line of inputs can be optional or mandatory);
3) For each item on the line, create the item with
"PLUTO_add_dataset" for a dataset chooser,
"PLUTO_add_string" for a string chooser,
"PLUTO_add_number" for a number chooser.
************************************************************************/
DEFINE_PLUGIN_PROTOTYPE
PLUGIN_interface * PLUGIN_init( int ncall )
{
PLUGIN_interface * plint ; /* will be the output of this routine */
if( ncall > 1 ) return NULL ; /* two interfaces */
#ifdef ALLOW_TESTING
if( ncall == 1 ) return TEST_init() ;
#else
if( ncall == 1 ) return NULL ;
#endif
/*---------------- set titles and call point ----------------*/
plint = PLUTO_new_interface( "Power Spectrum" ,
"Power Spectrum of a 3D+time Dataset" ,
helpstring ,
PLUGIN_CALL_VIA_MENU , POWER_main ) ;
PLUTO_add_hint( plint , "Power Spectrum of a 3D+time Dataset" ) ;
PLUTO_set_sequence( plint , "A:newdset:statistics" ) ;
/*--------- 1st line: Input dataset ---------*/
PLUTO_add_option( plint ,
"Input" , /* label at left of input line */
"Input" , /* tag to return to plugin */
TRUE /* is this mandatory? */
) ;
PLUTO_add_dataset( plint ,
"---->>" , /* label next to button */
ANAT_ALL_MASK , /* take any anat datasets */
FUNC_FIM_MASK , /* only allow fim funcs */
DIMEN_4D_MASK | /* need 3D+time datasets */
BRICK_ALLREAL_MASK /* need real-valued datasets */
) ;
/*---------- 2nd line: Output dataset ----------*/
PLUTO_add_option( plint ,
"Output" , /* label at left of input line */
"Output" , /* tag to return to plugin */
TRUE /* is this mandatory? */
) ;
PLUTO_add_string( plint ,
"Prefix" , /* label next to textfield */
0,NULL , /* no fixed strings to choose among */
19 /* 19 spaces for typing in value */
) ;
PLUTO_add_string( plint ,
"Datum" , /* label next to chooser button */
NUM_TYPE_STRINGS , /* number of strings to choose among */
type_strings , /* list of strings to choose among */
0 /* index of default string */
) ;
/*--------- Other lines: Parameters ---------*/
PLUTO_add_option( plint , "Ignore" , "Ignore" , TRUE ) ;
PLUTO_add_number( plint ,
"Count" , /* label next to chooser */
0 , /* smallest possible value */
999 , /* largest possible value */
0 , /* decimal shift (none in this case) */
4 , /* default value */
TRUE /* allow user to edit value? */
) ;
PLUTO_add_option( plint , "Taper" , "Taper" , TRUE ) ;
PLUTO_add_number( plint ,
"Percent" , /* label next to chooser */
0 , /* smallest possible value */
10 , /* largest possible value */
-1 , /* decimal shift (1 right == 0 to 100) */
0 , /* default value (with shift == 0) */
FALSE /* allow user to edit value? */
) ;
PLUTO_add_option( plint , "FFT" , "FFT" , TRUE ) ;
PLUTO_add_string( plint ,
"Length" , /* label next to chooser */
NUM_FFT_STRINGS , /* number of strings to choose among */
fft_strings , /* list of strings to choose among */
0 /* index of default string */
) ;
/*--------- done with interface setup ---------*/
return plint ;
}
/***************************************************************************
Main routine for this plugin (will be called from AFNI).
If the return string is not NULL, some error transpired, and
AFNI will popup the return string in a message box.
****************************************************************************/
/*------------------ macros to return workspace at exit -------------------*/
#undef FREEUP
#define FREEUP(x) do{ if((x) != NULL){free((x)); (x)=NULL;} } while(0)
#define FREE_WORKSPACE \
do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ; \
FREEUP(fout) ; FREEUP(cxar) ; FREEUP(tar) ; \
FREEUP(fxar) ; FREEUP(fyar) ; FREEUP(dtr) ; \
} while(0)
/*-------------------------------------------------------------------------*/
char * POWER_main( PLUGIN_interface * plint )
{
MCW_idcode * idc ; /* input dataset idcode */
THD_3dim_dataset * old_dset , * new_dset ; /* input and output datasets */
char * new_prefix , * str ; /* strings from user */
int new_datum , ignore , nfft , ninp , /* control parameters */
old_datum , nuse , ntaper , ktbot ;
float taper ;
byte ** bptr = NULL ; /* one of these will be the array of */
short ** sptr = NULL ; /* pointers to input dataset sub-bricks */
float ** fptr = NULL ; /* (depending on input datum type) */
complex * cxar = NULL ; /* will be array of data to FFT */
float * fxar = NULL ; /* array loaded from input dataset */
float * fyar = NULL ; /* array loaded from input dataset */
float ** fout = NULL ; /* will be array of output floats */
float * tar = NULL ; /* will be array of taper coefficients */
float * dtr = NULL ; /* will be array of detrending coeff */
float dfreq , pfact , phi , xr,xi , yr,yi ;
float x0,x1 , y0,y1 , d0fac,d1fac ;
int nfreq , nvox , perc , new_units ;
int istr , ii,iip , ibot,itop , kk , icx ; /* temp variables */
/*--------------------------------------------------------------------*/
/*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
/*--------- go to first input line ---------*/
PLUTO_next_option(plint) ;
idc = PLUTO_get_idcode(plint) ; /* get dataset item */
old_dset = PLUTO_find_dset(idc) ; /* get ptr to dataset */
if( old_dset == NULL )
return "*************************\n"
"Cannot find Input Dataset\n"
"*************************" ;
/*--------- go to second input line ---------*/
PLUTO_next_option(plint) ;
new_prefix = PLUTO_get_string(plint) ; /* get string item (the output prefix) */
if( ! PLUTO_prefix_ok(new_prefix) ) /* check if it is OK */
return "************************\n"
"Output Prefix is illegal\n"
"************************" ;
str = PLUTO_get_string(plint) ; /* get string item (the datum type) */
istr = PLUTO_string_index( str , /* find it in the list it came from */
NUM_TYPE_STRINGS ,
type_strings ) ;
switch( istr ){
default:
case 0:
new_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* use old dataset type */
break ;
case 1: new_datum = MRI_byte ; break ; /* assign type of user's choice */
case 2: new_datum = MRI_short ; break ;
case 3: new_datum = MRI_float ; break ;
}
/*--------- go to next input lines ---------*/
PLUTO_next_option(plint) ; /* skip to next line */
ignore = PLUTO_get_number(plint) ; /* get number item (ignore) */
PLUTO_next_option(plint) ; /* skip to next line */
taper = PLUTO_get_number(plint) * 0.01 ; /* get number item (taper %) */
/* compute FFT length to use */
PLUTO_next_option(plint) ; /* skip to next line */
str = PLUTO_get_string(plint) ; /* get string item for FFT count */
ninp = DSET_NUM_TIMES(old_dset) ; /* number of values in input */
nuse = ninp - ignore ; /* number of values to actually use */
if( nuse < 4 )
return "*****************************\n"
"Not enough time points to FFT\n"
"*****************************" ;
if( strcmp(str,fft_strings[0]) == 0 ){
/*-- get next larger power-of-2 --*/
#if 0
for( nfft=32 ; nfft < nuse ; nfft *= 2 ) ; /* loop until nfft >= nuse */
#else
nfft = csfft_nextup_one35(nuse) ;
#endif
} else {
nfft = strtol( str , NULL , 10 ) ; /* just convert string to integer */
}
/* if the input FFT length is less than the data length,
tell the user and truncate the amount of data to use */
if( nfft < nuse ){
char str[256] ;
sprintf( str , "******************************\n"
"Warning:\n"
" Number of points in FFT =%4d\n"
" is less than available data\n"
" in time series = %d\n"
"******************************" ,
nfft , nuse ) ;
PLUTO_popup_transient( plint , str ) ;
nuse = nfft ; /* can't use more data than the FFT length */
}
/* compute the number of output points and the output grid spacing */
nfreq = nfft / 2 ; /* # frequencies */
dfreq = 1.0 / (nfft * DSET_TIMESTEP(old_dset) ) ; /* frequency grid */
switch( DSET_TIMEUNITS(old_dset) ){
case UNITS_MSEC_TYPE: dfreq *= 1000.0 ; new_units = UNITS_HZ_TYPE ; break;
case UNITS_SEC_TYPE: new_units = UNITS_HZ_TYPE ; break;
case UNITS_HZ_TYPE: new_units = UNITS_SEC_TYPE; break;
default: new_units = DSET_TIMEUNITS(old_dset) ; break ; /* shouldn't happen */
}
/*------------------------------------------------------*/
/*---------- At this point, the inputs are OK ----------*/
PLUTO_popup_meter( plint ) ; /* popup a progress meter */
/*--------- set up pointers to each sub-brick in the input dataset ---------*/
DSET_load( old_dset ) ; /* must be in memory before we get pointers to it */
old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* get old dataset datum type */
switch( old_datum ){ /* pointer type depends on input datum type */
default:
return "******************************\n"
"Illegal datum in Input Dataset\n"
"******************************" ;
/** create array of pointers into old dataset sub-bricks **/
/** Note that we skip the first 'ignore' sub-bricks here **/
/*--------- input is bytes ----------*/
/* voxel #i at time #k is bptr[k][i] */
/* for i=0..nvox-1 and k=0..nuse-1. */
case MRI_byte:
bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;
if( bptr == NULL ) return "Malloc\nFailure!\n [bptr]" ;
for( kk=0 ; kk < nuse ; kk++ )
bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk+ignore) ;
break ;
/*--------- input is shorts ---------*/
/* voxel #i at time #k is sptr[k][i] */
/* for i=0..nvox-1 and k=0..nuse-1. */
case MRI_short:
sptr = (short **) malloc( sizeof(short *) * nuse ) ;
if( sptr == NULL ) return "Malloc\nFailure!\n [sptr]" ;
for( kk=0 ; kk < nuse ; kk++ )
sptr[kk] = (short *) DSET_ARRAY(old_dset,kk+ignore) ;
break ;
/*--------- input is floats ---------*/
/* voxel #i at time #k is fptr[k][i] */
/* for i=0..nvox-1 and k=0..nuse-1. */
case MRI_float:
fptr = (float **) malloc( sizeof(float *) * nuse ) ;
if( fptr == NULL ) return "Malloc\nFailure!\n [fptr]" ;
for( kk=0 ; kk < nuse ; kk++ )
fptr[kk] = (float *) DSET_ARRAY(old_dset,kk+ignore) ;
break ;
} /* end of switch on input type */
/*---- allocate space for 2 voxel timeseries and 1 FFT ----*/
cxar = (complex *) malloc( sizeof(complex) * nfft ) ; /* FFT */
fxar = (float *) malloc( sizeof(float) * nuse ) ; /* input */
fyar = (float *) malloc( sizeof(float) * nuse ) ; /* input */
if( cxar == NULL || fxar == NULL || fyar == NULL ){
FREE_WORKSPACE ;
return "Malloc\nFailure!\n [cxar]" ;
}
/*--------- make space for taper coefficient array ---------*/
tar = (float *) malloc( sizeof(float) * MAX(nuse,nfreq) ) ;
dtr = (float *) malloc( sizeof(float) * nuse ) ;
if( tar == NULL || dtr == NULL ){
FREE_WORKSPACE ;
return "Malloc\nFailure!\n [tar]" ;
}
ntaper = (int)(0.5 * taper * nuse + 0.49) ; /* will taper data over */
phi = PI / MAX(ntaper,1) ; /* kk=0..ntaper-1 on left */
ktbot = nuse - ntaper ; /* kk=ktbot..nuse-1 on right */
pfact = 0.0 ; /* sum of taper**2 */
for( kk=0 ; kk < nuse ; kk++ ){ /* Hamming-ize */
if( kk < ntaper )
tar[kk] = 0.54 - 0.46 * cos(kk*phi) ; /* ramp up */
else if( kk >= ktbot )
tar[kk] = 0.54 + 0.46 * cos((kk-ktbot+1)*phi) ; /* ramp down */
else
tar[kk] = 1.0 ; /* in the middle */
pfact += tar[kk] * tar[kk] ;
dtr[kk] = kk - 0.5 * (nuse-1) ; /* factors for linear detrending */
}
d0fac = 1.0 / nuse ;
d1fac = 12.0 / nuse / (nuse*nuse - 1.0) ;
/*--- compute factor to go from |FFT|**2 to PSD;
includes the scaling needed for loss of energy with tapering ---*/
pfact = DSET_TIMESTEP(old_dset) / pfact ;
/*--- include scaling factors for sub-bricks, if any ---*/
for( kk=0 ; kk < nuse ; kk++ )
if( DSET_BRICK_FACTOR(old_dset,kk+ignore) > 0.0 )
tar[kk] *= DSET_BRICK_FACTOR(old_dset,kk+ignore) ;
/*---------------------- make a new dataset ----------------------*/
new_dset = EDIT_empty_copy( old_dset ) ; /* start with copy of old one */
{ char * his = PLUTO_commandstring(plint) ;
tross_Copy_History( old_dset , new_dset ) ;
tross_Append_History( new_dset , his ) ; free(his) ;
}
/*-- edit some of its internal parameters --*/
ii = EDIT_dset_items(
new_dset ,
ADN_prefix , new_prefix , /* filename prefix */
ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
ADN_datum_all , new_datum , /* atomic datum */
ADN_nvals , nfreq , /* # sub-bricks */
ADN_ntt , nfreq , /* # time points */
ADN_ttorg , dfreq , /* time origin */
ADN_ttdel , dfreq , /* time step */
ADN_ttdur , dfreq , /* time duration */
ADN_nsl , 0 , /* z-axis time slicing */
ADN_tunits , new_units , /* time units */
ADN_none ) ;
if( ii != 0 ){
THD_delete_3dim_dataset( new_dset , False ) ;
FREE_WORKSPACE ;
return "***********************************\n"
"Error while creating output dataset\n"
"***********************************" ;
}
/*------ make floating point output sub-bricks
(only at the end will scale to byte or shorts)
Output #ii at freq #kk will go into fout[kk][ii],
for kk=0..nfreq-1, and for ii=0..nvox-1. ------*/
nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;
fout = (float **) malloc( sizeof(float *) * nfreq ) ; /* ptrs to sub-bricks */
if( fout == NULL ){
THD_delete_3dim_dataset( new_dset , False ) ;
FREE_WORKSPACE ;
return "Malloc\nFailure!\n [fout]" ;
}
for( kk=0 ; kk < nfreq ; kk++ ){
fout[kk] = (float *) malloc( sizeof(float) * nvox ) ; /* sub-brick # kk */
if( fout[kk] == NULL ) break ;
}
if( kk < nfreq ){
for( ; kk >= 0 ; kk-- ) FREEUP(fout[kk]) ; /* free all we did get */
THD_delete_3dim_dataset( new_dset , False ) ;
FREE_WORKSPACE ;
return "Malloc\nFailure!\n [arrays]" ;
}
{ char buf[128] ;
ii = (nfreq * nvox * sizeof(float)) / (1024*1024) ;
sprintf( buf , " \n"
"*** 3D+time Power Spectral Density:\n"
"*** Using %d MBytes of workspace,\n "
"*** with FFT length = %d\n" , ii,nfft ) ;
PLUTO_popup_transient( plint , buf ) ;
}
/*----------------------------------------------------*/
/*----- Setup has ended. Now do some real work. -----*/
/***** loop over voxels *****/
for( ii=0 ; ii < nvox ; ii += 2 ){ /* 2 time series at a time */
iip = (ii+1) % nvox ; /* voxel after ii */
/*** load data from input dataset, depending on type ***/
switch( old_datum ){
/*** input = bytes ***/
case MRI_byte:
for( kk=0 ; kk < nuse ; kk++ ){
fxar[kk] = bptr[kk][ii] ;
fyar[kk] = bptr[kk][iip] ;
}
break ;
/*** input = shorts ***/
case MRI_short:
for( kk=0 ; kk < nuse ; kk++ ){
fxar[kk] = sptr[kk][ii] ;
fyar[kk] = sptr[kk][iip] ;
}
break ;
/*** input = floats ***/
case MRI_float:
for( kk=0 ; kk < nuse ; kk++ ){
fxar[kk] = fptr[kk][ii] ;
fyar[kk] = fptr[kk][iip] ;
}
break ;
} /* end of switch over input type */
/*** detrend:
x0 = sum( fxar[kk] )
x1 = sum( fxar[kk] * (kk-0.5*(N-1)) )
x0 is used to remove the mean of fxar
x1 is used to remove the linear trend of fxar ***/
x0 = x1 = y0 = y1 = 0.0 ;
for( kk=0 ; kk < nuse ; kk++ ){
x0 += fxar[kk] ; x1 += fxar[kk] * dtr[kk] ;
y0 += fyar[kk] ; y1 += fyar[kk] * dtr[kk] ;
}
x0 *= d0fac ; x1 *= d1fac ; /* factors to remove mean and trend */
y0 *= d0fac ; y1 *= d1fac ;
for( kk=0 ; kk < nuse ; kk++ ){
fxar[kk] -= (x0 + x1 * dtr[kk]) ; /* remove mean and trend here! */
fyar[kk] -= (y0 + y1 * dtr[kk]) ;
}
/*** taper, scale, and put into cxar array ***/
for( kk=0 ; kk < nuse ; kk++ ){
cxar[kk].r = fxar[kk] * tar[kk] ;
cxar[kk].i = fyar[kk] * tar[kk] ;
}
/*** load zeros after where data was put ***/
for( kk=nuse ; kk < nfft ; kk++ ) cxar[kk].r = cxar[kk].i = 0.0 ;
/***** do the FFT (at long last) *****/
csfft_cox( -1 , nfft , cxar ) ;
/***** now compute output into corresponding voxels in fout *****/
/*--- Let x = fxar (1st real time series)
y = fyar (2nd real time series)
z = cxar (complex time series) = x + i y
N = nfft (length of FFT)
Then after FFT, since x and y are real, we have
zhat[k] = xhat[k] + i yhat[k] > for k=1..N/2
zhat[N-k]* = xhat[k] - i yhat[k]
so we can untangle the FFTs of x and y by
xhat[k] = 0.5 ( zhat[k] + zhat[N-k]* )
yhat[k] = 0.5 ( zhat[k] - zhat[N-k]* ) / i
This is the basis for doing 2 time series at once. ---*/
for( kk=1 ; kk <= nfreq ; kk++ ){
xr = 0.5 * ( cxar[kk].r + cxar[nfft-kk].r ) ; /* Re xhat[kk] */
xi = 0.5 * ( cxar[kk].i - cxar[nfft-kk].i ) ; /* Im xhat[kk] */
yr = 0.5 * ( cxar[kk].i + cxar[nfft-kk].i ) ; /* Re yhat[kk] */
yi = 0.5 * ( cxar[kk].r - cxar[nfft-kk].r ) ; /*-Im yhat[kk] */
fout[kk-1][ii] = pfact * (xr*xr + xi*xi) ;
fout[kk-1][iip] = pfact * (yr*yr + yi*yi) ;
}
perc = (100 * ii) / nvox ; /* display percentage done */
PLUTO_set_meter( plint , perc ) ; /* on the progress meter */
} /* end of outer loop over 2 voxels at a time */
DSET_unload( old_dset ) ; /* don't need this no more */
/*------------------------------------------------------------*/
/*------- The output is now in fout[kk][ii],
for kk=0..nfreq-1 , ii=0..nvox-1.
We must now put this into the output dataset -------*/
switch( new_datum ){
/*** output is floats is the simplest:
we just have to attach the fout bricks to the dataset ***/
case MRI_float:
for( kk=0 ; kk < nfreq ; kk++ )
EDIT_substitute_brick( new_dset , kk , MRI_float , fout[kk] ) ;
break ;
/*** output is shorts:
we have to create a scaled sub-brick from fout ***/
case MRI_short:{
short * bout ;
float fac ;
for( kk=0 ; kk < nfreq ; kk++ ){ /* loop over sub-bricks */
/*-- get output sub-brick --*/
bout = (short *) malloc( sizeof(short) * nvox ) ;
if( bout == NULL ){
fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
return("\nFinal malloc error in plug_power!\n") ;
/* EXIT(1) ;*/
}
/*-- find scaling and then scale --*/
fac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[kk] ) ;
if( fac > 0.0 ){
fac = 32767.0 / fac ;
EDIT_coerce_scale_type( nvox,fac ,
MRI_float,fout[kk] , MRI_short,bout ) ;
fac = 1.0 / fac ;
}
free( fout[kk] ) ; /* don't need this anymore */
/*-- put output brick into dataset, and store scale factor --*/
EDIT_substitute_brick( new_dset , kk , MRI_short , bout ) ;
tar[kk] = fac ;
perc = (100 * kk) / nfreq ;
PLUTO_set_meter( plint , perc ) ; /* on the progress meter */
}
/*-- save scale factor array into dataset --*/
EDIT_dset_items( new_dset , ADN_brick_fac , tar , ADN_none ) ;
}
break ;
/*** output is bytes (byte = unsigned char)
we have to create a scaled sub-brick from fout ***/
case MRI_byte:{
byte * bout ;
float fac ;
for( kk=0 ; kk < nfreq ; kk++ ){ /* loop over sub-bricks */
/*-- get output sub-brick --*/
bout = (byte *) malloc( sizeof(byte) * nvox ) ;
if( bout == NULL ){
fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
return("\nFinal malloc error in plug_power!\n\a") ;
/* EXIT(1) ;*/
}
/*-- find scaling and then scale --*/
fac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[kk] ) ;
if( fac > 0.0 ){
fac = 255.0 / fac ;
EDIT_coerce_scale_type( nvox,fac ,
MRI_float,fout[kk] , MRI_byte,bout ) ;
fac = 1.0 / fac ;
}
free( fout[kk] ) ; /* don't need this anymore */
/*-- put output brick into dataset, and store scale factor --*/
EDIT_substitute_brick( new_dset , kk , MRI_byte , bout ) ;
tar[kk] = fac ;
perc = (100 * kk) / nfreq ;
PLUTO_set_meter( plint , perc ) ; /* on the progress meter */
}
/*-- save scale factor array into dataset --*/
EDIT_dset_items( new_dset , ADN_brick_fac , tar , ADN_none ) ;
}
break ;
} /* end of switch on output data type */
/*-------------- Cleanup and go home ----------------*/
PLUTO_set_meter( plint , 100 ) ; /* set progress meter to 100% */
PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
FREE_WORKSPACE ;
return NULL ; /* null string returned means all was OK */
}
#ifdef ALLOW_TESTING
/*****************************************************************************
-----------------------------------------------------------------------------
Create the second interface within this plugin.
-----------------------------------------------------------------------------*/
PLUGIN_interface * TEST_init( void )
{
PLUGIN_interface * plint ; /* will be the output of this routine */
/*---------------- set titles and call point ----------------*/
plint = PLUTO_new_interface( "Testing" ,
"Testing, Testing, 1-2-3 ..." ,
NULL ,
PLUGIN_CALL_VIA_MENU , TEST_main ) ;
PLUTO_add_hint( plint , "1-2-3, 1-2-3, ..." ) ;
PLUTO_add_option( plint ,
"Input" , /* label at left of input line */
"Input" , /* tag to return to plugin */
TRUE /* is this mandatory? */
) ;
PLUTO_add_dataset_list( plint ,
"Datasets" , /* label next to button */
ANAT_ALL_MASK , /* take any anat datasets */
FUNC_FIM_MASK , /* only allow fim funcs */
DIMEN_4D_MASK | /* need 3D+time datasets */
BRICK_ALLREAL_MASK /* need real-valued datasets */
) ;
return plint ;
}
char * TEST_main( PLUGIN_interface * plint )
{
MRI_IMAGE * tsim ;
MCW_idclist * idclist ;
MCW_idcode * idc ;
THD_3dim_dataset * dset ;
char str[256] ;
int id ;
/*--------- go to first input line ---------*/
PLUTO_next_option(plint) ;
idclist = PLUTO_get_idclist(plint) ;
if( PLUTO_idclist_count(idclist) == 0 )
return " \nNo input dataset list!\n " ;
id = 0 ;
do {
idc = PLUTO_idclist_next(idclist) ;
dset = PLUTO_find_dset(idc) ;
if( dset == NULL ) return NULL ;
id++ ;
sprintf(str, " \nDataset %d = %s\n nx = %d\n ny = %d\n nz = %d\n " ,
id , DSET_FILECODE(dset) , dset->daxes->nxx,dset->daxes->nyy,dset->daxes->nzz ) ;
PLUTO_popup_transient( plint , str ) ;
} while(1) ;
return NULL ;
}
#endif /* ALLOW_TESTING */
syntax highlighted by Code2HTML, v. 0.9.1