Doxygen Source Code Documentation
3dRegAna.c File Reference
#include <stdio.h>#include <math.h>#include "mrilib.h"#include "matrix.h"#include "RegAna.c"Go to the source code of this file.
Data Structures | |
| struct | model |
| struct | RA_options |
Defines | |
| #define | PROGRAM_NAME "3dRegAna" |
| #define | PROGRAM_AUTHOR "B. Douglas Ward" |
| #define | PROGRAM_INITIAL "10 Oct 1997" |
| #define | PROGRAM_LATEST "02 Dec 2002" |
| #define | SUFFIX ".3dregana" |
| #define | DF "df -k ." |
| #define | MAX_XVARS 101 |
| #define | MAX_OBSERVATIONS 1000 |
| #define | MAX_NAME_LENGTH THD_MAX_NAME |
| #define | MEGA 1048576 |
| #define | DOPEN(ds, name) |
| #define | SUB_POINTER(ds, vv, ind, ptr) |
Typedefs | |
| typedef model | model |
| typedef RA_options | RA_options |
Functions | |
| void | RA_error (char *message) |
| void | display_help_menu () |
| void | get_dimensions (RA_options *option_data) |
| void | read_afni_data (RA_options *option_data, char *filename, int piece_len, int fim_offset, float *ffim) |
| void | check_piece (char *filename) |
| void | read_piece (char *filename, float *fin, int size) |
| void | write_piece (char *filename, float *fout, int size) |
| void | delete_piece (char *filename) |
| void | allocate_pieces (matrix xdata, model *regmodel, RA_options *option_data, float ***yfimar, float **freg_piece, float **rsqr_piece, float ***coef_piece, float ***tcoef_piece) |
| void | save_pieces (int piece, int piece_len, float *freg_piece, float *rsqr_piece, float **coef_piece, float **tcoef_piece) |
| void | deallocate_pieces (int n, float ***yfimar, float **freg_piece, float **rsqr_piece, float ***coef_piece, float ***tcoef_piece) |
| void | check_volume (char *filename, int num_pieces) |
| void | read_volume (char *filename, float *volume, int nxyz, int piece_size, int num_pieces) |
| void | delete_volume (char *filename, int nxyz, int piece_size, int num_pieces) |
| void | initialize_options (model *regmodel, RA_options *option_data) |
| void | sort_model_indices (model *regmodel) |
| void | get_inputs (int argc, char **argv, matrix *xdata, model *regmodel, RA_options *option_data) |
| void | count_volumes_and_files (model *regmodel, RA_options *option_data) |
| void | break_into_pieces (int num_datasets, RA_options *option_data) |
| void | identify_repeats (matrix *xdata, model *regmodel, RA_options *option_data) |
| void | check_for_valid_inputs (matrix *xdata, model *regmodel, RA_options *option_data) |
| void | check_one_output_file (RA_options *option_data, char *filename) |
| void | check_output_files (RA_options *option_data) |
| void | check_temporary_files (matrix *xdata, model *regmodel, RA_options *option_data) |
| void | check_disk_space (RA_options *option_data) |
| void | initialize_program (int argc, char **argv, matrix *xdata, model *regmodel, RA_options *option_data) |
| void | init_regression_analysis (int p, int q, int *flist, int *rlist, matrix xdata, matrix *x_full, matrix *xtxinv_full, matrix *xtxinvxt_full, matrix *x_base, matrix *xtxinvxt_base, matrix *x_rdcd, matrix *xtxinvxt_rdcd) |
| void | regression_analysis (int N, int p, int q, matrix x_full, matrix xtxinv_full, matrix xtxinvxt_full, matrix x_base, matrix xtxinvxt_base, matrix x_rdcd, matrix xtxinvxt_rdcd, vector y, float rms_min, int *levels, int *counts, int c, float flofmax, float *flof, vector *coef_full, vector *scoef_full, vector *tcoef_full, float *freg, float *rsqr) |
| void | save_voxel (int iv, vector y, float fdisp, model *regmodel, float flof, vector coef, vector tcoef, float freg, float rsqr, float *freg_piece, float *rsqr_piece, float **coef_piece, float **tcoef_piece) |
| void | calculate_results (matrix xdata, model *regmodel, RA_options *option_data) |
| float | EDIT_coerce_autoscale_new (int nxyz, int itype, void *ivol, int otype, void *ovol) |
| void | write_afni_data (RA_options *option_data, char *filename, float *ffim, float *ftr, int func_type, int numdof, int dendof) |
| void | write_bucket_data (matrix xdata, model *regmodel, RA_options *option_data) |
| void | output_results (matrix xdata, model *regmodel, RA_options *option_data) |
| void | terminate_program (matrix *xdata, model *regmodel, RA_options *option_data) |
| int | main (int argc, char **argv) |
Variables | |
| char * | commandline = NULL |
Define Documentation
|
|
Definition at line 104 of file 3dRegAna.c. Referenced by check_disk_space(). |
|
|
Value: do{ int pv ; (ds) = THD_open_dataset((name)) ; \ if( !ISVALID_3DIM_DATASET((ds)) ){ \ fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; } \ if( (ds)->daxes->nxx!=nx || (ds)->daxes->nyy!=ny || \ (ds)->daxes->nzz!=nz ){ \ fprintf(stderr,"*** Axes mismatch: %s\n",(name)) ; exit(1) ; } \ if( DSET_NVALS((ds)) != 1 ){ \ fprintf(stderr,"*** Must specify 1 sub-brick: %s\n",(name));exit(1);}\ THD_load_datablock( (ds)->dblk ) ; \ pv = DSET_PRINCIPAL_VALUE((ds)) ; \ if( DSET_ARRAY((ds),pv) == NULL ){ \ fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); } \ if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){ \ fprintf(stderr,"*** Can't use complex data: %s\n",(name)); exit(1); \ } \ break ; } while (0) Definition at line 269 of file 3dRegAna.c. |
|
|
Definition at line 117 of file 3dRegAna.c. Referenced by check_disk_space(), check_piece(), check_temporary_files(), check_volume(), delete_piece(), delete_volume(), get_inputs(), output_results(), read_piece(), read_volume(), save_pieces(), terminate_program(), write_bucket_data(), and write_piece(). |
|
|
Definition at line 116 of file 3dRegAna.c. Referenced by get_inputs(), initialize_options(), and terminate_program(). |
|
|
Definition at line 115 of file 3dRegAna.c. Referenced by allocate_pieces(), check_output_files(), deallocate_pieces(), get_inputs(), init_regression_analysis(), initialize_options(), save_pieces(), and terminate_program(). |
|
|
Definition at line 118 of file 3dRegAna.c. Referenced by break_into_pieces(). |
|
|
Definition at line 60 of file 3dRegAna.c. Referenced by main(). |
|
|
Definition at line 61 of file 3dRegAna.c. Referenced by main(). |
|
|
Definition at line 62 of file 3dRegAna.c. Referenced by main(). |
|
|
Definition at line 59 of file 3dRegAna.c. Referenced by get_inputs(), initialize_program(), main(), and RA_error(). |
|
|
Value: do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){ \ default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1); \ case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ; \ (ptr) = (void *)( fim + (ind) ) ; \ } break ; \ case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ; \ (ptr) = (void *)( fim + (ind) ) ; \ } break ; \ case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ; \ (ptr) = (void *)( fim + (ind) ) ; \ } break ; } break ; } while(0) Definition at line 292 of file 3dRegAna.c. |
|
|
Definition at line 66 of file 3dRegAna.c. Referenced by check_piece(), delete_piece(), read_piece(), and write_piece(). |
Typedef Documentation
|
|
|
|
|
|
Function Documentation
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 537 of file 3dRegAna.c. References FUNC_FIM_TYPE, FUNC_THR_TYPE, i, malloc, MAX_XVARS, MTEST, and p. Referenced by calculate_results().
00549 {
00550 int n; /* number of datasets */
00551 int p; /* number of parameters in the full model */
00552 int * flist = NULL; /* list of parameters in the full model */
00553 int i; /* dataset index */
00554 int ip; /* parameter index */
00555 int ix; /* x-variable index */
00556 int piece_size; /* number of voxels in datasset piece */
00557 int ib; /* sub-brick index */
00558 int nbricks; /* number of sub-bricks in bucket dataset */
00559
00560
00561 /*----- initialize local variables -----*/
00562 n = xdata.rows;
00563 p = regmodel->p;
00564 flist = regmodel->flist;
00565 piece_size = option_data->piece_size;
00566 nbricks = option_data->numbricks;
00567
00568
00569 /*----- allocate memory space for input data -----*/
00570 *yfimar = (float **) malloc (sizeof(float *) * n);
00571 MTEST (*yfimar);
00572 for (i = 0; i < n; i++)
00573 {
00574 (*yfimar)[i] = (float *) malloc(sizeof(float) * piece_size);
00575 MTEST ((*yfimar)[i]);
00576 }
00577
00578
00579 /*----- allocate memory space for F-statistics data -----*/
00580 for (ip = 0; ip < p; ip++)
00581 {
00582 ix = flist[ip];
00583
00584 if (option_data->fcoef_filename[ix] != NULL)
00585 {
00586 if (*freg_piece == NULL)
00587 {
00588 *freg_piece = (float *) malloc (sizeof(float) * piece_size);
00589 MTEST (*freg_piece);
00590 }
00591 }
00592 }
00593
00594 for (ib = 0; ib < nbricks; ib++)
00595 {
00596 if (option_data->brick_type[ib] == FUNC_FT_TYPE)
00597 {
00598 if (*freg_piece == NULL)
00599 {
00600 *freg_piece = (float *) malloc (sizeof(float) * piece_size);
00601 MTEST (*freg_piece);
00602 }
00603 }
00604 }
00605
00606
00607 /*----- allocate memory space for coef. of mult. determination R^2 -----*/
00608 for (ip = 0; ip < p; ip++)
00609 {
00610 ix = flist[ip];
00611
00612 if (option_data->rcoef_filename[ix] != NULL)
00613 {
00614 if (*rsqr_piece == NULL)
00615 {
00616 *rsqr_piece = (float *) malloc (sizeof(float) * piece_size);
00617 MTEST (*rsqr_piece);
00618 }
00619 }
00620 }
00621
00622 for (ib = 0; ib < nbricks; ib++)
00623 {
00624 if (option_data->brick_type[ib] == FUNC_THR_TYPE)
00625 {
00626 if (*rsqr_piece == NULL)
00627 {
00628 *rsqr_piece = (float *) malloc (sizeof(float) * piece_size);
00629 MTEST (*rsqr_piece);
00630 }
00631 }
00632 }
00633
00634
00635 /*----- allocate memory space for parameters and t-stats -----*/
00636 *coef_piece = (float **) malloc (sizeof(float *) * MAX_XVARS);
00637 MTEST (*coef_piece);
00638
00639 *tcoef_piece = (float **) malloc (sizeof(float *) * MAX_XVARS);
00640 MTEST (*tcoef_piece);
00641
00642 for (ix = 0; ix < MAX_XVARS; ix++)
00643 {
00644 (*coef_piece)[ix] = NULL;
00645 (*tcoef_piece)[ix] = NULL;
00646 }
00647
00648 for (ip = 0; ip < p; ip++)
00649 {
00650 ix = flist[ip];
00651
00652 if ( (option_data->fcoef_filename[ix] != NULL)
00653 || (option_data->rcoef_filename[ix] != NULL)
00654 || (option_data->tcoef_filename[ix] != NULL))
00655 {
00656 (*coef_piece)[ix] =
00657 (float *) malloc (sizeof(float) * piece_size);
00658 MTEST ((*coef_piece)[ix]);
00659 }
00660
00661 if (option_data->tcoef_filename[ix] != NULL)
00662 {
00663 (*tcoef_piece)[ix] =
00664 (float *) malloc (sizeof(float) * piece_size);
00665 MTEST ((*tcoef_piece)[ix]);
00666 }
00667 }
00668
00669 for (ib = 0; ib < nbricks; ib++)
00670 {
00671 if (option_data->brick_type[ib] == FUNC_FIM_TYPE)
00672 {
00673 ix = option_data->brick_coef[ib];
00674 if ((*coef_piece)[ix] == NULL)
00675 {
00676 (*coef_piece)[ix] = (float *) malloc (sizeof(float)*piece_size);
00677 MTEST ((*coef_piece)[ix]);
00678 }
00679 }
00680 }
00681
00682 for (ib = 0; ib < nbricks; ib++)
00683 {
00684 if (option_data->brick_type[ib] == FUNC_TT_TYPE)
00685 {
00686 ix = option_data->brick_coef[ib];
00687 if ((*tcoef_piece)[ix] == NULL)
00688 {
00689 (*tcoef_piece)[ix] = (float *) malloc (sizeof(float)*piece_size);
00690 MTEST ((*tcoef_piece)[ix]);
00691 }
00692 }
00693 }
00694
00695 }
|
|
||||||||||||
|
Definition at line 1636 of file 3dRegAna.c. References MEGA. Referenced by initialize_program().
01641 {
01642 int num_vols; /* number of output volumes */
01643
01644
01645 /*----- count number of distinct floating point volumes required -----*/
01646 num_vols = option_data->numf + option_data->numr
01647 + option_data->numt + option_data->numc;
01648
01649 /*----- calculate number of voxels per piece -----*/
01650 option_data->piece_size = option_data->workmem * MEGA
01651 / ((num_datasets + num_vols) * sizeof(float));
01652 if (option_data->piece_size > option_data->nxyz)
01653 option_data->piece_size = option_data->nxyz;
01654
01655 /*----- calculate number of pieces per dataset -----*/
01656 option_data->num_pieces = (option_data->nxyz + option_data->piece_size - 1)
01657 / option_data->piece_size;
01658
01659 printf ("num_pieces = %d piece_size = %d \n",
01660 option_data->num_pieces, option_data->piece_size);
01661
01662 }
|
|
||||||||||||||||
|
Definition at line 2321 of file 3dRegAna.c. References allocate_pieces(), deallocate_pieces(), vector::elts, i, init_regression_analysis(), matrix_destroy(), matrix_initialize(), matrix_print(), p, q, read_afni_data(), regression_analysis(), save_pieces(), save_voxel(), vector_create(), vector_destroy(), and vector_initialize().
02327 {
02328 int * flist = NULL; /* list of parameters in the full model */
02329 int * rlist = NULL; /* list of parameters in the rdcd model */
02330 int n; /* number of data points */
02331 int p; /* number of parameters in the full model */
02332 int q; /* number of parameters in the rdcd model */
02333 float flof; /* F-statistic for lack of fit */
02334 float freg; /* F-statistic for the full regression model */
02335 float rsqr; /* coeff. of multiple determination R^2 */
02336 vector coef; /* regression parameters */
02337 vector scoef; /* std. devs. for regression parameters */
02338 vector tcoef; /* t-statistics for regression parameters */
02339
02340 matrix x_full; /* extracted X matrix for full model */
02341 matrix xtxinv_full; /* matrix: 1/(X'X) for full model */
02342 matrix xtxinvxt_full; /* matrix: (1/(X'X))X' for full model */
02343 matrix x_base; /* extracted X matrix for baseline model */
02344 matrix xtxinvxt_base; /* matrix: (1/(X'X))X' for baseline model */
02345 matrix x_rdcd; /* extracted X matrix for reduced model */
02346 matrix xtxinvxt_rdcd; /* matrix: (1/(X'X))X' for reduced model */
02347 vector y; /* vector of measured data */
02348
02349 int i; /* dataset index */
02350 int nxyz; /* number of voxels per dataset */
02351 int num_datasets; /* total number of datasets */
02352 int piece_size; /* number of voxels in dataset piece */
02353 int num_pieces; /* dataset is divided into this many pieces */
02354 int piece; /* piece index */
02355 int piece_len; /* number of voxels in current piece */
02356 int fim_offset; /* array offset to current piece */
02357 int ivox; /* index of voxel in current piece */
02358 int nvox; /* index of voxel within entire volume */
02359
02360 float ** yfimar = NULL; /* array of pieces of Y-datasets */
02361 float * freg_piece = NULL; /* piece for F-statistics */
02362 float * rsqr_piece = NULL; /* piece for R^2 */
02363 float ** coef_piece = NULL; /* pieces for regression coefficients */
02364 float ** tcoef_piece = NULL; /* pieces for t-statistics */
02365
02366
02367 /*----- Initialize matrices and vectors -----*/
02368 matrix_initialize (&x_full);
02369 matrix_initialize (&xtxinv_full);
02370 matrix_initialize (&xtxinvxt_full);
02371 matrix_initialize (&x_base);
02372 matrix_initialize (&xtxinvxt_base);
02373 matrix_initialize (&x_rdcd);
02374 matrix_initialize (&xtxinvxt_rdcd);
02375 vector_initialize (&coef);
02376 vector_initialize (&scoef);
02377 vector_initialize (&tcoef);
02378 vector_initialize (&y);
02379
02380
02381 /*----- initialize local variables -----*/
02382 n = xdata.rows;
02383 p = regmodel->p;
02384 flist = regmodel->flist;
02385 q = regmodel->q;
02386 rlist = regmodel->rlist;
02387 nxyz = option_data->nxyz;
02388 piece_size = option_data->piece_size;
02389 num_pieces = option_data->num_pieces;
02390
02391
02392 /*----- allocate memory space for pieces -----*/
02393 allocate_pieces (xdata, regmodel, option_data, &yfimar,
02394 &freg_piece, &rsqr_piece, &coef_piece, &tcoef_piece);
02395
02396
02397 /*----- initialization for the regression analysis -----*/
02398 init_regression_analysis (p, q, flist, rlist, xdata,
02399 &x_full, &xtxinv_full, &xtxinvxt_full,
02400 &x_base, &xtxinvxt_base, &x_rdcd, &xtxinvxt_rdcd);
02401
02402
02403 vector_create (n, &y);
02404
02405
02406 if (option_data->fdisp >= 0)
02407 {
02408 printf ("\n");
02409 printf ("X matrix: \n");
02410 matrix_print (xdata);
02411 }
02412
02413
02414 /*----- loop over the pieces of the input datasets -----*/
02415 nvox = -1;
02416 for (piece = 0; piece < num_pieces; piece++)
02417 {
02418 printf ("piece = %d \n", piece);
02419 fim_offset = piece * piece_size;
02420 if (piece < num_pieces-1)
02421 piece_len = piece_size;
02422 else
02423 piece_len = nxyz - fim_offset;
02424
02425 /*----- read in the Y-data -----*/
02426 for (i = 0; i < n; i++)
02427 read_afni_data (option_data, option_data->yname[i],
02428 piece_len, fim_offset, yfimar[i]);
02429
02430
02431 /*----- loop over voxels in this piece -----*/
02432 for (ivox = 0; ivox < piece_len; ivox++)
02433 {
02434 nvox += 1;
02435
02436
02437 /*----- extract Y-data for this voxel -----*/
02438 for (i = 0; i < n; i++)
02439 y.elts[i] = yfimar[i][ivox];
02440
02441
02442 /*----- calculate results for this voxel -----*/
02443 regression_analysis (n, p, q,
02444 x_full, xtxinv_full, xtxinvxt_full, x_base,
02445 xtxinvxt_base, x_rdcd, xtxinvxt_rdcd,
02446 y, option_data->rms_min, option_data->levels,
02447 option_data->counts, option_data->c,
02448 option_data->flofmax, &flof,
02449 &coef, &scoef, &tcoef, &freg, &rsqr);
02450
02451
02452 /*----- save results for this voxel -----*/
02453 save_voxel (ivox, y, option_data->fdisp,
02454 regmodel, flof, coef, tcoef, freg, rsqr,
02455 freg_piece, rsqr_piece, coef_piece, tcoef_piece);
02456
02457
02458 } /* loop over voxels within this piece */
02459
02460
02461 /*----- save piece data into external files -----*/
02462 save_pieces (piece, piece_len,
02463 freg_piece, rsqr_piece, coef_piece, tcoef_piece);
02464
02465
02466 } /* loop over pieces */
02467
02468
02469 /*----- deallocate memory for pieces -----*/
02470 deallocate_pieces (n, &yfimar, &freg_piece, &rsqr_piece,
02471 &coef_piece, &tcoef_piece);
02472
02473
02474 /*----- Dispose of matrices -----*/
02475 vector_destroy (&y);
02476 vector_destroy (&tcoef);
02477 vector_destroy (&scoef);
02478 vector_destroy (&coef);
02479 matrix_destroy (&xtxinvxt_rdcd);
02480 matrix_destroy (&x_rdcd);
02481 matrix_destroy (&xtxinvxt_base);
02482 matrix_destroy (&x_base);
02483 matrix_destroy (&xtxinvxt_full);
02484 matrix_destroy (&xtxinv_full);
02485 matrix_destroy (&x_full);
02486
02487 }
|
|
|
Definition at line 1973 of file 3dRegAna.c. References DF, and MAX_NAME_LENGTH. Referenced by initialize(), and initialize_program().
01977 {
01978 char ch; /* user response */
01979 int nxyz; /* number of voxels per image */
01980 int nmax; /* maximum number of disk files */
01981 char filename[MAX_NAME_LENGTH]; /* output file name */
01982
01983
01984 /*----- initialize local variables -----*/
01985 nxyz = option_data->nxyz;
01986
01987 /*----- first, determine space required for temporary volume data -----*/
01988 nmax = 4 * nxyz * (option_data->numf + option_data->numr + option_data->numt
01989 + option_data->numc);
01990
01991 /*----- determine additional space required for output files -----*/
01992 nmax += 4 * nxyz * option_data->numfiles;
01993
01994 printf("\nThis problem requires approximately %d MB of free disk space.\n",
01995 (nmax/1000000) + 1);
01996
01997
01998 /*----- Determine how much disk space is available. -----*/
01999 printf ("Summary of available disk space: \n\n");
02000 system (DF);
02001 printf ("\nDo you wish to proceed? ");
02002 ch = getchar();
02003 if ((ch != 'Y') && (ch != 'y')) exit(0);
02004 printf ("\n");
02005 }
|
|
||||||||||||||||
|
Definition at line 1755 of file 3dRegAna.c. References RA_error.
01761 {
01762 int ib; /* sub-brick index */
01763
01764
01765 /*----- check data set dimensions -----*/
01766 if (xdata->cols < 2)
01767 RA_error ("Too few X variables ");
01768 if (xdata->rows < 3)
01769 RA_error ("Too few data sets for Y-observations ");
01770
01771 /*----- check model dimensions -----*/
01772 if (regmodel->q < 1)
01773 RA_error ("Reduced regression model is too small");
01774 if (regmodel->p <= regmodel->q)
01775 RA_error ("Full regression model is too small");
01776 if (xdata->rows <= regmodel->p)
01777 RA_error ("Too few data sets for fitting full regression model ");
01778 if (regmodel->rlist[0] != 0)
01779 RA_error ("Reduced model must include variable index 0");
01780
01781
01782 /*----- check for repeat observations -----*/
01783 if (option_data->flofmax >= 0.0)
01784 {
01785 if (xdata->rows <= option_data->c)
01786 RA_error ("Cannot conduct F-test for lack of fit (no repeat obs.) ");
01787 if (option_data->c <= regmodel->p)
01788 RA_error ("Cannot conduct F-test for lack of fit (c <= p) ");
01789 }
01790
01791 /*----- check bucket dataset parameters -----*/
01792 if (option_data->numbricks > 0)
01793 for (ib = 0; ib < option_data->numbricks; ib++)
01794 {
01795 if (option_data->brick_type[ib] < 0)
01796 RA_error ("Must specify contents of each brick in the bucket");
01797 }
01798 }
|
|
||||||||||||
|
Definition at line 1807 of file 3dRegAna.c. References ADN_directory_name, ADN_label1, ADN_none, ADN_prefix, ADN_self_name, ADN_type, THD_3dim_dataset::dblk, THD_datablock::diskptr, EDIT_dset_items(), EDIT_empty_copy(), GEN_FUNC_TYPE, HEAD_FUNC_TYPE, THD_diskptr::header_name, ISHEAD, ISVALID_3DIM_DATASET, THD_delete_3dim_dataset(), THD_is_file(), and THD_open_dataset().
01812 {
01813 THD_3dim_dataset * dset=NULL; /* input afni data set pointer */
01814 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
01815 int ierror; /* number of errors in editing data */
01816
01817
01818 /*----- read first dataset -----*/
01819 dset = THD_open_dataset (option_data->first_dataset ) ;
01820 if( ! ISVALID_3DIM_DATASET(dset) ){
01821 fprintf(stderr,"*** Unable to open dataset file %s\n",
01822 option_data->first_dataset);
01823 exit(1) ;
01824 }
01825
01826 /*-- make an empty copy of this dataset, for eventual output --*/
01827 new_dset = EDIT_empty_copy( dset ) ;
01828
01829
01830 ierror = EDIT_dset_items( new_dset ,
01831 ADN_prefix , filename ,
01832 ADN_label1 , filename ,
01833 ADN_directory_name , option_data->session ,
01834 ADN_self_name , filename ,
01835 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE :
01836 GEN_FUNC_TYPE ,
01837 ADN_none ) ;
01838
01839 if( ierror > 0 ){
01840 fprintf(stderr,
01841 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
01842 exit(1) ;
01843 }
01844
01845 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
01846 fprintf(stderr,
01847 "*** Output dataset file %s already exists--cannot continue!\a\n",
01848 new_dset->dblk->diskptr->header_name ) ;
01849 exit(1) ;
01850 }
01851
01852 /*----- deallocate memory -----*/
01853 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
01854 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01855
01856 }
|
|
|
Definition at line 1865 of file 3dRegAna.c. References check_one_output_file(), and MAX_XVARS.
01869 {
01870 int ix; /* x-variable index */
01871
01872
01873 for (ix = 0; ix < MAX_XVARS; ix++)
01874 {
01875 if (option_data->fcoef_filename[ix] != NULL)
01876 check_one_output_file (option_data, option_data->fcoef_filename[ix]);
01877 if (option_data->rcoef_filename[ix] != NULL)
01878 check_one_output_file (option_data, option_data->rcoef_filename[ix]);
01879 if (option_data->tcoef_filename[ix] != NULL)
01880 check_one_output_file (option_data, option_data->tcoef_filename[ix]);
01881 }
01882
01883 if (option_data->bucket_filename != NULL)
01884 check_one_output_file (option_data, option_data->bucket_filename);
01885
01886 }
|
|
|
Definition at line 387 of file 3dRegAna.c. References far, MAX_NAME_LENGTH, RA_error, and SUFFIX. Referenced by check_volume().
00391 {
00392 FILE * far = NULL; /* pointer to temporary file */
00393 char sfilename[MAX_NAME_LENGTH]; /* name of temporary file */
00394 char message[MAX_NAME_LENGTH]; /* error message */
00395
00396
00397 /*----- see if piece file already exists -----*/
00398 strcpy (sfilename, filename);
00399 strcat (sfilename, SUFFIX);
00400 far = fopen (sfilename, "r");
00401 if (far != NULL)
00402 {
00403 sprintf (message, "temporary file %s already exists. ", sfilename);
00404 RA_error (message);
00405 }
00406
00407 }
|
|
||||||||||||||||
|
Definition at line 1895 of file 3dRegAna.c. References check_volume(), MAX_NAME_LENGTH, and p. Referenced by initialize(), and initialize_program().
01901 {
01902 int p; /* number of parameters in full model */
01903 int ip, ix; /* parameter index */
01904 int num_pieces; /* dataset is divided into this many pieces */
01905 char filename[MAX_NAME_LENGTH]; /* name for temporary data file */
01906
01907
01908 /*----- initialize local variables -----*/
01909 p = regmodel->p;
01910 num_pieces = option_data->num_pieces;
01911
01912
01913 /*----- check for F-statistics data file -----*/
01914 if (option_data->numf > 0)
01915 {
01916 strcpy (filename, "freg");
01917 check_volume (filename, num_pieces);
01918 }
01919
01920
01921 /*----- check for R^2 data file -----*/
01922 if (option_data->numr > 0)
01923 {
01924 strcpy (filename, "rsqr");
01925 check_volume (filename, num_pieces);
01926 }
01927
01928
01929 /*----- check for t-statistics data files -----*/
01930 if (option_data->numt > 0)
01931 {
01932 for (ip = 0; ip < p; ip++)
01933 {
01934 ix = regmodel->flist[ip];
01935
01936 if (option_data->tcoef_filename[ix] != NULL)
01937 {
01938 sprintf (filename, "tcoef.%d", ix);
01939 check_volume (filename, num_pieces);
01940 }
01941 }
01942 }
01943
01944
01945 /*----- check for regression coefficients data files -----*/
01946 if (option_data->numc > 0)
01947 {
01948 for (ip = 0; ip < p; ip++)
01949 {
01950 ix = regmodel->flist[ip];
01951
01952 if ( (option_data->fcoef_filename[ix] != NULL)
01953 || (option_data->rcoef_filename[ix] != NULL)
01954 || (option_data->tcoef_filename[ix] != NULL) )
01955 {
01956 sprintf (filename, "coef.%d", ix);
01957 check_volume (filename, num_pieces);
01958 }
01959 }
01960 }
01961
01962
01963 }
|
|
||||||||||||
|
Definition at line 841 of file 3dRegAna.c. References check_piece(), and MAX_NAME_LENGTH. Referenced by check_temporary_files().
00846 {
00847 int piece; /* piece index */
00848 char sfilename[MAX_NAME_LENGTH]; /* name for temporary data file */
00849
00850
00851 /*----- loop over the temporary data file pieces -----*/
00852 for (piece = 0; piece < num_pieces; piece++)
00853 {
00854 /*----- check for piece data file -----*/
00855 sprintf (sfilename, "%s.p%d", filename, piece);
00856 check_piece (sfilename);
00857
00858 } /* loop over pieces */
00859
00860 }
|
|
||||||||||||
|
Definition at line 1539 of file 3dRegAna.c. References FUNC_FIM_TYPE, FUNC_THR_TYPE, and p. Referenced by initialize_program().
01544 {
01545 int ip; /* parameter index */
01546 int ix; /* x-variable index */
01547 int ib; /* sub-brick index */
01548 int p; /* number of parameters in the model */
01549 int nbricks; /* number of bricks in bucket dataset */
01550
01551
01552 /*----- initialize local variables -----*/
01553 nbricks = option_data->numbricks;
01554 p = regmodel->p;
01555
01556
01557 /*----- count number of volumes and files for F-statistics -----*/
01558 for (ip = 0; ip < p; ip++)
01559 {
01560 ix = regmodel->flist[ip];
01561
01562 if (option_data->fcoef_filename[ix] != NULL)
01563 {
01564 option_data->numf = 1;
01565 option_data->numfiles += 1;
01566 }
01567 }
01568
01569 for (ib = 0; ib < nbricks; ib++)
01570 if (option_data->brick_type[ib] == FUNC_FT_TYPE)
01571 option_data->numf = 1;
01572
01573
01574 /*----- count number of volumes and files for R^2 -----*/
01575 for (ip = 0; ip < p; ip++)
01576 {
01577 ix = regmodel->flist[ip];
01578
01579 if (option_data->rcoef_filename[ix] != NULL)
01580 {
01581 option_data->numr = 1;
01582 option_data->numfiles += 1;
01583 }
01584 }
01585
01586 for (ib = 0; ib < nbricks; ib++)
01587 if (option_data->brick_type[ib] == FUNC_THR_TYPE)
01588 option_data->numr = 1;
01589
01590
01591 /*----- count number of volumes and files for t-statistics -----*/
01592 for (ip = 0; ip < p; ip++)
01593 {
01594 ix = regmodel->flist[ip];
01595
01596 if (option_data->tcoef_filename[ix] != NULL)
01597 {
01598 option_data->numt += 1;
01599 option_data->numfiles += 1;
01600 }
01601
01602 else
01603 for (ib = 0; ib < nbricks; ib++)
01604 if ((option_data->brick_type[ib] == FUNC_TT_TYPE)
01605 && (option_data->brick_coef[ib] == ix))
01606 option_data->numt += 1;
01607 }
01608
01609
01610 /*----- count number of volumes for regression coefficients -----*/
01611 for (ip = 0; ip < p; ip++)
01612 {
01613 ix = regmodel->flist[ip];
01614
01615 if ( (option_data->fcoef_filename[ix] != NULL)
01616 || (option_data->rcoef_filename[ix] != NULL)
01617 || (option_data->tcoef_filename[ix] != NULL))
01618 option_data->numc += 1;
01619
01620 else
01621 for (ib = 0; ib < nbricks; ib++)
01622 if ((option_data->brick_type[ib] == FUNC_FIM_TYPE)
01623 && (option_data->brick_coef[ib] == ix))
01624 option_data->numc += 1;
01625 }
01626
01627 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 767 of file 3dRegAna.c. References free, i, and MAX_XVARS. Referenced by calculate_results().
00776 {
00777 int i; /* dataset index */
00778 int ip; /* parameter index */
00779
00780
00781 /*----- deallocate memory space for input data -----*/
00782 if (*yfimar != NULL)
00783 {
00784 for (i = 0; i < n; i++)
00785 {
00786 free ((*yfimar)[i]); (*yfimar)[i] = NULL;
00787 }
00788 free (*yfimar);
00789 *yfimar = NULL;
00790 }
00791
00792 /*----- deallocate memory space for F-statistics data -----*/
00793 if (*freg_piece != NULL)
00794 {
00795 free (*freg_piece);
00796 *freg_piece = NULL;
00797 }
00798
00799 /*----- deallocate memory space for coef. of mult. determination R^2 -----*/
00800 if (*rsqr_piece != NULL)
00801 {
00802 free (*rsqr_piece);
00803 *rsqr_piece = NULL;
00804 }
00805
00806 /*----- deallocate memory space for regression coefficients -----*/
00807 if (*coef_piece != NULL)
00808 {
00809 for (ip = 0; ip < MAX_XVARS; ip++)
00810 if ((*coef_piece)[ip] != NULL)
00811 {
00812 free ((*coef_piece)[ip]);
00813 (*coef_piece)[ip] = NULL;
00814 }
00815 free (*coef_piece);
00816 *coef_piece = NULL;
00817 }
00818
00819 /*----- deallocate memory space for t-statistics -----*/
00820 if (*tcoef_piece != NULL)
00821 {
00822 for (ip = 0; ip < MAX_XVARS; ip++)
00823 if ((*tcoef_piece)[ip] != NULL)
00824 {
00825 free ((*tcoef_piece)[ip]);
00826 (*tcoef_piece)[ip] = NULL;
00827 }
00828 free (*tcoef_piece);
00829 *tcoef_piece = NULL;
00830 }
00831
00832 }
|
|
|
Definition at line 514 of file 3dRegAna.c. References MAX_NAME_LENGTH, and SUFFIX. Referenced by delete_volume().
00518 {
00519 char sfilename[MAX_NAME_LENGTH]; /* file name */
00520
00521 /*----- build file name -----*/
00522 strcpy (sfilename, filename);
00523 strcat (sfilename, SUFFIX);
00524
00525 /*----- delete 3d data file -----*/
00526 remove (sfilename);
00527
00528 }
|
|
||||||||||||||||||||
|
Definition at line 912 of file 3dRegAna.c. References delete_piece(), and MAX_NAME_LENGTH. Referenced by terminate_program(), and write_bucket_data().
00919 {
00920 int piece; /* piece index */
00921 char sfilename[MAX_NAME_LENGTH]; /* name for temporary data file */
00922
00923
00924 /*----- loop over the temporary data file pieces -----*/
00925 for (piece = 0; piece < num_pieces; piece++)
00926 {
00927
00928 /*----- delete the piece data file -----*/
00929 sprintf (sfilename, "%s.p%d", filename, piece);
00930 delete_piece (sfilename);
00931
00932 } /* loop over pieces */
00933
00934 }
|
|
|
Definition at line 194 of file 3dRegAna.c. References MASTER_SHORTHELP_STRING.
00195 {
00196 printf
00197 (
00198 "This program performs multiple linear regression analysis. \n\n"
00199 "Usage: \n"
00200 "3dRegAna \n"
00201 "-rows n number of input datasets \n"
00202 "-cols m number of X variables \n"
00203 "-xydata X11 X12 ... X1m filename X variables and Y observations \n"
00204 " . . \n"
00205 " . . \n"
00206 " . . \n"
00207 "-xydata Xn1 Xn2 ... Xnm filename X variables and Y observations \n"
00208 " \n"
00209 "-model i1 ... iq : j1 ... jr definition of linear regression model; \n"
00210 " reduced model: \n"
00211 " Y = f(Xj1,...,Xjr) \n"
00212 " full model: \n"
00213 " Y = f(Xj1,...,Xjr,Xi1,...,Xiq) \n"
00214 " \n"
00215 "[-diskspace] print out disk space required for program execution\n"
00216 "[-workmem mega] number of megabytes of RAM to use for statistical \n"
00217 " workspace (default = 12) \n"
00218 "[-rmsmin r] r = minimum rms error to reject constant model \n"
00219 "[-fdisp fval] display (to screen) results for those voxels \n"
00220 " whose F-statistic is > fval \n"
00221 " \n"
00222 "[-flof alpha] alpha = minimum p value for F due to lack of fit \n"
00223 " \n"
00224 " \n"
00225 "The following commands generate individual AFNI 2 sub-brick datasets: \n"
00226 " \n"
00227 "[-fcoef k prefixname] estimate of kth regression coefficient \n"
00228 " along with F-test for the regression \n"
00229 " is written to AFNI `fift' dataset \n"
00230 "[-rcoef k prefixname] estimate of kth regression coefficient \n"
00231 " along with coef. of mult. deter. R^2 \n"
00232 " is written to AFNI `fith' dataset \n"
00233 "[-tcoef k prefixname] estimate of kth regression coefficient \n"
00234 " along with t-test for the coefficient \n"
00235 " is written to AFNI `fitt' dataset \n"
00236 " \n"
00237 " \n"
00238 "The following commands generate one AFNI 'bucket' type dataset: \n"
00239 " \n"
00240 "[-bucket n prefixname] create one AFNI 'bucket' dataset having \n"
00241 " n sub-bricks; n=0 creates default output;\n"
00242 " output 'bucket' is written to prefixname \n"
00243 "The mth sub-brick will contain: \n"
00244 "[-brick m coef k label] kth parameter regression coefficient \n"
00245 "[-brick m fstat label] F-stat for significance of regression \n"
00246 "[-brick m rstat label] coefficient of multiple determination R^2 \n"
00247 "[-brick m tstat k label] t-stat for kth regression coefficient \n"
00248 "\n" );
00249
00250 printf
00251 (
00252 "\n"
00253 "N.B.: For this program, the user must specify 1 and only 1 sub-brick \n"
00254 " with each -xydata command. That is, if an input dataset contains\n"
00255 " more than 1 sub-brick, a sub-brick selector must be used, e.g.: \n"
00256 " -xydata 2.17 4.59 7.18 'fred+orig[3]' \n"
00257 );
00258
00259 printf("\n" MASTER_SHORTHELP_STRING ) ;
00260
00261 exit(0);
00262 }
|
|
||||||||||||||||||||||||
|
compute start time of this timeseries * Definition at line 2497 of file 3dRegAna.c. References EDIT_coerce_scale_type(), MCW_vol_amax(), MRI_IS_INT_TYPE, and top.
02505 {
02506 float fac=0.0 , top ;
02507
02508 if( MRI_IS_INT_TYPE(otype) ){
02509 top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
02510 if (top == 0.0) fac = 0.0;
02511 else fac = MRI_TYPE_maxval[otype]/top;
02512 }
02513
02514 EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
02515 return ( fac );
02516 }
|
|
|
Definition at line 312 of file 3dRegAna.c. References THD_3dim_dataset::daxes, ISVALID_3DIM_DATASET, THD_dataxes::nxx, THD_dataxes::nyy, THD_dataxes::nzz, THD_delete_3dim_dataset(), and THD_open_dataset().
00316 {
00317
00318 THD_3dim_dataset * dset=NULL;
00319
00320 /*----- read first dataset to get dimensions, etc. -----*/
00321
00322 dset = THD_open_dataset( option_data->first_dataset ) ;
00323 if( ! ISVALID_3DIM_DATASET(dset) ){
00324 fprintf(stderr,"*** Unable to open dataset file %s\n",
00325 option_data->first_dataset);
00326 exit(1) ;
00327 }
00328
00329 /*----- data set dimensions in voxels -----*/
00330 option_data->nx = dset->daxes->nxx ;
00331 option_data->ny = dset->daxes->nyy ;
00332 option_data->nz = dset->daxes->nzz ;
00333 option_data->nxyz = option_data->nx * option_data->ny * option_data->nz ;
00334
00335 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00336
00337 }
|
|
||||||||||||||||||||||||
|
Definition at line 1077 of file 3dRegAna.c. References AFNI_logger(), argc, cols, display_help_menu(), FUNC_FIM_TYPE, FUNC_THR_TYPE, initialize_options(), ISVALID_3DIM_DATASET, malloc, matrix_create(), MAX_NAME_LENGTH, MAX_OBSERVATIONS, MAX_XVARS, p, PROGRAM_NAME, RA_error, rows, sort_model_indices(), THD_delete_3dim_dataset(), and THD_open_dataset(). Referenced by initialize_program().
01085 {
01086 const MAX_BRICKS = 100; /* max. number of bricks in the bucket */
01087 int nopt = 1; /* input option argument counter */
01088 int ival, index; /* integer input */
01089 float fval; /* float input */
01090 int rows, cols; /* number rows and columns for x matrix */
01091 int irows, jcols; /* data point counters */
01092 THD_3dim_dataset * dset=NULL; /* test whether data set exists */
01093 char message[MAX_NAME_LENGTH]; /* error message */
01094 char label[MAX_NAME_LENGTH]; /* sub-brick label */
01095 int ibrick; /* sub-brick index */
01096 int nbricks; /* number of sub-bricks in the bucket */
01097 int p; /* number of parameters */
01098 int ip, ix; /* parameter indices */
01099
01100
01101 /*----- does user request help menu? -----*/
01102 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
01103
01104
01105 /*----- add to program log -----*/
01106 AFNI_logger (PROGRAM_NAME,argc,argv);
01107
01108
01109 /*----- initialize the input options -----*/
01110 initialize_options (regmodel, option_data);
01111
01112
01113 /*----- main loop over input options -----*/
01114 while (nopt < argc && argv[nopt][0] == '-')
01115 {
01116
01117 /*----- -datum type -----*/
01118 if( strncmp(argv[nopt],"-datum",6) == 0 ){
01119 if( ++nopt >= argc ) RA_error("need an argument after -datum!") ;
01120
01121 if( strcmp(argv[nopt],"short") == 0 ){
01122 option_data->datum = MRI_short ;
01123 } else if( strcmp(argv[nopt],"float") == 0 ){
01124 option_data->datum = MRI_float ;
01125 } else {
01126 char buf[256] ;
01127 sprintf(buf,"-datum of type '%s' is not supported in 3dRegAna.",
01128 argv[nopt] ) ;
01129 RA_error(buf) ;
01130 }
01131 nopt++ ; continue ; /* go to next arg */
01132 }
01133
01134
01135 /*----- -session dirname -----*/
01136 if( strncmp(argv[nopt],"-session",8) == 0 ){
01137 nopt++ ;
01138 if( nopt >= argc ) RA_error("need argument after -session!") ;
01139 strcpy(option_data->session , argv[nopt++]) ;
01140 continue ;
01141 }
01142
01143
01144 /*----- -diskspace -----*/
01145 if( strncmp(argv[nopt],"-diskspace",10) == 0 )
01146 {
01147 option_data->diskspace = 1;
01148 nopt++ ; continue ; /* go to next arg */
01149 }
01150
01151
01152 /*----- -workmem megabytes -----*/
01153
01154 if( strncmp(argv[nopt],"-workmem",6) == 0 ){
01155 nopt++ ;
01156 if( nopt >= argc ) RA_error ("need argument after -workmem!") ;
01157 sscanf (argv[nopt], "%d", &ival);
01158 if( ival <= 0 ) RA_error ("illegal argument after -workmem!") ;
01159 option_data->workmem = ival ;
01160 nopt++ ; continue ;
01161 }
01162
01163
01164 /*----- -rmsmin r -----*/
01165 if (strncmp(argv[nopt], "-rmsmin", 7) == 0)
01166 {
01167 nopt++;
01168 if (nopt >= argc) RA_error ("need argument after -rmsmin ");
01169 sscanf (argv[nopt], "%f", &fval);
01170 if (fval < 0.0)
01171 RA_error ("illegal argument after -rmsmin ");
01172 option_data->rms_min = fval;
01173 nopt++;
01174 continue;
01175 }
01176
01177
01178 /*----- -flof alpha -----*/
01179 if (strncmp(argv[nopt], "-flof", 6) == 0)
01180 {
01181 nopt++;
01182 if (nopt >= argc) RA_error ("need argument after -flof ");
01183 sscanf (argv[nopt], "%f", &fval);
01184 if ((fval < 0.0) || (fval > 1.0))
01185 RA_error ("illegal argument after -flof ");
01186 option_data->flofmax = fval;
01187 nopt++;
01188 continue;
01189 }
01190
01191
01192 /*----- -fdisp fval -----*/
01193 if (strncmp(argv[nopt], "-fdisp", 6) == 0)
01194 {
01195 nopt++;
01196 if (nopt >= argc) RA_error ("need argument after -fdisp ");
01197 sscanf (argv[nopt], "%f", &fval);
01198 option_data->fdisp = fval;
01199 nopt++;
01200 continue;
01201 }
01202
01203
01204 /*----- -rows n -----*/
01205 if (strncmp(argv[nopt], "-rows", 5) == 0)
01206 {
01207 nopt++;
01208 if (nopt >= argc) RA_error ("need argument after -rows ");
01209 sscanf (argv[nopt], "%d", &ival);
01210 if ((ival <= 0) || (ival > MAX_OBSERVATIONS))
01211 RA_error ("illegal argument after -rows ");
01212 rows = ival;
01213 nopt++;
01214 continue;
01215 }
01216
01217
01218 /*----- -cols m -----*/
01219 if (strncmp(argv[nopt], "-cols", 5) == 0)
01220 {
01221 nopt++;
01222 if (nopt >= argc) RA_error ("need argument after -cols ");
01223 sscanf (argv[nopt], "%d", &ival);
01224 if ((ival <= 0) || (ival > MAX_XVARS))
01225 RA_error ("illegal argument after -cols ");
01226 cols = ival;
01227 nopt++;
01228
01229 matrix_create (rows, cols+1, xdata);
01230 irows = -1;
01231 continue;
01232 }
01233
01234
01235 /*----- -xydata x1 ... xm y -----*/
01236 if (strncmp(argv[nopt], "-xydata", 7) == 0)
01237 {
01238 nopt++;
01239 if (nopt + cols >= argc)
01240 RA_error ("need cols+1 arguments after -xydata ");
01241
01242 irows++;
01243 if (irows > rows-1) RA_error ("too many data files");
01244
01245 xdata->elts[irows][0] = 1.0;
01246 for (jcols = 1; jcols <= cols; jcols++)
01247 {
01248 sscanf (argv[nopt], "%f", &fval);
01249 xdata->elts[irows][jcols] = fval;
01250 nopt++;
01251 }
01252
01253 /*--- check whether input files exist ---*/
01254 dset = THD_open_dataset( argv[nopt] ) ;
01255 if( ! ISVALID_3DIM_DATASET(dset) )
01256 {
01257 sprintf(message,"Unable to open dataset file %s\n", argv[nopt]);
01258 RA_error (message);
01259 }
01260 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
01261
01262 option_data->yname[irows]
01263 = malloc (sizeof(char) * MAX_NAME_LENGTH);
01264 strcpy (option_data->yname[irows], argv[nopt]);
01265 nopt++;
01266 continue;
01267 }
01268
01269
01270 /*----- -model -----*/
01271 if (strncmp(argv[nopt], "-model", 6) == 0)
01272 {
01273 nopt++;
01274 if (nopt >= argc) RA_error ("need arguments after -model ");
01275
01276 while ((nopt < argc)
01277 && (strncmp(argv[nopt], "-", 1) != 0)
01278 && (strncmp(argv[nopt], ":", 1) != 0))
01279 {
01280 sscanf (argv[nopt], "%d", &ival);
01281 if ((ival <= 0) || (ival > cols))
01282 RA_error ("illegal argument after -model ");
01283 regmodel->flist[regmodel->p] = ival;
01284
01285 regmodel->p++;
01286 nopt++;
01287 }
01288
01289 if (strncmp(argv[nopt], ":", 1) == 0)
01290 {
01291 nopt++;
01292
01293 while ((nopt < argc)
01294 && (strncmp(argv[nopt], "-", 1) != 0))
01295 {
01296 sscanf (argv[nopt], "%d", &ival);
01297 if ((ival < 0) || (ival > cols))
01298 RA_error ("illegal argument after -model ");
01299 regmodel->flist[regmodel->p] = ival;
01300 regmodel->rlist[regmodel->q] = ival;
01301 regmodel->p++;
01302 regmodel->q++;
01303 nopt++;
01304 }
01305 }
01306
01307 /*----- sort model variable indices into ascending order -----*/
01308 sort_model_indices (regmodel);
01309
01310 continue;
01311 }
01312
01313
01314 /*----- -fcoef k filename -----*/
01315 if (strncmp(argv[nopt], "-fcoef", 6) == 0)
01316 {
01317 nopt++;
01318 if (nopt+1 >= argc) RA_error ("need 2 arguments after -fcoef ");
01319 sscanf (argv[nopt], "%d", &ival);
01320 if ((ival < 0) || (ival > cols))
01321 RA_error ("illegal argument after -fcoef ");
01322 index = ival;
01323 nopt++;
01324
01325 option_data->fcoef_filename[index] =
01326 malloc (sizeof(char) * MAX_NAME_LENGTH);
01327 if (option_data->fcoef_filename[index] == NULL)
01328 RA_error ("Unable to allocate memory for fcoef_filename");
01329 strcpy (option_data->fcoef_filename[index], argv[nopt]);
01330
01331 nopt++;
01332 continue;
01333 }
01334
01335
01336 /*----- -rcoef k filename -----*/
01337 if (strncmp(argv[nopt], "-rcoef", 6) == 0)
01338 {
01339 nopt++;
01340 if (nopt+1 >= argc) RA_error ("need 2 arguments after -rcoef ");
01341 sscanf (argv[nopt], "%d", &ival);
01342 if ((ival < 0) || (ival > cols))
01343 RA_error ("illegal argument after -rcoef ");
01344 index = ival;
01345 nopt++;
01346
01347 option_data->rcoef_filename[index] =
01348 malloc (sizeof(char) * MAX_NAME_LENGTH);
01349 if (option_data->rcoef_filename[index] == NULL)
01350 RA_error ("Unable to allocate memory for rcoef_filename");
01351 strcpy (option_data->rcoef_filename[index], argv[nopt]);
01352
01353 nopt++;
01354 continue;
01355 }
01356
01357
01358 /*----- -tcoef k filename -----*/
01359 if (strncmp(argv[nopt], "-tcoef", 6) == 0)
01360 {
01361 nopt++;
01362 if (nopt+1 >= argc) RA_error ("need 2 arguments after -tcoef ");
01363 sscanf (argv[nopt], "%d", &ival);
01364 if ((ival < 0) || (ival > cols))
01365 RA_error ("illegal argument after -tcoef ");
01366 index = ival;
01367 nopt++;
01368
01369 option_data->tcoef_filename[index] =
01370 malloc (sizeof(char) * MAX_NAME_LENGTH);
01371 if (option_data->tcoef_filename[index] == NULL)
01372 RA_error ("Unable to allocate memory for tcoef_filename");
01373 strcpy (option_data->tcoef_filename[index], argv[nopt]);
01374
01375 nopt++;
01376 continue;
01377 }
01378
01379
01380 /*----- -bucket n prefixname -----*/
01381 if (strncmp(argv[nopt], "-bucket", 7) == 0)
01382 {
01383 nopt++;
01384 if (nopt+1 >= argc) RA_error ("need 2 arguments after -bucket ");
01385 sscanf (argv[nopt], "%d", &ival);
01386 if ((ival < 0) || (ival > MAX_BRICKS))
01387 RA_error ("illegal argument after -bucket ");
01388 option_data->numbricks = ival;
01389 nopt++;
01390
01391 option_data->bucket_filename =
01392 malloc (sizeof(char) * MAX_NAME_LENGTH);
01393 if (option_data->bucket_filename == NULL)
01394 RA_error ("Unable to allocate memory for bucket_filename");
01395 strcpy (option_data->bucket_filename, argv[nopt]);
01396
01397 /*----- set number of sub-bricks in the bucket -----*/
01398 p = regmodel->p;
01399 if (ival == 0)
01400 nbricks = 2*p + 2;
01401 else
01402 nbricks = ival;
01403 option_data->numbricks = nbricks;
01404
01405 /*----- allocate memory and initialize bucket dataset options -----*/
01406 option_data->brick_type = malloc (sizeof(int) * nbricks);
01407 option_data->brick_coef = malloc (sizeof(int) * nbricks);
01408 option_data->brick_label = malloc (sizeof(char *) * nbricks);
01409 for (ibrick = 0; ibrick < nbricks; ibrick++)
01410 {
01411 option_data->brick_type[ibrick] = -1;
01412 option_data->brick_coef[ibrick] = -1;
01413 option_data->brick_label[ibrick] =
01414 malloc (sizeof(char) * MAX_NAME_LENGTH);
01415 }
01416
01417 if (ival == 0) /*----- throw everything into the bucket -----*/
01418 {
01419 for (ip = 0; ip < p; ip++)
01420 {
01421 ix = regmodel->flist[ip];
01422
01423 ibrick = 2*ip;
01424 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01425 option_data->brick_coef[ibrick] = ix;
01426 sprintf (label, "Coef #%d est.", ix);
01427 strcpy (option_data->brick_label[ibrick], label);
01428
01429 ibrick = 2*ip + 1;
01430 option_data->brick_type[ibrick] = FUNC_TT_TYPE;
01431 option_data->brick_coef[ibrick] = ix;
01432 sprintf (label, "Coef #%d t-stat", ix);
01433 strcpy (option_data->brick_label[ibrick], label);
01434 }
01435
01436 ibrick = 2*p;
01437 option_data->brick_type[ibrick] = FUNC_FT_TYPE;
01438 strcpy (option_data->brick_label[ibrick], "F-stat Regression");
01439
01440 ibrick = 2*p + 1;
01441 option_data->brick_type[ibrick] = FUNC_THR_TYPE;
01442 strcpy (option_data->brick_label[ibrick], "R^2 Regression");
01443
01444 }
01445
01446 nopt++;
01447 continue;
01448 }
01449
01450
01451 /*----- -brick m type k label -----*/
01452 if (strncmp(argv[nopt], "-brick", 6) == 0)
01453 {
01454 nopt++;
01455 if (nopt+2 >= argc) RA_error ("need more arguments after -brick ");
01456 sscanf (argv[nopt], "%d", &ibrick);
01457 if ((ibrick < 0) || (ibrick >= option_data->numbricks))
01458 RA_error ("illegal argument after -brick ");
01459 nopt++;
01460
01461 if (strncmp(argv[nopt], "coef", 4) == 0)
01462 {
01463 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01464
01465 nopt++;
01466 sscanf (argv[nopt], "%d", &ival);
01467 if ((ival < 0) || (ival > cols))
01468 RA_error ("illegal argument after coef ");
01469 option_data->brick_coef[ibrick] = ival;
01470
01471 nopt++;
01472 if (nopt >= argc) RA_error ("need more arguments after -brick ");
01473 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01474
01475 }
01476 else if (strncmp(argv[nopt], "tstat", 4) == 0)
01477 {
01478 option_data->brick_type[ibrick] = FUNC_TT_TYPE;
01479
01480 nopt++;
01481 sscanf (argv[nopt], "%d", &ival);
01482 if ((ival < 0) || (ival > cols))
01483 RA_error ("illegal argument after tstat ");
01484 option_data->brick_coef[ibrick] = ival;
01485
01486 nopt++;
01487 if (nopt >= argc) RA_error ("need more arguments after -brick ");
01488 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01489
01490 }
01491 else if (strncmp(argv[nopt], "fstat", 4) == 0)
01492 {
01493 option_data->brick_type[ibrick] = FUNC_FT_TYPE;
01494
01495 nopt++;
01496 if (nopt >= argc) RA_error ("need more arguments after -brick ");
01497 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01498
01499 }
01500 else if (strncmp(argv[nopt], "rstat", 4) == 0)
01501 {
01502 option_data->brick_type[ibrick] = FUNC_THR_TYPE;
01503
01504 nopt++;
01505 if (nopt >= argc) RA_error ("need more arguments after -brick ");
01506 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01507
01508 }
01509 else
01510 {
01511 sprintf(message,"Unrecognized option after -brick: %s\n",
01512 argv[nopt]);
01513 RA_error (message);
01514 }
01515
01516 nopt++;
01517 continue;
01518 }
01519
01520
01521 /*----- unknown command -----*/
01522 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
01523 RA_error (message);
01524
01525 }
01526
01527 /*----- check for fewer than expected datasets -----*/
01528 if (irows < rows-1)
01529 RA_error ("Fewer than expected datasets were entered");
01530 }
|
|
||||||||||||||||
|
Definition at line 1673 of file 3dRegAna.c. References cdff(), i, malloc, MTEST, p, q, and RA_error. Referenced by initialize_program().
01679 {
01680 int i, k; /* matrix row indices */
01681 int j; /* matrix column index */
01682 int match; /* boolean for repeat observation found */
01683
01684 int which; /* indicator for F-distribution calculation */
01685 double p, q; /* cumulative probabilities under F-dist. */
01686 double f; /* F value corresponding to probability p */
01687 double dfn, dfd; /* numerator and denominator dof */
01688 int status; /* calculation status */
01689 double bound; /* search bound */
01690
01691
01692 /*----- allocate memory -----*/
01693 option_data->levels = (int *) malloc (sizeof(int) * (xdata->rows));
01694 MTEST (option_data->levels);
01695 option_data->counts = (int *) malloc (sizeof(int) * (xdata->rows));
01696 MTEST (option_data->counts);
01697
01698
01699 /*----- initialization -----*/
01700 for (i = 0; i < xdata->rows; i++)
01701 option_data->counts[i] = 0;
01702 option_data->levels[0] = 0;
01703 option_data->counts[0] = 1;
01704 option_data->c = 1;
01705
01706
01707 /*----- loop over observations -----*/
01708 for (i = 1; i < xdata->rows; i++)
01709 {
01710 option_data->levels[i] = option_data->c;
01711
01712 /*----- determine if this is a repeat observation -----*/
01713 for (k = 0; k < i; k++)
01714 {
01715 match = 1;
01716 for (j = 1; j < xdata->cols; j++)
01717 if (xdata->elts[i][j] != xdata->elts[k][j]) match = 0;
01718
01719 if (match)
01720 {
01721 option_data->levels[i] = option_data->levels[k];
01722 break;
01723 }
01724 }
01725
01726 /*----- increment count of repeat observations -----*/
01727 k = option_data->levels[i];
01728 option_data->counts[k] ++;
01729
01730 /*----- increment count of distinct observation levels -----*/
01731 if (k == option_data->c) option_data->c++;
01732 }
01733
01734
01735 /*----- determine F value corresponding to alpha for lack of fit -----*/
01736 which = 2;
01737 q = option_data->flofmax;
01738 p = 1.0 - q;
01739 dfn = option_data->c - regmodel->p;
01740 dfd = xdata->rows - option_data->c;
01741 cdff (&which, &p, &q, &f, &dfn, &dfd, &status, &bound);
01742 if (status != 0) RA_error ("Error in calculating F - lack of fit ");
01743 option_data->flofmax = f;
01744
01745
01746 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 2080 of file 3dRegAna.c. References calc_matrices(), matrix_destroy(), matrix_initialize(), MAX_XVARS, p, and q.
02087 : 1/(X'X) for full model */ 02088 matrix * xtxinvxt_full, /* matrix: (1/(X'X))X' for full model */ 02089 matrix * x_base, /* extracted X matrix for baseline model */ 02090 matrix * xtxinvxt_base, /* matrix: (1/(X'X))X' for baseline model */ 02091 matrix * x_rdcd, /* extracted X matrices for reduced models */ 02092 matrix * xtxinvxt_rdcd /* matrix: (1/(X'X))X' for reduced models */ 02093 ) 02094 02095 { 02096 int zlist[MAX_XVARS]; /* list of parameters in constant model */ 02097 int ip; /* parameter index */ 02098 matrix xtxinv_temp; /* intermediate results */ 02099 02100 02101 /*----- Initialize matrix -----*/ 02102 matrix_initialize (&xtxinv_temp); 02103 02104 02105 /*----- Initialize list of parameters in the constant model -----*/ 02106 for (ip = 0; ip < MAX_XVARS; ip++) 02107 zlist[ip] = 0; 02108 02109 02110 /*----- Calculate constant matrices which will be needed later -----*/ 02111 calc_matrices (xdata, 1, zlist, x_base, &xtxinv_temp, xtxinvxt_base); 02112 calc_matrices (xdata, q, rlist, x_rdcd, &xtxinv_temp, xtxinvxt_rdcd); 02113 calc_matrices (xdata, p, flist, x_full, xtxinv_full, xtxinvxt_full); 02114 02115 02116 /*----- Destroy matrix -----*/ 02117 matrix_destroy (&xtxinv_temp); 02118 02119 } |
|
||||||||||||
|
Definition at line 943 of file 3dRegAna.c. References malloc, MAX_OBSERVATIONS, MAX_XVARS, and RA_error.
00948 {
00949 int ip; /* index */
00950
00951
00952 regmodel->p = 0;
00953 regmodel->q = 0;
00954
00955 option_data->datum = ILLEGAL_TYPE;
00956 strcpy (option_data->session, "./");
00957
00958 option_data->yname = NULL;
00959 option_data->first_dataset = NULL;
00960
00961 option_data->nx = 0;
00962 option_data->ny = 0;
00963 option_data->nz = 0;
00964 option_data->nxyz = 0;
00965
00966 option_data->diskspace = 0;
00967 option_data->workmem = 12;
00968 option_data->piece_size = 0;
00969 option_data->num_pieces = 0;
00970
00971 option_data->rms_min = 0.0;
00972 option_data->fdisp = -1.0;
00973
00974 option_data->levels = NULL;
00975 option_data->counts = NULL;
00976 option_data->c = 0;
00977 option_data->flofmax = -1;
00978
00979 option_data->numf = 0;
00980 option_data->numr = 0;
00981 option_data->numc = 0;
00982 option_data->numt = 0;
00983
00984 option_data->fcoef_filename = NULL;
00985 option_data->rcoef_filename = NULL;
00986 option_data->tcoef_filename = NULL;
00987
00988 option_data->numfiles = 0;
00989
00990 /*----- allocate memory for storing data file names -----*/
00991 option_data->yname
00992 = (char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
00993 for (ip = 0; ip < MAX_OBSERVATIONS; ip++)
00994 option_data->yname[ip] = NULL;
00995
00996
00997 /*----- allocate memory space and initialize pointers for filenames -----*/
00998 option_data->fcoef_filename =
00999 (char **) malloc (sizeof(char *) * MAX_XVARS);
01000 if (option_data->fcoef_filename == NULL)
01001 RA_error ("Unable to allocate memory for fcoef_filename");
01002
01003 option_data->rcoef_filename =
01004 (char **) malloc (sizeof(char *) * MAX_XVARS);
01005 if (option_data->rcoef_filename == NULL)
01006 RA_error ("Unable to allocate memory for rcoef_filename");
01007
01008 option_data->tcoef_filename =
01009 (char **) malloc (sizeof(char *) * MAX_XVARS);
01010 if (option_data->tcoef_filename == NULL)
01011 RA_error ("Unable to allocate memory for tcoef_filename");
01012
01013 for (ip = 0; ip < MAX_XVARS; ip++)
01014 {
01015 option_data->fcoef_filename[ip] = NULL;
01016 option_data->rcoef_filename[ip] = NULL;
01017 option_data->tcoef_filename[ip] = NULL;
01018 }
01019
01020
01021 /*----- initialize bucket dataset options -----*/
01022 option_data->bucket_filename = NULL;
01023 option_data->numbricks = -1;
01024 option_data->brick_type = NULL;
01025 option_data->brick_coef = NULL;
01026 option_data->brick_label = NULL;
01027
01028 }
|
|
||||||||||||||||||||||||
|
Definition at line 2014 of file 3dRegAna.c. References argc, break_into_pieces(), check_disk_space(), check_for_valid_inputs(), check_output_files(), check_temporary_files(), commandline, count_volumes_and_files(), get_dimensions(), get_inputs(), identify_repeats(), matrix_initialize(), PROGRAM_NAME, and tross_commandline().
02022 {
02023
02024 /*----- save command line for history notes -----*/
02025 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
02026
02027
02028 /*----- create independent variable data matrix -----*/
02029 matrix_initialize (xdata);
02030
02031
02032 /*----- get all operator inputs -----*/
02033 get_inputs (argc, argv, xdata, regmodel, option_data);
02034
02035
02036 /*----- use first data set to get data set dimensions -----*/
02037 option_data->first_dataset = option_data->yname[0];
02038 get_dimensions (option_data);
02039 printf ("Data set dimensions: nx = %d ny = %d nz = %d nxyz = %d \n",
02040 option_data->nx, option_data->ny, option_data->nz, option_data->nxyz);
02041
02042
02043 /*----- count the number of output volumes and files required -----*/
02044 count_volumes_and_files (regmodel, option_data);
02045
02046
02047 /*----- break problem into smaller pieces -----*/
02048 break_into_pieces (xdata->rows, option_data);
02049
02050
02051 /*----- identify repeat observations -----*/
02052 if (option_data->flofmax >= 0.0)
02053 identify_repeats (xdata, regmodel, option_data);
02054
02055
02056 /*----- check for valid inputs -----*/
02057 check_for_valid_inputs (xdata, regmodel, option_data);
02058
02059
02060 /*----- check whether any of the output files already exist -----*/
02061 check_output_files (option_data);
02062
02063
02064 /*----- check whether temporary files already exist -----*/
02065 check_temporary_files (xdata, regmodel, option_data);
02066
02067
02068 /*----- check whether there is sufficient disk space -----*/
02069 if (option_data->diskspace) check_disk_space (option_data);
02070
02071 }
|
|
||||||||||||
|
Definition at line 3184 of file 3dRegAna.c. References addto_args(), argc, calculate_results(), initialize_program(), machdep(), output_results(), PROGRAM_AUTHOR, PROGRAM_INITIAL, PROGRAM_LATEST, PROGRAM_NAME, and terminate_program().
03189 {
03190 RA_options option_data; /* user input options */
03191 matrix xdata; /* independent variable matrix */
03192 model regmodel; /* linear regression model */
03193 int piece_size; /* number of voxels in dataset piece */
03194
03195
03196 /*----- Identify software -----*/
03197 printf ("\n\n");
03198 printf ("Program: %s \n", PROGRAM_NAME);
03199 printf ("Author: %s \n", PROGRAM_AUTHOR);
03200 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
03201 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
03202 printf ("\n");
03203
03204 /*-- 22 Feb 1999: addto the arglist, if user wants to --*/
03205
03206 machdep() ;
03207 { int new_argc ; char ** new_argv ;
03208 addto_args( argc , argv , &new_argc , &new_argv ) ;
03209 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
03210 }
03211
03212 /*----- program initialization -----*/
03213 initialize_program (argc, argv, &xdata, ®model, &option_data);
03214
03215
03216 /*----- perform regression analysis -----*/
03217 calculate_results (xdata, ®model, &option_data);
03218
03219
03220 /*----- write requested output files -----*/
03221 output_results (xdata, ®model, &option_data);
03222
03223
03224 /*----- end of program -----*/
03225 terminate_program (&xdata, ®model, &option_data);
03226
03227 exit(0);
03228 }
|
|
||||||||||||||||
|
Definition at line 2873 of file 3dRegAna.c. References free, FUNC_THR_TYPE, malloc, MAX_NAME_LENGTH, MTEST, p, q, read_volume(), write_afni_data(), and write_bucket_data().
02879 {
02880 int p; /* number of parameters in full model */
02881 int q; /* number of parameters in reduced model */
02882 int n; /* number of data points */
02883 int nxyz; /* number of voxels */
02884 int ip, ix; /* parameter index */
02885 int numdof, dendof; /* numerator and denominator degrees of freedom */
02886 int piece_size; /* number of voxels in dataset piece */
02887 int num_pieces; /* dataset is divided into this many pieces */
02888 float * volume1 = NULL; /* volume data for 1st sub-brick */
02889 float * volume2 = NULL; /* volume data for 2nd sub-brick */
02890 char filename[MAX_NAME_LENGTH]; /* name for temporary data file */
02891
02892
02893 /*----- initialize local variables -----*/
02894 p = regmodel->p;
02895 q = regmodel->q;
02896 n = xdata.rows;
02897 nxyz = option_data->nxyz;
02898 piece_size = option_data->piece_size;
02899 num_pieces = option_data->num_pieces;
02900
02901
02902 /*----- allocate memory space for volume data -----*/
02903 volume1 = (float *) malloc (sizeof(float) * nxyz);
02904 MTEST (volume1);
02905 volume2 = (float *) malloc (sizeof(float) * nxyz);
02906 MTEST (volume2);
02907
02908
02909 /*----- write t-statistics data files -----*/
02910 if (option_data->numt > 0)
02911 {
02912 numdof = n - p;
02913 dendof = 0;
02914
02915 for (ip = 0; ip < p; ip++)
02916 {
02917 ix = regmodel->flist[ip];
02918
02919 if (option_data->tcoef_filename[ix] != NULL)
02920 {
02921 sprintf (filename, "coef.%d", ix);
02922 read_volume (filename, volume1, nxyz, piece_size, num_pieces);
02923
02924 sprintf (filename, "tcoef.%d", ix);
02925 read_volume (filename, volume2, nxyz, piece_size, num_pieces);
02926
02927 write_afni_data (option_data,
02928 option_data->tcoef_filename[ix],
02929 volume1, volume2,
02930 FUNC_TT_TYPE, numdof, dendof);
02931 }
02932 }
02933 }
02934
02935
02936 /*----- write R^2 data files -----*/
02937 if (option_data->numr > 0)
02938 {
02939 strcpy (filename, "rsqr");
02940 read_volume (filename, volume2, nxyz, piece_size, num_pieces);
02941
02942 for (ip = 0; ip < p; ip++)
02943 {
02944 ix = regmodel->flist[ip];
02945
02946 if (option_data->rcoef_filename[ix] != NULL)
02947 {
02948 sprintf (filename, "coef.%d", ix);
02949 read_volume (filename, volume1, nxyz, piece_size, num_pieces);
02950
02951 write_afni_data (option_data,
02952 option_data->rcoef_filename[ix],
02953 volume1, volume2,
02954 FUNC_THR_TYPE, 0, 0);
02955 }
02956 }
02957 }
02958
02959
02960 /*----- write F-statistics data files -----*/
02961 if (option_data->numf > 0)
02962 {
02963 strcpy (filename, "freg");
02964 read_volume (filename, volume2, nxyz, piece_size, num_pieces);
02965
02966 numdof = p - q;
02967 dendof = n - p;
02968
02969 for (ip = 0; ip < p; ip++)
02970 {
02971 ix = regmodel->flist[ip];
02972
02973 if (option_data->fcoef_filename[ix] != NULL)
02974 {
02975 sprintf (filename, "coef.%d", ix);
02976 read_volume (filename, volume1, nxyz, piece_size, num_pieces);
02977
02978 write_afni_data (option_data,
02979 option_data->fcoef_filename[ix],
02980 volume1, volume2,
02981 FUNC_FT_TYPE, numdof, dendof);
02982 }
02983 }
02984 }
02985
02986
02987 /*----- deallocate memory space for volume data -----*/
02988 free (volume1); volume1 = NULL;
02989 free (volume2); volume2 = NULL;
02990
02991
02992 /*----- write the bucket dataset -----*/
02993 if (option_data->numbricks > 0)
02994 write_bucket_data (xdata, regmodel, option_data);
02995
02996
02997 }
|
|
|
Definition at line 182 of file 3dRegAna.c. References PROGRAM_NAME.
00183 {
00184 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00185 exit(1);
00186 }
|
|
||||||||||||||||||||||||
|
Definition at line 347 of file 3dRegAna.c. References DOPEN, DSET_BRICK_FACTOR, DSET_BRICK_TYPE, DSET_PRINCIPAL_VALUE, EDIT_coerce_scale_type(), nz, SUB_POINTER, and THD_delete_3dim_dataset().
00355 {
00356 int iv; /* index number of intensity sub-brick */
00357 THD_3dim_dataset * dset=NULL; /* data set pointer */
00358 void * vfim = NULL; /* image data pointer */
00359 int nx, ny, nz, nxyz; /* data set dimensions in voxels */
00360
00361 nx = option_data->nx;
00362 ny = option_data->ny;
00363 nz = option_data->nz;
00364 nxyz = option_data->nxyz;
00365
00366
00367 /*----- read in the data -----*/
00368 DOPEN (dset, filename) ;
00369 iv = DSET_PRINCIPAL_VALUE(dset) ;
00370
00371 /*----- convert it to floats (in ffim) -----*/
00372 SUB_POINTER (dset, iv, fim_offset, vfim) ;
00373 EDIT_coerce_scale_type (piece_len, DSET_BRICK_FACTOR(dset,iv),
00374 DSET_BRICK_TYPE(dset,iv), vfim, /* input */
00375 MRI_float ,ffim ) ; /* output */
00376
00377 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00378 }
|
|
||||||||||||||||
|
Definition at line 416 of file 3dRegAna.c. References far, MAX_NAME_LENGTH, RA_error, and SUFFIX. Referenced by read_volume().
00422 {
00423 char sfilename[MAX_NAME_LENGTH]; /* piece file name */
00424 char message[MAX_NAME_LENGTH]; /* error message */
00425 FILE * far = NULL; /* floating point input file */
00426 int count; /* number of data items read from disk */
00427
00428
00429 /*----- input file name -----*/
00430 strcpy (sfilename, filename);
00431 strcat (sfilename, SUFFIX);
00432
00433 /*----- open temporary data file for input -----*/
00434 far = fopen (sfilename, "r");
00435 if (far == NULL)
00436 {
00437 sprintf (message, "Unable to open temporary file %s", sfilename);
00438 RA_error (message);
00439 }
00440
00441 /*----- read 3d data file -----*/
00442 count = fread (fin, sizeof(float), size, far);
00443 fclose (far);
00444
00445 /*----- error in reading file? -----*/
00446 if (count != size)
00447 {
00448 sprintf (message, "Error in reading temporary file %s", sfilename);
00449 RA_error (message);
00450 }
00451 }
|
|
||||||||||||||||||||||||
|
Definition at line 869 of file 3dRegAna.c. References MAX_NAME_LENGTH, and read_piece(). Referenced by output_results(), and write_bucket_data().
00877 {
00878 int piece; /* piece index */
00879 int piece_len; /* number of voxels in current piece */
00880 int fim_offset; /* array offset to current piece */
00881 char sfilename[MAX_NAME_LENGTH]; /* name for temporary data file */
00882
00883
00884 /*----- loop over the temporary data file pieces -----*/
00885 for (piece = 0; piece < num_pieces; piece++)
00886 {
00887
00888 /*----- current offset into volume -----*/
00889 fim_offset = piece * piece_size;
00890
00891 /*----- size of current piece -----*/
00892 if (piece < num_pieces-1)
00893 piece_len = piece_size;
00894 else
00895 piece_len = nxyz - fim_offset;
00896
00897 /*----- read in the piece data -----*/
00898 sprintf (sfilename, "%s.p%d", filename, piece);
00899 read_piece (sfilename, volume + fim_offset, piece_len);
00900
00901 } /* loop over pieces */
00902
00903 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 2128 of file 3dRegAna.c. References c, calc_coef(), calc_flof(), calc_freg(), calc_rsqr(), calc_sse(), calc_sspe(), calc_tcoef(), p, q, vector_create(), vector_destroy(), and vector_initialize().
02133 : 1/(X'X) for full model */ 02134 matrix xtxinvxt_full, /* matrix: (1/(X'X))X' for full model */ 02135 matrix x_base, /* extracted X matrix for baseline model */ 02136 matrix xtxinvxt_base, /* matrix: (1/(X'X))X' for baseline model */ 02137 matrix x_rdcd, /* extracted X matrix for reduced model */ 02138 matrix xtxinvxt_rdcd, /* matrix: (1/(X'X))X' for reduced model */ 02139 vector y, /* vector of measured data */ 02140 float rms_min, /* minimum rms error to reject zero model */ 02141 int * levels, /* indices for repeat observations */ 02142 int * counts, /* number of observations at each level */ 02143 int c, /* number of unique rows in ind. var. matrix */ 02144 float flofmax, /* max. allowed F-stat due to lack of fit */ 02145 float * flof, /* F-statistic for lack of fit */ 02146 vector * coef_full, /* regression parameters */ 02147 vector * scoef_full, /* std. devs. for regression parameters */ 02148 vector * tcoef_full, /* t-statistics for regression parameters */ 02149 float * freg, /* F-statistic for the full regression model */ 02150 float * rsqr /* coeff. of multiple determination R^2 */ 02151 ) 02152 02153 { 02154 float sse_base; /* error sum of squares, baseline model */ 02155 float sse_rdcd; /* error sum of squares, reduced model */ 02156 float sse_full; /* error sum of squares, full model */ 02157 float sspe; /* pure error sum of squares */ 02158 vector coef_temp; /* intermediate results */ 02159 02160 02161 /*----- Initialization -----*/ 02162 vector_initialize (&coef_temp); 02163 02164 02165 /*----- Calculate regression coefficients for baseline model -----*/ 02166 calc_coef (xtxinvxt_base, y, &coef_temp); 02167 02168 02169 /*----- Calculate the error sum of squares for the baseline model -----*/ 02170 sse_base = calc_sse (x_base, coef_temp, y); 02171 02172 02173 /*----- Stop here if variation about baseline is sufficiently low -----*/ 02174 if (sqrt(sse_base/N) < rms_min) 02175 { 02176 vector_create (p, coef_full); 02177 vector_create (p, scoef_full); 02178 vector_create (p, tcoef_full); 02179 *freg = 0.0; 02180 *rsqr = 0.0; 02181 vector_destroy (&coef_temp); 02182 return; 02183 } 02184 02185 02186 /*----- Calculate regression coefficients for the full model -----*/ 02187 calc_coef (xtxinvxt_full, y, coef_full); 02188 02189 02190 /*----- Calculate the error sum of squares for the full model -----*/ 02191 sse_full = calc_sse (x_full, *coef_full, y); 02192 02193 02194 /*----- Calculate t-statistics for the regression coefficients -----*/ 02195 calc_tcoef (N, p, sse_full, xtxinv_full, 02196 *coef_full, scoef_full, tcoef_full); 02197 02198 02199 /*----- Test for lack of fit -----*/ 02200 if (flofmax > 0.0) 02201 { 02202 /*----- Calculate the pure error sum of squares -----*/ 02203 sspe = calc_sspe (y, levels, counts, c); 02204 02205 /*----- Calculate F-statistic for lack of fit -----*/ 02206 *flof = calc_flof (N, p, c, sse_full, sspe); 02207 02208 if (*flof > flofmax) 02209 { 02210 vector_create (p, coef_full); 02211 vector_create (p, scoef_full); 02212 vector_create (p, tcoef_full); 02213 *freg = 0.0; 02214 *rsqr = 0.0; 02215 vector_destroy (&coef_temp); 02216 return; 02217 } 02218 } 02219 else 02220 *flof = -1.0; 02221 02222 02223 /*----- Calculate regression coefficients for reduced model -----*/ 02224 calc_coef (xtxinvxt_rdcd, y, &coef_temp); 02225 02226 02227 /*----- Calculate the error sum of squares for the reduced model -----*/ 02228 sse_rdcd = calc_sse (x_rdcd, coef_temp, y); 02229 02230 02231 /*----- Calculate F-statistic for significance of the regression -----*/ 02232 *freg = calc_freg (N, p, q, sse_full, sse_rdcd); 02233 02234 02235 /*----- Calculate coefficient of multiple determination R^2 -----*/ 02236 *rsqr = calc_rsqr (sse_full, sse_base); 02237 02238 02239 /*----- Dispose of vector -----*/ 02240 vector_destroy (&coef_temp); 02241 02242 } |
|
||||||||||||||||||||||||||||
|
Definition at line 704 of file 3dRegAna.c. References MAX_NAME_LENGTH, MAX_XVARS, and write_piece(). Referenced by calculate_results().
00714 {
00715 int ip; /* parameter index */
00716 char filename[MAX_NAME_LENGTH]; /* name for temporary data file */
00717
00718
00719 /*----- save piece containing F-statistics -----*/
00720 if (freg_piece != NULL)
00721 {
00722 sprintf (filename, "freg.p%d", piece);
00723 write_piece (filename, freg_piece, piece_len);
00724 }
00725
00726
00727 /*----- save piece containing R^2 -----*/
00728 if (rsqr_piece != NULL)
00729 {
00730 sprintf (filename, "rsqr.p%d", piece);
00731 write_piece (filename, rsqr_piece, piece_len);
00732 }
00733
00734
00735 /*----- save pieces containing regression coefficients -----*/
00736 if (coef_piece != NULL)
00737 {
00738 for (ip = 0; ip < MAX_XVARS; ip++)
00739 if (coef_piece[ip] != NULL)
00740 {
00741 sprintf (filename, "coef.%d.p%d", ip, piece);
00742 write_piece (filename, coef_piece[ip], piece_len);
00743 }
00744 }
00745
00746
00747 /*----- save pieces containing t-statistics -----*/
00748 if (tcoef_piece != NULL)
00749 {
00750 for (ip = 0; ip < MAX_XVARS; ip++)
00751 if (tcoef_piece[ip] != NULL)
00752 {
00753 sprintf (filename, "tcoef.%d.p%d", ip, piece);
00754 write_piece (filename, tcoef_piece[ip], piece_len);
00755 }
00756 }
00757
00758 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 2251 of file 3dRegAna.c.
02268 {
02269 int ip; /* parameter index */
02270 int ix; /* x-variable index */
02271
02272
02273 /*----- save regression results into piece data -----*/
02274 if (freg_piece != NULL) freg_piece[iv] = freg;
02275 if (rsqr_piece != NULL) rsqr_piece[iv] = rsqr;
02276
02277
02278 /*----- save regression parameter estimates -----*/
02279 for (ip = 0; ip < regmodel->p; ip++)
02280 {
02281 ix = regmodel->flist[ip];
02282
02283 if (coef_piece[ix] != NULL) coef_piece[ix][iv] = coef.elts[ip];
02284
02285 if (tcoef_piece[ix] != NULL) tcoef_piece[ix][iv] = tcoef.elts[ip];
02286
02287 }
02288
02289
02290 /*----- if so requested, display results for this voxel -----*/
02291 if ((fdisp >= 0.0) && (freg >= fdisp))
02292 {
02293 printf ("\n\nVoxel #%d: \n", iv);
02294 printf ("\nY data: \n");
02295 for (ip = 0; ip < y.dim; ip++)
02296 printf ("Y[%d] = %f \n", ip, y.elts[ip]);
02297
02298 if (flof >= 0.0) printf ("\nF lack of fit = %f \n", flof);
02299 printf ("\nF regression = %f \n", freg);
02300 printf ("R-squared = %f \n", rsqr);
02301
02302 printf ("\nRegression Coefficients: \n");
02303 for (ip = 0; ip < coef.dim; ip++)
02304 {
02305 ix = regmodel->flist[ip];
02306 printf ("b[%d] = %f ", ix, coef.elts[ip]);
02307 printf ("t[%d] = %f \n", ix, tcoef.elts[ip]);
02308 }
02309
02310 }
02311
02312 }
|
|
|
Definition at line 1037 of file 3dRegAna.c. Referenced by get_inputs().
01041 {
01042 int i, j, temp; /* model variable index numbers */
01043
01044
01045 /*----- sort full model indices into ascending order -----*/
01046 for (i = 0; i < regmodel->p - 1; i++)
01047 for (j = i+1; j < regmodel->p; j++)
01048 if (regmodel->flist[i] > regmodel->flist[j])
01049 {
01050 temp = regmodel->flist[i];
01051 regmodel->flist[i] = regmodel->flist[j];
01052 regmodel->flist[j] = temp;
01053 }
01054 else if (regmodel->flist[i] == regmodel->flist[j])
01055 RA_error ("Duplicate variable indices in model definition");
01056
01057
01058 /*----- sort reduced model indices into ascending order -----*/
01059 for (i = 0; i < regmodel->q - 1; i++)
01060 for (j = i+1; j < regmodel->q; j++)
01061 if (regmodel->rlist[i] > regmodel->rlist[j])
01062 {
01063 temp = regmodel->rlist[i];
01064 regmodel->rlist[i] = regmodel->rlist[j];
01065 regmodel->rlist[j] = temp;
01066 }
01067
01068 }
|
|
||||||||||||||||
|
Definition at line 3006 of file 3dRegAna.c. References delete_volume(), free, matrix_destroy(), MAX_NAME_LENGTH, MAX_OBSERVATIONS, MAX_XVARS, and p.
03012 {
03013 int p; /* number of parameters in full model */
03014 int ip, ix; /* parameter index */
03015 int ib; /* sub-brick index */
03016 int nxyz; /* number of voxels */
03017 int piece_size; /* number of voxels in dataset piece */
03018 int num_pieces; /* dataset is divided into this many pieces */
03019 char filename[MAX_NAME_LENGTH]; /* name for temporary data file */
03020
03021
03022 /*----- initialize local variables -----*/
03023 p = regmodel->p;
03024 nxyz = option_data->nxyz;
03025 piece_size = option_data->piece_size;
03026 num_pieces = option_data->num_pieces;
03027
03028
03029 /*----- delete F-statistics data files -----*/
03030 if (option_data->numf > 0)
03031 {
03032 strcpy (filename, "freg");
03033 delete_volume (filename, nxyz, piece_size, num_pieces);
03034 }
03035
03036
03037 /*----- delete R^2 data files -----*/
03038 if (option_data->numr > 0)
03039 {
03040 strcpy (filename, "rsqr");
03041 delete_volume (filename, nxyz, piece_size, num_pieces);
03042 }
03043
03044
03045 /*----- delete t-statistics data files -----*/
03046 if (option_data->numt > 0)
03047 {
03048 for (ip = 0; ip < p; ip++)
03049 {
03050 ix = regmodel->flist[ip];
03051 sprintf (filename, "tcoef.%d", ix);
03052 delete_volume (filename, nxyz, piece_size, num_pieces);
03053 }
03054 }
03055
03056
03057 /*----- delete regression coefficients data files -----*/
03058 if (option_data->numc > 0)
03059 {
03060 for (ip = 0; ip < p; ip++)
03061 {
03062 ix = regmodel->flist[ip];
03063 sprintf (filename, "coef.%d", ix);
03064 delete_volume (filename, nxyz, piece_size, num_pieces);
03065 }
03066 }
03067
03068
03069 /*----- dispose of data matrix -----*/
03070 matrix_destroy (xdata);
03071
03072
03073 /*----- deallocate memory -----*/
03074 if (option_data->yname != NULL)
03075 {
03076 for (ip = 0; ip < MAX_OBSERVATIONS; ip++)
03077 {
03078 if (option_data->yname[ip] != NULL)
03079 {
03080 free (option_data->yname[ip]);
03081 option_data->yname[ip] = NULL;
03082 }
03083 }
03084 free (option_data->yname);
03085 option_data->yname = NULL;
03086 }
03087
03088 if (option_data->fcoef_filename != NULL)
03089 {
03090 for (ip = 0; ip < MAX_XVARS; ip++)
03091 {
03092 if (option_data->fcoef_filename[ip] != NULL)
03093 {
03094 free (option_data->fcoef_filename[ip]);
03095 option_data->fcoef_filename[ip] = NULL;
03096 }
03097 }
03098 free (option_data->fcoef_filename);
03099 option_data->fcoef_filename = NULL;
03100 }
03101
03102 if (option_data->rcoef_filename != NULL)
03103 {
03104 for (ip = 0; ip < MAX_XVARS; ip++)
03105 {
03106 if (option_data->rcoef_filename[ip] != NULL)
03107 {
03108 free (option_data->rcoef_filename[ip]);
03109 option_data->rcoef_filename[ip] = NULL;
03110 }
03111 }
03112 free (option_data->rcoef_filename);
03113 option_data->rcoef_filename = NULL;
03114 }
03115
03116 if (option_data->tcoef_filename != NULL)
03117 {
03118 for (ip = 0; ip < MAX_XVARS; ip++)
03119 {
03120 if (option_data->tcoef_filename[ip] != NULL)
03121 {
03122 free (option_data->tcoef_filename[ip]);
03123 option_data->tcoef_filename[ip] = NULL;
03124 }
03125 }
03126 free (option_data->tcoef_filename);
03127 option_data->tcoef_filename = NULL;
03128 }
03129
03130 if (option_data->levels != NULL)
03131 {
03132 free (option_data->levels);
03133 option_data->levels = NULL;
03134 }
03135
03136 if (option_data->counts != NULL)
03137 {
03138 free (option_data->counts);
03139 option_data->counts = NULL;
03140 }
03141
03142
03143 if (option_data->bucket_filename != NULL)
03144 {
03145 free (option_data->bucket_filename);
03146 option_data->bucket_filename = NULL;
03147 }
03148
03149 if (option_data->brick_type != NULL)
03150 {
03151 free (option_data->brick_type);
03152 option_data->brick_type = NULL;
03153 }
03154
03155 if (option_data->brick_coef != NULL)
03156 {
03157 free (option_data->brick_coef);
03158 option_data->brick_coef = NULL;
03159 }
03160
03161 if (option_data->brick_label != NULL)
03162 {
03163 for (ib = 0; ib < option_data->numbricks; ib++)
03164 {
03165 if (option_data->brick_label[ib] != NULL)
03166 {
03167 free (option_data->brick_label[ib]);
03168 option_data->brick_label[ib] = NULL;
03169 }
03170 }
03171 free (option_data->brick_label);
03172 option_data->brick_label = NULL;
03173 }
03174
03175 }
|
|
||||||||||||||||||||||||||||||||
|
Definition at line 2532 of file 3dRegAna.c. References ADN_brick_fac, ADN_datum_array, ADN_directory_name, ADN_func_type, ADN_label1, ADN_malloc_type, ADN_none, ADN_nvals, ADN_prefix, ADN_self_name, ADN_stat_aux, ADN_type, commandline, DATABLOCK_MEM_MALLOC, THD_3dim_dataset::dblk, THD_datablock::diskptr, DSET_BRICK, DSET_BRIKNAME, EDIT_coerce_autoscale_new(), EDIT_dset_items(), EDIT_empty_copy(), FUNC_FT_SCALE_SHORT, FUNC_THR_SCALE_SHORT, FUNC_THR_TYPE, FUNC_TT_SCALE_SHORT, GEN_FUNC_TYPE, HEAD_FUNC_TYPE, THD_diskptr::header_name, ISHEAD, ISVALID_3DIM_DATASET, malloc, MAX_STAT_AUX, mri_datum_size(), mri_fix_data_pointer(), RA_error, THD_delete_3dim_dataset(), THD_is_file(), THD_load_statistics(), THD_open_dataset(), THD_write_3dim_dataset(), top, tross_Append_History(), and tross_multi_Append_History().
02542 {
02543 int nxyz; /* number of voxels */
02544 int ii; /* voxel index */
02545 THD_3dim_dataset * dset=NULL; /* input afni data set pointer */
02546 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
02547 int ierror; /* number of errors in editing data */
02548 int ibuf[32]; /* integer buffer */
02549 float fbuf[MAX_STAT_AUX]; /* float buffer */
02550 float fimfac; /* scale factor for short data */
02551 int output_datum; /* data type for output data */
02552 short * tsp = NULL; /* 2nd sub-brick data pointer */
02553 void * vdif = NULL; /* 1st sub-brick data pointer */
02554 float top, bot, func_scale_short; /* parameters for scaling data */
02555 int top_ss, bot_ss; /* 2nd sub-brick value limits */
02556 char label[80]; /* label for output file history */
02557
02558
02559 /*----- initialize local variables -----*/
02560 nxyz = option_data->nxyz;
02561
02562 /*----- read first dataset -----*/
02563 dset = THD_open_dataset (option_data->first_dataset) ;
02564 if( ! ISVALID_3DIM_DATASET(dset) ){
02565 fprintf(stderr,"*** Unable to open dataset file %s\n",
02566 option_data->first_dataset); exit(1) ;
02567 }
02568
02569
02570 /*-- make an empty copy of this dataset, for eventual output --*/
02571 new_dset = EDIT_empty_copy( dset ) ;
02572
02573
02574 /*----- Record history of dataset -----*/
02575
02576 sprintf (label, "Output prefix: %s", filename);
02577 if( commandline != NULL )
02578 tross_multi_Append_History( new_dset , commandline,label,NULL ) ;
02579 else
02580 tross_Append_History ( new_dset, label);
02581
02582
02583 output_datum = MRI_short ;
02584 ibuf[0] = output_datum ; ibuf[1] = MRI_short ;
02585
02586
02587 ierror = EDIT_dset_items( new_dset ,
02588 ADN_prefix , filename ,
02589 ADN_label1 , filename ,
02590 ADN_directory_name , option_data->session ,
02591 ADN_self_name , filename ,
02592 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE :
02593 GEN_FUNC_TYPE ,
02594 ADN_func_type , func_type ,
02595 ADN_nvals , FUNC_nvals[func_type] ,
02596 ADN_datum_array , ibuf ,
02597 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
02598 ADN_none ) ;
02599
02600 if( ierror > 0 ){
02601 fprintf(stderr,
02602 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
02603 exit(1) ;
02604 }
02605
02606 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
02607 fprintf(stderr,
02608 "*** Output dataset file %s already exists--cannot continue!\a\n",
02609 new_dset->dblk->diskptr->header_name ) ;
02610 exit(1) ;
02611 }
02612
02613 /*----- deleting exemplar dataset -----*/
02614 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
02615
02616
02617 /*----- allocate memory for output data -----*/
02618 vdif = (void *) malloc( mri_datum_size(output_datum) * nxyz ) ;
02619 tsp = (short *) malloc( sizeof(short) * nxyz ) ;
02620
02621 /*----- attach bricks to new data set -----*/
02622 mri_fix_data_pointer (vdif, DSET_BRICK(new_dset,0));
02623 mri_fix_data_pointer (tsp, DSET_BRICK(new_dset,1));
02624
02625
02626 /*----- convert data type to output specification -----*/
02627 fimfac = EDIT_coerce_autoscale_new (nxyz,
02628 MRI_float, ffim,
02629 output_datum, vdif);
02630 if (fimfac != 0.0) fimfac = 1.0 / fimfac;
02631
02632
02633 top_ss = 32700;
02634
02635 if (func_type == FUNC_THR_TYPE) /* threshold */
02636 {
02637 func_scale_short = FUNC_THR_SCALE_SHORT;
02638 bot_ss = 0;
02639 }
02640 else if (func_type == FUNC_TT_TYPE) /* t-statistic */
02641 {
02642 func_scale_short = FUNC_TT_SCALE_SHORT;
02643 bot_ss = -top_ss;
02644 }
02645 else if (func_type == FUNC_FT_TYPE) /* F-statistic */
02646 {
02647 func_scale_short = FUNC_FT_SCALE_SHORT;
02648 bot_ss = 0;
02649 }
02650 else
02651 RA_error ("Illegal ouput dataset function type");
02652
02653 top = top_ss / func_scale_short;
02654 bot = bot_ss / func_scale_short;
02655
02656
02657 for (ii = 0; ii < nxyz; ii++)
02658 {
02659 if (ftr[ii] > top)
02660 tsp[ii] = top_ss;
02661 else if (ftr[ii] < bot)
02662 tsp[ii] = bot_ss;
02663 else
02664 tsp[ii] = (short) (func_scale_short * ftr[ii] + 0.5);
02665 }
02666
02667
02668 /*----- write afni data set -----*/
02669 if (func_type == FUNC_THR_TYPE) /* threshold */
02670 printf("----- Writing `fith' dataset ");
02671 else if (func_type == FUNC_TT_TYPE) /* t-statistic */
02672 printf("----- Writing `fitt' dataset ");
02673 else if (func_type == FUNC_FT_TYPE) /* F-statistic */
02674 printf("----- Writing `fift' dataset ");
02675
02676 printf("into %s\n", DSET_BRIKNAME(new_dset) ) ;
02677
02678 fbuf[0] = numdof;
02679 fbuf[1] = dendof;
02680 for( ii=2 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
02681 (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
02682
02683 fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
02684 fbuf[1] = 1.0 / func_scale_short ;
02685 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
02686
02687 THD_load_statistics( new_dset ) ;
02688 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
02689
02690
02691 /*----- deallocate memory -----*/
02692 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
02693
02694 }
|
|
||||||||||||||||
|
Definition at line 2703 of file 3dRegAna.c. References ADN_directory_name, ADN_func_type, ADN_malloc_type, ADN_none, ADN_ntt, ADN_nvals, ADN_prefix, ADN_type, commandline, DATABLOCK_MEM_MALLOC, delete_volume(), DSET_BRIKNAME, DSET_HEADNAME, EDIT_BRICK_FACTOR, EDIT_BRICK_LABEL, EDIT_BRICK_TO_FIFT, EDIT_BRICK_TO_FITT, EDIT_coerce_autoscale_new(), EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), free, FUNC_BUCK_TYPE, FUNC_FIM_TYPE, FUNC_THR_TYPE, HEAD_FUNC_TYPE, malloc, MAX_NAME_LENGTH, MTEST, p, q, read_volume(), THD_delete_3dim_dataset(), THD_is_file(), THD_load_statistics(), THD_open_dataset(), THD_write_3dim_dataset(), and tross_Append_History().
02709 {
02710 const float EPSILON = 1.0e-10;
02711 int p; /* number of parameters in full model */
02712 int q; /* number of parameters in reduced model */
02713 int n; /* number of data points */
02714 int nxyz; /* number of voxels */
02715 THD_3dim_dataset * old_dset = NULL; /* prototype dataset */
02716 THD_3dim_dataset * new_dset = NULL; /* output bucket dataset */
02717 char * output_prefix; /* prefix name for bucket dataset */
02718 char * output_session; /* directory for bucket dataset */
02719 int nbricks, ib; /* number of sub-bricks in bucket dataset */
02720 short ** bar = NULL; /* bar[ib] points to data for sub-brick #ib */
02721 float factor; /* factor is new scale factor for sub-brick #ib */
02722 int brick_type; /* indicates statistical type of sub-brick */
02723 int brick_coef; /* regression coefficient index for sub-brick */
02724 char * brick_label; /* character string label for sub-brick */
02725 int ierror; /* number of errors in editing data */
02726 char filename[MAX_NAME_LENGTH]; /* name for temporary data file */
02727 int piece_size; /* number of voxels in dataset piece */
02728 int num_pieces; /* dataset is divided into this many pieces */
02729 float * volume = NULL; /* volume of floating point data */
02730 char label[80]; /* label for output file history */
02731
02732
02733 /*----- initialize local variables -----*/
02734 p = regmodel->p;
02735 q = regmodel->q;
02736 n = xdata.rows;
02737 nxyz = option_data->nxyz;
02738 piece_size = option_data->piece_size;
02739 num_pieces = option_data->num_pieces;
02740 nbricks = option_data->numbricks;
02741 output_prefix = option_data->bucket_filename;
02742 output_session = option_data->session;
02743
02744
02745 /*----- allocate memory -----*/
02746 volume = (float *) malloc (sizeof(float) * nxyz);
02747 MTEST (volume);
02748 bar = (short **) malloc (sizeof(short *) * nbricks);
02749 MTEST (bar);
02750
02751
02752 /*----- read first dataset -----*/
02753 old_dset = THD_open_dataset (option_data->first_dataset) ;
02754
02755
02756 /*-- make an empty copy of this dataset, for eventual output --*/
02757 new_dset = EDIT_empty_copy (old_dset);
02758
02759
02760 /*----- Record history of dataset -----*/
02761 if( commandline != NULL )
02762 tross_Append_History( new_dset , commandline ) ;
02763 sprintf (label, "Output prefix: %s", output_prefix);
02764 tross_Append_History ( new_dset, label);
02765
02766
02767 /*----- Modify some structural properties. Note that the nbricks
02768 just make empty sub-bricks, without any data attached. -----*/
02769 ierror = EDIT_dset_items (new_dset,
02770 ADN_prefix, output_prefix,
02771 ADN_directory_name, output_session,
02772 ADN_type, HEAD_FUNC_TYPE,
02773 ADN_func_type, FUNC_BUCK_TYPE,
02774 ADN_ntt, 0, /* no time */
02775 ADN_nvals, nbricks,
02776 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
02777 ADN_none ) ;
02778
02779 if( ierror > 0 )
02780 {
02781 fprintf(stderr,
02782 "*** %d errors in attempting to create output dataset!\n",
02783 ierror);
02784 exit(1);
02785 }
02786
02787 if (THD_is_file(DSET_HEADNAME(new_dset)))
02788 {
02789 fprintf(stderr,
02790 "*** Output dataset file %s already exists--cannot continue!\n",
02791 DSET_HEADNAME(new_dset));
02792 exit(1);
02793 }
02794
02795
02796 /*----- deleting exemplar dataset -----*/
02797 THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ;
02798
02799
02800 /*----- loop over new sub-brick index, attach data array with
02801 EDIT_substitute_brick then put some strings into the labels and
02802 keywords, and modify the sub-brick scaling factor -----*/
02803 for (ib = 0; ib < nbricks; ib++)
02804 {
02805 /*----- get information about this sub-brick -----*/
02806 brick_type = option_data->brick_type[ib];
02807 brick_coef = option_data->brick_coef[ib];
02808 brick_label = option_data->brick_label[ib];
02809
02810 if (brick_type == FUNC_FIM_TYPE)
02811 {
02812 sprintf (filename, "coef.%d", brick_coef);
02813 }
02814 else if (brick_type == FUNC_THR_TYPE)
02815 {
02816 sprintf (filename, "rsqr");
02817 }
02818 else if (brick_type == FUNC_TT_TYPE)
02819 {
02820 sprintf (filename, "tcoef.%d", brick_coef);
02821 EDIT_BRICK_TO_FITT (new_dset, ib, n-p);
02822 }
02823 else if (brick_type == FUNC_FT_TYPE)
02824 {
02825 sprintf (filename, "freg");
02826 EDIT_BRICK_TO_FIFT (new_dset, ib, p-q, n-p);
02827 }
02828
02829 /*----- allocate memory for output sub-brick -----*/
02830 bar[ib] = (short *) malloc (sizeof(short) * nxyz);
02831 MTEST (bar[ib]);
02832
02833 read_volume (filename, volume, nxyz, piece_size, num_pieces);
02834 delete_volume (filename, nxyz, piece_size, num_pieces);
02835
02836 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
02837 MRI_short, bar[ib]);
02838
02839 if (factor < EPSILON) factor = 0.0;
02840 else factor = 1.0 / factor;
02841
02842 /*----- edit the sub-brick -----*/
02843 EDIT_BRICK_LABEL (new_dset, ib, brick_label);
02844 EDIT_BRICK_FACTOR (new_dset, ib, factor);
02845
02846
02847 /*----- attach bar[ib] to be sub-brick #ib -----*/
02848 EDIT_substitute_brick (new_dset, ib, MRI_short, bar[ib]);
02849
02850 }
02851
02852
02853 /*----- write bucket data set -----*/
02854 THD_load_statistics (new_dset);
02855 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
02856 fprintf(stderr,"----- Wrote bucket dataset ");
02857 fprintf(stderr,"into %s\n", DSET_BRIKNAME(new_dset));
02858
02859
02860 /*----- deallocate memory -----*/
02861 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
02862 free (volume); volume = NULL;
02863
02864 }
|
|
||||||||||||||||
|
Definition at line 461 of file 3dRegAna.c. References far, fout, MAX_NAME_LENGTH, RA_error, and SUFFIX. Referenced by save_pieces().
00467 {
00468 char sfilename[MAX_NAME_LENGTH]; /* piece file name */
00469 char message[MAX_NAME_LENGTH]; /* error message */
00470 FILE * far = NULL; /* floating point output file */
00471 int count; /* number of data items written to disk */
00472
00473
00474 /*----- output file name -----*/
00475 strcpy (sfilename, filename);
00476 strcat (sfilename, SUFFIX);
00477
00478 /*----- first, see if file already exists -----*/
00479 far = fopen (sfilename, "r");
00480 if (far != NULL)
00481 {
00482 sprintf (message, "Temporary file %s already exists. ", sfilename);
00483 RA_error (message);
00484 }
00485
00486 /*----- open temporary data file for output -----*/
00487 far = fopen (sfilename, "w");
00488 if (far == NULL)
00489 {
00490 sprintf (message, "Unable to write temporary file %s ", sfilename);
00491 RA_error (message);
00492 }
00493
00494 /*----- write 3d data set -----*/
00495 count = fwrite (fout, sizeof(float), size, far);
00496 fclose (far);
00497
00498 /*----- error in writing file? -----*/
00499 if (count != size)
00500 {
00501 sprintf (message, "Error in writing temporary file %s ", sfilename);
00502 RA_error (message);
00503 }
00504
00505 }
|
Variable Documentation
|
|
Definition at line 120 of file 3dRegAna.c. Referenced by initialize_program(), write_afni_data(), and write_bucket_data(). |
