Doxygen Source Code Documentation
3dConvolve.c File Reference
#include "mrilib.h"
#include "matrix.h"
#include "Deconvolve.c"
#include "randgen.c"
Go to the source code of this file.
Data Structures | |
struct | DC_options |
Defines | |
#define | PROGRAM_NAME "3dConvolve" |
#define | PROGRAM_AUTHOR "B. Douglas Ward" |
#define | PROGRAM_INITIAL "28 June 2001" |
#define | PROGRAM_LATEST "28 Feb 2002" |
#define | RA_error DC_error |
Typedefs | |
typedef DC_options | DC_options |
Functions | |
void | DC_error (char *message) |
void | display_help_menu () |
void | initialize_options (DC_options *option_data) |
void | initialize_stim_options (DC_options *option_data, int num_stimts) |
void | get_options (int argc, char **argv, DC_options *option_data) |
float * | read_time_series (char *ts_filename, int *ts_length) |
void | read_input_data (DC_options *option_data, THD_3dim_dataset **dset_time, THD_3dim_dataset **mask_dset, THD_3dim_dataset **base_dset, THD_3dim_dataset **err_dset, THD_3dim_dataset ***irf_dset, float **base_data, int *base_length, float ***irf_data, int **irf_length, float **censor_array, int *censor_length, int **block_list, int *num_blocks, float ***stimulus, int **stim_length, float **errts_data, int *errts_length) |
void | check_one_output_file (THD_3dim_dataset *dset_time, char *filename) |
void | check_output_files (DC_options *option_data, THD_3dim_dataset *dset_time) |
void | check_for_valid_inputs (DC_options *option_data, THD_3dim_dataset *dset_time, THD_3dim_dataset *mask_dset, int base_length, int *irf_length, float *censor_array, int censor_length, int *block_list, int num_blocks, int *stim_length, int errts_length, int **good_list) |
void | zero_fill_volume (float **fvol, int nxyz) |
void | allocate_memory (DC_options *option_data, float ***outts_vol) |
void | initialize_program (int argc, char **argv, DC_options **option_data, THD_3dim_dataset **dset_time, THD_3dim_dataset **mask_dset, THD_3dim_dataset **base_dset, THD_3dim_dataset **err_dset, THD_3dim_dataset ***irf_dset, float **base_data, int *base_length, float ***irf_data, int **irf_length, float **censor_array, int *censor_length, int **good_list, int **block_list, int *num_blocks, float ***stimulus, int **stim_length, float **errts_data, int *errts_length, float ***predts_vol) |
void | extract_ts_array (THD_3dim_dataset *dset, int iv, float *ts_array) |
void | calc_response (int nt, float *ts_array, int *good_list, int N, matrix x, vector b, float *predts, float *errts, float sigma) |
void | save_voxel (DC_options *option_data, int iv, float *predts, float **predts_vol) |
void | calculate_results (DC_options *option_data, THD_3dim_dataset *dset, THD_3dim_dataset *mask, THD_3dim_dataset *base_dset, THD_3dim_dataset *err_dset, THD_3dim_dataset **irf_dset, float *base_data, int base_length, float **irf_data, int *irf_length, int *good_list, int *block_list, int num_blocks, float **stimulus, int *stim_length, float *errts_data, int errts_length, float **predts_vol) |
float | EDIT_coerce_autoscale_new (int nxyz, int itype, void *ivol, int otype, void *ovol) |
void | write_ts_array (int argc, char **argv, DC_options *option_data, int ts_length, float **vol_array, char *output_filename) |
void | write_one_ts (char *prefix, int ts_length, float **vol_array) |
void | output_results (int argc, char **argv, DC_options *option_data, float **predts_vol) |
void | terminate_program (DC_options **option_data, float **base_data, float ***irf_data, int **irf_length, float **censor_array, int **good_list, int **block_list, float ***stimulus, int **stim_length, float **errts_data, float ***predts_vol) |
int | main (int argc, char **argv) |
Define Documentation
|
Definition at line 37 of file 3dConvolve.c. Referenced by main(). |
|
Definition at line 38 of file 3dConvolve.c. Referenced by main(). |
|
Definition at line 39 of file 3dConvolve.c. Referenced by main(). |
|
Definition at line 36 of file 3dConvolve.c. Referenced by DC_error(), get_options(), main(), and write_ts_array(). |
|
Definition at line 43 of file 3dConvolve.c. Referenced by calc_glt_matrix(), calc_matrices(), check_for_valid_inputs(), check_piece(), get_inputs(), identify_repeats(), initialize_options(), read_piece(), sort_model_indices(), write_afni_data(), and write_piece(). |
Typedef Documentation
|
|
Function Documentation
|
Definition at line 1281 of file 3dConvolve.c. References malloc, MTEST, and zero_fill_volume().
01286 { 01287 int nxyz; /* total number of voxels */ 01288 int it; /* time point index */ 01289 int nt; /* number of images in input 3d+time dataset */ 01290 01291 01292 /*----- Initialize local variables -----*/ 01293 nxyz = option_data->nxyz; 01294 nt = option_data->nt; 01295 01296 01297 /*----- Allocate memory for fitted time series -----*/ 01298 if (option_data->output_filename != NULL) 01299 { 01300 *outts_vol = (float **) malloc (sizeof(float **) * nt); 01301 MTEST (*outts_vol); 01302 for (it = 0; it < nt; it++) 01303 { 01304 zero_fill_volume (&((*outts_vol)[it]), nxyz); 01305 } 01306 } 01307 01308 } |
|
Definition at line 1424 of file 3dConvolve.c. References vector::elts, rand_normal(), vector_destroy(), vector_initialize(), and vector_multiply(). Referenced by calculate_results().
01436 { 01437 vector yhat; /* product Xb */ 01438 int it; /* time point index */ 01439 01440 01441 /*----- initialize vectors -----*/ 01442 vector_initialize (&yhat); 01443 01444 01445 /*----- calculate the system reponse -----*/ 01446 vector_multiply (x, b, &yhat); 01447 01448 01449 /*----- copy to array -----*/ 01450 for (it = 0; it < nt; it++) 01451 predts[it] = ts_array[it]; 01452 for (it = 0; it < N; it++) 01453 predts[good_list[it]] = yhat.elts[it]; 01454 01455 01456 /*----- add measurement noise -----*/ 01457 if (errts != NULL) 01458 for (it = 0; it < nt; it++) 01459 { 01460 predts[it] += errts[it]; 01461 } 01462 01463 01464 /*----- add Gaussian noise -----*/ 01465 if (sigma > 0.0) 01466 { 01467 for (it = 0; it < N; it++) 01468 { 01469 predts[good_list[it]] += rand_normal (0.0, sigma*sigma); 01470 } 01471 } 01472 01473 01474 /*----- dispose of vectors -----*/ 01475 vector_destroy (&yhat); 01476 01477 01478 } |
|
Definition at line 1521 of file 3dConvolve.c. References calc_response(), vector::elts, extract_ts_array(), free, init_indep_var_matrix(), malloc, matrix_destroy(), matrix_initialize(), matrix_sprint(), MTEST, p, q, save_voxel(), vector_create(), vector_destroy(), and vector_initialize().
01542 { 01543 float * ts_array = NULL; /* array of measured data for one voxel */ 01544 float mask_val[1]; /* value of mask at current voxel */ 01545 01546 int p; /* number of parameters in the full model */ 01547 int q; /* number of parameters in the baseline model */ 01548 int polort; /* degree of polynomial for baseline model */ 01549 int m; /* parameter index */ 01550 int n; /* data point index */ 01551 01552 vector coef; /* regression parameters */ 01553 01554 matrix xdata; /* independent variable matrix */ 01555 01556 int ixyz; /* voxel index */ 01557 int nxyz; /* number of voxels per dataset */ 01558 01559 int nt; /* number of images in input 3d+time dataset */ 01560 int N; /* number of usable data points */ 01561 01562 int num_stimts; /* number of stimulus time series */ 01563 int * min_lag; /* minimum time delay for impulse response */ 01564 int * max_lag; /* maximum time delay for impulse response */ 01565 int * nptr; /* number of stim fn. time points per TR */ 01566 01567 int it; /* data point index */ 01568 int ip; /* paramter index */ 01569 int is; /* stimulus index */ 01570 int ilag; /* time lag index */ 01571 float stddev; /* normalized parameter standard deviation */ 01572 char * label; /* string containing stat. summary of results */ 01573 int novar; /* flag for insufficient variation in data */ 01574 01575 float * coefts = NULL; /* regression coefficients */ 01576 float * predts = NULL; /* predicted time series data */ 01577 float * errts = NULL; /* full model residual error time series */ 01578 float sigma; /* std. dev. of additive Gaussian noise */ 01579 01580 01581 /*----- Initialize matrices and vectors -----*/ 01582 matrix_initialize (&xdata); 01583 vector_initialize (&coef); 01584 01585 01586 /*----- Initialize local variables -----*/ 01587 nxyz = option_data->nxyz; 01588 nt = option_data->nt; 01589 num_stimts = option_data->num_stimts; 01590 min_lag = option_data->stim_minlag; 01591 max_lag = option_data->stim_maxlag; 01592 nptr = option_data->stim_nptr; 01593 01594 polort = option_data->polort; 01595 q = option_data->q; 01596 p = option_data->p; 01597 N = option_data->N; 01598 01599 sigma = option_data->sigma; 01600 01601 coefts = (float *) malloc (sizeof(float) * p); MTEST (coefts); 01602 if ((err_dset != NULL) || (errts_data != NULL)) 01603 { 01604 errts = (float *) malloc (sizeof(float) * nt); MTEST (errts); 01605 } 01606 predts = (float *) malloc (sizeof(float) * nt); MTEST (predts); 01607 ts_array = (float *) malloc (sizeof(float) * nt); MTEST (ts_array); 01608 if (option_data->input1D) 01609 for (it = 0; it < nt; it++) ts_array[it] = 0.0; 01610 01611 01612 /*----- Initialize the independent variable matrix -----*/ 01613 init_indep_var_matrix (p, q, polort, nt, N, good_list, block_list, 01614 num_blocks, num_stimts, stimulus, stim_length, 01615 min_lag, max_lag, nptr, NULL, &xdata); 01616 if (option_data->xout) matrix_sprint ("X matrix:", xdata); 01617 01618 01619 vector_create (p, &coef); 01620 01621 01622 /*----- Loop over all voxels -----*/ 01623 for (ixyz = 0; ixyz < nxyz; ixyz++) 01624 { 01625 /*----- Apply mask? -----*/ 01626 if (mask != NULL) 01627 { 01628 extract_ts_array (mask, ixyz, mask_val); 01629 if (mask_val[0] == 0.0) continue; 01630 } 01631 01632 01633 /*----- Extract model parameters for this voxel -----*/ 01634 if (option_data->input1D) 01635 { 01636 if (q > 0) 01637 { 01638 for (ip = 0; ip < q; ip++) 01639 coef.elts[ip] = base_data[ip]; 01640 } 01641 ip = q; 01642 for (is = 0; is < num_stimts; is++) 01643 for (ilag = min_lag[is]; ilag <= max_lag[is]; ilag++) 01644 { 01645 coef.elts[ip] = irf_data[is][ilag-min_lag[is]]; 01646 ip++; 01647 } 01648 } 01649 else 01650 { 01651 if (q > 0) 01652 { 01653 extract_ts_array (base_dset, ixyz, coefts); 01654 for (ip = 0; ip < q; ip++) 01655 coef.elts[ip] = coefts[ip]; 01656 } 01657 ip = q; 01658 for (is = 0; is < num_stimts; is++) 01659 { 01660 extract_ts_array (irf_dset[is], ixyz, coefts); 01661 for (ilag = min_lag[is]; ilag <= max_lag[is]; ilag++) 01662 { 01663 coef.elts[ip] = coefts[ilag-min_lag[is]]; 01664 ip++; 01665 } 01666 } 01667 } 01668 01669 01670 /*----- Extract error time series for this voxel -----*/ 01671 if ((option_data->input1D) && (errts_data != NULL)) 01672 for (it = 0; it < nt; it++) 01673 errts[it] = errts_data[it]; 01674 else if (err_dset != NULL) 01675 extract_ts_array (err_dset, ixyz, errts); 01676 01677 01678 /*----- Extract measured data for this voxel -----*/ 01679 if (option_data->input1D == 0) 01680 extract_ts_array (dset, ixyz, ts_array); 01681 01682 01683 /*----- Calculate the system response for this voxel -----*/ 01684 calc_response (nt, ts_array, good_list, N, xdata, coef, 01685 predts, errts, sigma); 01686 01687 01688 /*----- Save results for this voxel -----*/ 01689 save_voxel (option_data, ixyz, predts, predts_vol); 01690 01691 01692 /*----- Report results for this voxel -----*/ 01693 if (option_data->input1D) 01694 { 01695 printf ("\n\nResults for Voxel #%d: \n", ixyz); 01696 01697 for (it = 0; it < nt; it++) 01698 { 01699 printf (" %f \n", predts[it]); 01700 } 01701 01702 } 01703 01704 } /*----- Loop over voxels -----*/ 01705 01706 01707 /*----- Dispose of matrices and vectors -----*/ 01708 vector_destroy (&coef); 01709 matrix_destroy (&xdata); 01710 01711 if (coefts != NULL) { free (coefts); coefts = NULL; } 01712 if (ts_array != NULL) { free (ts_array); ts_array = NULL; } 01713 if (predts != NULL) { free (predts); predts = NULL; } 01714 if (errts != NULL) { free (errts); errts = NULL; } 01715 } |
|
Definition at line 1079 of file 3dConvolve.c. References check_output_files(), DC_error(), DSET_NVALS, DSET_NX, DSET_NY, DSET_NZ, malloc, MTEST, p, q, and THD_MAX_NAME.
01094 { 01095 char message[THD_MAX_NAME]; /* error message */ 01096 int is; /* stimulus index */ 01097 int num_stimts; /* number of stimulus time series */ 01098 int * min_lag; /* minimum time delay for impulse response */ 01099 int * max_lag; /* maximum time delay for impulse response */ 01100 int * nptr; /* number of stim fn. time points per TR */ 01101 int m; /* number of time delays for impulse response */ 01102 int q; /* number of baseline parameters */ 01103 int p; /* total number of parameters */ 01104 int it; /* time point index */ 01105 int nt; /* number of images in input 3d+time dataset */ 01106 int NFirst; /* first image from input 3d+time dataset to use */ 01107 int NLast; /* last image from input 3d+time dataset to use */ 01108 int N; /* number of usable time points */ 01109 int ib; /* block (run) index */ 01110 int irb; /* time index relative to start of block (run) */ 01111 int exp_length; /* expected length of IRF time series */ 01112 01113 01114 /*----- Initialize local variables -----*/ 01115 nt = option_data->nt; 01116 num_stimts = option_data->num_stimts; 01117 min_lag = option_data->stim_minlag; 01118 max_lag = option_data->stim_maxlag; 01119 nptr = option_data->stim_nptr; 01120 p = option_data->p; 01121 q = option_data->q; 01122 01123 01124 /*----- Check length of censor array -----*/ 01125 if (censor_length < nt) 01126 { 01127 sprintf (message, "Input censor time series file %s is too short", 01128 option_data->censor_filename); 01129 DC_error (message); 01130 } 01131 01132 01133 /*----- Create list of good (usable) time points -----*/ 01134 *good_list = (int *) malloc (sizeof(int) * nt); MTEST (*good_list); 01135 NFirst = option_data->NFirst; 01136 if (NFirst < 0) 01137 { 01138 for (is = 0; is < num_stimts; is++) 01139 if (NFirst < (max_lag[is]+nptr[is]-1)/nptr[is]) 01140 NFirst = (max_lag[is]+nptr[is]-1)/nptr[is]; 01141 } 01142 NLast = option_data->NLast; 01143 if (NLast < 0) NLast = nt; 01144 01145 N = 0; 01146 ib = 0; 01147 for (it = block_list[0]; it < nt; it++) 01148 { 01149 if (ib+1 < num_blocks) 01150 if (it >= block_list[ib+1]) ib++; 01151 01152 irb = it - block_list[ib]; 01153 01154 if ((irb >= NFirst) && (irb <= NLast) && (censor_array[it])) 01155 { 01156 (*good_list)[N] = it; 01157 N++; 01158 } 01159 } 01160 01161 01162 /*----- Check for sufficient data -----*/ 01163 if (N == 0) DC_error ("No usable time points?"); 01164 if (N <= p) 01165 { 01166 sprintf (message, "Insufficient data for estimating %d parameters", p); 01167 DC_error (message); 01168 } 01169 option_data->N = N; 01170 01171 01172 /*----- If mask is used, check for compatible dimensions -----*/ 01173 if (mask_dset != NULL) 01174 { 01175 if ( (DSET_NX(dset_time) != DSET_NX(mask_dset)) 01176 || (DSET_NY(dset_time) != DSET_NY(mask_dset)) 01177 || (DSET_NZ(dset_time) != DSET_NZ(mask_dset)) ) 01178 { 01179 sprintf (message, "%s and %s have incompatible dimensions", 01180 option_data->input_filename, option_data->mask_filename); 01181 DC_error (message); 01182 } 01183 01184 if (DSET_NVALS(mask_dset) != 1 ) 01185 DC_error ("Must specify 1 sub-brick from mask dataset"); 01186 } 01187 01188 01189 /*----- Check number of stimulus time series -----*/ 01190 if (num_stimts < 0) 01191 { 01192 DC_error ("Require: 0 <= num_stimts "); 01193 } 01194 01195 01196 /*----- Check lengths of stimulus time series -----*/ 01197 for (is = 0; is < num_stimts; is++) 01198 { 01199 if (stim_length[is] < nt*nptr[is]) 01200 { 01201 sprintf (message, "Input stimulus time series file %s is too short", 01202 option_data->stim_filename[is]); 01203 DC_error (message); 01204 } 01205 } 01206 01207 01208 /*----- Check length of error time series -----*/ 01209 if ((errts_length > 0) && (errts_length < nt)) 01210 { 01211 sprintf (message, "Input error time series file %s is too short", 01212 option_data->errts_filename); 01213 DC_error (message); 01214 } 01215 01216 01217 /*----- Check number of baseline parameters -----*/ 01218 if (base_length != (option_data->polort+1)*num_blocks) 01219 { 01220 sprintf (message, "Baseline parameters: Expected: %d Actual: %d", 01221 option_data->polort+1, base_length); 01222 DC_error (message); 01223 } 01224 01225 01226 /*----- Check lengths of IRF time series -----*/ 01227 for (is = 0; is < num_stimts; is++) 01228 { 01229 exp_length = max_lag[is] - min_lag[is] + 1; 01230 if (irf_length[is] != exp_length) 01231 { 01232 sprintf (message, "Length of IRF #%d: Expected: %d Actual: %d", 01233 is+1, exp_length, irf_length[is]); 01234 DC_error (message); 01235 } 01236 } 01237 01238 01239 /*----- Check for zero slice offsets with nptr option -----*/ 01240 if (dset_time != NULL) 01241 if (dset_time->taxis->nsl > 0) 01242 for (is = 0; is < num_stimts; is++) 01243 if (nptr[is] > 1) 01244 { 01245 sprintf (message, "Must align all slices to 0 offset time, \n "); 01246 strcat (message, "before using -stim_nptr option. "); 01247 strcat (message, "See program 3dTshift. "); 01248 DC_error (message); 01249 } 01250 01251 01252 /*----- Check whether any of the output files already exist -----*/ 01253 check_output_files (option_data, dset_time); 01254 01255 } |
|
Definition at line 1008 of file 3dConvolve.c. References ADN_label1, ADN_none, ADN_prefix, ADN_self_name, ADN_type, THD_3dim_dataset::dblk, DC_error(), THD_datablock::diskptr, EDIT_dset_items(), EDIT_empty_copy(), GEN_FUNC_TYPE, HEAD_FUNC_TYPE, THD_diskptr::header_name, ISHEAD, THD_delete_3dim_dataset(), THD_is_file(), and THD_MAX_NAME.
01013 { 01014 char message[THD_MAX_NAME]; /* error message */ 01015 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */ 01016 int ierror; /* number of errors in editing data */ 01017 01018 01019 /*----- make an empty copy of input dataset -----*/ 01020 new_dset = EDIT_empty_copy( dset_time ) ; 01021 01022 01023 ierror = EDIT_dset_items( new_dset , 01024 ADN_prefix , filename , 01025 ADN_label1 , filename , 01026 ADN_self_name , filename , 01027 ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE : 01028 GEN_FUNC_TYPE , 01029 ADN_none ) ; 01030 01031 if( ierror > 0 ) 01032 { 01033 sprintf (message, 01034 "*** %d errors in attempting to create output dataset!\n", 01035 ierror); 01036 DC_error (message); 01037 } 01038 01039 if( THD_is_file(new_dset->dblk->diskptr->header_name) ) 01040 { 01041 sprintf (message, 01042 "Output dataset file %s already exists " 01043 " -- cannot continue! ", 01044 new_dset->dblk->diskptr->header_name); 01045 DC_error (message); 01046 } 01047 01048 /*----- deallocate memory -----*/ 01049 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ; 01050 01051 } |
|
Definition at line 1060 of file 3dConvolve.c. References check_one_output_file().
01065 { 01066 01067 if ((option_data->output_filename != NULL) && (option_data->input1D != 1)) 01068 check_one_output_file (dset_time, option_data->output_filename); 01069 01070 } |
|
Definition at line 95 of file 3dConvolve.c. References PROGRAM_NAME.
00096 { 00097 fprintf (stderr, "%s Error: %s \a\n\n", PROGRAM_NAME, message); 00098 exit(1); 00099 } |
|
Definition at line 107 of file 3dConvolve.c.
00108 { 00109 printf ( 00110 "Program to calculate the voxelwise convolution of given impulse response \n" 00111 "function (IRF) time series contained in a 3d+time dataset with a specified \n" 00112 "input stimulus function time series. This program will also calculate \n" 00113 "convolutions involving multiple IRF's and multiple stimulus functions. \n" 00114 "Input options include addition of system noise to the estimated output. \n" 00115 "Output consists of an AFNI 3d+time dataset which contains the estimated \n" 00116 "system response. Alternatively, if all inputs are .1D time series files, \n" 00117 "then the output will be a single .1D time series file. \n" 00118 " \n" 00119 "Usage: \n" 00120 "3dConvolve \n" 00121 "-input fname fname = filename of 3d+time template dataset \n" 00122 "[-input1D] flag to indicate all inputs are .1D time series \n" 00123 "[-mask mname] mname = filename of 3d mask dataset \n" 00124 "[-censor cname] cname = filename of censor .1D time series \n" 00125 "[-concat rname] rname = filename for list of concatenated runs \n" 00126 "[-nfirst fnum] fnum = number of first time point to calculate by \n" 00127 " convolution procedure. (default = max maxlag) \n" 00128 "[-nlast lnum] lnum = number of last time point to calculate by \n" 00129 " convolution procedure. (default = last point) \n" 00130 "[-polort pnum] pnum = degree of polynomial corresponding to the \n" 00131 " baseline model (default: pnum = 1) \n" 00132 "[-base_file bname] bname = file containing baseline parameters \n" 00133 " \n" 00134 "-num_stimts num num = number of input stimulus time series \n" 00135 " (default: num = 0) \n" 00136 "-stim_file k sname sname = filename of kth time series input stimulus\n" 00137 "[-stim_minlag k m] m = minimum time lag for kth input stimulus \n" 00138 " (default: m = 0) \n" 00139 "[-stim_maxlag k n] n = maximum time lag for kth input stimulus \n" 00140 " (default: n = 0) \n" 00141 "[-stim_nptr k p] p = number of stimulus function points per TR \n" 00142 " Note: This option requires 0 slice offset times \n" 00143 " (default: p = 1) \n" 00144 " \n" 00145 "[-iresp k iprefix] iprefix = prefix of 3d+time input dataset which \n" 00146 " contains the kth impulse response function \n" 00147 " \n" 00148 "[-errts eprefix] eprefix = prefix of 3d+time input dataset which \n" 00149 " contains the residual error time series \n" 00150 " (i.e., noise which will be added to the output) \n" 00151 " \n" 00152 "[-sigma s] s = std. dev. of additive Gaussian noise \n" 00153 " (default: s = 0) \n" 00154 "[-seed d] d = seed for random number generator \n" 00155 " (default: d = 1234567) \n" 00156 " \n" 00157 "[-xout] flag to write X matrix to screen \n" 00158 "[-output tprefix] tprefix = prefix of 3d+time output dataset which \n" 00159 " will contain the convolved time series data \n" 00160 " (or tprefix = prefix of .1D output time series \n" 00161 " if the -input1D option is used) \n" 00162 " \n" 00163 ); 00164 00165 exit(0); 00166 } |
|
compute start time of this timeseries * Definition at line 1729 of file 3dConvolve.c. References EDIT_coerce_scale_type(), MCW_vol_amax(), MRI_IS_INT_TYPE, and top.
01731 { 01732 float fac=0.0 , top ; 01733 01734 if( MRI_IS_INT_TYPE(otype) ){ 01735 top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ; 01736 if (top == 0.0) fac = 0.0; 01737 else fac = MRI_TYPE_maxval[otype]/top; 01738 } 01739 01740 EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ; 01741 return ( fac ); 01742 } |
|
Definition at line 1382 of file 3dConvolve.c. References DC_error(), MRI_FLOAT_PTR, mri_free(), and THD_extract_series().
01388 { 01389 MRI_IMAGE * im; /* intermediate float data */ 01390 float * ar; /* pointer to float data */ 01391 int it; /* time index */ 01392 int nv; /* number of volumes */ 01393 01394 01395 /*----- Extract time series from 3d+time data set into MRI_IMAGE -----*/ 01396 im = THD_extract_series (iv, dset, 0); 01397 01398 01399 /*----- Verify extraction -----*/ 01400 if (im == NULL) DC_error ("Unable to extract data from 3d+time dataset"); 01401 01402 01403 /*----- Now extract time series from MRI_IMAGE -----*/ 01404 nv = dset->dblk->nvals ; 01405 ar = MRI_FLOAT_PTR (im); 01406 for (it = 0; it < nv; it++) 01407 { 01408 ts_array[it] = ar[it]; 01409 } 01410 01411 01412 /*----- Release memory -----*/ 01413 mri_free (im); im = NULL; 01414 01415 } |
|
Definition at line 275 of file 3dConvolve.c. References addto_args(), AFNI_logger(), argc, DC_error(), display_help_menu(), initialize_options(), initialize_stim_options(), machdep(), mainENTRY, malloc, MTEST, PRINT_VERSION, PROGRAM_NAME, and THD_MAX_NAME.
00281 { 00282 int nopt = 1; /* input option argument counter */ 00283 int ival, index; /* integer input */ 00284 float fval; /* float input */ 00285 long lval; /* long integer input */ 00286 char message[THD_MAX_NAME]; /* error message */ 00287 int k; /* stimulus time series index */ 00288 00289 00290 /*-- addto the arglist, if user wants to --*/ 00291 { int new_argc ; char ** new_argv ; 00292 addto_args( argc , argv , &new_argc , &new_argv ) ; 00293 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; } 00294 } 00295 00296 00297 /*----- does user request help menu? -----*/ 00298 if (argc < 2 || strcmp(argv[1], "-help") == 0) display_help_menu(); 00299 00300 00301 /*----- add to program log -----*/ 00302 mainENTRY("3dConvolve"); machdep(); PRINT_VERSION("3dConvolve"); 00303 AFNI_logger (PROGRAM_NAME,argc,argv); 00304 00305 00306 /*----- initialize the input options -----*/ 00307 initialize_options (option_data); 00308 00309 00310 /*----- main loop over input options -----*/ 00311 while (nopt < argc ) 00312 { 00313 /*----- -input filename -----*/ 00314 if (strcmp(argv[nopt], "-input") == 0) 00315 { 00316 nopt++; 00317 if (nopt >= argc) DC_error ("need argument after -input "); 00318 option_data->input_filename = malloc (sizeof(char)*THD_MAX_NAME); 00319 MTEST (option_data->input_filename); 00320 strcpy (option_data->input_filename, argv[nopt]); 00321 nopt++; 00322 continue; 00323 } 00324 00325 00326 /*----- -input1D -----*/ 00327 if (strcmp(argv[nopt], "-input1D") == 0) 00328 { 00329 option_data->input1D = 1; 00330 nopt++; 00331 continue; 00332 } 00333 00334 00335 /*----- -mask filename -----*/ 00336 if (strcmp(argv[nopt], "-mask") == 0) 00337 { 00338 nopt++; 00339 if (nopt >= argc) DC_error ("need argument after -mask "); 00340 option_data->mask_filename = malloc (sizeof(char)*THD_MAX_NAME); 00341 MTEST (option_data->mask_filename); 00342 strcpy (option_data->mask_filename, argv[nopt]); 00343 nopt++; 00344 continue; 00345 } 00346 00347 00348 /*----- -censor filename -----*/ 00349 if (strcmp(argv[nopt], "-censor") == 0) 00350 { 00351 nopt++; 00352 if (nopt >= argc) DC_error ("need argument after -censor "); 00353 option_data->censor_filename = 00354 malloc (sizeof(char)*THD_MAX_NAME); 00355 MTEST (option_data->censor_filename); 00356 strcpy (option_data->censor_filename, argv[nopt]); 00357 nopt++; 00358 continue; 00359 } 00360 00361 00362 /*----- -concat filename -----*/ 00363 if (strcmp(argv[nopt], "-concat") == 0) 00364 { 00365 nopt++; 00366 if (nopt >= argc) DC_error ("need argument after -concat "); 00367 option_data->concat_filename = 00368 malloc (sizeof(char)*THD_MAX_NAME); 00369 MTEST (option_data->concat_filename); 00370 strcpy (option_data->concat_filename, argv[nopt]); 00371 nopt++; 00372 continue; 00373 } 00374 00375 00376 /*----- -nfirst num -----*/ 00377 if (strcmp(argv[nopt], "-nfirst") == 0) 00378 { 00379 nopt++; 00380 if (nopt >= argc) DC_error ("need argument after -nfirst "); 00381 sscanf (argv[nopt], "%d", &ival); 00382 if (ival < 0) 00383 DC_error ("illegal argument after -nfirst "); 00384 option_data->NFirst = ival; 00385 nopt++; 00386 continue; 00387 } 00388 00389 00390 /*----- -nlast num -----*/ 00391 if (strcmp(argv[nopt], "-nlast") == 0) 00392 { 00393 nopt++; 00394 if (nopt >= argc) DC_error ("need argument after -nlast "); 00395 sscanf (argv[nopt], "%d", &ival); 00396 if (ival < 0) 00397 DC_error ("illegal argument after -nlast "); 00398 option_data->NLast = ival; 00399 nopt++; 00400 continue; 00401 } 00402 00403 00404 /*----- -polort num -----*/ 00405 if (strcmp(argv[nopt], "-polort") == 0) 00406 { 00407 nopt++; 00408 if (nopt >= argc) DC_error ("need argument after -polort "); 00409 sscanf (argv[nopt], "%d", &ival); 00410 if (ival < -1) 00411 DC_error ("illegal argument after -polort "); 00412 option_data->polort = ival; 00413 nopt++; 00414 continue; 00415 } 00416 00417 00418 /*----- -base_file filename -----*/ 00419 if (strcmp(argv[nopt], "-base_file") == 0) 00420 { 00421 nopt++; 00422 if (nopt >= argc) DC_error ("need argument after -base_file "); 00423 option_data->base_filename = malloc (sizeof(char)*THD_MAX_NAME); 00424 MTEST (option_data->base_filename); 00425 strcpy (option_data->base_filename, argv[nopt]); 00426 nopt++; 00427 continue; 00428 } 00429 00430 00431 /*----- -num_stimts num -----*/ 00432 if (strcmp(argv[nopt], "-num_stimts") == 0) 00433 { 00434 nopt++; 00435 if (nopt >= argc) DC_error ("need argument after -num_stimts "); 00436 sscanf (argv[nopt], "%d", &ival); 00437 if (ival < 0) 00438 { 00439 DC_error ("-num_stimts num Require: num >= 0 "); 00440 } 00441 00442 initialize_stim_options (option_data, ival); 00443 00444 nopt++; 00445 continue; 00446 } 00447 00448 00449 /*----- -stim_file k sname -----*/ 00450 if (strcmp(argv[nopt], "-stim_file") == 0) 00451 { 00452 nopt++; 00453 if (nopt+1 >= argc) DC_error ("need 2 arguments after -stim_file"); 00454 00455 sscanf (argv[nopt], "%d", &ival); 00456 if ((ival < 1) || (ival > option_data->num_stimts)) 00457 DC_error ("-stim_file k sname Require: 1 <= k <= num_stimts"); 00458 k = ival-1; 00459 nopt++; 00460 00461 option_data->stim_filename[k] 00462 = malloc (sizeof(char)*THD_MAX_NAME); 00463 MTEST (option_data->stim_filename[k]); 00464 strcpy (option_data->stim_filename[k], argv[nopt]); 00465 nopt++; 00466 continue; 00467 } 00468 00469 00470 /*----- -stim_minlag k lag -----*/ 00471 if (strcmp(argv[nopt], "-stim_minlag") == 0) 00472 { 00473 nopt++; 00474 if (nopt+1 >= argc) 00475 DC_error ("need 2 arguments after -stim_minlag"); 00476 00477 sscanf (argv[nopt], "%d", &ival); 00478 if ((ival < 1) || (ival > option_data->num_stimts)) 00479 DC_error ("-stim_minlag k lag Require: 1 <= k <= num_stimts"); 00480 k = ival-1; 00481 nopt++; 00482 00483 sscanf (argv[nopt], "%d", &ival); 00484 if (ival < 0) 00485 DC_error ("-stim_minlag k lag Require: 0 <= lag"); 00486 option_data->stim_minlag[k] = ival; 00487 nopt++; 00488 continue; 00489 } 00490 00491 00492 /*----- -stim_maxlag k lag -----*/ 00493 if (strcmp(argv[nopt], "-stim_maxlag") == 0) 00494 { 00495 nopt++; 00496 if (nopt+1 >= argc) 00497 DC_error ("need 2 arguments after -stim_maxlag"); 00498 00499 sscanf (argv[nopt], "%d", &ival); 00500 if ((ival < 1) || (ival > option_data->num_stimts)) 00501 DC_error ("-stim_maxlag k lag Require: 1 <= k <= num_stimts"); 00502 k = ival-1; 00503 nopt++; 00504 00505 sscanf (argv[nopt], "%d", &ival); 00506 if (ival < 0) 00507 DC_error ("-stim_maxlag k lag Require: 0 <= lag"); 00508 option_data->stim_maxlag[k] = ival; 00509 nopt++; 00510 continue; 00511 } 00512 00513 00514 /*----- -stim_nptr k p -----*/ 00515 if (strcmp(argv[nopt], "-stim_nptr") == 0) 00516 { 00517 nopt++; 00518 if (nopt+1 >= argc) 00519 DC_error ("need 2 arguments after -stim_nptr"); 00520 00521 sscanf (argv[nopt], "%d", &ival); 00522 if ((ival < 1) || (ival > option_data->num_stimts)) 00523 DC_error ("-stim_nptr k p Require: 1 <= k <= num_stimts"); 00524 k = ival-1; 00525 nopt++; 00526 00527 sscanf (argv[nopt], "%d", &ival); 00528 if (ival < 1) 00529 DC_error ("-stim_nptr k p Require: 1 <= p"); 00530 option_data->stim_nptr[k] = ival; 00531 nopt++; 00532 continue; 00533 } 00534 00535 00536 /*----- -iresp k iprefix -----*/ 00537 if (strcmp(argv[nopt], "-iresp") == 0) 00538 { 00539 nopt++; 00540 if (nopt+1 >= argc) DC_error ("need 2 arguments after -iresp"); 00541 00542 sscanf (argv[nopt], "%d", &ival); 00543 if ((ival < 1) || (ival > option_data->num_stimts)) 00544 DC_error ("-iresp k iprefix Require: 1 <= k <= num_stimts"); 00545 k = ival-1; 00546 nopt++; 00547 00548 option_data->iresp_filename[k] 00549 = malloc (sizeof(char)*THD_MAX_NAME); 00550 MTEST (option_data->iresp_filename[k]); 00551 strcpy (option_data->iresp_filename[k], argv[nopt]); 00552 nopt++; 00553 continue; 00554 } 00555 00556 00557 /*----- -errts filename -----*/ 00558 if (strcmp(argv[nopt], "-errts") == 0) 00559 { 00560 nopt++; 00561 if (nopt >= argc) DC_error ("need file prefixname after -errts "); 00562 option_data->errts_filename = malloc (sizeof(char)*THD_MAX_NAME); 00563 MTEST (option_data->errts_filename); 00564 strcpy (option_data->errts_filename, argv[nopt]); 00565 nopt++; 00566 continue; 00567 } 00568 00569 00570 /*----- -sigma s -----*/ 00571 if (strcmp(argv[nopt], "-sigma") == 0) 00572 { 00573 nopt++; 00574 if (nopt >= argc) DC_error ("need argument after -sigma "); 00575 sscanf (argv[nopt], "%f", &fval); 00576 if (fval < 0.0) 00577 DC_error ("illegal argument after -sigma "); 00578 option_data->sigma = fval; 00579 nopt++; 00580 continue; 00581 } 00582 00583 00584 /*----- -seed s -----*/ 00585 if (strcmp(argv[nopt], "-seed") == 0) 00586 { 00587 nopt++; 00588 if (nopt >= argc) DC_error ("need argument after -seed "); 00589 sscanf (argv[nopt], "%ld", &lval); 00590 if (lval <= 0) 00591 DC_error ("illegal argument after -seed "); 00592 option_data->seed = lval; 00593 nopt++; 00594 continue; 00595 } 00596 00597 00598 /*----- -xout -----*/ 00599 if (strcmp(argv[nopt], "-xout") == 0) 00600 { 00601 option_data->xout = 1; 00602 nopt++; 00603 continue; 00604 } 00605 00606 00607 /*----- -output filename -----*/ 00608 if (strcmp(argv[nopt], "-output") == 0) 00609 { 00610 nopt++; 00611 if (nopt >= argc) DC_error ("need argument after -output "); 00612 option_data->output_filename = 00613 malloc (sizeof(char)*THD_MAX_NAME); 00614 MTEST (option_data->output_filename); 00615 strcpy (option_data->output_filename, argv[nopt]); 00616 nopt++; 00617 continue; 00618 } 00619 00620 00621 /*----- unknown command -----*/ 00622 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]); 00623 DC_error (message); 00624 00625 } 00626 00627 00628 } |
|
Definition at line 175 of file 3dConvolve.c.
00179 { 00180 int is; /* input stimulus time series index */ 00181 00182 00183 00184 /*----- Initialize default values -----*/ 00185 option_data->nxyz = -1; 00186 option_data->nt = -1; 00187 option_data->NFirst = -1; 00188 option_data->NLast = -1; 00189 option_data->N = 0; 00190 option_data->polort = 1; 00191 option_data->p = 0; 00192 option_data->q = 0; 00193 option_data->input1D = 0; 00194 option_data->sigma = 0.0; 00195 option_data->seed = 1234567; 00196 option_data->xout = 0; 00197 00198 00199 /*----- Initialize stimulus options -----*/ 00200 option_data->num_stimts = 0; 00201 option_data->stim_filename = NULL; 00202 option_data->stim_minlag = NULL; 00203 option_data->stim_maxlag = NULL; 00204 option_data->stim_nptr = NULL; 00205 option_data->iresp_filename = NULL; 00206 00207 00208 /*----- Initialize character strings -----*/ 00209 option_data->input_filename = NULL; 00210 option_data->mask_filename = NULL; 00211 option_data->base_filename = NULL; 00212 option_data->censor_filename = NULL; 00213 option_data->concat_filename = NULL; 00214 option_data->errts_filename = NULL; 00215 option_data->output_filename = NULL; 00216 00217 00218 } |
|
Definition at line 1317 of file 3dConvolve.c. References allocate_memory(), argc, check_for_valid_inputs(), get_options(), malloc, and read_input_data().
01342 { 01343 01344 01345 /*----- Allocate memory -----*/ 01346 *option_data = (DC_options *) malloc (sizeof(DC_options)); 01347 01348 01349 /*----- Get command line inputs -----*/ 01350 get_options (argc, argv, *option_data); 01351 01352 01353 /*----- Read input data -----*/ 01354 read_input_data (*option_data, dset_time, mask_dset, base_dset, err_dset, 01355 irf_dset, base_data, base_length, irf_data, irf_length, 01356 censor_array, censor_length, block_list, num_blocks, 01357 stimulus, stim_length, errts_data, errts_length); 01358 01359 /*----- Check for valid inputs -----*/ 01360 check_for_valid_inputs (*option_data, *dset_time, *mask_dset, *base_length, 01361 *irf_length, *censor_array, *censor_length, 01362 *block_list, *num_blocks, *stim_length, 01363 *errts_length, good_list); 01364 01365 01366 /*----- Allocate memory for output volumes -----*/ 01367 allocate_memory (*option_data, predts_vol); 01368 01369 01370 /*----- Initialize random number generator -----*/ 01371 srand48 ((*option_data)->seed); 01372 01373 } |
|
Definition at line 227 of file 3dConvolve.c.
00232 { 00233 int is; /* input stimulus time series index */ 00234 00235 00236 /*----- Set number of input stimulus time series -----*/ 00237 if (num_stimts <= 0) return; 00238 else option_data->num_stimts = num_stimts; 00239 00240 00241 /*----- Allocate memory for stimulus options -----*/ 00242 option_data->stim_filename = (char **) malloc (sizeof(char *) * num_stimts); 00243 MTEST (option_data->stim_filename); 00244 option_data->stim_minlag = (int *) malloc (sizeof(int) * num_stimts); 00245 MTEST (option_data->stim_minlag); 00246 option_data->stim_maxlag = (int *) malloc (sizeof(int) * num_stimts); 00247 MTEST (option_data->stim_maxlag); 00248 option_data->stim_nptr = (int *) malloc (sizeof(int) * num_stimts); 00249 MTEST (option_data->stim_nptr); 00250 option_data->iresp_filename = (char **) malloc (sizeof(char *) * num_stimts); 00251 MTEST (option_data->iresp_filename); 00252 00253 00254 /*----- Initialize stimulus options -----*/ 00255 for (is = 0; is < num_stimts; is++) 00256 { 00257 option_data->stim_filename[is] = NULL; 00258 00259 option_data->stim_minlag[is] = 0; 00260 option_data->stim_maxlag[is] = 0; 00261 option_data->stim_nptr[is] = 1; 00262 00263 option_data->iresp_filename[is] = NULL; 00264 } 00265 00266 } |
|
Definition at line 2036 of file 3dConvolve.c. References argc, calculate_results(), free, initialize_program(), DC_options::num_stimts, output_results(), PROGRAM_AUTHOR, PROGRAM_INITIAL, PROGRAM_LATEST, PROGRAM_NAME, terminate_program(), and THD_delete_3dim_dataset().
02041 { 02042 DC_options * option_data; /* deconvolution algorithm options */ 02043 THD_3dim_dataset * dset_time = NULL; /* input 3d+time template data set */ 02044 THD_3dim_dataset * mask_dset = NULL; /* input mask data set */ 02045 THD_3dim_dataset * base_dset = NULL; /* input baseline parameter data set */ 02046 THD_3dim_dataset * err_dset = NULL; /* input error 3d+time data set */ 02047 THD_3dim_dataset ** irf_dset = NULL; /* IRF 3d+time dataset array */ 02048 float * base_data = NULL; /* input baseline parameters */ 02049 int base_length; /* number of input baseline parameters */ 02050 float ** irf_data = NULL; /* input IRF time series data */ 02051 int * irf_length = NULL; /* length of IRF time series */ 02052 float * censor_array = NULL; /* input censor time series array */ 02053 int censor_length; /* length of censor time series */ 02054 int * good_list = NULL; /* list of usable time points */ 02055 int * block_list = NULL; /* list of block starting points */ 02056 int num_blocks; /* number of blocks (runs) */ 02057 float ** stimulus = NULL; /* stimulus time series arrays */ 02058 int * stim_length = NULL; /* length of stimulus time series */ 02059 float * errts_data = NULL; /* input error time series data */ 02060 int errts_length; /* length of error time series */ 02061 02062 float ** predts_vol = NULL; /* volumes for estimated time series data */ 02063 int is; /* stimulus index */ 02064 02065 02066 02067 /*----- Identify software -----*/ 02068 printf ("\n\n"); 02069 printf ("Program: %s \n", PROGRAM_NAME); 02070 printf ("Author: %s \n", PROGRAM_AUTHOR); 02071 printf ("Initial Release: %s \n", PROGRAM_INITIAL); 02072 printf ("Latest Revision: %s \n", PROGRAM_LATEST); 02073 printf ("\n"); 02074 02075 02076 /*----- Program initialization -----*/ 02077 initialize_program (argc, argv, 02078 &option_data, &dset_time, &mask_dset, &base_dset, &err_dset, 02079 &irf_dset, &base_data, &base_length, &irf_data, &irf_length, 02080 &censor_array, &censor_length, 02081 &good_list, &block_list, &num_blocks, &stimulus, &stim_length, 02082 &errts_data, &errts_length, &predts_vol); 02083 02084 02085 /*----- Perform convolution -----*/ 02086 calculate_results (option_data, dset_time, mask_dset, base_dset, err_dset, 02087 irf_dset, base_data, base_length, irf_data, irf_length, 02088 good_list, block_list, num_blocks, stimulus, stim_length, 02089 errts_data, errts_length, predts_vol); 02090 02091 02092 /*----- Deallocate memory for input datasets -----*/ 02093 if (dset_time != NULL) 02094 { THD_delete_3dim_dataset (dset_time, False); dset_time = NULL; } 02095 if (mask_dset != NULL) 02096 { THD_delete_3dim_dataset (mask_dset, False); mask_dset = NULL; } 02097 if (base_dset != NULL) 02098 { THD_delete_3dim_dataset (base_dset, False); base_dset = NULL; } 02099 if (err_dset != NULL) 02100 { THD_delete_3dim_dataset (err_dset, False); err_dset = NULL; } 02101 if (irf_dset != NULL) 02102 { 02103 for (is = 0; is < option_data->num_stimts; is++) 02104 if (irf_dset[is] != NULL) 02105 { 02106 THD_delete_3dim_dataset (irf_dset[is], False); 02107 irf_dset[is] = NULL; 02108 } 02109 free(irf_dset); irf_dset = NULL; 02110 } 02111 02112 02113 /*----- Write requested output files -----*/ 02114 output_results (argc, argv, option_data, predts_vol); 02115 02116 02117 /*----- Terminate program -----*/ 02118 terminate_program (&option_data, &base_data, &irf_data, &irf_length, 02119 &censor_array, &good_list, &block_list, 02120 &stimulus, &stim_length, &errts_data, &predts_vol); 02121 02122 exit(0); 02123 } |
|
Definition at line 1911 of file 3dConvolve.c. References argc, write_one_ts(), and write_ts_array().
01919 { 01920 int * nptr; /* number of stim fn. time points per TR */ 01921 int nt; /* number of images in input 3d+time dataset */ 01922 01923 01924 /*----- Initialize local variables -----*/ 01925 nptr = option_data->stim_nptr; 01926 nt = option_data->nt; 01927 01928 01929 /*----- Write the predicted time series data -----*/ 01930 if (option_data->output_filename != NULL) 01931 01932 if (option_data->input1D) 01933 write_one_ts (option_data->output_filename, nt, predts_vol); 01934 else 01935 write_ts_array (argc, argv, option_data, nt, predts_vol, 01936 option_data->output_filename); 01937 01938 01939 } |
|
Definition at line 701 of file 3dConvolve.c. References DC_error(), DSET_NUM_TIMES, DSET_NVALS, DSET_NVOX, ISVALID_3DIM_DATASET, malloc, MTEST, p, q, read_time_series(), THD_load_datablock(), THD_MAX_NAME, THD_open_dataset(), and THD_open_one_dataset().
00722 { 00723 char message[THD_MAX_NAME]; /* error message */ 00724 int it; /* time point index */ 00725 int nt; /* number of input data time points */ 00726 int nxyz; /* number of voxels */ 00727 int num_stimts; /* number of stimulus time series arrays */ 00728 int * min_lag; /* minimum time delay for impulse response */ 00729 int * max_lag; /* maximum time delay for impulse response */ 00730 int q; /* number of baseline parameters */ 00731 int p; /* total number of parameters */ 00732 int is; /* stimulus time series index */ 00733 00734 00735 /*----- Initialize local variables -----*/ 00736 num_stimts = option_data->num_stimts; 00737 min_lag = option_data->stim_minlag; 00738 max_lag = option_data->stim_maxlag; 00739 00740 00741 /*----- Read the block list -----*/ 00742 if (option_data->concat_filename == NULL) 00743 { 00744 *num_blocks = 1; 00745 *block_list = (int *) malloc (sizeof(int) * 1); 00746 (*block_list)[0] = 0; 00747 } 00748 else 00749 { 00750 float * f = NULL; 00751 00752 f = read_time_series (option_data->concat_filename, num_blocks); 00753 if (*num_blocks < 1) 00754 { 00755 sprintf (message, "Problem reading concat file: %s ", 00756 option_data->concat_filename); 00757 DC_error (message); 00758 } 00759 else 00760 { 00761 *block_list = (int *) malloc (sizeof(int) * (*num_blocks)); 00762 for (it = 0; it < *num_blocks; it++) 00763 (*block_list)[it] = floor (f[it]+0.5); 00764 } 00765 } 00766 00767 00768 /*----- Determine total number of parameters in the model -----*/ 00769 q = (option_data->polort + 1) * (*num_blocks); 00770 p = q; 00771 for (is = 0; is < num_stimts; is++) 00772 { 00773 if (max_lag[is] < min_lag[is]) 00774 DC_error ("Require min lag <= max lag for all stimuli"); 00775 p += max_lag[is] - min_lag[is] + 1; 00776 } 00777 option_data->p = p; 00778 option_data->q = q; 00779 00780 00781 /*----- Read the input baseline, IRF, and error time series data -----*/ 00782 00783 00784 /*----- Read .1D files -----*/ 00785 if (option_data->input1D) 00786 { 00787 *dset_time = NULL; 00788 nt = option_data->NLast + 1; 00789 nxyz = 1; 00790 00791 /*----- Read the baseline parameter .1D time series -----*/ 00792 if (option_data->base_filename != NULL) 00793 { 00794 *base_data = read_time_series (option_data->base_filename, 00795 base_length); 00796 if (*base_data == NULL) 00797 { 00798 sprintf (message, "Unable to read baseline .1D file: %s", 00799 option_data->base_filename); 00800 DC_error (message); 00801 } 00802 } 00803 else 00804 { 00805 *base_data = NULL; 00806 *base_length = 0; 00807 } 00808 00809 /*----- Read the input IRF .1D time series -----*/ 00810 if (num_stimts > 0) 00811 { 00812 *irf_data = (float **) malloc (sizeof(float *) * num_stimts); 00813 MTEST (*irf_data); 00814 *irf_length = (int *) malloc (sizeof(int) * num_stimts); 00815 MTEST (*irf_length); 00816 for (is = 0; is < num_stimts; is++) 00817 { 00818 (*irf_data)[is] 00819 = read_time_series (option_data->iresp_filename[is], 00820 &(*irf_length)[is]); 00821 if ((*irf_data)[is] == NULL) 00822 { 00823 sprintf (message, "Unable to read IRF .1D file: %s", 00824 option_data->iresp_filename[is]); 00825 DC_error (message); 00826 } 00827 } 00828 } 00829 00830 /*----- Read the error .1D time series -----*/ 00831 if (option_data->errts_filename != NULL) 00832 { 00833 *errts_data = read_time_series (option_data->errts_filename, 00834 errts_length); 00835 if (*errts_data == NULL) 00836 { 00837 sprintf (message, "Unable to read residual error .1D file: %s", 00838 option_data->errts_filename); 00839 DC_error (message); 00840 } 00841 } 00842 else 00843 { 00844 *errts_data = NULL; 00845 *errts_length = 0; 00846 } 00847 } 00848 00849 00850 /*----- Read 3d+time datasets -----*/ 00851 else if (option_data->input_filename != NULL) 00852 { 00853 /*----- Read the input 3d+time (template) dataset -----*/ 00854 *dset_time = THD_open_one_dataset (option_data->input_filename); 00855 if (!ISVALID_3DIM_DATASET(*dset_time)) 00856 { 00857 sprintf (message, "Unable to open template dataset file: %s", 00858 option_data->input_filename); 00859 DC_error (message); 00860 } 00861 THD_load_datablock ((*dset_time)->dblk); 00862 nt = DSET_NUM_TIMES (*dset_time); 00863 nxyz = DSET_NVOX (*dset_time); 00864 00865 /*----- Read the input mask dataset -----*/ 00866 if (option_data->mask_filename != NULL) 00867 { 00868 *mask_dset = THD_open_dataset (option_data->mask_filename); 00869 if (!ISVALID_3DIM_DATASET(*mask_dset)) 00870 { 00871 sprintf (message, "Unable to open mask file: %s", 00872 option_data->mask_filename); 00873 DC_error (message); 00874 } 00875 THD_load_datablock ((*mask_dset)->dblk); 00876 } 00877 00878 /*----- Read the input baseline parameter dataset -----*/ 00879 if (option_data->base_filename != NULL) 00880 { 00881 *base_dset = THD_open_dataset (option_data->base_filename); 00882 if (!ISVALID_3DIM_DATASET(*base_dset)) 00883 { 00884 sprintf (message, "Unable to open baseline parameter file: %s", 00885 option_data->base_filename); 00886 DC_error (message); 00887 } 00888 THD_load_datablock ((*base_dset)->dblk); 00889 *base_length = DSET_NVALS (*base_dset); 00890 } 00891 else 00892 { 00893 *base_dset = NULL; 00894 *base_length = 0; 00895 } 00896 00897 /*----- Read the IRF 3d+time datasets -----*/ 00898 if (num_stimts > 0) 00899 { 00900 *irf_length = (int *) malloc (sizeof(int) * num_stimts); 00901 MTEST (*irf_length); 00902 00903 *irf_dset = (THD_3dim_dataset **) 00904 malloc (sizeof(THD_3dim_dataset *) * num_stimts); 00905 MTEST (*irf_dset); 00906 00907 for (is = 0; is < num_stimts; is++) 00908 { 00909 (*irf_dset)[is] 00910 = THD_open_dataset (option_data->iresp_filename[is]); 00911 if (!ISVALID_3DIM_DATASET((*irf_dset)[is])) 00912 { 00913 sprintf (message, "Unable to open IRF file: %s", 00914 option_data->iresp_filename[is]); 00915 DC_error (message); 00916 } 00917 THD_load_datablock ((*irf_dset)[is]->dblk); 00918 (*irf_length)[is] = DSET_NVALS ((*irf_dset)[is]); 00919 } 00920 } 00921 00922 /*----- Read the input error time series dataset -----*/ 00923 if (option_data->errts_filename != NULL) 00924 { 00925 *err_dset = THD_open_dataset (option_data->errts_filename); 00926 if (!ISVALID_3DIM_DATASET(*err_dset)) 00927 { 00928 sprintf (message, "Unable to open error time series file: %s", 00929 option_data->errts_filename); 00930 DC_error (message); 00931 } 00932 THD_load_datablock ((*err_dset)->dblk); 00933 *errts_length = DSET_NVALS (*err_dset); 00934 } 00935 else 00936 { 00937 *err_dset = NULL; 00938 *errts_length = 0; 00939 } 00940 } 00941 00942 else 00943 DC_error ("Must specify input 3d+time template dataset"); 00944 00945 00946 /*----- Check number of data points -----*/ 00947 if (nt <= 0) DC_error ("No time points? Please use -nlast option."); 00948 option_data->nt = nt; 00949 if (nxyz < 0) DC_error ("Program initialization error: nxyz < 0"); 00950 option_data->nxyz = nxyz; 00951 00952 00953 /*----- Read the censorship file -----*/ 00954 if (option_data->censor_filename != NULL) 00955 { 00956 /*----- Read the input censor time series array -----*/ 00957 *censor_array = read_time_series (option_data->censor_filename, 00958 censor_length); 00959 if (*censor_array == NULL) 00960 { 00961 sprintf (message, "Unable to read censor time series file: %s", 00962 option_data->censor_filename); 00963 DC_error (message); 00964 } 00965 } 00966 else 00967 { 00968 /*----- Create censor time series array -----*/ 00969 *censor_array = (float *) malloc (sizeof(float) * nt); 00970 MTEST (*censor_array); 00971 *censor_length = nt; 00972 for (it = 0; it < nt; it++) 00973 (*censor_array)[it] = 1.0; 00974 } 00975 00976 00977 /*----- Read the input stimulus time series -----*/ 00978 if (num_stimts > 0) 00979 { 00980 *stimulus = (float **) malloc (sizeof(float *) * num_stimts); 00981 MTEST (*stimulus); 00982 *stim_length = (int *) malloc (sizeof(int) * num_stimts); 00983 MTEST (*stim_length); 00984 00985 for (is = 0; is < num_stimts; is++) 00986 { 00987 (*stimulus)[is] = read_time_series (option_data->stim_filename[is], 00988 &((*stim_length)[is])); 00989 00990 if ((*stimulus)[is] == NULL) 00991 { 00992 sprintf (message, "Unable to read stimulus time series: %s", 00993 option_data->stim_filename[is]); 00994 DC_error (message); 00995 } 00996 } 00997 } 00998 00999 } |
|
--- Test various combinations for legality [11 Aug 2004] ---* Definition at line 638 of file 3dConvolve.c. References DC_error(), far, malloc, MRI_FLOAT_PTR, mri_free(), mri_read_1D(), MTEST, MRI_IMAGE::nx, MRI_IMAGE::ny, and THD_MAX_NAME.
00643 { 00644 char message[THD_MAX_NAME]; /* error message */ 00645 char * cpt; /* pointer to column suffix */ 00646 char filename[THD_MAX_NAME]; /* time series file name w/o column index */ 00647 char subv[THD_MAX_NAME]; /* string containing column index */ 00648 MRI_IMAGE * im, * flim; /* pointers to image structures 00649 -- used to read 1D ASCII */ 00650 float * far; /* pointer to MRI_IMAGE floating point data */ 00651 int nx; /* number of time points in time series */ 00652 int ny; /* number of columns in time series file */ 00653 int iy; /* time series file column index */ 00654 int ipt; /* time point index */ 00655 float * ts_data = NULL; /* input time series data */ 00656 00657 00658 /*----- First, check for empty filename -----*/ 00659 if (ts_filename == NULL) 00660 DC_error ("Missing input time series file name"); 00661 00662 00663 /*----- Read the time series file -----*/ 00664 00665 flim = mri_read_1D(ts_filename) ; 00666 if( flim == NULL ) 00667 { 00668 sprintf (message, "Unable to read time series file: %s", ts_filename); 00669 DC_error (message); 00670 } 00671 00672 00673 /*----- Set pointer to data, and set dimensions -----*/ 00674 far = MRI_FLOAT_PTR(flim); 00675 nx = flim->nx; 00676 ny = flim->ny; iy = 0 ; 00677 if( ny > 1 ){ 00678 fprintf(stderr,"WARNING: time series %s has %d columns\n",ts_filename,ny); 00679 } 00680 00681 /*----- Save the time series data -----*/ 00682 *ts_length = nx; 00683 ts_data = (float *) malloc (sizeof(float) * nx); 00684 MTEST (ts_data); 00685 for (ipt = 0; ipt < nx; ipt++) 00686 ts_data[ipt] = far[ipt + iy*nx]; 00687 00688 00689 mri_free (flim); flim = NULL; 00690 00691 return (ts_data); 00692 } |
|
Definition at line 1487 of file 3dConvolve.c.
01496 { 01497 int it; /* time point index */ 01498 int nt; /* total number of time points */ 01499 01500 01501 if (predts_vol == NULL) return; 01502 01503 01504 /*----- Initialize local variables -----*/ 01505 nt = option_data->nt; 01506 01507 01508 /*----- Save the predicted time series -----*/ 01509 for (it = 0; it < nt; it++) 01510 predts_vol[it][iv] = predts[it]; 01511 01512 } |
|
Definition at line 1945 of file 3dConvolve.c. References free.
01959 { 01960 int num_stimts; /* number of stimulus time series */ 01961 int is; /* stimulus index */ 01962 int it; /* time index */ 01963 int nt; /* number of images in input 3d+time dataset */ 01964 01965 01966 /*----- Initialize local variables -----*/ 01967 num_stimts = (*option_data)->num_stimts; 01968 nt = (*option_data)->nt; 01969 01970 01971 /*----- Deallocate memory for option data -----*/ 01972 if (*option_data != NULL) 01973 { free (*option_data); *option_data = NULL; } 01974 01975 /*----- Deallocate base parameter memory -----*/ 01976 if (*base_data != NULL) 01977 { free (*base_data); *base_data = NULL; } 01978 01979 /*----- Deallocate memory for IRF time series -----*/ 01980 if (*irf_data != NULL) 01981 { 01982 for (is = 0; is < num_stimts; is++) 01983 if ((*irf_data)[is] != NULL) 01984 { free ((*irf_data)[is]); (*irf_data)[is] = NULL; } 01985 free (*irf_data); *irf_data = NULL; 01986 } 01987 01988 /*----- Deallocate memory for IRF length -----*/ 01989 if (*irf_length != NULL) 01990 { free (*irf_length); *irf_length = NULL; } 01991 01992 /*----- Deallocate memory for censor array -----*/ 01993 if (*censor_array != NULL) 01994 { free (*censor_array); *censor_array = NULL; } 01995 01996 /*----- Deallocate memory for list of usable time points -----*/ 01997 if (*good_list != NULL) 01998 { free (*good_list); *good_list = NULL; } 01999 02000 /*----- Deallocate memory for list of block starting points -----*/ 02001 if (*block_list != NULL) 02002 { free (*block_list); *block_list = NULL; } 02003 02004 /*----- Deallocate memory for stimulus time series -----*/ 02005 if (*stimulus != NULL) 02006 { 02007 for (is = 0; is < num_stimts; is++) 02008 if ((*stimulus)[is] != NULL) 02009 { free ((*stimulus)[is]); (*stimulus)[is] = NULL; } 02010 free (*stimulus); *stimulus = NULL; 02011 } 02012 02013 /*----- Deallocate memory for length of stimulus time series -----*/ 02014 if (*stim_length != NULL) 02015 { free (*stim_length); *stim_length = NULL; } 02016 02017 /*----- Deallocate memory for residual error time series -----*/ 02018 if (*errts_data != NULL) 02019 { free (*errts_data); *errts_data = NULL; } 02020 02021 /*----- Deallocate space for predicted time series -----*/ 02022 if (*predts_vol != NULL) 02023 { 02024 for (it = 0; it < nt; it++) 02025 if ((*predts_vol)[it] != NULL) 02026 { free ((*predts_vol)[it]); (*predts_vol)[it] = NULL; } 02027 free (*predts_vol); *predts_vol = NULL; 02028 } 02029 02030 } |
|
Definition at line 1876 of file 3dConvolve.c.
01882 { 01883 char filename[80]; /* output time series file name */ 01884 int it; /* time index */ 01885 FILE * outfile = NULL; /* file pointer */ 01886 01887 01888 /*----- Create output filename by appending ".1D" -----*/ 01889 sprintf (filename, "%s.1D", prefix); 01890 outfile = fopen (filename, "w"); 01891 01892 01893 /*----- 'Volume' data consists of just one voxel -----*/ 01894 for (it = 0; it < ts_length; it++) 01895 { 01896 fprintf (outfile, "%f", vol_array[it][0]); 01897 fprintf (outfile, " \n"); 01898 } 01899 01900 01901 fclose (outfile); 01902 } |
|
Definition at line 1752 of file 3dConvolve.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().
01761 { 01762 const float EPSILON = 1.0e-10; 01763 01764 THD_3dim_dataset * dset = NULL; /* input afni data set pointer */ 01765 THD_3dim_dataset * new_dset = NULL; /* output afni data set pointer */ 01766 int ib; /* sub-brick index */ 01767 int ierror; /* number of errors in editing data */ 01768 int nxyz; /* total number of voxels */ 01769 float factor; /* factor is new scale factor for sub-brick #ib */ 01770 char * input_filename; /* input afni data set file name */ 01771 short ** bar = NULL; /* bar[ib] points to data for sub-brick #ib */ 01772 float * fbuf; /* float buffer */ 01773 float * volume; /* pointer to volume of data */ 01774 char label[80]; /* label for output file */ 01775 01776 01777 /*----- Initialize local variables -----*/ 01778 input_filename = option_data->input_filename; 01779 dset = THD_open_one_dataset (input_filename); 01780 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz; 01781 01782 01783 /*----- allocate memory -----*/ 01784 bar = (short **) malloc (sizeof(short *) * ts_length); MTEST (bar); 01785 fbuf = (float *) malloc (sizeof(float) * ts_length); MTEST (fbuf); 01786 01787 01788 /*-- make an empty copy of the prototype dataset, for eventual output --*/ 01789 new_dset = EDIT_empty_copy (dset); 01790 01791 01792 /*----- Record history of dataset -----*/ 01793 tross_Copy_History( dset , new_dset ) ; 01794 01795 { char * commandline = tross_commandline( PROGRAM_NAME , argc , argv ) ; 01796 sprintf (label, "Output prefix: %s", output_filename); 01797 if( commandline != NULL ) 01798 tross_multi_Append_History( new_dset , commandline,label,NULL ) ; 01799 else 01800 tross_Append_History ( new_dset, label); 01801 free(commandline) ; 01802 } 01803 01804 /*----- Delete prototype dataset -----*/ 01805 THD_delete_3dim_dataset (dset, False); dset = NULL ; 01806 01807 01808 ierror = EDIT_dset_items (new_dset, 01809 ADN_prefix, output_filename, 01810 ADN_label1, output_filename, 01811 ADN_self_name, output_filename, 01812 ADN_malloc_type, DATABLOCK_MEM_MALLOC, 01813 ADN_datum_all, MRI_short, 01814 ADN_nvals, ts_length, 01815 ADN_ntt, ts_length, 01816 ADN_none); 01817 01818 if( ierror > 0 ){ 01819 fprintf(stderr, 01820 "*** %d errors in attempting to create output dataset!\n", ierror ) ; 01821 exit(1) ; 01822 } 01823 01824 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){ 01825 fprintf(stderr, 01826 "*** Output dataset file %s already exists--cannot continue!\a\n", 01827 new_dset->dblk->diskptr->header_name ) ; 01828 exit(1) ; 01829 } 01830 01831 01832 /*----- attach bricks to new data set -----*/ 01833 for (ib = 0; ib < ts_length; ib++) 01834 { 01835 01836 /*----- Set pointer to appropriate volume -----*/ 01837 volume = vol_array[ib]; 01838 01839 /*----- Allocate memory for output sub-brick -----*/ 01840 bar[ib] = (short *) malloc (sizeof(short) * nxyz); 01841 MTEST (bar[ib]); 01842 01843 /*----- Convert data type to short for this sub-brick -----*/ 01844 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume, 01845 MRI_short, bar[ib]); 01846 if (factor < EPSILON) factor = 0.0; 01847 else factor = 1.0 / factor; 01848 fbuf[ib] = factor; 01849 01850 /*----- attach bar[ib] to be sub-brick #ib -----*/ 01851 mri_fix_data_pointer (bar[ib], DSET_BRICK(new_dset,ib)); 01852 } 01853 01854 01855 /*----- write afni data set -----*/ 01856 01857 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ; 01858 01859 THD_load_statistics (new_dset); 01860 THD_write_3dim_dataset (NULL, NULL, new_dset, True); 01861 fprintf (stderr,"-- Output 3D+time dataset into %s\n",DSET_BRIKNAME(new_dset)) ; 01862 01863 /*----- deallocate memory -----*/ 01864 THD_delete_3dim_dataset (new_dset, False); new_dset = NULL ; 01865 free (fbuf); fbuf = NULL; 01866 01867 } |
|
Definition at line 1263 of file 3dConvolve.c.
|