00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #define PROGRAM_NAME "3dWavelets"
00028 #define PROGRAM_AUTHOR "B. Douglas Ward"
00029 #define PROGRAM_INITIAL "28 March 2000"
00030 #define PROGRAM_LATEST "02 December 2002"
00031
00032
00033
00034
00035
00036
00037 #include "mrilib.h"
00038 #include "Wavelets.h"
00039 #include "Wavelets.c"
00040
00041
00042
00043
00044
00045
00046
00047 #define MAX_NAME_LENGTH THD_MAX_NAME
00048 #define MAX_FILTERS 20
00049
00050
00051
00052
00053
00054
00055
00056 typedef struct WA_options
00057 {
00058 int NFirst;
00059 int NLast;
00060 int N;
00061 int f;
00062 int p;
00063 int q;
00064
00065 float fdisp;
00066 char * input_filename;
00067 char * mask_filename;
00068 char * input1D_filename;
00069
00070 int wavelet_type;
00071
00072 int num_stop_filters;
00073 int stop_band[MAX_FILTERS];
00074 int stop_mintr[MAX_FILTERS];
00075 int stop_maxtr[MAX_FILTERS];
00076 float * stop_filter;
00077
00078 int num_base_filters;
00079 int base_band[MAX_FILTERS];
00080 int base_mintr[MAX_FILTERS];
00081 int base_maxtr[MAX_FILTERS];
00082 float * base_filter;
00083
00084 int num_sgnl_filters;
00085 int sgnl_band[MAX_FILTERS];
00086 int sgnl_mintr[MAX_FILTERS];
00087 int sgnl_maxtr[MAX_FILTERS];
00088 float * sgnl_filter;
00089
00090 char * bucket_filename;
00091 char * coefts_filename;
00092 char * fitts_filename;
00093 char * sgnlts_filename;
00094 char * errts_filename;
00095
00096 int fout;
00097 int rout;
00098 int cout;
00099 int vout;
00100 int stat_first;
00101
00102 } WA_options;
00103
00104
00105
00106
00107
00108
00109
00110 void display_help_menu()
00111 {
00112 printf (
00113 "Program to perform wavelet analysis of an FMRI 3d+time dataset. \n"
00114 " \n"
00115 "Usage: \n"
00116 "3dWavelets \n"
00117 "-type wname wname = name of wavelet to use for the analysis \n"
00118 " At present, there are only two choices for wname: \n"
00119 " Haar --> Haar wavelets \n"
00120 " Daub --> Daubechies wavelets \n"
00121 "-input fname fname = filename of 3d+time input dataset \n"
00122 "[-input1D dname] dname = filename of single (fMRI) .1D time series \n"
00123 "[-mask mname] mname = filename of 3d mask dataset \n"
00124 "[-nfirst fnum] fnum = number of first dataset image to use in \n"
00125 " the wavelet analysis. (default = 0) \n"
00126 "[-nlast lnum] lnum = number of last dataset image to use in \n"
00127 " the wavelet analysis. (default = last) \n"
00128 "[-fdisp fval] Write (to screen) results for those voxels \n"
00129 " whose F-statistic is >= fval \n"
00130 " \n"
00131 "Filter options: \n"
00132 "[-filt_stop band mintr maxtr] Specify wavelet coefs. to set to zero \n"
00133 "[-filt_base band mintr maxtr] Specify wavelet coefs. for baseline model\n"
00134 "[-filt_sgnl band mintr maxtr] Specify wavelet coefs. for signal model \n"
00135 " where band = frequency band \n"
00136 " mintr = min. value for time window (in TR) \n"
00137 " maxtr = max. value for time window (in TR) \n"
00138 " \n"
00139 "Output options: \n"
00140 "[-coefts cprefix] cprefix = prefix of 3d+time output dataset which \n"
00141 " will contain the forward wavelet transform \n"
00142 " coefficients \n"
00143 " \n"
00144 "[-fitts fprefix] fprefix = prefix of 3d+time output dataset which \n"
00145 " will contain the full model time series fit \n"
00146 " to the input data \n"
00147 " \n"
00148 "[-sgnlts sprefix] sprefix = prefix of 3d+time output dataset which \n"
00149 " will contain the signal model time series fit \n"
00150 " to the input data \n"
00151 " \n"
00152 "[-errts eprefix] eprefix = prefix of 3d+time output dataset which \n"
00153 " will contain the residual error time series \n"
00154 " from the full model fit to the input data \n"
00155 " \n"
00156 "The following options control the contents of the bucket dataset: \n"
00157 "[-fout] Flag to output the F-statistics \n"
00158 "[-rout] Flag to output the R^2 statistics \n"
00159 "[-cout] Flag to output the full model wavelet coefficients \n"
00160 "[-vout] Flag to output the sample variance (MSE) map \n"
00161 " \n"
00162 "[-stat_first] Flag to specify that the full model statistics will \n"
00163 " appear prior to the wavelet coefficients in the \n"
00164 " bucket dataset output \n"
00165 " \n"
00166 "[-bucket bprefix] bprefix = prefix of AFNI 'bucket' dataset containing\n"
00167 " parameters of interest, such as the F-statistic \n"
00168 " for significance of the wavelet signal model. \n"
00169 );
00170
00171 exit(0);
00172 }
00173
00174
00175
00176
00177
00178
00179
00180 void initialize_options
00181 (
00182 WA_options * option_data
00183 )
00184
00185 {
00186 int is;
00187
00188
00189
00190 option_data->NFirst = 0;
00191 option_data->NLast = 32767;
00192 option_data->N = 0;
00193 option_data->f = 0;
00194 option_data->p = 0;
00195 option_data->q = 0;
00196 option_data->fdisp = -1.0;
00197 option_data->wavelet_type = 0;
00198
00199
00200
00201 option_data->num_stop_filters = 0;
00202 option_data->num_base_filters = 0;
00203 option_data->num_sgnl_filters = 0;
00204 for (is = 0; is < MAX_FILTERS; is++)
00205 {
00206 option_data->stop_band[is] = 0;
00207 option_data->stop_mintr[is] = 0;
00208 option_data->stop_maxtr[is] = 0;
00209 option_data->base_band[is] = 0;
00210 option_data->base_mintr[is] = 0;
00211 option_data->base_maxtr[is] = 0;
00212 option_data->sgnl_band[is] = 0;
00213 option_data->sgnl_mintr[is] = 0;
00214 option_data->sgnl_maxtr[is] = 0;
00215 }
00216 option_data->stop_filter = NULL;
00217 option_data->base_filter = NULL;
00218 option_data->sgnl_filter = NULL;
00219
00220
00221
00222 option_data->fout = 0;
00223 option_data->rout = 0;
00224 option_data->cout = 0;
00225 option_data->vout = 0;
00226 option_data->stat_first = 0;
00227
00228
00229
00230 option_data->input_filename = NULL;
00231 option_data->mask_filename = NULL;
00232 option_data->input1D_filename = NULL;
00233 option_data->bucket_filename = NULL;
00234 option_data->coefts_filename = NULL;
00235 option_data->fitts_filename = NULL;
00236 option_data->sgnlts_filename = NULL;
00237 option_data->errts_filename = NULL;
00238
00239 }
00240
00241
00242
00243
00244
00245
00246
00247 void get_options
00248 (
00249 int argc,
00250 char ** argv,
00251 WA_options * option_data
00252 )
00253
00254 {
00255 int nopt = 1;
00256 int ival, index;
00257 float fval;
00258 char message[MAX_NAME_LENGTH];
00259 int ifilter;
00260
00261
00262
00263 machdep() ;
00264 { int new_argc ; char ** new_argv ;
00265 addto_args( argc , argv , &new_argc , &new_argv ) ;
00266 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00267 }
00268
00269
00270
00271 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
00272
00273
00274
00275 AFNI_logger (PROGRAM_NAME,argc,argv);
00276
00277
00278
00279 initialize_options (option_data);
00280
00281
00282
00283 while (nopt < argc )
00284 {
00285
00286
00287 if (strncmp(argv[nopt], "-type", 5) == 0)
00288 {
00289 nopt++;
00290 if (nopt >= argc) WA_error ("need argument after -type ");
00291 for (ival = 0; ival < MAX_WAVELET_TYPE; ival++)
00292 if (strncmp(argv[nopt], WAVELET_TYPE_name[ival], 4) == 0) break;
00293 if (ival >= MAX_WAVELET_TYPE)
00294 {
00295 sprintf(message,"Unrecognized wavelet type: %s\n", argv[nopt]);
00296 WA_error (message);
00297 }
00298 option_data->wavelet_type = ival;
00299 nopt++;
00300 continue;
00301 }
00302
00303
00304
00305 if (strncmp(argv[nopt], "-input", 8) == 0)
00306 {
00307 nopt++;
00308 if (nopt >= argc) WA_error ("need argument after -input ");
00309 option_data->input_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00310 MTEST (option_data->input_filename);
00311 strcpy (option_data->input_filename, argv[nopt]);
00312 nopt++;
00313 continue;
00314 }
00315
00316
00317
00318 if (strncmp(argv[nopt], "-input1D", 8) == 0)
00319 {
00320 nopt++;
00321 if (nopt >= argc) WA_error ("need argument after -input1D ");
00322 option_data->input1D_filename =
00323 malloc (sizeof(char)*MAX_NAME_LENGTH);
00324 MTEST (option_data->input1D_filename);
00325 strcpy (option_data->input1D_filename, argv[nopt]);
00326 nopt++;
00327 continue;
00328 }
00329
00330
00331
00332 if (strncmp(argv[nopt], "-mask", 5) == 0)
00333 {
00334 nopt++;
00335 if (nopt >= argc) WA_error ("need argument after -mask ");
00336 option_data->mask_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00337 MTEST (option_data->mask_filename);
00338 strcpy (option_data->mask_filename, argv[nopt]);
00339 nopt++;
00340 continue;
00341 }
00342
00343
00344
00345 if (strncmp(argv[nopt], "-nfirst", 7) == 0)
00346 {
00347 nopt++;
00348 if (nopt >= argc) WA_error ("need argument after -nfirst ");
00349 sscanf (argv[nopt], "%d", &ival);
00350 if (ival < 0)
00351 WA_error ("illegal argument after -nfirst ");
00352 option_data->NFirst = ival;
00353 nopt++;
00354 continue;
00355 }
00356
00357
00358
00359 if (strncmp(argv[nopt], "-nlast", 6) == 0)
00360 {
00361 nopt++;
00362 if (nopt >= argc) WA_error ("need argument after -nlast ");
00363 sscanf (argv[nopt], "%d", &ival);
00364 if (ival < 0)
00365 WA_error ("illegal argument after -nlast ");
00366 option_data->NLast = ival;
00367 nopt++;
00368 continue;
00369 }
00370
00371
00372
00373 if (strncmp(argv[nopt], "-fdisp", 6) == 0)
00374 {
00375 nopt++;
00376 if (nopt >= argc) WA_error ("need argument after -fdisp ");
00377 sscanf (argv[nopt], "%f", &fval);
00378 option_data->fdisp = fval;
00379 nopt++;
00380 continue;
00381 }
00382
00383
00384
00385 if (strncmp(argv[nopt], "-filt_stop", 10) == 0)
00386 {
00387 nopt++;
00388 if (nopt+2 >= argc)
00389 WA_error ("need 3 arguments after -filt_stop");
00390
00391 ifilter = option_data->num_stop_filters;
00392
00393 sscanf (argv[nopt], "%d", &ival);
00394 if (ival < -1)
00395 WA_error ("-filt_stop band mintr maxtr Require: -1 <= band");
00396 option_data->stop_band[ifilter] = ival;
00397 nopt++;
00398
00399 sscanf (argv[nopt], "%d", &ival);
00400 if (ival < 0)
00401 WA_error ("-filt_stop band mintr maxtr Require: 0 <= mintr");
00402 option_data->stop_mintr[ifilter] = ival;
00403 nopt++;
00404
00405 sscanf (argv[nopt], "%d", &ival);
00406 if (ival < 0)
00407 WA_error ("-filt_stop band mintr maxtr Require: 0 <= maxtr");
00408 option_data->stop_maxtr[ifilter] = ival;
00409 nopt++;
00410
00411 option_data->num_stop_filters++;
00412 continue;
00413 }
00414
00415
00416
00417 if (strncmp(argv[nopt], "-filt_base", 10) == 0)
00418 {
00419 nopt++;
00420 if (nopt+2 >= argc)
00421 WA_error ("need 3 arguments after -filt_base");
00422
00423 ifilter = option_data->num_base_filters;
00424
00425 sscanf (argv[nopt], "%d", &ival);
00426 if (ival < -1)
00427 WA_error ("-filt_base band mintr maxtr Require: -1 <= band");
00428 option_data->base_band[ifilter] = ival;
00429 nopt++;
00430
00431 sscanf (argv[nopt], "%d", &ival);
00432 if (ival < 0)
00433 WA_error ("-filt_base band mintr maxtr Require: 0 <= mintr");
00434 option_data->base_mintr[ifilter] = ival;
00435 nopt++;
00436
00437 sscanf (argv[nopt], "%d", &ival);
00438 if (ival < 0)
00439 WA_error ("-filt_base band mintr maxtr Require: 0 <= maxtr");
00440 option_data->base_maxtr[ifilter] = ival;
00441 nopt++;
00442
00443 option_data->num_base_filters++;
00444 continue;
00445 }
00446
00447
00448
00449 if (strncmp(argv[nopt], "-filt_sgnl", 10) == 0)
00450 {
00451 nopt++;
00452 if (nopt+2 >= argc)
00453 WA_error ("need 3 arguments after -filt_sgnl");
00454
00455 ifilter = option_data->num_sgnl_filters;
00456
00457 sscanf (argv[nopt], "%d", &ival);
00458 if (ival < -1)
00459 WA_error ("-filt_sgnl band mintr maxtr Require: -1 <= band");
00460 option_data->sgnl_band[ifilter] = ival;
00461 nopt++;
00462
00463 sscanf (argv[nopt], "%d", &ival);
00464 if (ival < 0)
00465 WA_error ("-filt_sgnl band mintr maxtr Require: 0 <= mintr");
00466 option_data->sgnl_mintr[ifilter] = ival;
00467 nopt++;
00468
00469 sscanf (argv[nopt], "%d", &ival);
00470 if (ival < 0)
00471 WA_error ("-filt_sgnl band mintr maxtr Require: 0 <= maxtr");
00472 option_data->sgnl_maxtr[ifilter] = ival;
00473 nopt++;
00474
00475 option_data->num_sgnl_filters++;
00476 continue;
00477 }
00478
00479
00480
00481 if (strncmp(argv[nopt], "-fout", 5) == 0)
00482 {
00483 option_data->fout = 1;
00484 nopt++;
00485 continue;
00486 }
00487
00488
00489
00490 if (strncmp(argv[nopt], "-rout", 5) == 0)
00491 {
00492 option_data->rout = 1;
00493 nopt++;
00494 continue;
00495 }
00496
00497
00498
00499 if (strncmp(argv[nopt], "-cout", 5) == 0)
00500 {
00501 option_data->cout = 1;
00502 nopt++;
00503 continue;
00504 }
00505
00506
00507
00508 if (strncmp(argv[nopt], "-vout", 5) == 0)
00509 {
00510 option_data->vout = 1;
00511 nopt++;
00512 continue;
00513 }
00514
00515
00516
00517 if (strncmp(argv[nopt], "-stat_first", 11) == 0)
00518 {
00519 option_data->stat_first = 1;
00520 nopt++;
00521 continue;
00522 }
00523
00524
00525
00526 if (strncmp(argv[nopt], "-bucket", 6) == 0)
00527 {
00528 nopt++;
00529 if (nopt >= argc) WA_error ("need file prefixname after -bucket ");
00530 option_data->bucket_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00531 MTEST (option_data->bucket_filename);
00532 strcpy (option_data->bucket_filename, argv[nopt]);
00533 nopt++;
00534 continue;
00535 }
00536
00537
00538
00539 if (strncmp(argv[nopt], "-coefts", 8) == 0)
00540 {
00541 nopt++;
00542 if (nopt >= argc) WA_error ("need file prefixname after -coefts ");
00543 option_data->coefts_filename =
00544 malloc (sizeof(char)*MAX_NAME_LENGTH);
00545 MTEST (option_data->coefts_filename);
00546 strcpy (option_data->coefts_filename, argv[nopt]);
00547 nopt++;
00548 continue;
00549 }
00550
00551
00552
00553 if (strncmp(argv[nopt], "-fitts", 7) == 0)
00554 {
00555 nopt++;
00556 if (nopt >= argc) WA_error ("need file prefixname after -fitts ");
00557 option_data->fitts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00558 MTEST (option_data->fitts_filename);
00559 strcpy (option_data->fitts_filename, argv[nopt]);
00560 nopt++;
00561 continue;
00562 }
00563
00564
00565
00566 if (strncmp(argv[nopt], "-sgnlts", 7) == 0)
00567 {
00568 nopt++;
00569 if (nopt >= argc) WA_error ("need file prefixname after -sgnlts ");
00570 option_data->sgnlts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00571 MTEST (option_data->sgnlts_filename);
00572 strcpy (option_data->sgnlts_filename, argv[nopt]);
00573 nopt++;
00574 continue;
00575 }
00576
00577
00578
00579 if (strncmp(argv[nopt], "-errts", 6) == 0)
00580 {
00581 nopt++;
00582 if (nopt >= argc) WA_error ("need file prefixname after -errts ");
00583 option_data->errts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00584 MTEST (option_data->errts_filename);
00585 strcpy (option_data->errts_filename, argv[nopt]);
00586 nopt++;
00587 continue;
00588 }
00589
00590
00591
00592 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
00593 WA_error (message);
00594
00595 }
00596
00597
00598 }
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608 float * read_time_series
00609 (
00610 char * ts_filename,
00611 int * ts_length
00612 )
00613
00614 {
00615 char message[MAX_NAME_LENGTH];
00616 char * cpt;
00617 char filename[THD_MAX_NAME];
00618 char subv[THD_MAX_NAME];
00619 MRI_IMAGE * im, * flim;
00620
00621 float * far;
00622 int nx;
00623 int ny;
00624 int iy;
00625 int ipt;
00626 float * ts_data = NULL;
00627
00628
00629
00630 flim = mri_read_1D(ts_filename) ;
00631 if (flim == NULL)
00632 {
00633 sprintf (message, "Unable to read time series file: %s", ts_filename);
00634 WA_error (message);
00635 }
00636
00637
00638
00639 far = MRI_FLOAT_PTR(flim);
00640 nx = flim->nx;
00641 ny = flim->ny; iy = 0 ;
00642 if( ny > 1 ){
00643 fprintf(stderr,"WARNING: time series %s has %d columns\n",ts_filename,ny);
00644 }
00645
00646
00647
00648 *ts_length = nx;
00649 ts_data = (float *) malloc (sizeof(float) * nx);
00650 MTEST (ts_data);
00651 for (ipt = 0; ipt < nx; ipt++)
00652 ts_data[ipt] = far[ipt + iy*nx];
00653
00654
00655 mri_free (flim); flim = NULL;
00656
00657 return (ts_data);
00658 }
00659
00660
00661
00662
00663
00664
00665
00666 void read_input_data
00667 (
00668 WA_options * option_data,
00669 THD_3dim_dataset ** dset_time,
00670 THD_3dim_dataset ** mask_dset,
00671 float ** fmri_data,
00672 int * fmri_length
00673 )
00674
00675 {
00676 char message[MAX_NAME_LENGTH];
00677
00678
00679
00680 if (option_data->input1D_filename != NULL)
00681 {
00682
00683 *fmri_data = read_time_series (option_data->input1D_filename,
00684 fmri_length);
00685 if (*fmri_data == NULL)
00686 {
00687 sprintf (message, "Unable to read time series file: %s",
00688 option_data->input1D_filename);
00689 WA_error (message);
00690 }
00691 *dset_time = NULL;
00692 }
00693
00694 else if (option_data->input_filename != NULL)
00695 {
00696
00697 *dset_time = THD_open_one_dataset (option_data->input_filename);
00698 if (!ISVALID_3DIM_DATASET(*dset_time))
00699 {
00700 sprintf (message, "Unable to open data file: %s",
00701 option_data->input_filename);
00702 WA_error (message);
00703 }
00704 THD_load_datablock ((*dset_time)->dblk);
00705
00706 if (option_data->mask_filename != NULL)
00707 {
00708
00709 *mask_dset = THD_open_dataset (option_data->mask_filename);
00710 if (!ISVALID_3DIM_DATASET(*mask_dset))
00711 {
00712 sprintf (message, "Unable to open mask file: %s",
00713 option_data->mask_filename);
00714 WA_error (message);
00715 }
00716 THD_load_datablock ((*mask_dset)->dblk);
00717 }
00718 }
00719 else
00720 WA_error ("Must specify input data");
00721
00722 }
00723
00724
00725
00726
00727
00728
00729
00730 void check_one_output_file
00731 (
00732 THD_3dim_dataset * dset_time,
00733 char * filename
00734 )
00735
00736 {
00737 char message[MAX_NAME_LENGTH];
00738 THD_3dim_dataset * new_dset=NULL;
00739 int ierror;
00740
00741
00742
00743 new_dset = EDIT_empty_copy( dset_time ) ;
00744
00745
00746 ierror = EDIT_dset_items( new_dset ,
00747 ADN_prefix , filename ,
00748 ADN_label1 , filename ,
00749 ADN_self_name , filename ,
00750 ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE :
00751 GEN_FUNC_TYPE ,
00752 ADN_none ) ;
00753
00754 if( ierror > 0 )
00755 {
00756 sprintf (message,
00757 "*** %d errors in attempting to create output dataset!\n",
00758 ierror);
00759 WA_error (message);
00760 }
00761
00762 if( THD_is_file(new_dset->dblk->diskptr->header_name) )
00763 {
00764 sprintf (message,
00765 "Output dataset file %s already exists "
00766 " -- cannot continue!\a\n",
00767 new_dset->dblk->diskptr->header_name);
00768 WA_error (message);
00769 }
00770
00771
00772 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00773
00774 }
00775
00776
00777
00778
00779
00780
00781
00782 void check_output_files
00783 (
00784 WA_options * option_data,
00785 THD_3dim_dataset * dset_time
00786 )
00787
00788 {
00789
00790 if (option_data->bucket_filename != NULL)
00791 check_one_output_file (dset_time, option_data->bucket_filename);
00792
00793 if (option_data->coefts_filename != NULL)
00794 check_one_output_file (dset_time, option_data->coefts_filename);
00795
00796 if (option_data->fitts_filename != NULL)
00797 check_one_output_file (dset_time, option_data->fitts_filename);
00798
00799 if (option_data->sgnlts_filename != NULL)
00800 check_one_output_file (dset_time, option_data->sgnlts_filename);
00801
00802 if (option_data->errts_filename != NULL)
00803 check_one_output_file (dset_time, option_data->errts_filename);
00804
00805 }
00806
00807
00808
00809
00810
00811
00812
00813 void check_for_valid_inputs
00814 (
00815 WA_options * option_data,
00816 THD_3dim_dataset * dset_time,
00817 THD_3dim_dataset * mask_dset
00818 )
00819
00820 {
00821 char message[MAX_NAME_LENGTH];
00822
00823
00824
00825 if (mask_dset != NULL)
00826 {
00827 if ( (DSET_NX(dset_time) != DSET_NX(mask_dset))
00828 || (DSET_NY(dset_time) != DSET_NY(mask_dset))
00829 || (DSET_NZ(dset_time) != DSET_NZ(mask_dset)) )
00830 {
00831 sprintf (message, "%s and %s have incompatible dimensions",
00832 option_data->input_filename, option_data->mask_filename);
00833 WA_error (message);
00834 }
00835
00836 if (DSET_NVALS(mask_dset) != 1 )
00837 WA_error ("Must specify 1 sub-brick from mask dataset");
00838 }
00839
00840
00841
00842 check_output_files (option_data, dset_time);
00843
00844 }
00845
00846
00847
00848
00849
00850
00851
00852 void initialize_filters
00853 (
00854 WA_options * option_data,
00855 THD_3dim_dataset * dset_time,
00856 THD_3dim_dataset * mask_dset,
00857 int fmri_length
00858 )
00859
00860 {
00861 char message[MAX_NAME_LENGTH];
00862 int nt;
00863 int NFirst;
00864 int NLast;
00865 int N;
00866 int ifilter;
00867 int i;
00868 int f, q, p;
00869
00870
00871
00872 if (option_data->input1D_filename != NULL)
00873 nt = fmri_length;
00874 else
00875 nt = DSET_NUM_TIMES (dset_time);
00876
00877
00878
00879 NFirst = option_data->NFirst;
00880
00881 NLast = option_data->NLast;
00882 if (NLast > nt-1) NLast = nt-1;
00883
00884 N = NLast - NFirst + 1;
00885 N = powerof2(my_log2(N));
00886 NLast = N + NFirst + 1;
00887
00888 option_data->N = N;
00889 option_data->NLast = NLast;
00890
00891
00892
00893 option_data->stop_filter =
00894 FWT_1d_stop_filter (option_data->num_stop_filters, option_data->stop_band,
00895 option_data->stop_mintr, option_data->stop_maxtr,
00896 option_data->NFirst, option_data->N);
00897
00898 option_data->base_filter =
00899 FWT_1d_pass_filter (option_data->num_base_filters, option_data->base_band,
00900 option_data->base_mintr, option_data->base_maxtr,
00901 option_data->NFirst, option_data->N);
00902
00903 option_data->sgnl_filter =
00904 FWT_1d_pass_filter (option_data->num_sgnl_filters, option_data->sgnl_band,
00905 option_data->sgnl_mintr, option_data->sgnl_maxtr,
00906 option_data->NFirst, option_data->N);
00907
00908
00909
00910 f = 0;
00911 for (i = 0; i < N; i++)
00912 if (option_data->stop_filter[i] == 0.0)
00913 {
00914 f++;
00915 option_data->base_filter[i] = 0.0;
00916 option_data->sgnl_filter[i] = 0.0;
00917 }
00918 option_data->f = f;
00919
00920 q = 0;
00921 for (i = 0; i < N; i++)
00922 if (option_data->base_filter[i] == 1.0)
00923 {
00924 q++;
00925 option_data->sgnl_filter[i] = 1.0;
00926 }
00927 option_data->q = q;
00928
00929 p = 0;
00930 for (i = 0; i < N; i++)
00931 if (option_data->sgnl_filter[i] == 1.0)
00932 {
00933 p++;
00934 }
00935 option_data->p = p;
00936
00937
00938 return;
00939 }
00940
00941
00942
00943
00944
00945
00946
00947 void allocate_memory
00948 (
00949 WA_options * option_data,
00950 THD_3dim_dataset * dset,
00951
00952 float *** coef_vol,
00953 float ** mse_vol,
00954 float ** ffull_vol,
00955 float ** rfull_vol,
00956
00957 float *** coefts_vol,
00958 float *** fitts_vol,
00959 float *** sgnlts_vol,
00960 float *** errts_vol
00961 )
00962
00963 {
00964 int nxyz;
00965 int ixyz;
00966 int it;
00967 int N;
00968 int ip;
00969 int p;
00970
00971
00972
00973 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
00974 N = option_data->N;
00975 p = option_data->p;
00976
00977
00978
00979 *coef_vol = (float **) malloc (sizeof(float *) * p); MTEST(*coef_vol);
00980
00981 for (ip = 0; ip < p; ip++)
00982 {
00983 (*coef_vol)[ip] = (float *) malloc (sizeof(float) * nxyz);
00984 MTEST((*coef_vol)[ip]);
00985 for (ixyz = 0; ixyz < nxyz; ixyz++)
00986 {
00987 (*coef_vol)[ip][ixyz] = 0.0;
00988 }
00989 }
00990
00991
00992 *mse_vol = (float *) malloc (sizeof(float) * nxyz); MTEST(*mse_vol);
00993 *ffull_vol = (float *) malloc (sizeof(float) * nxyz); MTEST(*ffull_vol);
00994 *rfull_vol = (float *) malloc (sizeof(float) * nxyz); MTEST(*rfull_vol);
00995 for (ixyz = 0; ixyz < nxyz; ixyz++)
00996 {
00997 (*mse_vol)[ixyz] = 0.0;
00998 (*ffull_vol)[ixyz] = 0.0;
00999 (*rfull_vol)[ixyz] = 0.0;
01000 }
01001
01002
01003
01004 if (option_data->coefts_filename != NULL)
01005 {
01006 *coefts_vol = (float **) malloc (sizeof(float **) * N);
01007 MTEST (*coefts_vol);
01008 for (it = 0; it < N; it++)
01009 {
01010 (*coefts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz);
01011 MTEST ((*coefts_vol)[it]);
01012 for (ixyz = 0; ixyz < nxyz; ixyz++)
01013 (*coefts_vol)[it][ixyz] = 0.0;
01014 }
01015 }
01016
01017
01018 if (option_data->fitts_filename != NULL)
01019 {
01020 *fitts_vol = (float **) malloc (sizeof(float **) * N);
01021 MTEST (*fitts_vol);
01022 for (it = 0; it < N; it++)
01023 {
01024 (*fitts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz);
01025 MTEST ((*fitts_vol)[it]);
01026 for (ixyz = 0; ixyz < nxyz; ixyz++)
01027 (*fitts_vol)[it][ixyz] = 0.0;
01028 }
01029 }
01030
01031
01032 if (option_data->sgnlts_filename != NULL)
01033 {
01034 *sgnlts_vol = (float **) malloc (sizeof(float **) * N);
01035 MTEST (*sgnlts_vol);
01036 for (it = 0; it < N; it++)
01037 {
01038 (*sgnlts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz);
01039 MTEST ((*sgnlts_vol)[it]);
01040 for (ixyz = 0; ixyz < nxyz; ixyz++)
01041 (*sgnlts_vol)[it][ixyz] = 0.0;
01042 }
01043 }
01044
01045
01046 if (option_data->errts_filename != NULL)
01047 {
01048 *errts_vol = (float **) malloc (sizeof(float **) * N);
01049 MTEST (*errts_vol);
01050 for (it = 0; it < N; it++)
01051 {
01052 (*errts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz);
01053 MTEST ((*errts_vol)[it]);
01054 for (ixyz = 0; ixyz < nxyz; ixyz++)
01055 (*errts_vol)[it][ixyz] = 0.0;
01056 }
01057 }
01058 }
01059
01060
01061
01062
01063
01064
01065
01066 void initialize_program
01067 (
01068 int argc,
01069 char ** argv,
01070 WA_options ** option_data,
01071 THD_3dim_dataset ** dset_time,
01072 THD_3dim_dataset ** mask_dset,
01073 float ** fmri_data,
01074 int * fmri_length,
01075
01076 float *** coef_vol,
01077 float ** mse_vol,
01078 float ** ffull_vol,
01079 float ** rfull_vol,
01080
01081 float *** coefts_vol,
01082 float *** fitts_vol,
01083 float *** sgnlts_vol,
01084 float *** errts_vol
01085 )
01086
01087 {
01088
01089
01090 *option_data = (WA_options *) malloc (sizeof(WA_options));
01091
01092
01093
01094 get_options (argc, argv, *option_data);
01095
01096
01097
01098 read_input_data (*option_data, dset_time, mask_dset, fmri_data, fmri_length);
01099
01100
01101
01102 check_for_valid_inputs (*option_data, *dset_time, *mask_dset);
01103
01104
01105
01106 initialize_filters (*option_data, *dset_time, *mask_dset, *fmri_length);
01107
01108
01109
01110 if ((*option_data)->input1D_filename == NULL)
01111 allocate_memory (*option_data, *dset_time,
01112 coef_vol, mse_vol, ffull_vol, rfull_vol,
01113 coefts_vol, fitts_vol, sgnlts_vol, errts_vol);
01114
01115 }
01116
01117
01118
01119
01120
01121
01122
01123 void extract_ts_array
01124 (
01125 THD_3dim_dataset * dset_time,
01126 int iv,
01127 float * ts_array
01128 )
01129
01130 {
01131 MRI_IMAGE * im = NULL;
01132 float * ar = NULL;
01133 int ts_length;
01134 int it;
01135
01136
01137
01138 im = THD_extract_series (iv, dset_time, 0);
01139
01140
01141
01142 if (im == NULL) WA_error ("Unable to extract data from 3d+time dataset");
01143
01144
01145
01146 ts_length = DSET_NUM_TIMES (dset_time);
01147 ar = MRI_FLOAT_PTR (im);
01148 for (it = 0; it < ts_length; it++)
01149 {
01150 ts_array[it] = ar[it];
01151 }
01152
01153
01154
01155 mri_free (im); im = NULL;
01156
01157 }
01158
01159
01160
01161
01162
01163
01164
01165 void save_voxel
01166 (
01167 WA_options * option_data,
01168 int iv,
01169 float * coef,
01170 float mse,
01171 float ffull,
01172 float rfull,
01173
01174 float * coefts,
01175 float * fitts,
01176 float * sgnlts,
01177 float * errts,
01178
01179 float ** coef_vol,
01180 float * mse_vol,
01181 float * ffull_vol,
01182 float * rfull_vol,
01183
01184 float ** coefts_vol,
01185 float ** fitts_vol,
01186 float ** sgnlts_vol,
01187 float ** errts_vol
01188
01189 )
01190
01191 {
01192 int ip;
01193 int p;
01194 int it;
01195 int N;
01196
01197
01198
01199 p = option_data->p;
01200 N = option_data->N;
01201
01202
01203
01204 for (ip = 0; ip < p; ip++)
01205 {
01206 coef_vol[ip][iv] = coef[ip];
01207 }
01208
01209
01210
01211 mse_vol[iv] = mse;
01212
01213
01214
01215 ffull_vol[iv] = ffull;
01216
01217
01218
01219 rfull_vol[iv] = rfull;
01220
01221
01222
01223 if (coefts_vol != NULL)
01224 for (it = 0; it < N; it++)
01225 coefts_vol[it][iv] = coefts[it];
01226
01227
01228
01229 if (fitts_vol != NULL)
01230 for (it = 0; it < N; it++)
01231 fitts_vol[it][iv] = fitts[it];
01232
01233
01234
01235 if (sgnlts_vol != NULL)
01236 for (it = 0; it < N; it++)
01237 sgnlts_vol[it][iv] = sgnlts[it];
01238
01239
01240
01241 if (errts_vol != NULL)
01242 for (it = 0; it < N; it++)
01243 errts_vol[it][iv] = errts[it];
01244
01245
01246 }
01247
01248
01249
01250
01251
01252
01253
01254 void calculate_results
01255 (
01256 WA_options * option_data,
01257 THD_3dim_dataset * dset,
01258 THD_3dim_dataset * mask,
01259 float * fmri_data,
01260 int fmri_length,
01261
01262 float ** coef_vol,
01263 float * mse_vol,
01264 float * ffull_vol,
01265 float * rfull_vol,
01266
01267 float ** coefts_vol,
01268 float ** fitts_vol,
01269 float ** sgnlts_vol,
01270 float ** errts_vol
01271 )
01272
01273 {
01274 float * ts_array = NULL;
01275 float mask_val[1];
01276
01277 int f;
01278 int p;
01279 int q;
01280 int m;
01281 int n;
01282
01283 float * coef = NULL;
01284 float sse_base;
01285 float sse_full;
01286 float ffull;
01287 float rfull;
01288
01289 int ixyz;
01290 int nxyz;
01291
01292 int nt;
01293 int NFirst;
01294 int NLast;
01295 int N;
01296
01297 int i;
01298
01299 float * coefts = NULL;
01300 float * fitts = NULL;
01301 float * sgnlts = NULL;
01302 float * errts = NULL;
01303
01304 char * label = NULL;
01305
01306
01307
01308 if (option_data->input1D_filename != NULL)
01309 {
01310 nxyz = 1;
01311 nt = fmri_length;
01312 }
01313 else
01314 {
01315 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
01316 nt = DSET_NUM_TIMES (dset);
01317 }
01318
01319 NFirst = option_data->NFirst;
01320 NLast = option_data->NLast;
01321 N = option_data->N;
01322 f = option_data->f;
01323 q = option_data->q;
01324 p = option_data->p;
01325
01326
01327
01328 ts_array = (float *) malloc (sizeof(float) * nt); MTEST (ts_array);
01329 coef = (float *) malloc (sizeof(float) * p); MTEST (coef);
01330 coefts = (float *) malloc (sizeof(float) * N); MTEST (coefts);
01331 fitts = (float *) malloc (sizeof(float) * N); MTEST (fitts);
01332 sgnlts = (float *) malloc (sizeof(float) * N); MTEST (sgnlts);
01333 errts = (float *) malloc (sizeof(float) * N); MTEST (errts);
01334
01335
01336
01337 for (ixyz = 0; ixyz < nxyz; ixyz++)
01338 {
01339
01340 if (mask != NULL)
01341 {
01342 extract_ts_array (mask, ixyz, mask_val);
01343 if (mask_val[0] == 0.0) continue;
01344 }
01345
01346
01347
01348 if (option_data->input1D_filename != NULL)
01349 {
01350 for (i = 0; i < N; i++)
01351 ts_array[i] = fmri_data[i+NFirst];
01352 }
01353 else
01354 extract_ts_array (dset, ixyz, ts_array);
01355
01356
01357
01358 wavelet_analysis (option_data->wavelet_type,
01359 f, option_data->stop_filter,
01360 q, option_data->base_filter,
01361 p, option_data->sgnl_filter,
01362 N, ts_array+NFirst, coef,
01363 &sse_base, &sse_full, &ffull, &rfull,
01364 coefts, fitts, sgnlts, errts);
01365
01366
01367
01368 if (option_data->input1D_filename == NULL)
01369 save_voxel (option_data, ixyz, coef, sse_full/(N-f-p), ffull, rfull,
01370 coefts, fitts, sgnlts, errts,
01371 coef_vol, mse_vol, ffull_vol, rfull_vol,
01372 coefts_vol, fitts_vol, sgnlts_vol, errts_vol);
01373 else
01374 {
01375
01376 ts_fprint ("WA.coefts.1D", N, coefts);
01377 ts_fprint ("WA.fitts.1D", N, fitts);
01378 ts_fprint ("WA.sgnlts.1D", N, sgnlts);
01379 ts_fprint ("WA.errts.1D", N, errts);
01380 }
01381
01382
01383
01384 if ( ((ffull > option_data->fdisp) && (option_data->fdisp >= 0.0))
01385 || (option_data->input1D_filename != NULL) )
01386 {
01387 printf ("\n\nResults for Voxel #%d: \n", ixyz);
01388 report_results (N, NFirst, f, p, q,
01389 option_data->base_filter, option_data->sgnl_filter,
01390 coef, sse_base, sse_full, ffull, rfull,
01391 &label);
01392 printf ("%s \n", label);
01393 }
01394
01395
01396 }
01397
01398
01399
01400
01401 free (ts_array); ts_array = NULL;
01402 free (coef); coef = NULL;
01403 free (coefts); coefts = NULL;
01404 free (fitts); fitts = NULL;
01405 free (sgnlts); sgnlts = NULL;
01406 free (errts); errts = NULL;
01407
01408 }
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422 float EDIT_coerce_autoscale_new( int nxyz ,
01423 int itype,void *ivol , int otype,void *ovol )
01424 {
01425 float fac=0.0 , top ;
01426
01427 if( MRI_IS_INT_TYPE(otype) ){
01428 top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
01429 if (top == 0.0) fac = 0.0;
01430 else fac = MRI_TYPE_maxval[otype]/top;
01431 }
01432
01433 EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
01434 return ( fac );
01435 }
01436
01437
01438
01439
01440
01441
01442
01443
01444 void write_ts_array
01445 (
01446 int argc,
01447 char ** argv,
01448 WA_options * option_data,
01449 int ts_length,
01450 float ** vol_array,
01451 char * output_filename
01452 )
01453
01454 {
01455 const float EPSILON = 1.0e-10;
01456
01457 THD_3dim_dataset * dset = NULL;
01458 THD_3dim_dataset * new_dset = NULL;
01459 int ib;
01460 int ierror;
01461 int nxyz;
01462 float factor;
01463 char * input_filename;
01464 short ** bar = NULL;
01465 float * fbuf;
01466 float * volume;
01467 char label[80];
01468
01469
01470
01471 input_filename = option_data->input_filename;
01472 dset = THD_open_one_dataset (input_filename);
01473 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
01474
01475
01476
01477 bar = (short **) malloc (sizeof(short *) * ts_length); MTEST (bar);
01478 fbuf = (float *) malloc (sizeof(float) * ts_length); MTEST (fbuf);
01479
01480
01481
01482 new_dset = EDIT_empty_copy (dset);
01483
01484
01485
01486 tross_Copy_History( dset , new_dset ) ;
01487
01488 { char * commandline = tross_commandline( PROGRAM_NAME , argc , argv ) ;
01489 sprintf (label, "Output prefix: %s", output_filename);
01490 if( commandline != NULL )
01491 tross_multi_Append_History( new_dset , commandline,label,NULL ) ;
01492 else
01493 tross_Append_History ( new_dset, label);
01494 free(commandline) ;
01495 }
01496
01497
01498 THD_delete_3dim_dataset (dset, False); dset = NULL ;
01499
01500
01501 ierror = EDIT_dset_items (new_dset,
01502 ADN_prefix, output_filename,
01503 ADN_label1, output_filename,
01504 ADN_self_name, output_filename,
01505 ADN_malloc_type, DATABLOCK_MEM_MALLOC,
01506 ADN_datum_all, MRI_short,
01507 ADN_nvals, ts_length,
01508 ADN_ntt, ts_length,
01509 ADN_none);
01510
01511 if( ierror > 0 ){
01512 fprintf(stderr,
01513 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
01514 exit(1) ;
01515 }
01516
01517 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
01518 fprintf(stderr,
01519 "*** Output dataset file %s already exists--cannot continue!\a\n",
01520 new_dset->dblk->diskptr->header_name ) ;
01521 exit(1) ;
01522 }
01523
01524
01525
01526 for (ib = 0; ib < ts_length; ib++)
01527 {
01528
01529
01530 volume = vol_array[ib];
01531
01532
01533 bar[ib] = (short *) malloc (sizeof(short) * nxyz);
01534 MTEST (bar[ib]);
01535
01536
01537 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
01538 MRI_short, bar[ib]);
01539 if (factor < EPSILON) factor = 0.0;
01540 else factor = 1.0 / factor;
01541 fbuf[ib] = factor;
01542
01543
01544 mri_fix_data_pointer (bar[ib], DSET_BRICK(new_dset,ib));
01545 }
01546
01547
01548
01549
01550 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
01551
01552 THD_load_statistics (new_dset);
01553 THD_write_3dim_dataset (NULL, NULL, new_dset, True);
01554 fprintf(stderr,"--- Wrote 3d+time dataset into %s\n",DSET_BRIKNAME(new_dset)) ;
01555
01556
01557
01558 THD_delete_3dim_dataset (new_dset, False); new_dset = NULL ;
01559 free (fbuf); fbuf = NULL;
01560
01561 }
01562
01563
01564
01565
01566
01567
01568
01569 void attach_sub_brick
01570 (
01571 THD_3dim_dataset * new_dset,
01572 int ibrick,
01573 float * volume,
01574 int nxyz,
01575 int brick_type,
01576 char * brick_label,
01577 int dof,
01578 int ndof,
01579 int ddof,
01580 short ** bar
01581 )
01582
01583 {
01584 const float EPSILON = 1.0e-10;
01585 float factor;
01586
01587
01588
01589 bar[ibrick] = (short *) malloc (sizeof(short) * nxyz);
01590 MTEST (bar[ibrick]);
01591 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
01592 MRI_short, bar[ibrick]);
01593
01594 if (factor < EPSILON) factor = 0.0;
01595 else factor = 1.0 / factor;
01596
01597
01598
01599 EDIT_BRICK_LABEL (new_dset, ibrick, brick_label);
01600 EDIT_BRICK_FACTOR (new_dset, ibrick, factor);
01601
01602 if (brick_type == FUNC_TT_TYPE)
01603 EDIT_BRICK_TO_FITT (new_dset, ibrick, dof);
01604 else if (brick_type == FUNC_FT_TYPE)
01605 EDIT_BRICK_TO_FIFT (new_dset, ibrick, ndof, ddof);
01606
01607
01608
01609 EDIT_substitute_brick (new_dset, ibrick, MRI_short, bar[ibrick]);
01610
01611 }
01612
01613
01614
01615
01616
01617
01618 void write_bucket_data
01619 (
01620 int argc,
01621 char ** argv,
01622 WA_options * option_data,
01623
01624 float ** coef_vol,
01625 float * mse_vol,
01626 float * ffull_vol,
01627 float * rfull_vol
01628 )
01629
01630 {
01631 THD_3dim_dataset * old_dset = NULL;
01632 THD_3dim_dataset * new_dset = NULL;
01633 char output_prefix[MAX_NAME_LENGTH];
01634 char output_session[MAX_NAME_LENGTH];
01635 int nbricks;
01636 short ** bar = NULL;
01637
01638 int brick_type;
01639 int brick_coef;
01640 char brick_label[MAX_NAME_LENGTH];
01641
01642 int ierror;
01643 float * volume;
01644
01645 int N;
01646 int NFirst;
01647 int f;
01648 int p;
01649 int q;
01650 int nxyz;
01651 int nt;
01652 int it;
01653 int ilag;
01654 int icoef;
01655 int ibrick;
01656 int dof, ndof, ddof;
01657 char label[MAX_NAME_LENGTH];
01658
01659
01660
01661 old_dset = THD_open_one_dataset (option_data->input_filename);
01662
01663
01664
01665 nxyz = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;
01666 nt = DSET_NUM_TIMES (old_dset);
01667
01668
01669 f = option_data->f;
01670 q = option_data->q;
01671 p = option_data->p;
01672 N = option_data->N;
01673 NFirst = option_data->NFirst;
01674 printf ("f = %d q = %d p = %d N = %d \n", f, q, p, N);
01675
01676
01677
01678 nbricks = p * option_data->cout
01679 + option_data->rout + option_data->fout + option_data->vout;
01680 printf ("nbricks = %d \n", nbricks);
01681
01682
01683 strcpy (output_prefix, option_data->bucket_filename);
01684 strcpy (output_session, "./");
01685
01686
01687
01688 bar = (short **) malloc (sizeof(short *) * nbricks);
01689 MTEST (bar);
01690
01691
01692
01693 new_dset = EDIT_empty_copy (old_dset);
01694
01695
01696
01697 tross_Copy_History( old_dset , new_dset ) ;
01698 tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
01699 sprintf (label, "Output prefix: %s", output_prefix);
01700 tross_Append_History ( new_dset, label);
01701
01702
01703
01704 THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ;
01705
01706
01707
01708
01709 ierror = EDIT_dset_items (new_dset,
01710 ADN_prefix, output_prefix,
01711 ADN_directory_name, output_session,
01712 ADN_type, HEAD_FUNC_TYPE,
01713 ADN_func_type, FUNC_BUCK_TYPE,
01714 ADN_datum_all, MRI_short ,
01715 ADN_ntt, 0,
01716 ADN_nvals, nbricks,
01717 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
01718 ADN_none ) ;
01719
01720 if( ierror > 0 )
01721 {
01722 fprintf(stderr,
01723 "*** %d errors in attempting to create bucket dataset!\n",
01724 ierror);
01725 exit(1);
01726 }
01727
01728 if (THD_is_file(DSET_HEADNAME(new_dset)))
01729 {
01730 fprintf(stderr,
01731 "*** Output dataset file %s already exists--cannot continue!\n",
01732 DSET_HEADNAME(new_dset));
01733 exit(1);
01734 }
01735
01736
01737
01738 ibrick = 0;
01739
01740
01741
01742 if (option_data->cout)
01743 {
01744
01745 if (option_data->stat_first)
01746 ibrick += option_data->vout + option_data->rout + option_data->fout;
01747
01748 icoef = 0;
01749 for (it = 0; it < N; it++)
01750 {
01751 if (option_data->sgnl_filter[it])
01752 {
01753
01754 brick_type = FUNC_FIM_TYPE;
01755
01756 {
01757 int band, mintr, maxtr;
01758
01759 if (it == 0)
01760 {
01761 band = -1;
01762 mintr = 0;
01763 maxtr = N-1;
01764 }
01765 else
01766 {
01767 band = my_log2(it);
01768 mintr = (it - powerof2(band)) * powerof2(my_log2(N)-band);
01769 maxtr = mintr + powerof2(my_log2(N)-band) - 1;
01770 }
01771
01772 mintr += NFirst;
01773 maxtr += NFirst;
01774
01775 if (option_data->base_filter[it])
01776 sprintf (brick_label, "B(%2d)[%3d,%3d]", band, mintr, maxtr);
01777 else
01778 sprintf (brick_label, "S(%2d)[%3d,%3d]", band, mintr, maxtr);
01779 }
01780
01781 volume = coef_vol[icoef];
01782 attach_sub_brick (new_dset, ibrick, volume, nxyz,
01783 brick_type, brick_label, 0, 0, 0, bar);
01784
01785 ibrick++;
01786 icoef++;
01787 }
01788
01789 }
01790 }
01791
01792
01793
01794
01795 if (option_data->stat_first) ibrick = 0;
01796
01797
01798 if (option_data->vout)
01799 {
01800 brick_type = FUNC_FIM_TYPE;
01801 sprintf (brick_label, "Full MSE");
01802 volume = mse_vol;
01803 attach_sub_brick (new_dset, ibrick, volume, nxyz,
01804 brick_type, brick_label, 0, 0, 0, bar);
01805 ibrick++;
01806 }
01807
01808
01809 if (option_data->rout)
01810 {
01811 brick_type = FUNC_THR_TYPE;
01812 sprintf (brick_label, "Full R^2");
01813 volume = rfull_vol;
01814 attach_sub_brick (new_dset, ibrick, volume, nxyz,
01815 brick_type, brick_label, 0, 0, 0, bar);
01816 ibrick++;
01817 }
01818
01819
01820 if (option_data->fout)
01821 {
01822 brick_type = FUNC_FT_TYPE;
01823 ndof = p - q;
01824 ddof = N - f - p;
01825 sprintf (brick_label, "Full F-stat");
01826 volume = ffull_vol;
01827 attach_sub_brick (new_dset, ibrick, volume, nxyz,
01828 brick_type, brick_label, 0, ndof, ddof, bar);
01829 ibrick++;
01830 }
01831
01832
01833
01834 THD_load_statistics (new_dset);
01835 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01836 fprintf(stderr,"Writing `bucket' dataset ");
01837 fprintf(stderr,"into %s\n", DSET_BRIKNAME(new_dset));
01838
01839
01840
01841 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01842
01843 }
01844
01845
01846
01847
01848
01849
01850
01851 void output_results
01852 (
01853 int argc,
01854 char ** argv,
01855 WA_options * option_data,
01856
01857 float ** coef_vol,
01858 float * mse_vol,
01859 float * ffull_vol,
01860 float * rfull_vol,
01861 float ** coefts_vol,
01862 float ** fitts_vol,
01863 float ** sgnlts_vol,
01864 float ** errts_vol
01865 )
01866
01867 {
01868 int N;
01869
01870
01871
01872 N = option_data->N;
01873
01874
01875
01876 if (option_data->bucket_filename != NULL)
01877 write_bucket_data (argc, argv, option_data, coef_vol,
01878 mse_vol, ffull_vol, rfull_vol);
01879
01880
01881
01882
01883 if (option_data->coefts_filename != NULL)
01884 write_ts_array (argc, argv, option_data, N, coefts_vol,
01885 option_data->coefts_filename);
01886
01887
01888
01889 if (option_data->fitts_filename != NULL)
01890 write_ts_array (argc, argv, option_data, N, fitts_vol,
01891 option_data->fitts_filename);
01892
01893
01894
01895 if (option_data->sgnlts_filename != NULL)
01896 write_ts_array (argc, argv, option_data, N, sgnlts_vol,
01897 option_data->sgnlts_filename);
01898
01899
01900
01901 if (option_data->errts_filename != NULL)
01902 write_ts_array (argc, argv, option_data, N, errts_vol,
01903 option_data->errts_filename);
01904
01905
01906
01907 }
01908
01909
01910
01911
01912 void terminate_program
01913 (
01914 WA_options ** option_data,
01915 float ** fmri_data,
01916
01917 float *** coef_vol,
01918 float ** mse_vol,
01919 float ** ffull_vol,
01920 float ** rfull_vol,
01921
01922 float *** coefts_vol,
01923 float *** fitts_vol,
01924 float *** sgnlts_vol,
01925 float *** errts_vol
01926 )
01927
01928 {
01929 int p;
01930 int ip;
01931 int is;
01932 int it;
01933 int N;
01934
01935
01936
01937 p = (*option_data)->p;
01938 N = (*option_data)->N;
01939
01940
01941
01942 if ((*option_data)->stop_filter != NULL) free ((*option_data)->stop_filter);
01943 if ((*option_data)->base_filter != NULL) free ((*option_data)->base_filter);
01944 if ((*option_data)->sgnl_filter != NULL) free ((*option_data)->sgnl_filter);
01945 if ((*option_data)->input_filename != NULL)
01946 free ((*option_data)->input_filename);
01947 if ((*option_data)->mask_filename != NULL)
01948 free ((*option_data)->mask_filename);
01949 if ((*option_data)->input1D_filename != NULL)
01950 free ((*option_data)->input1D_filename);
01951 if ((*option_data)->bucket_filename != NULL)
01952 free ((*option_data)->bucket_filename);
01953 if ((*option_data)->coefts_filename != NULL)
01954 free ((*option_data)->coefts_filename);
01955 if ((*option_data)->fitts_filename != NULL)
01956 free ((*option_data)->fitts_filename);
01957 if ((*option_data)->sgnlts_filename != NULL)
01958 free ((*option_data)->sgnlts_filename);
01959 if ((*option_data)->errts_filename != NULL)
01960 free ((*option_data)->errts_filename);
01961 free (*option_data); *option_data = NULL;
01962
01963
01964
01965 if (*fmri_data != NULL) { free (*fmri_data); *fmri_data = NULL; }
01966
01967
01968
01969 if (*coef_vol != NULL)
01970 {
01971 for (ip = 0; ip < p; ip++)
01972 {
01973 if ((*coef_vol)[ip] != NULL)
01974 { free ((*coef_vol)[ip]); (*coef_vol)[ip] = NULL; }
01975 }
01976 free (*coef_vol); *coef_vol = NULL;
01977 }
01978
01979 if (*mse_vol != NULL) { free (*mse_vol); *mse_vol = NULL; }
01980 if (*ffull_vol != NULL) { free (*ffull_vol); *ffull_vol = NULL; }
01981 if (*rfull_vol != NULL) { free (*rfull_vol); *rfull_vol = NULL; }
01982
01983
01984
01985 if (*coefts_vol != NULL)
01986 {
01987 for (it = 0; it < N; it++)
01988 { free ((*coefts_vol)[it]); (*coefts_vol)[it] = NULL; }
01989 free (*coefts_vol); *coefts_vol = NULL;
01990 }
01991
01992
01993
01994 if (*fitts_vol != NULL)
01995 {
01996 for (it = 0; it < N; it++)
01997 { free ((*fitts_vol)[it]); (*fitts_vol)[it] = NULL; }
01998 free (*fitts_vol); *fitts_vol = NULL;
01999 }
02000
02001
02002
02003 if (*sgnlts_vol != NULL)
02004 {
02005 for (it = 0; it < N; it++)
02006 { free ((*sgnlts_vol)[it]); (*sgnlts_vol)[it] = NULL; }
02007 free (*sgnlts_vol); *sgnlts_vol = NULL;
02008 }
02009
02010
02011
02012 if (*errts_vol != NULL)
02013 {
02014 for (it = 0; it < N; it++)
02015 { free ((*errts_vol)[it]); (*errts_vol)[it] = NULL; }
02016 free (*errts_vol); *errts_vol = NULL;
02017 }
02018
02019
02020 }
02021
02022
02023
02024
02025 int main
02026 (
02027 int argc,
02028 char ** argv
02029 )
02030
02031 {
02032 WA_options * option_data;
02033 THD_3dim_dataset * dset_time = NULL;
02034 THD_3dim_dataset * mask_dset = NULL;
02035 float * fmri_data = NULL;
02036 int fmri_length;
02037
02038 float ** coef_vol = NULL;
02039 float * mse_vol = NULL;
02040 float * ffull_vol = NULL;
02041 float * rfull_vol = NULL;
02042
02043 float ** coefts_vol = NULL;
02044 float ** fitts_vol = NULL;
02045 float ** sgnlts_vol = NULL;
02046 float ** errts_vol = NULL;
02047
02048
02049
02050 printf ("\n\n");
02051 printf ("Program: %s \n", PROGRAM_NAME);
02052 printf ("Author: %s \n", PROGRAM_AUTHOR);
02053 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
02054 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
02055 printf ("\n");
02056
02057
02058
02059 initialize_program (argc, argv, &option_data,
02060 &dset_time, &mask_dset, &fmri_data, &fmri_length,
02061 &coef_vol, &mse_vol, &ffull_vol, &rfull_vol,
02062 &coefts_vol, &fitts_vol, &sgnlts_vol, &errts_vol);
02063
02064
02065
02066 calculate_results (option_data, dset_time, mask_dset, fmri_data, fmri_length,
02067 coef_vol, mse_vol, ffull_vol, rfull_vol,
02068 coefts_vol, fitts_vol, sgnlts_vol, errts_vol);
02069
02070
02071
02072 if (dset_time != NULL)
02073 { THD_delete_3dim_dataset (dset_time, False); dset_time = NULL; }
02074 if (mask_dset != NULL)
02075 { THD_delete_3dim_dataset (mask_dset, False); mask_dset = NULL; }
02076
02077
02078
02079 if (option_data->input1D_filename == NULL)
02080 output_results (argc, argv, option_data, coef_vol, mse_vol,
02081 ffull_vol, rfull_vol,
02082 coefts_vol, fitts_vol, sgnlts_vol, errts_vol);
02083
02084
02085
02086 terminate_program (&option_data, &fmri_data,
02087 &coef_vol, &mse_vol, &ffull_vol, &rfull_vol,
02088 &coefts_vol, &fitts_vol, &sgnlts_vol, &errts_vol);
02089
02090 exit(0);
02091 }
02092
02093
02094
02095
02096