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 #define PROGRAM_NAME "3dTSgen"
00037 #define PROGRAM_AUTHOR "B. Douglas Ward"
00038 #define PROGRAM_DATE "09 September 1999"
00039
00040
00041 #include <stdio.h>
00042 #include <math.h>
00043 #include <stdlib.h>
00044
00045 #include "mrilib.h"
00046
00047 #include "matrix.h"
00048 #include "simplex.h"
00049 #include "NLfit.h"
00050
00051 #include "matrix.c"
00052 #include "simplex.c"
00053 #include "NLfit.c"
00054
00055
00056 typedef struct NL_options
00057 {
00058 char * bucket_filename;
00059 int numbricks;
00060 int * brick_type;
00061 int * brick_coef;
00062 char ** brick_label;
00063
00064 } NL_options;
00065
00066
00067
00068
00069 static float DELT = 1.0;
00070 static int inTR = 0 ;
00071 static float dsTR = 0.0 ;
00072
00073
00074 static char * commandline = NULL ;
00075
00076
00077
00078
00079
00080
00081
00082 void display_help_menu()
00083 {
00084 printf
00085 (
00086 "This program generates an AFNI 3d+time data set. The time series for \n"
00087 "each voxel is generated according to a user specified signal + noise \n"
00088 "model. \n\n"
00089 "Usage: \n"
00090 "3dTSgen \n"
00091 "-input fname fname = filename of prototype 3d + time data file \n"
00092 "[-inTR] set the TR of the created timeseries to be the TR \n"
00093 " of the prototype dataset \n"
00094 " [The default is to compute with TR = 1.] \n"
00095 " [The model functions are called for a ] \n"
00096 " [time grid of 0, TR, 2*TR, 3*TR, .... ] \n"
00097 "-signal slabel slabel = name of (non-linear) signal model \n"
00098 "-noise nlabel nlabel = name of (linear) noise model \n"
00099 "-sconstr k c d constraints for kth signal parameter: \n"
00100 " c <= gs[k] <= d \n"
00101 "-nconstr k c d constraints for kth noise parameter: \n"
00102 " c+b[k] <= gn[k] <= d+b[k] \n"
00103 "-sigma s s = std. dev. of additive Gaussian noise \n"
00104 "[-voxel num] screen output for voxel #num \n"
00105 "-output fname fname = filename of output 3d + time data file \n"
00106 " \n"
00107 " \n"
00108 "The following commands generate individual AFNI 1 sub-brick datasets: \n"
00109 " \n"
00110 "[-scoef k fname] write kth signal parameter gs[k]; \n"
00111 " output 'fim' is written to prefix filename fname \n"
00112 "[-ncoef k fname] write kth noise parameter gn[k]; \n"
00113 " output 'fim' is written to prefix filename fname \n"
00114 " \n"
00115 " \n"
00116 "The following commands generate one AFNI 'bucket' type dataset: \n"
00117 " \n"
00118 "[-bucket n prefixname] create one AFNI 'bucket' dataset containing \n"
00119 " n sub-bricks; n=0 creates default output; \n"
00120 " output 'bucket' is written to prefixname \n"
00121 "The mth sub-brick will contain: \n"
00122 "[-brick m scoef k label] kth signal parameter regression coefficient\n"
00123 "[-brick m ncoef k label] kth noise parameter regression coefficient \n"
00124 );
00125
00126 exit(0);
00127 }
00128
00129
00130
00131
00132
00133
00134 #define MTEST(ptr) \
00135 if((ptr)==NULL) \
00136 ( fprintf(stderr,"*** Cannot allocate memory for statistics!\n" \
00137 "*** Try using the -workmem option to reduce memory needs,\n" \
00138 "*** or create more swap space in the operating system.\n" \
00139 ), exit(0) )
00140
00141
00142
00143
00144
00145
00146
00147 void initialize_options
00148 (
00149 vfp * nmodel,
00150 vfp * smodel,
00151 char *** npname,
00152 char *** spname,
00153 float ** min_nconstr,
00154 float ** max_nconstr,
00155 float ** min_sconstr,
00156 float ** max_sconstr,
00157 float * sigma,
00158 int * nvoxel,
00159 char ** input_filename,
00160 char ** output_filename,
00161 char *** ncoef_filename,
00162 char *** scoef_filename,
00163 NL_options * option_data
00164 )
00165
00166 {
00167 int ip;
00168
00169
00170
00171 *sigma = 0.0;
00172 *nvoxel = -1;
00173 *smodel = NULL;
00174 *nmodel = NULL;
00175
00176
00177
00178 *npname = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00179 if (*npname == NULL)
00180 NLfit_error ("Unable to allocate memory for noise parameter names");
00181 for (ip = 0; ip < MAX_PARAMETERS; ip++)
00182 {
00183 (*npname)[ip] = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
00184 if ((*npname)[ip] == NULL)
00185 NLfit_error ("Unable to allocate memory for noise parameter names");
00186 }
00187
00188 *spname = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00189 if (*spname == NULL)
00190 NLfit_error ("Unable to allocate memory for signal parameter names");
00191 for (ip = 0; ip < MAX_PARAMETERS; ip++)
00192 {
00193 (*spname)[ip] = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
00194 if ((*spname)[ip] == NULL)
00195 NLfit_error ("Unable to allocate memory for signal parameter names");
00196 }
00197
00198
00199
00200 *min_nconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
00201 if (*min_nconstr == NULL)
00202 NLfit_error ("Unable to allocate memory for min_nconstr");
00203 *max_nconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
00204 if (*max_nconstr == NULL)
00205 NLfit_error ("Unable to allocate memory for max_nconstr");
00206 *min_sconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
00207 if (*min_sconstr == NULL)
00208 NLfit_error ("Unable to allocate memory for min_sconstr");
00209 *max_sconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
00210 if (*max_sconstr == NULL)
00211 NLfit_error ("Unable to allocate memory for max_sconstr");
00212
00213
00214
00215 *input_filename = NULL;
00216 *output_filename = NULL;
00217 *ncoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00218 if (*ncoef_filename == NULL)
00219 NLfit_error ("Unable to allocate memory for ncoef_filename");
00220 *scoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00221 if (*scoef_filename == NULL)
00222 NLfit_error ("Unable to allocate memory for scoef_filename");
00223 for (ip = 0; ip < MAX_PARAMETERS; ip++)
00224 {
00225 (*ncoef_filename)[ip] = NULL;
00226 (*scoef_filename)[ip] = NULL;
00227 }
00228
00229
00230
00231 option_data->bucket_filename = NULL;
00232 option_data->numbricks = -1;
00233 option_data->brick_type = NULL;
00234 option_data->brick_coef = NULL;
00235 option_data->brick_label = NULL;
00236
00237 }
00238
00239
00240
00241
00242
00243
00244
00245 void get_options
00246 (
00247 int argc,
00248 char ** argv,
00249 char ** nname,
00250 char ** sname,
00251 vfp * nmodel,
00252 vfp * smodel,
00253 int * r,
00254 int * p,
00255 char *** npname,
00256 char *** spname,
00257 float ** min_nconstr,
00258 float ** max_nconstr,
00259 float ** min_sconstr,
00260 float ** max_sconstr,
00261 float * sigma,
00262 int * nvoxel,
00263 char ** input_filename,
00264 char ** output_filename,
00265 char *** ncoef_filename,
00266 char *** scoef_filename,
00267
00268 int * nxyz,
00269 int * ts_length,
00270 NL_options * option_data
00271 )
00272
00273 {
00274 const MAX_BRICKS = 100;
00275 int nopt = 1;
00276 int ival, index;
00277 float fval;
00278 char message[MAX_NAME_LENGTH];
00279 int ok;
00280 THD_3dim_dataset * dset_time;
00281
00282 NLFIT_MODEL_array * model_array = NULL;
00283 int im;
00284 int ibrick;
00285 int nbricks;
00286
00287
00288
00289 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
00290
00291
00292
00293 model_array = NLFIT_get_many_MODELs ();
00294 if ((model_array == NULL) || (model_array->num == 0))
00295 NLfit_error ("Unable to locate any models");
00296
00297
00298
00299 initialize_options (nmodel, smodel, npname, spname,
00300 min_nconstr, max_nconstr, min_sconstr, max_sconstr,
00301 sigma, nvoxel, input_filename,
00302 output_filename, ncoef_filename, scoef_filename,
00303 option_data);
00304
00305
00306 while (nopt < argc )
00307 {
00308
00309
00310 if (strncmp(argv[nopt], "-input", 6) == 0)
00311 {
00312 nopt++;
00313 if (nopt >= argc) NLfit_error ("need argument after -input ");
00314 *input_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00315 if (*input_filename == NULL)
00316 NLfit_error ("Unable to allocate memory for input_filename");
00317 strcpy (*input_filename, argv[nopt]);
00318
00319
00320 dset_time = THD_open_one_dataset (*input_filename);
00321 if (dset_time == NULL)
00322 {
00323 sprintf (message,
00324 "Unable to open data file: %s", *input_filename);
00325 NLfit_error (message);
00326 }
00327 *nxyz = dset_time->dblk->diskptr->dimsizes[0]
00328 * dset_time->dblk->diskptr->dimsizes[1]
00329 * dset_time->dblk->diskptr->dimsizes[2] ;
00330 *ts_length = DSET_NUM_TIMES(dset_time);
00331
00332 dsTR = DSET_TIMESTEP(dset_time) ;
00333 if( DSET_TIMEUNITS(dset_time) == UNITS_MSEC_TYPE ) dsTR *= 0.001 ;
00334
00335 THD_delete_3dim_dataset(dset_time, False); dset_time = NULL ;
00336
00337 nopt++;
00338 continue;
00339 }
00340
00341
00342
00343 if( strncmp(argv[nopt],"-inTR",5) == 0 ){
00344 inTR = 1 ;
00345 nopt++ ; continue ;
00346 }
00347
00348
00349 if (strncmp(argv[nopt], "-signal", 7) == 0)
00350 {
00351 nopt++;
00352 if (nopt >= argc) NLfit_error ("need argument after -signal ");
00353 *sname = malloc (sizeof(char) * MAX_NAME_LENGTH);
00354 if (*sname == NULL)
00355 NLfit_error ("Unable to allocate memory for signal model name");
00356 strcpy (*sname, argv[nopt]);
00357 initialize_signal_model (model_array, *sname,
00358 smodel, p, *spname,
00359 *min_sconstr, *max_sconstr);
00360 nopt++;
00361 continue;
00362 }
00363
00364
00365
00366 if (strncmp(argv[nopt], "-noise", 6) == 0)
00367 {
00368 nopt++;
00369 if (nopt >= argc) NLfit_error ("need argument after -noise ");
00370 *nname = malloc (sizeof(char) * MAX_NAME_LENGTH);
00371 if (*nname == NULL)
00372 NLfit_error ("Unable to allocate memory for noise model name");
00373 strcpy (*nname, argv[nopt]);
00374 initialize_noise_model (model_array, *nname,
00375 nmodel, r, *npname,
00376 *min_nconstr, *max_nconstr);
00377 nopt++;
00378 continue;
00379 }
00380
00381
00382
00383 if (strncmp(argv[nopt], "-sconstr", 8) == 0 || strncmp(argv[nopt],"-scnstr",8) == 0 )
00384 {
00385 nopt++;
00386 if (nopt+2 >= argc) NLfit_error("need 3 arguments after -sconstr ");
00387
00388 sscanf (argv[nopt], "%d", &ival);
00389 if ((ival < 0) || (ival >= *p))
00390 NLfit_error ("illegal argument after -sconstr ");
00391 index = ival;
00392 nopt++;
00393
00394 sscanf (argv[nopt], "%f", &fval);
00395 (*min_sconstr)[index] = fval;
00396 nopt++;
00397
00398 sscanf (argv[nopt], "%f", &fval);
00399 (*max_sconstr)[index] = fval;
00400 nopt++;
00401 continue;
00402 }
00403
00404
00405
00406 if (strncmp(argv[nopt], "-nconstr", 8) == 0 || strncmp(argv[nopt],"-ncnstr",8) == 0 )
00407 {
00408 nopt++;
00409 if (nopt+2 >= argc) NLfit_error("need 3 arguments after -nconstr ");
00410
00411 sscanf (argv[nopt], "%d", &ival);
00412 if ((ival < 0) || (ival >= *r))
00413 NLfit_error ("illegal argument after -nconstr ");
00414 index = ival;
00415 nopt++;
00416
00417 sscanf (argv[nopt], "%f", &fval);
00418 (*min_nconstr)[index] = fval;
00419 nopt++;
00420
00421 sscanf (argv[nopt], "%f", &fval);
00422 (*max_nconstr)[index] = fval;
00423 nopt++;
00424 continue;
00425 }
00426
00427
00428
00429 if (strncmp(argv[nopt], "-sigma", 6) == 0)
00430 {
00431 nopt++;
00432 if (nopt >= argc) NLfit_error ("need argument after -sigma ");
00433 sscanf (argv[nopt], "%f", &fval);
00434 if (fval < 0.0)
00435 NLfit_error ("illegal argument after -sigma ");
00436 *sigma = fval;
00437 nopt++;
00438 continue;
00439 }
00440
00441
00442
00443 if (strncmp(argv[nopt], "-voxel", 6) == 0)
00444 {
00445 nopt++;
00446 if (nopt >= argc) NLfit_error ("need argument after -voxel ");
00447 sscanf (argv[nopt], "%d", &ival);
00448 if (ival <= 0)
00449 NLfit_error ("illegal argument after -voxel ");
00450 *nvoxel = ival;
00451 nopt++;
00452 continue;
00453 }
00454
00455
00456
00457 if (strncmp(argv[nopt], "-output", 7) == 0)
00458 {
00459 nopt++;
00460 if (nopt >= argc) NLfit_error ("need argument after -output ");
00461 *output_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00462 if (*output_filename == NULL)
00463 NLfit_error ("Unable to allocate memory for output_filename");
00464 strcpy (*output_filename, argv[nopt]);
00465 nopt++;
00466 continue;
00467 }
00468
00469
00470
00471 if (strncmp(argv[nopt], "-scoef", 6) == 0)
00472 {
00473 nopt++;
00474 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -scoef ");
00475 sscanf (argv[nopt], "%d", &ival);
00476 if ((ival < 0) || (ival >= *p))
00477 NLfit_error ("illegal argument after -scoef ");
00478 index = ival;
00479 nopt++;
00480
00481 (*scoef_filename)[index] = malloc (sizeof(char) * MAX_NAME_LENGTH);
00482 if ((*scoef_filename)[index] == NULL)
00483 NLfit_error ("Unable to allocate memory for scoef_filename");
00484 strcpy ((*scoef_filename)[index], argv[nopt]);
00485
00486 nopt++;
00487 continue;
00488 }
00489
00490
00491
00492 if (strncmp(argv[nopt], "-ncoef", 7) == 0)
00493 {
00494 nopt++;
00495 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -ncoef ");
00496 sscanf (argv[nopt], "%d", &ival);
00497 if ((ival < 0) || (ival >= *r))
00498 NLfit_error ("illegal argument after -ncoef ");
00499 index = ival;
00500 nopt++;
00501
00502 (*ncoef_filename)[index] = malloc (sizeof(char) * MAX_NAME_LENGTH);
00503 if ((*ncoef_filename)[index] == NULL)
00504 NLfit_error ("Unable to allocate memory for ncoef_filename");
00505 strcpy ((*ncoef_filename)[index], argv[nopt]);
00506
00507 nopt++;
00508 continue;
00509 }
00510
00511
00512
00513 if (strncmp(argv[nopt], "-bucket", 7) == 0)
00514 {
00515 nopt++;
00516 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -bucket ");
00517 sscanf (argv[nopt], "%d", &ival);
00518 if ((ival < 0) || (ival > MAX_BRICKS))
00519 NLfit_error ("illegal argument after -bucket ");
00520 nopt++;
00521
00522 option_data->bucket_filename =
00523 malloc (sizeof(char) * MAX_NAME_LENGTH);
00524 if (option_data->bucket_filename == NULL)
00525 NLfit_error ("Unable to allocate memory for bucket_filename");
00526 strcpy (option_data->bucket_filename, argv[nopt]);
00527
00528
00529 if (ival == 0)
00530 nbricks = (*p) + (*r);
00531 else
00532 nbricks = ival;
00533 option_data->numbricks = nbricks;
00534
00535
00536 option_data->brick_type = malloc (sizeof(int) * nbricks);
00537 option_data->brick_coef = malloc (sizeof(int) * nbricks);
00538 option_data->brick_label = malloc (sizeof(char *) * nbricks);
00539 for (ibrick = 0; ibrick < nbricks; ibrick++)
00540 {
00541 option_data->brick_type[ibrick] = -1;
00542 option_data->brick_coef[ibrick] = -1;
00543 option_data->brick_label[ibrick] =
00544 malloc (sizeof(char) * MAX_NAME_LENGTH);
00545 }
00546
00547
00548 if (ival == 0)
00549
00550 {
00551 for (ibrick = 0; ibrick < (*r); ibrick++)
00552 {
00553 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00554 option_data->brick_coef[ibrick] = ibrick;
00555 strcpy (option_data->brick_label[ibrick], (*npname)[ibrick]);
00556 }
00557
00558 for (ibrick = (*r); ibrick < (*p) + (*r); ibrick++)
00559 {
00560 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00561 option_data->brick_coef[ibrick] = ibrick;
00562 strcpy (option_data->brick_label[ibrick],
00563 (*spname)[ibrick-(*r)]);
00564 }
00565 }
00566
00567 nopt++;
00568 continue;
00569 }
00570
00571
00572
00573 if (strncmp(argv[nopt], "-brick", 6) == 0)
00574 {
00575 nopt++;
00576 if (nopt+2 >= argc)
00577 NLfit_error ("need more arguments after -brick ");
00578 sscanf (argv[nopt], "%d", &ibrick);
00579 if ((ibrick < 0) || (ibrick >= option_data->numbricks))
00580 NLfit_error ("illegal argument after -brick ");
00581 nopt++;
00582
00583 if (strncmp(argv[nopt], "scoef", 4) == 0)
00584 {
00585 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00586
00587 nopt++;
00588 sscanf (argv[nopt], "%d", &ival);
00589 if ((ival < 0) || (ival > (*p)))
00590 NLfit_error ("illegal argument after scoef ");
00591 option_data->brick_coef[ibrick] = ival + (*r);
00592
00593 nopt++;
00594 if (nopt >= argc)
00595 NLfit_error ("need more arguments after -brick ");
00596 strcpy (option_data->brick_label[ibrick], argv[nopt]);
00597 }
00598
00599 else if (strncmp(argv[nopt], "ncoef", 4) == 0)
00600 {
00601 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00602
00603 nopt++;
00604 sscanf (argv[nopt], "%d", &ival);
00605 if ((ival < 0) || (ival > (*r)))
00606 NLfit_error ("illegal argument after ncoef ");
00607 option_data->brick_coef[ibrick] = ival;
00608
00609 nopt++;
00610 if (nopt >= argc)
00611 NLfit_error ("need more arguments after -brick ");
00612 strcpy (option_data->brick_label[ibrick], argv[nopt]);
00613 }
00614
00615 else NLfit_error ("unable to interpret options after -brick ");
00616
00617 printf ("ibrick = %d \n", ibrick);
00618 printf ("brick_type = %d \n", option_data->brick_type[ibrick]);
00619 printf ("brick_coef = %d \n", option_data->brick_coef[ibrick]);
00620 printf ("brick_label = %s \n", option_data->brick_label[ibrick]);
00621
00622 nopt++;
00623 continue;
00624 }
00625
00626
00627
00628
00629 { char buf[256] ;
00630 sprintf(buf,"unrecognized command line option: %s",argv[nopt]) ;
00631 NLfit_error (buf);
00632 }
00633 }
00634
00635
00636
00637 DESTROY_MODEL_ARRAY (model_array) ;
00638
00639 }
00640
00641
00642
00643
00644
00645
00646
00647 void check_one_output_file
00648 (
00649 THD_3dim_dataset * dset,
00650 char * filename
00651 )
00652
00653 {
00654 THD_3dim_dataset * new_dset = NULL;
00655 int ierror;
00656
00657
00658
00659 new_dset = EDIT_empty_copy (dset);
00660
00661 ierror = EDIT_dset_items( new_dset,
00662 ADN_prefix, filename ,
00663 ADN_label1, filename ,
00664 ADN_self_name, filename ,
00665 ADN_type, ISHEAD(dset) ? HEAD_FUNC_TYPE :
00666 GEN_FUNC_TYPE ,
00667 ADN_none ) ;
00668
00669 if( ierror > 0 )
00670 {
00671 fprintf(stderr,
00672 "*** %d errors in attempting to create output dataset!\n",
00673 ierror ) ;
00674 exit(1) ;
00675 }
00676
00677 if( THD_is_file(new_dset->dblk->diskptr->header_name) )
00678 {
00679 fprintf(stderr,
00680 "*** Output dataset file %s already exists"
00681 "--cannot continue!\a\n",
00682 new_dset->dblk->diskptr->header_name ) ;
00683 exit(1) ;
00684 }
00685
00686
00687 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00688
00689 }
00690
00691
00692
00693
00694
00695
00696
00697 void check_output_files
00698 (
00699 char * input_filename,
00700 char * output_filename,
00701 char ** ncoef_filename,
00702 char ** scoef_filename,
00703 char * bucket_filename
00704 )
00705
00706 {
00707 THD_3dim_dataset * dset_time;
00708 int ip;
00709
00710
00711 dset_time = THD_open_one_dataset (input_filename);
00712
00713
00714 if (output_filename != NULL)
00715 check_one_output_file (dset_time, output_filename);
00716
00717 for (ip = 0; ip < MAX_PARAMETERS; ip++)
00718 {
00719 if (ncoef_filename[ip] != NULL)
00720 check_one_output_file (dset_time, ncoef_filename[ip]);
00721 if (scoef_filename[ip] != NULL)
00722 check_one_output_file (dset_time, scoef_filename[ip]);
00723 }
00724
00725
00726 if (bucket_filename != NULL)
00727 check_one_output_file (dset_time, bucket_filename);
00728
00729 THD_delete_3dim_dataset (dset_time, False); dset_time = NULL ;
00730 }
00731
00732
00733
00734
00735
00736
00737
00738 void check_for_valid_inputs
00739 (
00740 int r,
00741 int p,
00742 float * min_nconstr,
00743 float * max_nconstr,
00744 float * min_sconstr,
00745 float * max_sconstr,
00746
00747 char * input_filename,
00748 char * output_filename,
00749 char ** ncoef_filename,
00750 char ** scoef_filename,
00751 char * bucket_filename
00752 )
00753
00754 {
00755 int ip;
00756
00757
00758
00759 for (ip = 0; ip < r; ip++)
00760 if (min_nconstr[ip] > max_nconstr[ip])
00761 NLfit_error ("Must have minimum constraints <= maximum constraints");
00762 for (ip = 0; ip < p; ip++)
00763 if (min_sconstr[ip] > max_sconstr[ip])
00764 NLfit_error("Must have mininum constraints <= maximum constraints");
00765
00766
00767
00768 check_output_files (input_filename, output_filename,
00769 ncoef_filename, scoef_filename, bucket_filename);
00770
00771 }
00772
00773
00774
00775
00776
00777
00778
00779 void initialize_program
00780 (
00781 int argc,
00782 char ** argv,
00783 char ** nname,
00784 char ** sname,
00785 vfp * nmodel,
00786 vfp * smodel,
00787 int *r,
00788 int *p,
00789 char *** npname,
00790 char *** spname,
00791 float ** min_nconstr,
00792 float ** max_nconstr,
00793 float ** min_sconstr,
00794 float ** max_sconstr,
00795 float * sigma,
00796 int * nvoxel,
00797
00798 char ** input_filename,
00799 char ** output_filename,
00800 char *** ncoef_filename,
00801 char *** scoef_filename,
00802
00803 int * nxyz,
00804 int * ts_length,
00805
00806 float *** x_array,
00807 float ** ts_array,
00808 short *** d_array,
00809
00810 float ** par_full,
00811
00812 float *** ncoef_vol,
00813 float *** scoef_vol,
00814
00815 NL_options * option_data
00816
00817 )
00818
00819 {
00820 int dimension;
00821 int ip;
00822 int it;
00823
00824
00825
00826 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
00827
00828
00829
00830 get_options (argc, argv, nname, sname, nmodel, smodel,
00831 r, p, npname, spname,
00832 min_nconstr, max_nconstr, min_sconstr, max_sconstr,
00833 sigma, nvoxel,
00834 input_filename, output_filename, ncoef_filename, scoef_filename,
00835 nxyz, ts_length, option_data);
00836
00837
00838 check_for_valid_inputs (*r, *p, *min_nconstr, *max_nconstr,
00839 *min_sconstr, *max_sconstr, *input_filename,
00840 *output_filename, *ncoef_filename, *scoef_filename,
00841 option_data->bucket_filename);
00842
00843
00844 *ts_array = (float *) malloc (sizeof(float) * (*ts_length));
00845 if (*ts_array == NULL)
00846 NLfit_error ("Unable to allocate memory for ts_array");
00847
00848
00849
00850 *x_array = (float **) malloc (sizeof(float *) * (*ts_length));
00851 if (*x_array == NULL)
00852 NLfit_error ("Unable to allocate memory for x_array");
00853 for (it = 0; it < (*ts_length); it++)
00854 {
00855 (*x_array)[it] = (float *) malloc (sizeof(float) * 3);
00856 if ((*x_array)[it] == NULL)
00857 NLfit_error ("Unable to allocate memory for x_array[it]");
00858 }
00859
00860
00861
00862 if( inTR && dsTR > 0.0 ){
00863 DELT = dsTR ;
00864 fprintf(stderr,"--- computing with TR = %g\n",DELT) ;
00865 }
00866
00867 for (it = 0; it < (*ts_length); it++)
00868 {
00869 (*x_array)[it][0] = 1.0;
00870 (*x_array)[it][1] = it * DELT;
00871 (*x_array)[it][2] = (it * DELT) * (it * DELT);
00872 }
00873
00874
00875 dimension = (*r) + (*p);
00876 *par_full = (float *) malloc (sizeof(float) * dimension);
00877 if (*par_full == NULL)
00878 NLfit_error ("Unable to allocate memory for par_full");
00879
00880 *ncoef_vol = (float **) malloc (sizeof(float *) * (*r));
00881 if (*ncoef_vol == NULL)
00882 NLfit_error ("Unable to allocate memory for ncoef_vol");
00883 for (ip = 0; ip < (*r); ip++)
00884 {
00885 (*ncoef_vol)[ip] = (float *) malloc (sizeof(float) * (*nxyz));
00886 if ((*ncoef_vol)[ip] == NULL)
00887 NLfit_error ("Unable to allocate memory for ncoef_vol[ip]");
00888 }
00889
00890 *scoef_vol = (float **) malloc (sizeof(float *) * (*p));
00891 if (*scoef_vol == NULL)
00892 NLfit_error ("Unable to allocate memory for scoef_vol");
00893 for (ip = 0; ip < (*p); ip++)
00894 {
00895 (*scoef_vol)[ip] = (float *) malloc (sizeof(float) * (*nxyz));
00896 if ((*scoef_vol)[ip] == NULL)
00897 NLfit_error ("Unable to allocate memory for scoef_vol[ip]");
00898 }
00899
00900
00901
00902 *d_array = (short **) malloc (sizeof(short *) * (*ts_length));
00903 if (*d_array == NULL)
00904 NLfit_error ("Unable to allocate memory for d_array");
00905 for (it = 0; it < *ts_length; it++)
00906 {
00907 (*d_array)[it] = (short *) malloc (sizeof(short) * (*nxyz));
00908 if ((*d_array)[it] == NULL)
00909 NLfit_error ("Unable to allocate memory for d_array[it]");
00910 }
00911
00912
00913
00914 srand48 (1234567);
00915
00916 }
00917
00918
00919
00920
00921
00922
00923
00924 void generate_ts_array
00925 (
00926 vfp nmodel,
00927 vfp smodel,
00928 int r,
00929 int p,
00930 float * min_nconstr,
00931 float * max_nconstr,
00932 float * min_sconstr,
00933 float * max_sconstr,
00934 int ts_length,
00935 float ** x_array,
00936
00937 float sigma,
00938
00939 float * par_true,
00940 float * ts_array
00941 )
00942
00943 {
00944 int ip;
00945 int it;
00946 float n1, n2;
00947
00948
00949
00950 for (ip = 0; ip < r; ip++)
00951 par_true[ip] = get_random_value (min_nconstr[ip], max_nconstr[ip]);
00952
00953
00954 for (ip = 0; ip < p; ip++)
00955 par_true[ip+r] = get_random_value (min_sconstr[ip], max_sconstr[ip]);
00956
00957
00958 full_model (nmodel, smodel, par_true, par_true + r,
00959 ts_length, x_array, ts_array);
00960
00961
00962 for (it = 0; it < ts_length; it++)
00963 {
00964 normal (&n1, &n2);
00965 ts_array[it] += n1*sigma;
00966 }
00967
00968 }
00969
00970
00971
00972
00973
00974
00975
00976 void save_ts_array
00977 (
00978 int iv,
00979 int ts_length,
00980 float * ts_array,
00981 short ** d_array
00982 )
00983
00984 {
00985 int it;
00986
00987
00988 for (it = 0; it < ts_length; it++)
00989 {
00990 d_array[it][iv] = (short) ts_array[it];
00991 }
00992
00993 }
00994
00995
00996
00997
00998
00999
01000
01001 void write_ts_array
01002 (
01003 int iv,
01004 int ts_length,
01005 float * ts_array,
01006 short ** d_array
01007 )
01008
01009 {
01010 int it;
01011
01012 for (it = 0; it < ts_length; it++)
01013 printf ("ts_array[%d] = %9.3f d_array[%d][%d] = %d \n",
01014 it, ts_array[it], it, iv+1, d_array[it][iv]);
01015
01016 }
01017
01018
01019
01020
01021
01022
01023
01024 void write_parameters
01025 (
01026 int iv,
01027 char * nname,
01028 char * sname,
01029 int r,
01030 int p,
01031 char ** npname,
01032 char ** spname,
01033 float * par_full
01034 )
01035
01036 {
01037 int ip;
01038
01039
01040 if (iv < 0)
01041 printf ("\n\nVoxel Results: \n");
01042 else
01043 printf ("\n\nVoxel #%d\n", iv+1);
01044
01045
01046
01047 printf ("\nFull (%s + %s) Model: \n", nname, sname);
01048
01049 for (ip = 0; ip < r; ip++)
01050 printf ("gn[%d]=%12.6f %s \n", ip, par_full[ip], npname[ip]);
01051
01052 for (ip = 0; ip < p; ip++)
01053 printf ("gs[%d]=%12.6f %s \n", ip, par_full[ip+r], spname[ip]);
01054 }
01055
01056
01057
01058
01059
01060
01061
01062 void save_parameters
01063 (
01064 int iv,
01065 int r,
01066 int p,
01067 float * par_full,
01068
01069 float ** ncoef_vol,
01070 float ** scoef_vol
01071 )
01072
01073 {
01074 int ip;
01075
01076
01077
01078 for (ip = 0; ip < r; ip++)
01079 {
01080 ncoef_vol[ip][iv] = par_full[ip];
01081 }
01082
01083
01084 for (ip = 0; ip < p; ip++)
01085 {
01086 scoef_vol[ip][iv] = par_full[ip+r];
01087 }
01088
01089 }
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104 float EDIT_coerce_autoscale_new( int nxyz ,
01105 int itype,void *ivol , int otype,void *ovol )
01106 {
01107 float fac=0.0 , top ;
01108
01109 if( MRI_IS_INT_TYPE(otype) ){
01110 top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
01111 if (top == 0.0) fac = 0.0;
01112 else fac = MRI_TYPE_maxval[otype]/top;
01113 }
01114
01115 EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
01116 return ( fac );
01117 }
01118
01119
01120
01121
01122
01123
01124
01125
01126 void output_ts_array
01127 (
01128 short ** d_array,
01129 int ts_length,
01130 char * input_filename,
01131 char * filename
01132 )
01133
01134 {
01135 THD_3dim_dataset * dset = NULL;
01136 THD_3dim_dataset * new_dset = NULL;
01137 int ib;
01138 int ierror;
01139 float * fbuf;
01140 float fimfac;
01141 char label[80];
01142
01143
01144
01145 fbuf = (float *) malloc (sizeof(float) * ts_length); MTEST (fbuf);
01146 for (ib = 0; ib < ts_length; ib++) fbuf[ib] = 0.0;
01147
01148
01149
01150 dset = THD_open_one_dataset (input_filename);
01151 new_dset = EDIT_empty_copy (dset);
01152
01153
01154
01155
01156 sprintf (label, "Output prefix: %s", filename);
01157 if( commandline != NULL )
01158 tross_multi_Append_History( new_dset , commandline,label,NULL ) ;
01159 else
01160 tross_Append_History ( new_dset, label);
01161
01162
01163
01164 THD_delete_3dim_dataset (dset, False); dset = NULL ;
01165
01166
01167 ierror = EDIT_dset_items( new_dset ,
01168 ADN_prefix , filename ,
01169 ADN_label1 , filename ,
01170 ADN_self_name , filename ,
01171 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
01172 ADN_nvals , ts_length ,
01173 ADN_none ) ;
01174
01175 if( ierror > 0 ){
01176 fprintf(stderr,
01177 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
01178 exit(1) ;
01179 }
01180
01181 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
01182 fprintf(stderr,
01183 "*** Output dataset file %s already exists--cannot continue!\a\n",
01184 new_dset->dblk->diskptr->header_name ) ;
01185 exit(1) ;
01186 }
01187
01188
01189
01190 for (ib = 0; ib < ts_length; ib++)
01191 mri_fix_data_pointer (d_array[ib], DSET_BRICK(new_dset,ib));
01192
01193
01194
01195
01196 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
01197
01198 THD_load_statistics( new_dset ) ;
01199 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01200 fprintf(stderr,"--- Wrote combined dataset into %s\n",DSET_BRIKNAME(new_dset)) ;
01201
01202
01203
01204 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01205 free (fbuf); fbuf = NULL;
01206
01207 }
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217 void write_afni_data (char * input_filename, int nxyz, char * filename,
01218 float * ffim)
01219 {
01220 int ii;
01221 THD_3dim_dataset * dset=NULL;
01222 THD_3dim_dataset * new_dset=NULL;
01223 int iv;
01224 int ierror;
01225 int ibuf[32];
01226 float fbuf[MAX_STAT_AUX];
01227 float fimfac;
01228 int output_datum;
01229 void * vdif;
01230 int func_type;
01231 float top, func_scale_short;
01232 char label[80];
01233
01234
01235
01236 dset = THD_open_one_dataset (input_filename) ;
01237 if( ! ISVALID_3DIM_DATASET(dset) ){
01238 fprintf(stderr,"*** Unable to open dataset file %s\n",
01239 input_filename);
01240 exit(1) ;
01241 }
01242
01243
01244 new_dset = EDIT_empty_copy( dset ) ;
01245
01246
01247
01248 if( commandline != NULL )
01249 tross_Append_History( new_dset , commandline ) ;
01250 sprintf (label, "Output prefix: %s", filename);
01251 tross_Append_History ( new_dset, label);
01252
01253
01254 iv = DSET_PRINCIPAL_VALUE(dset) ;
01255 output_datum = DSET_BRICK_TYPE(dset,iv) ;
01256 if( output_datum == MRI_byte ) output_datum = MRI_short ;
01257
01258
01259 ibuf[0] = output_datum ;
01260
01261 func_type = FUNC_FIM_TYPE;
01262
01263 ierror = EDIT_dset_items( new_dset ,
01264 ADN_prefix , filename ,
01265 ADN_label1 , filename ,
01266 ADN_self_name , filename ,
01267 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE :
01268 GEN_FUNC_TYPE ,
01269 ADN_func_type , func_type ,
01270 ADN_nvals , FUNC_nvals[func_type] ,
01271 ADN_datum_array , ibuf ,
01272 ADN_ntt , 0 ,
01273 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
01274 ADN_none ) ;
01275
01276 if( ierror > 0 ){
01277 fprintf(stderr,
01278 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
01279 exit(1) ;
01280 }
01281
01282 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
01283 fprintf(stderr,
01284 "*** Output dataset file %s already exists--cannot continue!\a\n",
01285 new_dset->dblk->diskptr->header_name ) ;
01286 exit(1) ;
01287 }
01288
01289
01290
01291 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
01292
01293
01294
01295 vdif = (void *) malloc( mri_datum_size(output_datum) * nxyz );
01296 if (vdif == NULL) NLfit_error ("Unable to allocate memory for vdif");
01297
01298
01299
01300 mri_fix_data_pointer (vdif, DSET_BRICK(new_dset,0));
01301
01302
01303
01304 fimfac = EDIT_coerce_autoscale_new (nxyz,
01305 MRI_float, ffim,
01306 output_datum, vdif);
01307 if (fimfac != 0.0) fimfac = 1.0 / fimfac;
01308
01309
01310
01311 printf("--- Writing combined dataset into %s\n",DSET_BRIKNAME(new_dset)) ;
01312
01313 for( ii=0 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
01314 (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
01315
01316 fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
01317 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
01318
01319 THD_load_statistics( new_dset ) ;
01320 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01321
01322
01323 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01324
01325 }
01326
01327
01328
01329
01330
01331
01332
01333 void write_bucket_data
01334 (
01335 int q,
01336 int p,
01337 int nxyz,
01338 int n,
01339
01340 float ** ncoef_vol,
01341 float ** scoef_vol,
01342
01343 char * input_filename,
01344
01345 NL_options * option_data
01346 )
01347
01348 {
01349 const float EPSILON = 1.0e-10;
01350
01351 THD_3dim_dataset * old_dset = NULL;
01352 THD_3dim_dataset * new_dset = NULL;
01353 char * output_prefix;
01354 char * output_session;
01355 int nbricks, ib;
01356 short ** bar = NULL;
01357 float factor;
01358 int brick_type;
01359 int brick_coef;
01360 char * brick_label;
01361 int ierror;
01362 float * volume;
01363 int dimension;
01364 char label[80];
01365
01366
01367
01368 nbricks = option_data->numbricks;
01369 output_prefix = option_data->bucket_filename;
01370 output_session = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
01371 strcpy (output_session, "./");
01372 dimension = p + q;
01373
01374
01375
01376 bar = (short **) malloc (sizeof(short *) * nbricks);
01377 MTEST (bar);
01378
01379
01380
01381 old_dset = THD_open_one_dataset (input_filename);
01382
01383
01384
01385 new_dset = EDIT_empty_copy (old_dset);
01386
01387
01388
01389 if( commandline != NULL )
01390 tross_Append_History( new_dset , commandline ) ;
01391 sprintf (label, "Output prefix: %s", output_prefix);
01392 tross_Append_History ( new_dset, label);
01393
01394
01395
01396
01397 ierror = EDIT_dset_items (new_dset,
01398 ADN_prefix, output_prefix,
01399 ADN_directory_name, output_session,
01400 ADN_type, HEAD_FUNC_TYPE,
01401 ADN_func_type, FUNC_BUCK_TYPE,
01402 ADN_ntt, 0,
01403 ADN_nvals, nbricks,
01404 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
01405 ADN_none ) ;
01406
01407 if( ierror > 0 )
01408 {
01409 fprintf(stderr,
01410 "*** %d errors in attempting to create output dataset!\n",
01411 ierror);
01412 exit(1);
01413 }
01414
01415 if (THD_is_file(DSET_HEADNAME(new_dset)))
01416 {
01417 fprintf(stderr,
01418 "*** Output dataset file %s already exists--cannot continue!\n",
01419 DSET_HEADNAME(new_dset));
01420 exit(1);
01421 }
01422
01423
01424
01425 THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ;
01426
01427
01428
01429
01430
01431 for (ib = 0; ib < nbricks; ib++)
01432 {
01433
01434 brick_type = option_data->brick_type[ib];
01435 brick_coef = option_data->brick_coef[ib];
01436 brick_label = option_data->brick_label[ib];
01437
01438 if (brick_type == FUNC_FIM_TYPE)
01439 {
01440 if (brick_coef < q)
01441 volume = ncoef_vol[brick_coef];
01442 else if (brick_coef < p+q)
01443 volume = scoef_vol[brick_coef-q];
01444 }
01445
01446
01447 bar[ib] = (short *) malloc (sizeof(short) * nxyz);
01448 MTEST (bar[ib]);
01449 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
01450 MRI_short, bar[ib]);
01451
01452 if (factor < EPSILON) factor = 0.0;
01453 else factor = 1.0 / factor;
01454
01455
01456 EDIT_BRICK_LABEL (new_dset, ib, brick_label);
01457 EDIT_BRICK_FACTOR (new_dset, ib, factor);
01458
01459
01460
01461 EDIT_substitute_brick (new_dset, ib, MRI_short, bar[ib]);
01462
01463 }
01464
01465
01466
01467 THD_load_statistics (new_dset);
01468 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01469 fprintf(stderr,"Wrote bucket dataset ");
01470 fprintf(stderr,"into %s\n", DSET_BRIKNAME(new_dset));
01471
01472
01473
01474 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01475
01476 }
01477
01478
01479
01480
01481
01482
01483
01484 void output_parameters
01485 (
01486 int r,
01487 int p,
01488 int nxyz,
01489 int ts_length,
01490
01491 float ** ncoef_vol,
01492 float ** scoef_vol,
01493
01494 char * input_filename,
01495 char ** ncoef_filename,
01496 char ** scoef_filename,
01497
01498 NL_options * option_data
01499 )
01500
01501 {
01502 int ip;
01503 int dimension;
01504 int numdof, dendof;
01505
01506
01507 dimension = r + p;
01508
01509
01510
01511 if (option_data->numbricks > 0)
01512 write_bucket_data (r, p, nxyz, ts_length, ncoef_vol, scoef_vol,
01513 input_filename, option_data);
01514
01515
01516
01517 for (ip = 0; ip < r; ip++)
01518 {
01519 if (ncoef_filename[ip] != NULL)
01520 {
01521 write_afni_data (input_filename, nxyz, ncoef_filename[ip],
01522 ncoef_vol[ip]);
01523 }
01524 }
01525
01526
01527 for (ip = 0; ip < p; ip++)
01528 {
01529 if (scoef_filename[ip] != NULL)
01530 {
01531 write_afni_data (input_filename, nxyz, scoef_filename[ip],
01532 scoef_vol[ip]);
01533 }
01534 }
01535 }
01536
01537
01538
01539
01540
01541
01542
01543 void terminate_program
01544 (
01545 int r,
01546 int p,
01547 int ts_length,
01548 float *** x_array,
01549 float ** ts_array,
01550 short *** d_array,
01551 char ** nname,
01552 char *** npname,
01553 float ** min_nconstr,
01554 float ** max_nconstr,
01555 char ** sname,
01556 char *** spname,
01557 float ** par_full,
01558 float ** min_sconstr,
01559 float ** max_sconstr,
01560 float *** ncoef_vol,
01561 float *** scoef_vol,
01562 char ** input_filename,
01563 char ** output_filename,
01564 char *** ncoef_filename,
01565 char *** scoef_filename
01566 )
01567
01568 {
01569 int ip;
01570 int it;
01571
01572
01573
01574 if (*nname != NULL)
01575 { free (*nname); *nname = NULL; }
01576 if (*sname != NULL)
01577 { free (*sname); *sname = NULL; }
01578 for (ip = 0; ip < MAX_PARAMETERS; ip++)
01579 {
01580 if ((*npname)[ip] != NULL)
01581 { free ((*npname)[ip]); (*npname)[ip] = NULL; }
01582 if ((*spname)[ip] != NULL)
01583 { free ((*spname)[ip]); (*spname)[ip] = NULL; }
01584 }
01585
01586 if (*npname != NULL)
01587 { free (*npname); *npname = NULL; }
01588 if (*spname != NULL)
01589 { free (*spname); *spname = NULL; }
01590
01591
01592
01593 if (*min_nconstr != NULL) { free (*min_nconstr); *min_nconstr = NULL; }
01594 if (*max_nconstr != NULL) { free (*max_nconstr); *max_nconstr = NULL; }
01595 if (*min_sconstr != NULL) { free (*min_sconstr); *min_sconstr = NULL; }
01596 if (*max_sconstr != NULL) { free (*max_sconstr); *max_sconstr = NULL; }
01597
01598
01599
01600 if (*input_filename != NULL)
01601 { free (*input_filename); *input_filename = NULL; }
01602 if (*output_filename != NULL)
01603 { free (*output_filename); *output_filename = NULL; }
01604
01605 if (*ncoef_filename != NULL)
01606 {
01607 for (ip = 0; ip < MAX_PARAMETERS; ip++)
01608 {
01609 if ((*ncoef_filename)[ip] != NULL)
01610 { free ((*ncoef_filename)[ip]); (*ncoef_filename)[ip] = NULL; }
01611 }
01612 free (*ncoef_filename); *ncoef_filename = NULL;
01613 }
01614
01615 if (*scoef_filename != NULL)
01616 {
01617 for (ip = 0; ip < MAX_PARAMETERS; ip++)
01618 {
01619 if ((*scoef_filename)[ip] != NULL)
01620 { free ((*scoef_filename)[ip]); (*scoef_filename)[ip] = NULL; }
01621 }
01622 free (*scoef_filename); *scoef_filename = NULL;
01623 }
01624
01625
01626
01627 if (*ts_array != NULL) { free (*ts_array); *ts_array = NULL; }
01628
01629
01630
01631 if (*x_array != NULL)
01632 {
01633 for (it = 0; it < ts_length; it++)
01634 if ((*x_array)[it] != NULL)
01635 { free ((*x_array)[it]); (*x_array)[it] = NULL; }
01636 free (*x_array); *x_array = NULL;
01637 }
01638
01639
01640
01641 if (*par_full != NULL) { free (*par_full); *par_full = NULL; }
01642
01643
01644
01645 if (*ncoef_vol != NULL)
01646 {
01647 for (ip = 0; ip < r; ip++)
01648 {
01649 if ((*ncoef_vol)[ip] != NULL)
01650 { free ((*ncoef_vol)[ip]); (*ncoef_vol)[ip] = NULL; }
01651 }
01652 free (*ncoef_vol); *ncoef_vol = NULL;
01653 }
01654
01655 if (*scoef_vol != NULL)
01656 {
01657 for (ip = 0; ip < p; ip++)
01658 {
01659 if ((*scoef_vol)[ip] != NULL)
01660 { free ((*scoef_vol)[ip]); (*scoef_vol)[ip] = NULL; }
01661 }
01662 free (*scoef_vol); *scoef_vol = NULL;
01663 }
01664
01665 }
01666
01667
01668
01669
01670 int main
01671 (
01672 int argc,
01673 char ** argv
01674 )
01675
01676 {
01677 NL_options option_data;
01678
01679
01680 int ts_length;
01681 float ** x_array = NULL;
01682 float * ts_array = NULL;
01683 int nxyz;
01684 int iv;
01685 int nvoxel;
01686 short ** d_array = NULL;
01687
01688
01689 char * nname = NULL;
01690 vfp nmodel;
01691 int r;
01692 char ** npname = NULL;
01693 float * min_nconstr = NULL;
01694 float * max_nconstr = NULL;
01695
01696
01697 char * sname = NULL;
01698 vfp smodel;
01699 int p;
01700 float * par_full = NULL;
01701 char ** spname = NULL;
01702 float * min_sconstr = NULL;
01703 float * max_sconstr = NULL;
01704
01705
01706 float sigma;
01707
01708
01709 float ** ncoef_vol = NULL;
01710 float ** scoef_vol = NULL;
01711
01712
01713 char * input_filename = NULL;
01714 char * output_filename = NULL;
01715 char ** ncoef_filename = NULL;
01716 char ** scoef_filename = NULL;
01717
01718
01719 printf ("\n\n");
01720 printf ("Program: %s \n", PROGRAM_NAME);
01721 printf ("Author: %s \n", PROGRAM_AUTHOR);
01722 printf ("Date: %s \n", PROGRAM_DATE);
01723 printf ("\n");
01724
01725
01726
01727 initialize_program (argc, argv,
01728 &nname, &sname, &nmodel, &smodel,
01729 &r, &p, &npname, &spname,
01730 &min_nconstr, &max_nconstr, &min_sconstr, &max_sconstr,
01731 &sigma, &nvoxel,
01732 &input_filename, &output_filename, &ncoef_filename,
01733 &scoef_filename,
01734 &nxyz, &ts_length, &x_array, &ts_array, &d_array,
01735 &par_full,
01736 &ncoef_vol, &scoef_vol, &option_data);
01737
01738
01739
01740 for (iv = 0; iv < nxyz; iv++)
01741 {
01742
01743
01744 generate_ts_array (nmodel, smodel, r, p,
01745 min_nconstr, max_nconstr, min_sconstr, max_sconstr,
01746 ts_length, x_array, sigma,
01747 par_full, ts_array);
01748
01749
01750
01751 save_ts_array (iv, ts_length, ts_array, d_array);
01752
01753
01754
01755 if (iv == nvoxel-1)
01756 write_ts_array (iv, ts_length, ts_array, d_array);
01757
01758
01759
01760 save_parameters (iv, r, p, par_full, ncoef_vol, scoef_vol);
01761
01762
01763
01764 if (iv == nvoxel-1)
01765 write_parameters (iv, nname, sname, r, p, npname, spname, par_full);
01766
01767 }
01768
01769
01770
01771 output_ts_array (d_array, ts_length, input_filename, output_filename);
01772
01773
01774
01775 output_parameters (r, p, nxyz, ts_length, ncoef_vol, scoef_vol,
01776 input_filename,
01777 ncoef_filename, scoef_filename, &option_data);
01778
01779
01780
01781 terminate_program (r, p, ts_length, &x_array, &ts_array, &d_array,
01782 &nname, &npname, &min_nconstr, &max_nconstr,
01783 &sname, &spname, &par_full, &min_sconstr, &max_sconstr,
01784 &ncoef_vol, &scoef_vol,
01785 &input_filename, &output_filename,
01786 &ncoef_filename, &scoef_filename);
01787 }
01788
01789
01790
01791