00001
00002
00003
00004
00005
00006
00007 #ifdef USE_SUNPERF
00008 # include <sunperf.h>
00009 #endif
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282 #ifndef FLOATIZE
00283 # define PROGRAM_NAME "3dDeconvolve"
00284 #else
00285 # define PROGRAM_NAME "3dDeconvolve_f"
00286 #endif
00287
00288 #define PROGRAM_AUTHOR "B. Douglas Ward, et al."
00289 #define PROGRAM_INITIAL "02 September 1998"
00290 #define PROGRAM_LATEST "04 March 2005 - RWCox"
00291
00292
00293
00294 #define RA_error DC_error
00295
00296 #define USE_GET
00297
00298
00299
00300 #include "mrilib.h"
00301
00302
00303
00304
00305
00306 #if !defined(DONT_USE_SHM) && !defined(DONT_USE_FORK) && !defined(CYGWIN)
00307
00308 # include "thd_iochan.h"
00309
00310 # define PROC_MAX 32
00311
00312 static int proc_numjob = 1 ;
00313 static pid_t proc_pid[PROC_MAX] ;
00314 static int proc_shmid = 0 ;
00315 static char *proc_shmptr = NULL;
00316 static int proc_shmsize = 0 ;
00317
00318 static int proc_shm_arnum = 0 ;
00319 static float ***proc_shm_ar = NULL;
00320 static int *proc_shm_arsiz = NULL;
00321
00322 static int proc_vox_bot[PROC_MAX] ;
00323 static int proc_vox_top[PROC_MAX] ;
00324
00325 static int proc_ind ;
00326
00327 #else
00328
00329 # define proc_numjob 1
00330 # define proc_ind 0
00331
00332 #endif
00333
00334 static int proc_use_jobs = 0 ;
00335
00336
00337
00338 #ifndef FLOATIZE
00339 # include "matrix.h"
00340 # define MTYPE double
00341 #else
00342 # include "matrix_f.h"
00343 # define MTYPE float
00344 #endif
00345
00346
00347
00348 void JPEG_matrix_gray( matrix X , char *fname ) ;
00349
00350 void XSAVE_output( char * ) ;
00351
00352 static int xsave=0 ;
00353 static int nGoodList=0 , *GoodList=NULL ;
00354 static int nParam=0 , *ParamIndex=NULL , *ParamStim=NULL ;
00355 static char **ParamLabel=NULL ;
00356 static char *InputFilename=NULL, *BucketFilename=NULL, *CoefFilename=NULL ;
00357
00358 static matrix X , XtXinv , XtXinvXt ;
00359
00360 static int *Xcol_inbase ;
00361 static float *Xcol_mean ;
00362 static int voxel_num ;
00363 static float *voxel_base = NULL ;
00364
00365 float baseline_mean( vector coef ) ;
00366
00367 static int xrestore = 0 ;
00368 static char *xrestore_filename = NULL ;
00369 static int NumTimePoints=0 , NumRegressors=0 ;
00370
00371 static int verb = 1 ;
00372
00373 struct DC_options ;
00374
00375 void do_xrestore_stuff( int, char **, struct DC_options * ) ;
00376
00377 #define XSAVE_version "0.5"
00378
00379 static int nSymStim = 0 ;
00380 static SYM_irange *SymStim = NULL ;
00381
00382 void read_glt_matrix( char *fname, int *nrows, int ncol, matrix *cmat ) ;
00383 static void vstep_print(void) ;
00384
00385 static int show_singvals = 0 ;
00386
00387
00388
00389 #include "parser.h"
00390
00391
00392
00393 #define USE_BASIS
00394
00395
00396
00397 typedef struct {
00398 float a , b , c ;
00399 float ffac ;
00400 void *q ;
00401 float (*f)(float,float,float,float,void *) ;
00402 } basis_func ;
00403
00404
00405
00406 #define basis_funceval(bf,x) ((bf).f( (x), (bf).a,(bf).b,(bf).c,(bf).q )*(bf).ffac)
00407
00408
00409
00410 typedef struct {
00411 int nfunc , pbot ;
00412 float tbot,ttop ;
00413 basis_func *bfunc ;
00414 char *name ;
00415 } basis_expansion ;
00416
00417
00418
00419 basis_expansion * basis_parser( char *sym ) ;
00420 float basis_evaluation( basis_expansion *be , float *wt , float x ) ;
00421 void basis_write_iresp( int argc , char *argv[] ,
00422 struct DC_options *option_data ,
00423 basis_expansion *be , float dt ,
00424 float **wtar , char *output_filename ) ;
00425 void basis_write_sresp( int argc , char *argv[] ,
00426 struct DC_options *option_data ,
00427 basis_expansion *be , float dt ,
00428 float *mse ,
00429 int pbot, matrix cvar, char *output_filename ) ;
00430
00431
00432
00433 static basis_expansion **basis_stim = NULL ;
00434 static MRI_IMAGE **basis_times = NULL ;
00435 static MRI_IMAGE **basis_vect = NULL ;
00436 static float basis_TR = 1.0f ;
00437 static int basis_count = 0 ;
00438 static float basis_dtout = 0.0f ;
00439 static float irc_dt = 0.0f ;
00440
00441 static int basis_need_mse = 0 ;
00442
00443 static float basis_normall = 0.0f ;
00444
00445 #define basis_filler 3.e+33
00446
00447
00448
00449 typedef struct {
00450 int npar , pbot ;
00451 float *ww ;
00452 float scale_fac ;
00453 int denom_flag ;
00454 char *name ;
00455 } basis_irc ;
00456
00457 typedef struct { float a,b ; } floatpair ;
00458
00459 #define denom_BASELINE (1)
00460
00461 static int num_irc = 0 ;
00462 static basis_irc **irc = NULL ;
00463
00464 floatpair evaluate_irc( basis_irc *birc , vector coef ,
00465 float base , float mse , matrix cvar ) ;
00466
00467
00468
00469 #include "Deconvolve.c"
00470
00471
00472
00473 typedef struct DC_options
00474 {
00475 int nxyz;
00476 int nt;
00477 int NFirst;
00478 int NLast;
00479 int N;
00480 int polort;
00481 float rms_min;
00482 int quiet;
00483 int progress;
00484 float fdisp;
00485 char * input_filename;
00486 char * mask_filename;
00487 char * input1D_filename;
00488 char * censor_filename;
00489 char * concat_filename;
00490 int nodata;
00491 int p;
00492 int q;
00493 int qp;
00494
00495 int nbricks;
00496
00497 int num_stimts;
00498 char ** stim_filename;
00499 char ** stim_label;
00500 int * stim_base;
00501 int * stim_minlag;
00502 int * stim_maxlag;
00503 int * stim_nptr;
00504
00505 int num_glt;
00506 int * glt_rows;
00507 char ** glt_filename;
00508 char ** glt_label;
00509
00510 char * bucket_filename;
00511 char ** iresp_filename;
00512 char ** sresp_filename;
00513 char * fitts_filename;
00514 char * errts_filename;
00515
00516 int tshift;
00517 int fout;
00518 int rout;
00519 int tout;
00520 int vout;
00521 int nobout;
00522 int nocout;
00523 int xout;
00524 int full_first;
00525
00526 int nocond ;
00527
00528 char *xjpeg_filename;
00529
00530 int automask ;
00531
00532 int nodata_NT ;
00533 float nodata_TR ;
00534 } DC_options;
00535
00536
00537
00538
00539
00540
00541
00542 void DC_error (char * message)
00543 {
00544 fprintf (stderr, "%s Error: %s \a\n\n", PROGRAM_NAME, message);
00545 exit(1);
00546 }
00547
00548
00549
00550
00551 void identify_software ()
00552 {
00553
00554
00555 #if 1
00556 PRINT_VERSION("3dDeconvolve") ;
00557 #else
00558 printf ("\n\n");
00559 printf ("Program: %s \n", PROGRAM_NAME);
00560 printf ("Author: %s \n", PROGRAM_AUTHOR);
00561 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
00562 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
00563 #endif
00564
00565 #ifdef USE_ALTIVEC
00566 INFO_message ("Compiled with Altivec acceleration for Mac OS X\n") ;
00567 #elif defined(USE_SUNPERF) && !defined(FLOATIZE)
00568 INFO_message ("Compiled with BLAS-1 acceleration for Solaris\n") ;
00569 #elif defined(USE_SCSLBLAS)
00570 INFO_message ("Compiled with BLAS-1 acceleration for SGI Altix\n") ;
00571 #endif
00572 }
00573
00574
00575
00576
00577
00578
00579 void display_help_menu()
00580 {
00581 identify_software();
00582
00583 printf (
00584 "Program to calculate the deconvolution of a measurement 3d+time dataset \n"
00585 "with a specified input stimulus time series. This program will also \n"
00586 "perform multiple linear regression using multiple input stimulus time \n"
00587 "series. Output consists of an AFNI 'bucket' type dataset containing the \n"
00588 "least squares estimates of the linear regression coefficients, t-statistics\n"
00589 "for significance of the coefficients, partial F-statistics for significance\n"
00590 "of the individual input stimuli, and the F-statistic for significance of \n"
00591 "the overall regression. Additional output consists of a 3d+time dataset \n"
00592 "containing the estimated system impulse response function. \n"
00593 " \n"
00594 "Usage: \n"
00595 PROGRAM_NAME "\n"
00596 " \n"
00597 "**** Input data and control options: \n"
00598 "-input fname fname = filename of 3d+time input dataset \n"
00599 " (more than one filename can be given) \n"
00600 " (here, and these datasets will be) \n"
00601 " (catenated in time; if you do this, ) \n"
00602 " ('-concat' is not needed and is ignored) \n"
00603 "[-input1D dname] dname = filename of single (fMRI) .1D time series \n"
00604 "[-nodata [NT [TR]] Evaluate experimental design only (no input data) \n"
00605 "[-mask mname] mname = filename of 3d mask dataset \n"
00606 "[-automask] build a mask automatically from input data \n"
00607 " (will be slow for long time series datasets) \n"
00608 "[-censor cname] cname = filename of censor .1D time series \n"
00609 "[-concat rname] rname = filename for list of concatenated runs \n"
00610 "[-nfirst fnum] fnum = number of first dataset image to use in the\n"
00611 " deconvolution procedure. (default = max maxlag) \n"
00612 "[-nlast lnum] lnum = number of last dataset image to use in the \n"
00613 " deconvolution procedure. (default = last point) \n"
00614 "[-polort pnum] pnum = degree of polynomial corresponding to the \n"
00615 " null hypothesis (default: pnum = 1) \n"
00616 "[-legendre] use Legendre polynomials for null hypothesis \n"
00617 "[-nolegendre] use power polynomials for null hypotheses \n"
00618 " (default is -legendre) \n"
00619 "[-nodmbase] don't de-mean baseline time series \n"
00620 " (i.e., polort>1 and -stim_base inputs) \n"
00621 "[-dmbase] de-mean baseline time series (default if polort>0)\n"
00622 "[-nocond] don't calculate matrix condition number \n"
00623 "[-svd] Use SVD instead of Gaussian elimination (default) \n"
00624 "[-nosvd] Use Gaussian elimination instead of SVD \n"
00625 "[-rmsmin r] r = minimum rms error to reject reduced model \n"
00626 " \n"
00627 "**** Input stimulus options: \n"
00628 "-num_stimts num num = number of input stimulus time series \n"
00629 " (0 <= num) (default: num = 0) \n"
00630 "-stim_file k sname sname = filename of kth time series input stimulus\n"
00631 "[-stim_label k slabel] slabel = label for kth input stimulus \n"
00632 "[-stim_base k] kth input stimulus is part of the baseline model \n"
00633 "[-stim_minlag k m] m = minimum time lag for kth input stimulus \n"
00634 " (default: m = 0) \n"
00635 "[-stim_maxlag k n] n = maximum time lag for kth input stimulus \n"
00636 " (default: n = 0) \n"
00637 "[-stim_nptr k p] p = number of stimulus function points per TR \n"
00638 " Note: This option requires 0 slice offset times \n"
00639 " (default: p = 1) \n"
00640 " \n"
00641 "[-stim_times k tname Rmodel] \n"
00642 " Generate the k-th response model from a set of stimulus times \n"
00643 " given in file 'tname'. The response model is specified by the \n"
00644 " 'Rmodel' argument, which can be one of \n"
00645 " 'GAM(p,q)' = 1 parameter gamma variate \n"
00646 " 'SPMG' = 2 parameter SPM gamma variate + derivative \n"
00647 " 'POLY(b,c,n)' = n parameter polynomial expansion \n"
00648 " 'SIN(b,c,n)' = n parameter sine series expansion \n"
00649 " 'TENT(b,c,n)' = n parameter tent function expansion \n"
00650 " 'BLOCK(d,p)' = 1 parameter block stimulus of duration 'd' \n"
00651 " (can also be called 'IGFUN' which stands) \n"
00652 " (for 'incomplete gamma function' ) \n"
00653 " 'EXPR(b,c) exp1 ... expn' = n parameter; arbitrary expressions \n"
00654 " \n"
00655 "[-basis_normall a] \n"
00656 " Normalize all basis functions for '-stim_times' to have \n"
00657 " amplitude 'a' (must have a > 0). The peak absolute value \n"
00658 " of each basis function will be scaled to be 'a'. \n"
00659 " NOTE: -basis_normall only affect -stim_times options that \n"
00660 " appear LATER on the command line \n"
00661 " \n"
00662 "**** General linear test (GLT) options: \n"
00663 "-num_glt num num = number of general linear tests (GLTs) \n"
00664 " (0 <= num) (default: num = 0) \n"
00665 "[-glt s gltname] Perform s simultaneous linear tests, as specified \n"
00666 " by the matrix contained in file gltname \n"
00667 "[-glt_label k glabel] glabel = label for kth general linear test \n"
00668 "[-gltsym gltname] Read the GLT with symbolic names from the file \n"
00669 " \n"
00670 "[-TR_irc dt] \n"
00671 " Use 'dt' as the stepsize for computation of integrals in -IRC_times \n"
00672 " options. Default is to use value given in '-TR_times'. \n"
00673 " \n"
00674 "**** Options for output 3d+time datasets: \n"
00675 "[-iresp k iprefix] iprefix = prefix of 3d+time output dataset which \n"
00676 " will contain the kth estimated impulse response \n"
00677 "[-tshift] Use cubic spline interpolation to time shift the \n"
00678 " estimated impulse response function, in order to\n"
00679 " correct for differences in slice acquisition \n"
00680 " times. Note that this effects only the 3d+time \n"
00681 " output dataset generated by the -iresp option. \n"
00682 "[-sresp k sprefix] sprefix = prefix of 3d+time output dataset which \n"
00683 " will contain the standard deviations of the \n"
00684 " kth impulse response function parameters \n"
00685 "[-fitts fprefix] fprefix = prefix of 3d+time output dataset which \n"
00686 " will contain the (full model) time series fit \n"
00687 " to the input data \n"
00688 "[-errts eprefix] eprefix = prefix of 3d+time output dataset which \n"
00689 " will contain the residual error time series \n"
00690 " from the full model fit to the input data \n"
00691 "[-TR_times dt] \n"
00692 " Use 'dt' as the stepsize for output of -iresp and -sresp file \n"
00693 " for response models generated by '-stim_times' options. \n"
00694 " Default is same as time spacing in the '-input' 3D+time dataset. \n"
00695 " \n"
00696 "**** Options to control the contents of the output bucket dataset: \n"
00697 "[-fout] Flag to output the F-statistics \n"
00698 "[-rout] Flag to output the R^2 statistics \n"
00699 "[-tout] Flag to output the t-statistics \n"
00700 "[-vout] Flag to output the sample variance (MSE) map \n"
00701 "[-nobout] Flag to suppress output of baseline coefficients \n"
00702 " (and associated statistics) \n"
00703 "[-nocout] Flag to suppress output of regression coefficients \n"
00704 " (and associated statistics) \n"
00705 "[-full_first] Flag to specify that the full model statistics will \n"
00706 " appear first in the bucket dataset output \n"
00707 "[-bucket bprefix] Create one AFNI 'bucket' dataset containing various \n"
00708 " parameters of interest, such as the estimated IRF \n"
00709 " coefficients, and full model fit statistics. \n"
00710 " Output 'bucket' dataset is written to bprefix. \n"
00711 " \n"
00712 "[-xsave] Flag to save X matrix into file bprefix.xsave \n"
00713 " (only works if -bucket option is also given) \n"
00714 "[-noxsave] Don't save X matrix (this is the default) \n"
00715 "[-cbucket cprefix] Save the regression coefficients (no statistics) \n"
00716 " into a dataset named 'cprefix'. This dataset \n"
00717 " will be used in a -xrestore run instead of the \n"
00718 " bucket dataset, if possible. \n"
00719 " \n"
00720 "[-xrestore f.xsave] Restore the X matrix, etc. from a previous run \n"
00721 " that was saved into file 'f.xsave'. You can \n"
00722 " then carry out new -glt tests. When -xrestore \n"
00723 " is used, most other command line options are \n"
00724 " ignored. \n"
00725 " \n"
00726 "**** The following options control the screen output only: \n"
00727 "[-quiet] Flag to suppress most screen output \n"
00728 "[-xout] Flag to write X and inv(X'X) matrices to screen \n"
00729 "[-xjpeg filename] Write a JPEG file graphing the X matrix \n"
00730 "[-progress n] Write statistical results for every nth voxel \n"
00731 "[-fdisp fval] Write statistical results for those voxels \n"
00732 " whose full model F-statistic is > fval \n"
00733 );
00734
00735 #ifdef PROC_MAX
00736 printf( "\n"
00737 " -jobs J Run the program with 'J' jobs (sub-processes).\n"
00738 " On a multi-CPU machine, this can speed the\n"
00739 " program up considerably. On a single CPU\n"
00740 " machine, using this option is silly.\n"
00741 " J should be a number from 1 up to the\n"
00742 " number of CPU sharing memory on the system.\n"
00743 " J=1 is normal (single process) operation.\n"
00744 " The maximum allowed value of J is %d.\n"
00745 " * For more information on parallelizing, see\n"
00746 " http://afni.nimh.nih.gov/afni/doc/misc/parallize.html \n"
00747 " * Use -mask to get more speed; cf. 3dAutomask.\n"
00748 , PROC_MAX ) ;
00749 #endif
00750
00751 printf("\n"
00752 "** NOTE **\n"
00753 "This version of the program has been compiled to use\n"
00754 #ifdef FLOATIZE
00755 "single precision arithmetic for most internal calculations.\n"
00756 #else
00757 "double precision arithmetic for most internal calculations.\n"
00758 #endif
00759 ) ;
00760
00761 exit(0);
00762 }
00763
00764
00765
00766
00767
00768
00769
00770 void initialize_options
00771 (
00772 DC_options * option_data
00773 )
00774
00775 {
00776 int is;
00777
00778
00779 option_data->nxyz = -1;
00780 option_data->nt = -1;
00781 option_data->NFirst = -1;
00782 option_data->NLast = -1;
00783 option_data->N = 0;
00784 option_data->polort = 1;
00785 option_data->rms_min = 0.0;
00786 option_data->quiet = 0;
00787 option_data->progress = 0;
00788 option_data->fdisp = -1.0;
00789 option_data->nodata = 0;
00790 option_data->p = 0;
00791 option_data->q = 0;
00792 option_data->qp = 0;
00793 option_data->nbricks = 0;
00794 option_data->nocond = 0;
00795 option_data->nodata_NT= 0;
00796 option_data->nodata_TR= 0.0;
00797
00798 option_data->xjpeg_filename = NULL ;
00799
00800
00801 option_data->num_stimts = 0;
00802 option_data->stim_filename = NULL;
00803 option_data->stim_label = NULL;
00804 option_data->stim_base = NULL;
00805 option_data->stim_minlag = NULL;
00806 option_data->stim_maxlag = NULL;
00807 option_data->stim_nptr = NULL;
00808 option_data->iresp_filename = NULL;
00809 option_data->sresp_filename = NULL;
00810
00811
00812 option_data->num_glt = 0;
00813 option_data->glt_filename = NULL;
00814 option_data->glt_label = NULL;
00815 option_data->glt_rows = NULL;
00816
00817
00818 option_data->tshift = 0;
00819 option_data->fout = 0;
00820 option_data->rout = 0;
00821 option_data->tout = 0;
00822 option_data->vout = 0;
00823 option_data->xout = 0;
00824 option_data->nobout = 0;
00825 option_data->nocout = 0;
00826 option_data->full_first = 0;
00827
00828
00829 option_data->input_filename = NULL;
00830 option_data->mask_filename = NULL;
00831 option_data->input1D_filename = NULL;
00832 option_data->censor_filename = NULL;
00833 option_data->concat_filename = NULL;
00834 option_data->bucket_filename = NULL;
00835 option_data->fitts_filename = NULL;
00836 option_data->errts_filename = NULL;
00837
00838 option_data->automask = 0 ;
00839 }
00840
00841
00842
00843
00844
00845
00846
00847 void initialize_stim_options
00848 (
00849 DC_options * option_data,
00850 int num_stimts
00851 )
00852
00853 {
00854 int is;
00855
00856 ENTRY("initialize_stim_options") ;
00857
00858
00859 if (num_stimts <= 0) EXRETURN ;
00860 else option_data->num_stimts = num_stimts;
00861
00862
00863
00864 option_data->stim_filename = (char **) malloc (sizeof(char *) * num_stimts);
00865 MTEST (option_data->stim_filename);
00866 option_data->stim_label = (char **) malloc (sizeof(char *) * num_stimts);
00867 MTEST (option_data->stim_label);
00868 option_data->stim_base = (int *) malloc (sizeof(int) * num_stimts);
00869 MTEST (option_data->stim_base);
00870 option_data->stim_minlag = (int *) malloc (sizeof(int) * num_stimts);
00871 MTEST (option_data->stim_minlag);
00872 option_data->stim_maxlag = (int *) malloc (sizeof(int) * num_stimts);
00873 MTEST (option_data->stim_maxlag);
00874 option_data->stim_nptr = (int *) malloc (sizeof(int) * num_stimts);
00875 MTEST (option_data->stim_nptr);
00876 option_data->iresp_filename = (char **) malloc (sizeof(char *) * num_stimts);
00877 MTEST (option_data->iresp_filename);
00878 option_data->sresp_filename = (char **) malloc (sizeof(char *) * num_stimts);
00879 MTEST (option_data->sresp_filename);
00880
00881
00882
00883 basis_stim = (basis_expansion **)malloc(sizeof(basis_expansion *)*num_stimts);
00884 basis_times = (MRI_IMAGE **) malloc(sizeof(MRI_IMAGE *) *num_stimts);
00885 basis_vect = (MRI_IMAGE **) malloc(sizeof(MRI_IMAGE *) *num_stimts);
00886
00887
00888 for (is = 0; is < num_stimts; is++)
00889 {
00890 option_data->stim_filename[is] = NULL;
00891 option_data->stim_label[is] = malloc (sizeof(char)*THD_MAX_NAME);
00892 MTEST (option_data->stim_label[is]);
00893 sprintf (option_data->stim_label[is], "Stim#%d", is+1);
00894
00895 option_data->stim_base[is] = 0;
00896 option_data->stim_minlag[is] = 0;
00897 option_data->stim_maxlag[is] = 0;
00898 option_data->stim_nptr[is] = 1;
00899
00900 option_data->iresp_filename[is] = NULL;
00901 option_data->sresp_filename[is] = NULL;
00902
00903 basis_stim [is] = NULL ;
00904 basis_times[is] = NULL ;
00905 basis_vect [is] = NULL ;
00906 }
00907
00908 EXRETURN ;
00909 }
00910
00911
00912
00913
00914
00915
00916
00917 void initialize_glt_options
00918 (
00919 DC_options * option_data,
00920 int num_glt
00921 )
00922
00923 {
00924 int iglt;
00925
00926 ENTRY("initialize_glt_options") ;
00927
00928
00929
00930 if (num_glt <= 0) EXRETURN ;
00931 else option_data->num_glt = num_glt;
00932
00933
00934
00935 option_data->glt_filename = (char **) malloc (sizeof(char *) * num_glt);
00936 MTEST (option_data->glt_filename);
00937 option_data->glt_label = (char **) malloc (sizeof(char *) * num_glt);
00938 MTEST (option_data->glt_label);
00939 option_data->glt_rows = (int *) malloc (sizeof(int) * num_glt);
00940 MTEST (option_data->glt_rows);
00941
00942
00943
00944 for (iglt = 0; iglt < num_glt; iglt++)
00945 {
00946 option_data->glt_filename[iglt] = NULL;
00947 option_data->glt_label[iglt] = malloc (sizeof(char)*THD_MAX_NAME);
00948 MTEST (option_data->glt_label[iglt]);
00949 sprintf (option_data->glt_label[iglt], "GLT #%d ", iglt+1);
00950 option_data->glt_rows[iglt] = 0;
00951 }
00952
00953
00954 EXRETURN ;
00955 }
00956
00957
00958
00959
00960
00961
00962
00963 void get_options
00964 (
00965 int argc,
00966 char ** argv,
00967 DC_options * option_data
00968 )
00969
00970 {
00971 int nopt = 1;
00972 int ival, index;
00973 float fval;
00974 char message[THD_MAX_NAME];
00975 int k;
00976 int s;
00977 int iglt = 0;
00978 int nerr ;
00979
00980
00981 mainENTRY("3dDeconvolve"); machdep() ;
00982 { int new_argc ; char ** new_argv ;
00983 addto_args( argc , argv , &new_argc , &new_argv ) ;
00984 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00985 }
00986
00987
00988
00989 if (argc < 2 || strcmp(argv[1], "-help") == 0) display_help_menu();
00990
00991
00992
00993 AFNI_logger (PROGRAM_NAME,argc,argv);
00994
00995
00996
00997 initialize_options (option_data);
00998
00999
01000
01001 while (nopt < argc )
01002 {
01003
01004 if( strcmp(argv[nopt],"-OK") == 0 ){ nopt++; continue; }
01005
01006
01007
01008 if( strcmp(argv[nopt],"-nocond") == 0 ){
01009 #ifndef FLOATIZE
01010 option_data->nocond = 1 ;
01011 #else
01012 fprintf(stderr,"** WARNING: -nocond is ignored in 3dDeconvolve_f!\n") ;
01013 #endif
01014 nopt++ ; continue ;
01015 }
01016
01017 if( strcmp(argv[nopt],"-singvals") == 0 ){
01018 show_singvals = 1 ; option_data->nocond = 0 ;
01019 nopt++ ; continue ;
01020 }
01021
01022
01023 if (strcmp(argv[nopt], "-xjpeg") == 0)
01024 {
01025 nopt++;
01026 if (nopt >= argc) DC_error ("need argument after -xjpeg ");
01027 option_data->xjpeg_filename = malloc (sizeof(char)*THD_MAX_NAME);
01028 MTEST (option_data->xjpeg_filename);
01029 strcpy (option_data->xjpeg_filename, argv[nopt]);
01030 if( strstr(option_data->xjpeg_filename,".jpg") == NULL &&
01031 strstr(option_data->xjpeg_filename,".JPG") == NULL )
01032 strcat( option_data->xjpeg_filename , ".jpg" ) ;
01033 nopt++; continue;
01034 }
01035
01036
01037
01038 if (strcmp(argv[nopt], "-input") == 0)
01039 {
01040 int iopt , slen ;
01041 nopt++;
01042 if (nopt >= argc) DC_error ("need argument after -input ");
01043 #if 0
01044 option_data->input_filename = malloc (sizeof(char)*THD_MAX_NAME);
01045 MTEST (option_data->input_filename);
01046 strcpy (option_data->input_filename, argv[nopt]);
01047 nopt++;
01048 #else
01049 slen = 0 ;
01050 for( iopt=nopt ; iopt < argc && argv[iopt][0] != '-' ; iopt++ )
01051 slen += strlen(argv[iopt])+8 ;
01052 option_data->input_filename = calloc(sizeof(char),slen) ;
01053 MTEST (option_data->input_filename);
01054 for( iopt=nopt ; iopt < argc && argv[iopt][0] != '-' ; iopt++ ){
01055 strcat( option_data->input_filename , argv[iopt] ) ;
01056 strcat( option_data->input_filename , " " ) ;
01057 }
01058 slen = strlen(option_data->input_filename) ;
01059 option_data->input_filename[slen-1] = '\0' ;
01060 nopt = iopt ;
01061 #endif
01062 continue;
01063 }
01064
01065
01066
01067 if (strcmp(argv[nopt], "-mask") == 0)
01068 {
01069 nopt++;
01070 if (nopt >= argc) DC_error ("need argument after -mask ");
01071 if( option_data->automask ) DC_error("can't use -mask AND -automask!") ;
01072 option_data->mask_filename = malloc (sizeof(char)*THD_MAX_NAME);
01073 MTEST (option_data->mask_filename);
01074 strcpy (option_data->mask_filename, argv[nopt]);
01075 nopt++;
01076 continue;
01077 }
01078
01079
01080 if( strcmp(argv[nopt],"-automask") == 0 ){
01081 if( option_data->mask_filename != NULL )
01082 DC_error("can't use -automask AND -mask!") ;
01083 option_data->automask = 1 ;
01084 nopt++ ; continue ;
01085 }
01086
01087
01088
01089 if (strcmp(argv[nopt], "-input1D") == 0)
01090 {
01091 nopt++;
01092 if (nopt >= argc) DC_error ("need argument after -input1D ");
01093 option_data->input1D_filename =
01094 malloc (sizeof(char)*THD_MAX_NAME);
01095 MTEST (option_data->input1D_filename);
01096 strcpy (option_data->input1D_filename, argv[nopt]);
01097 nopt++;
01098 continue;
01099 }
01100
01101
01102
01103 if (strcmp(argv[nopt], "-censor") == 0)
01104 {
01105 nopt++;
01106 if (nopt >= argc) DC_error ("need argument after -censor ");
01107 option_data->censor_filename =
01108 malloc (sizeof(char)*THD_MAX_NAME);
01109 MTEST (option_data->censor_filename);
01110 strcpy (option_data->censor_filename, argv[nopt]);
01111 nopt++;
01112 continue;
01113 }
01114
01115
01116
01117 if (strcmp(argv[nopt], "-concat") == 0)
01118 {
01119 nopt++;
01120 if (nopt >= argc) DC_error ("need argument after -concat ");
01121 option_data->concat_filename =
01122 malloc (sizeof(char)*THD_MAX_NAME);
01123 MTEST (option_data->concat_filename);
01124 strcpy (option_data->concat_filename, argv[nopt]);
01125 nopt++;
01126 continue;
01127 }
01128
01129
01130
01131 if (strcmp(argv[nopt], "-nodata") == 0)
01132 {
01133 option_data->nodata = 1;
01134 nopt++;
01135
01136
01137
01138 if( isdigit(argv[nopt][0]) ){
01139 option_data->nodata_NT = (int)strtol(argv[nopt++],NULL,10) ;
01140 if( isdigit(argv[nopt][0]) ){
01141 option_data->nodata_TR = (float)strtod(argv[nopt++],NULL) ;
01142 }
01143 }
01144 continue;
01145 }
01146
01147
01148
01149 if (strcmp(argv[nopt], "-nfirst") == 0)
01150 {
01151 nopt++;
01152 if (nopt >= argc) DC_error ("need argument after -nfirst ");
01153 sscanf (argv[nopt], "%d", &ival);
01154 if (ival < 0)
01155 DC_error ("illegal argument after -nfirst ");
01156 option_data->NFirst = ival;
01157 nopt++;
01158 continue;
01159 }
01160
01161
01162
01163 if (strcmp(argv[nopt], "-nlast") == 0)
01164 {
01165 nopt++;
01166 if (nopt >= argc) DC_error ("need argument after -nlast ");
01167 sscanf (argv[nopt], "%d", &ival);
01168 if (ival < 0)
01169 DC_error ("illegal argument after -nlast ");
01170 option_data->NLast = ival;
01171 nopt++;
01172 continue;
01173 }
01174
01175
01176
01177 if (strcmp(argv[nopt], "-polort") == 0)
01178 {
01179 nopt++;
01180 if (nopt >= argc) DC_error ("need argument after -polort ");
01181 sscanf (argv[nopt], "%d", &ival);
01182 if (ival < -1)
01183 DC_error ("illegal argument after -polort ");
01184 option_data->polort = ival;
01185 nopt++;
01186 continue;
01187 }
01188
01189
01190
01191 if( strstr(argv[nopt],"legendre") != NULL ){
01192 legendre_polort = (strncmp(argv[nopt],"-leg",4) == 0) ;
01193 nopt++ ; continue ;
01194 }
01195
01196
01197
01198 if( strstr(argv[nopt],"svd") != NULL ){
01199 use_psinv = (strncmp(argv[nopt],"-svd",4) == 0) ;
01200 nopt++ ; continue ;
01201 }
01202
01203
01204
01205 if( strstr(argv[nopt],"dmbase") != NULL ){
01206 demean_base = (strncmp(argv[nopt],"-dmb",4) == 0) ;
01207 nopt++ ; continue ;
01208 }
01209
01210
01211
01212 if( strstr(argv[nopt],"xsave") != NULL ){
01213 xsave = (strncmp(argv[nopt],"-xsave",5) == 0) ;
01214 nopt++ ; continue ;
01215 }
01216
01217
01218
01219 if( strcmp(argv[nopt],"-xrestore") == 0 ){
01220 nopt++;
01221 if( nopt >= argc) DC_error ("need argument after -xrestore ");
01222 xrestore_filename = strdup( argv[nopt] ) ;
01223 if( !THD_is_file(xrestore_filename) )
01224 DC_error("file named after -xrestore doesn't exist") ;
01225 xrestore = 1 ; nopt++ ; continue ;
01226 }
01227
01228
01229 if (strcmp(argv[nopt], "-quiet") == 0)
01230 {
01231 option_data->quiet = 1; verb = 0 ;
01232 nopt++;
01233 continue;
01234 }
01235
01236
01237 if (strcmp(argv[nopt], "-progress") == 0)
01238 {
01239 nopt++;
01240 if (nopt >= argc) DC_error ("need argument after -progress ");
01241 sscanf (argv[nopt], "%d", &ival);
01242 if (ival < 0)
01243 DC_error ("illegal argument after -progress ");
01244 option_data->progress = ival;
01245 nopt++;
01246 continue;
01247 }
01248
01249
01250
01251 if (strcmp(argv[nopt], "-rmsmin") == 0)
01252 {
01253 nopt++;
01254 if (nopt >= argc) DC_error ("need argument after -rmsmin ");
01255 sscanf (argv[nopt], "%f", &fval);
01256 if (fval < 0.0)
01257 DC_error ("illegal argument after -rmsmin ");
01258 option_data->rms_min = fval;
01259 nopt++;
01260 continue;
01261 }
01262
01263
01264
01265 if (strcmp(argv[nopt], "-fdisp") == 0)
01266 {
01267 nopt++;
01268 if (nopt >= argc) DC_error ("need argument after -fdisp ");
01269 sscanf (argv[nopt], "%f", &fval);
01270 option_data->fdisp = fval;
01271 nopt++;
01272 continue;
01273 }
01274
01275
01276
01277 if (strcmp(argv[nopt], "-num_stimts") == 0)
01278 {
01279 nopt++;
01280 if (nopt >= argc) DC_error ("need argument after -num_stimts ");
01281 sscanf (argv[nopt], "%d", &ival);
01282 if (ival < 0)
01283 {
01284 DC_error ("-num_stimts num Require: num >= 0 ");
01285 }
01286
01287 initialize_stim_options (option_data, ival);
01288
01289 nopt++;
01290 continue;
01291 }
01292
01293
01294 if( strcmp(argv[nopt],"-TR_times") == 0 ){
01295 nopt++ ;
01296 if( nopt >= argc ) DC_error("need argument after -TR_times") ;
01297 sscanf( argv[nopt] , "%f" , &basis_dtout ) ;
01298 if( basis_dtout <= 0.0f ){
01299 fprintf(stderr,"** ERROR: -TR_times '%s' is illegal\n",argv[nopt]) ;
01300 exit(1) ;
01301 }
01302 nopt++ ; continue ;
01303 }
01304
01305
01306 if( strcmp(argv[nopt],"-TR_irc") == 0 ){
01307 nopt++ ;
01308 if( nopt >= argc ) DC_error("need argument after -TR_irc") ;
01309 sscanf( argv[nopt] , "%f" , &irc_dt ) ;
01310 if( irc_dt <= 0.0f ){
01311 fprintf(stderr,"** ERROR: -TR_irc '%s' is illegal\n",argv[nopt]) ;
01312 exit(1) ;
01313 }
01314 nopt++ ; continue ;
01315 }
01316
01317
01318
01319 if( strcmp(argv[nopt],"-basis_normall") == 0 ){
01320 nopt++ ;
01321 if( nopt >= argc ) DC_error("need argument after -basis_normall") ;
01322 basis_normall = strtod( argv[nopt] , NULL ) ;
01323 if( basis_normall <= 0.0f )
01324 DC_error("value after -basis_normall is illegal!") ;
01325 nopt++ ; continue ;
01326 }
01327
01328
01329 if( strcmp(argv[nopt],"-stim_times") == 0 ){
01330 nopt++ ;
01331 if( nopt+2 >= argc ) DC_error("need 3 arguments after -stim_times");
01332 sscanf( argv[nopt] , "%d" , &ival ) ;
01333 if( (ival < 1) || (ival > option_data->num_stimts) ){
01334 fprintf(stderr,
01335 "** ERROR: '-stim_times %d' value out of range 1..%d\n",
01336 ival , option_data->num_stimts ) ;
01337 exit(1) ;
01338 }
01339 k = ival-1 ; nopt++ ;
01340 if( option_data->stim_filename[k] != NULL ){
01341 fprintf(stderr,
01342 "** ERROR: '-stim_times %d' trying to overwrite previous stimulus\n",
01343 ival ) ;
01344 exit(1) ;
01345 }
01346 option_data->stim_filename[k] = strdup( argv[nopt] ) ;
01347 basis_times[k] = mri_read_ascii_ragged( argv[nopt] , basis_filler ) ;
01348 if( basis_times[k] == NULL ){
01349 fprintf(stderr,
01350 "** ERROR: '-stim_times %d' can't read file '%s'\n",
01351 ival , argv[nopt] ) ;
01352 exit(1) ;
01353 }
01354 nopt++ ;
01355 basis_stim[k] = basis_parser( argv[nopt] ) ;
01356 if( basis_stim[k] == NULL ){
01357 fprintf(stderr,
01358 "** ERROR: '-stim_times %d %s' can't parse '%s'\n",
01359 ival , argv[nopt-1] , argv[nopt] ) ;
01360 }
01361 basis_count++ ;
01362 nopt++ ; continue ;
01363 }
01364
01365
01366 if (strcmp(argv[nopt], "-stim_file") == 0)
01367 {
01368 nopt++;
01369 if (nopt+1 >= argc) DC_error ("need 2 arguments after -stim_file");
01370
01371 sscanf (argv[nopt], "%d", &ival);
01372 if ((ival < 1) || (ival > option_data->num_stimts))
01373 DC_error ("-stim_file k sname Require: 1 <= k <= num_stimts");
01374 k = ival-1;
01375 nopt++;
01376 if( option_data->stim_filename[k] != NULL ){
01377 fprintf(stderr,
01378 "** ERROR: '-stim_file %d' trying to overwrite previous stimulus\n",
01379 ival ) ;
01380 exit(1) ;
01381 }
01382
01383 option_data->stim_filename[k] = malloc (sizeof(char)*THD_MAX_NAME);
01384 MTEST (option_data->stim_filename[k]);
01385 strcpy (option_data->stim_filename[k], argv[nopt]);
01386 nopt++;
01387 continue;
01388 }
01389
01390
01391
01392 if (strcmp(argv[nopt], "-stim_label") == 0)
01393 {
01394 nopt++;
01395 if (nopt+1 >= argc) DC_error ("need 2 arguments after -stim_label");
01396
01397 sscanf (argv[nopt], "%d", &ival);
01398 if ((ival < 1) || (ival > option_data->num_stimts))
01399 DC_error ("-stim_label k slabel Require: 1 <= k <= num_stimts");
01400 k = ival-1;
01401 nopt++;
01402
01403 strcpy (option_data->stim_label[k], argv[nopt]);
01404 nopt++;
01405 continue;
01406 }
01407
01408
01409
01410 if (strcmp(argv[nopt], "-stim_base") == 0)
01411 {
01412 nopt++;
01413 if (nopt >= argc)
01414 DC_error ("need 1 argument after -stim_base");
01415
01416 sscanf (argv[nopt], "%d", &ival);
01417 if ((ival < 1) || (ival > option_data->num_stimts))
01418 DC_error ("-stim_base k Require: 1 <= k <= num_stimts");
01419 k = ival-1;
01420 option_data->stim_base[k] = 1;
01421 nopt++;
01422 continue;
01423 }
01424
01425
01426
01427 if (strcmp(argv[nopt], "-stim_minlag") == 0)
01428 {
01429 nopt++;
01430 if (nopt+1 >= argc)
01431 DC_error ("need 2 arguments after -stim_minlag");
01432
01433 sscanf (argv[nopt], "%d", &ival);
01434 if ((ival < 1) || (ival > option_data->num_stimts))
01435 DC_error ("-stim_minlag k lag Require: 1 <= k <= num_stimts");
01436 k = ival-1;
01437 nopt++;
01438
01439 sscanf (argv[nopt], "%d", &ival);
01440 if (ival < 0)
01441 DC_error ("-stim_minlag k lag Require: 0 <= lag");
01442 option_data->stim_minlag[k] = ival;
01443 nopt++;
01444 continue;
01445 }
01446
01447
01448
01449 if (strcmp(argv[nopt], "-stim_maxlag") == 0)
01450 {
01451 nopt++;
01452 if (nopt+1 >= argc)
01453 DC_error ("need 2 arguments after -stim_maxlag");
01454
01455 sscanf (argv[nopt], "%d", &ival);
01456 if ((ival < 1) || (ival > option_data->num_stimts))
01457 DC_error ("-stim_maxlag k lag Require: 1 <= k <= num_stimts");
01458 k = ival-1;
01459 nopt++;
01460
01461 sscanf (argv[nopt], "%d", &ival);
01462 if (ival < 0)
01463 DC_error ("-stim_maxlag k lag Require: 0 <= lag");
01464 option_data->stim_maxlag[k] = ival;
01465 nopt++;
01466 continue;
01467 }
01468
01469
01470
01471 if (strcmp(argv[nopt], "-stim_nptr") == 0)
01472 {
01473 nopt++;
01474 if (nopt+1 >= argc)
01475 DC_error ("need 2 arguments after -stim_nptr");
01476
01477 sscanf (argv[nopt], "%d", &ival);
01478 if ((ival < 1) || (ival > option_data->num_stimts))
01479 DC_error ("-stim_nptr k p Require: 1 <= k <= num_stimts");
01480 k = ival-1;
01481 nopt++;
01482
01483 sscanf (argv[nopt], "%d", &ival);
01484 if (ival < 1)
01485 DC_error ("-stim_nptr k p Require: 1 <= p");
01486 option_data->stim_nptr[k] = ival;
01487 nopt++;
01488 continue;
01489 }
01490
01491
01492
01493 if (strcmp(argv[nopt], "-num_glt") == 0)
01494 {
01495 nopt++;
01496 if (nopt >= argc) DC_error ("need argument after -num_glt ");
01497 sscanf (argv[nopt], "%d", &ival);
01498 if (ival < 0)
01499 {
01500 DC_error ("-num_glt num Require: num >= 0 ");
01501 }
01502
01503 if (option_data->num_glt > 0)
01504 DC_error ("-num_glt option must precede any -glt options ");
01505
01506 initialize_glt_options (option_data, ival);
01507
01508 nopt++;
01509 continue;
01510 }
01511
01512
01513
01514 if (strcmp(argv[nopt], "-glt") == 0)
01515 {
01516 nopt++;
01517 if (nopt+1 >= argc) DC_error ("need 2 arguments after -glt");
01518
01519 sscanf (argv[nopt], "%d", &ival);
01520 if (ival < 1)
01521 {
01522 DC_error ("-glt s gltname Require: s >= 1 (s = #rows in GLT)");
01523 }
01524 s = ival;
01525
01526 if (option_data->num_glt == 0)
01527 initialize_glt_options (option_data, 10);
01528
01529 if (iglt+1 > option_data->num_glt)
01530 DC_error ("Use -num_glt option to specify number of GLTs");
01531
01532 option_data->glt_rows[iglt] = s;
01533 nopt++;
01534
01535 option_data->glt_filename[iglt]
01536 = malloc (sizeof(char)*THD_MAX_NAME);
01537 MTEST (option_data->glt_filename[iglt]);
01538 strcpy (option_data->glt_filename[iglt],
01539 argv[nopt]);
01540 iglt++;
01541
01542 nopt++;
01543 continue;
01544 }
01545
01546
01547
01548 if( strcmp(argv[nopt],"-gltsym") == 0 ){
01549 nopt++ ;
01550 if( nopt >= argc ) DC_error("need 1 argument after -gltsym") ;
01551
01552 if( option_data->num_glt == 0 )
01553 initialize_glt_options(option_data,10) ;
01554
01555 if( iglt+1 > option_data->num_glt )
01556 DC_error("Use -num_glt option to specify number of GLTs when more than 10") ;
01557
01558 option_data->glt_rows[iglt] = -1 ;
01559
01560 option_data->glt_filename[iglt] = strdup( argv[nopt] ) ;
01561 iglt++ ; nopt++ ; continue ;
01562 }
01563
01564
01565
01566 if (strcmp(argv[nopt], "-glt_label") == 0)
01567 {
01568 nopt++;
01569 if (nopt+1 >= argc) DC_error ("need 2 arguments after -glt_label");
01570
01571 sscanf (argv[nopt], "%d", &ival);
01572 if ((ival < 1) || (ival > option_data->num_glt))
01573 DC_error ("-stim_label k slabel Require: 1 <= k <= num_glt");
01574 k = ival-1;
01575 nopt++;
01576
01577 strcpy (option_data->glt_label[k], argv[nopt]);
01578 nopt++;
01579 continue;
01580 }
01581
01582
01583
01584 if (strcmp(argv[nopt], "-iresp") == 0)
01585 {
01586 nopt++;
01587 if (nopt+1 >= argc) DC_error ("need 2 arguments after -iresp");
01588
01589 sscanf (argv[nopt], "%d", &ival);
01590 if ((ival < 1) || (ival > option_data->num_stimts))
01591 DC_error ("-iresp k iprefix Require: 1 <= k <= num_stimts");
01592 k = ival-1;
01593 nopt++;
01594
01595 option_data->iresp_filename[k]
01596 = malloc (sizeof(char)*THD_MAX_NAME);
01597 MTEST (option_data->iresp_filename[k]);
01598 strcpy (option_data->iresp_filename[k], argv[nopt]);
01599 nopt++;
01600 continue;
01601 }
01602
01603
01604
01605 if (strcmp(argv[nopt], "-tshift") == 0)
01606 {
01607 option_data->tshift = 1;
01608 nopt++;
01609 continue;
01610 }
01611
01612
01613
01614 if (strcmp(argv[nopt], "-sresp") == 0)
01615 {
01616 nopt++;
01617 if (nopt+1 >= argc) DC_error ("need 2 arguments after -sresp");
01618
01619 sscanf (argv[nopt], "%d", &ival);
01620 if ((ival < 1) || (ival > option_data->num_stimts))
01621 DC_error ("-sresp k iprefix Require: 1 <= k <= num_stimts");
01622 k = ival-1;
01623 nopt++;
01624
01625 option_data->sresp_filename[k]
01626 = malloc (sizeof(char)*THD_MAX_NAME);
01627 MTEST (option_data->sresp_filename[k]);
01628 strcpy (option_data->sresp_filename[k], argv[nopt]);
01629 nopt++;
01630 continue;
01631 }
01632
01633
01634
01635 if (strcmp(argv[nopt], "-fout") == 0)
01636 {
01637 option_data->fout = 1;
01638 nopt++;
01639 continue;
01640 }
01641
01642
01643
01644 if (strcmp(argv[nopt], "-rout") == 0)
01645 {
01646 option_data->rout = 1;
01647 nopt++;
01648 continue;
01649 }
01650
01651
01652
01653 if (strcmp(argv[nopt], "-tout") == 0)
01654 {
01655 option_data->tout = 1;
01656 nopt++;
01657 continue;
01658 }
01659
01660
01661
01662 if (strcmp(argv[nopt], "-vout") == 0)
01663 {
01664 option_data->vout = 1;
01665 nopt++;
01666 continue;
01667 }
01668
01669
01670
01671 if (strcmp(argv[nopt], "-xout") == 0)
01672 {
01673 option_data->xout = 1;
01674 nopt++;
01675 continue;
01676 }
01677
01678
01679
01680 if (strcmp(argv[nopt], "-nobout") == 0)
01681 {
01682 option_data->nobout = 1;
01683 nopt++;
01684 continue;
01685 }
01686
01687
01688
01689 if (strcmp(argv[nopt], "-nocout") == 0)
01690 {
01691 option_data->nocout = 1;
01692 nopt++;
01693 continue;
01694 }
01695
01696
01697
01698 if (strcmp(argv[nopt], "-full_first") == 0)
01699 {
01700 option_data->full_first = 1;
01701 nopt++;
01702 continue;
01703 }
01704
01705
01706
01707 if (strcmp(argv[nopt], "-bucket") == 0 || strcmp(argv[nopt],"-prefix") == 0 )
01708 {
01709 nopt++;
01710 if (nopt >= argc) DC_error ("need file prefixname after -bucket ");
01711 option_data->bucket_filename = malloc (sizeof(char)*THD_MAX_NAME);
01712 MTEST (option_data->bucket_filename);
01713 strcpy (option_data->bucket_filename, argv[nopt]);
01714 nopt++;
01715 continue;
01716 }
01717
01718
01719 if (strcmp(argv[nopt], "-cbucket") == 0 || strcmp(argv[nopt],"-cprefix") == 0 )
01720 {
01721 nopt++;
01722 if (nopt >= argc) DC_error ("need file prefixname after -cbucket ");
01723 CoefFilename = strdup( argv[nopt] ) ;
01724 nopt++; continue;
01725 }
01726
01727
01728
01729 if (strcmp(argv[nopt], "-fitts") == 0)
01730 {
01731 nopt++;
01732 if (nopt >= argc) DC_error ("need file prefixname after -fitts ");
01733 option_data->fitts_filename = malloc (sizeof(char)*THD_MAX_NAME);
01734 MTEST (option_data->fitts_filename);
01735 strcpy (option_data->fitts_filename, argv[nopt]);
01736 nopt++;
01737 continue;
01738 }
01739
01740
01741
01742 if (strcmp(argv[nopt], "-errts") == 0)
01743 {
01744 nopt++;
01745 if (nopt >= argc) DC_error ("need file prefixname after -errts ");
01746 option_data->errts_filename = malloc (sizeof(char)*THD_MAX_NAME);
01747 MTEST (option_data->errts_filename);
01748 strcpy (option_data->errts_filename, argv[nopt]);
01749 nopt++;
01750 continue;
01751 }
01752
01753
01754 if( strcmp(argv[nopt],"-jobs") == 0 ){
01755 nopt++ ;
01756 if (nopt >= argc) DC_error ("need J parameter after -jobs ");
01757 #ifdef PROC_MAX
01758 proc_numjob = strtol(argv[nopt],NULL,10) ;
01759 if( proc_numjob < 1 ){
01760 fprintf(stderr,"** NOTICE: setting number of processes to 1!\n") ;
01761 proc_numjob = 1 ;
01762 } else if( proc_numjob > PROC_MAX ){
01763 fprintf(stderr,"** NOTICE: setting number of processes to %d!\n",PROC_MAX);
01764 proc_numjob = PROC_MAX ;
01765 }
01766 #else
01767 fprintf(stderr,"** WARNING: -jobs not supported in this version\n") ;
01768 #endif
01769 proc_use_jobs = 1 ;
01770 nopt++; continue;
01771 }
01772
01773
01774
01775 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
01776 DC_error (message);
01777
01778 }
01779
01780
01781
01782 option_data->num_glt = iglt;
01783
01784
01785
01786 #ifdef PROC_MAX
01787 if( (xrestore || option_data->input1D_filename != NULL) && proc_numjob > 1 ){
01788 proc_numjob = 1 ;
01789 if( verb ) fprintf(stderr,"** WARNING: -jobs reset to 1\n") ;
01790 }
01791 #endif
01792
01793
01794
01795 if( option_data->polort < 0 ) demean_base = 0 ;
01796
01797 nerr = 0 ;
01798 for( k=0 ; k < option_data->num_stimts ; k++ ){
01799
01800 if( basis_stim[k] != NULL ){
01801
01802 basis_stim[k]->name = strdup( option_data->stim_label[k] ) ;
01803
01804 if( option_data->sresp_filename[k] != NULL ) basis_need_mse = 1 ;
01805
01806 if( option_data->stim_nptr[k] != 1 ){
01807 fprintf(stderr,
01808 "** ERROR: '-stim_nptr %d %d' illegal with '-stim_times'\n",
01809 k+1 , option_data->stim_nptr[k] ) ;
01810 nerr++ ;
01811 }
01812 if( option_data->stim_base[k] != 0 ){
01813 fprintf(stderr,
01814 "** ERROR: '-stim_base %d' illegal with '-stim_times'\n",
01815 k+1 ) ;
01816 nerr++ ;
01817 }
01818 if( option_data->stim_minlag[k] != 0 ){
01819 fprintf(stderr,
01820 "** ERROR: '-stim_minlag %d %d' illegal with '-stim_times'\n",
01821 k+1 , option_data->stim_minlag[k] ) ;
01822 nerr++ ;
01823 }
01824 if( option_data->stim_maxlag[k] != 0 ){
01825 fprintf(stderr,
01826 "** ERROR: '-stim_maxlag %d %d' illegal with '-stim_times'\n",
01827 k+1 , option_data->stim_maxlag[k] ) ;
01828 nerr++ ;
01829 }
01830
01831 } else {
01832
01833 if( option_data->stim_filename[k] == NULL ){
01834 fprintf(stderr,"** ERROR: no stimulus #%d given\n",k+1) ;
01835 nerr++ ;
01836 }
01837
01838 }
01839 }
01840 if( nerr > 0 ) exit(1) ;
01841
01842 }
01843
01844
01845
01846
01847
01848
01849
01850
01851 float * read_time_series
01852 (
01853 char * ts_filename,
01854 int * ts_length
01855 )
01856
01857 {
01858 char message[THD_MAX_NAME];
01859 char * cpt;
01860 char filename[THD_MAX_NAME];
01861 char subv[THD_MAX_NAME];
01862 MRI_IMAGE * im, * flim;
01863
01864 float * far;
01865 int nx;
01866 int ny;
01867 int iy;
01868 int ipt;
01869 float * ts_data = NULL;
01870
01871
01872 ENTRY("read_time_series") ;
01873
01874
01875 if (ts_filename == NULL)
01876 DC_error ("Missing input time series file name");
01877
01878
01879
01880 flim = mri_read_1D( ts_filename ) ;
01881 if (flim == NULL)
01882 {
01883 sprintf (message, "Unable to read time series file: %s", ts_filename);
01884 DC_error (message);
01885 }
01886 far = MRI_FLOAT_PTR(flim);
01887 nx = flim->nx;
01888 ny = flim->ny; iy = 0 ;
01889 if( ny > 1 ){
01890 fprintf(stderr,"** WARNING: time series %s has %d columns\n",ts_filename,ny);
01891 }
01892
01893
01894
01895 *ts_length = nx;
01896 ts_data = (float *) malloc (sizeof(float) * nx);
01897 MTEST (ts_data);
01898 for (ipt = 0; ipt < nx; ipt++)
01899 ts_data[ipt] = far[ipt + iy*nx];
01900
01901
01902 mri_free (flim); flim = NULL;
01903
01904 RETURN (ts_data);
01905 }
01906
01907
01908
01909
01910
01911
01912
01913 void read_input_data
01914 (
01915 DC_options * option_data,
01916 THD_3dim_dataset ** dset_time,
01917 byte ** mask_vol,
01918 float ** fmri_data,
01919 int * fmri_length,
01920 float ** censor_array,
01921 int * censor_length,
01922 int ** block_list,
01923 int * num_blocks,
01924 float *** stimulus,
01925 int ** stim_length,
01926 matrix ** glt_cmat
01927 )
01928
01929 {
01930 char message[THD_MAX_NAME];
01931 int it;
01932 int nt;
01933 int nxyz;
01934 int num_stimts;
01935 int * baseline;
01936 int * min_lag;
01937 int * max_lag;
01938 int qp;
01939 int q;
01940 int p;
01941 int is;
01942 int num_glt;
01943 int iglt;
01944
01945
01946 ENTRY("read_input_data") ;
01947
01948
01949 num_stimts = option_data->num_stimts;
01950 num_glt = option_data->num_glt;
01951 baseline = option_data->stim_base;
01952 min_lag = option_data->stim_minlag;
01953 max_lag = option_data->stim_maxlag;
01954
01955
01956
01957 if (num_stimts > 0)
01958 {
01959 int nerr=0 ;
01960
01961 *stimulus = (float **) calloc (sizeof(float *) , num_stimts);
01962 MTEST (*stimulus);
01963 *stim_length = (int *) calloc (sizeof(int) , num_stimts);
01964 MTEST (*stim_length);
01965
01966
01967
01968 { int js ;
01969 for( is=0 ; is < num_stimts ; is++ ){
01970 for( js=is+1 ; js < num_stimts ; js++ ){
01971 if( strcmp( option_data->stim_filename[is] ,
01972 option_data->stim_filename[js] ) == 0 )
01973 fprintf(stderr,"** WARNING: stimulus filename "
01974 "#%d:'%s' same as #%d:'%s'\n" ,
01975 is+1,option_data->stim_filename[is] ,
01976 js+1,option_data->stim_filename[js] ) ; nerr++ ;
01977 }
01978 }
01979 }
01980
01981 if( nerr > 0 && AFNI_yesenv("AFNI_3dDeconvolve_nodup") ){
01982 fprintf(stderr,"** ERROR: can't continue after above warnings!\n");
01983 exit(1) ;
01984 }
01985
01986 for (is = 0; is < num_stimts; is++)
01987 {
01988 if( basis_stim[is] != NULL ) continue ;
01989
01990 (*stimulus)[is] = read_time_series (option_data->stim_filename[is],
01991 &((*stim_length)[is]));
01992
01993 if ((*stimulus)[is] == NULL)
01994 {
01995 sprintf (message, "Unable to read stimulus time series: %s",
01996 option_data->stim_filename[is]);
01997 DC_error (message);
01998 }
01999
02000 }
02001 }
02002
02003
02004 *num_blocks = 0 ;
02005
02006
02007 if (option_data->nodata)
02008 {
02009
02010 if (num_stimts <= 0)
02011 DC_error ("Must have num_stimts > 0 for -nodata option");
02012
02013 if( basis_count > 0 ){
02014 if( option_data->nodata_TR > 0.0 ) basis_TR = option_data->nodata_TR ;
02015 else if( basis_dtout > 0.0 ) basis_TR = basis_dtout ;
02016 fprintf(stderr,
02017 "** NOTICE: using TR=%.3f for -stim_times and -nodata\n",basis_TR);
02018 }
02019
02020 *dset_time = NULL;
02021 if( option_data->nodata_NT > 0 )
02022 nt = option_data->nodata_NT ;
02023 else
02024 nt = (*stim_length)[0] / option_data->stim_nptr[0];
02025
02026 nxyz = 0;
02027 }
02028
02029 else if (option_data->input1D_filename != NULL)
02030 {
02031
02032 *fmri_data = read_time_series (option_data->input1D_filename,
02033 fmri_length);
02034 if (*fmri_data == NULL)
02035 {
02036 sprintf (message, "Unable to read time series file: %s",
02037 option_data->input1D_filename);
02038 DC_error (message);
02039 }
02040 *dset_time = NULL;
02041 nt = *fmri_length;
02042 nxyz = 1;
02043 }
02044
02045 else if (option_data->input_filename != NULL)
02046 {
02047
02048 *dset_time = THD_open_dataset (option_data->input_filename);
02049 if (!ISVALID_3DIM_DATASET(*dset_time))
02050 {
02051 sprintf (message, "Unable to open dataset '%s'",
02052 option_data->input_filename);
02053 DC_error (message);
02054 }
02055 THD_load_datablock ((*dset_time)->dblk);
02056 if( !DSET_LOADED(*dset_time) ){
02057 sprintf (message, "Unable to load dataset '%s'",
02058 option_data->input_filename);
02059 DC_error (message);
02060 }
02061 if( (*dset_time)->taxis == NULL ){
02062 fprintf(stderr,"** WARNING: dataset '%s' has no time axis!!\n",
02063 option_data->input_filename) ;
02064 }
02065 nt = DSET_NUM_TIMES (*dset_time);
02066 nxyz = DSET_NVOX (*dset_time);
02067
02068 basis_TR = DSET_TR(*dset_time) ;
02069 if( basis_TR <= 0.0f ) basis_TR = 1.0f ;
02070 if( DSET_TIMEUNITS(*dset_time) == UNITS_MSEC_TYPE ) basis_TR *= 0.001 ;
02071
02072 if( DSET_IS_TCAT(*dset_time) ){
02073 if( option_data->concat_filename != NULL ){
02074 fprintf(stderr,
02075 "** WARNING: '-concat %s' ignored: input dataset is auto-catenated\n" ,
02076 option_data->concat_filename ) ;
02077 option_data->concat_filename = NULL ;
02078 }
02079 *num_blocks = (*dset_time)->tcat_num ;
02080 *block_list = (int *) malloc (sizeof(int) * (*num_blocks));
02081 (*block_list)[0] = 0;
02082 for( it=0 ; it < (*num_blocks)-1 ; it++ )
02083 (*block_list)[it+1] = (*block_list)[it] + (*dset_time)->tcat_len[it] ;
02084 if( verb ){
02085 fprintf(stderr,"++ Auto-catenated datasets start at:") ;
02086 for( it=0 ; it < (*num_blocks) ; it++ )
02087 fprintf(stderr," %d",(*block_list)[it]) ;
02088 fprintf(stderr,"\n") ;
02089 }
02090 }
02091
02092 if( option_data->automask ){
02093 MRI_IMAGE *qim ; int mc ;
02094 qim = THD_rms_brick( *dset_time ) ;
02095 *mask_vol = mri_automask_image( qim ) ;
02096 mri_free( qim ) ;
02097 if( *mask_vol == NULL ){
02098 fprintf(stderr,"** WARNING: unable to generate automask?!\n") ;
02099 } else {
02100 mc = THD_countmask( DSET_NVOX(*dset_time) , *mask_vol ) ;
02101 if( mc <= 1 ){
02102 fprintf(stderr,"** WARNING: automask is empty!?\n") ;
02103 free(*mask_vol) ; *mask_vol = NULL ;
02104 } else {
02105 fprintf(stderr,"++ %d voxels in automask (out of %d)\n",
02106 mc,DSET_NVOX(*dset_time)) ;
02107 }
02108 }
02109 }
02110
02111 if (option_data->mask_filename != NULL)
02112 {
02113 THD_3dim_dataset * mask_dset = NULL;
02114
02115
02116 mask_dset = THD_open_dataset (option_data->mask_filename);
02117 if (!ISVALID_3DIM_DATASET(mask_dset))
02118 {
02119 sprintf (message, "Unable to open mask file: %s",
02120 option_data->mask_filename);
02121 DC_error (message);
02122 }
02123
02124
02125 if ( (DSET_NX(*dset_time) != DSET_NX(mask_dset))
02126 || (DSET_NY(*dset_time) != DSET_NY(mask_dset))
02127 || (DSET_NZ(*dset_time) != DSET_NZ(mask_dset)) )
02128 {
02129 sprintf (message, "datasets '%s' and '%s' have incompatible dimensions",
02130 option_data->input_filename, option_data->mask_filename);
02131 DC_error (message);
02132 }
02133
02134 if (DSET_NVALS(mask_dset) != 1 )
02135 DC_error ("Must specify exactly 1 sub-brick from mask dataset");
02136
02137 *mask_vol = THD_makemask( mask_dset , 0 , 1.0,0.0 ) ;
02138 if (*mask_vol == NULL) DC_error ("Unable to read mask dataset");
02139
02140 DSET_delete(mask_dset) ;
02141 }
02142
02143 } else {
02144 DC_error ("Must specify input data");
02145 }
02146
02147
02148
02149 if (nt <= 0) DC_error ("No time points?");
02150 option_data->nt = nt;
02151 if (nxyz < 0) DC_error ("Program initialization error: nxyz < 0");
02152 option_data->nxyz = nxyz;
02153
02154 voxel_num = nxyz ;
02155
02156
02157
02158 if (option_data->concat_filename == NULL)
02159 {
02160 if( *num_blocks <= 0 ){
02161 *num_blocks = 1;
02162 *block_list = (int *) malloc (sizeof(int) * 1);
02163 (*block_list)[0] = 0;
02164 }
02165 }
02166 else
02167 {
02168 float * f = NULL;
02169
02170 f = read_time_series (option_data->concat_filename, num_blocks);
02171 if (*num_blocks < 1)
02172 {
02173 sprintf (message, "Problem reading concat file: %s ",
02174 option_data->concat_filename);
02175 DC_error (message);
02176 }
02177 else
02178 {
02179 *block_list = (int *) malloc (sizeof(int) * (*num_blocks));
02180 for (it = 0; it < *num_blocks; it++)
02181 (*block_list)[it] = floor (f[it]+0.5);
02182 }
02183 }
02184
02185
02186
02187 if( basis_count > 0 ){
02188 MRI_IMAGE *tim , *qim ;
02189 float *tar , *qar , toff , tmax,tt,ss ;
02190 int ii,jj , kk , ngood , nx,ny , nf,dd ;
02191 int nbl=*num_blocks , *bst=*block_list ;
02192 basis_expansion *be ;
02193
02194 if( basis_TR <= 0.0f ) basis_TR = 1.0f ;
02195 if( basis_dtout <= 0.0f ) basis_dtout = basis_TR ;
02196
02197 for( is=0 ; is < num_stimts ; is++ ){
02198 be = basis_stim[is] ;
02199 if( be == NULL ) continue ;
02200
02201
02202
02203 tim = basis_times[is] ; nx = tim->nx ; ny = tim->ny ;
02204 ngood = 0 ; tar = MRI_FLOAT_PTR(tim) ;
02205 qar = (float *)calloc(sizeof(float),nx*ny) ;
02206
02207 if( nx == 1 ){
02208 tmax = (nt-1)*basis_TR ;
02209 for( ii=0 ; ii < nx*ny ; ii++ ){
02210 tt = tar[ii] ;
02211 if( tt >= 0.0f && tt <= tmax ) qar[ngood++] = tt/basis_TR ;
02212 }
02213 } else {
02214 if( ny != nbl ){
02215 fprintf(stderr,
02216 "** WARNING: '-stim_times %d' file '%s' has %d rows,"
02217 " but data has %d time blocks\n" ,
02218 is+1 , option_data->stim_filename[is] , ny,nbl ) ;
02219 if( ny > nbl ) ny = nbl ;
02220 }
02221 for( jj=0 ; jj < ny ; jj++ ){
02222 if( jj < nbl-1 ) tmax = (bst[jj+1]-1-bst[jj])*basis_TR ;
02223 else tmax = (nt -1-bst[jj])*basis_TR ;
02224 for( ii=0 ; ii < nx ; ii++ ){
02225 tt = tar[ii+jj*nx] ;
02226 if( tt >= 0.0f && tt <= tmax ) qar[ngood++] = tt/basis_TR+bst[jj] ;
02227 }
02228 }
02229 }
02230
02231 if( ngood == 0 ){
02232 fprintf(stderr,
02233 "** WARNING: '-stim_times %d' file '%s' has no good entries\n",
02234 is+1 , option_data->stim_filename[is] ) ;
02235 free((void *)qar) ; qim = NULL ;
02236 } else {
02237 qim = mri_new_vol_empty(ngood,1,1,MRI_float) ;
02238 qar = (float *)realloc((void *)qar,sizeof(float)*ngood) ;
02239 mri_fix_data_pointer( qar , qim ) ;
02240 }
02241
02242 mri_free(tim) ; basis_times[is] = qim ;
02243
02244
02245
02246 nf = basis_stim[is]->nfunc ;
02247 basis_vect[is] = mri_new( nt , nf , MRI_float ) ;
02248 if( qim != NULL ){
02249 int imin,imax , ibot,itop ;
02250 float *bv = MRI_FLOAT_PTR(basis_vect[is]) , dbot,dtop ;
02251
02252 #if 0
02253 fprintf(stderr,"-stim_times %d: adjusted time indexes follow:\n",is+1) ;
02254 for( kk=0 ; kk < ngood ; kk++ ) fprintf(stderr," %g",qar[kk]) ;
02255 fprintf(stderr,"\n") ;
02256 #endif
02257
02258 dbot = be->tbot / basis_TR ;
02259 dtop = be->ttop / basis_TR ;
02260 imin = 0 ; imax = nt-1 ;
02261
02262 for( kk=0 ; kk < ngood ; kk++ ){
02263 tt = qar[kk] ; if( tt < 0.0f || tt >= (float)nt ) continue ;
02264
02265 if( nbl > 1 ){
02266 for( jj=1 ; jj < nbl ; jj++ ) if( tt < bst[jj] ) break ;
02267 jj-- ;
02268 imin = bst[jj] ;
02269 if( jj < nbl-1 ) imax = bst[jj+1] - 1 ;
02270 else imax = nt - 1 ;
02271 }
02272
02273
02274
02275 ibot = (int)ceil ( tt + dbot ) ; if( ibot < imin ) ibot = imin ;
02276 itop = (int)floor( tt + dtop ) ; if( itop > imax ) itop = imax ;
02277
02278 for( ii=ibot ; ii <= itop ; ii++ ){
02279 ss = basis_TR*(ii-tt) ;
02280 if( ss < be->tbot || ss > be->ttop ) continue ;
02281 for( jj=0 ; jj < nf ; jj++ )
02282 bv[ii+jj*nt] += basis_funceval( be->bfunc[jj] , basis_TR*(ii-tt) );
02283 }
02284 }
02285 #if 0
02286 fprintf(stderr,"-stim_times %d: basis vectors follow:\n",is+1) ;
02287 for( ii=0 ; ii < nt ; ii++ ){
02288 fprintf(stderr,"%d:",ii) ;
02289 for( jj=0 ; jj < nf ; jj++ ) fprintf(stderr," %g",bv[ii+jj*nt]) ;
02290 fprintf(stderr,"\n") ;
02291 }
02292 #endif
02293 }
02294
02295 }
02296 }
02297
02298
02299
02300 qp = (option_data->polort + 1) * (*num_blocks);
02301 q = qp;
02302 p = qp;
02303 for (is = 0; is < num_stimts; is++){
02304 if( basis_stim[is] != NULL ){
02305 basis_stim[is]->pbot = p ;
02306 p += basis_stim[is]->nfunc ;
02307 } else {
02308 if (max_lag[is] < min_lag[is])
02309 DC_error ("Require min lag <= max lag for all stimuli");
02310 p += max_lag[is] - min_lag[is] + 1;
02311 if (baseline[is]) q += max_lag[is] - min_lag[is] + 1;
02312 }
02313 }
02314 option_data->p = p;
02315 option_data->q = q;
02316 option_data->qp = qp;
02317
02318
02319 if (option_data->censor_filename != NULL)
02320 {
02321
02322 *censor_array = read_time_series (option_data->censor_filename,
02323 censor_length);
02324 if (*censor_array == NULL)
02325 {
02326 sprintf (message, "Unable to read censor time series file: %s",
02327 option_data->censor_filename);
02328 DC_error (message);
02329 }
02330 }
02331 else
02332 {
02333
02334 *censor_array = (float *) malloc (sizeof(float) * nt);
02335 MTEST (*censor_array);
02336 *censor_length = nt;
02337 for (it = 0; it < nt; it++)
02338 (*censor_array)[it] = 1.0;
02339 }
02340
02341
02342
02343
02344 nSymStim = num_stimts+1 ;
02345 SymStim = (SYM_irange *)calloc(sizeof(SYM_irange),nSymStim) ;
02346
02347 strcpy( SymStim[num_stimts].name , "Ort" ) ;
02348 SymStim[num_stimts].nbot = 0 ;
02349 SymStim[num_stimts].ntop = qp-1 ;
02350 SymStim[num_stimts].gbot = 0 ;
02351 it = qp ;
02352 for( is=0 ; is < num_stimts ; is++ ){
02353 MCW_strncpy( SymStim[is].name , option_data->stim_label[is] , 64 ) ;
02354 if( basis_stim[is] != NULL ){
02355 SymStim[is].nbot = 0 ;
02356 SymStim[is].ntop = basis_stim[is]->nfunc-1 ;
02357 } else {
02358 SymStim[is].nbot = min_lag[is] ;
02359 SymStim[is].ntop = max_lag[is] ;
02360 }
02361 SymStim[is].gbot = it ;
02362 it += SymStim[is].ntop - SymStim[is].nbot + 1 ;
02363 if( strchr(SymStim[is].name,' ') != NULL ||
02364 strchr(SymStim[is].name,'*') != NULL ||
02365 strchr(SymStim[is].name,';') != NULL ){
02366 fprintf(stderr,
02367 "** WARNING: -stim_label #%d '%s' has characters bad for -gltsym\n",
02368 is+1 , SymStim[is].name ) ;
02369 }
02370 }
02371
02372
02373 if (num_glt > 0)
02374 {
02375 *glt_cmat = (matrix *) malloc (sizeof(matrix) * num_glt);
02376
02377
02378 for (iglt = 0; iglt < num_glt; iglt++)
02379 matrix_initialize (&((*glt_cmat)[iglt]));
02380
02381
02382 for (iglt = 0; iglt < num_glt; iglt++)
02383 {
02384 #if 1
02385 read_glt_matrix( option_data->glt_filename[iglt] ,
02386 option_data->glt_rows + iglt ,
02387 p , *glt_cmat + iglt ) ;
02388 #else
02389 matrix_file_read (option_data->glt_filename[iglt],
02390 option_data->glt_rows[iglt],
02391 p, &((*glt_cmat)[iglt]), 1);
02392 if ((*glt_cmat)[iglt].elts == NULL)
02393 {
02394 sprintf (message, "Unable to read GLT matrix from file: %s",
02395 option_data->glt_filename[iglt]);
02396 DC_error (message);
02397 }
02398 #endif
02399 }
02400 }
02401
02402 EXRETURN ;
02403 }
02404
02405
02406
02407
02408
02409
02410
02411
02412 void remove_zero_stimfns
02413 (
02414 DC_options * option_data,
02415 float ** stimulus,
02416 int * stim_length,
02417 matrix * glt_cmat
02418 )
02419
02420 {
02421 int num_stimts;
02422 int is, isp;
02423 int it;
02424 int num_glt;
02425 int iglt;
02426 int all_zero;
02427
02428
02429 ENTRY("remove_zero_stimfns") ;
02430
02431
02432 num_stimts = option_data->num_stimts;
02433 num_glt = option_data->num_glt;
02434
02435
02436
02437 is = 0;
02438 while (is < num_stimts)
02439 {
02440 if( basis_stim[is] != NULL ){ is++ ; continue ; }
02441
02442
02443 all_zero = TRUE;
02444 for (it = 0; it < stim_length[is]; it++)
02445 {
02446 if (stimulus[is][it] != 0.0)
02447 {
02448 all_zero = FALSE;
02449 break;
02450 }
02451 }
02452
02453 if (all_zero)
02454 {
02455 printf("** WARNING! -stim_file function %s comprises all zeros!\n",
02456 option_data->stim_filename[is]); fflush(stdout);
02457 if (option_data->num_glt > 0)
02458 DC_error
02459 ("Cannot process -glt option(s) when -stim_file function is all zero");
02460 if( basis_count > 0 )
02461 DC_error
02462 ("Cannot process -basis_stim option(s) when -stim_file function is all zero");
02463
02464 option_data->p -=
02465 option_data->stim_maxlag[is] - option_data->stim_minlag[is] + 1;
02466 if (option_data->stim_base[is])
02467 option_data->q -=
02468 option_data->stim_maxlag[is] - option_data->stim_minlag[is] + 1;
02469 for (isp = is; isp < num_stimts-1; isp++)
02470 {
02471 stimulus[isp] = stimulus[isp+1];
02472 stim_length[isp] = stim_length[isp+1];
02473 option_data->stim_filename[isp]
02474 = option_data->stim_filename[isp+1];
02475 option_data->stim_label[isp] = option_data->stim_label[isp+1];
02476 option_data->stim_base[isp]
02477 = option_data->stim_base[isp+1];
02478 option_data->stim_minlag[isp] = option_data->stim_minlag[isp+1];
02479 option_data->stim_maxlag[isp] = option_data->stim_maxlag[isp+1];
02480 option_data->iresp_filename[isp]
02481 = option_data->iresp_filename[isp+1];
02482 option_data->sresp_filename[isp]
02483 = option_data->sresp_filename[isp+1];
02484 }
02485
02486 num_stimts--;
02487 option_data->num_stimts = num_stimts;
02488 }
02489 else
02490 is++;
02491 }
02492
02493 EXRETURN ;
02494 }
02495
02496
02497
02498
02499
02500
02501
02502 void check_one_output_file
02503 (
02504 THD_3dim_dataset * dset_time,
02505 char * filename
02506 )
02507
02508 {
02509 char message[THD_MAX_NAME];
02510 THD_3dim_dataset * new_dset=NULL;
02511 int ierror;
02512
02513
02514
02515 new_dset = EDIT_empty_copy( dset_time ) ;
02516
02517
02518 ierror = EDIT_dset_items( new_dset ,
02519 ADN_prefix , filename ,
02520 ADN_label1 , filename ,
02521 ADN_self_name , filename ,
02522 ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE :
02523 GEN_FUNC_TYPE ,
02524 ADN_none ) ;
02525
02526 if( ierror > 0 )
02527 {
02528 sprintf (message,
02529 "*** %d errors in attempting to create output dataset!\n",
02530 ierror);
02531 DC_error (message);
02532 }
02533
02534 if( THD_is_file(new_dset->dblk->diskptr->header_name) )
02535 {
02536 sprintf (message,
02537 "Output dataset file %s already exists "
02538 " -- cannot continue! ",
02539 new_dset->dblk->diskptr->header_name);
02540 DC_error (message);
02541 }
02542
02543
02544 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
02545
02546 }
02547
02548
02549
02550
02551
02552
02553
02554 void check_output_files
02555 (
02556 DC_options * option_data,
02557 THD_3dim_dataset * dset_time
02558 )
02559
02560 {
02561 int is;
02562
02563
02564 if (option_data->bucket_filename != NULL)
02565 check_one_output_file (dset_time, option_data->bucket_filename);
02566
02567 if (option_data->fitts_filename != NULL)
02568 check_one_output_file (dset_time, option_data->fitts_filename);
02569
02570 if (option_data->errts_filename != NULL)
02571 check_one_output_file (dset_time, option_data->errts_filename);
02572
02573 for (is = 0; is < option_data->num_stimts; is++)
02574 {
02575 if (option_data->iresp_filename[is] != NULL)
02576 check_one_output_file (dset_time, option_data->iresp_filename[is]);
02577
02578 if (option_data->sresp_filename[is] != NULL)
02579 check_one_output_file (dset_time, option_data->sresp_filename[is]);
02580 }
02581 }
02582
02583
02584
02585
02586
02587
02588
02589 void check_for_valid_inputs
02590 (
02591 DC_options * option_data,
02592 THD_3dim_dataset * dset_time,
02593 int fmri_length,
02594 float * censor_array,
02595 int censor_length,
02596 int * block_list,
02597 int num_blocks,
02598 int * stim_length,
02599 float **stimulus ,
02600 int ** good_list
02601 )
02602
02603 {
02604 char message[THD_MAX_NAME];
02605 int is;
02606 int num_stimts;
02607 int * min_lag;
02608 int * max_lag;
02609 int * nptr;
02610 int m;
02611 int qp;
02612 int q;
02613 int p;
02614 int nbricks;
02615 int it;
02616 int nt;
02617 int NFirst;
02618 int NLast;
02619 int N;
02620 int ib;
02621 int irb;
02622 int num_glt;
02623 int * glt_rows;
02624 int iglt;
02625 int nerr=0 ;
02626
02627
02628
02629 nt = option_data->nt;
02630 num_stimts = option_data->num_stimts;
02631 min_lag = option_data->stim_minlag;
02632 max_lag = option_data->stim_maxlag;
02633 num_glt = option_data->num_glt;
02634 glt_rows = option_data->glt_rows;
02635 nptr = option_data->stim_nptr;
02636 p = option_data->p;
02637 q = option_data->q;
02638 qp = option_data->qp;
02639
02640
02641 if( xsave && option_data->bucket_filename == NULL ){
02642 fprintf(stderr,"** WARNING: -xsave given without -bucket; -xsave is disabled!\n") ;
02643 xsave = 0 ;
02644 }
02645
02646
02647 if (censor_length < nt)
02648 {
02649 sprintf (message, "Input censor time series file %s is too short",
02650 option_data->censor_filename);
02651 DC_error (message);
02652 }
02653
02654
02655
02656 for (ib = 0; ib < num_blocks; ib++)
02657 if ((block_list[ib] < 0) || (block_list[ib] >= nt))
02658 {
02659 sprintf (message, "Invalid -concat input: %d ", block_list[ib]);
02660 DC_error (message);
02661 }
02662 if (num_blocks > 1)
02663 for (ib = 1; ib < num_blocks; ib++)
02664 if (block_list[ib] <= block_list[ib-1])
02665 DC_error ("Invalid concatenated runs list");
02666
02667
02668
02669 *good_list = (int *) malloc (sizeof(int) * nt); MTEST (*good_list);
02670 NFirst = option_data->NFirst;
02671 if (NFirst < 0){
02672 for (is = 0; is < num_stimts; is++)
02673 if( basis_stim[is] == NULL && NFirst < (max_lag[is]+nptr[is]-1)/nptr[is] )
02674 NFirst = (max_lag[is]+nptr[is]-1)/nptr[is];
02675 }
02676 NLast = option_data->NLast;
02677 if (NLast < 0) NLast = nt;
02678
02679 N = 0;
02680 ib = 0;
02681 for (it = block_list[0]; it < nt; it++)
02682 {
02683 if (ib+1 < num_blocks)
02684 if (it >= block_list[ib+1]) ib++;
02685
02686 irb = it - block_list[ib];
02687
02688 if ((irb >= NFirst) && (irb <= NLast) && (censor_array[it]))
02689 {
02690 (*good_list)[N] = it;
02691 N++;
02692 }
02693 }
02694
02695
02696
02697 if (N == 0) DC_error ("No usable time points?");
02698 if (N <= p)
02699 {
02700 sprintf (message, "Insufficient data for estimating %d parameters", p);
02701 DC_error (message);
02702 }
02703 option_data->N = N;
02704
02705 GoodList = *good_list ;
02706 nGoodList = N ;
02707 nParam = p ;
02708
02709
02710 if (num_stimts < 0)
02711 {
02712 DC_error ("Require: 0 <= num_stimts (DUH! --or-- D'oh!)");
02713 }
02714
02715
02716
02717
02718 #define ALLOW_EXTEND
02719 for (is = 0; is < num_stimts; is++)
02720 {
02721 if( basis_stim[is] != NULL ) continue ;
02722
02723 if (stim_length[is] < nt*nptr[is])
02724 {
02725 #ifndef ALLOW_EXTEND
02726 sprintf (message, "Input stimulus time series file %s is too short",
02727 option_data->stim_filename[is]);
02728 DC_error (message);
02729 #else
02730 int nlen=nt*nptr[is], qq ;
02731 fprintf(stderr,
02732 "** WARNING: input stimulus time series file %s is too short:\n"
02733 " length = %d, but should be at least %d.\n" ,
02734 option_data->stim_filename[is] , stim_length[is] , nlen ) ;
02735 stimulus[is] = (float *) realloc( stimulus[is] , sizeof(float)*nlen ) ;
02736 for( qq=stim_length[is] ; qq < nlen ; qq++ ) stimulus[is][qq] = 0.0 ;
02737 stim_length[is] = nlen ; nerr++ ;
02738 #endif
02739 }
02740 }
02741 #ifdef ALLOW_EXTEND
02742 if( nerr > 0 ){
02743 char *eee = getenv("AFNI_3dDeconvolve_extend") ;
02744 if( eee != NULL && (*eee=='n' || *eee=='N') ){
02745 fprintf(stderr,"** ERROR: Can't continue with too short files!\n"
02746 " AFNI_3dDeconvolve_extend = %s\n",eee ) ;
02747 exit(1) ;
02748 }
02749 fprintf(stderr,"++ EXTENDING short files with zero values\n"
02750 " (to stop this behavior, setenv AFNI_3dDeconvolve_extend NO)\n") ;
02751 }
02752 #endif
02753
02754
02755
02756 for (is = 0; is < num_stimts; is++)
02757 {
02758 if( basis_stim[is] != NULL ) continue ;
02759
02760 m = max_lag[is] - min_lag[is] + 1;
02761 if (m < 2)
02762 {
02763 if (option_data->iresp_filename[is] != NULL)
02764 {
02765 sprintf (message, "Only %d time point for output dataset %s ",
02766 m, option_data->iresp_filename[is]);
02767 DC_error (message);
02768 }
02769
02770 if (option_data->sresp_filename[is] != NULL)
02771 {
02772 sprintf (message, "Only %d time point for output dataset %s",
02773 m, option_data->sresp_filename[is]);
02774 DC_error (message);
02775 }
02776 }
02777 if ((m < 4) && (option_data->tshift))
02778 {
02779 if (option_data->iresp_filename[is] != NULL)
02780 {
02781 sprintf (message, "Only %d time points for 3d+time dataset %s\n",
02782 m, option_data->iresp_filename[is]);
02783 strcat (message, "Require >= 4 data points for -tshift option");
02784 DC_error (message);
02785 }
02786 }
02787 }
02788
02789
02790
02791
02792 nbricks = 0;
02793 if (option_data->bucket_filename != NULL)
02794 {
02795 if (! option_data->nocout)
02796 {
02797 if (! option_data->nobout)
02798 nbricks += qp * (1 + option_data->tout);
02799
02800 for (is = 0; is < num_stimts; is++)
02801 {
02802 if ((!option_data->stim_base[is]) || (!option_data->nobout))
02803 {
02804 if( basis_stim[is] != NULL ) m = basis_stim[is]->nfunc ;
02805 else m = max_lag[is] - min_lag[is] + 1;
02806 nbricks += m * (1 + option_data->tout);
02807 nbricks += option_data->rout + option_data->fout;
02808 }
02809 }
02810 }
02811
02812 nbricks += option_data->rout + option_data->fout + option_data->vout;
02813
02814 if (num_glt > 0)
02815 for (iglt = 0; iglt < num_glt; iglt++)
02816 {
02817 nbricks += glt_rows[iglt] * (1 + option_data->tout);
02818 nbricks += option_data->rout + option_data->fout;
02819 }
02820
02821 if (nbricks <= 0)
02822 {
02823 sprintf (message,
02824 "User requested bucket dataset with only %d sub-bricks",
02825 nbricks);
02826 DC_error (message);
02827 }
02828
02829 }
02830 option_data->nbricks = nbricks;
02831
02832
02833
02834 if (dset_time != NULL)
02835 if (dset_time->taxis->nsl > 0)
02836 for (is = 0; is < num_stimts; is++)
02837 if (nptr[is] > 1)
02838 {
02839 sprintf (message, "Must align all slices to 0 offset time, \n ");
02840 strcat (message, "before using -stim_nptr option. ");
02841 strcat (message, "See program 3dTshift. ");
02842 DC_error (message);
02843 }
02844
02845
02846
02847 if ((option_data->tshift) && (option_data->input1D_filename != NULL))
02848 DC_error ("-tshift option is not compatible with -input1D option");
02849
02850
02851
02852 if (option_data->input_filename != NULL)
02853 check_output_files (option_data, dset_time);
02854
02855 }
02856
02857
02858
02859
02860
02861
02862
02863 void zero_fill_volume (float ** fvol, int nxyz)
02864 {
02865 int ixyz;
02866
02867 if( proc_numjob == 1 ){
02868
02869 *fvol = (float *) malloc (sizeof(float) * nxyz); MTEST(*fvol);
02870 for (ixyz = 0; ixyz < nxyz; ixyz++)
02871 (*fvol)[ixyz] = 0.0;
02872
02873 }
02874 #ifdef PROC_MAX
02875 else {
02876
02877
02878 proc_shm_arnum++ ;
02879 proc_shm_arsiz = (int *) realloc( proc_shm_arsiz ,
02880 sizeof(int) *proc_shm_arnum ) ;
02881 proc_shm_ar = (float ***) realloc( proc_shm_ar ,
02882 sizeof(float **)*proc_shm_arnum ) ;
02883 proc_shm_arsiz[proc_shm_arnum-1] = nxyz ;
02884 proc_shm_ar[proc_shm_arnum-1] = fvol ;
02885
02886
02887
02888 }
02889 #endif
02890 }
02891
02892
02893
02894 #ifdef PROC_MAX
02895
02896
02897
02898 #include <signal.h>
02899
02900 void proc_sigfunc(int sig)
02901 {
02902 char * sname ; int ii ;
02903 static volatile int fff=0 ;
02904 if( fff ) _exit(1); else fff=1 ;
02905 switch(sig){
02906 default: sname = "unknown" ; break ;
02907 case SIGHUP: sname = "SIGHUP" ; break ;
02908 case SIGTERM: sname = "SIGTERM" ; break ;
02909 case SIGILL: sname = "SIGILL" ; break ;
02910 case SIGKILL: sname = "SIGKILL" ; break ;
02911 case SIGPIPE: sname = "SIGPIPE" ; break ;
02912 case SIGSEGV: sname = "SIGSEGV" ; break ;
02913 case SIGBUS: sname = "SIGBUS" ; break ;
02914 case SIGINT: sname = "SIGINT" ; break ;
02915 case SIGFPE: sname = "SIGFPE" ; break ;
02916 }
02917 if( proc_shmid > 0 ){
02918 shmctl( proc_shmid , IPC_RMID , NULL ) ; proc_shmid = 0 ;
02919 }
02920 fprintf(stderr,"Fatal Signal %d (%s) received in job #%d\n",
02921 sig,sname,proc_ind) ;
02922 exit(1) ;
02923 }
02924
02925 void proc_atexit( void )
02926 {
02927 if( proc_shmid > 0 )
02928 shmctl( proc_shmid , IPC_RMID , NULL ) ;
02929 }
02930
02931
02932
02933
02934
02935 void proc_finalize_shm_volumes(void)
02936 {
02937 char kstr[32] ; int ii ;
02938
02939 if( proc_shm_arnum == 0 ) return ;
02940
02941 proc_shmsize = 0 ;
02942 for( ii=0 ; ii < proc_shm_arnum ; ii++ )
02943 proc_shmsize += proc_shm_arsiz[ii] ;
02944
02945 proc_shmsize *= sizeof(float) ;
02946
02947
02948
02949 UNIQ_idcode_fill( kstr ) ;
02950 proc_shmid = shm_create( kstr , proc_shmsize ) ;
02951 if( proc_shmid < 0 ){
02952 fprintf(stderr,"\n** Can't create shared memory of size %d!\n"
02953 "** Try re-running without -jobs option!\n" ,
02954 proc_shmsize ) ;
02955
02956
02957
02958 { char *cmd=NULL ;
02959 #if defined(LINUX)
02960 cmd = "cat /proc/sys/kernel/shmmax" ;
02961 #elif defined(SOLARIS)
02962 cmd = "/usr/sbin/sysdef | grep SHMMAX" ;
02963 #elif defined(DARWIN)
02964 cmd = "sysctl -n kern.sysv.shmmax" ;
02965 #endif
02966 if( cmd != NULL ){
02967 FILE *fp = popen( cmd , "r" ) ;
02968 if( fp != NULL ){
02969 unsigned int smax=0 ;
02970 fscanf(fp,"%u",&smax) ; pclose(fp) ;
02971 if( smax > 0 )
02972 fprintf(stderr ,
02973 "\n"
02974 "** POSSIBLY USEFUL ADVICE:\n"
02975 "** Current max shared memory size = %u bytes.\n"
02976 "** For information on how to change this, see\n"
02977 "** http://afni.nimh.nih.gov/afni/doc/misc/parallize.html \n"
02978 "** and also contact your system administrator.\n"
02979 , smax ) ;
02980 }
02981 }
02982 }
02983 exit(1) ;
02984 }
02985
02986
02987
02988
02989 signal(SIGPIPE,proc_sigfunc) ; signal(SIGSEGV,proc_sigfunc) ;
02990 signal(SIGINT ,proc_sigfunc) ; signal(SIGFPE ,proc_sigfunc) ;
02991 signal(SIGBUS ,proc_sigfunc) ; signal(SIGHUP ,proc_sigfunc) ;
02992 signal(SIGTERM,proc_sigfunc) ; signal(SIGILL ,proc_sigfunc) ;
02993 signal(SIGKILL,proc_sigfunc) ; signal(SIGPIPE,proc_sigfunc) ;
02994 atexit(proc_atexit) ;
02995
02996 fprintf(stderr , "++ Shared memory: %d bytes at id=%d\n" ,
02997 proc_shmsize , proc_shmid ) ;
02998
02999
03000
03001 proc_shmptr = shm_attach( proc_shmid ) ;
03002 if( proc_shmptr == NULL ){
03003 fprintf(stderr,"\n** Can't attach to shared memory!?\n"
03004 "** This is bizarre.\n" ) ;
03005 shmctl( proc_shmid , IPC_RMID , NULL ) ;
03006 exit(1) ;
03007 }
03008
03009
03010
03011 memset( proc_shmptr , 0 , proc_shmsize ) ;
03012
03013
03014
03015 *proc_shm_ar[0] = (float *) proc_shmptr ;
03016 for( ii=1 ; ii < proc_shm_arnum ; ii++ )
03017 *proc_shm_ar[ii] = *proc_shm_ar[ii-1] + proc_shm_arsiz[ii] ;
03018 }
03019
03020
03021
03022
03023
03024 void proc_free( void *ptr )
03025 {
03026 int ii ;
03027
03028 if( ptr == NULL ) return ;
03029 if( proc_shmid == 0 ){ free(ptr); return; }
03030 for( ii=0 ; ii < proc_shm_arnum ; ii++ )
03031 if( ((float *)ptr) == *proc_shm_ar[ii] ) return;
03032 free(ptr); return;
03033 }
03034
03035 #undef free
03036 #define free proc_free
03037
03038 #endif
03039
03040
03041
03042
03043
03044
03045 void allocate_memory
03046 (
03047 DC_options * option_data,
03048
03049 float *** coef_vol,
03050 float *** scoef_vol,
03051 float *** tcoef_vol,
03052 float *** fpart_vol,
03053 float *** rpart_vol,
03054 float ** mse_vol,
03055 float ** ffull_vol,
03056 float ** rfull_vol,
03057
03058 float **** glt_coef_vol,
03059 float **** glt_tcoef_vol,
03060 float *** glt_fstat_vol,
03061 float *** glt_rstat_vol,
03062
03063 float *** fitts_vol,
03064 float *** errts_vol
03065 )
03066
03067 {
03068 int ip;
03069 int qp;
03070 int q;
03071 int p;
03072 int nxyz;
03073 int ixyz;
03074 int is;
03075 int num_stimts;
03076 int num_glt;
03077 int iglt;
03078 int nlc;
03079 int ilc;
03080 int it;
03081 int nt;
03082 int * min_lag;
03083 int * max_lag;
03084
03085 int fout;
03086 int rout;
03087 int tout;
03088 int vout;
03089 int bout;
03090 int cout;
03091 int ibot,itop ;
03092
03093
03094 ENTRY("allocate_memory") ;
03095
03096
03097 nxyz = option_data->nxyz;
03098 num_stimts = option_data->num_stimts;
03099 num_glt = option_data->num_glt;
03100 qp = option_data->qp;
03101 q = option_data->q;
03102 p = option_data->p;
03103 nt = option_data->nt;
03104 min_lag = option_data->stim_minlag;
03105 max_lag = option_data->stim_maxlag;
03106
03107 fout = option_data->fout;
03108 rout = option_data->rout;
03109 tout = option_data->tout;
03110 vout = option_data->vout;
03111 bout = 1 - (option_data->nobout || option_data->nocout);
03112 cout = 1 - option_data->nocout;
03113
03114 if( CoefFilename != NULL ){ bout = cout = 1 ; }
03115
03116
03117 *coef_vol = (float **) malloc (sizeof(float *) * p); MTEST(*coef_vol);
03118 *scoef_vol = (float **) malloc (sizeof(float *) * p); MTEST(*scoef_vol);
03119 *tcoef_vol = (float **) malloc (sizeof(float *) * p); MTEST(*tcoef_vol);
03120 for (ip = 0; ip < p; ip++)
03121 {
03122 (*coef_vol)[ip] = NULL;
03123 (*scoef_vol)[ip] = NULL;
03124 (*tcoef_vol)[ip] = NULL;
03125 }
03126
03127 if (bout)
03128 for (ip = 0; ip < qp; ip++)
03129 {
03130 zero_fill_volume (&((*coef_vol)[ip]), nxyz);
03131 if (tout) zero_fill_volume (&((*tcoef_vol)[ip]), nxyz);
03132 }
03133
03134 ip = qp - 1;
03135 for (is = 0; is < num_stimts; is++)
03136 {
03137 if( basis_stim[is] != NULL ){ ibot=0 ; itop=basis_stim[is]->nfunc-1 ; }
03138 else { ibot=min_lag[is] ; itop=max_lag[is] ; }
03139 for (it = ibot; it <= itop; it++)
03140 {
03141 ip++;
03142 if (option_data->stim_base[is])
03143 {
03144 if (bout || (option_data->iresp_filename[is] != NULL))
03145 zero_fill_volume (&((*coef_vol)[ip]), nxyz);
03146 if (bout && tout)
03147 zero_fill_volume (&((*tcoef_vol)[ip]), nxyz);
03148 }
03149 else
03150 {
03151 if (cout || (option_data->iresp_filename[is] != NULL))
03152 zero_fill_volume (&((*coef_vol)[ip]), nxyz);
03153 if (cout && tout)
03154 zero_fill_volume (&((*tcoef_vol)[ip]), nxyz);
03155 }
03156 if (option_data->sresp_filename[is] != NULL)
03157 zero_fill_volume (&((*scoef_vol)[ip]), nxyz);
03158 }
03159 }
03160
03161
03162 if (fout)
03163 {
03164 *fpart_vol = (float **) malloc (sizeof(float *) * num_stimts);
03165 MTEST(*fpart_vol);
03166 for (is = 0; is < num_stimts; is++)
03167 if ((! option_data->stim_base[is]) || (! option_data->nobout))
03168 zero_fill_volume (&((*fpart_vol)[is]), nxyz);
03169 else
03170 (*fpart_vol)[is] = NULL;
03171 }
03172
03173
03174 if (rout)
03175 {
03176 *rpart_vol = (float **) malloc (sizeof(float *) * num_stimts);
03177 MTEST(*rpart_vol);
03178 for (is = 0; is < num_stimts; is++)
03179 if ((! option_data->stim_base[is]) || (! option_data->nobout))
03180 zero_fill_volume (&((*rpart_vol)[is]), nxyz);
03181 else
03182 (*rpart_vol)[is] = NULL;
03183 }
03184
03185
03186 if (vout || basis_need_mse) zero_fill_volume (&(*mse_vol), nxyz);
03187 if (fout) zero_fill_volume (&(*ffull_vol), nxyz);
03188 if (rout) zero_fill_volume (&(*rfull_vol), nxyz);
03189
03190
03191
03192 if (num_glt > 0)
03193 {
03194 *glt_coef_vol = (float ***) malloc (sizeof(float **) * num_glt);
03195 MTEST(*glt_coef_vol);
03196
03197 if (tout)
03198 {
03199 *glt_tcoef_vol = (float ***) malloc (sizeof(float **) * num_glt);
03200 MTEST(*glt_tcoef_vol);
03201 }
03202
03203 if (fout)
03204 {
03205 *glt_fstat_vol = (float **) malloc (sizeof(float *) * num_glt);
03206 MTEST(*glt_fstat_vol);
03207 }
03208
03209 if (rout)
03210 {
03211 *glt_rstat_vol = (float **) malloc (sizeof(float *) * num_glt);
03212 MTEST(*glt_rstat_vol);
03213 }
03214
03215 for (iglt = 0; iglt < num_glt; iglt++)
03216 {
03217 nlc = option_data->glt_rows[iglt];
03218
03219 (*glt_coef_vol)[iglt] = (float **) malloc (sizeof(float *) * nlc);
03220 MTEST((*glt_coef_vol)[iglt]);
03221
03222 if (tout)
03223 {
03224 (*glt_tcoef_vol)[iglt]
03225 = (float **) malloc (sizeof(float *) * nlc);
03226 MTEST((*glt_tcoef_vol)[iglt]);
03227 }
03228
03229 for (ilc = 0; ilc < nlc; ilc++)
03230 {
03231 zero_fill_volume (&((*glt_coef_vol)[iglt][ilc]), nxyz);
03232 if (tout)
03233 zero_fill_volume (&((*glt_tcoef_vol)[iglt][ilc]), nxyz);
03234 }
03235
03236 if (fout) zero_fill_volume (&((*glt_fstat_vol)[iglt]), nxyz);
03237 if (rout) zero_fill_volume (&((*glt_rstat_vol)[iglt]), nxyz);
03238 }
03239 }
03240
03241
03242
03243 if (option_data->fitts_filename != NULL)
03244 {
03245 *fitts_vol = (float **) malloc (sizeof(float **) * nt);
03246 MTEST (*fitts_vol);
03247 for (it = 0; it < nt; it++)
03248 {
03249 zero_fill_volume (&((*fitts_vol)[it]), nxyz);
03250 }
03251 }
03252
03253
03254 if (option_data->errts_filename != NULL)
03255 {
03256 *errts_vol = (float **) malloc (sizeof(float **) * nt);
03257 MTEST (*errts_vol);
03258 for (it = 0; it < nt; it++)
03259 {
03260 zero_fill_volume (&((*errts_vol)[it]), nxyz);
03261 }
03262 }
03263
03264 #ifdef PROC_MAX
03265 if( proc_numjob > 1 ) proc_finalize_shm_volumes() ;
03266 #endif
03267
03268 EXRETURN ;
03269 }
03270
03271
03272
03273
03274
03275
03276
03277 void initialize_program
03278 (
03279 int argc,
03280 char ** argv,
03281 DC_options ** option_data,
03282 THD_3dim_dataset ** dset_time,
03283 byte ** mask_vol,
03284 float ** fmri_data,
03285 int * fmri_length,
03286 float ** censor_array,
03287 int * censor_length,
03288 int ** good_list,
03289 int ** block_list,
03290 int * num_blocks,
03291 float *** stimulus,
03292 int ** stim_length,
03293 matrix ** glt_cmat,
03294
03295 float *** coef_vol,
03296 float *** scoef_vol,
03297 float *** tcoef_vol,
03298 float *** fpart_vol,
03299 float *** rpart_vol,
03300 float ** mse_vol,
03301 float ** ffull_vol,
03302 float ** rfull_vol,
03303
03304 float **** glt_coef_vol,
03305 float **** glt_tcoef_vol,
03306 float *** glt_fstat_vol,
03307 float *** glt_rstat_vol,
03308
03309 float *** fitts_vol,
03310 float *** errts_vol
03311 )
03312
03313 {
03314 int iglt;
03315
03316 ENTRY("initialize_program") ;
03317
03318
03319
03320 *option_data = (DC_options *) malloc (sizeof(DC_options));
03321
03322
03323
03324 get_options (argc, argv, *option_data);
03325
03326
03327
03328 if (!(*option_data)->quiet) identify_software();
03329
03330 if( xrestore ) EXRETURN ;
03331
03332
03333 if( !(*option_data)->quiet &&
03334 !legendre_polort &&
03335 (*option_data)->polort > 1 ){
03336 fprintf(stderr,"** WARNING: you have polynomials of order %d for the baseline\n"
03337 "** but disabled use of the Legendre polynomials!\n"
03338 "** Check the matrix condition and accuracy of results!\n" ,
03339 (*option_data)->polort ) ;
03340 (*option_data)->nocond = 0 ;
03341 }
03342
03343
03344 read_input_data (*option_data, dset_time, mask_vol, fmri_data, fmri_length,
03345 censor_array, censor_length, block_list, num_blocks,
03346 stimulus, stim_length, glt_cmat);
03347
03348
03349
03350 if( !use_psinv )
03351 remove_zero_stimfns (*option_data, *stimulus, *stim_length, *glt_cmat);
03352
03353
03354
03355 check_for_valid_inputs (*option_data, *dset_time,
03356 *fmri_length, *censor_array, *censor_length,
03357 *block_list, *num_blocks, *stim_length, *stimulus, good_list);
03358
03359
03360
03361 if (!(*option_data)->nodata)
03362 allocate_memory (*option_data, coef_vol, scoef_vol, tcoef_vol,
03363 fpart_vol, rpart_vol, mse_vol, ffull_vol, rfull_vol,
03364 glt_coef_vol, glt_tcoef_vol, glt_fstat_vol, glt_rstat_vol,
03365 fitts_vol, errts_vol);
03366
03367 EXRETURN ;
03368 }
03369
03370
03371 #ifndef USE_GET
03372
03373
03374
03375
03376
03377 void extract_ts_array
03378 (
03379 THD_3dim_dataset * dset_time,
03380 int iv,
03381 float * ts_array
03382 )
03383
03384 {
03385 MRI_IMAGE * im;
03386 float * ar;
03387 int ts_length;
03388 int it;
03389
03390
03391
03392 im = THD_extract_series (iv, dset_time, 0);
03393
03394
03395
03396 if (im == NULL) DC_error ("Unable to extract data from 3d+time dataset");
03397
03398
03399
03400 ts_length = DSET_NUM_TIMES (dset_time);
03401 ar = MRI_FLOAT_PTR (im);
03402 #if 0
03403 for (it = 0; it < ts_length; it++)
03404 {
03405 ts_array[it] = ar[it];
03406 }
03407 #else
03408 memcpy( ts_array , ar , sizeof(float)*ts_length ) ;
03409 #endif
03410
03411
03412
03413 mri_free (im); im = NULL;
03414
03415 }
03416 #endif
03417
03418
03419
03420
03421
03422
03423
03424 void save_voxel
03425 (
03426 DC_options * option_data,
03427 int iv,
03428 vector coef,
03429 vector scoef,
03430 vector tcoef,
03431 float * fpart,
03432 float * rpart,
03433 float mse,
03434 float ffull,
03435 float rfull,
03436 vector * glt_coef,
03437 vector * glt_tcoef,
03438 float * fglt,
03439 float * rglt,
03440 int nt,
03441 float * ts_array,
03442 int * good_list,
03443 float * fitts,
03444 float * errts,
03445
03446 float ** coef_vol,
03447 float ** scoef_vol,
03448 float ** tcoef_vol,
03449 float ** fpart_vol,
03450 float ** rpart_vol,
03451 float * mse_vol,
03452 float * ffull_vol,
03453 float * rfull_vol,
03454 float *** glt_coef_vol,
03455 float *** glt_tcoef_vol,
03456 float ** glt_fstat_vol,
03457 float ** glt_rstat_vol,
03458 float ** fitts_vol,
03459 float ** errts_vol
03460
03461 )
03462
03463 {
03464 int ip;
03465 int p;
03466 int is;
03467 int num_stimts;
03468 int num_glt;
03469 int * glt_rows;
03470 int iglt;
03471 int ilc;
03472 int it;
03473 int N;
03474
03475
03476
03477 num_stimts = option_data->num_stimts;
03478 p = option_data->p;
03479 num_glt = option_data->num_glt;
03480 glt_rows = option_data->glt_rows;
03481 N = option_data->N;
03482
03483
03484
03485 if (coef_vol != NULL)
03486 for (ip = 0; ip < p; ip++)
03487 if (coef_vol[ip] != NULL)
03488 coef_vol[ip][iv] = coef.elts[ip];
03489 if (scoef_vol != NULL)
03490 for (ip = 0; ip < p; ip++)
03491 if (scoef_vol[ip] != NULL)
03492 scoef_vol[ip][iv] = scoef.elts[ip];
03493 if (tcoef_vol != NULL)
03494 for (ip = 0; ip < p; ip++)
03495 if (tcoef_vol[ip] != NULL)
03496 tcoef_vol[ip][iv] = tcoef.elts[ip];
03497
03498
03499
03500 if (fpart_vol != NULL)
03501 for (is = 0; is < num_stimts; is++)
03502 if (fpart_vol[is] != NULL)
03503 fpart_vol[is][iv] = fpart[is];
03504 if (rpart_vol != NULL)
03505 for (is = 0; is < num_stimts; is++)
03506 if (rpart_vol[is] != NULL)
03507 rpart_vol[is][iv] = rpart[is];
03508
03509
03510
03511 if (mse_vol != NULL) mse_vol[iv] = mse;
03512
03513
03514
03515 if (ffull_vol != NULL) ffull_vol[iv] = ffull;
03516
03517
03518
03519 if (rfull_vol != NULL) rfull_vol[iv] = rfull;
03520
03521
03522
03523 if (num_glt > 0)
03524 {
03525
03526 if (glt_coef_vol != NULL)
03527 for (iglt = 0; iglt < num_glt; iglt++)
03528 if (glt_coef_vol[iglt] != NULL)
03529 for (ilc = 0; ilc < glt_rows[iglt]; ilc++)
03530 glt_coef_vol[iglt][ilc][iv] = glt_coef[iglt].elts[ilc];
03531
03532
03533 if (glt_tcoef_vol != NULL)
03534 for (iglt = 0; iglt < num_glt; iglt++)
03535 if (glt_tcoef_vol[iglt] != NULL)
03536 for (ilc = 0; ilc < glt_rows[iglt]; ilc++)
03537 glt_tcoef_vol[iglt][ilc][iv] = glt_tcoef[iglt].elts[ilc];
03538
03539
03540 if (glt_fstat_vol != NULL)
03541 for (iglt = 0; iglt < num_glt; iglt++)
03542 if (glt_fstat_vol[iglt] != NULL)
03543 glt_fstat_vol[iglt][iv] = fglt[iglt];
03544
03545
03546 if (glt_rstat_vol != NULL)
03547 for (iglt = 0; iglt < num_glt; iglt++)
03548 if (glt_rstat_vol[iglt] != NULL)
03549 glt_rstat_vol[iglt][iv] = rglt[iglt];
03550 }
03551
03552
03553
03554 if (fitts_vol != NULL)
03555 {
03556 for (it = 0; it < nt; it++)
03557 fitts_vol[it][iv] = ts_array[it];
03558
03559 for (it = 0; it < N; it++)
03560 fitts_vol[good_list[it]][iv] = fitts[it];
03561 }
03562
03563 if (errts_vol != NULL)
03564 {
03565 for (it = 0; it < nt; it++)
03566 errts_vol[it][iv] = 0.0;
03567
03568 for (it = 0; it < N; it++)
03569 errts_vol[good_list[it]][iv] = errts[it];
03570 }
03571
03572 }
03573
03574
03575
03576
03577
03578
03579
03580 void report_evaluation
03581 (
03582 int qp,
03583 int num_stimts,
03584 char ** stim_label,
03585 int * min_lag,
03586 int * max_lag,
03587 matrix x_full ,
03588 matrix xtxinv_full,
03589 int num_glt,
03590 char ** glt_label,
03591 int * glt_rows,
03592 matrix * cxtxinvct
03593 )
03594
03595 {
03596 int m;
03597 int is;
03598 int ilag;
03599 int iglt;
03600 int ilc;
03601 float stddev;
03602 int ibot,itop ;
03603 int j , do_extras ; float jnorm ;
03604
03605
03606
03607 do_extras = AFNI_yesenv("AFNI_3dDeconvolve_nodata_extras") ;
03608 if( do_extras ){
03609 printf("\nBaseline -polort parameters:\n") ;
03610 for( m=0 ; m < qp ; m++ ){
03611 jnorm = 0.0 ;
03612 for( j=0 ; j < x_full.rows ; j++ ) jnorm += SQR(x_full.elts[j][m]) ;
03613 jnorm = sqrt(jnorm) ;
03614 stddev = sqrt ( xtxinv_full.elts[m][m] );
03615 printf (" h[%2d] norm. std. dev. = %8.4f X col. norm = %8.4f\n",
03616 ilag, stddev, jnorm );
03617 }
03618 }
03619
03620
03621 m = qp;
03622 for (is = 0; is < num_stimts; is++)
03623 {
03624 printf ("\nStimulus: %s \n", stim_label[is]);
03625 if( basis_stim[is] != NULL ){ ibot = 0 ; itop = basis_stim[is]->nfunc-1 ; }
03626 else { ibot = min_lag[is] ; itop = max_lag[is] ; }
03627 for (ilag = ibot; ilag <= itop; ilag++)
03628 {
03629 jnorm = 0.0 ;
03630 for( j=0 ; j < x_full.rows ; j++ ) jnorm += SQR(x_full.elts[j][m]) ;
03631 jnorm = sqrt(jnorm) ;
03632 stddev = sqrt ( xtxinv_full.elts[m][m] );
03633 if( do_extras )
03634 printf (" h[%2d] norm. std. dev. = %8.4f X col. norm = %8.4f\n",
03635 ilag, stddev, jnorm );
03636 else
03637 printf (" h[%2d] norm. std. dev. = %8.4f\n",
03638 ilag, stddev );
03639 m++;
03640 }
03641 }
03642
03643
03644 if (num_glt > 0)
03645 {
03646 for (iglt = 0; iglt < num_glt; iglt++)
03647 {
03648 printf ("\nGeneral Linear Test: %s \n", glt_label[iglt]);
03649
03650 for (ilc = 0; ilc < glt_rows[iglt]; ilc++)
03651 {
03652 stddev = sqrt ( 1.0 * cxtxinvct[iglt].elts[ilc][ilc] );
03653 printf (" LC[%d] norm. std. dev. = %8.4f \n",
03654 ilc, stddev);
03655 }
03656 }
03657 }
03658
03659 fflush(stdout);
03660 }
03661
03662
03663
03664
03665
03666
03667
03668 void calculate_results
03669 (
03670 DC_options * option_data,
03671 THD_3dim_dataset * dset,
03672 byte * mask_vol,
03673 float * fmri_data,
03674 int fmri_length,
03675 int * good_list,
03676 int * block_list,
03677 int num_blocks,
03678 float ** stimulus,
03679 int * stim_length,
03680 matrix * glt_cmat,
03681
03682 float ** coef_vol,
03683 float ** scoef_vol,
03684 float ** tcoef_vol,
03685 float ** fpart_vol,
03686 float ** rpart_vol,
03687 float * mse_vol,
03688 float * ffull_vol,
03689 float * rfull_vol,
03690 float *** glt_coef_vol,
03691 float *** glt_tcoef_vol,
03692 float ** glt_fstat_vol,
03693 float ** glt_rstat_vol,
03694 float ** fitts_vol,
03695 float ** errts_vol
03696 )
03697
03698 {
03699 float * ts_array = NULL;
03700
03701 #ifdef USE_GET
03702 int do_get = 0 ;
03703 #endif
03704
03705 int qp;
03706 int q;
03707 int p;
03708 int polort;
03709 int m;
03710 int n;
03711
03712 vector coef;
03713 vector scoef;
03714 vector tcoef;
03715 float * fpart = NULL;
03716 float * rpart = NULL;
03717 float ffull;
03718 float rfull;
03719 float mse;
03720
03721 matrix xdata;
03722 matrix x_full;
03723 matrix xtxinv_full;
03724 matrix xtxinvxt_full;
03725 matrix x_base;
03726 matrix xtxinvxt_base;
03727 matrix * x_rdcd = NULL;
03728 matrix * xtxinvxt_rdcd = NULL;
03729
03730 vector y;
03731
03732 int ixyz;
03733 int nxyz;
03734
03735 int nt;
03736 int N;
03737
03738 int num_stimts;
03739 int * baseline;
03740 int * min_lag;
03741 int * max_lag;
03742 int * nptr;
03743 char ** stim_label;
03744
03745 int i;
03746 int is;
03747 int ilag;
03748 float stddev;
03749 float rms_min;
03750 char * label;
03751 int nodata;
03752 int novar;
03753
03754 int iglt;
03755 int ilc;
03756 int num_glt;
03757 char ** glt_label;
03758 int * glt_rows;
03759 matrix * cxtxinvct = NULL;
03760 matrix * glt_amat = NULL;
03761 vector * glt_coef = NULL;
03762 vector * glt_tcoef = NULL;
03763 float * fglt = NULL;
03764 float * rglt = NULL;
03765
03766 float * fitts = NULL;
03767 float * errts = NULL;
03768
03769 int vstep ;
03770 double ct ;
03771
03772 ENTRY("calculate_results") ;
03773
03774
03775 nodata = option_data->nodata;
03776 nxyz = option_data->nxyz;
03777 nt = option_data->nt;
03778 rms_min = option_data->rms_min;
03779 num_stimts = option_data->num_stimts;
03780 stim_label = option_data->stim_label;
03781 baseline = option_data->stim_base;
03782 min_lag = option_data->stim_minlag;
03783 max_lag = option_data->stim_maxlag;
03784 nptr = option_data->stim_nptr;
03785 num_glt = option_data->num_glt;
03786 glt_label = option_data->glt_label;
03787 glt_rows = option_data->glt_rows;
03788
03789 polort = option_data->polort;
03790 qp = option_data->qp;
03791 q = option_data->q;
03792 p = option_data->p;
03793
03794 N = option_data->N;
03795
03796
03797
03798 matrix_initialize (&xdata);
03799 matrix_initialize (&x_full);
03800 matrix_initialize (&xtxinv_full);
03801 matrix_initialize (&xtxinvxt_full);
03802 matrix_initialize (&x_base);
03803 matrix_initialize (&xtxinvxt_base);
03804 vector_initialize (&coef) ;
03805 vector_initialize (&scoef);
03806 vector_initialize (&tcoef);
03807 vector_initialize (&y);
03808
03809 if (num_stimts > 0)
03810 {
03811 x_rdcd = (matrix *) malloc (sizeof(matrix) * num_stimts);
03812 MTEST (x_rdcd);
03813 xtxinvxt_rdcd = (matrix *) malloc (sizeof(matrix) * num_stimts);
03814 MTEST (xtxinvxt_rdcd);
03815
03816 for (is =0; is < num_stimts; is++)
03817 {
03818 matrix_initialize (&x_rdcd[is]);
03819 matrix_initialize (&xtxinvxt_rdcd[is]);
03820 }
03821 }
03822
03823 if (num_glt > 0)
03824 {
03825 cxtxinvct = (matrix *) malloc (sizeof(matrix) * num_glt);
03826 glt_amat = (matrix *) malloc (sizeof(matrix) * num_glt);
03827 glt_coef = (vector *) malloc (sizeof(vector) * num_glt);
03828 glt_tcoef = (vector *) malloc (sizeof(vector) * num_glt);
03829
03830 for (iglt =0; iglt < num_glt; iglt++)
03831 {
03832 matrix_initialize (&cxtxinvct[iglt]);
03833 matrix_initialize (&glt_amat[iglt]);
03834 vector_initialize (&glt_coef[iglt]);
03835 vector_initialize (&glt_tcoef[iglt]);
03836 }
03837 }
03838
03839
03840
03841 if (num_stimts > 0)
03842 {
03843 fpart = (float *) malloc (sizeof(float) * num_stimts); MTEST (fpart);
03844 rpart = (float *) malloc (sizeof(float) * num_stimts); MTEST (rpart);
03845 }
03846 if (num_glt > 0)
03847 {
03848 fglt = (float *) malloc (sizeof(float) * num_glt); MTEST (fglt);
03849 rglt = (float *) malloc (sizeof(float) * num_glt); MTEST (rglt);
03850 }
03851
03852 if (option_data->input1D_filename == NULL)
03853 {
03854 #ifdef USE_GET
03855 do_get = 1 ;
03856 #else
03857 ts_array = (float *) malloc (sizeof(float) * nt); MTEST (ts_array);
03858 #endif
03859 }
03860 fitts = (float *) malloc (sizeof(float) * nt); MTEST (fitts);
03861 errts = (float *) malloc (sizeof(float) * nt); MTEST (errts);
03862
03863
03864
03865
03866 ct = COX_clock_time() ;
03867
03868 init_indep_var_matrix (p, qp, polort, nt, N, good_list, block_list,
03869 num_blocks, num_stimts, stimulus, stim_length,
03870 min_lag, max_lag, nptr, option_data->stim_base , &xdata);
03871 if (option_data->xout) matrix_sprint ("X matrix:", xdata);
03872
03873 if( option_data->xjpeg_filename != NULL )
03874 JPEG_matrix_gray( xdata , option_data->xjpeg_filename ) ;
03875
03876
03877
03878
03879 { int *iar , k , nerr=0 ;
03880 iar = matrix_check_columns( xdata , 1.e-7 ) ;
03881 if( iar != NULL ){
03882 fprintf(stderr,"** WARNING: Problems with the X matrix columns:\n") ;
03883 for( k=0 ; iar[2*k] >= 0 ; k++ ){
03884 if( iar[2*k+1] >= 0 ){
03885 fprintf(stderr," * Columns %d and %d are nearly collinear!\n",
03886 iar[2*k],iar[2*k+1] ) ; nerr++ ;
03887 } else {
03888 fprintf(stderr," * Column %d is all zeros: %s\n",
03889 iar[2*k] ,
03890 use_psinv ? "SVD on => will be kept"
03891 : "SVD off => will be excised" ) ;
03892 }
03893 }
03894 if( nerr > 0 && AFNI_yesenv("AFNI_3dDeconvolve_nodup") ){
03895 fprintf(stderr,"** ERROR: can't continue after above warnings!\n");
03896 exit(1) ;
03897 }
03898 free(iar) ;
03899 }
03900 }
03901
03902
03903
03904 #ifndef PSINV_EPS
03905 #define PSINV_EPS 1.e-12
03906 #endif
03907 if( !option_data->nocond ){
03908 double *ev , ebot,emin,emax ; int i,nsmall=0 ;
03909 ev = matrix_singvals( xdata ) ;
03910 emax = 1.e-38 ;
03911 for( i=0 ; i < xdata.cols ; i++ ){
03912 if( ev[i] > emax ) emax = ev[i] ;
03913 }
03914 ebot = sqrt(PSINV_EPS)*emax ; emin = 1.e+38 ;
03915 for( i=0 ; i < xdata.cols ; i++ ){
03916 if( ev[i] >= ebot && ev[i] < emin ) emin = ev[i] ;
03917 if( ev[i] < ebot ) nsmall++ ;
03918 }
03919 if( nsmall > 0 ){
03920 fprintf(stderr,
03921 "** WARNING: Largest singular value=%g;"
03922 " %d %s less than cutoff=%g\n"
03923 " Implies strong collinearity in the input regressors\n",
03924 emax , nsmall , (nsmall==1)?"is":"are" , ebot ) ;
03925 show_singvals = 1 ;
03926 }
03927 if( show_singvals ){
03928 fprintf(stderr,"++ Matrix singular values:") ;
03929 for( i=0; i < xdata.cols ; i++ ) fprintf(stderr," %g",ev[i]) ;
03930 fprintf(stderr,"\n") ;
03931 }
03932 free((void *)ev) ;
03933 if( emin <= 0.0 || emax <= 0.0 ){
03934 fprintf(stderr,"** Matrix condition: UNDEFINED: "
03935 "min ev=%g max ev=%g ** VERY BAD **\n",emin,emax ) ;
03936 } else {
03937 double cond = sqrt(emax/emin) ;
03938 if( !use_psinv ) cond = cond*cond ;
03939 fprintf(stderr,"++ (%dx%d) Matrix condition [%s]: %g",
03940 xdata.rows,xdata.cols , (use_psinv) ? "X" : "XtX" , cond ) ;
03941 #ifdef FLOATIZE
03942 if( cond > 100.0 ) fprintf(stderr," ** BEWARE **") ;
03943 #else
03944 if( cond > 1.e6 ) fprintf(stderr," ** BEWARE **") ;
03945 #endif
03946 fprintf(stderr,"\n") ;
03947 }
03948 }
03949
03950
03951 init_regression_analysis (p, qp, num_stimts, baseline, min_lag, max_lag,
03952 xdata, &x_full, &xtxinv_full, &xtxinvxt_full,
03953 &x_base, &xtxinvxt_base, x_rdcd, xtxinvxt_rdcd);
03954 if ((option_data->xout) || nodata)
03955 matrix_sprint ("(X'X) inverse matrix:", xtxinv_full);
03956
03957
03958 X = x_full ;
03959 XtXinv = xtxinv_full ;
03960 XtXinvXt = xtxinvxt_full ;
03961
03962 { int m , npar , j ;
03963 register double sum ;
03964 float mmax ;
03965
03966 Xcol_inbase = (int *) calloc(sizeof(int) ,p) ;
03967 Xcol_mean = (float *)calloc(sizeof(float),p) ;
03968
03969 for( is=0 ; is < qp ; is++ ) Xcol_inbase[is] = 1 ;
03970 m = qp ;
03971 for( is=0 ; is < num_stimts ; is++ ){
03972 npar = (basis_stim[is] != NULL)
03973 ? basis_stim[is]->nfunc
03974 : option_data->stim_maxlag[is] - option_data->stim_minlag[is] + 1 ;
03975
03976 if( baseline[is] ) for( j=0 ; j < npar ; j++ ) Xcol_inbase[m+j] = 1 ;
03977 m += npar ;
03978 }
03979
03980 mmax = 0.0f ;
03981 for( j=0 ; j < p ; j++ ){
03982 sum = 0.0 ;
03983 for( i=0 ; i < X.rows ; i++ ) sum += X.elts[i][j] ;
03984 Xcol_mean[j] = (float)(sum/X.rows) ;
03985 if( Xcol_inbase[j] && fabs(Xcol_mean[j]) > mmax )
03986 mmax = fabs(Xcol_mean[j]) ;
03987 }
03988
03989 if( mmax > 0.0f ){
03990 mmax *= 9.99e-6 ;
03991 for( j=0 ; j < p ; j++ )
03992 if( Xcol_inbase[j] && fabs(Xcol_mean[j]) > mmax ) Xcol_inbase[j] = 2 ;
03993 }
03994 }
03995
03996
03997 if( !option_data->nocond ){
03998 double esum , sum ;
03999 int nn=xdata.rows , mm=xdata.cols , ii,jj,kk ;
04000 char *www = "\0" ;
04001 fprintf(stderr,"++ Matrix inverse average error = ") ;
04002 esum = 0.0 ;
04003 for( ii=0 ; ii < mm ; ii++ ){
04004 for( jj=0 ; jj < mm ; jj++ ){
04005 sum = (ii==jj) ? -1.0 : 0.0 ;
04006 for( kk=0 ; kk < nn ; kk++ )
04007 sum += xtxinvxt_full.elts[ii][kk]*xdata.elts[kk][jj] ;
04008 esum += fabs(sum) ;
04009 }
04010 }
04011 esum /= (mm*mm) ;
04012 if( esum > 1.e-3 ) www = " ** WARNING!!! **" ;
04013 fprintf(stderr,"%g %s\n",esum,www) ;
04014 }
04015
04016
04017 if( option_data->xjpeg_filename != NULL ){
04018 char *jpt , *jsuf=".jpg" ;
04019 char *fn = calloc( sizeof(char) , strlen(option_data->xjpeg_filename)+16 ) ;
04020 matrix xpsinv ;
04021
04022 matrix_initialize( &xpsinv ) ;
04023 matrix_transpose( xtxinvxt_full , &xpsinv ) ;
04024
04025 strcpy(fn,option_data->xjpeg_filename) ;
04026 jpt = strstr(fn,".jpg") ;
04027 if( jpt == NULL ){ jpt = strstr(fn,".JPG") ; jsuf = ".JPG" ; }
04028 if( jpt == NULL ) jpt = fn + strlen(fn) ;
04029 strcpy(jpt,"_psinv") ; strcat(fn,jsuf) ;
04030
04031 JPEG_matrix_gray( xpsinv , fn ) ;
04032 free((void *)fn) ; matrix_destroy( &xpsinv ) ;
04033 }
04034
04035
04036 if (num_glt > 0)
04037 init_glt_analysis (xtxinv_full, num_glt, glt_cmat, glt_amat, cxtxinvct);
04038
04039 ct = COX_clock_time() - ct ;
04040 fprintf(stderr,"++ Matrix setup time = %.2f s\n",ct) ;
04041
04042
04043 vector_create (N, &y);
04044
04045
04046 if (nodata)
04047 {
04048 report_evaluation(qp, num_stimts, stim_label, min_lag, max_lag,
04049 x_full, xtxinv_full,
04050 num_glt, glt_label, glt_rows, cxtxinvct);
04051 }
04052
04053 else
04054 {
04055
04056 int ixyz_bot=0 , ixyz_top=nxyz ;
04057
04058 #ifdef USE_GET
04059 #define NGET 128
04060 int nget=0 ,
04061 cget=0 ,
04062 jget ,
04063 iget[NGET] ;
04064 MRI_IMARR *imget=NULL ;
04065 #endif
04066
04067 #ifdef PROC_MAX
04068 if( proc_numjob > 1 ){
04069 int vv , nvox=nxyz , nper , pp , nv ;
04070 pid_t newpid ;
04071
04072
04073
04074 if( mask_vol != NULL ){
04075 for( vv=nvox=0 ; vv < nxyz ; vv++ )
04076 if( mask_vol[vv] != 0 ) nvox++ ;
04077 }
04078
04079 if( nvox < proc_numjob ){
04080
04081 proc_numjob = 1 ;
04082
04083 } else {
04084
04085
04086
04087 nper = nvox / proc_numjob ;
04088 if( mask_vol == NULL ){
04089 proc_vox_bot[0] = 0 ;
04090 for( pp=0 ; pp < proc_numjob ; pp++ ){
04091 proc_vox_top[pp] = proc_vox_bot[pp] + nper ;
04092 if( pp < proc_numjob-1 ) proc_vox_bot[pp+1] = proc_vox_top[pp] ;
04093 }
04094 proc_vox_top[proc_numjob-1] = nxyz ;
04095 } else {
04096 proc_vox_bot[0] = 0 ;
04097 for( pp=0 ; pp < proc_numjob ; pp++ ){
04098 for( nv=0,vv=proc_vox_bot[pp] ;
04099 nv < nper && vv < nxyz ; vv++ ){
04100 if( mask_vol[vv] != 0 ) nv++ ;
04101 }
04102 proc_vox_top[pp] = vv ;
04103 if( pp < proc_numjob-1 ) proc_vox_bot[pp+1] = proc_vox_top[pp] ;
04104 }
04105 proc_vox_top[proc_numjob-1] = nxyz ;
04106 }
04107
04108
04109
04110 DSET_load(dset) ;
04111
04112
04113
04114 if( !option_data->quiet ) fprintf(stderr,"++ Voxels in dataset: %d\n",nxyz) ;
04115 if( nvox < nxyz )
04116 if( !option_data->quiet ) fprintf(stderr,"++ Voxels in mask: %d\n",nvox) ;
04117 if( !option_data->quiet ) fprintf(stderr,"++ Voxels per job: %d\n",nper) ;
04118
04119 for( pp=1 ; pp < proc_numjob ; pp++ ){
04120 ixyz_bot = proc_vox_bot[pp] ;
04121 ixyz_top = proc_vox_top[pp] ;
04122 proc_ind = pp ;
04123 newpid = fork() ;
04124 if( newpid == -1 ){
04125 fprintf(stderr,"** Can't fork job #%d! Error exit!\n",pp);
04126 exit(1) ;
04127 }
04128 if( newpid == 0 ) break ;
04129 proc_pid[pp] = newpid ;
04130 iochan_sleep(10) ;
04131 }
04132 if( pp == proc_numjob ){
04133 ixyz_bot = proc_vox_bot[0] ;
04134 ixyz_top = proc_vox_top[0] ;
04135 proc_ind = 0 ;
04136 }
04137 if( !option_data->quiet )
04138 fprintf(stderr,"++ Job #%d: processing voxels %d to %d; elapsed time=%.3f\n",
04139 proc_ind,ixyz_bot,ixyz_top-1,COX_clock_time()) ;
04140 }
04141 }
04142 #endif
04143
04144 if( proc_numjob == 1 && !option_data->quiet )
04145 fprintf(stderr,"++ Calculations starting; elapsed time=%.3f\n",COX_clock_time()) ;
04146
04147 vstep = nxyz / 50 ;
04148 if( option_data->quiet ||
04149 option_data->fdisp >= 0.0 ||
04150 option_data->progress > 0 ||
04151 proc_numjob > 1 ) vstep = 0 ;
04152
04153 if( vstep > 0 ) fprintf(stderr,"++ voxel loop:") ;
04154
04155
04156 for (ixyz = ixyz_bot; ixyz < ixyz_top; ixyz++)
04157 {
04158
04159 if( vstep > 0 && ixyz%vstep==vstep-1 ) vstep_print() ;
04160
04161
04162 if (mask_vol != NULL)
04163 if (mask_vol[ixyz] == 0) continue;
04164
04165 #ifdef USE_GET
04166
04167
04168 if( do_get && cget == nget ){
04169 if( imget != NULL ) DESTROY_IMARR(imget) ;
04170 iget[0] = ixyz ; nget = 1 ;
04171 for( jget=ixyz+1 ; jget < nxyz && nget < NGET ; jget++ ){
04172 if( mask_vol == NULL || mask_vol[jget] != 0 )
04173 iget[nget++] = jget ;
04174 }
04175 imget = THD_extract_many_series( nget, iget, dset ) ;
04176 cget = 0 ;
04177 }
04178 #endif
04179
04180
04181 if (option_data->input1D_filename != NULL)
04182 ts_array = fmri_data;
04183 else {
04184 #ifdef USE_GET
04185 ts_array = MRI_FLOAT_PTR(IMARR_SUBIM(imget,cget));
04186 cget++ ;
04187 #else
04188 extract_ts_array (dset, ixyz, ts_array);
04189 #endif
04190 }
04191
04192 for (i = 0; i < N; i++)
04193 y.elts[i] = ts_array[good_list[i]];
04194
04195
04196
04197 regression_analysis (N, p, q, num_stimts, min_lag, max_lag,
04198 x_full, xtxinv_full, xtxinvxt_full, x_base,
04199 xtxinvxt_base, x_rdcd, xtxinvxt_rdcd,
04200 y, rms_min, &mse, &coef, &scoef, &tcoef,
04201 fpart, rpart, &ffull, &rfull, &novar,
04202 fitts, errts);
04203
04204 if( voxel_base != NULL )
04205 voxel_base[ixyz] = baseline_mean( coef ) ;
04206
04207
04208
04209 if (num_glt > 0)
04210 glt_analysis (N, p, x_full, y, mse*(N-p), coef, novar, cxtxinvct,
04211 num_glt, glt_rows, glt_cmat, glt_amat,
04212 glt_coef, glt_tcoef, fglt, rglt);
04213
04214
04215
04216 save_voxel (option_data, ixyz, coef, scoef, tcoef, fpart, rpart, mse,
04217 ffull, rfull, glt_coef, glt_tcoef, fglt, rglt,
04218 nt, ts_array, good_list, fitts, errts,
04219 coef_vol, scoef_vol, tcoef_vol, fpart_vol, rpart_vol,
04220 mse_vol, ffull_vol, rfull_vol, glt_coef_vol,
04221 glt_tcoef_vol, glt_fstat_vol, glt_rstat_vol,
04222 fitts_vol, errts_vol);
04223
04224
04225
04226 if ( ((ffull > option_data->fdisp) && (option_data->fdisp >= 0.0))
04227 || ((option_data->progress > 0)
04228 && (ixyz % option_data->progress == 0))
04229 || (option_data->input1D_filename != NULL) )
04230 {
04231
04232 if( proc_ind == 0 ){
04233 printf ("\n\nResults for Voxel #%d: \n", ixyz);
04234 report_results (N, qp, q, p, polort, block_list, num_blocks,
04235 num_stimts, stim_label, baseline, min_lag, max_lag,
04236 coef, tcoef, fpart, rpart, ffull, rfull, mse,
04237 num_glt, glt_label, glt_rows, glt_coef,
04238 glt_tcoef, fglt, rglt, &label);
04239 printf ("%s \n", label); fflush(stdout);
04240 }
04241 }
04242
04243 }
04244
04245 if( vstep > 0 ) fprintf(stderr,"\n") ;
04246
04247 #ifdef USE_GET
04248 if( do_get ){
04249 if( imget != NULL ) DESTROY_IMARR(imget) ;
04250 ts_array = NULL ;
04251 }
04252 #endif
04253
04254
04255
04256
04257 #ifdef PROC_MAX
04258 if( proc_numjob > 1 ){
04259 if( proc_ind > 0 ){
04260 if( !option_data->quiet )
04261 fprintf(stderr,"++ Job #%d finished; elapsed time=%.3f\n",proc_ind,COX_clock_time()) ;
04262 _exit(0) ;
04263
04264 } else {
04265 int pp ;
04266 if( !option_data->quiet )
04267 fprintf(stderr,"++ Job #0 waiting for children to finish; elapsed time=%.3f\n",COX_clock_time()) ;
04268 for( pp=1 ; pp < proc_numjob ; pp++ )
04269 waitpid( proc_pid[pp] , NULL , 0 ) ;
04270 if( !option_data->quiet )
04271 fprintf(stderr,"++ Job #0 now finishing up; elapsed time=%.3f\n",COX_clock_time()) ;
04272 }
04273
04274
04275
04276 }
04277 #endif
04278 if( proc_numjob == 1 && !option_data->quiet )
04279 fprintf(stderr,"++ Calculations finished; elapsed time=%.3f\n",COX_clock_time()) ;
04280
04281 if( option_data->input1D_filename == NULL)
04282 DSET_unload(dset) ;
04283
04284 }
04285
04286
04287
04288 vector_destroy (&y);
04289 vector_destroy (&tcoef);
04290 vector_destroy (&scoef);
04291 vector_destroy (&coef);
04292
04293 if (num_stimts > 0)
04294 {
04295 for (is = 0; is < num_stimts; is++)
04296 {
04297 matrix_destroy (&x_rdcd[is]);
04298 matrix_destroy (&xtxinvxt_rdcd[is]);
04299 }
04300 free (x_rdcd); x_rdcd = NULL;
04301 free (xtxinvxt_rdcd); xtxinvxt_rdcd = NULL;
04302 }
04303
04304 matrix_destroy (&xtxinvxt_base);
04305 matrix_destroy (&x_base);
04306 matrix_destroy (&xdata);
04307
04308 if (num_glt > 0)
04309 {
04310 for (iglt = 0; iglt < num_glt; iglt++)
04311 {
04312 matrix_destroy (&cxtxinvct[iglt]);
04313 matrix_destroy (&glt_amat[iglt]);
04314 vector_destroy (&glt_coef[iglt]);
04315 vector_destroy (&glt_tcoef[iglt]);
04316 }
04317 free (cxtxinvct); cxtxinvct = NULL;
04318 free (glt_amat); glt_amat = NULL;
04319 free (glt_coef); glt_coef = NULL;
04320 free (glt_tcoef); glt_tcoef = NULL;
04321 }
04322
04323
04324 if (fpart != NULL) { free (fpart); fpart = NULL; }
04325 if (rpart != NULL) { free (rpart); rpart = NULL; }
04326 if (fglt != NULL) { free (fglt); fglt = NULL; }
04327 if (rglt != NULL) { free (rglt); rglt = NULL; }
04328 if (ts_array != NULL) { free (ts_array); ts_array = NULL; }
04329 if (fitts != NULL) { free (fitts); fitts = NULL; }
04330 if (errts != NULL) { free (errts); errts = NULL; }
04331
04332 EXRETURN ;
04333 }
04334
04335
04336
04337
04338
04339
04340
04341
04342
04343
04344
04345
04346
04347 float EDIT_coerce_autoscale_new( int nxyz ,
04348 int itype,void *ivol , int otype,void *ovol )
04349 {
04350 float fac=0.0 , top ;
04351
04352 ENTRY("EDIT_coerce_autoscale_new") ;
04353
04354 if( MRI_IS_INT_TYPE(otype) ){
04355 top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
04356 if (top == 0.0) fac = 0.0;
04357 else fac = MRI_TYPE_maxval[otype]/top;
04358 }
04359
04360 EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
04361 RETURN ( fac );
04362 }
04363
04364
04365
04366
04367
04368
04369
04370
04371
04372 void cubic_spline
04373 (
04374 DC_options * option_data,
04375 int ts_length,
04376 float ** vol_array
04377 )
04378
04379 {
04380 THD_3dim_dataset * dset = NULL;
04381 int nx, ny, nz, nxyz;
04382 int ixyz;
04383 int isl;
04384 int i;
04385 float * yarray = NULL;
04386 float * sarray = NULL;
04387 matrix m, minv;
04388 vector v, sv;
04389 int n;
04390 float * a = NULL,
04391 * b = NULL,
04392 * c = NULL,
04393 * d = NULL;
04394 float tslice;
04395 float tdelta;
04396 float frac;
04397 int k;
04398 float t;
04399 float delt;
04400 float y;
04401
04402
04403
04404 matrix_initialize (&m);
04405 matrix_initialize (&minv);
04406 vector_initialize (&v);
04407 vector_initialize (&sv);
04408
04409
04410
04411 dset = THD_open_dataset (option_data->input_filename);
04412 n = ts_length - 1;
04413 tdelta = dset->taxis->ttdel;
04414 nx = dset->daxes->nxx; ny = dset->daxes->nyy; nz = dset->daxes->nzz;
04415 nxyz = nx * ny * nz;
04416
04417
04418
04419 yarray = (float *) malloc (sizeof(float) * ts_length);
04420 sarray = (float *) malloc (sizeof(float) * (n+1));
04421 a = (float *) malloc (sizeof(float) * n);
04422 b = (float *) malloc (sizeof(float) * n);
04423 c = (float *) malloc (sizeof(float) * n);
04424 d = (float *) malloc (sizeof(float) * n);
04425
04426
04427
04428 matrix_create (n-1, n-1, &m);
04429 m.elts[0][0] = 4.0;
04430 m.elts[0][1] = 1.0;
04431 m.elts[n-2][n-3] = 1.0;
04432 m.elts[n-2][n-2] = 4.0;
04433 for (i = 1; i < n-2; i++)
04434 {
04435 m.elts[i][i] = 4.0;
04436 m.elts[i][i-1] = 1.0;
04437 m.elts[i][i+1] = 1.0;
04438 }
04439 matrix_inverse (m, &minv);
04440
04441
04442 vector_create (n-1, &v);
04443
04444
04445
04446 for (ixyz = 0; ixyz < nxyz; ixyz++)
04447 {
04448
04449
04450 isl = ixyz / (nx*ny);
04451 tslice = THD_timeof_slice (0, isl, dset);
04452 frac = -((tslice/tdelta) - 0.5);
04453
04454
04455 for (i = 0; i < ts_length; i++)
04456 yarray[i] = vol_array[i][ixyz];
04457
04458
04459
04460 for (i = 1; i < n; i++)
04461 v.elts[i-1] = 6.0 * (yarray[i+1] - 2.0 * yarray[i] + yarray[i-1]);
04462 vector_multiply (minv, v, &sv);
04463
04464
04465
04466 for (i = 1; i < n; i++)
04467 {
04468 sarray[i] = sv.elts[i-1];
04469 }
04470 sarray[0] = 0.0;
04471 sarray[n] = 0.0;
04472
04473
04474
04475 for (i = 0; i < n; i++)
04476 {
04477 a[i] = (sarray[i+1] - sarray[i]) / 6.0;
04478 b[i] = sarray[i] / 2;
04479 c[i] = (yarray[i+1]-yarray[i]) - (2.0*sarray[i]+sarray[i+1]) / 6.0;
04480 d[i] = yarray[i];
04481 }
04482
04483
04484
04485 for (i = 0; i < ts_length; i++)
04486 {
04487 t = i + frac;
04488
04489 if (frac < 0.0) k = i-1;
04490 else k = i;
04491
04492 if (k < 0) k = 0;
04493 if (k > n-1) k = n-1;
04494
04495 delt = t - k;
04496
04497 yarray[i] = a[k]*delt*delt*delt + b[k]*delt*delt + c[k]*delt + d[k];
04498 }
04499
04500
04501
04502 for (i = 0; i < ts_length; i++)
04503 vol_array[i][ixyz] = yarray[i];
04504
04505
04506
04507 }
04508
04509
04510
04511 THD_delete_3dim_dataset (dset, False); dset = NULL ;
04512
04513 matrix_destroy (&m);
04514 matrix_destroy (&minv);
04515 vector_destroy (&v);
04516 vector_destroy (&sv);
04517
04518 free (sarray); sarray = NULL;
04519 free (yarray); yarray = NULL;
04520 free (a); a = NULL;
04521 free (b); b = NULL;
04522 free (c); c = NULL;
04523 free (d); d = NULL;
04524 }
04525
04526
04527
04528
04529
04530
04531
04532
04533 void write_ts_array
04534 (
04535 int argc,
04536 char ** argv,
04537 DC_options * option_data,
04538 int ts_length,
04539 int nptr,
04540 int tshift,
04541 float ** vol_array,
04542 char * output_filename
04543 )
04544
04545 {
04546 const float EPSILON = 1.0e-10;
04547
04548 THD_3dim_dataset * dset = NULL;
04549 THD_3dim_dataset * new_dset = NULL;
04550 int ib;
04551 int ierror;
04552 int nxyz;
04553 float factor;
04554 char * input_filename;
04555 short ** bar = NULL;
04556 float * fbuf = NULL;
04557 float * volume;
04558 char label[THD_MAX_NAME];
04559 float newtr;
04560
04561
04562
04563 input_filename = option_data->input_filename;
04564 dset = THD_open_dataset (input_filename);
04565 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
04566 newtr = DSET_TIMESTEP(dset) / nptr;
04567
04568
04569 bar = (short **) malloc (sizeof(short *) * ts_length); MTEST (bar);
04570 fbuf = (float *) malloc (sizeof(float) * ts_length); MTEST (fbuf);
04571
04572
04573 new_dset = EDIT_empty_copy (dset);
04574
04575
04576 tross_Copy_History( dset , new_dset ) ;
04577
04578 { char * commandline = tross_commandline( PROGRAM_NAME , argc , argv ) ;
04579 sprintf (label, "Output prefix: %s", output_filename);
04580 if( commandline != NULL )
04581 tross_multi_Append_History( new_dset , commandline,label,NULL ) ;
04582 else
04583 tross_Append_History ( new_dset, label);
04584 free(commandline) ;
04585 }
04586
04587
04588 THD_delete_3dim_dataset (dset, False); dset = NULL ;
04589
04590
04591 ierror = EDIT_dset_items (new_dset,
04592 ADN_prefix, output_filename,
04593 ADN_label1, output_filename,
04594 ADN_self_name, output_filename,
04595 ADN_malloc_type, DATABLOCK_MEM_MALLOC,
04596 ADN_datum_all, MRI_short,
04597 ADN_nvals, ts_length,
04598 ADN_ntt, ts_length,
04599 ADN_ttdel, newtr,
04600 ADN_none);
04601
04602 if( ierror > 0 ){
04603 fprintf(stderr,
04604 "** %d errors in attempting to create output dataset!\n", ierror ) ;
04605 exit(1) ;
04606 }
04607
04608 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
04609 fprintf(stderr,
04610 "** Output dataset file %s already exists--cannot continue!\a\n",
04611 new_dset->dblk->diskptr->header_name ) ;
04612 exit(1) ;
04613 }
04614
04615
04616
04617 if (tshift)
04618 EDIT_dset_items (new_dset,
04619 ADN_nsl, 0,
04620 ADN_ttorg, 0.0,
04621 ADN_ttdur, 0.0,
04622 ADN_none);
04623
04624
04625
04626 for (ib = 0; ib < ts_length; ib++)
04627 {
04628
04629
04630 volume = vol_array[ib];
04631
04632
04633 bar[ib] = (short *) malloc (sizeof(short) * nxyz);
04634 MTEST (bar[ib]);
04635
04636
04637 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
04638 MRI_short, bar[ib]);
04639 if (factor < EPSILON) factor = 0.0;
04640 else factor = 1.0 / factor;
04641 fbuf[ib] = factor;
04642
04643
04644 mri_fix_data_pointer (bar[ib], DSET_BRICK(new_dset,ib));
04645 }
04646
04647
04648
04649
04650 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
04651
04652 THD_load_statistics (new_dset);
04653 THD_write_3dim_dataset (NULL, NULL, new_dset, True);
04654 if (! option_data->quiet)
04655 fprintf(stderr,"++ Wrote 3D+time dataset into %s\n",DSET_BRIKNAME(new_dset)) ;
04656
04657
04658
04659 THD_delete_3dim_dataset (new_dset, False); new_dset = NULL ;
04660 free (fbuf); fbuf = NULL;
04661
04662 }
04663
04664
04665
04666
04667
04668
04669
04670 void attach_sub_brick
04671 (
04672 THD_3dim_dataset * new_dset,
04673 int ibrick,
04674 float * volume,
04675 int nxyz,
04676 int brick_type,
04677 char * brick_label,
04678 int dof,
04679 int ndof,
04680 int ddof,
04681 short ** bar
04682 )
04683
04684 {
04685 const float EPSILON = 1.0e-10;
04686 float factor;
04687 short *sbr ;
04688
04689 ENTRY("attach_sub_brick") ;
04690
04691
04692
04693 sbr = (short *) malloc (sizeof(short) * nxyz);
04694 MTEST (sbr);
04695 factor = EDIT_coerce_autoscale_new(nxyz, MRI_float,volume, MRI_short,sbr);
04696
04697 if (factor < EPSILON) factor = 0.0;
04698 else factor = 1.0 / factor;
04699
04700
04701 EDIT_BRICK_LABEL (new_dset, ibrick, brick_label);
04702 EDIT_BRICK_FACTOR(new_dset, ibrick, factor);
04703
04704 if (brick_type == FUNC_TT_TYPE)
04705 EDIT_BRICK_TO_FITT (new_dset, ibrick, dof);
04706 else if (brick_type == FUNC_FT_TYPE)
04707 EDIT_BRICK_TO_FIFT (new_dset, ibrick, ndof, ddof);
04708
04709
04710
04711 EDIT_substitute_brick (new_dset, ibrick, MRI_short, sbr);
04712 if( bar != NULL ) bar[ibrick] = sbr ;
04713
04714 EXRETURN ;
04715 }
04716
04717
04718
04719
04720
04721
04722 void write_bucket_data
04723 (
04724 int argc,
04725 char ** argv,
04726 DC_options * option_data,
04727
04728 float ** coef_vol,
04729 float ** tcoef_vol,
04730 float ** fpart_vol,
04731 float ** rpart_vol,
04732 float * mse_vol,
04733 float * ffull_vol,
04734 float * rfull_vol,
04735 float *** glt_coef_vol,
04736 float *** glt_tcoef_vol,
04737 float ** glt_fstat_vol,
04738 float ** glt_rstat_vol
04739 )
04740
04741 {
04742 THD_3dim_dataset * old_dset = NULL;
04743 THD_3dim_dataset * new_dset = NULL;
04744 char output_prefix[THD_MAX_NAME];
04745 char output_session[THD_MAX_NAME];
04746 int nbricks;
04747 short ** bar = NULL;
04748
04749 int brick_type;
04750 int brick_coef;
04751 char brick_label[THD_MAX_NAME];
04752
04753 int ierror;
04754 float * volume;
04755
04756 int N;
04757 int qp;
04758 int q;
04759 int p;
04760 int polort;
04761 int num_stimts;
04762 int istim;
04763 int nxyz;
04764 int nt;
04765 int ilag;
04766 int icoef;
04767 int ibrick;
04768 int dof, ndof, ddof;
04769 char label[THD_MAX_NAME];
04770 char blab[THD_MAX_NAME] ;
04771 int num_glt;
04772 int * glt_rows;
04773 int iglt;
04774 int ilc;
04775
04776 THD_3dim_dataset *coef_dset = NULL ;
04777 int cbuck , bout,cout , ibot,itop ;
04778
04779
04780 old_dset = THD_open_dataset (option_data->input_filename);
04781
04782 bout = !option_data->nobout ;
04783 cout = !option_data->nocout ;
04784
04785
04786 nxyz = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;
04787 num_stimts = option_data->num_stimts;
04788 nt = DSET_NUM_TIMES (old_dset);
04789 num_glt = option_data->num_glt;
04790 glt_rows = option_data->glt_rows;
04791
04792 polort = option_data->polort;
04793 qp = option_data->qp;
04794 q = option_data->q;
04795 p = option_data->p;
04796 N = option_data->N;
04797 nbricks = option_data->nbricks;
04798
04799
04800
04801 strcpy (output_prefix, option_data->bucket_filename);
04802 strcpy (output_session, "./");
04803
04804
04805
04806 bar = (short **) malloc (sizeof(short *) * nbricks);
04807 MTEST (bar);
04808
04809
04810
04811 new_dset = EDIT_empty_copy (old_dset);
04812
04813
04814
04815 tross_Copy_History( old_dset , new_dset ) ;
04816 tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
04817 sprintf (label, "Output prefix: %s", output_prefix);
04818 tross_Append_History ( new_dset, label);
04819
04820
04821
04822 ierror = EDIT_dset_items (new_dset,
04823 ADN_prefix, output_prefix,
04824 ADN_type, HEAD_FUNC_TYPE,
04825 ADN_func_type, FUNC_BUCK_TYPE,
04826 ADN_datum_all, MRI_short ,
04827 ADN_ntt, 0,
04828 ADN_nvals, nbricks,
04829 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
04830 ADN_none ) ;
04831
04832 if( ierror > 0 )
04833 {
04834 fprintf(stderr,
04835 "** %d errors in attempting to create bucket dataset!\n",
04836 ierror);
04837 exit(1);
04838 }
04839
04840 if( strstr(output_prefix,"/") == NULL )
04841 (void) EDIT_dset_items (new_dset,
04842 ADN_directory_name, output_session,
04843 ADN_none ) ;
04844
04845 if (THD_is_file(DSET_HEADNAME(new_dset)))
04846 {
04847 fprintf(stderr,
04848 "** Bucket dataset file %s already exists--cannot continue!\n",
04849 DSET_HEADNAME(new_dset));
04850 exit(1);
04851 }
04852
04853 if( CoefFilename != NULL ){
04854 coef_dset = EDIT_empty_copy( new_dset ) ;
04855 tross_Copy_History( old_dset , coef_dset ) ;
04856 tross_Make_History( PROGRAM_NAME , argc , argv , coef_dset ) ;
04857 (void) EDIT_dset_items( coef_dset,
04858 ADN_prefix, CoefFilename ,
04859 ADN_type, HEAD_FUNC_TYPE,
04860 ADN_func_type, FUNC_BUCK_TYPE,
04861 ADN_datum_all, MRI_short ,
04862 ADN_ntt, 0,
04863 ADN_nvals, p ,
04864 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
04865 ADN_none ) ;
04866 if( strstr(CoefFilename,"/") == NULL )
04867 (void) EDIT_dset_items( coef_dset ,
04868 ADN_directory_name, output_session,
04869 ADN_none ) ;
04870 if( THD_is_file(DSET_HEADNAME(coef_dset)) ){
04871 fprintf(stderr,
04872 "** Coefficient dataset file %s already exists--cannot continue!\n",
04873 DSET_HEADNAME(coef_dset));
04874 exit(1);
04875 }
04876 }
04877
04878
04879 if( DSET_IS_TCAT(old_dset) )
04880 InputFilename = strdup(old_dset->tcat_list) ;
04881 else
04882 InputFilename = strdup(THD_trailname(DSET_HEADNAME(old_dset),0)) ;
04883
04884 BucketFilename = strdup(THD_trailname(DSET_HEADNAME(new_dset),0)) ;
04885
04886 if( coef_dset != NULL )
04887 CoefFilename = strdup(THD_trailname(DSET_HEADNAME(coef_dset),0)) ;
04888
04889
04890 THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ;
04891
04892
04893
04894
04895
04896
04897 ibrick = -1;
04898 if (option_data->full_first)
04899 ibrick += option_data->vout + option_data->rout + option_data->fout;
04900
04901 cbuck = -1 ;
04902
04903
04904
04905 if( cout || coef_dset != NULL )
04906 {
04907
04908 { int ii ;
04909 ParamStim = (int *) calloc(sizeof(int) ,nParam) ;
04910 ParamLabel = (char **)calloc(sizeof(char *),nParam) ;
04911 for( ii=0 ; ii < qp ; ii++ ) ParamLabel[ii] = strdup("ort") ;
04912 for( ; ii < nParam ; ii++ ) ParamLabel[ii] = strdup("coef") ;
04913 if( cout && bout )
04914 ParamIndex = (int *)calloc(sizeof(int) ,nParam) ;
04915 else
04916 ParamIndex = NULL ;
04917 }
04918
04919
04920 if( bout || coef_dset != NULL )
04921 {
04922 strcpy (label, "Base");
04923
04924 if( legendre_polort ) strcpy(blab,"P_") ;
04925 else strcpy(blab,"t^") ;
04926
04927 for (icoef = 0; icoef < qp; icoef++)
04928 {
04929 if (qp == polort+1)
04930 strcpy (label, "Base");
04931 else
04932 sprintf (label, "Run#%d", icoef/(polort+1) + 1);
04933
04934
04935 brick_type = FUNC_FIM_TYPE;
04936 sprintf (brick_label, "%s %s%d Coef", label,blab, icoef % (polort+1));
04937 volume = coef_vol[icoef];
04938
04939 if( cout && bout )
04940 attach_sub_brick (new_dset, ++ibrick, volume, nxyz,
04941 brick_type, brick_label, 0, 0, 0, bar);
04942
04943 sprintf(brick_label,"%s:%s%d" , label,blab,icoef%(polort+1)) ;
04944
04945 if( ParamIndex != NULL ) ParamIndex[icoef] = ibrick ;
04946 ParamStim [icoef] = 0 ;
04947 ParamLabel[icoef] = strdup(brick_label) ;
04948
04949 if( coef_dset != NULL ){
04950 cbuck++ ;
04951 attach_sub_brick( coef_dset , cbuck , volume , nxyz ,
04952 brick_type, brick_label, 0, 0, 0, NULL);
04953 }
04954
04955
04956 if ( cout && bout && option_data->tout)
04957 {
04958 ibrick++;
04959 brick_type = FUNC_TT_TYPE;
04960 dof = N - p;
04961 sprintf (brick_label, "%s %s%d t-st",
04962 label,blab, icoef % (polort+1));
04963 volume = tcoef_vol[icoef];
04964 attach_sub_brick (new_dset, ibrick, volume, nxyz,
04965 brick_type, brick_label, dof, 0, 0, bar);
04966 }
04967 }
04968 }
04969
04970
04971
04972 icoef = qp;
04973 for (istim = 0; istim < num_stimts; istim++)
04974 {
04975 strcpy (label, option_data->stim_label[istim]);
04976
04977 if( basis_stim[istim] != NULL ){
04978 ibot = 0 ; itop = basis_stim[istim]->nfunc-1 ;
04979 } else {
04980 ibot = option_data->stim_minlag[istim] ;
04981 itop = option_data->stim_maxlag[istim] ;
04982 }
04983
04984
04985 for( ilag=ibot ; ilag <= itop ; ilag++ )
04986 {
04987 if( !option_data->stim_base[istim] ||
04988 bout ||
04989 coef_dset != NULL )
04990 {
04991
04992 brick_type = FUNC_FIM_TYPE;
04993 sprintf (brick_label, "%s[%d] Coef", label, ilag);
04994 volume = coef_vol[icoef];
04995 if( cout && (!option_data->stim_base[istim] || bout) )
04996 attach_sub_brick (new_dset, ++ibrick, volume, nxyz,
04997 brick_type, brick_label, 0, 0, 0, bar);
04998
04999 sprintf(brick_label,"%s:%d",label,ilag) ;
05000
05001 if( ParamIndex != NULL ) ParamIndex[icoef] = ibrick ;
05002 ParamStim [icoef] = option_data->stim_base[istim] ? 0
05003 : istim+1 ;
05004 ParamLabel[icoef] = strdup(brick_label) ;
05005
05006 if( coef_dset != NULL ){
05007 cbuck++ ;
05008 attach_sub_brick( coef_dset , cbuck , volume , nxyz ,
05009 brick_type, brick_label, 0, 0, 0, NULL);
05010 }
05011
05012
05013 if (cout && option_data->tout && (!option_data->stim_base[istim] || bout) )
05014 {
05015 ibrick++;
05016 brick_type = FUNC_TT_TYPE;
05017 dof = N - p;
05018 sprintf (brick_label, "%s[%d] t-st", label, ilag);
05019 volume = tcoef_vol[icoef];
05020 attach_sub_brick (new_dset, ibrick, volume, nxyz,
05021 brick_type, brick_label, dof, 0, 0, bar);
05022 }
05023 }
05024
05025 icoef++;
05026 }
05027
05028
05029 if( cout && option_data->rout
05030 && (!option_data->stim_base[istim] || bout) )
05031 {
05032 ibrick++;
05033 brick_type = FUNC_THR_TYPE;
05034 sprintf (brick_label, "%s R^2", label);
05035 volume = rpart_vol[istim];
05036 attach_sub_brick (new_dset, ibrick, volume, nxyz,
05037 brick_type, brick_label, 0, 0, 0, bar);
05038 }
05039
05040
05041 if( cout && option_data->fout
05042 && (!option_data->stim_base[istim] || bout) )
05043 {
05044 ibrick++;
05045 brick_type = FUNC_FT_TYPE;
05046 ndof = itop - ibot + 1 ;
05047 ddof = N - p;
05048 sprintf (brick_label, "%s F-stat", label);
05049 volume = fpart_vol[istim];
05050 attach_sub_brick (new_dset, ibrick, volume, nxyz,
05051 brick_type, brick_label, 0, ndof, ddof, bar);
05052 }
05053
05054 }
05055
05056 }
05057
05058 if( coef_dset != NULL ){
05059 DSET_write(coef_dset) ;
05060 if( !option_data->quiet )
05061 fprintf(stderr,"++ Wrote cbucket to %s\n",DSET_BRIKNAME(coef_dset)) ;
05062 DSET_delete(coef_dset) ; coef_dset = NULL ;
05063 }
05064
05065
05066
05067
05068 if( xsave ) XSAVE_output( DSET_PREFIX(new_dset) ) ;
05069
05070
05071 for (iglt = 0; iglt < num_glt; iglt++)
05072 {
05073 strcpy (label, option_data->glt_label[iglt]);
05074
05075
05076 for (ilc = 0; ilc < glt_rows[iglt]; ilc++)
05077 {
05078
05079 ibrick++;
05080 brick_type = FUNC_FIM_TYPE;
05081 sprintf (brick_label, "%s LC[%d] coef", label, ilc);
05082 volume = glt_coef_vol[iglt][ilc];
05083 attach_sub_brick (new_dset, ibrick, volume, nxyz,
05084 brick_type, brick_label, 0, 0, 0, bar);
05085
05086
05087 if (option_data->tout)
05088 {
05089 ibrick++;
05090 brick_type = FUNC_TT_TYPE;
05091 dof = N - p;
05092 sprintf (brick_label, "%s LC[%d] t-st", label, ilc);
05093 volume = glt_tcoef_vol[iglt][ilc];
05094 attach_sub_brick (new_dset, ibrick, volume, nxyz,
05095 brick_type, brick_label, dof, 0, 0, bar);
05096 }
05097 }
05098
05099
05100 if (option_data->rout)
05101 {
05102 ibrick++;
05103 brick_type = FUNC_THR_TYPE;
05104 sprintf (brick_label, "%s R^2", label);
05105 volume = glt_rstat_vol[iglt];
05106 attach_sub_brick (new_dset, ibrick, volume, nxyz,
05107 brick_type, brick_label, 0, 0, 0, bar);
05108 }
05109
05110
05111 if (option_data->fout)
05112 {
05113 ibrick++;
05114 brick_type = FUNC_FT_TYPE;
05115 ndof = glt_rows[iglt];
05116 ddof = N - p;
05117 sprintf (brick_label, "%s F-stat", label);
05118 volume = glt_fstat_vol[iglt];
05119 attach_sub_brick (new_dset, ibrick, volume, nxyz,
05120 brick_type, brick_label, 0, ndof, ddof, bar);
05121 }
05122
05123 }
05124
05125
05126
05127 if (option_data->full_first) ibrick = -1;
05128
05129
05130 if (option_data->vout)
05131 {
05132 ibrick++;
05133 brick_type = FUNC_FIM_TYPE;
05134 sprintf (brick_label, "Full MSE");
05135 volume = mse_vol;
05136 attach_sub_brick (new_dset, ibrick, volume, nxyz,
05137 brick_type, brick_label, 0, 0, 0, bar);
05138 }
05139
05140
05141 if (option_data->rout)
05142 {
05143 ibrick++;
05144 brick_type = FUNC_THR_TYPE;
05145 sprintf (brick_label, "Full R^2");
05146 volume = rfull_vol;
05147 attach_sub_brick (new_dset, ibrick, volume, nxyz,
05148 brick_type, brick_label, 0, 0, 0, bar);
05149 }
05150
05151
05152 if (option_data->fout)
05153 {
05154 ibrick++;
05155 brick_type = FUNC_FT_TYPE;
05156 ndof = p - q;
05157 ddof = N - p;
05158 sprintf (brick_label, "Full F-stat");
05159 volume = ffull_vol;
05160 attach_sub_brick (new_dset, ibrick, volume, nxyz,
05161 brick_type, brick_label, 0, ndof, ddof, bar);
05162 }
05163
05164
05165
05166
05167 THD_load_statistics (new_dset);
05168 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
05169 if (! option_data->quiet)
05170 printf("++ Wrote bucket dataset into %s\n", DSET_BRIKNAME(new_dset));
05171
05172
05173
05174 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
05175
05176 }
05177
05178
05179
05180
05181
05182
05183
05184 void write_one_ts
05185 (
05186 char * prefix,
05187 int ts_length,
05188 float ** vol_array
05189 )
05190
05191 {
05192 char filename[80];
05193 int it;
05194 FILE * outfile = NULL;
05195
05196
05197
05198 sprintf (filename, "%s.1D", prefix);
05199 outfile = fopen (filename, "w");
05200
05201
05202
05203 for (it = 0; it < ts_length; it++)
05204 {
05205 fprintf (outfile, "%f", vol_array[it][0]);
05206 fprintf (outfile, " \n");
05207 }
05208
05209
05210 fclose (outfile);
05211 }
05212
05213
05214
05215
05216
05217
05218
05219 void output_results
05220 (
05221 int argc,
05222 char ** argv,
05223 DC_options * option_data,
05224
05225 float ** coef_vol,
05226 float ** scoef_vol,
05227 float ** tcoef_vol,
05228 float ** fpart_vol,
05229 float ** rpart_vol,
05230 float * mse_vol,
05231 float * ffull_vol,
05232 float * rfull_vol,
05233 float *** glt_coef_vol,
05234 float *** glt_tcoef_vol,
05235 float ** glt_fstat_vol,
05236 float ** glt_rstat_vol,
05237 float ** fitts_vol,
05238 float ** errts_vol
05239 )
05240
05241 {
05242 int qp;
05243 int q;
05244 int p;
05245 int polort;
05246 int num_stimts;
05247 int * min_lag;
05248 int * max_lag;
05249 int * nptr;
05250 int ib;
05251 int is;
05252 int ts_length;
05253 int nt;
05254 int nxyz;
05255
05256
05257
05258 polort = option_data->polort;
05259 qp = option_data->qp;
05260 q = option_data->q;
05261 p = option_data->p;
05262 num_stimts = option_data->num_stimts;
05263 min_lag = option_data->stim_minlag;
05264 max_lag = option_data->stim_maxlag;
05265 nptr = option_data->stim_nptr;
05266 nt = option_data->nt;
05267 nxyz = option_data->nxyz;
05268
05269
05270
05271 if (option_data->bucket_filename != NULL)
05272 if (nxyz > 1)
05273 write_bucket_data (argc, argv, option_data, coef_vol, tcoef_vol,
05274 fpart_vol, rpart_vol, mse_vol, ffull_vol, rfull_vol,
05275 glt_coef_vol, glt_tcoef_vol, glt_fstat_vol, glt_rstat_vol);
05276 else
05277 write_one_ts (option_data->bucket_filename, p, coef_vol);
05278
05279
05280
05281 ib = qp;
05282 for (is = 0; is < num_stimts; is++)
05283 {
05284 if( basis_stim[is] != NULL ){
05285 if( option_data->iresp_filename[is] != NULL )
05286 #if 0
05287 fprintf(stderr,
05288 "** WARNING: Ignoring -iresp %d '%s'\n",
05289 is+1,option_data->iresp_filename[is]) ;
05290 #else
05291 basis_write_iresp( argc , argv , option_data ,
05292 basis_stim[is] , basis_dtout ,
05293 coef_vol+ib ,
05294 option_data->iresp_filename[is] ) ;
05295 ib += basis_stim[is]->nfunc ;
05296 #endif
05297 continue ;
05298 }
05299
05300
05301
05302 ts_length = max_lag[is] - min_lag[is] + 1;
05303 if (option_data->iresp_filename[is] != NULL)
05304 {
05305
05306 if ((option_data->tshift) && (nptr[is] == 1) && (nxyz > 1))
05307 cubic_spline (option_data, ts_length, coef_vol+ib);
05308
05309 if (nxyz > 1)
05310 write_ts_array (argc, argv, option_data, ts_length,
05311 nptr[is], option_data->tshift, coef_vol+ib,
05312 option_data->iresp_filename[is]);
05313 else
05314 write_one_ts (option_data->iresp_filename[is],
05315 ts_length, coef_vol+ib);
05316 }
05317 ib += ts_length;
05318 }
05319
05320
05321
05322 ib = qp;
05323 for (is = 0; is < num_stimts; is++)
05324 {
05325 if( basis_stim[is] != NULL ){
05326 if( option_data->sresp_filename[is] != NULL )
05327 #if 0
05328 fprintf(stderr,
05329 "** WARNING: Ignoring -sresp %d '%s'\n",
05330 is+1 , option_data->sresp_filename[is]) ;
05331 #else
05332 basis_write_sresp( argc , argv , option_data ,
05333 basis_stim[is] , basis_dtout ,
05334 mse_vol , ib , XtXinv ,
05335 option_data->sresp_filename[is] ) ;
05336 ib += basis_stim[is]->nfunc ;
05337 #endif
05338 continue ;
05339 }
05340
05341
05342
05343 ts_length = max_lag[is] - min_lag[is] + 1;
05344 if (option_data->sresp_filename[is] != NULL)
05345 if (nxyz > 1)
05346 write_ts_array (argc, argv, option_data, ts_length,
05347 nptr[is], 0, scoef_vol+ib,
05348 option_data->sresp_filename[is]);
05349 else
05350 write_one_ts (option_data->sresp_filename[is],
05351 ts_length, scoef_vol+ib);
05352
05353 ib += ts_length;
05354 }
05355
05356
05357
05358 if (option_data->fitts_filename != NULL)
05359 if (nxyz > 1)
05360 write_ts_array (argc, argv, option_data, nt, 1, 0, fitts_vol,
05361 option_data->fitts_filename);
05362 else
05363 write_one_ts (option_data->fitts_filename, nt, fitts_vol);
05364
05365
05366
05367
05368 if (option_data->errts_filename != NULL)
05369 if (nxyz > 1)
05370 write_ts_array (argc, argv, option_data, nt, 1, 0, errts_vol,
05371 option_data->errts_filename);
05372 else
05373 write_one_ts (option_data->errts_filename, nt, errts_vol);
05374
05375 }
05376
05377
05378
05379 #if 0
05380 void terminate_program
05381 (
05382 DC_options ** option_data,
05383 float *** stimulus,
05384 matrix ** glt_cmat,
05385
05386 float *** coef_vol,
05387 float *** scoef_vol,
05388 float *** tcoef_vol,
05389 float *** fpart_vol,
05390 float *** rpart_vol,
05391 float ** mse_vol,
05392 float ** ffull_vol,
05393 float ** rfull_vol,
05394
05395 float **** glt_coef_vol,
05396 float **** glt_tcoef_vol,
05397 float *** glt_fstat_vol,
05398 float *** glt_rstat_vol,
05399
05400 float *** fitts_vol,
05401 float *** errts_vol
05402 )
05403
05404 {
05405 int p;
05406 int num_stimts;
05407 int ip;
05408 int is;
05409 int num_glt;
05410 int iglt;
05411 int * glt_rows;
05412 int ilc;
05413 int it;
05414 int nt;
05415
05416
05417
05418 if (*option_data == NULL) return;
05419 p = (*option_data)->p;
05420 num_stimts = (*option_data)->num_stimts;
05421 num_glt = (*option_data)->num_glt;
05422 glt_rows = (*option_data)->glt_rows;
05423 nt = (*option_data)->nt;
05424
05425
05426
05427 free (*option_data); *option_data = NULL;
05428
05429
05430
05431 if (*stimulus != NULL)
05432 {
05433 for (is = 0; is < num_stimts; is++)
05434 if ((*stimulus)[is] != NULL)
05435 { free ((*stimulus)[is]); (*stimulus)[is] = NULL; }
05436 free (*stimulus); *stimulus = NULL;
05437 }
05438
05439
05440
05441 if (*glt_cmat != NULL)
05442 {
05443 for (iglt = 0; iglt < num_glt; iglt++)
05444 matrix_destroy (&((*glt_cmat)[iglt]));
05445 free (*glt_cmat); *glt_cmat = NULL;
05446 }
05447
05448
05449
05450 if (*coef_vol != NULL)
05451 {
05452 for (ip = 0; ip < p; ip++)
05453 if ((*coef_vol)[ip] != NULL)
05454 { free ((*coef_vol)[ip]); (*coef_vol)[ip] = NULL; }
05455 free (*coef_vol); *coef_vol = NULL;
05456 }
05457
05458 if (*scoef_vol != NULL)
05459 {
05460 for (ip = 0; ip < p; ip++)
05461 if ((*scoef_vol)[ip] != NULL)
05462 { free ((*scoef_vol)[ip]); (*scoef_vol)[ip] = NULL; }
05463 free (*scoef_vol); *scoef_vol = NULL;
05464 }
05465
05466 if (*tcoef_vol != NULL)
05467 {
05468 for (ip = 0; ip < p; ip++)
05469 if ((*tcoef_vol)[ip] != NULL)
05470 { free ((*tcoef_vol)[ip]); (*tcoef_vol)[ip] = NULL; }
05471 free (*tcoef_vol); *tcoef_vol = NULL;
05472 }
05473
05474 if (*fpart_vol != NULL)
05475 {
05476 for (is = 0; is < num_stimts; is++)
05477 if ((*fpart_vol)[is] != NULL)
05478 { free ((*fpart_vol)[is]); (*fpart_vol)[is] = NULL; }
05479 free (*fpart_vol); *fpart_vol = NULL;
05480 }
05481
05482 if (*rpart_vol != NULL)
05483 {
05484 for (is = 0; is < num_stimts; is++)
05485 if ((*rpart_vol)[is] != NULL)
05486 { free ((*rpart_vol)[is]); (*rpart_vol)[is] = NULL; }
05487 free (*rpart_vol); *rpart_vol = NULL;
05488 }
05489
05490 if (*mse_vol != NULL) { free (*mse_vol); *mse_vol = NULL; }
05491 if (*ffull_vol != NULL) { free (*ffull_vol); *ffull_vol = NULL; }
05492 if (*rfull_vol != NULL) { free (*rfull_vol); *rfull_vol = NULL; }
05493
05494
05495
05496 if (*glt_coef_vol != NULL)
05497 {
05498 for (iglt = 0; iglt < num_glt; iglt++)
05499 if ((*glt_coef_vol)[iglt] != NULL)
05500 {
05501 for (ilc = 0; ilc < glt_rows[iglt]; ilc++)
05502 if ((*glt_coef_vol)[iglt][ilc] != NULL)
05503 {
05504 free ((*glt_coef_vol)[iglt][ilc]);
05505 (*glt_coef_vol)[iglt][ilc] = NULL;
05506 }
05507 free ((*glt_coef_vol)[iglt]); (*glt_coef_vol)[iglt] = NULL;
05508 }
05509 free (*glt_coef_vol); *glt_coef_vol = NULL;
05510 }
05511
05512 if (*glt_tcoef_vol != NULL)
05513 {
05514 for (iglt = 0; iglt < num_glt; iglt++)
05515 if ((*glt_tcoef_vol)[iglt] != NULL)
05516 {
05517 for (ilc = 0; ilc < glt_rows[iglt]; ilc++)
05518 if ((*glt_tcoef_vol)[iglt][ilc] != NULL)
05519 {
05520 free ((*glt_tcoef_vol)[iglt][ilc]);
05521 (*glt_tcoef_vol)[iglt][ilc] = NULL;
05522 }
05523 free ((*glt_tcoef_vol)[iglt]); (*glt_tcoef_vol)[iglt] = NULL;
05524 }
05525 free (*glt_tcoef_vol); *glt_tcoef_vol = NULL;
05526 }
05527
05528 if (*glt_fstat_vol != NULL)
05529 {
05530 for (iglt = 0; iglt < num_glt; iglt++)
05531 if ((*glt_fstat_vol)[iglt] != NULL)
05532 { free ((*glt_fstat_vol)[iglt]); (*glt_fstat_vol)[iglt] = NULL; }
05533 free (*glt_fstat_vol); *glt_fstat_vol = NULL;
05534 }
05535
05536 if (*glt_rstat_vol != NULL)
05537 {
05538 for (iglt = 0; iglt < num_glt; iglt++)
05539 if ((*glt_rstat_vol)[iglt] != NULL)
05540 { free ((*glt_rstat_vol)[iglt]); (*glt_rstat_vol)[iglt] = NULL; }
05541 free (*glt_rstat_vol); *glt_rstat_vol = NULL;
05542 }
05543
05544
05545
05546 if (*fitts_vol != NULL)
05547 {
05548 for (it = 0; it < nt; it++)
05549 { free ((*fitts_vol)[it]); (*fitts_vol)[it] = NULL; }
05550 free (*fitts_vol); *fitts_vol = NULL;
05551 }
05552
05553 if (*errts_vol != NULL)
05554 {
05555 for (it = 0; it < nt; it++)
05556 { free ((*errts_vol)[it]); (*errts_vol)[it] = NULL; }
05557 free (*errts_vol); *errts_vol = NULL;
05558 }
05559
05560 }
05561 #endif
05562
05563
05564
05565
05566
05567 int main
05568 (
05569 int argc,
05570 char ** argv
05571 )
05572
05573 {
05574 DC_options * option_data;
05575 THD_3dim_dataset * dset_time = NULL;
05576 byte * mask_vol = NULL;
05577 float * fmri_data = NULL;
05578 int fmri_length;
05579 float * censor_array = NULL;
05580 int censor_length;
05581 int * good_list = NULL;
05582 int * block_list = NULL;
05583 int num_blocks;
05584 float ** stimulus = NULL;
05585 int * stim_length = NULL;
05586 matrix * glt_cmat = NULL;
05587
05588 float ** coef_vol = NULL;
05589 float ** scoef_vol = NULL;
05590 float ** tcoef_vol = NULL;
05591 float ** fpart_vol = NULL;
05592 float ** rpart_vol = NULL;
05593
05594 float * mse_vol = NULL;
05595 float * ffull_vol = NULL;
05596 float * rfull_vol = NULL;
05597
05598 float *** glt_coef_vol = NULL;
05599 float *** glt_tcoef_vol = NULL;
05600 float ** glt_fstat_vol = NULL;
05601 float ** glt_rstat_vol = NULL;
05602
05603 float ** fitts_vol = NULL;
05604 float ** errts_vol = NULL;
05605
05606
05607 #ifdef FLOATIZE
05608
05609 if( argc < 2 || strcmp(argv[1],"-OK") != 0 ){
05610 fprintf(stderr,"**\n"
05611 "** 3dDeconvolve_f is now disabled by default.\n"
05612 "** It is dangerous, due to roundoff problems.\n"
05613 "** Please use 3dDeconvolve from now on!\n"
05614 "**\n"
05615 "** HOWEVER, if you insist on using 3dDeconvolve_f, then:\n"
05616 "** + Use '-OK' as the first command line option.\n"
05617 "** + Check the matrix condition number;\n"
05618 "** if it is greater than 100, BEWARE!\n"
05619 "**\n"
05620 "** RWCox - July 2004\n"
05621 "**\n"
05622 ) ;
05623 exit(0) ;
05624 }
05625 #endif
05626
05627 #ifdef USING_MCW_MALLOC
05628 enable_mcw_malloc() ;
05629 #endif
05630 mainENTRY("3dDeconvolve main") ;
05631
05632
05633 (void) COX_clock_time() ;
05634
05635
05636 initialize_program (argc, argv, &option_data, &dset_time, &mask_vol,
05637 &fmri_data, &fmri_length, &censor_array, &censor_length, &good_list,
05638 &block_list, &num_blocks, &stimulus, &stim_length, &glt_cmat,
05639 &coef_vol, &scoef_vol, &tcoef_vol, &fpart_vol, &rpart_vol,
05640 &mse_vol, &ffull_vol, &rfull_vol, &glt_coef_vol, &glt_tcoef_vol,
05641 &glt_fstat_vol, &glt_rstat_vol, &fitts_vol, &errts_vol);
05642
05643 if( xrestore ){
05644 do_xrestore_stuff( argc,argv , option_data ) ;
05645 exit(0) ;
05646 }
05647
05648
05649
05650 calculate_results (option_data, dset_time, mask_vol, fmri_data, fmri_length,
05651 good_list, block_list, num_blocks, stimulus, stim_length,
05652 glt_cmat, coef_vol, scoef_vol, tcoef_vol, fpart_vol,
05653 rpart_vol, mse_vol, ffull_vol, rfull_vol, glt_coef_vol, glt_tcoef_vol,
05654 glt_fstat_vol, glt_rstat_vol, fitts_vol, errts_vol);
05655
05656
05657
05658 if (dset_time != NULL)
05659 { THD_delete_3dim_dataset (dset_time, False); dset_time = NULL; }
05660 if (mask_vol != NULL)
05661 { free (mask_vol); mask_vol = NULL; }
05662
05663
05664
05665 if (!option_data->nodata)
05666 output_results (argc, argv, option_data, coef_vol, scoef_vol, tcoef_vol,
05667 fpart_vol, rpart_vol, mse_vol, ffull_vol, rfull_vol,
05668 glt_coef_vol, glt_tcoef_vol, glt_fstat_vol, glt_rstat_vol,
05669 fitts_vol, errts_vol);
05670
05671
05672 #if 0
05673
05674 terminate_program (&option_data, &stimulus, &glt_cmat, &coef_vol, &scoef_vol,
05675 &tcoef_vol, &fpart_vol, &rpart_vol, & mse_vol, &ffull_vol,
05676 &rfull_vol, &glt_coef_vol, &glt_tcoef_vol, &glt_fstat_vol,
05677 &glt_rstat_vol, &fitts_vol, &errts_vol);
05678 #endif
05679
05680 if( proc_use_jobs == 1 && verb ){
05681 fprintf(stderr,"++ Program finished; elapsed time=%.3f\n",COX_clock_time());
05682 }
05683 #ifndef FLOATIZE
05684 if( proc_numjob == 1 && verb ){
05685 double fv=get_matrix_flops() , fd=get_matrix_dotlen() ;
05686 if( fv > 0.0 && fd > 0.0 )
05687 fprintf(stderr,"++ #Flops=%g Average Dot Product=%g\n",fv,fd) ;
05688 }
05689 #endif
05690
05691 exit(0);
05692 }
05693
05694
05695
05696
05697 #include "coxplot.h"
05698
05699 #define TSGRAY_SEPARATE_YSCALE (1<<0)
05700 #define TSGRAY_FLIP_XY (1<<1)
05701
05702
05703
05704
05705
05706
05707
05708
05709
05710
05711
05712
05713 MEM_plotdata * PLOT_tsgray( int npt , int nts , int ymask , float **y )
05714 {
05715 MEM_plotdata *mp ;
05716 float ybot,ytop , yfac , dx,dy , val ;
05717 int ii,jj , flipxy ;
05718 char str[32] ;
05719 int sepscl ;
05720 float boxr=1.0,boxg=0.9,boxb=0.0 ;
05721 char *eee ;
05722 int dobox=1 ;
05723
05724 if( npt < 2 || nts < 1 || y == NULL ) return NULL ;
05725
05726 eee = my_getenv( "AFNI_XJPEG_COLOR" ) ;
05727 if( eee != NULL ){
05728 float rf=-1.0, gf=-1.0 , bf=-1.0 ;
05729 sscanf( eee , "rgbi:%f/%f/%f" , &rf,&gf,&bf ) ;
05730 if( rf >= 0.0 && rf <= 1.0 && gf >= 0.0 && gf <= 1.0 && bf >= 0.0 && bf <= 1.0 ){
05731 boxr = rf ; boxg = gf ; boxb = bf ;
05732 }
05733 if( NOISH(eee) ) dobox = 0 ;
05734 }
05735
05736
05737
05738 ybot = ytop = y[0][0] ;
05739 for( jj=0 ; jj < nts ; jj++ ){
05740 for( ii=0 ; ii < npt ; ii++ ){
05741 val = y[jj][ii] ;
05742 if( ybot > val ) ybot = val ;
05743 else if( ytop < val ) ytop = val ;
05744 }
05745 }
05746 if( ybot >= ytop ) return NULL ;
05747 yfac = 1.0/(ytop-ybot) ;
05748 dx = 0.9999/npt ;
05749 dy = 0.9999/nts ;
05750
05751 create_memplot_surely( "Gplot" , 1.0 ) ;
05752 set_color_memplot( 0.0 , 0.0 , 0.0 ) ;
05753 set_thick_memplot( 0.0 ) ;
05754 set_opacity_memplot( 1.0 ) ;
05755
05756 flipxy = (ymask & TSGRAY_FLIP_XY) != 0 ;
05757 sepscl = (ymask & TSGRAY_SEPARATE_YSCALE) != 0 ;
05758
05759 for( jj=0 ; jj < nts ; jj++ ){
05760
05761 if( sepscl ){
05762 ybot = ytop = y[jj][0] ;
05763 for( ii=1 ; ii < npt ; ii++ ){
05764 val = y[jj][ii] ;
05765 if( ybot > val ) ybot = val ;
05766 else if( ytop < val ) ytop = val ;
05767 }
05768 if( ybot >= ytop ) yfac = 1.0 ;
05769 else yfac = 1.0/(ytop-ybot) ;
05770 }
05771
05772 for( ii=0 ; ii < npt ; ii++ ){
05773 val = yfac*(ytop-y[jj][ii]) ;
05774 set_color_memplot( val,val,val ) ;
05775 if( flipxy )
05776 plotrect_memplot( ii*dx,jj*dy , (ii+1)*dx,(jj+1)*dy ) ;
05777 else
05778 plotrect_memplot( jj*dy,1.0-ii*dx , (jj+1)*dy,1.0-(ii+1)*dy ) ;
05779 }
05780 }
05781
05782 if( dobox ){
05783 set_color_memplot( boxr, boxg, boxb ) ;
05784 set_opacity_memplot( 0.789 ) ;
05785 for( jj=0 ; jj <= nts ; jj++ ){
05786 if( flipxy ){
05787 plotline_memplot( 1.0,jj*dy , 0.0,jj*dy ) ;
05788 } else {
05789 plotline_memplot( jj*dy,1.0 , jj*dy,0.0 ) ;
05790 }
05791 }
05792 }
05793
05794 set_opacity_memplot( 1.0 ) ;
05795 set_color_memplot( 0.0 , 0.0 , 0.0 ) ;
05796 mp = get_active_memplot() ;
05797 return mp ;
05798 }
05799
05800
05801
05802 MRI_IMAGE * PLOT_matrix_gray( matrix X )
05803 {
05804 int nts=X.cols , npt=X.rows , ii,jj , nxim=768 , nyim=1024 ;
05805 MEM_plotdata *mp ;
05806 float **xar ;
05807 MRI_IMAGE *im ;
05808 char *eee ;
05809
05810 if( nts < 1 || npt < 2 ) return NULL ;
05811
05812 xar = (float **)malloc( sizeof(float *)*nts ) ;
05813 for( jj=0 ; jj < nts ; jj++ ){
05814 xar[jj] = (float *)malloc( sizeof(float *)*npt ) ;
05815 for( ii=0 ; ii < npt ; ii++ ) xar[jj][ii] = X.elts[ii][jj] ;
05816 }
05817
05818 mp = PLOT_tsgray( npt , nts , TSGRAY_SEPARATE_YSCALE , xar ) ;
05819
05820 for( jj=0 ; jj < nts ; jj++ ) free((void *)xar[jj]) ;
05821 free((void *)xar) ;
05822
05823 if( mp == NULL ) return NULL ;
05824
05825 eee = my_getenv( "AFNI_XJPEG_IMXY") ;
05826 if( eee != NULL ){
05827 int a=-1, b=-1 ;
05828 sscanf( eee , "%dx%d" , &a,&b ) ;
05829 if( a > 99 && b > 99 ){ nxim = a ; nyim = b ; }
05830 if( a < -1 ) nxim = -a*nts ;
05831 if( b < 0 ) nyim = -b*npt ;
05832 }
05833
05834 im = mri_new( nxim , nyim , MRI_rgb ) ;
05835 memset( MRI_RGB_PTR(im) , 255 , 3*im->nvox ) ;
05836 memplot_to_RGB_sef( im , mp , 0,0,1 ) ;
05837 delete_memplot( mp ) ;
05838 return im ;
05839 }
05840
05841
05842
05843 void JPEG_matrix_gray( matrix X , char *fname )
05844 {
05845 char *pg , *jpfilt ;
05846 MRI_IMAGE *im ;
05847 FILE *fp ;
05848
05849 if( fname == NULL || *fname == '\0' ) return ;
05850
05851 pg = THD_find_executable( "cjpeg" ) ;
05852 if( pg == NULL ){
05853 fprintf(stderr,
05854 "** WARNING: can't save %s because program 'cjpeg' not in path!\n",
05855 fname) ;
05856 return ;
05857 }
05858
05859 im = PLOT_matrix_gray( X ) ;
05860 if( im == NULL ){
05861 fprintf(stderr,
05862 "** WARNING: can't save %s because of internal error!\n",fname) ;
05863 return ;
05864 }
05865
05866 jpfilt = (char *)malloc( sizeof(char)*(strlen(pg)+strlen(fname)+32) ) ;
05867 sprintf( jpfilt , "%s -quality 95 > %s" , pg , fname ) ;
05868 #ifndef CYGWIN
05869 signal( SIGPIPE , SIG_IGN ) ; errno = 0 ;
05870 #endif
05871 fp = popen( jpfilt , "w" ) ;
05872 if( fp == NULL ){
05873 mri_free(im) ; free((void *)jpfilt) ;
05874 fprintf(stderr,"** WARNING: can't save %s because cjpeg filter fails!\n",fname) ;
05875 return ;
05876 }
05877 fprintf(fp,"P6\n%d %d\n255\n" , im->nx,im->ny ) ;
05878 fwrite( MRI_RGB_PTR(im), sizeof(byte), 3*im->nvox, fp ) ;
05879 (void) pclose(fp) ;
05880 if( verb ) fprintf(stderr,"++ Wrote matrix image to file %s\n",fname) ;
05881
05882 mri_free(im) ; free((void *)jpfilt) ; return ;
05883 }
05884
05885
05886
05887
05888
05889
05890
05891
05892
05893 #include "niml.h"
05894
05895 #if 0
05896
05897
05898
05899
05900 char * niml_element_to_string( NI_element *nel )
05901 {
05902 NI_stream ns ;
05903 char *stout ;
05904 int ii,jj ;
05905
05906 if( nel == NULL ) return NULL ;
05907
05908 ns = NI_stream_open( "str:" , "w" ) ;
05909 (void) NI_write_element( ns , nel , NI_TEXT_MODE ) ;
05910 stout = strdup( NI_stream_getbuf(ns) ) ;
05911 NI_stream_close( ns ) ;
05912 jj = strlen(stout) ;
05913 for( ii=jj-1 ; ii > 0 && isspace(stout[ii]) ; ii-- ) ;
05914 stout[ii+1] = '\0' ;
05915 return stout ;
05916 }
05917 #endif
05918
05919
05920
05921
05922
05923
05924
05925 NI_element * matrix_to_niml( matrix a , char *ename )
05926 {
05927 int m=a.rows , n=a.cols , ii,jj ;
05928 MTYPE **aar = a.elts ;
05929 NI_element *nel ;
05930 double *ecol ;
05931
05932 if( m < 1 || n < 1 || aar == NULL ) return NULL ;
05933
05934 if( ename == NULL || *ename == '\0' ) ename = "matrix" ;
05935
05936 nel = NI_new_data_element( ename , m ) ;
05937 ecol = (double *) malloc(sizeof(double)*m) ;
05938
05939 for( jj=0 ; jj < n ; jj++ ){
05940 for( ii=0 ; ii < m ; ii++ ) ecol[ii] = aar[ii][jj] ;
05941 NI_add_column( nel , NI_DOUBLE , ecol ) ;
05942 }
05943
05944 free((void *)ecol) ;
05945 return nel ;
05946 }
05947
05948
05949
05950
05951
05952 void niml_to_matrix( NI_element *nel , matrix *a )
05953 {
05954 int m , n , ii , jj ;
05955 double *ecol ;
05956 MTYPE **aar ;
05957 char *ename ;
05958
05959 if( nel == NULL || a == NULL ) return ;
05960
05961 m = nel->vec_len ;
05962 n = nel->vec_num ;
05963
05964 if( m < 1 || n < 1 ) return ;
05965 for( jj=0 ; jj < n ; jj++ )
05966 if( nel->vec_typ[jj] != NI_DOUBLE ) return ;
05967
05968 matrix_create( m , n , a ) ;
05969 aar = a->elts ;
05970
05971 for( jj=0 ; jj < n ; jj++ ){
05972 ecol = (double *)nel->vec[jj] ;
05973 for( ii=0 ; ii < m ; ii++ ) aar[ii][jj] = (MTYPE)ecol[ii] ;;
05974 }
05975
05976 return ;
05977 }
05978
05979
05980
05981
05982
05983
05984
05985 NI_element * intvec_to_niml( int nvec , int *vec , char *ename )
05986 {
05987 NI_element *nel ;
05988
05989 if( nvec < 1 || vec == NULL ) return NULL ;
05990
05991 if( ename == NULL || *ename == '\0' ) ename = "intvec" ;
05992
05993 nel = NI_new_data_element( ename , nvec ) ;
05994 NI_add_column( nel , NI_INT , vec ) ;
05995 return nel ;
05996 }
05997
05998
05999
06000 void niml_to_intvec( NI_element *nel , int *nvec , int **vec )
06001 {
06002 if( nel == NULL || nvec == NULL || vec == NULL ) return ;
06003
06004 if( nel->vec_len < 1 || nel->vec_num < 1 ) return ;
06005 if( nel->vec_typ[0] != NI_INT ) return ;
06006
06007 *nvec = nel->vec_len ;
06008 *vec = (int *)malloc(*nvec*sizeof(int)) ;
06009 memcpy(*vec,nel->vec[0],*nvec*sizeof(int)) ;
06010 return ;
06011 }
06012
06013
06014
06015 NI_element * stringvec_to_niml( int nstr , char **str , char *ename )
06016 {
06017 NI_element *nel ;
06018
06019 if( nstr < 1 || str == NULL ) return NULL ;
06020
06021 if( ename == NULL || *ename == '\0' ) ename = "stringvec" ;
06022
06023 nel = NI_new_data_element( ename , nstr ) ;
06024 NI_add_column( nel , NI_STRING , str ) ;
06025 return nel ;
06026 }
06027
06028
06029
06030 void niml_to_stringvec( NI_element *nel , int *nstr , char ***str )
06031 {
06032 int ii ;
06033 char **sv ;
06034
06035 if( nel == NULL || nstr == NULL || str == NULL ) return ;
06036
06037 if( nel->vec_len < 1 || nel->vec_num < 1 ) return ;
06038 if( nel->vec_typ[0] != NI_STRING ) return ;
06039 sv = (char **) nel->vec[0] ;
06040
06041 *nstr = nel->vec_len ;
06042 *str = (char **)calloc(sizeof(char *),*nstr) ;
06043 for( ii=0 ; ii < *nstr ; ii++ ) (*str)[ii] = strdup(sv[ii]) ;
06044 return ;
06045 }
06046
06047
06048
06049 NI_element * symvec_to_niml( int ns , SYM_irange *sv , char *ename )
06050 {
06051 NI_element *nel ;
06052 int ii , *bc,*tc,*gc ; char **sc ;
06053
06054 if( ns < 1 || sv == NULL ) return NULL ;
06055
06056 if( ename == NULL || *ename == '\0' ) ename = "symvec" ;
06057
06058 nel = NI_new_data_element( ename , ns ) ;
06059
06060 bc = (int *) malloc(sizeof(int) *ns) ;
06061 tc = (int *) malloc(sizeof(int) *ns) ;
06062 gc = (int *) malloc(sizeof(int) *ns) ;
06063 sc = (char **)malloc(sizeof(char *)*ns) ;
06064
06065 for( ii=0 ; ii < ns ; ii++ ){
06066 bc[ii] = sv[ii].nbot ; tc[ii] = sv[ii].ntop ;
06067 gc[ii] = sv[ii].gbot ; sc[ii] = sv[ii].name ;
06068 }
06069
06070 NI_add_column( nel , NI_INT , bc ); free((void *)bc) ;
06071 NI_add_column( nel , NI_INT , tc ); free((void *)tc) ;
06072 NI_add_column( nel , NI_INT , gc ); free((void *)gc) ;
06073 NI_add_column( nel , NI_STRING , sc ); free((void *)sc) ;
06074
06075 return nel ;
06076 }
06077
06078
06079
06080 void niml_to_symvec( NI_element *nel , int *ns , SYM_irange **sv )
06081 {
06082 int ii , *bc,*tc,*gc ; char **sc ;
06083
06084 if( nel == NULL || ns == NULL || sv == NULL ) return ;
06085
06086 if( nel->vec_len < 1 || nel->vec_num < 4 ) return ;
06087 if( nel->vec_typ[0] != NI_INT ) return ; bc = (int *) nel->vec[0] ;
06088 if( nel->vec_typ[1] != NI_INT ) return ; tc = (int *) nel->vec[1] ;
06089 if( nel->vec_typ[2] != NI_INT ) return ; gc = (int *) nel->vec[2] ;
06090 if( nel->vec_typ[3] != NI_STRING ) return ; sc = (char **)nel->vec[3] ;
06091
06092 *ns = nel->vec_len ;
06093 *sv = (SYM_irange *)calloc(sizeof(SYM_irange),*ns) ;
06094 for( ii=0 ; ii < *ns ; ii++ ){
06095 (*sv)[ii].nbot = bc[ii] ;
06096 (*sv)[ii].ntop = tc[ii] ;
06097 (*sv)[ii].gbot = gc[ii] ;
06098 NI_strncpy( (*sv)[ii].name , sc[ii] , 64 ) ;
06099 }
06100 return ;
06101 }
06102
06103
06104
06105 void XSAVE_output( char *prefix )
06106 {
06107 char *fname , *cpt ;
06108 NI_stream ns ;
06109 NI_element *nel ;
06110 int nimode = NI_BINARY_MODE ;
06111
06112 if( !xsave ) return ;
06113
06114 if( AFNI_yesenv("AFNI_XSAVE_TEXT") ) nimode = NI_TEXT_MODE ;
06115
06116
06117
06118 if( prefix == NULL || *prefix == '\0' ) prefix = "X" ;
06119
06120 fname = malloc( sizeof(char) * (strlen(prefix)+32) ) ;
06121 strcpy(fname,"file:") ; strcat(fname,prefix) ; strcat(fname,".xsave") ;
06122 if( THD_is_ondisk(fname+5) && verb )
06123 fprintf(stderr,
06124 "** WARNING: -xsave output file %s will be overwritten!\n",fname+5) ;
06125 ns = NI_stream_open( fname , "w" ) ;
06126 if( ns == (NI_stream)NULL ){
06127 fprintf(stderr,
06128 "** ERROR: -xsave output file %s can't be opened!\n",fname+5) ;
06129 free((void *)fname) ;
06130 return ;
06131 }
06132
06133
06134
06135 nel = NI_new_data_element( "XSAVE", 0 ) ;
06136 if( InputFilename != NULL )
06137 NI_set_attribute( nel , "InputFilename" , InputFilename ) ;
06138 if( BucketFilename != NULL )
06139 NI_set_attribute( nel , "BucketFilename" , BucketFilename ) ;
06140 if( CoefFilename != NULL )
06141 NI_set_attribute( nel , "CoefFilename" , CoefFilename ) ;
06142
06143 sprintf(fname,"%d",X.rows) ;
06144 NI_set_attribute( nel , "NumTimePoints" , fname ) ;
06145 sprintf(fname,"%d",X.cols) ;
06146 NI_set_attribute( nel , "NumRegressors" , fname ) ;
06147
06148 NI_set_attribute( nel , "Deconvolveries" , XSAVE_version ) ;
06149
06150 cpt = tross_datetime() ;
06151 NI_set_attribute( nel , "DateTime" , cpt ) ;
06152 free((void *)cpt) ;
06153 cpt = tross_hostname() ;
06154 NI_set_attribute( nel , "Hostname" , cpt ) ;
06155 free((void *)cpt) ;
06156 cpt = tross_username() ;
06157 NI_set_attribute( nel , "Username" , cpt ) ;
06158 free((void *)cpt) ;
06159
06160 (void) NI_write_element( ns , nel , nimode ) ;
06161 NI_free_element( nel ) ;
06162
06163
06164
06165 nel = matrix_to_niml( X , "matrix" ) ;
06166 NI_set_attribute( nel , "Xname" , "X" ) ;
06167 (void) NI_write_element( ns , nel , nimode ) ;
06168 NI_free_element( nel ) ;
06169
06170 nel = matrix_to_niml( XtXinv , "matrix" ) ;
06171 NI_set_attribute( nel , "Xname" , "XtXinv" ) ;
06172 (void) NI_write_element( ns , nel , nimode ) ;
06173 NI_free_element( nel ) ;
06174
06175 nel = matrix_to_niml( XtXinvXt , "matrix" ) ;
06176 NI_set_attribute( nel , "Xname" , "XtXinvXt" ) ;
06177 (void) NI_write_element( ns , nel , nimode ) ;
06178 NI_free_element( nel ) ;
06179
06180
06181
06182 nel = intvec_to_niml( nGoodList , GoodList , "intvec" ) ;
06183 NI_set_attribute( nel , "Xname" , "GoodList" ) ;
06184 (void) NI_write_element( ns , nel , nimode ) ;
06185 NI_free_element( nel ) ;
06186
06187
06188
06189 if( ParamIndex != NULL ){
06190 nel = intvec_to_niml( nParam, ParamIndex , "intvec" ) ;
06191 NI_set_attribute( nel , "Xname" , "ParamIndex" ) ;
06192 (void) NI_write_element( ns , nel , nimode ) ;
06193 NI_free_element( nel ) ;
06194 }
06195
06196
06197
06198 nel = intvec_to_niml( nParam, ParamStim , "intvec" ) ;
06199 NI_set_attribute( nel , "Xname" , "ParamStim" ) ;
06200 (void) NI_write_element( ns , nel , nimode ) ;
06201 NI_free_element( nel ) ;
06202
06203
06204
06205 nel = stringvec_to_niml( nParam , ParamLabel , "stringvec" ) ;
06206 NI_set_attribute( nel , "Xname" , "ParamLabel" ) ;
06207 (void) NI_write_element( ns , nel , nimode ) ;
06208 NI_free_element( nel ) ;
06209
06210
06211
06212 if( nSymStim > 0 ){
06213 nel = symvec_to_niml( nSymStim , SymStim , "symvec" ) ;
06214 NI_set_attribute( nel , "Xname" , "SymStim" ) ;
06215 (void) NI_write_element( ns , nel , nimode ) ;
06216 NI_free_element( nel ) ;
06217 }
06218
06219
06220
06221 NI_stream_close(ns) ; free((void *)fname) ; return ;
06222 }
06223
06224
06225
06226 #define DPR(s) fprintf(stderr,"%s\n",(s))
06227
06228 void XSAVE_input( char *xname )
06229 {
06230 char *fname , *cpt ;
06231 NI_stream ns ;
06232 NI_element *nel ;
06233
06234 if( xname == NULL || *xname == '\0' ) return ;
06235
06236
06237
06238 fname = malloc( sizeof(char) * (strlen(xname)+32) ) ;
06239 strcpy(fname,"file:") ; strcat(fname,xname) ;
06240 ns = NI_stream_open( fname , "r" ) ;
06241 free((void *)fname) ;
06242 if( ns == (NI_stream)NULL ){
06243 fprintf(stderr,"** ERROR: can't open file %s for -xrestore!\n",xname) ;
06244 return ;
06245 }
06246
06247
06248
06249 nel = NI_read_element( ns , 1 ) ;
06250 if( nel == NULL ){
06251 fprintf(stderr,"** ERROR: can't read header in file %s for -xrestore!\n",xname) ;
06252 NI_stream_close( ns ) ; return ;
06253 }
06254
06255
06256
06257 cpt = NI_get_attribute( nel , "InputFilename" ) ;
06258 if( cpt != NULL ) InputFilename = strdup(cpt) ;
06259 cpt = NI_get_attribute( nel , "BucketFilename" ) ;
06260 if( cpt != NULL ) BucketFilename = strdup(cpt) ;
06261 cpt = NI_get_attribute( nel , "CoefFilename" ) ;
06262 if( cpt != NULL ) CoefFilename = strdup(cpt) ;
06263
06264
06265
06266 cpt = NI_get_attribute( nel , "NumTimePoints" ) ;
06267 if( cpt != NULL ) NumTimePoints = strtol(cpt,NULL,10) ;
06268 cpt = NI_get_attribute( nel , "NumRegressors" ) ;
06269 if( cpt != NULL ) NumRegressors = strtol(cpt,NULL,10) ;
06270
06271 NI_free_element( nel ) ;
06272
06273
06274
06275 while(1){
06276
06277 nel = NI_read_element( ns , 1 ) ;
06278 if( nel == NULL ) break ;
06279
06280 cpt = NI_get_attribute( nel , "Xname" ) ;
06281 if( cpt == NULL ){
06282 NI_free_element( nel ) ; continue ;
06283 }
06284
06285 if( strcmp(nel->name,"matrix") == 0 ){
06286
06287 if( strcmp(cpt,"X" )==0 ) niml_to_matrix( nel , &X );
06288 else if( strcmp(cpt,"XtXinv" )==0 ) niml_to_matrix( nel , &XtXinv );
06289 else if( strcmp(cpt,"XtXinvXt")==0 ) niml_to_matrix( nel , &XtXinvXt );
06290
06291 } else if( strcmp(nel->name,"intvec") == 0 ){
06292
06293 if( strcmp(cpt,"GoodList") == 0 )
06294 niml_to_intvec( nel , &nGoodList, &GoodList );
06295 else if( strcmp(cpt,"ParamIndex") == 0 )
06296 niml_to_intvec( nel , &nParam , &ParamIndex );
06297 else if( strcmp(cpt,"ParamStim" ) == 0 )
06298 niml_to_intvec( nel , &nParam , &ParamStim );
06299
06300 } else if( strcmp(nel->name,"stringvec") == 0 ){
06301
06302 if( strcmp(cpt,"ParamLabel") == 0 )
06303 niml_to_stringvec( nel , &nParam , &ParamLabel );
06304
06305 } else if( strcmp(nel->name,"symvec") == 0 ){
06306
06307 if( strcmp(cpt,"SymStim") == 0 )
06308 niml_to_symvec( nel , &nSymStim , &SymStim );
06309
06310 } else {
06311
06312 }
06313
06314 NI_free_element( nel ) ;
06315
06316 }
06317
06318 NI_stream_close( ns ) ; return ;
06319 }
06320
06321
06322
06323
06324
06325 void check_xrestore_data(void)
06326 {
06327 int nerr = 0 ;
06328
06329 if( X.rows < 1 || X.cols < 1 ){
06330 fprintf(stderr,
06331 "** ERROR: -xrestore %s has bad X matrix\n",xrestore_filename) ;
06332 nerr++ ;
06333 }
06334 if( XtXinv.rows < 1 || XtXinv.cols < 1 ){
06335 fprintf(stderr,
06336 "** ERROR: -xrestore %s has bad XtXinv matrix\n",xrestore_filename) ;
06337 nerr++ ;
06338 }
06339 if( XtXinvXt.rows < 1 || XtXinvXt.cols < 1 ){
06340 fprintf(stderr,
06341 "** ERROR: -xrestore %s has bad XtXinvXt matrix\n",xrestore_filename) ;
06342 nerr++ ;
06343 }
06344 if( nParam > 0 && X.cols != nParam ){
06345 fprintf(stderr,
06346 "** ERROR: -xrestore %s X matrix cols mismatch: %d != %d\n",
06347 xrestore_filename, X.cols , nParam ) ;
06348 nerr++ ;
06349 }
06350 if( GoodList == NULL ){
06351 fprintf(stderr,
06352 "** ERROR: -xrestore %s missing GoodList field\n",xrestore_filename) ;
06353 nerr++ ;
06354 } else if( nGoodList != X.rows ){
06355 fprintf(stderr,
06356 "** ERROR: -xrestore %s X matrix rows mismatch: %d != %d\n",
06357 xrestore_filename, X.cols , nGoodList ) ;
06358 nerr++ ;
06359 }
06360 #if 0
06361 if( ParamStim == NULL ){
06362 fprintf(stderr,
06363 "** ERROR: -xrestore %s missing ParamStim field\n",xrestore_filename) ;
06364 nerr++ ;
06365 }
06366 if( ParamLabel == NULL ){
06367 fprintf(stderr,
06368 "** ERROR: -xrestore %s missing ParamLabel field\n",xrestore_filename) ;
06369 nerr++ ;
06370 }
06371 #endif
06372 if( InputFilename == NULL ){
06373 fprintf(stderr,
06374 "** ERROR: -xrestore %s missing InputFilename field\n",xrestore_filename) ;
06375 nerr++ ;
06376 }
06377
06378 if( nerr > 0 ) exit(1) ;
06379 return ;
06380 }
06381
06382
06383
06384 void do_xrestore_stuff( int argc , char **argv , DC_options *option_data )
06385 {
06386 THD_3dim_dataset *dset_time , *dset_coef, *dset_buck ;
06387 int nt , np , ii, nxyz, tout,rout,fout, ixyz,novar,view ;
06388 char *buck_prefix , *buck_name ;
06389 char brick_label[THD_MAX_NAME] ;
06390
06391 float *ts_array = NULL;
06392 #undef NGET
06393 #define NGET 32
06394 int nget=0 ,
06395 cget=0 ,
06396 jget ,
06397 iget[NGET] ;
06398 MRI_IMARR *imget=NULL ;
06399
06400 int num_glt , iglt , ilc,nlc ;
06401 matrix *glt_cmat , *glt_amat , *cxtxinvct ;
06402 vector *glt_coef , *glt_tcoef, y , coef ;
06403 float *fglt , *rglt ;
06404 char **glt_label ;
06405 int *glt_rows ;
06406 float ***glt_coef_vol ,
06407 ***glt_tcoef_vol=NULL ,
06408 ** glt_fstat_vol=NULL ,
06409 ** glt_rstat_vol=NULL ;
06410 float *cdar=NULL , ssef , *volume ;
06411 int ivol , nvol , nbuck , vstep ;
06412
06413
06414
06415 num_glt = option_data->num_glt ;
06416 glt_label = option_data->glt_label ;
06417 glt_rows = option_data->glt_rows ;
06418
06419 tout = option_data->tout ;
06420 rout = option_data->rout ;
06421 fout = option_data->fout ;
06422
06423 if( num_glt < 1 ){
06424 fprintf(stderr,"** ERROR: -xrestore with no new GLTs???\n"); exit(1);
06425 }
06426
06427
06428
06429 matrix_initialize( &X ) ;
06430 matrix_initialize( &XtXinv ) ;
06431 matrix_initialize( &XtXinvXt ) ;
06432 nGoodList = 0 ; GoodList = NULL ;
06433 nParam = 0 ; ParamIndex = NULL ; ParamStim = NULL ; ParamLabel = NULL ;
06434
06435
06436
06437 if( verb ) fprintf(stderr,"++ Starting -xrestore %s\n",xrestore_filename) ;
06438
06439 XSAVE_input( xrestore_filename ) ;
06440 check_xrestore_data() ;
06441
06442 nt = X.rows ; np = X.cols ;
06443
06444
06445
06446 if( verb ) fprintf(stderr,"++ loading time series dataset '%s'\n",InputFilename);
06447 dset_time = THD_open_dataset( InputFilename ) ;
06448 if( dset_time == NULL ){
06449 fprintf(stderr,
06450 "** ERROR: -xrestore can't open time series dataset '%s'\n" ,
06451 InputFilename ) ;
06452 exit(1) ;
06453 }
06454 DSET_load( dset_time ) ;
06455 if( !DSET_LOADED(dset_time) ){
06456 fprintf(stderr,
06457 "** ERROR: -xrestore can't load time series dataset '%s'\n" ,
06458 InputFilename ) ;
06459 exit(1) ;
06460 }
06461 nxyz = DSET_NVOX(dset_time) ;
06462 view = dset_time->view_type ;
06463
06464
06465
06466 dset_coef = NULL ;
06467 if( CoefFilename != NULL ){
06468 if( verb) fprintf(stderr,"++ loading coefficient dataset %s\n",CoefFilename);
06469 dset_coef = THD_open_dataset( CoefFilename ) ;
06470 if( dset_coef == NULL ){
06471 fprintf(stderr,
06472 "** WARNING: -xrestore can't open coefficient dataset %s\n",
06473 CoefFilename);
06474 } else {
06475 DSET_load(dset_coef) ;
06476 if( !DSET_LOADED(dset_coef) ){
06477 fprintf(stderr,
06478 "** WARNING: -xrestore can't load coefficient dataset %s\n",
06479 CoefFilename);
06480 DSET_delete(dset_coef) ; dset_coef = NULL ;
06481 }
06482 }
06483 if( dset_coef != NULL && DSET_NVALS(dset_coef) < np ){
06484 fprintf(stderr,
06485 "** WARNING: -xrestore coefficient dataset %s too short\n",
06486 CoefFilename);
06487 DSET_delete(dset_coef) ; dset_coef = NULL ;
06488 }
06489 if( dset_coef != NULL ){
06490 if( ParamIndex != NULL ) free((void *)ParamIndex) ;
06491 ParamIndex = (int *)malloc(sizeof(int)*np) ;
06492 for( ii=0 ; ii < np ; ii++ ) ParamIndex[ii] = ii ;
06493 }
06494 }
06495
06496
06497
06498 if( dset_coef == NULL && BucketFilename != NULL && ParamIndex != NULL ){
06499 if( verb ) fprintf(stderr,"++ loading original bucket dataset %s\n",BucketFilename) ;
06500 dset_coef = THD_open_dataset( BucketFilename ) ;
06501 if( dset_coef == NULL ){
06502 fprintf(stderr,
06503 "** WARNING: -xrestore can't open old bucket dataset %s\n",
06504 BucketFilename);
06505 } else {
06506 DSET_load(dset_coef) ;
06507 if( !DSET_LOADED(dset_coef) ){
06508 fprintf(stderr,
06509 "** WARNING: -xrestore can't load old bucket dataset %s\n",
06510 BucketFilename);
06511 DSET_delete(dset_coef) ; dset_coef = NULL ;
06512 }
06513 }
06514 if( dset_coef != NULL && DSET_NVALS(dset_coef) < ParamIndex[np-1] ){
06515 fprintf(stderr,
06516 "** WARNING: -xrestore old bucket dataset %s too short\n",
06517 BucketFilename);
06518 DSET_delete(dset_coef) ; dset_coef = NULL ;
06519 }
06520 }
06521
06522 if( ISVALID_DSET(dset_coef) && DSET_NVOX(dset_coef) != nxyz ){
06523 fprintf(stderr,
06524 "** ERROR: dataset mismatch between time series and coef!\n") ;
06525 exit(1) ;
06526 }
06527
06528
06529
06530 if( dset_coef == NULL && verb )
06531 fprintf(stderr,
06532 "++ -xrestore recomputing coefficients from time series\n");
06533
06534
06535
06536 glt_cmat = (matrix *) malloc( sizeof(matrix) * num_glt ) ;
06537
06538 for( iglt=0; iglt < num_glt ; iglt++){
06539 matrix_initialize( glt_cmat + iglt ) ;
06540 #if 1
06541 read_glt_matrix( option_data->glt_filename[iglt] ,
06542 option_data->glt_rows + iglt ,
06543 np , glt_cmat + iglt ) ;
06544 #else
06545 matrix_file_read ( option_data->glt_filename[iglt] ,
06546 option_data->glt_rows[iglt] ,
06547 np , glt_cmat+iglt , 1 ) ;
06548 if( glt_cmat[iglt].elts == NULL ){
06549 fprintf(stderr,
06550 "** ERROR: Can't read GLT matrix from file %s\n",
06551 option_data->glt_filename[iglt] ) ;
06552 exit(1) ;
06553 }
06554 #endif
06555 }
06556
06557
06558
06559
06560
06561
06562 cxtxinvct = (matrix *) malloc( sizeof(matrix) * num_glt ) ;
06563 glt_amat = (matrix *) malloc( sizeof(matrix) * num_glt ) ;
06564 glt_coef = (vector *) malloc( sizeof(vector) * num_glt ) ;
06565 glt_tcoef = (vector *) malloc( sizeof(vector) * num_glt ) ;
06566
06567 for( iglt=0 ; iglt < num_glt ; iglt++ ){
06568 matrix_initialize( cxtxinvct+iglt ) ;
06569 matrix_initialize( glt_amat +iglt ) ;
06570 vector_initialize( glt_coef +iglt ) ;
06571 vector_initialize( glt_tcoef+iglt ) ;
06572 }
06573
06574 fglt = (float *) malloc( sizeof(float) * num_glt ) ;
06575 rglt = (float *) malloc( sizeof(float) * num_glt ) ;
06576
06577 vector_initialize( &y ) ; vector_create( nt, &y );
06578 vector_initialize( &coef ); vector_create( np, &coef );
06579
06580 if( dset_coef != NULL )
06581 cdar = (float *)malloc( sizeof(float) * DSET_NVALS(dset_coef) ) ;
06582
06583 init_glt_analysis( XtXinv , num_glt , glt_cmat , glt_amat , cxtxinvct ) ;
06584
06585
06586
06587 glt_coef_vol = (float ***) malloc( sizeof(float **) * num_glt ) ;
06588 if( tout )
06589 glt_tcoef_vol = (float ***) malloc( sizeof(float **) * num_glt ) ;
06590 if( fout )
06591 glt_fstat_vol = (float **) malloc( sizeof(float *) * num_glt ) ;
06592 if( rout )
06593 glt_rstat_vol = (float **) malloc( sizeof(float *) * num_glt ) ;
06594
06595 nvol = 0 ;
06596 for( iglt=0 ; iglt < num_glt ; iglt++ ){
06597 nlc = glt_rows[iglt];
06598 glt_coef_vol[iglt] = (float **) malloc( sizeof(float *) * nlc ) ;
06599 if( tout )
06600 glt_tcoef_vol[iglt] = (float **) malloc( sizeof(float *) * nlc ) ;
06601 for( ilc=0 ; ilc < nlc ; ilc++ ){
06602 glt_coef_vol[iglt][ilc] = (float *)calloc(sizeof(float),nxyz); nvol++;
06603 if( tout ){
06604 glt_tcoef_vol[iglt][ilc] = (float *)calloc(sizeof(float),nxyz); nvol++;
06605 }
06606 }
06607 if( fout ){
06608 glt_fstat_vol[iglt] = (float *)calloc( sizeof(float), nxyz ); nvol++;
06609 }
06610 if( rout ){
06611 glt_rstat_vol[iglt] = (float *)calloc( sizeof(float), nxyz ); nvol++;
06612 }
06613 }
06614
06615
06616
06617
06618
06619
06620
06621 vstep = nxyz / 50 ; if( !verb ) vstep = 0 ;
06622 if( vstep > 0 ) fprintf(stderr,"++ voxel loop:") ;
06623 for( ixyz=0 ; ixyz < nxyz ; ixyz++ ){
06624
06625 if( vstep > 0 && ixyz%vstep == vstep-1 ) vstep_print() ;
06626
06627
06628
06629 if( cget == nget ){
06630 if( imget != NULL ) DESTROY_IMARR(imget) ;
06631 iget[0] = ixyz ; nget = 1 ;
06632 for( jget=ixyz+1 ; jget < nxyz && nget < NGET ; jget++ )
06633 iget[nget++] = jget ;
06634 imget = THD_extract_many_series( nget, iget, dset_time ) ;
06635 cget = 0 ;
06636 }
06637
06638
06639
06640 ts_array = MRI_FLOAT_PTR(IMARR_SUBIM(imget,cget)) ; cget++ ;
06641 for( ii=0 ; ii < nt ; ii++ ) y.elts[ii] = ts_array[GoodList[ii]];
06642
06643
06644
06645 if( dset_coef != NULL ){
06646 (void) THD_extract_array( ixyz , dset_coef , 0 , cdar ) ;
06647 for( ii=0 ; ii < np ; ii++ ) coef.elts[ii] = cdar[ParamIndex[ii]] ;
06648 } else {
06649 vector_multiply( XtXinvXt , y , &coef ) ;
06650 }
06651
06652
06653
06654 for( ii=0 ; ii < np ; ii++ ) if( coef.elts[ii] != 0.0 ) break ;
06655 novar = (ii==np) ;
06656
06657 if( !novar ) ssef = calc_sse( X , coef , y ) ;
06658
06659
06660
06661
06662
06663 glt_analysis( nt , np ,
06664 X , y , ssef , coef , novar ,
06665 cxtxinvct , num_glt , glt_rows , glt_cmat , glt_amat ,
06666 glt_coef , glt_tcoef , fglt , rglt ) ;
06667
06668
06669
06670 if (glt_coef_vol != NULL)
06671 for (iglt = 0; iglt < num_glt; iglt++)
06672 if (glt_coef_vol[iglt] != NULL)
06673 for (ilc = 0; ilc < glt_rows[iglt]; ilc++)
06674 glt_coef_vol[iglt][ilc][ixyz] = glt_coef[iglt].elts[ilc];
06675
06676 if (glt_tcoef_vol != NULL)
06677 for (iglt = 0; iglt < num_glt; iglt++)
06678 if (glt_tcoef_vol[iglt] != NULL)
06679 for (ilc = 0; ilc < glt_rows[iglt]; ilc++)
06680 glt_tcoef_vol[iglt][ilc][ixyz] = glt_tcoef[iglt].elts[ilc];
06681
06682 if (glt_fstat_vol != NULL)
06683 for (iglt = 0; iglt < num_glt; iglt++)
06684 if (glt_fstat_vol[iglt] != NULL)
06685 glt_fstat_vol[iglt][ixyz] = fglt[iglt];
06686
06687 if (glt_rstat_vol != NULL)
06688 for (iglt = 0; iglt < num_glt; iglt++)
06689 if (glt_rstat_vol[iglt] != NULL)
06690 glt_rstat_vol[iglt][ixyz] = rglt[iglt];
06691
06692 }
06693
06694 if( vstep > 0 ) fprintf(stderr,"\n") ;
06695
06696
06697
06698 DSET_unload( dset_time ) ;
06699 if( dset_coef != NULL ) DSET_delete( dset_coef ) ;
06700
06701
06702
06703
06704
06705
06706 if( option_data->bucket_filename != NULL ){
06707 buck_prefix = strdup(option_data->bucket_filename) ;
06708 } else if( BucketFilename != NULL ){
06709 buck_prefix = (char *)malloc(strlen(BucketFilename)+4) ;
06710 FILENAME_TO_PREFIX( BucketFilename , buck_prefix ) ;
06711 } else {
06712 buck_prefix = strdup("X") ;
06713 }
06714
06715
06716
06717 buck_name = (char *)malloc(strlen(buck_prefix)+32) ;
06718 strcpy(buck_name,buck_prefix) ;
06719 strcat(buck_name,"+") ; strcat(buck_name,VIEW_codestr[view]) ;
06720
06721 dset_buck = THD_open_dataset( buck_name ) ;
06722
06723 if( dset_buck != NULL ){
06724
06725 if( verb) fprintf(stderr,"++ -xrestore appending to dataset %s\n",buck_name) ;
06726 DSET_mallocize( dset_buck ) ;
06727 if( DSET_NVOX(dset_buck) != nxyz ){
06728 fprintf(stderr,
06729 "** ERROR: dataset %s mismatch with time series '%s'\n" ,
06730 buck_name , DSET_HEADNAME(dset_time) );
06731 exit(0) ;
06732 }
06733 DSET_load( dset_buck ) ;
06734 if( !DSET_LOADED(dset_buck) ){
06735 fprintf(stderr,
06736 "** ERROR: can't load dataset %s from disk\n" , buck_name ) ;
06737 exit(1) ;
06738 }
06739
06740
06741
06742 ivol = DSET_NVALS(dset_buck) ;
06743
06744 EDIT_add_bricklist( dset_buck , nvol, NULL, NULL, NULL ) ;
06745
06746 } else {
06747
06748 if( verb ) fprintf(stderr,"++ -xrestore creating new dataset %s\n",buck_name) ;
06749
06750 dset_buck = EDIT_empty_copy( dset_time ) ;
06751 (void) EDIT_dset_items( dset_buck ,
06752 ADN_prefix, buck_prefix ,
06753 ADN_type, HEAD_FUNC_TYPE,
06754 ADN_func_type, FUNC_BUCK_TYPE,
06755 ADN_datum_all, MRI_short ,
06756 ADN_ntt, 0,
06757 ADN_nvals, nvol ,
06758 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
06759 ADN_none ) ;
06760 ivol = 0 ;
06761 }
06762
06763 tross_Make_History( PROGRAM_NAME , argc , argv , dset_buck ) ;
06764
06765 nbuck = DSET_NVALS(dset_buck) ;
06766
06767
06768
06769 for( iglt=0 ; iglt < num_glt ; iglt++ ){
06770
06771 for( ilc=0 ; ilc < glt_rows[iglt] ; ilc++ ){
06772 sprintf( brick_label , "%s LC[%d] coef" , glt_label[iglt], ilc ) ;
06773 volume = glt_coef_vol[iglt][ilc];
06774 attach_sub_brick( dset_buck, ivol, volume, nxyz,
06775 FUNC_FIM_TYPE, brick_label, 0, 0, 0, NULL);
06776 free((void *)volume) ; ivol++ ;
06777
06778 if( tout ){
06779 sprintf( brick_label , "%s LC[%d] t-st" , glt_label[iglt], ilc ) ;
06780 volume = glt_tcoef_vol[iglt][ilc];
06781 attach_sub_brick( dset_buck, ivol, volume, nxyz,
06782 FUNC_TT_TYPE, brick_label, nt-np, 0, 0, NULL);
06783 free((void *)volume) ; ivol++ ;
06784 }
06785 }
06786
06787 if( rout ){
06788 sprintf( brick_label , "%s R^2" , glt_label[iglt] ) ;
06789 volume = glt_rstat_vol[iglt];
06790 attach_sub_brick( dset_buck, ivol, volume, nxyz,
06791 FUNC_THR_TYPE, brick_label, 0, 0, 0, NULL);
06792 free((void *)volume) ; ivol++ ;
06793 }
06794
06795 if( fout ){
06796 sprintf( brick_label , "%s F-stat" , glt_label[iglt] ) ;
06797 volume = glt_fstat_vol[iglt];
06798 attach_sub_brick( dset_buck, ivol, volume, nxyz,
06799 FUNC_FT_TYPE, brick_label, 0, glt_rows[iglt], nt-np, NULL);
06800 free((void *)volume) ; ivol++ ;
06801 }
06802 }
06803
06804
06805
06806 DSET_write( dset_buck ) ;
06807 return ;
06808 }
06809
06810
06811
06812
06813
06814 #define GLT_ERR \
06815 do{ fprintf(stderr,"** ERROR: Can't read GLT matrix from file %s\n",fname); \
06816 exit(1) ; } while(0)
06817
06818 void read_glt_matrix( char *fname, int *nrows, int ncol, matrix *cmat )
06819 {
06820 int ii,jj ;
06821
06822 ENTRY("read_glt_matrix") ;
06823
06824 if( *nrows > 0 ){
06825
06826 matrix_file_read( fname , *nrows , ncol , cmat , 1 ) ;
06827 if( cmat->elts == NULL ) GLT_ERR ;
06828
06829 } else {
06830 floatvecvec *fvv ;
06831 float **far=NULL ;
06832 int nr=0 , iv ;
06833
06834 if( nSymStim < 1 ){
06835 fprintf(stderr,"** ERROR: use of -gltsym without SymStim being defined\n");
06836 exit(1) ;
06837 }
06838
06839 if( strncmp(fname,"SYM:",4) == 0 ){
06840 char *fdup=strdup(fname+4) , *fpt , *buf ;
06841 int ss , ns ;
06842
06843 buf = fdup ;
06844 while(1){
06845 fpt = strchr(buf,'\\') ;
06846 if( fpt != NULL ) *fpt = '\0' ;
06847 fvv = SYM_expand_ranges( ncol-1 , nSymStim,SymStim , buf ) ;
06848 if( fvv == NULL || fvv->nvec < 1 ) continue ;
06849 far = (float **)realloc((void *)far , sizeof(float *)*(nr+fvv->nvec)) ;
06850 for( iv=0 ; iv < fvv->nvec ; iv++ ) far[nr++] = fvv->fvar[iv].ar ;
06851 free((void *)fvv->fvar) ; free((void *)fvv) ;
06852 if( fpt == NULL ) break ;
06853 buf = fpt+1 ;
06854 }
06855 free((void *)fdup) ;
06856
06857 } else {
06858 char buf[8192] , *cpt ;
06859 FILE *fp = fopen( fname , "r" ) ;
06860 if( fp == NULL ) GLT_ERR ;
06861 while(1){
06862 cpt = fgets( buf , 8192 , fp ) ;
06863 if( cpt == NULL ) break ;
06864 fvv = SYM_expand_ranges( ncol-1 , nSymStim,SymStim , buf ) ;
06865 if( fvv == NULL || fvv->nvec < 1 ) continue ;
06866 far = (float **)realloc((void *)far , sizeof(float *)*(nr+fvv->nvec)) ;
06867 for( iv=0 ; iv < fvv->nvec ; iv++ ) far[nr++] = fvv->fvar[iv].ar ;
06868 free((void *)fvv->fvar) ; free((void *)fvv) ;
06869 }
06870 fclose(fp) ;
06871 }
06872 if( nr == 0 ) GLT_ERR ;
06873 *nrows = nr ;
06874 array_to_matrix( nr , ncol , far , cmat ) ;
06875
06876 for( ii=0 ; ii < nr ; ii++ ) free((void *)far[ii]) ;
06877 free((void *)far) ;
06878
06879 if( !AFNI_noenv("AFNI_GLTSYM_PRINT") ){
06880 printf("GLT matrix from '%s':\n",fname) ;
06881 matrix_print( *cmat ) ;
06882 }
06883 }
06884
06885
06886
06887 for( ii=0 ; ii < *nrows ; ii++ ){
06888 for( jj=0 ; jj < ncol && cmat->elts[ii][jj] == 0.0 ; jj++ ) ;
06889 if( jj == ncol )
06890 fprintf(stderr,"** ERROR: row #%d of matrix '%s' is all zero!\n",
06891 ii+1 , fname ) ;
06892 }
06893
06894 EXRETURN ;
06895 }
06896
06897
06898
06899 static void vstep_print(void)
06900 {
06901 static int nn=0 ;
06902 static char xx[10] = "0123456789" ;
06903 fprintf(stderr , "%c" , xx[nn%10] ) ;
06904 if( nn%10 == 9) fprintf(stderr,".") ;
06905 nn++ ;
06906 }
06907
06908
06909
06910
06911
06912
06913
06914
06915
06916
06917 float basis_evaluation( basis_expansion *be , float *wt , float x )
06918 {
06919 float sum=0.0 ;
06920 int ii ;
06921
06922 if( x >= be->tbot && x <= be->ttop ){
06923 for( ii=0 ; ii < be->nfunc ; ii++ )
06924 sum += wt[ii] * basis_funceval( be->bfunc[ii] , x ) ;
06925 }
06926
06927 return sum ;
06928 }
06929
06930
06931
06932
06933
06934
06935
06936 static float basis_tent( float x, float bot, float mid, float top, void *q )
06937 {
06938 if( x <= bot || x >= top ) return 0.0f ;
06939 if( x <= mid ) return (x-bot)/(mid-bot) ;
06940 return (top-x)/(top-mid) ;
06941 }
06942
06943
06944
06945
06946
06947 static float basis_one( float x, float bot, float top, float junk, void *q )
06948 {
06949 if( x < bot || x > top ) return 0.0f ;
06950 return 1.0f ;
06951 }
06952
06953
06954
06955
06956
06957
06958
06959 static float basis_cos( float x, float bot, float top, float freq, void *q )
06960 {
06961 if( x < bot || x > top ) return 0.0f ;
06962 return (float)cos(freq*(x-bot)) ;
06963 }
06964
06965
06966
06967
06968
06969
06970
06971 static float basis_sin( float x, float bot, float top, float freq, void *q )
06972 {
06973 if( x <= bot || x >= top ) return 0.0f ;
06974 return (float)sin(freq*(x-bot)) ;
06975 }
06976
06977
06978
06979
06980
06981
06982
06983 static float basis_gam( float x, float b, float c, float top, void *q )
06984 {
06985 if( x <= 0.0f || x > top ) return 0.0f ;
06986 return (float)(pow(x/(b*c),b)*exp(b-x/c)) ;
06987 }
06988
06989
06990
06991 #undef SPM_A1
06992 #undef SPM_A2
06993 #undef SPM_P1
06994 #undef SPM_P2
06995
06996 #define SPM_A1 0.0083333333
06997 #define SPM_P1 4.0
06998 #define SPM_A2 1.274527e-13
06999 #define SPM_P2 15.0
07000
07001 static float basis_spmg1( float x, float a, float b, float c, void *q )
07002 {
07003 if( x <= 0.0f || x >= 25.0f ) return 0.0f ;
07004 return (float)(exp(-x)*( SPM_A1*pow(x,SPM_P1)
07005 -SPM_A2*pow(x,SPM_P2) )) ;
07006 }
07007
07008 static float basis_spmg2( float x, float a, float b, float c, void *q )
07009 {
07010 if( x <= 0.0f || x >= 25.0f ) return 0.0f ;
07011 return (float)(exp(-x)*( SPM_A1*pow(x,SPM_P1-1.0)*(SPM_P1-x)
07012 -SPM_A2*pow(x,SPM_P2-1.0)*(SPM_P2-x) )) ;
07013 }
07014
07015
07016
07017
07018
07019
07020
07021 #undef TPEAK4
07022 #define TPEAK4(TT) ((TT)/(1.0-exp(-0.25*(TT))))
07023
07024 static float basis_block_hrf4( float tt , float TT )
07025 {
07026 register double t26, t2, t4, t1, t42, t12, t34, t35, t16, t46, t,L ;
07027 double w ;
07028
07029 if( tt <= 0.0f || tt >= (TT+15.0f) ) return 0.0f ;
07030
07031 t = tt ; L = TT ; t4 = exp(0.4e1 - t);
07032 if( t < L ){ L = t ; t16 = 54.5982 ; }
07033 else { t16 = exp(4.0-t+L) ; }
07034
07035 t1 = t * t;
07036 t2 = t1 * t1;
07037 t4 = exp(4.0 - t);
07038 t12 = t1 * t;
07039 t26 = t16 * L;
07040 t34 = L * L;
07041 t35 = t16 * t34;
07042 t42 = t16 * t34 * L;
07043 t46 = t34 * t34;
07044
07045 w = -t2 * t4 / 0.256e3 - 0.3e1 / 0.32e2 * t4 - 0.3e1 / 0.32e2 * t4 * t
07046 - 0.3e1 / 0.64e2 * t4 * t1 - t4 * t12 / 0.64e2 + t16 * t2 / 0.256e3
07047 + 0.3e1 / 0.32e2 * t16 + 0.3e1 / 0.32e2 * t16 * t
07048 + 0.3e1 / 0.64e2 * t1 * t16 + t16 * t12 / 0.64e2 - 0.3e1 / 0.32e2 * t26
07049 - 0.3e1 / 0.32e2 * t26 * t - 0.3e1 / 0.64e2 * t1 * t26
07050 - t26 * t12 / 0.64e2 + 0.3e1 / 0.64e2 * t35 + 0.3e1 / 0.64e2 * t35 * t
07051 + 0.3e1 / 0.128e3 * t1 * t35 - t42 / 0.64e2 - t42 * t / 0.64e2
07052 + t16 * t46 / 0.256e3 ;
07053 return (float)w ;
07054 }
07055
07056 static float basis_block4( float t, float T, float peak, float junk, void *q )
07057 {
07058 float w , tp , pp ;
07059
07060 w = basis_block_hrf4(t,T) ;
07061 if( w > 0.0f && peak > 0.0f ){
07062 tp = TPEAK4(T) ; pp = basis_block_hrf4(tp,T) ;
07063 if( pp > 0.0f ) w *= peak / pp ;
07064 }
07065 return w ;
07066 }
07067
07068
07069
07070
07071
07072
07073
07074 #undef TPEAK5
07075 #define TPEAK5(TT) ((TT)/(1.0-exp(-0.2*(TT))))
07076
07077 static float basis_block_hrf5( float tt, float TT )
07078 {
07079 register double t , T ;
07080 register double t2,t3,t4,t5,t6,t7,t9,t10,t11,t14,t20,t25,t28,t37,t57 ;
07081 double w ;
07082
07083 if( tt <= 0.0f || tt >= (TT+15.0f) ) return 0.0f ;
07084
07085 t = tt ; T = TT ;
07086
07087 #if 1
07088 t2 = exp(-t) ;
07089 if( t <= T ){ t3 = t ; t4 = 1.0/t2 ; }
07090 else { t3 = T ; t4 = exp(T) ; }
07091 t2 *= 148.413 ;
07092 #else
07093 t2 = exp(0.5e1 - t);
07094 t3 = (t <= T ? t : T);
07095 t4 = exp(t3);
07096 #endif
07097 t5 = t * t;
07098 t6 = t5 * t5;
07099 t7 = t6 * t;
07100 t9 = t3 * t3;
07101 t10 = t9 * t9;
07102 t11 = t4 * t10;
07103 t14 = t4 * t9 * t3;
07104 t20 = t4 * t3;
07105 t25 = t4 * t9;
07106 t28 = t5 * t;
07107 t37 = -0.120e3 + t4 * t7 + 0.5e1 * t11 - 0.20e2 * t14 - t4 * t10 * t3
07108 - 0.10e2 * t14 * t5 - 0.120e3 * t20 * t - 0.20e2 * t14 * t
07109 + 0.30e2 * t25 * t5 + 0.10e2 * t25 * t28 + 0.5e1 * t11 * t
07110 + 0.20e2 * t4 * t28 + 0.60e2 * t25 * t;
07111 t57 = -0.5e1 * t20 * t6 - 0.20e2 * t20 * t28 - 0.60e2 * t20 * t5
07112 - 0.5e1 * t6 - 0.20e2 * t28 + 0.120e3 * t4 - 0.120e3 * t
07113 - 0.120e3 * t20 + 0.60e2 * t25 - t7 - 0.60e2 * t5 + 0.120e3 * t4 * t
07114 + 0.60e2 * t4 * t5 + 0.5e1 * t4 * t6;
07115 w = t2 * (t37 + t57) / 0.3125e4;
07116
07117 return (float)w ;
07118 }
07119
07120 static float basis_block5( float t, float T, float peak, float junk, void *q )
07121 {
07122 float w , tp , pp ;
07123
07124 w = basis_block_hrf5(t,T) ;
07125 if( w > 0.0f && peak > 0.0f ){
07126 tp = TPEAK5(T) ; pp = basis_block_hrf5(tp,T) ;
07127 if( pp > 0.0f ) w *= peak / pp ;
07128 }
07129 return w ;
07130 }
07131
07132
07133
07134
07135
07136
07137
07138 static float basis_legendre( float x, float bot, float top, float n, void *q )
07139 {
07140 float xq ;
07141
07142 x = 2.0f*(x-bot)/(top-bot) - 1.0f ;
07143
07144 if( x < -1.0f || x > 1.0f ) return 0.0f ;
07145
07146 xq = x*x ;
07147
07148 switch( (int)n ){
07149 case 0: return 1.0f ;
07150 case 1: return x ;
07151 case 2: return (3.0f*xq-1.0f)/2.0f ;
07152 case 3: return (5.0f*xq-3.0f)*x/2.0f ;
07153 case 4: return ((35.0f*xq-30.0f)*xq+3.0f)/8.0f ;
07154 case 5: return ((63.0f*xq-70.0f)*xq+15.0f)*x/8.0f ;
07155 case 6: return (((231.0f*xq-315.0f)*xq+105.0f)*xq-5.0f)/16.0f ;
07156 case 7: return (((429.0f*xq-693.0f)*xq+315.0f)*xq-35.0f)*x/16.0f ;
07157
07158 case 8: return ((((6435.0f*xq-12012.0f)*xq+6930.0f)*xq-1260.0f)*xq+35.0f)
07159 /128.0f;
07160
07161 case 9: return ((((12155.0f*xq-25740.0f)*xq+18018.0f)*xq-4620.0f)*xq+315.0f)
07162 *x/128.0f ;
07163 }
07164
07165 return 0.0f ;
07166 }
07167
07168 #undef POLY_MAX
07169 #define POLY_MAX 9
07170
07171
07172
07173 #define ITT 19
07174 #define IXX 23
07175 #define IZZ 25
07176
07177
07178
07179
07180 static float basis_expr( float x, float bot, float top, float dtinv, void *q )
07181 {
07182 PARSER_code *pc = (PARSER_code *)q ;
07183 double atoz[26] , val ;
07184
07185 if( x < bot || x > top ) return 0.0f ;
07186 atoz[ITT] = x ;
07187 atoz[IXX] = (x-bot)*dtinv ;
07188 atoz[IZZ] = 2.0*atoz[IXX] - 1.0 ;
07189 val = PARSER_evaluate_one( pc , atoz ) ;
07190 return (float)val ;
07191 }
07192
07193
07194
07195
07196
07197 basis_expansion * basis_parser( char *sym )
07198 {
07199 basis_expansion *be ;
07200 char *cpt , *scp ;
07201 float bot=0.0f, top=0.0f ;
07202 int nn , nord=0 ;
07203
07204 if( sym == NULL ) return NULL ;
07205
07206 scp = strdup(sym) ;
07207 cpt = strchr(scp,'(') ;
07208 if( cpt != NULL ){ *cpt = '\0' ; cpt++ ; }
07209
07210 be = (basis_expansion *)malloc(sizeof(basis_expansion)) ;
07211 be->name = NULL ;
07212
07213
07214
07215 if( strcmp(scp,"GAM") == 0 ){
07216
07217 be->nfunc = 1 ;
07218 be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07219 be->bfunc[0].f = basis_gam ;
07220 if( cpt == NULL ){
07221 be->bfunc[0].a = 8.6f ;
07222 be->bfunc[0].b = 0.547f ;
07223 be->bfunc[0].c = 11.1f ;
07224 } else {
07225 sscanf(cpt,"%f,%f",&bot,&top) ;
07226 if( bot <= 0.0f || top <= 0.0f ){
07227 fprintf(stderr,"** ERROR: 'GAM(%s' is illegal\n",cpt) ;
07228 fprintf(stderr,
07229 " Correct format: 'GAM(b,c)' with b > 0 and c > 0.\n");
07230 free((void *)be->bfunc); free((void *)be); free(scp); return NULL;
07231 }
07232 be->bfunc[0].a = bot ;
07233 be->bfunc[0].b = top ;
07234 be->bfunc[0].c = bot*top + 4.0f*sqrt(bot)*top ;
07235 }
07236 be->tbot = 0.0f ; be->ttop = be->bfunc[0].c ;
07237
07238
07239
07240 } else if( strcmp(scp,"TENT") == 0 ){
07241 float dx ;
07242
07243 if( cpt == NULL ){
07244 fprintf(stderr,"** ERROR: 'TENT' by itself is illegal\n") ;
07245 fprintf(stderr,
07246 " Correct format: 'TENT(bot,top,n)' with bot < top and n > 0.\n") ;
07247 free((void *)be); free(scp); return NULL;
07248 }
07249 sscanf(cpt,"%f,%f,%d",&bot,&top,&nord) ;
07250 if( bot >= top || nord < 2 ){
07251 fprintf(stderr,"** ERROR: 'TENT(%s' is illegal\n",cpt) ;
07252 fprintf(stderr,
07253 " Correct format: 'TENT(bot,top,n)' with bot < top and n > 1.\n") ;
07254 free((void *)be); free(scp); return NULL;
07255 }
07256 be->nfunc = nord ;
07257 be->tbot = bot ; be->ttop = top ;
07258 be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07259 dx = (top-bot) / (nord-1) ;
07260
07261 be->bfunc[0].f = basis_tent ;
07262 be->bfunc[0].a = bot-0.001*dx ;
07263 be->bfunc[0].b = bot ;
07264 be->bfunc[0].c = bot+dx ;
07265 for( nn=1 ; nn < nord-1 ; nn++ ){
07266 be->bfunc[nn].f = basis_tent ;
07267 be->bfunc[nn].a = bot + (nn-1)*dx ;
07268 be->bfunc[nn].b = bot + nn *dx ;
07269 be->bfunc[nn].c = bot + (nn+1)*dx ;
07270 }
07271 be->bfunc[nord-1].f = basis_tent ;
07272 be->bfunc[nord-1].a = bot + (nord-2)*dx ;
07273 be->bfunc[nord-1].b = top ;
07274 be->bfunc[nord-1].c = top + 0.001*dx ;
07275
07276
07277
07278 } else if( strcmp(scp,"TRIG") == 0 ){
07279
07280 if( cpt == NULL ){
07281 fprintf(stderr,"** ERROR: 'TRIG' by itself is illegal\n") ;
07282 fprintf(stderr,
07283 " Correct format: 'TRIG(bot,top,n)' with bot < top and n > 2.\n") ;
07284 free((void *)be); free(scp); return NULL;
07285 }
07286 sscanf(cpt,"%f,%f,%d",&bot,&top,&nord) ;
07287 if( bot >= top || nord < 3 ){
07288 fprintf(stderr,"** ERROR: 'TRIG(%s' is illegal\n",cpt) ;
07289 fprintf(stderr,
07290 " Correct format: 'TRIG(bot,top,n)' with bot < top and n > 2.\n") ;
07291 free((void *)be); free(scp); return NULL;
07292 }
07293 be->nfunc = nord ;
07294 be->tbot = bot ; be->ttop = top ;
07295 be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07296
07297 be->bfunc[0].f = basis_one ;
07298 be->bfunc[0].a = bot ;
07299 be->bfunc[0].b = top ;
07300 be->bfunc[0].c = 0.0f ;
07301 for( nn=1 ; nn < nord ; nn++ ){
07302 be->bfunc[nn].f = basis_cos ;
07303 be->bfunc[nn].a = bot ;
07304 be->bfunc[nn].b = top ;
07305 be->bfunc[nn].c = (2.0*PI)*((nn+1)/2)/(top-bot) ;
07306 nn++ ; if( nn >= nord ) break ;
07307 be->bfunc[nn].f = basis_sin ;
07308 be->bfunc[nn].a = bot ;
07309 be->bfunc[nn].b = top ;
07310 be->bfunc[nn].c = (2.0*PI)*((nn+1)/2)/(top-bot) ;
07311 }
07312
07313
07314
07315 } else if( strcmp(scp,"SIN") == 0 ){
07316
07317 if( cpt == NULL ){
07318 fprintf(stderr,"** ERROR: 'SIN' by itself is illegal\n") ;
07319 fprintf(stderr,
07320 " Correct format: 'SIN(bot,top,n)' with bot < top and n > 0.\n") ;
07321 free((void *)be); free(scp); return NULL;
07322 }
07323 sscanf(cpt,"%f,%f,%d",&bot,&top,&nord) ;
07324 if( bot >= top || nord < 1 ){
07325 fprintf(stderr,"** ERROR: 'SIN(%s' is illegal\n",cpt) ;
07326 fprintf(stderr,
07327 " Correct format: 'SIN(bot,top,n)' with bot < top and n > 0.\n") ;
07328 free((void *)be); free(scp); return NULL;
07329 }
07330 be->nfunc = nord ;
07331 be->tbot = bot ; be->ttop = top ;
07332 be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07333
07334 for( nn=0 ; nn < nord ; nn++ ){
07335 be->bfunc[nn].f = basis_sin ;
07336 be->bfunc[nn].a = bot ;
07337 be->bfunc[nn].b = top ;
07338 be->bfunc[nn].c = PI*(nn+1)/(top-bot) ;
07339 }
07340
07341
07342
07343 } else if( strcmp(scp,"POLY") == 0 ){
07344
07345 if( cpt == NULL ){
07346 fprintf(stderr,"** ERROR: 'POLY' by itself is illegal\n") ;
07347 fprintf(stderr,
07348 " Correct format: 'POLY(bot,top,n)' with bot < top and n > 0.\n") ;
07349 free((void *)be); free(scp); return NULL;
07350 }
07351 sscanf(cpt,"%f,%f,%d",&bot,&top,&nord) ;
07352 if( bot >= top || nord < 1 || nord > POLY_MAX ){
07353 fprintf(stderr,"** ERROR: 'POLY(%s' is illegal\n",cpt) ;
07354 fprintf(stderr,
07355 " Correct format: 'POLY(bot,top,n)' with bot < top and n > 0.\n") ;
07356 free((void *)be); free(scp); return NULL;
07357 }
07358 be->nfunc = nord ;
07359 be->tbot = bot ; be->ttop = top ;
07360 be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07361 be->bfunc[0].f = basis_one ;
07362 be->bfunc[0].a = bot ;
07363 be->bfunc[0].b = top ;
07364 be->bfunc[0].c = 0.0f ;
07365 for( nn=1 ; nn < nord ; nn++ ){
07366 be->bfunc[nn].f = basis_legendre ;
07367 be->bfunc[nn].a = bot ;
07368 be->bfunc[nn].b = top ;
07369 be->bfunc[nn].c = (float)nn ;
07370 }
07371
07372
07373
07374 } else if( strncmp(scp,"SPMG",4) == 0 ){
07375
07376 be->nfunc = 2 ;
07377 be->tbot = 0.0f ; be->ttop = 25.0f ;
07378 be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07379 be->bfunc[0].f = basis_spmg1 ;
07380 be->bfunc[0].a = 0.0f ;
07381 be->bfunc[0].b = 0.0f ;
07382 be->bfunc[0].c = 0.0f ;
07383 be->bfunc[1].f = basis_spmg2 ;
07384 be->bfunc[1].a = 0.0f ;
07385 be->bfunc[1].b = 0.0f ;
07386 be->bfunc[1].c = 0.0f ;
07387
07388
07389
07390 } else if( strncmp(scp,"BLOCK",5) == 0 || strncmp(scp,"IGFUN",5) == 0 ){
07391 int nb=4 ;
07392
07393 if( scp[5] != '\0' ){
07394 nb = strtol( scp+5 , NULL , 10 ) ;
07395 if( nb != 4 && nb != 5 ){
07396 fprintf(stderr,
07397 "** ERROR: '%s' has illegal power: only 4 or 5 allowed\n",scp) ;
07398 free((void *)be); free(scp); return NULL;
07399 }
07400 }
07401
07402 if( cpt == NULL ){
07403 fprintf(stderr,"** ERROR: '%s' by itself is illegal\n",scp) ;
07404 fprintf(stderr,
07405 " Correct format: 'BLOCKn(dur)' with dur > 0, n=4 or 5\n"
07406 " *OR* 'BLOCKn(dur,peak)' with peak > 0 too.\n") ;
07407 free((void *)be); free(scp); return NULL;
07408 }
07409 sscanf(cpt,"%f,%f",&top,&bot) ;
07410 if( top <= 0.0f ){
07411 fprintf(stderr,"** ERROR: '%s(%s' is illegal\n",scp,cpt) ;
07412 fprintf(stderr,
07413 " Correct format: 'BLOCKn(dur)' with dur > 0, n=4 or 4\n"
07414 " or: 'BLOCKn(dur,peak)' with peak > 0 too.\n") ;
07415 free((void *)be); free(scp); return NULL;
07416 }
07417 if( bot < 0.0f ) bot = 0.0f ;
07418
07419 be->nfunc = 1 ;
07420 be->tbot = 0.0f ; be->ttop = top+15.0f ;
07421 be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07422 be->bfunc[0].f = (nb==5) ? basis_block5 : basis_block4 ;
07423 be->bfunc[0].a = top ;
07424 be->bfunc[0].b = bot ;
07425 be->bfunc[0].c = 0.0f ;
07426
07427
07428
07429 } else if( strcmp(scp,"EXPR") == 0 ){
07430 char *ept ;
07431 NI_str_array *sar ;
07432 int nexpr , ie ;
07433 PARSER_code *pc ;
07434
07435 if( cpt == NULL ){
07436 fprintf(stderr,"** ERROR: 'EXPR' by itself is illegal\n") ;
07437 free((void *)be); free(scp); return NULL;
07438 }
07439 sscanf(cpt,"%f,%f",&bot,&top) ;
07440 if( top <= bot ){
07441 fprintf(stderr,"** ERROR: 'EXPR(%f,%f)' has illegal time range\n",bot,top) ;
07442 free((void *)be); free(scp); return NULL;
07443 }
07444 ept = strchr( cpt , ')' ) ;
07445 if( ept == NULL ){
07446 fprintf(stderr,"** ERROR: 'EXPR(%f,%f)' has no expressions!?\n",bot,top);
07447 free((void *)be); free(scp); return NULL;
07448 }
07449 sar = NI_decode_string_list( ept+1 , "~" ) ;
07450 if( sar == NULL || sar->num == 0 ){
07451 fprintf(stderr,"** ERROR: 'EXPR(%f,%f)' has no expressions!?\n",bot,top);
07452 free((void *)be); free(scp); return NULL;
07453 }
07454 be->nfunc = nexpr = sar->num ;
07455 be->tbot = bot ; be->ttop = top ;
07456 be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07457 PARSER_set_printout(1) ;
07458 for( ie=0 ; ie < nexpr ; ie++ ){
07459 pc = PARSER_generate_code( sar->str[ie] ) ;
07460 if( pc == NULL ){
07461 fprintf(stderr,"** ERROR: unparsable EXPRession '%s'\n",sar->str[ie]) ;
07462 free((void *)be); free(scp); return NULL;
07463 }
07464 if( strcmp(sar->str[ie],"1") != 0 &&
07465 !PARSER_has_symbol("t",pc) &&
07466 !PARSER_has_symbol("x",pc) &&
07467 !PARSER_has_symbol("z",pc) ){
07468 fprintf(stderr,
07469 "** ERROR: illegal EXPRession '%s' lacks symbol t, x, or z\n",
07470 sar->str[ie]) ;
07471 free((void *)be); free(scp); return NULL;
07472 }
07473 be->bfunc[ie].f = basis_expr ;
07474 be->bfunc[ie].a = bot ;
07475 be->bfunc[ie].b = top ;
07476 be->bfunc[ie].c = 1.0f/(top-bot) ;
07477 be->bfunc[ie].q = (void *)pc ;
07478 }
07479 PARSER_set_printout(0) ;
07480
07481 NI_delete_str_array(sar) ;
07482
07483
07484
07485 } else {
07486 fprintf(stderr,"** ERROR: '%s' is unknown response function type\n",scp) ;
07487 free((void *)be); free(scp); return NULL;
07488 }
07489
07490
07491
07492 for( nn=0 ; nn < be->nfunc ; nn++ )
07493 be->bfunc[nn].ffac = 1.0 ;
07494
07495
07496
07497
07498 #undef BNSUB
07499 #define BNSUB 999
07500 if( basis_normall > 0.0f ){
07501 int jj ; float dt , ftop , val ;
07502 bot = be->tbot ; top = be->ttop ; dt = (top-bot)/BNSUB ;
07503 for( nn=0 ; nn < be->nfunc ; nn++ ){
07504 ftop = 0.0f ;
07505 for( jj=0 ; jj <= BNSUB ; jj++ ){
07506 val = basis_funceval( be->bfunc[nn] , bot+jj*dt ) ;
07507 val = fabs(val) ; if( val > ftop ) ftop = val ;
07508 }
07509 if( ftop > 0.0f ) be->bfunc[nn].ffac = basis_normall / ftop ;
07510 }
07511 }
07512
07513 free(scp); return be;
07514 }
07515
07516
07517
07518
07519
07520
07521 void basis_write_iresp( int argc , char *argv[] ,
07522 DC_options *option_data ,
07523 basis_expansion *be , float dt ,
07524 float **wtar , char *output_filename )
07525 {
07526 int nvox, ii, nf, allz, ts_length ;
07527 register int pp, ib ;
07528 float *wt , *tt , **hout , factor , **bb ;
07529 short *bar ;
07530 THD_3dim_dataset *in_dset = NULL;
07531 THD_3dim_dataset *out_dset = NULL;
07532 char *commandline , label[512] ;
07533 const float EPSILON = 1.0e-10 ;
07534
07535
07536
07537 in_dset = THD_open_dataset(option_data->input_filename);
07538 nvox = in_dset->daxes->nxx * in_dset->daxes->nyy * in_dset->daxes->nzz;
07539
07540 if( dt <= 0.0f ) dt = DSET_TR(in_dset) ;
07541 if( dt <= 0.0f ) dt = 1.0f ;
07542
07543
07544
07545 out_dset = EDIT_empty_copy( in_dset ) ;
07546 tross_Copy_History( in_dset , out_dset ) ;
07547 DSET_delete( in_dset ) ;
07548
07549
07550
07551 commandline = tross_commandline( PROGRAM_NAME , argc , argv ) ;
07552 sprintf( label , "Impulse response: %s" , output_filename ) ;
07553 if( commandline != NULL )
07554 tross_multi_Append_History( out_dset , commandline,label,NULL ) ;
07555 else
07556 tross_Append_History ( out_dset, label);
07557 free((void *)commandline) ;
07558
07559 ts_length = 1 + (int)ceil( (be->ttop - be->tbot)/dt ) ;
07560
07561
07562
07563 (void ) EDIT_dset_items( out_dset,
07564 ADN_prefix, output_filename,
07565 ADN_label1, output_filename,
07566 ADN_self_name, output_filename,
07567 ADN_malloc_type, DATABLOCK_MEM_MALLOC,
07568 ADN_datum_all, MRI_short,
07569 ADN_nvals, ts_length,
07570 ADN_ntt, ts_length,
07571 ADN_ttdel, dt ,
07572 ADN_ttorg, be->tbot,
07573 ADN_none ) ;
07574
07575 if( THD_is_file(out_dset->dblk->diskptr->header_name) ){
07576 fprintf(stderr,
07577 "** ERROR: Output dataset file %s already exists - won't overwrite\n",
07578 out_dset->dblk->diskptr->header_name ) ;
07579 DSET_delete(out_dset) ; return ;
07580 }
07581
07582
07583
07584 hout = (float **) malloc( sizeof(float *) * ts_length ) ;
07585 for( ib=0 ; ib < ts_length ; ib++ )
07586 hout[ib] = (float *)calloc(sizeof(float),nvox) ;
07587
07588
07589
07590 nf = be->nfunc ;
07591 wt = (float *) malloc( sizeof(float) * nf ) ;
07592 tt = (float *) malloc( sizeof(float) * ts_length ) ;
07593
07594 for( ib=0 ; ib < ts_length ; ib++ )
07595 tt[ib] = be->tbot + ib*dt ;
07596
07597 bb = (float **) malloc( sizeof(float *) * nf ) ;
07598 for( pp=0 ; pp < nf ; pp++ ){
07599 bb[pp] = (float *) malloc( sizeof(float) * ts_length ) ;
07600 for( ib=0 ; ib < ts_length ; ib++ )
07601 bb[pp][ib] = basis_funceval( be->bfunc[pp] , tt[ib] ) ;
07602 }
07603
07604
07605
07606
07607
07608 for( ii=0 ; ii < nvox ; ii++ ){
07609 allz = 1 ;
07610 for( pp=0 ; pp < nf ; pp++ ){
07611 wt[pp] = wtar[pp][ii] ; allz = ( allz && (wt[pp] == 0.0f) );
07612 }
07613
07614 if( allz ){
07615 for( ib=0 ; ib < ts_length ; ib++ ) hout[ib][ii] = 0.0f ;
07616 } else {
07617 register float sum ;
07618 for( ib=0 ; ib < ts_length ; ib++ ){
07619 sum = 0.0f ;
07620 for( pp=0 ; pp < nf ; pp++ ) sum += wt[pp] * bb[pp][ib] ;
07621 hout[ib][ii] = sum ;
07622 }
07623 }
07624 }
07625
07626
07627
07628 for( pp=0 ; pp < nf ; pp++ ) free((void *)bb[pp]) ;
07629 free((void *)bb) ; free((void *)tt) ; free((void *)wt) ;
07630
07631
07632
07633 for( ib=0 ; ib < ts_length ; ib++ ){
07634 bar = (short *) malloc( sizeof(short) * nvox ) ;
07635 factor = EDIT_coerce_autoscale_new(nvox, MRI_float,hout[ib], MRI_short,bar) ;
07636 if( factor < EPSILON ) factor = 0.0f ;
07637 else factor = 1.0f / factor ;
07638 EDIT_BRICK_FACTOR( out_dset , ib , factor ) ;
07639 EDIT_substitute_brick( out_dset , ib , MRI_short , bar ) ;
07640 free((void *)hout[ib]) ;
07641 }
07642 free((void *)hout) ;
07643
07644
07645
07646 DSET_write( out_dset ) ;
07647 if( verb )
07648 fprintf(stderr,"++ Wrote iresp 3D+time dataset into %s\n",DSET_BRIKNAME(out_dset)) ;
07649
07650 DSET_delete( out_dset ) ;
07651 return ;
07652 }
07653
07654
07655
07656 void basis_write_sresp( int argc , char *argv[] ,
07657 struct DC_options *option_data ,
07658 basis_expansion *be , float dt ,
07659 float *mse ,
07660 int pbot, matrix cvar, char *output_filename )
07661 {
07662 int nvox, ii, nf, allz, ts_length ;
07663 register int pp,qq, ib ;
07664 register float sum ;
07665 float *vv , *tt , **hout , factor , **bb ;
07666 short *bar ;
07667 THD_3dim_dataset *in_dset = NULL;
07668 THD_3dim_dataset *out_dset = NULL;
07669 char *commandline , label[512] ;
07670 const float EPSILON = 1.0e-10 ;
07671
07672
07673
07674 in_dset = THD_open_dataset(option_data->input_filename);
07675 nvox = in_dset->daxes->nxx * in_dset->daxes->nyy * in_dset->daxes->nzz;
07676
07677 if( dt <= 0.0f ) dt = DSET_TR(in_dset) ;
07678 if( dt <= 0.0f ) dt = 1.0f ;
07679
07680
07681
07682 out_dset = EDIT_empty_copy( in_dset ) ;
07683 tross_Copy_History( in_dset , out_dset ) ;
07684 DSET_delete( in_dset ) ;
07685
07686
07687
07688 commandline = tross_commandline( PROGRAM_NAME , argc , argv ) ;
07689 sprintf( label , "Sigma response: %s" , output_filename ) ;
07690 if( commandline != NULL )
07691 tross_multi_Append_History( out_dset , commandline,label,NULL ) ;
07692 else
07693 tross_Append_History ( out_dset, label);
07694 free((void *)commandline) ;
07695
07696 ts_length = 1 + (int)ceil( (be->ttop - be->tbot)/dt ) ;
07697
07698
07699
07700 (void ) EDIT_dset_items( out_dset,
07701 ADN_prefix, output_filename,
07702 ADN_label1, output_filename,
07703 ADN_self_name, output_filename,
07704 ADN_malloc_type, DATABLOCK_MEM_MALLOC,
07705 ADN_datum_all, MRI_short,
07706 ADN_nvals, ts_length,
07707 ADN_ntt, ts_length,
07708 ADN_ttdel, dt ,
07709 ADN_ttorg, be->tbot,
07710 ADN_none ) ;
07711
07712 if( THD_is_file(out_dset->dblk->diskptr->header_name) ){
07713 fprintf(stderr,
07714 "** ERROR: Output dataset file %s already exists - won't overwrite\n",
07715 out_dset->dblk->diskptr->header_name ) ;
07716 DSET_delete(out_dset) ; return ;
07717 }
07718
07719
07720
07721 hout = (float **) malloc( sizeof(float *) * ts_length ) ;
07722 for( ib=0 ; ib < ts_length ; ib++ )
07723 hout[ib] = (float *)calloc(sizeof(float),nvox) ;
07724
07725 nf = be->nfunc ;
07726 tt = (float *) malloc( sizeof(float) * ts_length ) ;
07727
07728 for( ib=0 ; ib < ts_length ; ib++ )
07729 tt[ib] = be->tbot + ib*dt ;
07730
07731
07732
07733 bb = (float ** ) malloc( sizeof(float *) * nf ) ;
07734 for( pp=0 ; pp < nf ; pp++ ){
07735 bb[pp] = (float * ) malloc( sizeof(float ) * ts_length ) ;
07736 for( ib=0 ; ib < ts_length ; ib++ )
07737 bb[pp][ib] = basis_funceval( be->bfunc[pp] , tt[ib] ) ;
07738 }
07739 free((void *)tt) ;
07740
07741
07742
07743 vv = (float *) malloc( sizeof(float) * ts_length ) ;
07744 for( ib=0 ; ib < ts_length ; ib++ ){
07745 sum = 0.0f ;
07746 for( pp=0 ; pp < nf ; pp++ ){
07747 for( qq=0 ; qq < nf ; qq++ )
07748 sum += cvar.elts[pbot+pp][pbot+qq] * bb[pp][ib] * bb[qq][ib] ;
07749 }
07750 vv[ib] = (sum >= 0.0f) ? sum : 0.0f ;
07751 }
07752 for( pp=0 ; pp < nf ; pp++ ) free((void *)bb[pp]) ;
07753 free((void *)bb) ;
07754
07755
07756
07757 for( ii=0 ; ii < nvox ; ii++ ){
07758 for( ib=0 ; ib < ts_length ; ib++ )
07759 hout[ib][ii] = sqrt( vv[ib] * mse[ii] ) ;
07760 }
07761 free((void *)vv) ;
07762
07763
07764
07765 for( ib=0 ; ib < ts_length ; ib++ ){
07766 bar = (short *) malloc( sizeof(short) * nvox ) ;
07767 factor = EDIT_coerce_autoscale_new(nvox, MRI_float,hout[ib], MRI_short,bar) ;
07768 if( factor < EPSILON ) factor = 0.0f ;
07769 else factor = 1.0f / factor ;
07770 EDIT_BRICK_FACTOR( out_dset , ib , factor ) ;
07771 EDIT_substitute_brick( out_dset , ib , MRI_short , bar ) ;
07772 free((void *)hout[ib]) ;
07773 }
07774 free((void *)hout) ;
07775
07776
07777
07778 DSET_write( out_dset ) ;
07779 if( verb )
07780 fprintf(stderr,"++ Wrote sresp 3D+time dataset into %s\n",DSET_BRIKNAME(out_dset)) ;
07781
07782 DSET_delete( out_dset ) ;
07783 return ;
07784 }
07785
07786
07787
07788 float baseline_mean( vector coef )
07789 {
07790 register int jj ;
07791 register double sum ;
07792
07793 sum = 0.0 ;
07794 for( jj=0 ; jj < nParam ; jj++ )
07795 if( Xcol_inbase[jj] == 2 ) sum += coef.elts[jj] * Xcol_mean[jj] ;
07796
07797 return (float)sum ;
07798 }
07799
07800
07801
07802 #if 0
07803 typedef struct {
07804 int npar , pbot ;
07805 float *ww ;
07806 float scale_fac ;
07807 int denom_flag ;
07808 char *name ;
07809 } basis_irc ;
07810 #endif
07811
07812 floatpair evaluate_irc( basis_irc *birc , vector coef ,
07813 float base , float mse , matrix cvar )
07814 {
07815 floatpair vt={0.0f,0.0f} ;
07816 int ii,jj, np, pb ;
07817 double asum , bsum ;
07818 float *ww ;
07819
07820 np = birc->npar ;
07821 pb = birc->pbot ;
07822 ww = birc->ww ;
07823
07824 asum = 0.0 ;
07825 for( ii=0 ; ii < np ; ii++ )
07826 asum += coef.elts[pb+ii] * ww[ii] ;
07827
07828 bsum = 0.0 ;
07829 for( ii=0 ; ii < np ; ii++ )
07830 for( jj=0 ; jj < np ; jj++ )
07831 bsum += cvar.elts[pb+ii][pb+jj] * ww[ii] * ww[jj] ;
07832
07833 bsum *= mse ;
07834
07835 if( bsum > 0.0 ) vt.b = asum / sqrt(bsum) ;
07836
07837 vt.a = (float)(asum * birc->scale_fac) ;
07838 if( birc->denom_flag && denom_BASELINE ){
07839 if( base == 0.0f ) vt.a = 0.0f ;
07840 else vt.a /= base ;
07841 }
07842 return vt ;
07843 }