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
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 #define PROGRAM_NAME "3dNLfim"
00083 #define PROGRAM_AUTHOR "B. Douglas Ward"
00084 #define PROGRAM_INITIAL "19 June 1997"
00085 #define PROGRAM_LATEST "07 May 2003"
00086
00087
00088
00089 #include <stdio.h>
00090 #include <math.h>
00091 #include <stdlib.h>
00092
00093 #include "mrilib.h"
00094
00095
00096
00097
00098
00099 #if !defined(DONT_USE_SHM) && !defined(DONT_USE_FORK) && !defined(CYGWIN)
00100
00101 # include "thd_iochan.h"
00102
00103 # define PROC_MAX 32
00104
00105 static int proc_numjob = 1 ;
00106 static pid_t proc_pid[PROC_MAX] ;
00107 static int proc_shmid = 0 ;
00108 static char *proc_shmptr = NULL;
00109 static int proc_shmsize = 0 ;
00110
00111 static int proc_shm_arnum = 0 ;
00112 static float ***proc_shm_ar = NULL;
00113 static int *proc_shm_arsiz = NULL;
00114
00115 static int proc_vox_bot[PROC_MAX] ;
00116 static int proc_vox_top[PROC_MAX] ;
00117
00118 static int proc_ind ;
00119
00120 #else
00121
00122 # define proc_numjob 1
00123 # define proc_ind 0
00124
00125 #endif
00126
00127 #include "matrix.h"
00128 #include "simplex.h"
00129 #include "NLfit.h"
00130
00131 #include "matrix.c"
00132 #include "simplex.c"
00133 #include "NLfit.c"
00134
00135 typedef struct NL_options
00136 {
00137 char * bucket_filename;
00138 int numbricks;
00139 int * brick_type;
00140 int * brick_coef;
00141 char ** brick_label;
00142
00143 } NL_options;
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 static float DELT = 1.0;
00155 static int inTR = 0 ;
00156 static float dsTR = 0.0 ;
00157
00158 static char * commandline = NULL ;
00159
00160 static byte * mask_vol = NULL;
00161 static int mask_nvox = 0;
00162
00163
00164
00165
00166
00167
00168
00169 void display_help_menu()
00170 {
00171 printf
00172 (
00173 "This program calculates a nonlinear regression for each voxel of the \n"
00174 "input AFNI 3d+time data set. The nonlinear regression is calculated \n"
00175 "by means of a least squares fit to the signal plus noise models which \n"
00176 "are specified by the user. \n"
00177 " \n"
00178 "Usage: \n"
00179 "3dNLfim \n"
00180 "-input fname fname = filename of 3d + time data file for input \n"
00181 "[-mask mset] Use the 0 sub-brick of dataset 'mset' as a mask \n"
00182 " to indicate which voxels to analyze (a sub-brick \n"
00183 " selector is allowed) [default = use all voxels] \n"
00184 "[-ignore num] num = skip this number of initial images in the \n"
00185 " time series for regresion analysis; default = 3 \n"
00186 "[-inTR] set delt = TR of the input 3d+time dataset \n"
00187 " [The default is to compute with delt = 1.0 ] \n"
00188 " [The model functions are calculated using a \n"
00189 " time grid of: 0, delt, 2*delt, 3*delt, ... ] \n"
00190 "[-time fname] fname = ASCII file containing each time point \n"
00191 " in the time series. Defaults to even spacing \n"
00192 " given by TR (this option overrides -inTR). \n"
00193 "-signal slabel slabel = name of (non-linear) signal model \n"
00194 "-noise nlabel nlabel = name of (linear) noise model \n"
00195 "-sconstr k c d constraints for kth signal parameter: \n"
00196 " c <= gs[k] <= d \n"
00197 "-nconstr k c d constraints for kth noise parameter: \n"
00198 " c+b[k] <= gn[k] <= d+b[k] \n"
00199 "[-nabs] use absolute constraints for noise parameters: \n"
00200 " c <= gn[k] <= d \n"
00201 "[-nrand n] n = number of random test points \n"
00202 "[-nbest b] b = find opt. soln. for b best test points \n"
00203 "[-rmsmin r] r = minimum rms error to reject reduced model \n"
00204 "[-fdisp fval] display (to screen) results for those voxels \n"
00205 " whose f-statistic is > fval \n"
00206 " \n"
00207 " \n"
00208 "The following commands generate individual AFNI 2 sub-brick datasets: \n"
00209 " \n"
00210 "[-freg fname] perform f-test for significance of the regression; \n"
00211 " output 'fift' is written to prefix filename fname\n"
00212 "[-frsqr fname] calculate R^2 (coef. of multiple determination); \n"
00213 " store along with f-test for regression; \n"
00214 " output 'fift' is written to prefix filename fname\n"
00215 "[-fsmax fname] estimate signed maximum of signal; store along \n"
00216 " with f-test for regression; output 'fift' is \n"
00217 " written to prefix filename fname \n"
00218 "[-ftmax fname] estimate time of signed maximum; store along \n"
00219 " with f-test for regression; output 'fift' is \n"
00220 " written to prefix filename fname \n"
00221 "[-fpsmax fname] calculate (signed) maximum percentage change of \n"
00222 " signal from baseline; output 'fift' is \n"
00223 " written to prefix filename fname \n"
00224 "[-farea fname] calculate area between signal and baseline; store \n"
00225 " with f-test for regression; output 'fift' is \n"
00226 " written to prefix filename fname \n"
00227 "[-fparea fname] percentage area of signal relative to baseline; \n"
00228 " store with f-test for regression; output 'fift' \n"
00229 " is written to prefix filename fname \n"
00230 "[-fscoef k fname] estimate kth signal parameter gs[k]; store along \n"
00231 " with f-test for regression; output 'fift' is \n"
00232 " written to prefix filename fname \n"
00233 "[-fncoef k fname] estimate kth noise parameter gn[k]; store along \n"
00234 " with f-test for regression; output 'fift' is \n"
00235 " written to prefix filename fname \n"
00236 "[-tscoef k fname] perform t-test for significance of the kth signal \n"
00237 " parameter gs[k]; output 'fitt' is written \n"
00238 " to prefix filename fname \n"
00239 "[-tncoef k fname] perform t-test for significance of the kth noise \n"
00240 " parameter gn[k]; output 'fitt' is written \n"
00241 " to prefix filename fname \n"
00242 " \n"
00243 " \n"
00244 "The following commands generate one AFNI 'bucket' type dataset: \n"
00245 " \n"
00246 "[-bucket n prefixname] create one AFNI 'bucket' dataset containing \n"
00247 " n sub-bricks; n=0 creates default output; \n"
00248 " output 'bucket' is written to prefixname \n"
00249 "The mth sub-brick will contain: \n"
00250 "[-brick m scoef k label] kth signal parameter regression coefficient\n"
00251 "[-brick m ncoef k label] kth noise parameter regression coefficient \n"
00252 "[-brick m tmax label] time at max. abs. value of signal \n"
00253 "[-brick m smax label] signed max. value of signal \n"
00254 "[-brick m psmax label] signed max. value of signal as percent \n"
00255 " above baseline level \n"
00256 "[-brick m area label] area between signal and baseline \n"
00257 "[-brick m parea label] signed area between signal and baseline \n"
00258 " as percent of baseline area \n"
00259 "[-brick m tscoef k label] t-stat for kth signal parameter coefficient\n"
00260 "[-brick m tncoef k label] t-stat for kth noise parameter coefficient \n"
00261 "[-brick m resid label] std. dev. of the full model fit residuals \n"
00262 "[-brick m rsqr label] R^2 (coefficient of multiple determination)\n"
00263 "[-brick m fstat label] F-stat for significance of the regression \n"
00264 " \n"
00265 " \n"
00266 "The following commands write the time series fit for each voxel \n"
00267 "to an AFNI 3d+time dataset: \n"
00268 "[-sfit fname] fname = prefix for output 3d+time signal model fit \n"
00269 "[-snfit fname] fname = prefix for output 3d+time signal+noise fit \n"
00270 " \n"
00271 );
00272
00273
00274 #ifdef PROC_MAX
00275 printf( "\n"
00276 " -jobs J Run the program with 'J' jobs (sub-processes).\n"
00277 " On a multi-CPU machine, this can speed the\n"
00278 " program up considerably. On a single CPU\n"
00279 " machine, using this option is silly.\n"
00280 " J should be a number from 1 up to the\n"
00281 " number of CPU sharing memory on the system.\n"
00282 " J=1 is normal (single process) operation.\n"
00283 " The maximum allowed value of J is %d.\n"
00284 " * For more information on parallelizing, see\n"
00285 " http://afni.nimh.nih.gov/afni/doc/misc/parallize.html\n"
00286 " * Use -mask to get more speed; cf. 3dAutomask.\n"
00287 , PROC_MAX ) ;
00288 #endif
00289
00290 exit(0);
00291 }
00292
00293
00294
00295
00296
00297
00298 #define MTEST(ptr) \
00299 if((ptr)==NULL) \
00300 ( NLfit_error ("Cannot allocate memory") )
00301
00302
00303
00304
00305
00306
00307
00308 void initialize_options
00309 (
00310 int * ignore,
00311 vfp * nmodel,
00312 vfp * smodel,
00313 int * r,
00314 int * p,
00315 char *** npname,
00316 char *** spname,
00317 float ** min_nconstr,
00318 float ** max_nconstr,
00319 float ** min_sconstr,
00320 float ** max_sconstr,
00321 int * nabs,
00322 int * nrand,
00323 int * nbest,
00324 float * rms_min,
00325 float * fdisp,
00326 char ** input_filename,
00327 char ** tfilename,
00328 char ** freg_filename,
00329 char ** frsqr_filename,
00330 char *** fncoef_filename,
00331 char *** fscoef_filename,
00332 char *** tncoef_filename,
00333 char *** tscoef_filename,
00334 char ** sfit_filename,
00335 char ** snfit_filename,
00336 NL_options * option_data
00337 )
00338
00339 {
00340 int ip;
00341
00342
00343
00344 *ignore = 3;
00345 *nabs = 0;
00346 *nrand = 100;
00347 *nbest = 5;
00348 *rms_min = 0.0;
00349 *fdisp = 10.0;
00350 *smodel = NULL;
00351 *nmodel = NULL;
00352 *r = -1;
00353 *p = -1;
00354
00355
00356
00357 *npname = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00358 MTEST (*npname);
00359 for (ip = 0; ip < MAX_PARAMETERS; ip++)
00360 {
00361 (*npname)[ip] = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
00362 MTEST ((*npname)[ip]);
00363 }
00364
00365
00366
00367 *spname = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00368 MTEST (*spname);
00369 for (ip = 0; ip < MAX_PARAMETERS; ip++)
00370 {
00371 (*spname)[ip] = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
00372 MTEST ((*spname)[ip]);
00373 }
00374
00375
00376
00377 *min_nconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
00378 MTEST (*min_nconstr);
00379 *max_nconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
00380 MTEST (*max_nconstr);
00381 *min_sconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
00382 MTEST (*min_sconstr);
00383 *max_sconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
00384 MTEST (*max_sconstr);
00385
00386
00387
00388 *input_filename = NULL;
00389 *tfilename = NULL;
00390 *freg_filename = NULL;
00391 *frsqr_filename = NULL;
00392 *fncoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00393 MTEST (*fncoef_filename);
00394 *fscoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00395 MTEST (*fscoef_filename);
00396 *tncoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00397 MTEST (*tncoef_filename);
00398 *tscoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00399 MTEST (*tscoef_filename);
00400 for (ip = 0; ip < MAX_PARAMETERS; ip++)
00401 {
00402 (*fncoef_filename)[ip] = NULL;
00403 (*fscoef_filename)[ip] = NULL;
00404 (*tncoef_filename)[ip] = NULL;
00405 (*tscoef_filename)[ip] = NULL;
00406 }
00407 *sfit_filename = NULL;
00408 *snfit_filename = NULL;
00409
00410
00411
00412 option_data->bucket_filename = NULL;
00413 option_data->numbricks = -1;
00414 option_data->brick_type = NULL;
00415 option_data->brick_coef = NULL;
00416 option_data->brick_label = NULL;
00417
00418 }
00419
00420
00421
00422
00423
00424
00425
00426 void get_options
00427 (
00428 int argc,
00429 char ** argv,
00430 int * ignore,
00431 char ** nname,
00432 char ** sname,
00433 vfp * nmodel,
00434 vfp * smodel,
00435 int * r,
00436 int * p,
00437 char *** npname,
00438 char *** spname,
00439 float ** min_nconstr,
00440 float ** max_nconstr,
00441 float ** min_sconstr,
00442 float ** max_sconstr,
00443 int * nabs,
00444 int * nrand,
00445 int * nbest,
00446 float * rms_min,
00447 float * fdisp,
00448 char ** input_filename,
00449 char ** tfilename,
00450 char ** freg_filename,
00451 char ** frsqr_filename,
00452 char ** fsmax_filename,
00453 char ** ftmax_filename,
00454 char ** fpmax_filename,
00455 char ** farea_filename,
00456 char ** fparea_filename,
00457 char *** fncoef_filename,
00458 char *** fscoef_filename,
00459 char *** tncoef_filename,
00460 char *** tscoef_filename,
00461 char ** sfit_filename,
00462 char ** snfit_filename,
00463
00464 THD_3dim_dataset ** dset_time,
00465 int * nxyz,
00466 int * ts_length,
00467 NL_options * option_data
00468 )
00469
00470 {
00471 const MAX_BRICKS = 100;
00472 int nopt = 1;
00473 int ival, index;
00474 float fval;
00475 char message[MAX_NAME_LENGTH];
00476 int ok;
00477
00478 NLFIT_MODEL_array * model_array = NULL;
00479 int im;
00480 int ibrick;
00481 int nbricks;
00482
00483
00484
00485 if (argc < 2 || strcmp(argv[1], "-help") == 0) display_help_menu();
00486
00487
00488
00489 AFNI_logger (PROGRAM_NAME,argc,argv);
00490
00491
00492
00493 model_array = NLFIT_get_many_MODELs ();
00494 if ((model_array == NULL) || (model_array->num == 0))
00495 NLfit_error ("Unable to locate any models");
00496
00497
00498 initialize_options (ignore, nmodel, smodel, r, p, npname, spname,
00499 min_nconstr, max_nconstr, min_sconstr, max_sconstr, nabs,
00500 nrand, nbest, rms_min, fdisp,
00501 input_filename, tfilename, freg_filename,
00502 frsqr_filename, fncoef_filename, fscoef_filename,
00503 tncoef_filename, tscoef_filename,
00504 sfit_filename, snfit_filename, option_data);
00505
00506
00507
00508 while (nopt < argc )
00509 {
00510
00511
00512 if (strcmp(argv[nopt], "-input") == 0)
00513 {
00514 nopt++;
00515 if (nopt >= argc) NLfit_error ("need argument after -input ");
00516 *input_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00517 MTEST (*input_filename);
00518 strcpy (*input_filename, argv[nopt]);
00519 *dset_time = THD_open_one_dataset (*input_filename);
00520 if ((*dset_time) == NULL)
00521 {
00522 sprintf (message,
00523 "Unable to open data file: %s", *input_filename);
00524 NLfit_error (message);
00525 }
00526
00527 THD_load_datablock ((*dset_time)->dblk);
00528
00529 *nxyz = (*dset_time)->dblk->diskptr->dimsizes[0]
00530 * (*dset_time)->dblk->diskptr->dimsizes[1]
00531 * (*dset_time)->dblk->diskptr->dimsizes[2] ;
00532 *ts_length = DSET_NUM_TIMES(*dset_time);
00533
00534 dsTR = DSET_TIMESTEP(*dset_time) ;
00535 if( DSET_TIMEUNITS(*dset_time) == UNITS_MSEC_TYPE ) dsTR *= 0.001 ;
00536
00537 nopt++;
00538 continue;
00539 }
00540
00541
00542
00543
00544 if( strcmp(argv[nopt],"-mask") == 0 ){
00545 THD_3dim_dataset * mset ; int ii,mc ;
00546 nopt++ ;
00547 if (nopt >= argc) NLfit_error ("need argument after -mask!") ;
00548 mset = THD_open_dataset( argv[nopt] ) ;
00549 if (mset == NULL) NLfit_error ("can't open -mask dataset!") ;
00550
00551 mask_vol = THD_makemask( mset , 0 , 1.0,0.0 ) ;
00552 mask_nvox = DSET_NVOX(mset) ;
00553 DSET_delete(mset) ;
00554
00555 if (mask_vol == NULL ) NLfit_error ("can't use -mask dataset!") ;
00556 for( ii=mc=0 ; ii < mask_nvox ; ii++ ) if (mask_vol[ii]) mc++ ;
00557 if (mc == 0) NLfit_error ("mask is all zeros!") ;
00558 printf("++ %d voxels in mask %s\n",mc,argv[nopt]) ;
00559 nopt++ ; continue ;
00560 }
00561
00562
00563
00564
00565 if( strcmp(argv[nopt],"-inTR") == 0 ){
00566 inTR = 1 ;
00567 nopt++ ; continue ;
00568 }
00569
00570
00571 if (strcmp(argv[nopt], "-ignore") == 0)
00572 {
00573 nopt++;
00574 if (nopt >= argc) NLfit_error ("need argument after -ignore ");
00575 sscanf (argv[nopt], "%d", &ival);
00576 if (ival < 0)
00577 NLfit_error ("illegal argument after -ignore ");
00578 *ignore = ival;
00579 nopt++;
00580 continue;
00581 }
00582
00583
00584 if (strcmp(argv[nopt], "-time") == 0)
00585 {
00586 nopt++;
00587 if (nopt >= argc) NLfit_error ("need argument after -time ");
00588 *tfilename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00589 MTEST (*tfilename);
00590 strcpy (*tfilename, argv[nopt]);
00591 nopt++;
00592 continue;
00593 }
00594
00595
00596
00597 if (strcmp(argv[nopt], "-signal") == 0)
00598 {
00599 nopt++;
00600 if (nopt >= argc) NLfit_error ("need argument after -signal ");
00601 *sname = malloc (sizeof(char) * MAX_NAME_LENGTH);
00602 MTEST (*sname);
00603 strcpy (*sname, argv[nopt]);
00604 initialize_signal_model (model_array, *sname,
00605 smodel, p, *spname,
00606 *min_sconstr, *max_sconstr);
00607 nopt++;
00608 continue;
00609 }
00610
00611
00612
00613 if (strcmp(argv[nopt], "-noise") == 0)
00614 {
00615 nopt++;
00616 if (nopt >= argc) NLfit_error ("need argument after -noise ");
00617 *nname = malloc (sizeof(char) * MAX_NAME_LENGTH);
00618 MTEST (*nname);
00619 strcpy (*nname, argv[nopt]);
00620 initialize_noise_model (model_array, *nname,
00621 nmodel, r, *npname,
00622 *min_nconstr, *max_nconstr);
00623 nopt++;
00624 continue;
00625 }
00626
00627
00628
00629 if ((*smodel) == NULL) NLfit_error ("Must specify signal model");
00630 if ((*nmodel) == NULL) NLfit_error ("Must specify noise model");
00631
00632
00633
00634 if (strcmp(argv[nopt], "-sconstr") == 0)
00635 {
00636 nopt++;
00637 if (nopt+2 >= argc) NLfit_error("need 3 arguments after -sconstr ");
00638
00639 sscanf (argv[nopt], "%d", &ival);
00640 if ((ival < 0) || (ival >= *p))
00641 NLfit_error ("illegal argument after -sconstr ");
00642 index = ival;
00643 nopt++;
00644
00645 sscanf (argv[nopt], "%f", &fval);
00646 (*min_sconstr)[index] = fval;
00647 nopt++;
00648
00649 sscanf (argv[nopt], "%f", &fval);
00650 (*max_sconstr)[index] = fval;
00651 nopt++;
00652 continue;
00653 }
00654
00655
00656
00657 if (strcmp(argv[nopt], "-nconstr") == 0)
00658 {
00659 nopt++;
00660 if (nopt+2 >= argc) NLfit_error("need 3 arguments after -nconstr ");
00661
00662 sscanf (argv[nopt], "%d", &ival);
00663 if ((ival < 0) || (ival >= *r))
00664 NLfit_error ("illegal argument after -nconstr ");
00665 index = ival;
00666 nopt++;
00667
00668 sscanf (argv[nopt], "%f", &fval);
00669 (*min_nconstr)[index] = fval;
00670 nopt++;
00671
00672 sscanf (argv[nopt], "%f", &fval);
00673 (*max_nconstr)[index] = fval;
00674 nopt++;
00675 continue;
00676 }
00677
00678
00679
00680 if (strcmp(argv[nopt], "-nabs") == 0)
00681 {
00682 *nabs = 1;
00683 nopt++;
00684 continue;
00685 }
00686
00687
00688
00689 if (strcmp(argv[nopt], "-nrand") == 0)
00690 {
00691 nopt++;
00692 if (nopt >= argc) NLfit_error ("need argument after -nrand ");
00693 sscanf (argv[nopt], "%d", &ival);
00694 if (ival <= 0)
00695 NLfit_error ("illegal argument after -nrand ");
00696 *nrand = ival;
00697 nopt++;
00698 continue;
00699 }
00700
00701
00702
00703 if (strcmp(argv[nopt], "-nbest") == 0)
00704 {
00705 nopt++;
00706 if (nopt >= argc) NLfit_error ("need argument after -nbest ");
00707 sscanf (argv[nopt], "%d", &ival);
00708 if (ival <= 0)
00709 NLfit_error ("illegal argument after -nbest ");
00710 *nbest = ival;
00711 nopt++;
00712 continue;
00713 }
00714
00715
00716
00717 if (strcmp(argv[nopt], "-rmsmin") == 0)
00718 {
00719 nopt++;
00720 if (nopt >= argc) NLfit_error ("need argument after -rmsmin ");
00721 sscanf (argv[nopt], "%f", &fval);
00722 if (fval < 0.0)
00723 NLfit_error ("illegal argument after -rmsmin ");
00724 *rms_min = fval;
00725 nopt++;
00726 continue;
00727 }
00728
00729
00730
00731 if (strcmp(argv[nopt], "-fdisp") == 0)
00732 {
00733 nopt++;
00734 if (nopt >= argc) NLfit_error ("need argument after -fdisp ");
00735 sscanf (argv[nopt], "%f", &fval);
00736 *fdisp = fval;
00737 nopt++;
00738 continue;
00739 }
00740
00741
00742
00743 if (strcmp(argv[nopt], "-freg") == 0)
00744 {
00745 nopt++;
00746 if (nopt >= argc) NLfit_error ("need argument after -freg ");
00747 *freg_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00748 MTEST (*freg_filename);
00749 strcpy (*freg_filename, argv[nopt]);
00750 nopt++;
00751 continue;
00752 }
00753
00754
00755
00756 if (strcmp(argv[nopt], "-frsqr") == 0)
00757 {
00758 nopt++;
00759 if (nopt >= argc) NLfit_error ("need argument after -frsqr ");
00760 *frsqr_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00761 MTEST (*frsqr_filename);
00762 strcpy (*frsqr_filename, argv[nopt]);
00763 nopt++;
00764 continue;
00765 }
00766
00767
00768
00769 if (strcmp(argv[nopt], "-fsmax") == 0)
00770 {
00771 nopt++;
00772 if (nopt >= argc) NLfit_error ("need argument after -fsmax ");
00773 *fsmax_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00774 MTEST (*fsmax_filename);
00775 strcpy (*fsmax_filename, argv[nopt]);
00776 nopt++;
00777 continue;
00778 }
00779
00780
00781
00782 if (strcmp(argv[nopt], "-ftmax") == 0)
00783 {
00784 nopt++;
00785 if (nopt >= argc) NLfit_error ("need argument after -ftmax ");
00786 *ftmax_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00787 MTEST (*ftmax_filename);
00788 strcpy (*ftmax_filename, argv[nopt]);
00789 nopt++;
00790 continue;
00791 }
00792
00793
00794
00795 if (strcmp(argv[nopt], "-fpmax") == 0)
00796 {
00797 nopt++;
00798 if (nopt >= argc) NLfit_error ("need argument after -fpmax ");
00799 *fpmax_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00800 MTEST (*fpmax_filename);
00801 strcpy (*fpmax_filename, argv[nopt]);
00802 nopt++;
00803 continue;
00804 }
00805
00806
00807
00808 if (strcmp(argv[nopt], "-farea") == 0)
00809 {
00810 nopt++;
00811 if (nopt >= argc) NLfit_error ("need argument after -farea ");
00812 *farea_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00813 MTEST (*farea_filename);
00814 strcpy (*farea_filename, argv[nopt]);
00815 nopt++;
00816 continue;
00817 }
00818
00819
00820
00821 if (strcmp(argv[nopt], "-fparea") == 0)
00822 {
00823 nopt++;
00824 if (nopt >= argc) NLfit_error ("need argument after -fparea ");
00825 *fparea_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00826 MTEST (*fparea_filename);
00827 strcpy (*fparea_filename, argv[nopt]);
00828 nopt++;
00829 continue;
00830 }
00831
00832
00833
00834 if (strcmp(argv[nopt], "-fscoef") == 0)
00835 {
00836 nopt++;
00837 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -fscoef ");
00838 sscanf (argv[nopt], "%d", &ival);
00839 if ((ival < 0) || (ival >= *p))
00840 NLfit_error ("illegal argument after -fscoef ");
00841 index = ival;
00842 nopt++;
00843
00844 (*fscoef_filename)[index] = malloc (sizeof(char) * MAX_NAME_LENGTH);
00845 MTEST ((*fscoef_filename)[index]);
00846 strcpy ((*fscoef_filename)[index], argv[nopt]);
00847
00848 nopt++;
00849 continue;
00850 }
00851
00852
00853
00854 if (strcmp(argv[nopt], "-fncoef") == 0)
00855 {
00856 nopt++;
00857 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -fncoef ");
00858 sscanf (argv[nopt], "%d", &ival);
00859 if ((ival < 0) || (ival >= *r))
00860 NLfit_error ("illegal argument after -fncoef ");
00861 index = ival;
00862 nopt++;
00863
00864 (*fncoef_filename)[index] = malloc (sizeof(char) * MAX_NAME_LENGTH);
00865 MTEST ((*fncoef_filename)[index]);
00866 strcpy ((*fncoef_filename)[index], argv[nopt]);
00867
00868 nopt++;
00869 continue;
00870 }
00871
00872
00873
00874 if (strcmp(argv[nopt], "-tscoef") == 0)
00875 {
00876 nopt++;
00877 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -tscoef ");
00878 sscanf (argv[nopt], "%d", &ival);
00879 if ((ival < 0) || (ival >= *p))
00880 NLfit_error ("illegal argument after -tscoef ");
00881 index = ival;
00882 nopt++;
00883
00884 (*tscoef_filename)[index] = malloc (sizeof(char) * MAX_NAME_LENGTH);
00885 MTEST ((*tscoef_filename)[index]);
00886 strcpy ((*tscoef_filename)[index], argv[nopt]);
00887
00888 calc_tstats = 1;
00889
00890 nopt++;
00891 continue;
00892 }
00893
00894
00895
00896 if (strcmp(argv[nopt], "-tncoef") == 0)
00897 {
00898 nopt++;
00899 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -tncoef ");
00900 sscanf (argv[nopt], "%d", &ival);
00901 if ((ival < 0) || (ival >= *r))
00902 NLfit_error ("illegal argument after -tncoef ");
00903 index = ival;
00904 nopt++;
00905
00906 (*tncoef_filename)[index] = malloc (sizeof(char) * MAX_NAME_LENGTH);
00907 MTEST ((*tncoef_filename)[index]);
00908 strcpy ((*tncoef_filename)[index], argv[nopt]);
00909
00910 calc_tstats = 1;
00911
00912 nopt++;
00913 continue;
00914 }
00915
00916
00917
00918 if (strcmp(argv[nopt], "-bucket") == 0)
00919 {
00920 nopt++;
00921 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -bucket ");
00922 sscanf (argv[nopt], "%d", &ival);
00923 if ((ival < 0) || (ival > MAX_BRICKS))
00924 NLfit_error ("illegal argument after -bucket ");
00925 nopt++;
00926
00927 option_data->bucket_filename =
00928 malloc (sizeof(char) * MAX_NAME_LENGTH);
00929 MTEST (option_data->bucket_filename);
00930 strcpy (option_data->bucket_filename, argv[nopt]);
00931
00932
00933 if (ival == 0)
00934 nbricks = (*p) + (*r) + 8;
00935 else
00936 nbricks = ival;
00937 option_data->numbricks = nbricks;
00938
00939
00940 option_data->brick_type = malloc (sizeof(int) * nbricks);
00941 MTEST (option_data->brick_type);
00942 option_data->brick_coef = malloc (sizeof(int) * nbricks);
00943 MTEST (option_data->brick_coef);
00944 option_data->brick_label = malloc (sizeof(char *) * nbricks);
00945 MTEST (option_data->brick_label);
00946 for (ibrick = 0; ibrick < nbricks; ibrick++)
00947 {
00948 option_data->brick_type[ibrick] = -1;
00949 option_data->brick_coef[ibrick] = -1;
00950 option_data->brick_label[ibrick] =
00951 malloc (sizeof(char) * MAX_NAME_LENGTH);
00952 MTEST (option_data->brick_label[ibrick]);
00953 }
00954
00955
00956 if (ival == 0)
00957
00958 {
00959 for (ibrick = 0; ibrick < (*r); ibrick++)
00960 {
00961 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00962 option_data->brick_coef[ibrick] = ibrick;
00963 strcpy (option_data->brick_label[ibrick], (*npname)[ibrick]);
00964 }
00965
00966 for (ibrick = (*r); ibrick < (*p) + (*r); ibrick++)
00967 {
00968 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00969 option_data->brick_coef[ibrick] = ibrick;
00970 strcpy (option_data->brick_label[ibrick],
00971 (*spname)[ibrick-(*r)]);
00972 }
00973
00974 ibrick = (*p) + (*r);
00975 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00976 option_data->brick_coef[ibrick] = ibrick;
00977 strcpy (option_data->brick_label[ibrick], "Signal TMax");
00978
00979 ibrick++;
00980 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00981 option_data->brick_coef[ibrick] = ibrick;
00982 strcpy (option_data->brick_label[ibrick], "Signal SMax");
00983
00984 ibrick++;
00985 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00986 option_data->brick_coef[ibrick] = ibrick;
00987 strcpy (option_data->brick_label[ibrick], "Signal % SMax");
00988
00989 ibrick++;
00990 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00991 option_data->brick_coef[ibrick] = ibrick;
00992 strcpy (option_data->brick_label[ibrick], "Signal Area");
00993
00994 ibrick++;
00995 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00996 option_data->brick_coef[ibrick] = ibrick;
00997 strcpy (option_data->brick_label[ibrick], "Signal % Area");
00998
00999 ibrick++;
01000 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01001 option_data->brick_coef[ibrick] = ibrick;
01002 strcpy (option_data->brick_label[ibrick], "Sigma Resid");
01003
01004 ibrick++;
01005 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01006 option_data->brick_coef[ibrick] = ibrick;
01007 strcpy (option_data->brick_label[ibrick], "R^2");
01008
01009 ibrick++;
01010 option_data->brick_type[ibrick] = FUNC_FT_TYPE;
01011 option_data->brick_coef[ibrick] = ibrick;
01012 strcpy (option_data->brick_label[ibrick], "F-stat Regression");
01013
01014 }
01015
01016 nopt++;
01017 continue;
01018 }
01019
01020
01021
01022 if (strcmp(argv[nopt], "-brick") == 0)
01023 {
01024 nopt++;
01025 if (nopt+2 >= argc)
01026 NLfit_error ("need more arguments after -brick ");
01027 sscanf (argv[nopt], "%d", &ibrick);
01028 if ((ibrick < 0) || (ibrick >= option_data->numbricks))
01029 NLfit_error ("illegal argument after -brick ");
01030 nopt++;
01031
01032 if (strcmp(argv[nopt], "scoef") == 0)
01033 {
01034 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01035
01036 nopt++;
01037 sscanf (argv[nopt], "%d", &ival);
01038 if ((ival < 0) || (ival > (*p)))
01039 NLfit_error ("illegal argument after scoef ");
01040 option_data->brick_coef[ibrick] = ival + (*r);
01041
01042 nopt++;
01043 if (nopt >= argc)
01044 NLfit_error ("need more arguments after -brick ");
01045 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01046 }
01047
01048 else if (strcmp(argv[nopt], "ncoef") == 0)
01049 {
01050 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01051
01052 nopt++;
01053 sscanf (argv[nopt], "%d", &ival);
01054 if ((ival < 0) || (ival > (*r)))
01055 NLfit_error ("illegal argument after ncoef ");
01056 option_data->brick_coef[ibrick] = ival;
01057
01058 nopt++;
01059 if (nopt >= argc)
01060 NLfit_error ("need more arguments after -brick ");
01061 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01062 }
01063
01064 else if (strcmp(argv[nopt], "tmax") == 0)
01065 {
01066 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01067 option_data->brick_coef[ibrick] = (*r) + (*p);
01068 nopt++;
01069 if (nopt >= argc)
01070 NLfit_error ("need more arguments after -brick ");
01071 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01072 }
01073
01074 else if (strcmp(argv[nopt], "smax") == 0)
01075 {
01076 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01077 option_data->brick_coef[ibrick] = (*r) + (*p) + 1;
01078 nopt++;
01079 if (nopt >= argc)
01080 NLfit_error ("need more arguments after -brick ");
01081 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01082 }
01083
01084 else if (strcmp(argv[nopt], "psmax") == 0)
01085 {
01086 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01087 option_data->brick_coef[ibrick] = (*r) + (*p) + 2;
01088 nopt++;
01089 if (nopt >= argc)
01090 NLfit_error ("need more arguments after -brick ");
01091 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01092 }
01093
01094 else if (strcmp(argv[nopt], "area") == 0)
01095 {
01096 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01097 option_data->brick_coef[ibrick] = (*r) + (*p) + 3;
01098 nopt++;
01099 if (nopt >= argc)
01100 NLfit_error ("need more arguments after -brick ");
01101 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01102 }
01103
01104 else if (strcmp(argv[nopt], "parea") == 0)
01105 {
01106 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01107 option_data->brick_coef[ibrick] = (*r) + (*p) + 4;
01108 nopt++;
01109 if (nopt >= argc)
01110 NLfit_error ("need more arguments after -brick ");
01111 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01112 }
01113
01114 else if (strcmp(argv[nopt], "resid") == 0)
01115 {
01116 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01117 option_data->brick_coef[ibrick] = (*r) + (*p) + 5;
01118 nopt++;
01119 if (nopt >= argc)
01120 NLfit_error ("need more arguments after -brick ");
01121 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01122 }
01123
01124 else if (strcmp(argv[nopt], "rsqr") == 0)
01125 {
01126 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01127 option_data->brick_coef[ibrick] = (*r) + (*p) + 6;
01128 nopt++;
01129 if (nopt >= argc)
01130 NLfit_error ("need more arguments after -brick ");
01131 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01132 }
01133
01134 else if (strcmp(argv[nopt], "fstat") == 0)
01135 {
01136 option_data->brick_type[ibrick] = FUNC_FT_TYPE;
01137 option_data->brick_coef[ibrick] = (*r) + (*p) + 7;
01138 nopt++;
01139 if (nopt >= argc)
01140 NLfit_error ("need more arguments after -brick ");
01141 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01142 }
01143
01144 else if (strcmp(argv[nopt], "tscoef") == 0)
01145 {
01146 option_data->brick_type[ibrick] = FUNC_TT_TYPE;
01147
01148 nopt++;
01149 sscanf (argv[nopt], "%d", &ival);
01150 if ((ival < 0) || (ival > (*p)))
01151 NLfit_error ("illegal argument after tscoef ");
01152 option_data->brick_coef[ibrick] = ival + (*r);
01153
01154 nopt++;
01155 if (nopt >= argc)
01156 NLfit_error ("need more arguments after -brick ");
01157 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01158
01159 calc_tstats = 1;
01160 }
01161
01162 else if (strcmp(argv[nopt], "tncoef") == 0)
01163 {
01164 option_data->brick_type[ibrick] = FUNC_TT_TYPE;
01165
01166 nopt++;
01167 sscanf (argv[nopt], "%d", &ival);
01168 if ((ival < 0) || (ival > (*r)))
01169 NLfit_error ("illegal argument after tncoef ");
01170 option_data->brick_coef[ibrick] = ival;
01171
01172 nopt++;
01173 if (nopt >= argc)
01174 NLfit_error ("need more arguments after -brick ");
01175 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01176
01177 calc_tstats = 1;
01178 }
01179
01180 else NLfit_error ("unable to interpret options after -brick ");
01181
01182 nopt++;
01183 continue;
01184 }
01185
01186
01187
01188 if (strcmp(argv[nopt], "-sfit") == 0)
01189 {
01190 nopt++;
01191 if (nopt >= argc) NLfit_error ("need argument after -sfit ");
01192 *sfit_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
01193 MTEST (*sfit_filename);
01194 strcpy (*sfit_filename, argv[nopt]);
01195 nopt++;
01196 continue;
01197 }
01198
01199
01200
01201 if (strcmp(argv[nopt], "-snfit") == 0)
01202 {
01203 nopt++;
01204 if (nopt >= argc) NLfit_error ("need argument after -snfit ");
01205 *snfit_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
01206 MTEST (*snfit_filename);
01207 strcpy (*snfit_filename, argv[nopt]);
01208 nopt++;
01209 continue;
01210 }
01211
01212
01213
01214 if( strcmp(argv[nopt],"-jobs") == 0 ){
01215 nopt++ ;
01216 if (nopt >= argc) NLfit_error ("need J parameter after -jobs ");
01217 #ifdef PROC_MAX
01218 proc_numjob = strtol(argv[nopt],NULL,10) ;
01219 if( proc_numjob < 1 ){
01220 fprintf(stderr,"** setting number of processes to 1!\n") ;
01221 proc_numjob = 1 ;
01222 } else if( proc_numjob > PROC_MAX ){
01223 fprintf(stderr,"** setting number of processes to %d!\n",PROC_MAX);
01224 proc_numjob = PROC_MAX ;
01225 }
01226 #else
01227 fprintf(stderr,"** -jobs not supported in this version\n") ;
01228 #endif
01229 nopt++; continue;
01230 }
01231
01232
01233
01234 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
01235 NLfit_error (message);
01236
01237 }
01238
01239
01240
01241 DESTROY_MODEL_ARRAY (model_array) ;
01242
01243 }
01244
01245
01246
01247
01248
01249
01250
01251 void check_one_output_file
01252 (
01253 THD_3dim_dataset * dset_time,
01254 char * filename
01255 )
01256
01257 {
01258 THD_3dim_dataset * new_dset=NULL;
01259 int ierror;
01260 char message[MAX_NAME_LENGTH];
01261
01262
01263
01264 new_dset = EDIT_empty_copy( dset_time ) ;
01265
01266
01267 ierror = EDIT_dset_items( new_dset ,
01268 ADN_prefix , filename ,
01269 ADN_label1 , filename ,
01270 ADN_self_name , filename ,
01271 ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE :
01272 GEN_FUNC_TYPE ,
01273 ADN_none ) ;
01274
01275 if( ierror > 0 )
01276 {
01277 fprintf(stderr,
01278 "*** %d errors in attempting to create output dataset!\n",
01279 ierror ) ;
01280 exit(1) ;
01281 }
01282
01283 if( THD_is_file(new_dset->dblk->diskptr->header_name) )
01284 {
01285 sprintf (message, "Output dataset file %s already exists",
01286 new_dset->dblk->diskptr->header_name ) ;
01287 NLfit_error (message);
01288 }
01289
01290
01291 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01292
01293 }
01294
01295
01296
01297
01298
01299
01300
01301 void check_output_files
01302 (
01303 char * freg_filename,
01304 char * frsqr_filename,
01305 char * fsmax_filename,
01306 char * ftmax_filename,
01307 char * fpmax_filename,
01308 char * farea_filename,
01309 char * fparea_filename,
01310 char ** fncoef_filename,
01311 char ** fscoef_filename,
01312 char ** tncoef_filename,
01313 char ** tscoef_filename,
01314 char * bucket_filename,
01315 char * sfit_filename,
01316 char * snfit_filename,
01317 THD_3dim_dataset * dset_time
01318 )
01319
01320 {
01321 int ip;
01322
01323 if (freg_filename != NULL)
01324 check_one_output_file (dset_time, freg_filename);
01325
01326 if (frsqr_filename != NULL)
01327 check_one_output_file (dset_time, frsqr_filename);
01328
01329 if (fsmax_filename != NULL)
01330 check_one_output_file (dset_time, fsmax_filename);
01331
01332 if (ftmax_filename != NULL)
01333 check_one_output_file (dset_time, ftmax_filename);
01334
01335 if (fpmax_filename != NULL)
01336 check_one_output_file (dset_time, fpmax_filename);
01337
01338 if (farea_filename != NULL)
01339 check_one_output_file (dset_time, farea_filename);
01340
01341 if (fparea_filename != NULL)
01342 check_one_output_file (dset_time, fparea_filename);
01343
01344 for (ip = 0; ip < MAX_PARAMETERS; ip++)
01345 {
01346 if (fncoef_filename[ip] != NULL)
01347 check_one_output_file (dset_time, fncoef_filename[ip]);
01348 if (fscoef_filename[ip] != NULL)
01349 check_one_output_file (dset_time, fscoef_filename[ip]);
01350 if (tncoef_filename[ip] != NULL)
01351 check_one_output_file (dset_time, tncoef_filename[ip]);
01352 if (tscoef_filename[ip] != NULL)
01353 check_one_output_file (dset_time, tscoef_filename[ip]);
01354 }
01355
01356 if (bucket_filename != NULL)
01357 check_one_output_file (dset_time, bucket_filename);
01358
01359
01360 if (sfit_filename != NULL)
01361 check_one_output_file (dset_time, sfit_filename);
01362 if (snfit_filename != NULL)
01363 check_one_output_file (dset_time, snfit_filename);
01364
01365
01366 }
01367
01368
01369
01370
01371
01372
01373
01374 void check_for_valid_inputs
01375 (
01376 int nxyz,
01377 int r,
01378 int p,
01379 float * min_nconstr,
01380 float * max_nconstr,
01381 float * min_sconstr,
01382 float * max_sconstr,
01383 int nrand,
01384 int nbest,
01385
01386 char * freg_filename,
01387 char * frsqr_filename,
01388 char * fsmax_filename,
01389 char * ftmax_filename,
01390 char * fpmax_filename,
01391 char * farea_filename,
01392 char * fparea_filename,
01393 char ** fncoef_filename,
01394 char ** fscoef_filename,
01395 char ** tncoef_filename,
01396 char ** tscoef_filename,
01397 char * bucket_filename,
01398 char * sfit_filename,
01399 char * snfit_filename,
01400 THD_3dim_dataset * dset_time
01401 )
01402
01403 {
01404 int ip;
01405
01406
01407
01408 if (mask_vol != NULL)
01409 if (mask_nvox != nxyz)
01410 NLfit_error ("mask and input dataset bricks don't match in size");
01411
01412
01413
01414 for (ip = 0; ip < r; ip++)
01415 if (min_nconstr[ip] > max_nconstr[ip])
01416 NLfit_error ("Must have minimum constraints <= maximum constraints");
01417 for (ip = 0; ip < p; ip++)
01418 if (min_sconstr[ip] > max_sconstr[ip])
01419 NLfit_error("Must have mininum constraints <= maximum constraints");
01420
01421
01422
01423 if (nrand < nbest)
01424 NLfit_error ("Must have nbest <= nrand");
01425
01426
01427
01428 check_output_files (freg_filename, frsqr_filename,
01429 fsmax_filename, ftmax_filename,
01430 fpmax_filename, farea_filename, fparea_filename,
01431 fncoef_filename, fscoef_filename,
01432 tncoef_filename, tscoef_filename, bucket_filename,
01433 sfit_filename, snfit_filename, dset_time);
01434
01435 }
01436
01437
01438
01439
01440
01441
01442 void zero_fill_volume (float ** fvol, int nxyz)
01443 {
01444 int ixyz;
01445
01446 if( proc_numjob == 1 ){
01447
01448 *fvol = (float *) malloc (sizeof(float) * nxyz); MTEST(*fvol);
01449 for (ixyz = 0; ixyz < nxyz; ixyz++)
01450 (*fvol)[ixyz] = 0.0;
01451
01452 }
01453 #ifdef PROC_MAX
01454 else {
01455
01456
01457 proc_shm_arnum++ ;
01458 proc_shm_arsiz = (int *) realloc( proc_shm_arsiz ,
01459 sizeof(int) *proc_shm_arnum ) ;
01460 proc_shm_ar = (float ***) realloc( proc_shm_ar ,
01461 sizeof(float **)*proc_shm_arnum ) ;
01462 proc_shm_arsiz[proc_shm_arnum-1] = nxyz ;
01463 proc_shm_ar[proc_shm_arnum-1] = fvol ;
01464
01465
01466
01467 }
01468 #endif
01469 }
01470
01471
01472
01473
01474 #ifdef PROC_MAX
01475
01476
01477
01478 #include <signal.h>
01479
01480 void proc_sigfunc(int sig)
01481 {
01482 char * sname ; int ii ;
01483 static volatile int fff=0 ;
01484 if( fff ) _exit(1); else fff=1 ;
01485 switch(sig){
01486 default: sname = "unknown" ; break ;
01487 case SIGHUP: sname = "SIGHUP" ; break ;
01488 case SIGTERM: sname = "SIGTERM" ; break ;
01489 case SIGILL: sname = "SIGILL" ; break ;
01490 case SIGKILL: sname = "SIGKILL" ; break ;
01491 case SIGPIPE: sname = "SIGPIPE" ; break ;
01492 case SIGSEGV: sname = "SIGSEGV" ; break ;
01493 case SIGBUS: sname = "SIGBUS" ; break ;
01494 case SIGINT: sname = "SIGINT" ; break ;
01495 case SIGFPE: sname = "SIGFPE" ; break ;
01496 }
01497 if( proc_shmid > 0 ){
01498 shmctl( proc_shmid , IPC_RMID , NULL ) ; proc_shmid = 0 ;
01499 }
01500 fprintf(stderr,"Fatal Signal %d (%s) received in job #%d\n",
01501 sig,sname,proc_ind) ;
01502 exit(1) ;
01503 }
01504
01505 void proc_atexit( void )
01506 {
01507 if( proc_shmid > 0 )
01508 shmctl( proc_shmid , IPC_RMID , NULL ) ;
01509 }
01510
01511
01512
01513
01514
01515 void proc_finalize_shm_volumes(void)
01516 {
01517 char kstr[32] ; int ii ;
01518
01519 if( proc_shm_arnum == 0 ) return ;
01520
01521 proc_shmsize = 0 ;
01522 for( ii=0 ; ii < proc_shm_arnum ; ii++ )
01523 proc_shmsize += proc_shm_arsiz[ii] ;
01524
01525 proc_shmsize *= sizeof(float) ;
01526
01527
01528
01529 UNIQ_idcode_fill( kstr ) ;
01530 proc_shmid = shm_create( kstr , proc_shmsize ) ;
01531 if( proc_shmid < 0 ){
01532 fprintf(stderr,"\n** Can't create shared memory of size %d!\n"
01533 "** Try re-running without -jobs option!\n" ,
01534 proc_shmsize ) ;
01535
01536
01537
01538 { char *cmd=NULL ;
01539 #if defined(LINUX)
01540 cmd = "cat /proc/sys/kernel/shmmax" ;
01541 #elif defined(SOLARIS)
01542 cmd = "/usr/sbin/sysdef | grep SHMMAX" ;
01543 #elif defined(DARWIN)
01544 cmd = "sysctl -n kern.sysv.shmmax" ;
01545 #endif
01546 if( cmd != NULL ){
01547 FILE *fp = popen( cmd , "r" ) ;
01548 if( fp != NULL ){
01549 unsigned int smax=0 ;
01550 fscanf(fp,"%u",&smax) ; pclose(fp) ;
01551 if( smax > 0 )
01552 fprintf(stderr ,
01553 "\n"
01554 "** POSSIBLY USEFUL ADVICE:\n"
01555 "** Current max shared memory size = %u bytes.\n"
01556 "** For information on how to change this, see\n"
01557 "** http://afni.nimh.nih.gov/afni/parallize.htm\n"
01558 "** and also contact your system administrator.\n"
01559 , smax ) ;
01560 }
01561 }
01562 }
01563 exit(1) ;
01564 }
01565
01566
01567
01568
01569 signal(SIGPIPE,proc_sigfunc) ; signal(SIGSEGV,proc_sigfunc) ;
01570 signal(SIGINT ,proc_sigfunc) ; signal(SIGFPE ,proc_sigfunc) ;
01571 signal(SIGBUS ,proc_sigfunc) ; signal(SIGHUP ,proc_sigfunc) ;
01572 signal(SIGTERM,proc_sigfunc) ; signal(SIGILL ,proc_sigfunc) ;
01573 signal(SIGKILL,proc_sigfunc) ; signal(SIGPIPE,proc_sigfunc) ;
01574 atexit(proc_atexit) ;
01575
01576 fprintf(stderr , "++ Shared memory: %d bytes at id=%d\n" ,
01577 proc_shmsize , proc_shmid ) ;
01578
01579
01580
01581 proc_shmptr = shm_attach( proc_shmid ) ;
01582 if( proc_shmptr == NULL ){
01583 fprintf(stderr,"\n** Can't attach to shared memory!?\n"
01584 "** This is bizarre.\n" ) ;
01585 shmctl( proc_shmid , IPC_RMID , NULL ) ;
01586 exit(1) ;
01587 }
01588
01589
01590
01591 memset( proc_shmptr , 0 , proc_shmsize ) ;
01592
01593
01594
01595 *proc_shm_ar[0] = (float *) proc_shmptr ;
01596 for( ii=1 ; ii < proc_shm_arnum ; ii++ )
01597 *proc_shm_ar[ii] = *proc_shm_ar[ii-1] + proc_shm_arsiz[ii] ;
01598 }
01599
01600
01601
01602
01603
01604 void proc_free( void *ptr )
01605 {
01606 int ii ;
01607
01608 if( ptr == NULL ) return ;
01609 if( proc_shmid == 0 ){ free(ptr); return; }
01610 for( ii=0 ; ii < proc_shm_arnum ; ii++ )
01611 if( ((float *)ptr) == *proc_shm_ar[ii] ) return;
01612 free(ptr); return;
01613 }
01614
01615 #undef free
01616 #define free proc_free
01617
01618 #endif
01619
01620
01621
01622
01623
01624
01625 void initialize_program
01626 (
01627 int argc,
01628 char ** argv,
01629 int * ignore,
01630 char ** nname,
01631 char ** sname,
01632 vfp * nmodel,
01633 vfp * smodel,
01634 int *r,
01635 int *p,
01636 char *** npname,
01637 char *** spname,
01638 float ** min_nconstr,
01639 float ** max_nconstr,
01640 float ** min_sconstr,
01641 float ** max_sconstr,
01642 int * nabs,
01643 int * nrand,
01644 int * nbest,
01645 float * rms_min,
01646 float * fdisp,
01647
01648 char ** input_filename,
01649 char ** tfilename,
01650 char ** freg_filename,
01651 char ** frsqr_filename,
01652 char ** fsmax_filename,
01653 char ** ftmax_filename,
01654 char ** fpmax_filename,
01655 char ** farea_filename,
01656 char ** fparea_filename,
01657 char *** fncoef_filename,
01658 char *** fscoef_filename,
01659 char *** tncoef_filename,
01660 char *** tscoef_filename,
01661 char ** sfit_filename,
01662 char ** snfit_filename,
01663
01664 THD_3dim_dataset ** dset_time,
01665 int * nxyz,
01666 int * ts_length,
01667
01668 float *** x_array,
01669 float ** ts_array,
01670
01671 float ** par_rdcd,
01672 float ** par_full,
01673 float ** tpar_full,
01674
01675 float ** rmsreg_vol,
01676 float ** freg_vol,
01677 float ** rsqr_vol,
01678 float ** smax_vol,
01679 float ** tmax_vol,
01680 float ** pmax_vol,
01681 float ** area_vol,
01682 float ** parea_vol,
01683 float *** ncoef_vol,
01684 float *** scoef_vol,
01685 float *** tncoef_vol,
01686 float *** tscoef_vol,
01687 float *** sfit_vol,
01688 float *** snfit_vol,
01689
01690 NL_options * option_data
01691
01692 )
01693
01694 {
01695 int dimension;
01696 int ip;
01697 int it;
01698 MRI_IMAGE * im, * flim;
01699
01700 int nt;
01701 float * tar;
01702 int ixyz;
01703
01704
01705
01706 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
01707
01708
01709
01710 get_options(argc, argv, ignore, nname, sname, nmodel, smodel,
01711 r, p, npname, spname,
01712 min_nconstr, max_nconstr, min_sconstr, max_sconstr, nabs,
01713 nrand, nbest, rms_min, fdisp, input_filename, tfilename,
01714 freg_filename, frsqr_filename, fsmax_filename,
01715 ftmax_filename, fpmax_filename, farea_filename, fparea_filename,
01716 fncoef_filename, fscoef_filename,
01717 tncoef_filename, tscoef_filename, sfit_filename, snfit_filename,
01718 dset_time, nxyz, ts_length, option_data);
01719
01720
01721
01722 check_for_valid_inputs (*nxyz, *r, *p, *min_nconstr, *max_nconstr,
01723 *min_sconstr, *max_sconstr, *nrand, *nbest,
01724 *freg_filename, *frsqr_filename,
01725 *fsmax_filename, *ftmax_filename,
01726 *fpmax_filename, *farea_filename, *fparea_filename,
01727 *fncoef_filename, *fscoef_filename,
01728 *tncoef_filename, *tscoef_filename,
01729 *sfit_filename, *snfit_filename,
01730 option_data->bucket_filename, *dset_time);
01731
01732
01733
01734 *ts_length = *ts_length - *ignore;
01735 *ts_array = (float *) malloc (sizeof(float) * (*ts_length));
01736 MTEST (*ts_array);
01737
01738
01739
01740 *x_array = (float **) malloc (sizeof(float *) * (*ts_length));
01741 MTEST (*x_array);
01742 for (it = 0; it < (*ts_length); it++)
01743 {
01744 (*x_array)[it] = (float *) malloc (sizeof(float) * 3);
01745 MTEST ((*x_array)[it]);
01746 }
01747
01748
01749
01750 if (*tfilename == NULL)
01751 {
01752 if( inTR && dsTR > 0.0 ){
01753 DELT = dsTR ;
01754 fprintf(stderr,"--- computing with TR = %g\n",DELT) ;
01755 }
01756 for (it = 0; it < (*ts_length); it++)
01757 {
01758 (*x_array)[it][0] = 1.0;
01759 (*x_array)[it][1] = it * DELT;
01760 (*x_array)[it][2] = (it * DELT) * (it * DELT);
01761 }
01762 }
01763 else
01764 {
01765 flim = mri_read_1D(*tfilename) ;
01766 if (flim == NULL)
01767 NLfit_error ("Unable to read time reference file \n");
01768 nt = flim -> nx;
01769 if (nt < (*ts_length))
01770 NLfit_error ("Time reference array is too short");
01771 tar = MRI_FLOAT_PTR(flim) ;
01772 for (it = 0; it < (*ts_length); it++)
01773 {
01774 (*x_array)[it][0] = 1.0;
01775 (*x_array)[it][1] = tar[it] ;
01776 (*x_array)[it][2] = tar[it] * tar[it];
01777 }
01778 mri_free (flim);
01779 }
01780
01781
01782 dimension = (*r) + (*p);
01783 *par_rdcd = (float *) malloc (sizeof(float) * dimension);
01784 MTEST (*par_rdcd);
01785 *par_full = (float *) malloc (sizeof(float) * dimension);
01786 MTEST (*par_full);
01787 *tpar_full = (float *) malloc (sizeof(float) * dimension);
01788 MTEST (*tpar_full);
01789
01790
01791
01792 zero_fill_volume( freg_vol , *nxyz ) ;
01793
01794 if ((*freg_filename != NULL) || (option_data->bucket_filename != NULL))
01795 zero_fill_volume( rmsreg_vol , *nxyz ) ;
01796
01797 if ((*frsqr_filename != NULL) || (option_data->bucket_filename != NULL))
01798 zero_fill_volume( rsqr_vol , *nxyz ) ;
01799
01800 if ((*fsmax_filename != NULL) || (option_data->bucket_filename != NULL))
01801 zero_fill_volume( smax_vol , *nxyz ) ;
01802
01803 if ((*ftmax_filename != NULL) || (option_data->bucket_filename != NULL))
01804 zero_fill_volume( tmax_vol , *nxyz ) ;
01805
01806
01807 if ((*fpmax_filename != NULL) || (option_data->bucket_filename != NULL))
01808 zero_fill_volume( pmax_vol , *nxyz ) ;
01809
01810 if ((*farea_filename != NULL) || (option_data->bucket_filename != NULL))
01811 zero_fill_volume( area_vol , *nxyz ) ;
01812
01813 if ((*fparea_filename != NULL) || (option_data->bucket_filename != NULL))
01814 zero_fill_volume( parea_vol , *nxyz ) ;
01815
01816
01817 *ncoef_vol = (float **) malloc (sizeof(float *) * (*r));
01818 MTEST (*ncoef_vol);
01819 *tncoef_vol = (float **) malloc (sizeof(float *) * (*r));
01820 MTEST (*tncoef_vol);
01821
01822 for (ip = 0; ip < (*r); ip++)
01823 {
01824 if (((*fncoef_filename)[ip] != NULL) || ((*tncoef_filename)[ip] != NULL)
01825 || (option_data->bucket_filename != NULL))
01826 zero_fill_volume( &((*ncoef_vol)[ip]) , *nxyz ) ;
01827 else
01828 (*ncoef_vol)[ip] = NULL;
01829
01830 if (((*tncoef_filename)[ip] != NULL)
01831 || (option_data->bucket_filename != NULL))
01832 zero_fill_volume( &((*tncoef_vol)[ip]) , *nxyz ) ;
01833 else
01834 (*tncoef_vol)[ip] = NULL;
01835 }
01836
01837
01838 *scoef_vol = (float **) malloc (sizeof(float *) * (*p));
01839 MTEST (*scoef_vol);
01840 *tscoef_vol = (float **) malloc (sizeof(float *) * (*p));
01841 MTEST (*tscoef_vol);
01842
01843 for (ip = 0; ip < (*p); ip++)
01844 {
01845 if (((*fscoef_filename)[ip] != NULL) || ((*tscoef_filename)[ip] != NULL)
01846 || (option_data->bucket_filename != NULL))
01847 zero_fill_volume( &((*scoef_vol)[ip]) , *nxyz ) ;
01848 else
01849 (*scoef_vol)[ip] = NULL;
01850
01851 if (((*tscoef_filename)[ip] != NULL)
01852 || (option_data->bucket_filename != NULL))
01853 zero_fill_volume( &((*tscoef_vol)[ip]) , *nxyz ) ;
01854 else
01855 (*tscoef_vol)[ip] = NULL;
01856 }
01857
01858
01859
01860 if (*sfit_filename != NULL)
01861 {
01862 *sfit_vol = (float **) malloc (sizeof(float *) * (*ts_length));
01863 MTEST(*sfit_vol);
01864
01865 for (it = 0; it < *ts_length; it++)
01866 zero_fill_volume( &((*sfit_vol)[it]) , *nxyz ) ;
01867 }
01868
01869
01870 if (*snfit_filename != NULL)
01871 {
01872 *snfit_vol = (float **) malloc (sizeof(float *) * (*ts_length));
01873 MTEST(*snfit_vol);
01874
01875 for (it = 0; it < *ts_length; it++)
01876 zero_fill_volume( &((*snfit_vol)[it]) , *nxyz ) ;
01877 }
01878
01879 #ifdef PROC_MAX
01880 if( proc_numjob > 1 ) proc_finalize_shm_volumes() ;
01881 #endif
01882
01883 }
01884
01885
01886
01887
01888
01889
01890
01891 void read_ts_array
01892 (
01893 THD_3dim_dataset * dset_time,
01894 int iv,
01895 int ts_length,
01896 int ignore,
01897 float * ts_array
01898 )
01899
01900 {
01901 MRI_IMAGE * im;
01902 float * ar;
01903 int it;
01904
01905
01906
01907 im = THD_extract_series (iv, dset_time, 0);
01908
01909
01910
01911 if (im == NULL) NLfit_error ("Unable to extract data from 3d+time dataset");
01912
01913
01914
01915 ar = MRI_FLOAT_PTR (im);
01916 for (it = 0; it < ts_length; it++)
01917 {
01918 ts_array[it] = ar[it+ignore];
01919 }
01920
01921
01922
01923 mri_free (im); im = NULL;
01924
01925 }
01926
01927
01928
01929
01930
01931
01932
01933 void write_ts_array
01934 (
01935 int ts_length,
01936 float * ts_array
01937 )
01938
01939 {
01940 int it;
01941
01942 for (it = 0; it < ts_length; it++)
01943 printf ("%d %f \n", it, ts_array[it]);
01944 }
01945
01946
01947
01948
01949
01950
01951
01952 void save_results
01953 (
01954 int iv,
01955 vfp nmodel,
01956 vfp smodel,
01957 int r,
01958 int p,
01959 int novar,
01960 int ts_length,
01961 float ** x_array,
01962 float * par_full,
01963 float * tpar_full,
01964 float rmsreg,
01965 float freg,
01966 float rsqr,
01967 float smax,
01968 float tmax,
01969 float pmax,
01970 float area,
01971 float parea,
01972
01973 float * rmsreg_vol,
01974 float * freg_vol,
01975 float * rsqr_vol,
01976 float * smax_vol,
01977 float * tmax_vol,
01978 float * pmax_vol,
01979 float * area_vol,
01980 float * parea_vol,
01981 float ** ncoef_vol,
01982 float ** scoef_vol,
01983 float ** tncoef_vol,
01984 float ** tscoef_vol,
01985 float ** sfit_vol,
01986 float ** snfit_vol
01987 )
01988
01989 {
01990 int ip;
01991 int it;
01992 float * s_array;
01993 float * n_array;
01994
01995
01996
01997 if (freg_vol != NULL) freg_vol[iv] = freg;
01998 if (rmsreg_vol != NULL) rmsreg_vol[iv] = rmsreg;
01999 if (rsqr_vol != NULL) rsqr_vol[iv] = rsqr;
02000
02001
02002 if (smax_vol != NULL) smax_vol[iv] = smax;
02003 if (tmax_vol != NULL) tmax_vol[iv] = tmax;
02004
02005
02006 if (pmax_vol != NULL) pmax_vol[iv] = pmax;
02007 if (area_vol != NULL) area_vol[iv] = area;
02008 if (parea_vol != NULL) parea_vol[iv] = parea;
02009
02010
02011 for (ip = 0; ip < r; ip++)
02012 {
02013 if (ncoef_vol[ip] != NULL) ncoef_vol[ip][iv] = par_full[ip];
02014 if (tncoef_vol[ip] != NULL) tncoef_vol[ip][iv] = tpar_full[ip];
02015 }
02016
02017
02018 for (ip = 0; ip < p; ip++)
02019 {
02020 if (scoef_vol[ip] != NULL) scoef_vol[ip][iv] = par_full[ip+r];
02021 if (tscoef_vol[ip] != NULL) tscoef_vol[ip][iv] = tpar_full[ip+r];
02022 }
02023
02024
02025
02026 if (sfit_vol != NULL)
02027 {
02028 if (novar)
02029 {
02030 for (it = 0; it < ts_length; it++)
02031 sfit_vol[it][iv] = 0.0;
02032 }
02033 else
02034 {
02035 s_array = (float *) malloc (sizeof(float) * (ts_length));
02036 MTEST (s_array);
02037
02038 smodel (par_full + r, ts_length, x_array, s_array);
02039
02040 for (it = 0; it < ts_length; it++)
02041 sfit_vol[it][iv] = s_array[it];
02042
02043 free (s_array); s_array = NULL;
02044 }
02045 }
02046
02047
02048
02049 if (snfit_vol != NULL)
02050 {
02051 n_array = (float *) malloc (sizeof(float) * (ts_length));
02052 MTEST (n_array);
02053 nmodel (par_full, ts_length, x_array, n_array);
02054
02055 for (it = 0; it < ts_length; it++)
02056 {
02057 snfit_vol[it][iv] = n_array[it];
02058 }
02059
02060 free (n_array); n_array = NULL;
02061
02062
02063 if (! novar)
02064 {
02065 s_array = (float *) malloc (sizeof(float) * (ts_length));
02066 MTEST (s_array);
02067 smodel (par_full + r, ts_length, x_array, s_array);
02068
02069 for (it = 0; it < ts_length; it++)
02070 {
02071 snfit_vol[it][iv] += s_array[it];
02072 }
02073
02074 free (s_array); s_array = NULL;
02075 }
02076 }
02077
02078 }
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092 float EDIT_coerce_autoscale_new( int nxyz ,
02093 int itype,void *ivol , int otype,void *ovol )
02094 {
02095 float fac=0.0 , top ;
02096
02097 if( MRI_IS_INT_TYPE(otype) ){
02098 top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
02099 if (top == 0.0) fac = 0.0;
02100 else fac = MRI_TYPE_maxval[otype]/top;
02101 }
02102
02103 EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
02104 return ( fac );
02105 }
02106
02107
02108
02109
02110
02111
02112
02113
02114
02115
02116
02117
02118
02119
02120
02121
02122
02123 void write_afni_data (char * input_filename, int nxyz, char * filename,
02124 float * ffim, float * ftr, int numdof, int dendof)
02125 {
02126 int ii;
02127 THD_3dim_dataset * dset=NULL;
02128 THD_3dim_dataset * new_dset=NULL;
02129 int iv;
02130 int ierror;
02131 int ibuf[32];
02132 float fbuf[MAX_STAT_AUX];
02133 float fimfac;
02134 int output_datum;
02135 short * tsp;
02136 void * vdif;
02137 int func_type;
02138 float top, func_scale_short;
02139 char label[80];
02140
02141
02142
02143 dset = THD_open_one_dataset (input_filename) ;
02144 if( ! ISVALID_3DIM_DATASET(dset) ){
02145 fprintf(stderr,"*** Unable to open dataset file %s\n",
02146 input_filename);
02147 exit(1) ;
02148 }
02149
02150
02151 new_dset = EDIT_empty_copy( dset ) ;
02152
02153
02154
02155 tross_Copy_History( dset , new_dset ) ;
02156 sprintf (label, "Output prefix: %s", filename);
02157 if( commandline != NULL )
02158 tross_multi_Append_History( new_dset , commandline,label,NULL ) ;
02159 else
02160 tross_Append_History ( new_dset, label);
02161
02162
02163 iv = DSET_PRINCIPAL_VALUE(dset) ;
02164 output_datum = DSET_BRICK_TYPE(dset,iv) ;
02165 if( output_datum == MRI_byte ) output_datum = MRI_short ;
02166
02167
02168 ibuf[0] = output_datum ; ibuf[1] = MRI_short ;
02169
02170 if (dendof == 0) func_type = FUNC_TT_TYPE;
02171 else func_type = FUNC_FT_TYPE;
02172
02173 ierror = EDIT_dset_items( new_dset ,
02174 ADN_prefix , filename ,
02175 ADN_label1 , filename ,
02176 ADN_self_name , filename ,
02177 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE :
02178 GEN_FUNC_TYPE ,
02179 ADN_func_type , func_type ,
02180 ADN_nvals , FUNC_nvals[func_type] ,
02181 ADN_datum_array , ibuf ,
02182 ADN_ntt , 0 ,
02183 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
02184 ADN_none ) ;
02185
02186 if( ierror > 0 ){
02187 fprintf(stderr,
02188 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
02189 exit(1) ;
02190 }
02191
02192 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
02193 fprintf(stderr,
02194 "*** Output dataset file %s already exists--cannot continue!\a\n",
02195 new_dset->dblk->diskptr->header_name ) ;
02196 exit(1) ;
02197 }
02198
02199
02200 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
02201
02202
02203
02204 vdif = (void *) malloc( mri_datum_size(output_datum) * nxyz );
02205 if (vdif == NULL) NLfit_error ("Unable to allocate memory for vdif");
02206 tsp = (short *) malloc( sizeof(short) * nxyz );
02207 if (tsp == NULL) NLfit_error ("Unable to allocate memory for tsp");
02208
02209
02210 mri_fix_data_pointer (vdif, DSET_BRICK(new_dset,0));
02211 mri_fix_data_pointer (tsp, DSET_BRICK(new_dset,1));
02212
02213
02214
02215 fimfac = EDIT_coerce_autoscale_new (nxyz,
02216 MRI_float, ffim,
02217 output_datum, vdif);
02218 if (fimfac != 0.0) fimfac = 1.0 / fimfac;
02219
02220 #define TOP_SS 32700
02221
02222 if (dendof == 0)
02223 {
02224 top = TOP_SS/FUNC_TT_SCALE_SHORT;
02225 func_scale_short = FUNC_TT_SCALE_SHORT;
02226 }
02227 else
02228 {
02229 top = TOP_SS/FUNC_FT_SCALE_SHORT;
02230 func_scale_short = FUNC_FT_SCALE_SHORT;
02231 }
02232
02233 for (ii = 0; ii < nxyz; ii++)
02234 {
02235 if (ftr[ii] > top)
02236 tsp[ii] = TOP_SS;
02237 else if (ftr[ii] < -top)
02238 tsp[ii] = -TOP_SS;
02239 else if (ftr[ii] >= 0.0)
02240 tsp[ii] = (short) (func_scale_short * ftr[ii] + 0.5);
02241 else
02242 tsp[ii] = (short) (func_scale_short * ftr[ii] - 0.5);
02243
02244 }
02245
02246
02247
02248 printf("++ Writing combined dataset into %s\n",DSET_BRIKNAME(new_dset)) ;
02249
02250 fbuf[0] = numdof;
02251 fbuf[1] = dendof;
02252 for( ii=2 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
02253 (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
02254
02255 fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
02256 fbuf[1] = 1.0 / func_scale_short ;
02257 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
02258
02259 THD_load_statistics( new_dset ) ;
02260 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
02261
02262
02263 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
02264
02265 }
02266
02267
02268
02269
02270
02271
02272
02273 void write_bucket_data
02274 (
02275 int q,
02276 int p,
02277 int numdof,
02278 int dendof,
02279 int nxyz,
02280 int n,
02281
02282 float * rmsreg_vol,
02283 float * freg_vol,
02284 float * rsqr_vol,
02285 float * smax_vol,
02286 float * tmax_vol,
02287 float * pmax_vol,
02288 float * area_vol,
02289 float * parea_vol,
02290 float ** ncoef_vol,
02291 float ** scoef_vol,
02292 float ** tncoef_vol,
02293 float ** tscoef_vol,
02294
02295 char * input_filename,
02296
02297 NL_options * option_data
02298 )
02299
02300 {
02301 const float EPSILON = 1.0e-10;
02302
02303 THD_3dim_dataset * old_dset = NULL;
02304 THD_3dim_dataset * new_dset = NULL;
02305 char * output_prefix;
02306 char * output_session;
02307 int nbricks, ib;
02308 short ** bar = NULL;
02309 float factor;
02310 int brick_type;
02311 int brick_coef;
02312 char * brick_label;
02313 int ierror;
02314 float * volume;
02315 int dimension;
02316 char label[80];
02317
02318
02319
02320 nbricks = option_data->numbricks;
02321 output_prefix = option_data->bucket_filename;
02322 output_session = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
02323 strcpy (output_session, "./");
02324 dimension = p + q;
02325
02326
02327
02328 bar = (short **) malloc (sizeof(short *) * nbricks);
02329 MTEST (bar);
02330
02331
02332
02333 old_dset = THD_open_one_dataset (input_filename);
02334
02335
02336
02337 new_dset = EDIT_empty_copy (old_dset);
02338
02339
02340
02341 tross_Copy_History( old_dset , new_dset ) ;
02342 if( commandline != NULL )
02343 tross_Append_History( new_dset , commandline ) ;
02344 sprintf (label, "Output prefix: %s", output_prefix);
02345 tross_Append_History ( new_dset, label);
02346
02347
02348
02349
02350 ierror = EDIT_dset_items (new_dset,
02351 ADN_prefix, output_prefix,
02352 ADN_directory_name, output_session,
02353 ADN_type, HEAD_FUNC_TYPE,
02354 ADN_func_type, FUNC_BUCK_TYPE,
02355 ADN_ntt, 0,
02356 ADN_nvals, nbricks,
02357 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
02358 ADN_none ) ;
02359
02360 if( ierror > 0 )
02361 {
02362 fprintf(stderr,
02363 "*** %d errors in attempting to create output dataset!\n",
02364 ierror);
02365 exit(1);
02366 }
02367
02368 if (THD_is_file(DSET_HEADNAME(new_dset)))
02369 {
02370 fprintf(stderr,
02371 "*** Output dataset file %s already exists--cannot continue!\n",
02372 DSET_HEADNAME(new_dset));
02373 exit(1);
02374 }
02375
02376
02377
02378 THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ;
02379
02380
02381
02382
02383
02384 for (ib = 0; ib < nbricks; ib++)
02385 {
02386
02387 brick_type = option_data->brick_type[ib];
02388 brick_coef = option_data->brick_coef[ib];
02389 brick_label = option_data->brick_label[ib];
02390
02391 if (brick_type == FUNC_FIM_TYPE)
02392 {
02393 if (brick_coef < q)
02394 volume = ncoef_vol[brick_coef];
02395 else if (brick_coef < p+q)
02396 volume = scoef_vol[brick_coef-q];
02397 else if (brick_coef == p+q)
02398 volume = tmax_vol;
02399 else if (brick_coef == p+q+1)
02400 volume = smax_vol;
02401 else if (brick_coef == p+q+2)
02402 volume = pmax_vol;
02403 else if (brick_coef == p+q+3)
02404 volume = area_vol;
02405 else if (brick_coef == p+q+4)
02406 volume = parea_vol;
02407 else if (brick_coef == p+q+5)
02408 volume = rmsreg_vol;
02409 else if (brick_coef == p+q+6)
02410 volume = rsqr_vol;
02411 }
02412 else if (brick_type == FUNC_TT_TYPE)
02413 {
02414 if (brick_coef < q)
02415 volume = tncoef_vol[brick_coef];
02416 else if (brick_coef < p+q)
02417 volume = tscoef_vol[brick_coef-q];
02418 EDIT_BRICK_TO_FITT (new_dset, ib, dendof);
02419 }
02420 else if (brick_type == FUNC_FT_TYPE)
02421 {
02422 volume = freg_vol;
02423 EDIT_BRICK_TO_FIFT (new_dset, ib, numdof, dendof);
02424 }
02425
02426
02427 bar[ib] = (short *) malloc (sizeof(short) * nxyz);
02428 MTEST (bar[ib]);
02429 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
02430 MRI_short, bar[ib]);
02431
02432 if (factor < EPSILON) factor = 0.0;
02433 else factor = 1.0 / factor;
02434
02435
02436 EDIT_BRICK_LABEL (new_dset, ib, brick_label);
02437 EDIT_BRICK_FACTOR (new_dset, ib, factor);
02438
02439
02440
02441 EDIT_substitute_brick (new_dset, ib, MRI_short, bar[ib]);
02442
02443 }
02444
02445
02446
02447 THD_load_statistics (new_dset);
02448 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
02449 fprintf(stderr,"++ Wrote bucket dataset: %s\n",DSET_BRIKNAME(new_dset)) ;
02450
02451
02452 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
02453
02454 }
02455
02456
02457
02458
02459
02460
02461
02462
02463 void write_3dtime
02464 (
02465 char * input_filename,
02466 int ts_length,
02467 float ** vol_array,
02468 char * output_filename
02469 )
02470
02471 {
02472 const float EPSILON = 1.0e-10;
02473
02474 THD_3dim_dataset * dset = NULL;
02475 THD_3dim_dataset * new_dset = NULL;
02476 int ib;
02477 int ierror;
02478 int nxyz;
02479 float factor;
02480 short ** bar = NULL;
02481 float * fbuf;
02482 float * volume;
02483 char label[80];
02484
02485
02486
02487 dset = THD_open_one_dataset (input_filename);
02488 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
02489
02490
02491
02492 bar = (short **) malloc (sizeof(short *) * ts_length); MTEST (bar);
02493 fbuf = (float *) malloc (sizeof(float) * ts_length); MTEST (fbuf);
02494 for (ib = 0; ib < ts_length; ib++) fbuf[ib] = 0.0;
02495
02496
02497
02498 new_dset = EDIT_empty_copy (dset);
02499
02500
02501
02502 tross_Copy_History( dset , new_dset ) ;
02503 if( commandline != NULL )
02504 tross_Append_History( new_dset , commandline ) ;
02505 sprintf (label, "Output prefix: %s", output_filename);
02506 tross_Append_History ( new_dset, label);
02507
02508
02509
02510 THD_delete_3dim_dataset (dset, False); dset = NULL ;
02511
02512
02513 ierror = EDIT_dset_items (new_dset,
02514 ADN_prefix, output_filename,
02515 ADN_label1, output_filename,
02516 ADN_self_name, output_filename,
02517 ADN_malloc_type, DATABLOCK_MEM_MALLOC,
02518 ADN_datum_all, MRI_short,
02519 ADN_nvals, ts_length,
02520 ADN_ntt, ts_length,
02521 ADN_none);
02522
02523 if( ierror > 0 ){
02524 fprintf(stderr,
02525 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
02526 exit(1) ;
02527 }
02528
02529 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
02530 fprintf(stderr,
02531 "*** Output dataset file %s already exists--cannot continue!\a\n",
02532 new_dset->dblk->diskptr->header_name ) ;
02533 exit(1) ;
02534 }
02535
02536
02537
02538 for (ib = 0; ib < ts_length; ib++)
02539 {
02540
02541
02542 volume = vol_array[ib];
02543
02544
02545 bar[ib] = (short *) malloc (sizeof(short) * nxyz);
02546 MTEST (bar[ib]);
02547
02548
02549 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
02550 MRI_short, bar[ib]);
02551 if (factor < EPSILON) factor = 0.0;
02552 else factor = 1.0 / factor;
02553 fbuf[ib] = factor;
02554
02555
02556 mri_fix_data_pointer (bar[ib], DSET_BRICK(new_dset,ib));
02557 }
02558
02559
02560
02561
02562 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
02563
02564 THD_load_statistics (new_dset);
02565 THD_write_3dim_dataset (NULL, NULL, new_dset, True);
02566 fprintf (stderr,"++ Wrote 3D+time dataset %s\n",DSET_BRIKNAME(new_dset)) ;
02567
02568
02569
02570 THD_delete_3dim_dataset (new_dset, False); new_dset = NULL ;
02571 free (fbuf); fbuf = NULL;
02572
02573 }
02574
02575
02576
02577
02578
02579
02580
02581 void output_results
02582 (
02583 int r,
02584 int p,
02585 float * min_nconstr,
02586 float * max_nconstr,
02587 float * min_sconstr,
02588 float * max_sconstr,
02589 int nxyz,
02590 int ts_length,
02591
02592 float * rmsreg_vol,
02593 float * freg_vol,
02594 float * rsqr_vol,
02595 float * smax_vol,
02596 float * tmax_vol,
02597 float * pmax_vol,
02598 float * area_vol,
02599 float * parea_vol,
02600 float ** ncoef_vol,
02601 float ** scoef_vol,
02602 float ** tncoef_vol,
02603 float ** tscoef_vol,
02604 float ** sfit_vol,
02605 float ** snfit_vol,
02606
02607 char * input_filename,
02608 char * freg_filename,
02609 char * frsqr_filename,
02610 char * fsmax_filename,
02611 char * ftmax_filename,
02612 char * fpmax_filename,
02613 char * farea_filename,
02614 char * fparea_filename,
02615 char ** fncoef_filename,
02616 char ** fscoef_filename,
02617 char ** tncoef_filename,
02618 char ** tscoef_filename,
02619 char * sfit_filename,
02620 char * snfit_filename,
02621
02622 NL_options * option_data
02623
02624 )
02625
02626 {
02627 int ip;
02628 int dimension;
02629 int numdof, dendof;
02630
02631
02632
02633 dimension = r + p;
02634 numdof = p;
02635 dendof = ts_length - dimension;
02636
02637
02638
02639 for (ip = 0; ip < r; ip++)
02640 if (min_nconstr[ip] == max_nconstr[ip])
02641 dendof++;
02642
02643 for (ip = 0; ip < p; ip++)
02644 if (min_sconstr[ip] == max_sconstr[ip])
02645 {
02646 numdof--;
02647 dendof++;
02648 }
02649
02650
02651
02652 if (option_data->numbricks > 0)
02653 write_bucket_data (r, p, numdof, dendof, nxyz, ts_length, rmsreg_vol,
02654 freg_vol, rsqr_vol, smax_vol, tmax_vol, pmax_vol,
02655 area_vol, parea_vol, ncoef_vol, scoef_vol,
02656 tncoef_vol, tscoef_vol, input_filename, option_data);
02657
02658
02659
02660 if (freg_filename != NULL)
02661 write_afni_data (input_filename, nxyz, freg_filename,
02662 rmsreg_vol, freg_vol, numdof, dendof);
02663
02664
02665
02666 if (frsqr_filename != NULL)
02667 write_afni_data (input_filename, nxyz, frsqr_filename,
02668 rsqr_vol, freg_vol, numdof, dendof);
02669
02670
02671
02672 if (fsmax_filename != NULL)
02673 write_afni_data (input_filename, nxyz, fsmax_filename,
02674 smax_vol, freg_vol, numdof, dendof);
02675
02676
02677
02678 if (ftmax_filename != NULL)
02679 write_afni_data (input_filename, nxyz, ftmax_filename,
02680 tmax_vol, freg_vol, numdof, dendof);
02681
02682
02683
02684 if (fpmax_filename != NULL)
02685 write_afni_data (input_filename, nxyz, fpmax_filename,
02686 pmax_vol, freg_vol, numdof, dendof);
02687
02688
02689
02690 if (farea_filename != NULL)
02691 write_afni_data (input_filename, nxyz, farea_filename,
02692 area_vol, freg_vol, numdof, dendof);
02693
02694
02695
02696 if (fparea_filename != NULL)
02697 write_afni_data (input_filename, nxyz, fparea_filename,
02698 parea_vol, freg_vol, numdof, dendof);
02699
02700
02701
02702 for (ip = 0; ip < r; ip++)
02703 {
02704 if (tncoef_filename[ip] != NULL)
02705 write_afni_data (input_filename, nxyz, tncoef_filename[ip],
02706 ncoef_vol[ip], tncoef_vol[ip], dendof, 0);
02707
02708 if (fncoef_filename[ip] != NULL)
02709 write_afni_data (input_filename, nxyz, fncoef_filename[ip],
02710 ncoef_vol[ip], freg_vol, numdof, dendof);
02711 }
02712
02713
02714
02715 for (ip = 0; ip < p; ip++)
02716 {
02717 if (tscoef_filename[ip] != NULL)
02718 write_afni_data (input_filename, nxyz, tscoef_filename[ip],
02719 scoef_vol[ip], tscoef_vol[ip], dendof, 0);
02720
02721 if (fscoef_filename[ip] != NULL)
02722 write_afni_data (input_filename, nxyz, fscoef_filename[ip],
02723 scoef_vol[ip], freg_vol, numdof, dendof);
02724 }
02725
02726
02727
02728 if (sfit_filename != NULL)
02729 {
02730 write_3dtime (input_filename, ts_length, sfit_vol, sfit_filename);
02731 }
02732
02733
02734
02735 if (snfit_filename != NULL)
02736 {
02737 write_3dtime (input_filename, ts_length, snfit_vol, snfit_filename);
02738 }
02739
02740 }
02741
02742
02743
02744
02745
02746
02747
02748 void terminate_program
02749 (
02750 int r,
02751 int p,
02752 int ts_length,
02753 float *** x_array,
02754 float ** ts_array,
02755 char ** nname,
02756 char *** npname,
02757 float ** par_rdcd,
02758 float ** min_nconstr,
02759 float ** max_nconstr,
02760 char ** sname,
02761 char *** spname,
02762 float ** par_full,
02763 float ** tpar_full,
02764 float ** min_sconstr,
02765 float ** max_sconstr,
02766
02767 float ** rmsreg_vol,
02768 float ** freg_vol,
02769 float ** rsqr_vol,
02770 float ** smax_vol,
02771 float ** tmax_vol,
02772 float ** pmax_vol,
02773 float ** area_vol,
02774 float ** parea_vol,
02775 float *** ncoef_vol,
02776 float *** scoef_vol,
02777 float *** tncoef_vol,
02778 float *** tscoef_vol,
02779 float *** sfit_vol,
02780 float *** snfit_vol,
02781
02782 char ** input_filename,
02783 char ** freg_filename,
02784 char ** frsqr_filename,
02785 char ** fsmax_filename,
02786 char ** ftmax_filename,
02787 char ** fpmax_filename,
02788 char ** farea_filename,
02789 char ** fparea_filename,
02790 char *** fncoef_filename,
02791 char *** fscoef_filename,
02792 char *** tncoef_filename,
02793 char *** tscoef_filename,
02794 char ** sfit_filename,
02795 char ** snfit_filename
02796 )
02797
02798 {
02799 int ip;
02800 int it;
02801
02802
02803
02804 if (*nname != NULL)
02805 { free (*nname); *nname = NULL; }
02806 if (*sname != NULL)
02807 { free (*sname); *sname = NULL; }
02808 for (ip = 0; ip < MAX_PARAMETERS; ip++)
02809 {
02810 if ((*npname)[ip] != NULL)
02811 { free ((*npname)[ip]); (*npname)[ip] = NULL; }
02812 if ((*spname)[ip] != NULL)
02813 { free ((*spname)[ip]); (*spname)[ip] = NULL; }
02814 }
02815
02816 if (*npname != NULL)
02817 { free (*npname); *npname = NULL; }
02818 if (*spname != NULL)
02819 { free (*spname); *spname = NULL; }
02820
02821
02822
02823 if (*min_nconstr != NULL) { free (*min_nconstr); *min_nconstr = NULL; }
02824 if (*max_nconstr != NULL) { free (*max_nconstr); *max_nconstr = NULL; }
02825 if (*min_sconstr != NULL) { free (*min_sconstr); *min_sconstr = NULL; }
02826 if (*max_sconstr != NULL) { free (*max_sconstr); *max_sconstr = NULL; }
02827
02828
02829
02830 if (*input_filename != NULL)
02831 { free (*input_filename); *input_filename = NULL; }
02832 if (*freg_filename != NULL)
02833 { free (*freg_filename); *freg_filename = NULL; }
02834 if (*frsqr_filename != NULL)
02835 { free (*frsqr_filename); *frsqr_filename = NULL; }
02836 if (*fsmax_filename != NULL)
02837 { free (*fsmax_filename); *fsmax_filename = NULL; }
02838 if (*ftmax_filename != NULL)
02839 { free (*ftmax_filename); *ftmax_filename = NULL; }
02840 if (*fpmax_filename != NULL)
02841 { free (*fpmax_filename); *fpmax_filename = NULL; }
02842 if (*farea_filename != NULL)
02843 { free (*farea_filename); *farea_filename = NULL; }
02844 if (*fparea_filename != NULL)
02845 { free (*fparea_filename); *fparea_filename = NULL; }
02846
02847 for (ip = 0; ip < MAX_PARAMETERS; ip++)
02848 {
02849 if ((*fncoef_filename)[ip] != NULL)
02850 { free ((*fncoef_filename)[ip]); (*fncoef_filename)[ip] = NULL; }
02851 if ((*fscoef_filename)[ip] != NULL)
02852 { free ((*fscoef_filename)[ip]); (*fscoef_filename)[ip] = NULL; }
02853 if ((*tncoef_filename)[ip] != NULL)
02854 { free ((*tncoef_filename)[ip]); (*tncoef_filename)[ip] = NULL; }
02855 if ((*tscoef_filename)[ip] != NULL)
02856 { free ((*tscoef_filename)[ip]); (*tscoef_filename)[ip] = NULL; }
02857 }
02858
02859 if (*fncoef_filename != NULL)
02860 { free (*fncoef_filename); *fncoef_filename = NULL; }
02861 if (*fscoef_filename != NULL)
02862 { free (*fscoef_filename); *fscoef_filename = NULL; }
02863 if (*tncoef_filename != NULL)
02864 { free (*tncoef_filename); *tncoef_filename = NULL; }
02865 if (*tscoef_filename != NULL)
02866 { free (*tscoef_filename); *tscoef_filename = NULL; }
02867
02868 if (*sfit_filename != NULL)
02869 { free (*sfit_filename); *sfit_filename = NULL; }
02870 if (*snfit_filename != NULL)
02871 { free (*snfit_filename); *snfit_filename = NULL; }
02872
02873
02874
02875 if (*ts_array != NULL) { free (*ts_array); *ts_array = NULL; }
02876
02877
02878
02879 for (it = 0; it < ts_length; it++)
02880 if ((*x_array)[it] != NULL)
02881 { free ((*x_array)[it]); (*x_array)[it] = NULL; }
02882 if (*x_array != NULL) { free (*x_array); *x_array = NULL; }
02883
02884
02885
02886 if (*par_rdcd != NULL) { free (*par_rdcd); *par_rdcd = NULL; }
02887 if (*par_full != NULL) { free (*par_full); *par_full = NULL; }
02888 if (*tpar_full != NULL) { free (*tpar_full); *tpar_full = NULL; }
02889
02890
02891
02892 if (*rmsreg_vol != NULL) { free (*rmsreg_vol); *rmsreg_vol = NULL; }
02893 if (*freg_vol != NULL) { free (*freg_vol); *freg_vol = NULL; }
02894 if (*rsqr_vol != NULL) { free (*rsqr_vol); *rsqr_vol = NULL; }
02895 if (*smax_vol != NULL) { free (*smax_vol); *smax_vol = NULL; }
02896 if (*tmax_vol != NULL) { free (*tmax_vol); *tmax_vol = NULL; }
02897 if (*pmax_vol != NULL) { free (*pmax_vol); *pmax_vol = NULL; }
02898 if (*area_vol != NULL) { free (*area_vol); *area_vol = NULL; }
02899 if (*parea_vol != NULL) { free (*parea_vol); *parea_vol = NULL; }
02900
02901 if (*ncoef_vol != NULL)
02902 {
02903 for (ip = 0; ip < r; ip++)
02904 if ((*ncoef_vol)[ip] != NULL)
02905 { free ((*ncoef_vol)[ip]); (*ncoef_vol)[ip] = NULL; }
02906 free (*ncoef_vol); *ncoef_vol = NULL;
02907 }
02908
02909 if (*tncoef_vol != NULL)
02910 {
02911 for (ip = 0; ip < r; ip++)
02912 if ((*tncoef_vol)[ip] != NULL)
02913 { free ((*tncoef_vol)[ip]); (*tncoef_vol)[ip] = NULL; }
02914 free (*tncoef_vol); *tncoef_vol = NULL;
02915 }
02916
02917 if (*scoef_vol != NULL)
02918 {
02919 for (ip = 0; ip < p; ip++)
02920 if ((*scoef_vol)[ip] != NULL)
02921 { free ((*scoef_vol)[ip]); (*scoef_vol)[ip] = NULL; }
02922 free (*scoef_vol); *scoef_vol = NULL;
02923 }
02924
02925 if (*tscoef_vol != NULL)
02926 {
02927 for (ip = 0; ip < p; ip++)
02928 if ((*tscoef_vol)[ip] != NULL)
02929 { free ((*tscoef_vol)[ip]); (*tscoef_vol)[ip] = NULL; }
02930 free (*tscoef_vol); *tscoef_vol = NULL;
02931 }
02932
02933 if (*sfit_vol != NULL)
02934 {
02935 for (it = 0; it < ts_length; it++)
02936 if ((*sfit_vol)[it] != NULL)
02937 { free ((*sfit_vol)[it]); (*sfit_vol)[it] = NULL; }
02938 free (*sfit_vol); *sfit_vol = NULL;
02939 }
02940
02941 if (*snfit_vol != NULL)
02942 {
02943 for (it = 0; it < ts_length; it++)
02944 if ((*snfit_vol)[it] != NULL)
02945 { free ((*snfit_vol)[it]); (*snfit_vol)[it] = NULL; }
02946 free (*snfit_vol); *snfit_vol = NULL;
02947 }
02948
02949 }
02950
02951
02952
02953
02954 int main
02955 (
02956 int argc,
02957 char ** argv
02958 )
02959
02960 {
02961
02962 NL_options option_data;
02963 int nabs;
02964 int nrand;
02965 int nbest;
02966 float rms_min;
02967 float fdisp;
02968
02969
02970 THD_3dim_dataset * dset_time = NULL;
02971 int ts_length;
02972 int ignore;
02973 float ** x_array = NULL;
02974 float * ts_array = NULL;
02975 int nxyz;
02976 int iv;
02977
02978
02979 char * nname = NULL;
02980 vfp nmodel;
02981 int r;
02982 char ** npname = NULL;
02983 float * par_rdcd = NULL;
02984 float sse_rdcd;
02985 float * min_nconstr = NULL;
02986 float * max_nconstr = NULL;
02987
02988
02989 char * sname = NULL;
02990 vfp smodel;
02991 int p;
02992 char ** spname = NULL;
02993 float * par_full = NULL;
02994 float sse_full;
02995 float * tpar_full = NULL;
02996 float freg;
02997 float rmsreg;
02998 float rsqr;
02999 float smax;
03000 float tmax;
03001 float pmax;
03002 float area;
03003 float parea;
03004 float * min_sconstr = NULL;
03005 float * max_sconstr = NULL;
03006
03007
03008 float * rmsreg_vol = NULL;
03009 float * freg_vol = NULL;
03010 float * rsqr_vol = NULL;
03011 float * smax_vol = NULL;
03012 float * tmax_vol = NULL;
03013 float * pmax_vol = NULL;
03014 float * area_vol = NULL;
03015 float * parea_vol = NULL;
03016 float ** ncoef_vol = NULL;
03017 float ** scoef_vol = NULL;
03018 float ** tncoef_vol = NULL;
03019 float ** tscoef_vol = NULL;
03020 float ** sfit_vol = NULL;
03021 float ** snfit_vol = NULL;
03022
03023
03024 char * input_filename = NULL;
03025 char * tfilename = NULL;
03026 char * freg_filename = NULL;
03027 char * frsqr_filename= NULL;
03028 char * fsmax_filename = NULL;
03029 char * ftmax_filename = NULL;
03030 char * fpmax_filename = NULL;
03031 char * farea_filename = NULL;
03032 char * fparea_filename = NULL;
03033 char ** fncoef_filename = NULL;
03034 char ** fscoef_filename = NULL;
03035 char ** tncoef_filename = NULL;
03036 char ** tscoef_filename = NULL;
03037 char * sfit_filename = NULL;
03038 char * snfit_filename = NULL;
03039
03040 char * label;
03041 int novar;
03042
03043 int ixyz_bot , ixyz_top ;
03044
03045
03046 (void) COX_clock_time() ;
03047
03048
03049 printf ("\n\n");
03050 printf ("Program: %s \n", PROGRAM_NAME);
03051 printf ("Author: %s \n", PROGRAM_AUTHOR);
03052 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
03053 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
03054 printf ("\n");
03055
03056
03057
03058 machdep() ;
03059 {
03060 int new_argc ; char ** new_argv ;
03061 addto_args( argc , argv , &new_argc , &new_argv ) ;
03062 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
03063 }
03064
03065
03066
03067 initialize_program (argc, argv, &ignore,
03068 &nname, &sname, &nmodel, &smodel,
03069 &r, &p, &npname, &spname,
03070 &min_nconstr, &max_nconstr, &min_sconstr, &max_sconstr,
03071 &nabs, &nrand, &nbest, &rms_min, &fdisp,
03072 &input_filename, &tfilename,
03073 &freg_filename, &frsqr_filename,
03074 &fsmax_filename, &ftmax_filename, &fpmax_filename,
03075 &farea_filename, &fparea_filename, &fncoef_filename,
03076 &fscoef_filename, &tncoef_filename, &tscoef_filename,
03077 &sfit_filename, &snfit_filename,
03078 &dset_time, &nxyz, &ts_length, &x_array, &ts_array,
03079 &par_rdcd, &par_full, &tpar_full,
03080 &rmsreg_vol, &freg_vol, &rsqr_vol,
03081 &smax_vol, &tmax_vol, &pmax_vol, &area_vol, &parea_vol,
03082 &ncoef_vol, &scoef_vol, &tncoef_vol, &tscoef_vol,
03083 &sfit_vol, &snfit_vol, &option_data);
03084
03085 #if 0
03086 #ifdef SAVE_RAN
03087 RAN_setup( nmodel , smodel , r , p , nabs ,
03088 min_nconstr, max_nconstr,
03089 min_sconstr, max_sconstr,
03090 ts_length, x_array, nrand ) ;
03091 #endif
03092 #endif
03093
03094 ixyz_bot = 0 ; ixyz_top = nxyz ;
03095
03096 #ifdef PROC_MAX
03097 if( proc_numjob > 1 ){
03098 int vv , nvox=nxyz , nper , pp , nv ;
03099 pid_t newpid ;
03100
03101
03102
03103 if( mask_vol != NULL ){
03104 for( vv=nvox=0 ; vv < nxyz ; vv++ )
03105 if( mask_vol[vv] != 0 ) nvox++ ;
03106 }
03107
03108 if( nvox < proc_numjob ){
03109
03110 proc_numjob = 1 ;
03111
03112 } else {
03113
03114
03115
03116 nper = nvox / proc_numjob ;
03117 if( mask_vol == NULL ){
03118 proc_vox_bot[0] = 0 ;
03119 for( pp=0 ; pp < proc_numjob ; pp++ ){
03120 proc_vox_top[pp] = proc_vox_bot[pp] + nper ;
03121 if( pp < proc_numjob-1 ) proc_vox_bot[pp+1] = proc_vox_top[pp] ;
03122 }
03123 proc_vox_top[proc_numjob-1] = nxyz ;
03124 } else {
03125 proc_vox_bot[0] = 0 ;
03126 for( pp=0 ; pp < proc_numjob ; pp++ ){
03127 for( nv=0,vv=proc_vox_bot[pp] ;
03128 nv < nper && vv < nxyz ; vv++ ){
03129 if( mask_vol[vv] != 0 ) nv++ ;
03130 }
03131 proc_vox_top[pp] = vv ;
03132 if( pp < proc_numjob-1 ) proc_vox_bot[pp+1] = proc_vox_top[pp] ;
03133 }
03134 proc_vox_top[proc_numjob-1] = nxyz ;
03135 }
03136
03137
03138
03139 DSET_load(dset_time) ;
03140
03141
03142
03143 fprintf(stderr,"++ Voxels in dataset: %d\n",nxyz) ;
03144 if( nvox < nxyz )
03145 fprintf(stderr,"++ Voxels in mask: %d\n",nvox) ;
03146 fprintf(stderr,"++ Voxels per job: %d\n",nper) ;
03147
03148 for( pp=1 ; pp < proc_numjob ; pp++ ){
03149 ixyz_bot = proc_vox_bot[pp] ;
03150 ixyz_top = proc_vox_top[pp] ;
03151 proc_ind = pp ;
03152 newpid = fork() ;
03153 if( newpid == -1 ){
03154 fprintf(stderr,"** Can't fork job #%d! Error exit!\n",pp);
03155 exit(1) ;
03156 }
03157 if( newpid == 0 ) break ;
03158 proc_pid[pp] = newpid ;
03159 iochan_sleep(10) ;
03160 }
03161 if( pp == proc_numjob ){
03162 ixyz_bot = proc_vox_bot[0] ;
03163 ixyz_top = proc_vox_top[0] ;
03164 proc_ind = 0 ;
03165 }
03166 fprintf(stderr,"++ Job #%d: processing voxels %d to %d; elapsed time=%.3f\n",
03167 proc_ind,ixyz_bot,ixyz_top-1,COX_clock_time()) ;
03168 }
03169 }
03170 #endif
03171
03172 if( proc_numjob == 1 )
03173 fprintf(stderr,"++ Calculations starting; elapsed time=%.3f\n",COX_clock_time()) ;
03174
03175
03176
03177 for (iv = ixyz_bot; iv < ixyz_top; iv++)
03178 {
03179
03180 if (mask_vol != NULL)
03181 if (mask_vol[iv] == 0) continue;
03182
03183
03184
03185 read_ts_array (dset_time, iv, ts_length, ignore, ts_array);
03186
03187
03188
03189 calc_reduced_model (ts_length, r, x_array, ts_array,
03190 par_rdcd, &sse_rdcd);
03191
03192
03193
03194 calc_full_model (nmodel, smodel, r, p,
03195 min_nconstr, max_nconstr, min_sconstr, max_sconstr,
03196 ts_length, x_array, ts_array, par_rdcd, sse_rdcd, nabs,
03197 nrand, nbest, rms_min, par_full, &sse_full, &novar);
03198
03199
03200
03201 analyze_results (nmodel, smodel, r, p, novar,
03202 min_nconstr, max_nconstr, min_sconstr, max_sconstr,
03203 ts_length, x_array,
03204 par_rdcd, sse_rdcd, par_full, sse_full,
03205 &rmsreg, &freg, &rsqr, &smax, &tmax, &pmax,
03206 &area, &parea, tpar_full);
03207
03208
03209
03210 if (freg >= fdisp && proc_ind == 0 )
03211 {
03212 report_results (nname, sname, r, p, npname, spname, ts_length,
03213 par_rdcd, sse_rdcd, par_full, sse_full, tpar_full,
03214 rmsreg, freg, rsqr, smax, tmax, pmax,
03215 area, parea, &label);
03216 printf ("\n\nVoxel #%d\n", iv);
03217 printf ("%s \n", label);
03218 }
03219
03220
03221
03222 save_results (iv, nmodel, smodel, r, p, novar, ts_length, x_array,
03223 par_full, tpar_full, rmsreg, freg, rsqr, smax,
03224 tmax, pmax, area, parea, rmsreg_vol,
03225 freg_vol, rsqr_vol, smax_vol,
03226 tmax_vol, pmax_vol, area_vol, parea_vol,ncoef_vol,
03227 scoef_vol, tncoef_vol, tscoef_vol, sfit_vol, snfit_vol);
03228 }
03229
03230
03231
03232
03233
03234 #ifdef PROC_MAX
03235 if( proc_numjob > 1 ){
03236 if( proc_ind > 0 ){
03237 fprintf(stderr,"++ Job #%d finished; elapsed time=%.3f\n",proc_ind,COX_clock_time()) ;
03238 _exit(0) ;
03239
03240 } else {
03241 int pp ;
03242 fprintf(stderr,"++ Job #0 waiting for children to finish; elapsed time=%.3f\n",COX_clock_time()) ;
03243 for( pp=1 ; pp < proc_numjob ; pp++ )
03244 waitpid( proc_pid[pp] , NULL , 0 ) ;
03245 fprintf(stderr,"++ Job #0 now finishing up; elapsed time=%.3f\n",COX_clock_time()) ;
03246 }
03247
03248
03249
03250 }
03251 #endif
03252 if( proc_numjob == 1 )
03253 fprintf(stderr,"++ Calculations finished; elapsed time=%.3f\n",COX_clock_time()) ;
03254
03255
03256
03257 THD_delete_3dim_dataset( dset_time , False ) ; dset_time = NULL ;
03258
03259
03260
03261 output_results (r, p, min_nconstr, max_nconstr, min_sconstr, max_sconstr,
03262 nxyz, ts_length, rmsreg_vol, freg_vol, rsqr_vol,
03263 smax_vol, tmax_vol, pmax_vol, area_vol, parea_vol,
03264 ncoef_vol, scoef_vol, tncoef_vol, tscoef_vol,
03265 sfit_vol, snfit_vol,
03266 input_filename, freg_filename, frsqr_filename,
03267 fsmax_filename, ftmax_filename,
03268 fpmax_filename, farea_filename, fparea_filename,
03269 fncoef_filename, fscoef_filename,
03270 tncoef_filename, tscoef_filename,
03271 sfit_filename, snfit_filename, &option_data);
03272
03273
03274
03275 terminate_program (r, p, ts_length, &x_array, &ts_array,
03276 &nname, &npname, &par_rdcd, &min_nconstr, &max_nconstr,
03277 &sname, &spname, &par_full, &tpar_full,
03278 &min_sconstr, &max_sconstr,
03279 &rmsreg_vol, &freg_vol, &rsqr_vol,
03280 &smax_vol, &tmax_vol, &pmax_vol, &area_vol, &parea_vol,
03281 &ncoef_vol, &scoef_vol, &tncoef_vol, &tscoef_vol,
03282 &sfit_vol, &snfit_vol, &input_filename,
03283 &freg_filename, &frsqr_filename,
03284 &fsmax_filename, &ftmax_filename,
03285 &fpmax_filename, &farea_filename, &fparea_filename,
03286 &fncoef_filename, &fscoef_filename,
03287 &tncoef_filename, &tscoef_filename,
03288 &sfit_filename, &snfit_filename);
03289
03290 fprintf(stderr,"++ Program finished; elapsed time=%.3f\n",COX_clock_time()) ;
03291 exit (0);
03292 }
