Doxygen Source Code Documentation
3dWavelets.c File Reference
#include "mrilib.h"
#include "Wavelets.h"
#include "Wavelets.c"
Go to the source code of this file.
Data Structures | |
struct | WA_options |
Defines | |
#define | PROGRAM_NAME "3dWavelets" |
#define | PROGRAM_AUTHOR "B. Douglas Ward" |
#define | PROGRAM_INITIAL "28 March 2000" |
#define | PROGRAM_LATEST "02 December 2002" |
#define | MAX_NAME_LENGTH THD_MAX_NAME |
#define | MAX_FILTERS 20 |
Typedefs | |
typedef WA_options | WA_options |
Functions | |
void | display_help_menu () |
void | initialize_options (WA_options *option_data) |
void | get_options (int argc, char **argv, WA_options *option_data) |
float * | read_time_series (char *ts_filename, int *ts_length) |
void | read_input_data (WA_options *option_data, THD_3dim_dataset **dset_time, THD_3dim_dataset **mask_dset, float **fmri_data, int *fmri_length) |
void | check_one_output_file (THD_3dim_dataset *dset_time, char *filename) |
void | check_output_files (WA_options *option_data, THD_3dim_dataset *dset_time) |
void | check_for_valid_inputs (WA_options *option_data, THD_3dim_dataset *dset_time, THD_3dim_dataset *mask_dset) |
void | initialize_filters (WA_options *option_data, THD_3dim_dataset *dset_time, THD_3dim_dataset *mask_dset, int fmri_length) |
void | allocate_memory (WA_options *option_data, THD_3dim_dataset *dset, float ***coef_vol, float **mse_vol, float **ffull_vol, float **rfull_vol, float ***coefts_vol, float ***fitts_vol, float ***sgnlts_vol, float ***errts_vol) |
void | initialize_program (int argc, char **argv, WA_options **option_data, THD_3dim_dataset **dset_time, THD_3dim_dataset **mask_dset, float **fmri_data, int *fmri_length, float ***coef_vol, float **mse_vol, float **ffull_vol, float **rfull_vol, float ***coefts_vol, float ***fitts_vol, float ***sgnlts_vol, float ***errts_vol) |
void | extract_ts_array (THD_3dim_dataset *dset_time, int iv, float *ts_array) |
void | save_voxel (WA_options *option_data, int iv, float *coef, float mse, float ffull, float rfull, float *coefts, float *fitts, float *sgnlts, float *errts, float **coef_vol, float *mse_vol, float *ffull_vol, float *rfull_vol, float **coefts_vol, float **fitts_vol, float **sgnlts_vol, float **errts_vol) |
void | calculate_results (WA_options *option_data, THD_3dim_dataset *dset, THD_3dim_dataset *mask, float *fmri_data, int fmri_length, float **coef_vol, float *mse_vol, float *ffull_vol, float *rfull_vol, float **coefts_vol, float **fitts_vol, float **sgnlts_vol, float **errts_vol) |
float | EDIT_coerce_autoscale_new (int nxyz, int itype, void *ivol, int otype, void *ovol) |
void | write_ts_array (int argc, char **argv, WA_options *option_data, int ts_length, float **vol_array, char *output_filename) |
void | attach_sub_brick (THD_3dim_dataset *new_dset, int ibrick, float *volume, int nxyz, int brick_type, char *brick_label, int dof, int ndof, int ddof, short **bar) |
void | write_bucket_data (int argc, char **argv, WA_options *option_data, float **coef_vol, float *mse_vol, float *ffull_vol, float *rfull_vol) |
void | output_results (int argc, char **argv, WA_options *option_data, float **coef_vol, float *mse_vol, float *ffull_vol, float *rfull_vol, float **coefts_vol, float **fitts_vol, float **sgnlts_vol, float **errts_vol) |
void | terminate_program (WA_options **option_data, float **fmri_data, float ***coef_vol, float **mse_vol, float **ffull_vol, float **rfull_vol, float ***coefts_vol, float ***fitts_vol, float ***sgnlts_vol, float ***errts_vol) |
int | main (int argc, char **argv) |
Define Documentation
|
Definition at line 48 of file 3dWavelets.c. Referenced by initialize_options(). |
|
Definition at line 47 of file 3dWavelets.c. Referenced by check_for_valid_inputs(), check_one_output_file(), get_options(), initialize_filters(), read_input_data(), read_time_series(), and write_bucket_data(). |
|
Definition at line 28 of file 3dWavelets.c. Referenced by main(). |
|
Definition at line 29 of file 3dWavelets.c. Referenced by main(). |
|
Definition at line 30 of file 3dWavelets.c. Referenced by main(). |
|
Definition at line 27 of file 3dWavelets.c. Referenced by get_options(), main(), write_bucket_data(), and write_ts_array(). |
Typedef Documentation
|
|
Function Documentation
|
Definition at line 948 of file 3dWavelets.c. References malloc, MTEST, and p. Referenced by initialize_program().
00963 { 00964 int nxyz; /* total number of voxels */ 00965 int ixyz; /* voxel index */ 00966 int it; /* time point index */ 00967 int N; /* number of usable data points */ 00968 int ip; /* parameter index */ 00969 int p; /* number of parameters in the full model */ 00970 00971 00972 /*----- Initialize local variables -----*/ 00973 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz; 00974 N = option_data->N; 00975 p = option_data->p; 00976 00977 00978 /*----- Allocate memory space for volume data -----*/ 00979 *coef_vol = (float **) malloc (sizeof(float *) * p); MTEST(*coef_vol); 00980 00981 for (ip = 0; ip < p; ip++) 00982 { 00983 (*coef_vol)[ip] = (float *) malloc (sizeof(float) * nxyz); 00984 MTEST((*coef_vol)[ip]); 00985 for (ixyz = 0; ixyz < nxyz; ixyz++) 00986 { 00987 (*coef_vol)[ip][ixyz] = 0.0; 00988 } 00989 } 00990 00991 00992 *mse_vol = (float *) malloc (sizeof(float) * nxyz); MTEST(*mse_vol); 00993 *ffull_vol = (float *) malloc (sizeof(float) * nxyz); MTEST(*ffull_vol); 00994 *rfull_vol = (float *) malloc (sizeof(float) * nxyz); MTEST(*rfull_vol); 00995 for (ixyz = 0; ixyz < nxyz; ixyz++) 00996 { 00997 (*mse_vol)[ixyz] = 0.0; 00998 (*ffull_vol)[ixyz] = 0.0; 00999 (*rfull_vol)[ixyz] = 0.0; 01000 } 01001 01002 01003 /*----- Allocate memory for forward wavelet transform coefficients -----*/ 01004 if (option_data->coefts_filename != NULL) 01005 { 01006 *coefts_vol = (float **) malloc (sizeof(float **) * N); 01007 MTEST (*coefts_vol); 01008 for (it = 0; it < N; it++) 01009 { 01010 (*coefts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz); 01011 MTEST ((*coefts_vol)[it]); 01012 for (ixyz = 0; ixyz < nxyz; ixyz++) 01013 (*coefts_vol)[it][ixyz] = 0.0; 01014 } 01015 } 01016 01017 /*----- Allocate memory for filtered time series -----*/ 01018 if (option_data->fitts_filename != NULL) 01019 { 01020 *fitts_vol = (float **) malloc (sizeof(float **) * N); 01021 MTEST (*fitts_vol); 01022 for (it = 0; it < N; it++) 01023 { 01024 (*fitts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz); 01025 MTEST ((*fitts_vol)[it]); 01026 for (ixyz = 0; ixyz < nxyz; ixyz++) 01027 (*fitts_vol)[it][ixyz] = 0.0; 01028 } 01029 } 01030 01031 /*----- Allocate memory for signal model time series -----*/ 01032 if (option_data->sgnlts_filename != NULL) 01033 { 01034 *sgnlts_vol = (float **) malloc (sizeof(float **) * N); 01035 MTEST (*sgnlts_vol); 01036 for (it = 0; it < N; it++) 01037 { 01038 (*sgnlts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz); 01039 MTEST ((*sgnlts_vol)[it]); 01040 for (ixyz = 0; ixyz < nxyz; ixyz++) 01041 (*sgnlts_vol)[it][ixyz] = 0.0; 01042 } 01043 } 01044 01045 /*----- Allocate memory for residual time series -----*/ 01046 if (option_data->errts_filename != NULL) 01047 { 01048 *errts_vol = (float **) malloc (sizeof(float **) * N); 01049 MTEST (*errts_vol); 01050 for (it = 0; it < N; it++) 01051 { 01052 (*errts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz); 01053 MTEST ((*errts_vol)[it]); 01054 for (ixyz = 0; ixyz < nxyz; ixyz++) 01055 (*errts_vol)[it][ixyz] = 0.0; 01056 } 01057 } 01058 } |
|
Definition at line 1570 of file 3dWavelets.c. References EDIT_BRICK_FACTOR, EDIT_BRICK_LABEL, EDIT_BRICK_TO_FIFT, EDIT_BRICK_TO_FITT, EDIT_coerce_autoscale_new(), EDIT_substitute_brick(), malloc, and MTEST. Referenced by do_xrestore_stuff(), and write_bucket_data().
01583 { 01584 const float EPSILON = 1.0e-10; 01585 float factor; /* factor is new scale factor for sub-brick #ib */ 01586 01587 01588 /*----- allocate memory for output sub-brick -----*/ 01589 bar[ibrick] = (short *) malloc (sizeof(short) * nxyz); 01590 MTEST (bar[ibrick]); 01591 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume, 01592 MRI_short, bar[ibrick]); 01593 01594 if (factor < EPSILON) factor = 0.0; 01595 else factor = 1.0 / factor; 01596 01597 01598 /*----- edit the sub-brick -----*/ 01599 EDIT_BRICK_LABEL (new_dset, ibrick, brick_label); 01600 EDIT_BRICK_FACTOR (new_dset, ibrick, factor); 01601 01602 if (brick_type == FUNC_TT_TYPE) 01603 EDIT_BRICK_TO_FITT (new_dset, ibrick, dof); 01604 else if (brick_type == FUNC_FT_TYPE) 01605 EDIT_BRICK_TO_FIFT (new_dset, ibrick, ndof, ddof); 01606 01607 01608 /*----- attach bar[ib] to be sub-brick #ibrick -----*/ 01609 EDIT_substitute_brick (new_dset, ibrick, MRI_short, bar[ibrick]); 01610 01611 } |
|
Definition at line 1255 of file 3dWavelets.c. References DSET_NUM_TIMES, extract_ts_array(), free, i, malloc, MTEST, p, q, report_results(), save_voxel(), ts_fprint(), and wavelet_analysis().
01273 { 01274 float * ts_array = NULL; /* array of measured data for one voxel */ 01275 float mask_val[1]; /* value of mask at current voxel */ 01276 01277 int f; 01278 int p; /* number of parameters in the full model */ 01279 int q; /* number of parameters in the baseline model */ 01280 int m; /* parameter index */ 01281 int n; /* data point index */ 01282 01283 float * coef = NULL; /* regression parameters */ 01284 float sse_base; /* baseline model error sum of squares */ 01285 float sse_full; /* full model error sum of squares */ 01286 float ffull; /* full model F-statistic */ 01287 float rfull; /* full model R^2 stat. */ 01288 01289 int ixyz; /* voxel index */ 01290 int nxyz; /* number of voxels per dataset */ 01291 01292 int nt; /* number of images in input 3d+time dataset */ 01293 int NFirst; /* first image from input 3d+time dataset to use */ 01294 int NLast; /* last image from input 3d+time dataset to use */ 01295 int N; /* number of usable data points */ 01296 01297 int i; /* data point index */ 01298 01299 float * coefts = NULL; /* filtered FWT coefficients */ 01300 float * fitts = NULL; /* filterd time series */ 01301 float * sgnlts = NULL; /* signal model fitted time series */ 01302 float * errts = NULL; /* residual error time series */ 01303 01304 char * label = NULL; /* string containing stat. summary of results */ 01305 01306 01307 /*----- Initialize local variables -----*/ 01308 if (option_data->input1D_filename != NULL) 01309 { 01310 nxyz = 1; 01311 nt = fmri_length; 01312 } 01313 else 01314 { 01315 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz; 01316 nt = DSET_NUM_TIMES (dset); 01317 } 01318 01319 NFirst = option_data->NFirst; 01320 NLast = option_data->NLast; 01321 N = option_data->N; 01322 f = option_data->f; 01323 q = option_data->q; 01324 p = option_data->p; 01325 01326 01327 /*----- Allocate memory for arrays -----*/ 01328 ts_array = (float *) malloc (sizeof(float) * nt); MTEST (ts_array); 01329 coef = (float *) malloc (sizeof(float) * p); MTEST (coef); 01330 coefts = (float *) malloc (sizeof(float) * N); MTEST (coefts); 01331 fitts = (float *) malloc (sizeof(float) * N); MTEST (fitts); 01332 sgnlts = (float *) malloc (sizeof(float) * N); MTEST (sgnlts); 01333 errts = (float *) malloc (sizeof(float) * N); MTEST (errts); 01334 01335 01336 /*----- Loop over all voxels -----*/ 01337 for (ixyz = 0; ixyz < nxyz; ixyz++) 01338 { 01339 /*----- Apply mask? -----*/ 01340 if (mask != NULL) 01341 { 01342 extract_ts_array (mask, ixyz, mask_val); 01343 if (mask_val[0] == 0.0) continue; 01344 } 01345 01346 01347 /*----- Extract Y-data for this voxel -----*/ 01348 if (option_data->input1D_filename != NULL) 01349 { 01350 for (i = 0; i < N; i++) 01351 ts_array[i] = fmri_data[i+NFirst]; 01352 } 01353 else 01354 extract_ts_array (dset, ixyz, ts_array); 01355 01356 01357 /*----- Perform the wavelet analysis for this voxel-----*/ 01358 wavelet_analysis (option_data->wavelet_type, 01359 f, option_data->stop_filter, 01360 q, option_data->base_filter, 01361 p, option_data->sgnl_filter, 01362 N, ts_array+NFirst, coef, 01363 &sse_base, &sse_full, &ffull, &rfull, 01364 coefts, fitts, sgnlts, errts); 01365 01366 01367 /*----- Save results for this voxel -----*/ 01368 if (option_data->input1D_filename == NULL) 01369 save_voxel (option_data, ixyz, coef, sse_full/(N-f-p), ffull, rfull, 01370 coefts, fitts, sgnlts, errts, 01371 coef_vol, mse_vol, ffull_vol, rfull_vol, 01372 coefts_vol, fitts_vol, sgnlts_vol, errts_vol); 01373 else 01374 { 01375 /*----- If only one voxel, save results as .1D files -----*/ 01376 ts_fprint ("WA.coefts.1D", N, coefts); 01377 ts_fprint ("WA.fitts.1D", N, fitts); 01378 ts_fprint ("WA.sgnlts.1D", N, sgnlts); 01379 ts_fprint ("WA.errts.1D", N, errts); 01380 } 01381 01382 01383 /*----- Report results for this voxel -----*/ 01384 if ( ((ffull > option_data->fdisp) && (option_data->fdisp >= 0.0)) 01385 || (option_data->input1D_filename != NULL) ) 01386 { 01387 printf ("\n\nResults for Voxel #%d: \n", ixyz); 01388 report_results (N, NFirst, f, p, q, 01389 option_data->base_filter, option_data->sgnl_filter, 01390 coef, sse_base, sse_full, ffull, rfull, 01391 &label); 01392 printf ("%s \n", label); 01393 } 01394 01395 01396 } /*----- Loop over voxels -----*/ 01397 01398 01399 01400 /*----- Release memory -----*/ 01401 free (ts_array); ts_array = NULL; 01402 free (coef); coef = NULL; 01403 free (coefts); coefts = NULL; 01404 free (fitts); fitts = NULL; 01405 free (sgnlts); sgnlts = NULL; 01406 free (errts); errts = NULL; 01407 01408 } |
|
Definition at line 814 of file 3dWavelets.c. References check_output_files(), DSET_NVALS, DSET_NX, DSET_NY, DSET_NZ, MAX_NAME_LENGTH, and WA_error().
00820 { 00821 char message[MAX_NAME_LENGTH]; /* error message */ 00822 00823 00824 /*----- If mask is used, check for compatible dimensions -----*/ 00825 if (mask_dset != NULL) 00826 { 00827 if ( (DSET_NX(dset_time) != DSET_NX(mask_dset)) 00828 || (DSET_NY(dset_time) != DSET_NY(mask_dset)) 00829 || (DSET_NZ(dset_time) != DSET_NZ(mask_dset)) ) 00830 { 00831 sprintf (message, "%s and %s have incompatible dimensions", 00832 option_data->input_filename, option_data->mask_filename); 00833 WA_error (message); 00834 } 00835 00836 if (DSET_NVALS(mask_dset) != 1 ) 00837 WA_error ("Must specify 1 sub-brick from mask dataset"); 00838 } 00839 00840 00841 /*----- Check whether any of the output files already exist -----*/ 00842 check_output_files (option_data, dset_time); 00843 00844 } |
|
Definition at line 731 of file 3dWavelets.c. References ADN_label1, ADN_none, ADN_prefix, ADN_self_name, ADN_type, THD_3dim_dataset::dblk, THD_datablock::diskptr, EDIT_dset_items(), EDIT_empty_copy(), GEN_FUNC_TYPE, HEAD_FUNC_TYPE, THD_diskptr::header_name, ISHEAD, MAX_NAME_LENGTH, THD_delete_3dim_dataset(), THD_is_file(), and WA_error().
00736 { 00737 char message[MAX_NAME_LENGTH]; /* error message */ 00738 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */ 00739 int ierror; /* number of errors in editing data */ 00740 00741 00742 /*----- make an empty copy of input dataset -----*/ 00743 new_dset = EDIT_empty_copy( dset_time ) ; 00744 00745 00746 ierror = EDIT_dset_items( new_dset , 00747 ADN_prefix , filename , 00748 ADN_label1 , filename , 00749 ADN_self_name , filename , 00750 ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE : 00751 GEN_FUNC_TYPE , 00752 ADN_none ) ; 00753 00754 if( ierror > 0 ) 00755 { 00756 sprintf (message, 00757 "*** %d errors in attempting to create output dataset!\n", 00758 ierror); 00759 WA_error (message); 00760 } 00761 00762 if( THD_is_file(new_dset->dblk->diskptr->header_name) ) 00763 { 00764 sprintf (message, 00765 "Output dataset file %s already exists " 00766 " -- cannot continue!\a\n", 00767 new_dset->dblk->diskptr->header_name); 00768 WA_error (message); 00769 } 00770 00771 /*----- deallocate memory -----*/ 00772 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ; 00773 00774 } |
|
Definition at line 783 of file 3dWavelets.c. References check_one_output_file(). Referenced by check_for_valid_inputs(), initialize(), and initialize_program().
00788 { 00789 00790 if (option_data->bucket_filename != NULL) 00791 check_one_output_file (dset_time, option_data->bucket_filename); 00792 00793 if (option_data->coefts_filename != NULL) 00794 check_one_output_file (dset_time, option_data->coefts_filename); 00795 00796 if (option_data->fitts_filename != NULL) 00797 check_one_output_file (dset_time, option_data->fitts_filename); 00798 00799 if (option_data->sgnlts_filename != NULL) 00800 check_one_output_file (dset_time, option_data->sgnlts_filename); 00801 00802 if (option_data->errts_filename != NULL) 00803 check_one_output_file (dset_time, option_data->errts_filename); 00804 00805 } |
|
Definition at line 110 of file 3dWavelets.c.
00111 { 00112 printf ( 00113 "Program to perform wavelet analysis of an FMRI 3d+time dataset. \n" 00114 " \n" 00115 "Usage: \n" 00116 "3dWavelets \n" 00117 "-type wname wname = name of wavelet to use for the analysis \n" 00118 " At present, there are only two choices for wname: \n" 00119 " Haar --> Haar wavelets \n" 00120 " Daub --> Daubechies wavelets \n" 00121 "-input fname fname = filename of 3d+time input dataset \n" 00122 "[-input1D dname] dname = filename of single (fMRI) .1D time series \n" 00123 "[-mask mname] mname = filename of 3d mask dataset \n" 00124 "[-nfirst fnum] fnum = number of first dataset image to use in \n" 00125 " the wavelet analysis. (default = 0) \n" 00126 "[-nlast lnum] lnum = number of last dataset image to use in \n" 00127 " the wavelet analysis. (default = last) \n" 00128 "[-fdisp fval] Write (to screen) results for those voxels \n" 00129 " whose F-statistic is >= fval \n" 00130 " \n" 00131 "Filter options: \n" 00132 "[-filt_stop band mintr maxtr] Specify wavelet coefs. to set to zero \n" 00133 "[-filt_base band mintr maxtr] Specify wavelet coefs. for baseline model\n" 00134 "[-filt_sgnl band mintr maxtr] Specify wavelet coefs. for signal model \n" 00135 " where band = frequency band \n" 00136 " mintr = min. value for time window (in TR) \n" 00137 " maxtr = max. value for time window (in TR) \n" 00138 " \n" 00139 "Output options: \n" 00140 "[-coefts cprefix] cprefix = prefix of 3d+time output dataset which \n" 00141 " will contain the forward wavelet transform \n" 00142 " coefficients \n" 00143 " \n" 00144 "[-fitts fprefix] fprefix = prefix of 3d+time output dataset which \n" 00145 " will contain the full model time series fit \n" 00146 " to the input data \n" 00147 " \n" 00148 "[-sgnlts sprefix] sprefix = prefix of 3d+time output dataset which \n" 00149 " will contain the signal model time series fit \n" 00150 " to the input data \n" 00151 " \n" 00152 "[-errts eprefix] eprefix = prefix of 3d+time output dataset which \n" 00153 " will contain the residual error time series \n" 00154 " from the full model fit to the input data \n" 00155 " \n" 00156 "The following options control the contents of the bucket dataset: \n" 00157 "[-fout] Flag to output the F-statistics \n" 00158 "[-rout] Flag to output the R^2 statistics \n" 00159 "[-cout] Flag to output the full model wavelet coefficients \n" 00160 "[-vout] Flag to output the sample variance (MSE) map \n" 00161 " \n" 00162 "[-stat_first] Flag to specify that the full model statistics will \n" 00163 " appear prior to the wavelet coefficients in the \n" 00164 " bucket dataset output \n" 00165 " \n" 00166 "[-bucket bprefix] bprefix = prefix of AFNI 'bucket' dataset containing\n" 00167 " parameters of interest, such as the F-statistic \n" 00168 " for significance of the wavelet signal model. \n" 00169 ); 00170 00171 exit(0); 00172 } |
|
compute start time of this timeseries * Definition at line 1422 of file 3dWavelets.c. References EDIT_coerce_scale_type(), MCW_vol_amax(), MRI_IS_INT_TYPE, and top.
01424 { 01425 float fac=0.0 , top ; 01426 01427 if( MRI_IS_INT_TYPE(otype) ){ 01428 top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ; 01429 if (top == 0.0) fac = 0.0; 01430 else fac = MRI_TYPE_maxval[otype]/top; 01431 } 01432 01433 EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ; 01434 return ( fac ); 01435 } |
|
Definition at line 1124 of file 3dWavelets.c. References DSET_NUM_TIMES, FIM_error(), MRI_FLOAT_PTR, mri_free(), THD_extract_series(), and WA_error(). Referenced by calculate_results(), and set_fim_thr_level().
01130 { 01131 MRI_IMAGE * im = NULL; /* intermediate float data */ 01132 float * ar = NULL; /* pointer to float data */ 01133 int ts_length; /* length of input 3d+time data set */ 01134 int it; /* time index */ 01135 01136 01137 /*----- Extract time series from 3d+time data set into MRI_IMAGE -----*/ 01138 im = THD_extract_series (iv, dset_time, 0); 01139 01140 01141 /*----- Verify extraction -----*/ 01142 if (im == NULL) WA_error ("Unable to extract data from 3d+time dataset"); 01143 01144 01145 /*----- Now extract time series from MRI_IMAGE -----*/ 01146 ts_length = DSET_NUM_TIMES (dset_time); 01147 ar = MRI_FLOAT_PTR (im); 01148 for (it = 0; it < ts_length; it++) 01149 { 01150 ts_array[it] = ar[it]; 01151 } 01152 01153 01154 /*----- Release memory -----*/ 01155 mri_free (im); im = NULL; 01156 01157 } |
|
Definition at line 248 of file 3dWavelets.c. References addto_args(), AFNI_logger(), argc, display_help_menu(), initialize_options(), machdep(), malloc, MAX_NAME_LENGTH, MAX_WAVELET_TYPE, MTEST, PROGRAM_NAME, and WA_error().
00254 { 00255 int nopt = 1; /* input option argument counter */ 00256 int ival, index; /* integer input */ 00257 float fval; /* float input */ 00258 char message[MAX_NAME_LENGTH]; /* error message */ 00259 int ifilter; /* filter index */ 00260 00261 00262 /*-- addto the arglist, if user wants to --*/ 00263 machdep() ; 00264 { int new_argc ; char ** new_argv ; 00265 addto_args( argc , argv , &new_argc , &new_argv ) ; 00266 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; } 00267 } 00268 00269 00270 /*----- does user request help menu? -----*/ 00271 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu(); 00272 00273 00274 /*----- Add to program log -----*/ 00275 AFNI_logger (PROGRAM_NAME,argc,argv); 00276 00277 00278 /*----- initialize the input options -----*/ 00279 initialize_options (option_data); 00280 00281 00282 /*----- main loop over input options -----*/ 00283 while (nopt < argc ) 00284 { 00285 00286 /*----- -type wname -----*/ 00287 if (strncmp(argv[nopt], "-type", 5) == 0) 00288 { 00289 nopt++; 00290 if (nopt >= argc) WA_error ("need argument after -type "); 00291 for (ival = 0; ival < MAX_WAVELET_TYPE; ival++) 00292 if (strncmp(argv[nopt], WAVELET_TYPE_name[ival], 4) == 0) break; 00293 if (ival >= MAX_WAVELET_TYPE) 00294 { 00295 sprintf(message,"Unrecognized wavelet type: %s\n", argv[nopt]); 00296 WA_error (message); 00297 } 00298 option_data->wavelet_type = ival; 00299 nopt++; 00300 continue; 00301 } 00302 00303 00304 /*----- -input filename -----*/ 00305 if (strncmp(argv[nopt], "-input", 8) == 0) 00306 { 00307 nopt++; 00308 if (nopt >= argc) WA_error ("need argument after -input "); 00309 option_data->input_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); 00310 MTEST (option_data->input_filename); 00311 strcpy (option_data->input_filename, argv[nopt]); 00312 nopt++; 00313 continue; 00314 } 00315 00316 00317 /*----- -input1D filename -----*/ 00318 if (strncmp(argv[nopt], "-input1D", 8) == 0) 00319 { 00320 nopt++; 00321 if (nopt >= argc) WA_error ("need argument after -input1D "); 00322 option_data->input1D_filename = 00323 malloc (sizeof(char)*MAX_NAME_LENGTH); 00324 MTEST (option_data->input1D_filename); 00325 strcpy (option_data->input1D_filename, argv[nopt]); 00326 nopt++; 00327 continue; 00328 } 00329 00330 00331 /*----- -mask filename -----*/ 00332 if (strncmp(argv[nopt], "-mask", 5) == 0) 00333 { 00334 nopt++; 00335 if (nopt >= argc) WA_error ("need argument after -mask "); 00336 option_data->mask_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); 00337 MTEST (option_data->mask_filename); 00338 strcpy (option_data->mask_filename, argv[nopt]); 00339 nopt++; 00340 continue; 00341 } 00342 00343 00344 /*----- -nfirst num -----*/ 00345 if (strncmp(argv[nopt], "-nfirst", 7) == 0) 00346 { 00347 nopt++; 00348 if (nopt >= argc) WA_error ("need argument after -nfirst "); 00349 sscanf (argv[nopt], "%d", &ival); 00350 if (ival < 0) 00351 WA_error ("illegal argument after -nfirst "); 00352 option_data->NFirst = ival; 00353 nopt++; 00354 continue; 00355 } 00356 00357 00358 /*----- -nlast num -----*/ 00359 if (strncmp(argv[nopt], "-nlast", 6) == 0) 00360 { 00361 nopt++; 00362 if (nopt >= argc) WA_error ("need argument after -nlast "); 00363 sscanf (argv[nopt], "%d", &ival); 00364 if (ival < 0) 00365 WA_error ("illegal argument after -nlast "); 00366 option_data->NLast = ival; 00367 nopt++; 00368 continue; 00369 } 00370 00371 00372 /*----- -fdisp fval -----*/ 00373 if (strncmp(argv[nopt], "-fdisp", 6) == 0) 00374 { 00375 nopt++; 00376 if (nopt >= argc) WA_error ("need argument after -fdisp "); 00377 sscanf (argv[nopt], "%f", &fval); 00378 option_data->fdisp = fval; 00379 nopt++; 00380 continue; 00381 } 00382 00383 00384 /*----- -filt_stop band mintr maxtr -----*/ 00385 if (strncmp(argv[nopt], "-filt_stop", 10) == 0) 00386 { 00387 nopt++; 00388 if (nopt+2 >= argc) 00389 WA_error ("need 3 arguments after -filt_stop"); 00390 00391 ifilter = option_data->num_stop_filters; 00392 00393 sscanf (argv[nopt], "%d", &ival); 00394 if (ival < -1) 00395 WA_error ("-filt_stop band mintr maxtr Require: -1 <= band"); 00396 option_data->stop_band[ifilter] = ival; 00397 nopt++; 00398 00399 sscanf (argv[nopt], "%d", &ival); 00400 if (ival < 0) 00401 WA_error ("-filt_stop band mintr maxtr Require: 0 <= mintr"); 00402 option_data->stop_mintr[ifilter] = ival; 00403 nopt++; 00404 00405 sscanf (argv[nopt], "%d", &ival); 00406 if (ival < 0) 00407 WA_error ("-filt_stop band mintr maxtr Require: 0 <= maxtr"); 00408 option_data->stop_maxtr[ifilter] = ival; 00409 nopt++; 00410 00411 option_data->num_stop_filters++; 00412 continue; 00413 } 00414 00415 00416 /*----- -filt_base band mintr maxtr -----*/ 00417 if (strncmp(argv[nopt], "-filt_base", 10) == 0) 00418 { 00419 nopt++; 00420 if (nopt+2 >= argc) 00421 WA_error ("need 3 arguments after -filt_base"); 00422 00423 ifilter = option_data->num_base_filters; 00424 00425 sscanf (argv[nopt], "%d", &ival); 00426 if (ival < -1) 00427 WA_error ("-filt_base band mintr maxtr Require: -1 <= band"); 00428 option_data->base_band[ifilter] = ival; 00429 nopt++; 00430 00431 sscanf (argv[nopt], "%d", &ival); 00432 if (ival < 0) 00433 WA_error ("-filt_base band mintr maxtr Require: 0 <= mintr"); 00434 option_data->base_mintr[ifilter] = ival; 00435 nopt++; 00436 00437 sscanf (argv[nopt], "%d", &ival); 00438 if (ival < 0) 00439 WA_error ("-filt_base band mintr maxtr Require: 0 <= maxtr"); 00440 option_data->base_maxtr[ifilter] = ival; 00441 nopt++; 00442 00443 option_data->num_base_filters++; 00444 continue; 00445 } 00446 00447 00448 /*----- -filt_sgnl band mintr maxtr -----*/ 00449 if (strncmp(argv[nopt], "-filt_sgnl", 10) == 0) 00450 { 00451 nopt++; 00452 if (nopt+2 >= argc) 00453 WA_error ("need 3 arguments after -filt_sgnl"); 00454 00455 ifilter = option_data->num_sgnl_filters; 00456 00457 sscanf (argv[nopt], "%d", &ival); 00458 if (ival < -1) 00459 WA_error ("-filt_sgnl band mintr maxtr Require: -1 <= band"); 00460 option_data->sgnl_band[ifilter] = ival; 00461 nopt++; 00462 00463 sscanf (argv[nopt], "%d", &ival); 00464 if (ival < 0) 00465 WA_error ("-filt_sgnl band mintr maxtr Require: 0 <= mintr"); 00466 option_data->sgnl_mintr[ifilter] = ival; 00467 nopt++; 00468 00469 sscanf (argv[nopt], "%d", &ival); 00470 if (ival < 0) 00471 WA_error ("-filt_sgnl band mintr maxtr Require: 0 <= maxtr"); 00472 option_data->sgnl_maxtr[ifilter] = ival; 00473 nopt++; 00474 00475 option_data->num_sgnl_filters++; 00476 continue; 00477 } 00478 00479 00480 /*----- -fout -----*/ 00481 if (strncmp(argv[nopt], "-fout", 5) == 0) 00482 { 00483 option_data->fout = 1; 00484 nopt++; 00485 continue; 00486 } 00487 00488 00489 /*----- -rout -----*/ 00490 if (strncmp(argv[nopt], "-rout", 5) == 0) 00491 { 00492 option_data->rout = 1; 00493 nopt++; 00494 continue; 00495 } 00496 00497 00498 /*----- -cout -----*/ 00499 if (strncmp(argv[nopt], "-cout", 5) == 0) 00500 { 00501 option_data->cout = 1; 00502 nopt++; 00503 continue; 00504 } 00505 00506 00507 /*----- -vout -----*/ 00508 if (strncmp(argv[nopt], "-vout", 5) == 0) 00509 { 00510 option_data->vout = 1; 00511 nopt++; 00512 continue; 00513 } 00514 00515 00516 /*----- -stat_first -----*/ 00517 if (strncmp(argv[nopt], "-stat_first", 11) == 0) 00518 { 00519 option_data->stat_first = 1; 00520 nopt++; 00521 continue; 00522 } 00523 00524 00525 /*----- -bucket filename -----*/ 00526 if (strncmp(argv[nopt], "-bucket", 6) == 0) 00527 { 00528 nopt++; 00529 if (nopt >= argc) WA_error ("need file prefixname after -bucket "); 00530 option_data->bucket_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); 00531 MTEST (option_data->bucket_filename); 00532 strcpy (option_data->bucket_filename, argv[nopt]); 00533 nopt++; 00534 continue; 00535 } 00536 00537 00538 /*----- -coefts filename -----*/ 00539 if (strncmp(argv[nopt], "-coefts", 8) == 0) 00540 { 00541 nopt++; 00542 if (nopt >= argc) WA_error ("need file prefixname after -coefts "); 00543 option_data->coefts_filename = 00544 malloc (sizeof(char)*MAX_NAME_LENGTH); 00545 MTEST (option_data->coefts_filename); 00546 strcpy (option_data->coefts_filename, argv[nopt]); 00547 nopt++; 00548 continue; 00549 } 00550 00551 00552 /*----- -fitts filename -----*/ 00553 if (strncmp(argv[nopt], "-fitts", 7) == 0) 00554 { 00555 nopt++; 00556 if (nopt >= argc) WA_error ("need file prefixname after -fitts "); 00557 option_data->fitts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); 00558 MTEST (option_data->fitts_filename); 00559 strcpy (option_data->fitts_filename, argv[nopt]); 00560 nopt++; 00561 continue; 00562 } 00563 00564 00565 /*----- -sgnlts filename -----*/ 00566 if (strncmp(argv[nopt], "-sgnlts", 7) == 0) 00567 { 00568 nopt++; 00569 if (nopt >= argc) WA_error ("need file prefixname after -sgnlts "); 00570 option_data->sgnlts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); 00571 MTEST (option_data->sgnlts_filename); 00572 strcpy (option_data->sgnlts_filename, argv[nopt]); 00573 nopt++; 00574 continue; 00575 } 00576 00577 00578 /*----- -errts filename -----*/ 00579 if (strncmp(argv[nopt], "-errts", 6) == 0) 00580 { 00581 nopt++; 00582 if (nopt >= argc) WA_error ("need file prefixname after -errts "); 00583 option_data->errts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); 00584 MTEST (option_data->errts_filename); 00585 strcpy (option_data->errts_filename, argv[nopt]); 00586 nopt++; 00587 continue; 00588 } 00589 00590 00591 /*----- unknown command -----*/ 00592 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]); 00593 WA_error (message); 00594 00595 } 00596 00597 00598 } |
|
Definition at line 853 of file 3dWavelets.c. References DSET_NUM_TIMES, FWT_1d_pass_filter(), FWT_1d_stop_filter(), i, MAX_NAME_LENGTH, my_log2(), p, powerof2(), and q. Referenced by initialize_program().
00860 { 00861 char message[MAX_NAME_LENGTH]; /* error message */ 00862 int nt; /* number of images in input 3d+time dataset */ 00863 int NFirst; /* first image from input 3d+time dataset to use */ 00864 int NLast; /* last image from input 3d+time dataset to use */ 00865 int N; /* number of usable time points */ 00866 int ifilter; /* filter number index */ 00867 int i; /* filter coefficient index */ 00868 int f, q, p; /* numbers of model parameters */ 00869 00870 00871 /*----- Initialize local variables -----*/ 00872 if (option_data->input1D_filename != NULL) 00873 nt = fmri_length; 00874 else 00875 nt = DSET_NUM_TIMES (dset_time); 00876 00877 00878 /*----- Determine the number of usable time points -----*/ 00879 NFirst = option_data->NFirst; 00880 00881 NLast = option_data->NLast; 00882 if (NLast > nt-1) NLast = nt-1; 00883 00884 N = NLast - NFirst + 1; 00885 N = powerof2(my_log2(N)); 00886 NLast = N + NFirst + 1; 00887 00888 option_data->N = N; 00889 option_data->NLast = NLast; 00890 00891 00892 /*----- Initialize the filters -----*/ 00893 option_data->stop_filter = 00894 FWT_1d_stop_filter (option_data->num_stop_filters, option_data->stop_band, 00895 option_data->stop_mintr, option_data->stop_maxtr, 00896 option_data->NFirst, option_data->N); 00897 00898 option_data->base_filter = 00899 FWT_1d_pass_filter (option_data->num_base_filters, option_data->base_band, 00900 option_data->base_mintr, option_data->base_maxtr, 00901 option_data->NFirst, option_data->N); 00902 00903 option_data->sgnl_filter = 00904 FWT_1d_pass_filter (option_data->num_sgnl_filters, option_data->sgnl_band, 00905 option_data->sgnl_mintr, option_data->sgnl_maxtr, 00906 option_data->NFirst, option_data->N); 00907 00908 00909 /*----- Count numbers of model parameters -----*/ 00910 f = 0; 00911 for (i = 0; i < N; i++) 00912 if (option_data->stop_filter[i] == 0.0) 00913 { 00914 f++; 00915 option_data->base_filter[i] = 0.0; 00916 option_data->sgnl_filter[i] = 0.0; 00917 } 00918 option_data->f = f; 00919 00920 q = 0; 00921 for (i = 0; i < N; i++) 00922 if (option_data->base_filter[i] == 1.0) 00923 { 00924 q++; 00925 option_data->sgnl_filter[i] = 1.0; 00926 } 00927 option_data->q = q; 00928 00929 p = 0; 00930 for (i = 0; i < N; i++) 00931 if (option_data->sgnl_filter[i] == 1.0) 00932 { 00933 p++; 00934 } 00935 option_data->p = p; 00936 00937 00938 return; 00939 } |
|
Definition at line 181 of file 3dWavelets.c. References MAX_FILTERS.
00185 { 00186 int is; /* filter index */ 00187 00188 00189 /*----- Initialize default values -----*/ 00190 option_data->NFirst = 0; 00191 option_data->NLast = 32767; 00192 option_data->N = 0; 00193 option_data->f = 0; 00194 option_data->p = 0; 00195 option_data->q = 0; 00196 option_data->fdisp = -1.0; 00197 option_data->wavelet_type = 0; 00198 00199 00200 /*----- Initialize filter specifications -----*/ 00201 option_data->num_stop_filters = 0; 00202 option_data->num_base_filters = 0; 00203 option_data->num_sgnl_filters = 0; 00204 for (is = 0; is < MAX_FILTERS; is++) 00205 { 00206 option_data->stop_band[is] = 0; 00207 option_data->stop_mintr[is] = 0; 00208 option_data->stop_maxtr[is] = 0; 00209 option_data->base_band[is] = 0; 00210 option_data->base_mintr[is] = 0; 00211 option_data->base_maxtr[is] = 0; 00212 option_data->sgnl_band[is] = 0; 00213 option_data->sgnl_mintr[is] = 0; 00214 option_data->sgnl_maxtr[is] = 0; 00215 } 00216 option_data->stop_filter = NULL; 00217 option_data->base_filter = NULL; 00218 option_data->sgnl_filter = NULL; 00219 00220 00221 /*----- Initialize output flags -----*/ 00222 option_data->fout = 0; 00223 option_data->rout = 0; 00224 option_data->cout = 0; 00225 option_data->vout = 0; 00226 option_data->stat_first = 0; 00227 00228 00229 /*----- Initialize file name character strings -----*/ 00230 option_data->input_filename = NULL; 00231 option_data->mask_filename = NULL; 00232 option_data->input1D_filename = NULL; 00233 option_data->bucket_filename = NULL; 00234 option_data->coefts_filename = NULL; 00235 option_data->fitts_filename = NULL; 00236 option_data->sgnlts_filename = NULL; 00237 option_data->errts_filename = NULL; 00238 00239 } |
|
Definition at line 1067 of file 3dWavelets.c. References allocate_memory(), argc, check_for_valid_inputs(), get_options(), initialize_filters(), malloc, and read_input_data().
01087 { 01088 01089 /*----- Allocate memory -----*/ 01090 *option_data = (WA_options *) malloc (sizeof(WA_options)); 01091 01092 01093 /*----- Get command line inputs -----*/ 01094 get_options (argc, argv, *option_data); 01095 01096 01097 /*----- Read input data -----*/ 01098 read_input_data (*option_data, dset_time, mask_dset, fmri_data, fmri_length); 01099 01100 01101 /*----- Check for valid inputs -----*/ 01102 check_for_valid_inputs (*option_data, *dset_time, *mask_dset); 01103 01104 01105 /*----- Initialize the filters and control variables -----*/ 01106 initialize_filters (*option_data, *dset_time, *mask_dset, *fmri_length); 01107 01108 01109 /*----- Allocate memory for output volumes -----*/ 01110 if ((*option_data)->input1D_filename == NULL) 01111 allocate_memory (*option_data, *dset_time, 01112 coef_vol, mse_vol, ffull_vol, rfull_vol, 01113 coefts_vol, fitts_vol, sgnlts_vol, errts_vol); 01114 01115 } |
|
Definition at line 2026 of file 3dWavelets.c. References argc, calculate_results(), initialize_program(), WA_options::input1D_filename, output_results(), PROGRAM_AUTHOR, PROGRAM_INITIAL, PROGRAM_LATEST, PROGRAM_NAME, terminate_program(), and THD_delete_3dim_dataset().
02031 { 02032 WA_options * option_data; /* wavelet analysis options */ 02033 THD_3dim_dataset * dset_time = NULL; /* input 3d+time data set */ 02034 THD_3dim_dataset * mask_dset = NULL; /* input mask data set */ 02035 float * fmri_data = NULL; /* input fMRI time series data */ 02036 int fmri_length; /* length of fMRI time series */ 02037 02038 float ** coef_vol = NULL; /* array of volumes for model parameters */ 02039 float * mse_vol = NULL; /* volume of full model mean square error */ 02040 float * ffull_vol = NULL; /* volume of full model F-statistics */ 02041 float * rfull_vol = NULL; /* volume of full model R^2 stats. */ 02042 02043 float ** coefts_vol = NULL; /* volumes for wavelet transform coefs.*/ 02044 float ** fitts_vol = NULL; /* volumes for wavelet filtered time series */ 02045 float ** sgnlts_vol = NULL; /* volumes for signal model fit to input data */ 02046 float ** errts_vol = NULL; /* volumes for residual errors */ 02047 02048 02049 /*----- Identify software -----*/ 02050 printf ("\n\n"); 02051 printf ("Program: %s \n", PROGRAM_NAME); 02052 printf ("Author: %s \n", PROGRAM_AUTHOR); 02053 printf ("Initial Release: %s \n", PROGRAM_INITIAL); 02054 printf ("Latest Revision: %s \n", PROGRAM_LATEST); 02055 printf ("\n"); 02056 02057 02058 /*----- Program initialization -----*/ 02059 initialize_program (argc, argv, &option_data, 02060 &dset_time, &mask_dset, &fmri_data, &fmri_length, 02061 &coef_vol, &mse_vol, &ffull_vol, &rfull_vol, 02062 &coefts_vol, &fitts_vol, &sgnlts_vol, &errts_vol); 02063 02064 02065 /*----- Perform wavelet analysis -----*/ 02066 calculate_results (option_data, dset_time, mask_dset, fmri_data, fmri_length, 02067 coef_vol, mse_vol, ffull_vol, rfull_vol, 02068 coefts_vol, fitts_vol, sgnlts_vol, errts_vol); 02069 02070 02071 /*----- Deallocate memory for input datasets -----*/ 02072 if (dset_time != NULL) 02073 { THD_delete_3dim_dataset (dset_time, False); dset_time = NULL; } 02074 if (mask_dset != NULL) 02075 { THD_delete_3dim_dataset (mask_dset, False); mask_dset = NULL; } 02076 02077 02078 /*----- Write requested output files -----*/ 02079 if (option_data->input1D_filename == NULL) 02080 output_results (argc, argv, option_data, coef_vol, mse_vol, 02081 ffull_vol, rfull_vol, 02082 coefts_vol, fitts_vol, sgnlts_vol, errts_vol); 02083 02084 02085 /*----- Terminate program -----*/ 02086 terminate_program (&option_data, &fmri_data, 02087 &coef_vol, &mse_vol, &ffull_vol, &rfull_vol, 02088 &coefts_vol, &fitts_vol, &sgnlts_vol, &errts_vol); 02089 02090 exit(0); 02091 } |
|
Definition at line 1852 of file 3dWavelets.c. References argc, write_bucket_data(), and write_ts_array().
01867 { 01868 int N; /* number of usable data points */ 01869 01870 01871 /*----- Initialize local variables -----*/ 01872 N = option_data->N; 01873 01874 01875 /*----- Write the bucket dataset -----*/ 01876 if (option_data->bucket_filename != NULL) 01877 write_bucket_data (argc, argv, option_data, coef_vol, 01878 mse_vol, ffull_vol, rfull_vol); 01879 01880 01881 01882 /*----- Write the forward wavelet transform coefficients -----*/ 01883 if (option_data->coefts_filename != NULL) 01884 write_ts_array (argc, argv, option_data, N, coefts_vol, 01885 option_data->coefts_filename); 01886 01887 01888 /*----- Write the wavelet filtered 3d+time dataset -----*/ 01889 if (option_data->fitts_filename != NULL) 01890 write_ts_array (argc, argv, option_data, N, fitts_vol, 01891 option_data->fitts_filename); 01892 01893 01894 /*----- Write the signal model 3d+time dataset -----*/ 01895 if (option_data->sgnlts_filename != NULL) 01896 write_ts_array (argc, argv, option_data, N, sgnlts_vol, 01897 option_data->sgnlts_filename); 01898 01899 01900 /*----- Write the residual errors 3d+time dataset -----*/ 01901 if (option_data->errts_filename != NULL) 01902 write_ts_array (argc, argv, option_data, N, errts_vol, 01903 option_data->errts_filename); 01904 01905 01906 01907 } |
|
Definition at line 667 of file 3dWavelets.c. References ISVALID_3DIM_DATASET, MAX_NAME_LENGTH, read_time_series(), THD_load_datablock(), THD_open_dataset(), THD_open_one_dataset(), and WA_error(). Referenced by initialize_program().
00675 { 00676 char message[MAX_NAME_LENGTH]; /* error message */ 00677 00678 00679 /*----- Read the input fMRI measurement data -----*/ 00680 if (option_data->input1D_filename != NULL) 00681 { 00682 /*----- Read the input fMRI 1D time series -----*/ 00683 *fmri_data = read_time_series (option_data->input1D_filename, 00684 fmri_length); 00685 if (*fmri_data == NULL) 00686 { 00687 sprintf (message, "Unable to read time series file: %s", 00688 option_data->input1D_filename); 00689 WA_error (message); 00690 } 00691 *dset_time = NULL; 00692 } 00693 00694 else if (option_data->input_filename != NULL) 00695 { 00696 /*----- Read the input 3d+time dataset -----*/ 00697 *dset_time = THD_open_one_dataset (option_data->input_filename); 00698 if (!ISVALID_3DIM_DATASET(*dset_time)) 00699 { 00700 sprintf (message, "Unable to open data file: %s", 00701 option_data->input_filename); 00702 WA_error (message); 00703 } 00704 THD_load_datablock ((*dset_time)->dblk); 00705 00706 if (option_data->mask_filename != NULL) 00707 { 00708 /*----- Read the input mask dataset -----*/ 00709 *mask_dset = THD_open_dataset (option_data->mask_filename); 00710 if (!ISVALID_3DIM_DATASET(*mask_dset)) 00711 { 00712 sprintf (message, "Unable to open mask file: %s", 00713 option_data->mask_filename); 00714 WA_error (message); 00715 } 00716 THD_load_datablock ((*mask_dset)->dblk); 00717 } 00718 } 00719 else 00720 WA_error ("Must specify input data"); 00721 00722 } |
|
Definition at line 609 of file 3dWavelets.c. References far, malloc, MAX_NAME_LENGTH, MRI_FLOAT_PTR, mri_free(), mri_read_1D(), MTEST, MRI_IMAGE::nx, MRI_IMAGE::ny, THD_MAX_NAME, and WA_error(). Referenced by initialize_program(), and read_input_data().
00614 { 00615 char message[MAX_NAME_LENGTH]; /* error message */ 00616 char * cpt; /* pointer to column suffix */ 00617 char filename[THD_MAX_NAME]; /* time series file name w/o column index */ 00618 char subv[THD_MAX_NAME]; /* string containing column index */ 00619 MRI_IMAGE * im, * flim; /* pointers to image structures 00620 -- used to read 1D ASCII */ 00621 float * far; /* pointer to MRI_IMAGE floating point data */ 00622 int nx; /* number of time points in time series */ 00623 int ny; /* number of columns in time series file */ 00624 int iy; /* time series file column index */ 00625 int ipt; /* time point index */ 00626 float * ts_data = NULL; /* input time series data */ 00627 00628 00629 /*----- Read the time series file -----*/ 00630 flim = mri_read_1D(ts_filename) ; 00631 if (flim == NULL) 00632 { 00633 sprintf (message, "Unable to read time series file: %s", ts_filename); 00634 WA_error (message); 00635 } 00636 00637 00638 /*----- Set pointer to data, and set dimensions -----*/ 00639 far = MRI_FLOAT_PTR(flim); 00640 nx = flim->nx; 00641 ny = flim->ny; iy = 0 ; 00642 if( ny > 1 ){ 00643 fprintf(stderr,"WARNING: time series %s has %d columns\n",ts_filename,ny); 00644 } 00645 00646 00647 /*----- Save the time series data -----*/ 00648 *ts_length = nx; 00649 ts_data = (float *) malloc (sizeof(float) * nx); 00650 MTEST (ts_data); 00651 for (ipt = 0; ipt < nx; ipt++) 00652 ts_data[ipt] = far[ipt + iy*nx]; 00653 00654 00655 mri_free (flim); flim = NULL; 00656 00657 return (ts_data); 00658 } |
|
Definition at line 1166 of file 3dWavelets.c. References p. Referenced by calculate_results().
01191 { 01192 int ip; /* parameter index */ 01193 int p; /* total number of parameters */ 01194 int it; /* time point index */ 01195 int N; /* number of usable data points */ 01196 01197 01198 /*----- Initialize local variables -----*/ 01199 p = option_data->p; 01200 N = option_data->N; 01201 01202 01203 /*----- Saved wavelet coefficients -----*/ 01204 for (ip = 0; ip < p; ip++) 01205 { 01206 coef_vol[ip][iv] = coef[ip]; 01207 } 01208 01209 01210 /*----- Save full model mean square error -----*/ 01211 mse_vol[iv] = mse; 01212 01213 01214 /*----- Save regression F-statistic -----*/ 01215 ffull_vol[iv] = ffull; 01216 01217 01218 /*----- Save R^2 values -----*/ 01219 rfull_vol[iv] = rfull; 01220 01221 01222 /*----- Save the forward wavelet transform coefficients -----*/ 01223 if (coefts_vol != NULL) 01224 for (it = 0; it < N; it++) 01225 coefts_vol[it][iv] = coefts[it]; 01226 01227 01228 /*----- Save the wavelet filtered time series -----*/ 01229 if (fitts_vol != NULL) 01230 for (it = 0; it < N; it++) 01231 fitts_vol[it][iv] = fitts[it]; 01232 01233 01234 /*----- Save the fitted signal model time series -----*/ 01235 if (sgnlts_vol != NULL) 01236 for (it = 0; it < N; it++) 01237 sgnlts_vol[it][iv] = sgnlts[it]; 01238 01239 01240 /*----- Save the residual errors time series -----*/ 01241 if (errts_vol != NULL) 01242 for (it = 0; it < N; it++) 01243 errts_vol[it][iv] = errts[it]; 01244 01245 01246 } |
|
Definition at line 1913 of file 3dWavelets.c.
01928 { 01929 int p; /* number of parameters in full model */ 01930 int ip; /* parameter index */ 01931 int is; /* stimulus index */ 01932 int it; /* time index */ 01933 int N; /* number of usable data points */ 01934 01935 01936 /*----- Initialize local variables -----*/ 01937 p = (*option_data)->p; 01938 N = (*option_data)->N; 01939 01940 01941 /*----- Deallocate memory for option data -----*/ 01942 if ((*option_data)->stop_filter != NULL) free ((*option_data)->stop_filter); 01943 if ((*option_data)->base_filter != NULL) free ((*option_data)->base_filter); 01944 if ((*option_data)->sgnl_filter != NULL) free ((*option_data)->sgnl_filter); 01945 if ((*option_data)->input_filename != NULL) 01946 free ((*option_data)->input_filename); 01947 if ((*option_data)->mask_filename != NULL) 01948 free ((*option_data)->mask_filename); 01949 if ((*option_data)->input1D_filename != NULL) 01950 free ((*option_data)->input1D_filename); 01951 if ((*option_data)->bucket_filename != NULL) 01952 free ((*option_data)->bucket_filename); 01953 if ((*option_data)->coefts_filename != NULL) 01954 free ((*option_data)->coefts_filename); 01955 if ((*option_data)->fitts_filename != NULL) 01956 free ((*option_data)->fitts_filename); 01957 if ((*option_data)->sgnlts_filename != NULL) 01958 free ((*option_data)->sgnlts_filename); 01959 if ((*option_data)->errts_filename != NULL) 01960 free ((*option_data)->errts_filename); 01961 free (*option_data); *option_data = NULL; 01962 01963 01964 /*----- Deallocate memory for fMRI time series data -----*/ 01965 if (*fmri_data != NULL) { free (*fmri_data); *fmri_data = NULL; } 01966 01967 01968 /*----- Deallocate space for volume data -----*/ 01969 if (*coef_vol != NULL) 01970 { 01971 for (ip = 0; ip < p; ip++) 01972 { 01973 if ((*coef_vol)[ip] != NULL) 01974 { free ((*coef_vol)[ip]); (*coef_vol)[ip] = NULL; } 01975 } 01976 free (*coef_vol); *coef_vol = NULL; 01977 } 01978 01979 if (*mse_vol != NULL) { free (*mse_vol); *mse_vol = NULL; } 01980 if (*ffull_vol != NULL) { free (*ffull_vol); *ffull_vol = NULL; } 01981 if (*rfull_vol != NULL) { free (*rfull_vol); *rfull_vol = NULL; } 01982 01983 01984 /*----- Deallocate space for forward wavelet transform coefficients -----*/ 01985 if (*coefts_vol != NULL) 01986 { 01987 for (it = 0; it < N; it++) 01988 { free ((*coefts_vol)[it]); (*coefts_vol)[it] = NULL; } 01989 free (*coefts_vol); *coefts_vol = NULL; 01990 } 01991 01992 01993 /*----- Deallocate space for wavelet filtered time series -----*/ 01994 if (*fitts_vol != NULL) 01995 { 01996 for (it = 0; it < N; it++) 01997 { free ((*fitts_vol)[it]); (*fitts_vol)[it] = NULL; } 01998 free (*fitts_vol); *fitts_vol = NULL; 01999 } 02000 02001 02002 /*----- Deallocate space for signal model time series -----*/ 02003 if (*sgnlts_vol != NULL) 02004 { 02005 for (it = 0; it < N; it++) 02006 { free ((*sgnlts_vol)[it]); (*sgnlts_vol)[it] = NULL; } 02007 free (*sgnlts_vol); *sgnlts_vol = NULL; 02008 } 02009 02010 02011 /*----- Deallocate space for residual errors time series -----*/ 02012 if (*errts_vol != NULL) 02013 { 02014 for (it = 0; it < N; it++) 02015 { free ((*errts_vol)[it]); (*errts_vol)[it] = NULL; } 02016 free (*errts_vol); *errts_vol = NULL; 02017 } 02018 02019 02020 } |
|
Definition at line 1619 of file 3dWavelets.c. References ADN_datum_all, ADN_directory_name, ADN_func_type, ADN_malloc_type, ADN_none, ADN_ntt, ADN_nvals, ADN_prefix, ADN_type, argc, attach_sub_brick(), DATABLOCK_MEM_MALLOC, THD_3dim_dataset::daxes, DSET_BRIKNAME, DSET_HEADNAME, DSET_NUM_TIMES, EDIT_dset_items(), EDIT_empty_copy(), FUNC_BUCK_TYPE, FUNC_FIM_TYPE, FUNC_THR_TYPE, HEAD_FUNC_TYPE, malloc, MAX_NAME_LENGTH, MTEST, my_log2(), THD_dataxes::nxx, THD_dataxes::nyy, THD_dataxes::nzz, p, powerof2(), PROGRAM_NAME, q, THD_delete_3dim_dataset(), THD_is_file(), THD_load_statistics(), THD_open_one_dataset(), THD_write_3dim_dataset(), tross_Append_History(), tross_Copy_History(), and tross_Make_History(). Referenced by output_parameters(), and output_results().
01630 { 01631 THD_3dim_dataset * old_dset = NULL; /* prototype dataset */ 01632 THD_3dim_dataset * new_dset = NULL; /* output bucket dataset */ 01633 char output_prefix[MAX_NAME_LENGTH]; /* prefix name for bucket dataset */ 01634 char output_session[MAX_NAME_LENGTH]; /* directory for bucket dataset */ 01635 int nbricks; /* number of sub-bricks in bucket dataset */ 01636 short ** bar = NULL; /* bar[ib] points to data for sub-brick #ib */ 01637 01638 int brick_type; /* indicates statistical type of sub-brick */ 01639 int brick_coef; /* regression coefficient index for sub-brick */ 01640 char brick_label[MAX_NAME_LENGTH]; /* character string label for sub-brick */ 01641 01642 int ierror; /* number of errors in editing data */ 01643 float * volume; /* volume of floating point data */ 01644 01645 int N; /* number of usable data points */ 01646 int NFirst; 01647 int f; /* number of parameters removed by stop filter */ 01648 int p; /* number of parameters in the signal model */ 01649 int q; /* number of parameters in the rdcd model */ 01650 int nxyz; /* total number of voxels */ 01651 int nt; /* number of images in input 3d+time dataset */ 01652 int it; 01653 int ilag; /* lag index */ 01654 int icoef; /* coefficient index */ 01655 int ibrick; /* sub-brick index */ 01656 int dof, ndof, ddof; /* degrees of freedom */ 01657 char label[MAX_NAME_LENGTH]; /* general label for sub-bricks */ 01658 01659 01660 /*----- read prototype dataset -----*/ 01661 old_dset = THD_open_one_dataset (option_data->input_filename); 01662 01663 01664 /*----- Initialize local variables -----*/ 01665 nxyz = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz; 01666 nt = DSET_NUM_TIMES (old_dset); 01667 01668 01669 f = option_data->f; 01670 q = option_data->q; 01671 p = option_data->p; 01672 N = option_data->N; 01673 NFirst = option_data->NFirst; 01674 printf ("f = %d q = %d p = %d N = %d \n", f, q, p, N); 01675 01676 01677 /*----- Calculate number of sub-bricks in the bucket -----*/ 01678 nbricks = p * option_data->cout 01679 + option_data->rout + option_data->fout + option_data->vout; 01680 printf ("nbricks = %d \n", nbricks); 01681 01682 01683 strcpy (output_prefix, option_data->bucket_filename); 01684 strcpy (output_session, "./"); 01685 01686 01687 /*----- allocate memory -----*/ 01688 bar = (short **) malloc (sizeof(short *) * nbricks); 01689 MTEST (bar); 01690 01691 01692 /*-- make an empty copy of prototype dataset, for eventual output --*/ 01693 new_dset = EDIT_empty_copy (old_dset); 01694 01695 01696 /*----- Record history of dataset -----*/ 01697 tross_Copy_History( old_dset , new_dset ) ; 01698 tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ; 01699 sprintf (label, "Output prefix: %s", output_prefix); 01700 tross_Append_History ( new_dset, label); 01701 01702 01703 /*----- delete prototype dataset -----*/ 01704 THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ; 01705 01706 01707 /*----- Modify some structural properties. Note that the nbricks 01708 just make empty sub-bricks, without any data attached. -----*/ 01709 ierror = EDIT_dset_items (new_dset, 01710 ADN_prefix, output_prefix, 01711 ADN_directory_name, output_session, 01712 ADN_type, HEAD_FUNC_TYPE, 01713 ADN_func_type, FUNC_BUCK_TYPE, 01714 ADN_datum_all, MRI_short , 01715 ADN_ntt, 0, /* no time */ 01716 ADN_nvals, nbricks, 01717 ADN_malloc_type, DATABLOCK_MEM_MALLOC , 01718 ADN_none ) ; 01719 01720 if( ierror > 0 ) 01721 { 01722 fprintf(stderr, 01723 "*** %d errors in attempting to create bucket dataset!\n", 01724 ierror); 01725 exit(1); 01726 } 01727 01728 if (THD_is_file(DSET_HEADNAME(new_dset))) 01729 { 01730 fprintf(stderr, 01731 "*** Output dataset file %s already exists--cannot continue!\n", 01732 DSET_HEADNAME(new_dset)); 01733 exit(1); 01734 } 01735 01736 01737 /*----- Attach individual sub-bricks to the bucket dataset -----*/ 01738 ibrick = 0; 01739 01740 01741 /*----- Full model wavelet coefficients -----*/ 01742 if (option_data->cout) 01743 { 01744 /*----- User can choose to place full model stats first -----*/ 01745 if (option_data->stat_first) 01746 ibrick += option_data->vout + option_data->rout + option_data->fout; 01747 01748 icoef = 0; 01749 for (it = 0; it < N; it++) 01750 { 01751 if (option_data->sgnl_filter[it]) 01752 { 01753 /*----- Wavelet coefficient -----*/ 01754 brick_type = FUNC_FIM_TYPE; 01755 01756 { 01757 int band, mintr, maxtr; 01758 01759 if (it == 0) 01760 { 01761 band = -1; 01762 mintr = 0; 01763 maxtr = N-1; 01764 } 01765 else 01766 { 01767 band = my_log2(it); 01768 mintr = (it - powerof2(band)) * powerof2(my_log2(N)-band); 01769 maxtr = mintr + powerof2(my_log2(N)-band) - 1; 01770 } 01771 01772 mintr += NFirst; 01773 maxtr += NFirst; 01774 01775 if (option_data->base_filter[it]) 01776 sprintf (brick_label, "B(%2d)[%3d,%3d]", band, mintr, maxtr); 01777 else 01778 sprintf (brick_label, "S(%2d)[%3d,%3d]", band, mintr, maxtr); 01779 } 01780 01781 volume = coef_vol[icoef]; 01782 attach_sub_brick (new_dset, ibrick, volume, nxyz, 01783 brick_type, brick_label, 0, 0, 0, bar); 01784 01785 ibrick++; 01786 icoef++; 01787 } 01788 01789 } /* End loop over Wavelet coefficients */ 01790 } /* if (option_data->cout) */ 01791 01792 01793 01794 /*----- Statistics for full model -----*/ 01795 if (option_data->stat_first) ibrick = 0; 01796 01797 /*----- Full model MSE -----*/ 01798 if (option_data->vout) 01799 { 01800 brick_type = FUNC_FIM_TYPE; 01801 sprintf (brick_label, "Full MSE"); 01802 volume = mse_vol; 01803 attach_sub_brick (new_dset, ibrick, volume, nxyz, 01804 brick_type, brick_label, 0, 0, 0, bar); 01805 ibrick++; 01806 } 01807 01808 /*----- Full model R^2 -----*/ 01809 if (option_data->rout) 01810 { 01811 brick_type = FUNC_THR_TYPE; 01812 sprintf (brick_label, "Full R^2"); 01813 volume = rfull_vol; 01814 attach_sub_brick (new_dset, ibrick, volume, nxyz, 01815 brick_type, brick_label, 0, 0, 0, bar); 01816 ibrick++; 01817 } 01818 01819 /*----- Full model F-stat -----*/ 01820 if (option_data->fout) 01821 { 01822 brick_type = FUNC_FT_TYPE; 01823 ndof = p - q; 01824 ddof = N - f - p; 01825 sprintf (brick_label, "Full F-stat"); 01826 volume = ffull_vol; 01827 attach_sub_brick (new_dset, ibrick, volume, nxyz, 01828 brick_type, brick_label, 0, ndof, ddof, bar); 01829 ibrick++; 01830 } 01831 01832 01833 /*----- write bucket data set -----*/ 01834 THD_load_statistics (new_dset); 01835 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ; 01836 fprintf(stderr,"Writing `bucket' dataset "); 01837 fprintf(stderr,"into %s\n", DSET_BRIKNAME(new_dset)); 01838 01839 01840 /*----- deallocate memory -----*/ 01841 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ; 01842 01843 } |
|
Definition at line 1445 of file 3dWavelets.c. References ADN_brick_fac, ADN_datum_all, ADN_label1, ADN_malloc_type, ADN_none, ADN_ntt, ADN_nvals, ADN_prefix, ADN_self_name, argc, DATABLOCK_MEM_MALLOC, THD_3dim_dataset::daxes, THD_3dim_dataset::dblk, THD_datablock::diskptr, DSET_BRICK, DSET_BRIKNAME, EDIT_coerce_autoscale_new(), EDIT_dset_items(), EDIT_empty_copy(), free, THD_diskptr::header_name, malloc, mri_fix_data_pointer(), MTEST, THD_dataxes::nxx, THD_dataxes::nyy, THD_dataxes::nzz, PROGRAM_NAME, THD_delete_3dim_dataset(), THD_is_file(), THD_load_statistics(), THD_open_one_dataset(), THD_write_3dim_dataset(), tross_Append_History(), tross_commandline(), tross_Copy_History(), and tross_multi_Append_History(). Referenced by main(), and output_results().
01454 { 01455 const float EPSILON = 1.0e-10; 01456 01457 THD_3dim_dataset * dset = NULL; /* input afni data set pointer */ 01458 THD_3dim_dataset * new_dset = NULL; /* output afni data set pointer */ 01459 int ib; /* sub-brick index */ 01460 int ierror; /* number of errors in editing data */ 01461 int nxyz; /* total number of voxels */ 01462 float factor; /* factor is new scale factor for sub-brick #ib */ 01463 char * input_filename; /* input afni data set file name */ 01464 short ** bar = NULL; /* bar[ib] points to data for sub-brick #ib */ 01465 float * fbuf; /* float buffer */ 01466 float * volume; /* pointer to volume of data */ 01467 char label[80]; /* label for output file */ 01468 01469 01470 /*----- Initialize local variables -----*/ 01471 input_filename = option_data->input_filename; 01472 dset = THD_open_one_dataset (input_filename); 01473 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz; 01474 01475 01476 /*----- allocate memory -----*/ 01477 bar = (short **) malloc (sizeof(short *) * ts_length); MTEST (bar); 01478 fbuf = (float *) malloc (sizeof(float) * ts_length); MTEST (fbuf); 01479 01480 01481 /*-- make an empty copy of the prototype dataset, for eventual output --*/ 01482 new_dset = EDIT_empty_copy (dset); 01483 01484 01485 /*----- Record history of dataset -----*/ 01486 tross_Copy_History( dset , new_dset ) ; 01487 01488 { char * commandline = tross_commandline( PROGRAM_NAME , argc , argv ) ; 01489 sprintf (label, "Output prefix: %s", output_filename); 01490 if( commandline != NULL ) 01491 tross_multi_Append_History( new_dset , commandline,label,NULL ) ; 01492 else 01493 tross_Append_History ( new_dset, label); 01494 free(commandline) ; 01495 } 01496 01497 /*----- Delete prototype dataset -----*/ 01498 THD_delete_3dim_dataset (dset, False); dset = NULL ; 01499 01500 01501 ierror = EDIT_dset_items (new_dset, 01502 ADN_prefix, output_filename, 01503 ADN_label1, output_filename, 01504 ADN_self_name, output_filename, 01505 ADN_malloc_type, DATABLOCK_MEM_MALLOC, 01506 ADN_datum_all, MRI_short, 01507 ADN_nvals, ts_length, 01508 ADN_ntt, ts_length, 01509 ADN_none); 01510 01511 if( ierror > 0 ){ 01512 fprintf(stderr, 01513 "*** %d errors in attempting to create output dataset!\n", ierror ) ; 01514 exit(1) ; 01515 } 01516 01517 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){ 01518 fprintf(stderr, 01519 "*** Output dataset file %s already exists--cannot continue!\a\n", 01520 new_dset->dblk->diskptr->header_name ) ; 01521 exit(1) ; 01522 } 01523 01524 01525 /*----- attach bricks to new data set -----*/ 01526 for (ib = 0; ib < ts_length; ib++) 01527 { 01528 01529 /*----- Set pointer to appropriate volume -----*/ 01530 volume = vol_array[ib]; 01531 01532 /*----- Allocate memory for output sub-brick -----*/ 01533 bar[ib] = (short *) malloc (sizeof(short) * nxyz); 01534 MTEST (bar[ib]); 01535 01536 /*----- Convert data type to short for this sub-brick -----*/ 01537 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume, 01538 MRI_short, bar[ib]); 01539 if (factor < EPSILON) factor = 0.0; 01540 else factor = 1.0 / factor; 01541 fbuf[ib] = factor; 01542 01543 /*----- attach bar[ib] to be sub-brick #ib -----*/ 01544 mri_fix_data_pointer (bar[ib], DSET_BRICK(new_dset,ib)); 01545 } 01546 01547 01548 /*----- write afni data set -----*/ 01549 01550 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ; 01551 01552 THD_load_statistics (new_dset); 01553 THD_write_3dim_dataset (NULL, NULL, new_dset, True); 01554 fprintf(stderr,"--- Wrote 3d+time dataset into %s\n",DSET_BRIKNAME(new_dset)) ; 01555 01556 01557 /*----- deallocate memory -----*/ 01558 THD_delete_3dim_dataset (new_dset, False); new_dset = NULL ; 01559 free (fbuf); fbuf = NULL; 01560 01561 } |