Doxygen Source Code Documentation
3dFDR.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Data Structures | |
struct | voxel |
Defines | |
#define | PROGRAM_NAME "3dFDR" |
#define | PROGRAM_AUTHOR "B. Douglas Ward" |
#define | PROGRAM_INITIAL "31 January 2002" |
#define | PROGRAM_LATEST "31 January 2002" |
#define | FDR_MAX_LL 10000 |
#define | DOPEN(ds, name) |
#define | SUB_POINTER(ds, vv, ind, ptr) |
#define | MTEST(ptr) |
Typedefs | |
typedef voxel | voxel |
Functions | |
void | FDR_error (char *message) |
void | FDR_Syntax (void) |
void | read_options (int argc, char *argv[]) |
float * | read_time_series (char *ts_filename, int *ts_length) |
void | check_one_output_file (THD_3dim_dataset *dset_time, char *filename) |
void * | initialize_program (int argc, char *argv[]) |
voxel * | create_voxel () |
voxel * | add_voxel (voxel *new_voxel, voxel *head_voxel) |
voxel * | new_voxel (int ixyz, float pvalue, voxel *head_voxel) |
void | delete_all_voxels () |
void | save_all_voxels (float *far) |
void | process_volume (float *ffim, int statcode, float *stataux) |
void | process_1ddata () |
float | EDIT_coerce_autoscale_new (int nxyz, int itype, void *ivol, int otype, void *ovol) |
void | process_subbrick (THD_3dim_dataset *dset, int ibrick) |
THD_3dim_dataset * | process_dataset () |
void | output_results (THD_3dim_dataset *new_dset) |
int | main (int argc, char *argv[]) |
Variables | |
int | FDR_quiet = 0 |
int | FDR_list = 0 |
float | FDR_mask_thr = 1.0 |
float | FDR_cn = 1.0 |
int | FDR_nxyz = 0 |
int | FDR_nthr = 0 |
char * | FDR_input_filename = NULL |
char * | FDR_input1D_filename = NULL |
char * | FDR_mask_filename = NULL |
char * | FDR_output_prefix = NULL |
byte * | FDR_mask = NULL |
float * | FDR_input1D_data = NULL |
voxel * | FDR_head_voxel [FDR_MAX_LL] |
char * | commandline = NULL |
THD_3dim_dataset * | FDR_dset = NULL |
Define Documentation
|
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) ; } \ 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 51 of file 3dFDR.c. Referenced by delete_all_voxels(), initialize_program(), process_volume(), and save_all_voxels(). |
|
Value: if((ptr)==NULL) \ ( FDR_error ("Cannot allocate memory") ) |
|
Definition at line 20 of file 3dFDR.c. Referenced by main(). |
|
Definition at line 21 of file 3dFDR.c. Referenced by main(). |
|
Definition at line 22 of file 3dFDR.c. Referenced by main(). |
|
Definition at line 19 of file 3dFDR.c. Referenced by FDR_error(), initialize_program(), and main(). |
|
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) |
Typedef Documentation
|
Definition at line 471 of file vp_extract.c. |
Function Documentation
|
Definition at line 658 of file 3dFDR.c. References new_voxel(), voxel::next_voxel, and voxel::pvalue. Referenced by new_voxel().
00659 { 00660 voxel * voxel_ptr = NULL; 00661 00662 if ((head_voxel == NULL) || (new_voxel->pvalue >= head_voxel->pvalue)) 00663 { 00664 new_voxel->next_voxel = head_voxel; 00665 head_voxel = new_voxel; 00666 } 00667 00668 else 00669 { 00670 voxel_ptr = head_voxel; 00671 00672 while ((voxel_ptr->next_voxel != NULL) && 00673 (new_voxel->pvalue < voxel_ptr->next_voxel->pvalue)) 00674 voxel_ptr = voxel_ptr->next_voxel; 00675 00676 new_voxel->next_voxel = voxel_ptr->next_voxel; 00677 voxel_ptr->next_voxel = new_voxel; 00678 } 00679 00680 return (head_voxel); 00681 } |
|
Definition at line 388 of file 3dFDR.c. References ADN_label1, ADN_none, ADN_prefix, ADN_self_name, ADN_type, THD_3dim_dataset::dblk, THD_datablock::diskptr, EDIT_dset_items(), EDIT_empty_copy(), FDR_error(), GEN_FUNC_TYPE, HEAD_FUNC_TYPE, THD_diskptr::header_name, ISHEAD, THD_delete_3dim_dataset(), THD_is_file(), and THD_MAX_NAME.
00393 { 00394 char message[THD_MAX_NAME]; /* error message */ 00395 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */ 00396 int ierror; /* number of errors in editing data */ 00397 00398 00399 /*----- make an empty copy of input dataset -----*/ 00400 new_dset = EDIT_empty_copy( dset_time ) ; 00401 00402 00403 ierror = EDIT_dset_items( new_dset , 00404 ADN_prefix , filename , 00405 ADN_label1 , filename , 00406 ADN_self_name , filename , 00407 ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE : 00408 GEN_FUNC_TYPE , 00409 ADN_none ) ; 00410 00411 if( ierror > 0 ) 00412 { 00413 sprintf (message, 00414 "*** %d errors in attempting to create output dataset!\n", 00415 ierror); 00416 FDR_error (message); 00417 } 00418 00419 if( THD_is_file(new_dset->dblk->diskptr->header_name) ) 00420 { 00421 sprintf (message, 00422 "Output dataset file %s already exists " 00423 " -- cannot continue! ", 00424 new_dset->dblk->diskptr->header_name); 00425 FDR_error (message); 00426 } 00427 00428 /*----- deallocate memory -----*/ 00429 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ; 00430 00431 } |
|
Definition at line 637 of file 3dFDR.c. References voxel::ixyz, malloc, MTEST, voxel::next_voxel, and voxel::pvalue. Referenced by new_voxel().
|
|
Definition at line 711 of file 3dFDR.c. References FDR_MAX_LL, free, and voxel::next_voxel. Referenced by process_volume().
00712 { 00713 int ibin; 00714 voxel * voxel_ptr = NULL; /* pointer to current voxel */ 00715 voxel * next_voxel = NULL; /* pointer to next voxel */ 00716 00717 00718 for (ibin = 0; ibin < FDR_MAX_LL; ibin++) 00719 { 00720 voxel_ptr = FDR_head_voxel[ibin]; 00721 while (voxel_ptr != NULL) 00722 { 00723 next_voxel = voxel_ptr->next_voxel; 00724 free (voxel_ptr); 00725 voxel_ptr = next_voxel; 00726 } 00727 FDR_head_voxel[ibin] = NULL; 00728 } 00729 00730 } |
|
compute start time of this timeseries * Definition at line 954 of file 3dFDR.c. References EDIT_coerce_scale_type(), MCW_vol_amax(), MRI_IS_INT_TYPE, and top.
00956 { 00957 float fac=0.0 , top ; 00958 00959 if( MRI_IS_INT_TYPE(otype) ){ 00960 top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ; 00961 if (top == 0.0) fac = 0.0; 00962 else fac = MRI_TYPE_maxval[otype]/top; 00963 } 00964 00965 EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ; 00966 return ( fac ); 00967 } |
|
Definition at line 115 of file 3dFDR.c. References PROGRAM_NAME. Referenced by check_one_output_file(), initialize_program(), output_results(), read_options(), and read_time_series().
00116 { 00117 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message); 00118 exit(1); 00119 } |
|
Definition at line 136 of file 3dFDR.c. References MASTER_SHORTHELP_STRING. Referenced by initialize_program().
00137 { 00138 printf( 00139 "This program implements the False Discovery Rate (FDR) algorithm for \n" 00140 "thresholding of voxelwise statistics. \n" 00141 " \n" 00142 "Program input consists of a functional dataset containing one (or more) \n" 00143 "statistical sub-bricks. Output consists of a bucket dataset with one \n" 00144 "sub-brick for each input sub-brick. For non-statistical input sub-bricks, \n" 00145 "the output is a copy of the input. However, statistical input sub-bricks \n" 00146 "are replaced by their corresponding FDR values, as follows: \n" 00147 " \n" 00148 "For each voxel, the minimum value of q is determined such that \n" 00149 " E(FDR) <= q \n" 00150 "leads to rejection of the null hypothesis in that voxel. Only voxels inside\n" 00151 "the user specified mask will be considered. These q-values are then mapped\n" 00152 "to z-scores for compatibility with the AFNI statistical threshold display: \n" 00153 " \n" 00154 " stat ==> p-value ==> FDR q-value ==> FDR z-score \n" 00155 " \n" 00156 ); 00157 00158 printf( 00159 "Usage: \n" 00160 " 3dFDR \n" 00161 " -input fname fname = filename of input 3d functional dataset \n" 00162 " OR \n" 00163 " -input1D dname dname = .1D file containing column of p-values \n" 00164 " \n" 00165 " -mask_file mname Use mask values from file mname. \n" 00166 " Note: If file mname contains more than 1 sub-brick, \n" 00167 " the mask sub-brick must be specified! \n" 00168 " Default: No mask \n" 00169 " \n" 00170 " -mask_thr m Only voxels whose corresponding mask value is \n" 00171 " greater than or equal to m in absolute value will \n" 00172 " be considered. Default: m=1 \n" 00173 " \n" 00174 " Constant c(N) depends on assumption about p-values: \n" 00175 " -cind c(N) = 1 p-values are independent across N voxels \n" 00176 " -cdep c(N) = sum(1/i), i=1,...,N any joint distribution \n" 00177 " Default: c(N) = 1 \n" 00178 " \n" 00179 " -quiet Flag to suppress screen output \n" 00180 " \n" 00181 " -list Write sorted list of voxel q-values to screen \n" 00182 " \n" 00183 " -prefix pname Use 'pname' for the output dataset prefix name. \n" 00184 " OR \n" 00185 " -output pname \n" 00186 " \n" 00187 "\n") ; 00188 00189 printf("\n" MASTER_SHORTHELP_STRING ) ; 00190 00191 exit(0) ; 00192 00193 } |
|
Definition at line 440 of file 3dFDR.c. References addto_args(), AFNI_logger(), argc, check_one_output_file(), commandline, DOPEN, DSET_BRICK_FACTOR, DSET_BRICK_TYPE, DSET_NVALS, DSET_NX, DSET_NY, DSET_NZ, DSET_PRINCIPAL_VALUE, EDIT_coerce_scale_type(), FDR_cn, FDR_error(), FDR_input1D_data, FDR_input1D_filename, FDR_input_filename, FDR_mask, FDR_mask_filename, FDR_mask_thr, FDR_MAX_LL, FDR_nthr, FDR_nxyz, FDR_output_prefix, FDR_quiet, FDR_Syntax(), free, machdep(), malloc, MTEST, nz, PROGRAM_NAME, read_options(), read_time_series(), SUB_POINTER, THD_delete_3dim_dataset(), THD_open_dataset(), and tross_commandline().
00441 { 00442 int iv; /* index number of sub-brick */ 00443 void * vfim = NULL; /* sub-brick data pointer */ 00444 float * ffim = NULL; /* sub-brick data in floating point format */ 00445 int ixyz; /* voxel index */ 00446 int nx, ny, nz, nxyz; /* numbers of voxels in input dataset */ 00447 int mx, my, mz, mxyz; /* numbers of voxels in mask dataset */ 00448 int nthr; /* number of voxels above mask threshold */ 00449 char message[80]; /* error message */ 00450 int ibin; /* p-value bin index */ 00451 00452 00453 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/ 00454 machdep() ; 00455 { int new_argc ; char ** new_argv ; 00456 addto_args( argc , argv , &new_argc , &new_argv ) ; 00457 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; } 00458 } 00459 00460 00461 /*----- Save command line for history notes -----*/ 00462 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ; 00463 00464 00465 /*----- Does user request help menu? -----*/ 00466 if( argc < 2 || strcmp(argv[1],"-help") == 0 ) FDR_Syntax() ; 00467 00468 00469 /*----- Add to program log -----*/ 00470 AFNI_logger (PROGRAM_NAME,argc,argv); 00471 00472 00473 /*----- Read input options -----*/ 00474 read_options( argc , argv ) ; 00475 00476 00477 /*----- Open the mask dataset -----*/ 00478 if (FDR_mask_filename != NULL) 00479 { 00480 if (!FDR_quiet) 00481 printf ("Reading mask dataset: %s \n", FDR_mask_filename); 00482 DOPEN (FDR_dset, FDR_mask_filename); 00483 00484 if (FDR_dset == NULL) 00485 { 00486 sprintf (message, "Cannot open mask dataset %s", FDR_mask_filename); 00487 FDR_error (message); 00488 } 00489 00490 if (DSET_NVALS(FDR_dset) != 1) 00491 FDR_error ("Must specify single sub-brick for mask data"); 00492 00493 00494 /*----- Get dimensions of mask dataset -----*/ 00495 mx = DSET_NX(FDR_dset); 00496 my = DSET_NY(FDR_dset); 00497 mz = DSET_NZ(FDR_dset); 00498 mxyz = mx*my*mz; 00499 00500 00501 /*----- Allocate memory for float data -----*/ 00502 ffim = (float *) malloc (sizeof(float) * mxyz); MTEST (ffim); 00503 00504 00505 /*----- Convert mask dataset sub-brick to floats (in ffim) -----*/ 00506 iv = DSET_PRINCIPAL_VALUE (FDR_dset); 00507 SUB_POINTER (FDR_dset, iv, 0, vfim); 00508 EDIT_coerce_scale_type (mxyz, DSET_BRICK_FACTOR(FDR_dset,iv), 00509 DSET_BRICK_TYPE(FDR_dset,iv), vfim, /* input */ 00510 MRI_float , ffim); /* output */ 00511 00512 00513 /*----- Allocate memory for mask volume -----*/ 00514 FDR_mask = (byte *) malloc (sizeof(byte) * mxyz); 00515 MTEST (FDR_mask); 00516 00517 00518 /*----- Create mask of voxels above mask threshold -----*/ 00519 nthr = 0; 00520 for (ixyz = 0; ixyz < mxyz; ixyz++) 00521 if (fabs(ffim[ixyz]) >= FDR_mask_thr) 00522 { 00523 FDR_mask[ixyz] = 1; 00524 nthr++; 00525 } 00526 else 00527 FDR_mask[ixyz] = 0; 00528 00529 if (!FDR_quiet) 00530 printf ("Number of voxels above mask threshold = %d \n", nthr); 00531 if (nthr < 1) 00532 FDR_error ("No voxels above mask threshold. Cannot continue."); 00533 00534 00535 /*----- Delete floating point sub-brick -----*/ 00536 if (ffim != NULL) { free (ffim); ffim = NULL; } 00537 00538 /*----- Delete mask dataset -----*/ 00539 THD_delete_3dim_dataset (FDR_dset, False); FDR_dset = NULL ; 00540 00541 } 00542 00543 00544 /*----- Get the input data -----*/ 00545 00546 if (FDR_input1D_filename != NULL) 00547 { 00548 /*----- Read the input .1D file -----*/ 00549 if (!FDR_quiet) printf ("Reading input data: %s \n", 00550 FDR_input1D_filename); 00551 FDR_input1D_data = read_time_series (FDR_input1D_filename, &nxyz); 00552 00553 if (FDR_input1D_data == NULL) 00554 { 00555 sprintf (message, "Unable to read input .1D data file: %s", 00556 FDR_input1D_filename); 00557 FDR_error (message); 00558 } 00559 00560 if (nxyz < 1) 00561 { 00562 sprintf (message, "No p-values in input .1D data file: %s", 00563 FDR_input1D_filename); 00564 FDR_error (message); 00565 } 00566 00567 FDR_nxyz = nxyz; 00568 FDR_nthr = nxyz; 00569 } 00570 00571 else 00572 { 00573 /*----- Open the input 3D dataset -----*/ 00574 if (!FDR_quiet) printf ("Reading input dataset: %s \n", 00575 FDR_input_filename); 00576 FDR_dset = THD_open_dataset(FDR_input_filename); 00577 00578 if (FDR_dset == NULL) 00579 { 00580 sprintf (message, "Cannot open input dataset %s", 00581 FDR_input_filename); 00582 FDR_error (message); 00583 } 00584 00585 /*----- Get dimensions of input dataset -----*/ 00586 nx = DSET_NX(FDR_dset); 00587 ny = DSET_NY(FDR_dset); 00588 nz = DSET_NZ(FDR_dset); 00589 nxyz = nx*ny*nz; 00590 00591 00592 /*----- Check for compatible dimensions -----*/ 00593 if (FDR_mask != NULL) 00594 { 00595 if ((nx != mx) || (ny != my) || (nz != mz)) 00596 FDR_error ("Mask and input dataset have incompatible dimensions"); 00597 FDR_nxyz = nxyz; 00598 FDR_nthr = nthr; 00599 } 00600 else 00601 { 00602 FDR_nxyz = nxyz; 00603 FDR_nthr = nxyz; 00604 } 00605 00606 00607 /*----- Check whether output dataset already exists -----*/ 00608 check_one_output_file (FDR_dset, FDR_output_prefix); 00609 } 00610 00611 00612 /*----- Initialize constant c(N) -----*/ 00613 if (FDR_cn < 0.0) 00614 { 00615 double cn; 00616 cn = 0.0; 00617 for (ixyz = 1; ixyz <= FDR_nthr; ixyz++) 00618 cn += 1.0 / ixyz; 00619 FDR_cn = cn; 00620 if (!FDR_quiet) 00621 printf ("c(N) = %f \n", FDR_cn); 00622 } 00623 00624 /*----- Initialize voxel pointers -----*/ 00625 for (ibin = 0; ibin < FDR_MAX_LL; ibin++) 00626 FDR_head_voxel[ibin] = NULL; 00627 00628 00629 } |
|
compute the overall minimum and maximum voxel values for a dataset Definition at line 1111 of file 3dFDR.c. References argc, FDR_input1D_filename, initialize_program(), output_results(), process_1ddata(), process_dataset(), PROGRAM_AUTHOR, PROGRAM_INITIAL, PROGRAM_LATEST, and PROGRAM_NAME.
01113 { 01114 THD_3dim_dataset * new_dset = NULL; /* output bucket dataset */ 01115 01116 01117 /*----- Identify software -----*/ 01118 printf ("\n\n"); 01119 printf ("Program: %s \n", PROGRAM_NAME); 01120 printf ("Author: %s \n", PROGRAM_AUTHOR); 01121 printf ("Initial Release: %s \n", PROGRAM_INITIAL); 01122 printf ("Latest Revision: %s \n", PROGRAM_LATEST); 01123 printf ("\n"); 01124 01125 01126 /*----- Initialize program: get all operator inputs; 01127 create mask for voxels above mask threshold -----*/ 01128 initialize_program (argc, argv); 01129 01130 01131 if (FDR_input1D_filename != NULL) 01132 { 01133 /*----- Process list of p-values -----*/ 01134 process_1ddata (); 01135 } 01136 else 01137 { 01138 /*----- Process 3d dataset -----*/ 01139 new_dset = process_dataset (); 01140 01141 /*----- Output the results as a bucket dataset -----*/ 01142 output_results (new_dset); 01143 } 01144 01145 exit(0) ; 01146 } |
|
Definition at line 689 of file 3dFDR.c. References add_voxel(), create_voxel(), voxel::ixyz, and voxel::pvalue. Referenced by add_voxel().
00691 { 00692 voxel * voxel_ptr = NULL; 00693 00694 voxel_ptr = create_voxel (); 00695 00696 voxel_ptr->ixyz = ixyz; 00697 voxel_ptr->pvalue = pvalue; 00698 00699 head_voxel = add_voxel (voxel_ptr, head_voxel); 00700 00701 return (head_voxel); 00702 00703 } |
|
Definition at line 1082 of file 3dFDR.c. References ADN_func_type, ADN_none, DSET_BRIKNAME, EDIT_dset_items(), FDR_error(), FDR_quiet, FUNC_BUCK_TYPE, THD_delete_3dim_dataset(), THD_load_statistics(), and THD_write_3dim_dataset().
01083 { 01084 int ierror; /* flag for errors in editing dataset */ 01085 01086 01087 /*----- Make sure that output is a bucket dataset -----*/ 01088 ierror = EDIT_dset_items( new_dset , 01089 ADN_func_type , FUNC_BUCK_TYPE, 01090 ADN_none ) ; 01091 if (ierror > 0) 01092 FDR_error ("Errors in attempting to create output dataset."); 01093 01094 01095 /*----- Output the FDR dataset -----*/ 01096 if( !FDR_quiet ) fprintf(stderr,"Computing sub-brick statistics\n") ; 01097 THD_load_statistics( new_dset ) ; 01098 01099 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ; 01100 if( !FDR_quiet ) fprintf(stderr,"Wrote output to %s\n", DSET_BRIKNAME(new_dset) ); 01101 01102 01103 /*----- Deallocate memory for output dataset -----*/ 01104 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ; 01105 01106 } |
|
Definition at line 916 of file 3dFDR.c. References FDR_input1D_data, FDR_nxyz, FDR_output_prefix, free, mri_fix_data_pointer(), mri_free(), mri_new_vol_empty(), mri_write_1D(), and process_volume(). Referenced by main().
00918 { 00919 float * ffim = NULL; /* input list of p-values */ 00920 00921 00922 /*----- Set pointer to input array of p-values -----*/ 00923 ffim = FDR_input1D_data; 00924 00925 00926 /*----- Calculate FDR z-scores for all input p-values -----*/ 00927 process_volume (ffim, -1, NULL); 00928 00929 if( FDR_output_prefix != NULL && ffim != NULL ){ 00930 MRI_IMAGE *im = mri_new_vol_empty( FDR_nxyz,1,1 , MRI_float ) ; 00931 mri_fix_data_pointer( ffim , im ) ; 00932 mri_write_1D( FDR_output_prefix , im ) ; 00933 mri_fix_data_pointer( NULL , im ) ; 00934 mri_free( im ) ; 00935 } 00936 00937 /*----- Deallocate memory -----*/ 00938 if (ffim != NULL) { free (ffim); ffim = NULL; } 00939 00940 } |
|
Definition at line 1031 of file 3dFDR.c. References commandline, DSET_BRICK_STATCODE, DSET_NVALS, EDIT_full_copy(), FDR_output_prefix, FDR_quiet, free, FUNC_IS_STAT, process_subbrick(), THD_delete_3dim_dataset(), tross_Append_History(), and tross_Copy_History(). Referenced by main().
01033 { 01034 THD_3dim_dataset * new_dset = NULL; /* output bucket dataset */ 01035 int ibrick, nbricks; /* sub-brick indices */ 01036 int statcode; /* type of stat. sub-brick */ 01037 01038 01039 /*----- Make full copy of input dataset -----*/ 01040 new_dset = EDIT_full_copy(FDR_dset, FDR_output_prefix); 01041 01042 01043 /*----- Record history of dataset -----*/ 01044 tross_Copy_History( FDR_dset , new_dset ) ; 01045 01046 if( commandline != NULL ) 01047 { 01048 tross_Append_History ( new_dset, commandline); 01049 free(commandline) ; 01050 } 01051 01052 01053 /*----- Deallocate memory for input dataset -----*/ 01054 THD_delete_3dim_dataset (FDR_dset , False ); FDR_dset = NULL ; 01055 01056 01057 /*----- Loop over sub-bricks in the dataset -----*/ 01058 nbricks = DSET_NVALS(new_dset); 01059 for (ibrick = 0; ibrick < nbricks; ibrick++) 01060 { 01061 statcode = DSET_BRICK_STATCODE(new_dset, ibrick); 01062 if (FUNC_IS_STAT(statcode)) 01063 { 01064 /*----- Process the statistical sub-bricks -----*/ 01065 if (!FDR_quiet) 01066 printf ("ibrick = %3d statcode = %5s \n", 01067 ibrick, FUNC_prefixstr[statcode]); 01068 process_subbrick (new_dset, ibrick); 01069 } 01070 } 01071 01072 01073 return (new_dset); 01074 } |
|
Definition at line 975 of file 3dFDR.c. References DSET_BRICK_FACTOR, DSET_BRICK_LABEL, DSET_BRICK_STATAUX, DSET_BRICK_STATCODE, DSET_BRICK_TYPE, EDIT_BRICK_FACTOR, EDIT_BRICK_LABEL, EDIT_BRICK_TO_FIZT, EDIT_coerce_autoscale_new(), EDIT_coerce_scale_type(), FDR_nxyz, FDR_quiet, free, malloc, MTEST, process_volume(), SUB_POINTER, and THD_MAX_NAME. Referenced by process_dataset().
00977 { 00978 const float EPSILON = 1.0e-10; 00979 float factor; /* factor is new scale factor for this sub-brick */ 00980 void * vfim = NULL; /* sub-brick data pointer */ 00981 float * ffim = NULL; /* sub-brick data in floating point format */ 00982 char brick_label[THD_MAX_NAME]; /* sub-brick label */ 00983 00984 00985 if (!FDR_quiet) printf ("Processing sub-brick #%d \n", ibrick); 00986 00987 00988 /*----- Allocate memory for float data -----*/ 00989 ffim = (float *) malloc (sizeof(float) * FDR_nxyz); MTEST (ffim); 00990 00991 00992 /*----- Convert sub-brick to float stats -----*/ 00993 SUB_POINTER (dset, ibrick, 0, vfim); 00994 EDIT_coerce_scale_type (FDR_nxyz, DSET_BRICK_FACTOR(dset,ibrick), 00995 DSET_BRICK_TYPE(dset,ibrick), vfim, /* input */ 00996 MRI_float , ffim); /* output */ 00997 00998 00999 /*----- Calculate FDR z-scores for all voxels within this volume -----*/ 01000 process_volume (ffim, DSET_BRICK_STATCODE(dset,ibrick), 01001 DSET_BRICK_STATAUX (dset,ibrick)); 01002 01003 01004 /*----- Replace old sub-brick with new z-scores -----*/ 01005 SUB_POINTER (dset, ibrick, 0, vfim); 01006 factor = EDIT_coerce_autoscale_new (FDR_nxyz, MRI_float, ffim, 01007 DSET_BRICK_TYPE(dset,ibrick), vfim); 01008 if (factor < EPSILON) factor = 0.0; 01009 else factor = 1.0 / factor; 01010 01011 01012 /*----- edit the sub-brick -----*/ 01013 strcpy (brick_label, "FDRz "); 01014 strcat (brick_label, DSET_BRICK_LABEL(dset, ibrick)); 01015 EDIT_BRICK_LABEL (dset, ibrick, brick_label); 01016 EDIT_BRICK_FACTOR (dset, ibrick, factor); 01017 EDIT_BRICK_TO_FIZT (dset, ibrick); 01018 01019 01020 /*----- Deallocate memory -----*/ 01021 if (ffim != NULL) { free (ffim); ffim = NULL; } 01022 01023 } |
|
Definition at line 769 of file 3dFDR.c. References delete_all_voxels(), FDR_cn, FDR_input1D_filename, FDR_mask, FDR_MAX_LL, FDR_nthr, FDR_nxyz, free, voxel::ixyz, malloc, MTEST, new_voxel(), voxel::next_voxel, normal_p2t(), voxel::pvalue, save_all_voxels(), and THD_stat_to_pval(). Referenced by process_1ddata(), and process_subbrick().
00771 { 00772 int ixyz; /* voxel index */ 00773 int icount; /* count of sorted p-values */ 00774 float fval; /* voxel input statistical value */ 00775 float pval; /* voxel input stat. p-value */ 00776 float qval; /* voxel FDR q-value */ 00777 float zval; /* voxel FDR z-score */ 00778 float qval_min; /* smallest previous q-value */ 00779 voxel * head_voxel = NULL; /* linked list of voxels */ 00780 voxel * voxel_ptr = NULL; /* pointer to current voxel */ 00781 int ibin; /* p-value bin */ 00782 int * iarray = NULL; /* output array of voxel indices */ 00783 float * parray = NULL; /* output array of voxel p-values */ 00784 float * qarray = NULL; /* output array of voxel FDR q-values */ 00785 float * zarray = NULL; /* output array of voxel FDR z-scores */ 00786 00787 00788 /*----- Allocate memory for screen output arrays -----*/ 00789 if (FDR_list) 00790 { 00791 iarray = (int *) malloc (sizeof(int) * FDR_nthr); MTEST(iarray); 00792 parray = (float *) malloc (sizeof(float) * FDR_nthr); MTEST(parray); 00793 qarray = (float *) malloc (sizeof(float) * FDR_nthr); MTEST(qarray); 00794 zarray = (float *) malloc (sizeof(float) * FDR_nthr); MTEST(zarray); 00795 } 00796 00797 00798 /*----- Loop over all voxels; sort p-values -----*/ 00799 icount = FDR_nthr; 00800 for (ixyz = 0; ixyz < FDR_nxyz; ixyz++) 00801 { 00802 00803 /*----- First, check if voxel is inside the mask -----*/ 00804 if (FDR_mask != NULL) 00805 if (!FDR_mask[ixyz]) continue; 00806 00807 00808 /*----- Convert stats to p-values -----*/ 00809 fval = fabs(ffim[ixyz]); 00810 if (statcode <= 0) 00811 pval = fval; 00812 else 00813 pval = THD_stat_to_pval (fval, statcode, stataux); 00814 00815 if (pval >= 1.0) 00816 { 00817 /*----- Count but don't sort voxels with p-value = 1 -----*/ 00818 icount--; 00819 if (FDR_list) 00820 { 00821 iarray[icount] = ixyz; 00822 parray[icount] = 1.0; 00823 qarray[icount] = 1.0; 00824 zarray[icount] = 0.0; 00825 } 00826 } 00827 else 00828 { 00829 /*----- Place voxel in p-value bin -----*/ 00830 ibin = (int) (pval * (FDR_MAX_LL)); 00831 if (ibin < 0) ibin = 0; 00832 if (ibin > FDR_MAX_LL-1) ibin = FDR_MAX_LL-1; 00833 head_voxel = new_voxel (ixyz, pval, FDR_head_voxel[ibin]); 00834 FDR_head_voxel[ibin] = head_voxel; 00835 } 00836 } 00837 00838 00839 /*----- Calculate FDR q-values -----*/ 00840 qval_min = 1.0; 00841 ibin = FDR_MAX_LL-1; 00842 while (ibin >= 0) 00843 { 00844 voxel_ptr = FDR_head_voxel[ibin]; 00845 00846 while (voxel_ptr != NULL) 00847 { 00848 /*----- Convert sorted p-values to FDR q-values -----*/ 00849 pval = voxel_ptr->pvalue; 00850 qval = FDR_cn * (pval*FDR_nthr) / icount; 00851 if (qval > qval_min) 00852 qval = qval_min; 00853 else 00854 qval_min = qval; 00855 00856 /*----- Convert FDR q-value to FDR z-score -----*/ 00857 if (qval < 1.0e-20) 00858 zval = 10.0; 00859 else 00860 zval = normal_p2t(qval); 00861 00862 icount--; 00863 00864 /*----- Save calculated values -----*/ 00865 if (FDR_list) 00866 { 00867 iarray[icount] = voxel_ptr->ixyz; 00868 parray[icount] = pval; 00869 qarray[icount] = qval; 00870 zarray[icount] = zval; 00871 } 00872 00873 voxel_ptr->pvalue = zval; 00874 voxel_ptr = voxel_ptr->next_voxel; 00875 } 00876 00877 ibin--; 00878 } 00879 00880 00881 /*----- Write out the calculated values -----*/ 00882 if (FDR_list) 00883 { 00884 printf ("%12s %12s %12s %12s \n", 00885 "Index", "p-value", "q-value", "z-score"); 00886 for (icount = 0; icount < FDR_nthr; icount++) 00887 { 00888 if (FDR_input1D_filename != NULL) 00889 ixyz = iarray[icount] + 1; 00890 else 00891 ixyz = iarray[icount]; 00892 printf ("%12d %12.6f %12.6f %12.6f \n", 00893 ixyz, parray[icount], qarray[icount], zarray[icount]); 00894 } 00895 00896 /*----- Deallocate memory for output arrays -----*/ 00897 free (iarray); free (parray); free (qarray); free (zarray); 00898 } 00899 00900 00901 /*----- Place FDR z-scores into float array -----*/ 00902 save_all_voxels (ffim); 00903 00904 00905 /*----- Deallocate linked-list memory -----*/ 00906 delete_all_voxels(); 00907 00908 } |
|
Definition at line 202 of file 3dFDR.c. References argc, FDR_cn, FDR_error(), FDR_input1D_filename, FDR_input_filename, FDR_list, FDR_mask_filename, FDR_mask_thr, FDR_output_prefix, FDR_quiet, malloc, MCW_strncpy, MTEST, THD_MAX_NAME, and THD_MAX_PREFIX. Referenced by initialize_program().
00203 { 00204 int nopt = 1 ; /* count of input arguments */ 00205 char message[80]; /* error message */ 00206 00207 00208 00209 /*----- main loop over input options -----*/ 00210 while( nopt < argc ) 00211 { 00212 00213 /*----- -input fname -----*/ 00214 if (strcmp(argv[nopt], "-input") == 0) 00215 { 00216 nopt++; 00217 if (nopt >= argc) FDR_error ("need argument after -input "); 00218 FDR_input_filename = (char *) malloc (sizeof(char)*THD_MAX_NAME); 00219 MTEST (FDR_input_filename); 00220 strcpy (FDR_input_filename, argv[nopt]); 00221 nopt++; 00222 continue; 00223 } 00224 00225 00226 /*----- -input1D dname -----*/ 00227 if (strcmp(argv[nopt], "-input1D") == 0) 00228 { 00229 nopt++; 00230 if (nopt >= argc) FDR_error ("need argument after -input1D "); 00231 FDR_input1D_filename = (char *) malloc (sizeof(char)*THD_MAX_NAME); 00232 MTEST (FDR_input1D_filename); 00233 strcpy (FDR_input1D_filename, argv[nopt]); 00234 nopt++; 00235 continue; 00236 } 00237 00238 00239 /*----- -mask_file mname -----*/ 00240 if (strcmp(argv[nopt], "-mask_file") == 0) 00241 { 00242 nopt++; 00243 if (nopt >= argc) FDR_error ("need argument after -mask_file "); 00244 FDR_mask_filename = (char *) malloc (sizeof(char)*THD_MAX_NAME); 00245 MTEST (FDR_mask_filename); 00246 strcpy (FDR_mask_filename, argv[nopt]); 00247 nopt++; 00248 continue; 00249 } 00250 00251 00252 /*----- -mask_thr m -----*/ 00253 if( strcmp(argv[nopt],"-mask_thr") == 0 ){ 00254 float fval; 00255 nopt++ ; 00256 if( nopt >= argc ){ 00257 FDR_error (" need 1 argument after -mask_thr"); 00258 } 00259 sscanf (argv[nopt], "%f", &fval); 00260 if (fval < 0.0){ 00261 FDR_error (" Require mask_thr >= 0.0 "); 00262 } 00263 FDR_mask_thr = fval; 00264 nopt++; continue; 00265 } 00266 00267 00268 /*----- -cind -----*/ 00269 if( strcmp(argv[nopt],"-cind") == 0 ){ 00270 FDR_cn = 1.0; 00271 nopt++ ; continue ; 00272 } 00273 00274 00275 /*----- -cdep -----*/ 00276 if( strcmp(argv[nopt],"-cdep") == 0 ){ 00277 FDR_cn = -1.0; 00278 nopt++ ; continue ; 00279 } 00280 00281 00282 /*----- -quiet -----*/ 00283 if( strcmp(argv[nopt],"-quiet") == 0 ){ 00284 FDR_quiet = 1; 00285 nopt++ ; continue ; 00286 } 00287 00288 00289 /*----- -list -----*/ 00290 if( strcmp(argv[nopt],"-list") == 0 ){ 00291 FDR_list = 1; 00292 nopt++ ; continue ; 00293 } 00294 00295 00296 /*----- -prefix prefix -----*/ 00297 if( strcmp(argv[nopt],"-prefix") == 0 || 00298 strcmp(argv[nopt],"-output") == 0 ){ 00299 nopt++ ; 00300 if( nopt >= argc ){ 00301 FDR_error (" need argument after -prefix!"); 00302 } 00303 FDR_output_prefix = (char *) malloc (sizeof(char) * THD_MAX_PREFIX); 00304 MCW_strncpy( FDR_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ; 00305 continue ; 00306 } 00307 00308 00309 /*----- unknown command -----*/ 00310 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]); 00311 FDR_error (message); 00312 00313 00314 } /*----- end of loop over command line arguments -----*/ 00315 00316 00317 } |
|
Definition at line 327 of file 3dFDR.c. References far, FDR_error(), malloc, MRI_FLOAT_PTR, mri_free(), mri_read_1D(), MTEST, MRI_IMAGE::nx, MRI_IMAGE::ny, and THD_MAX_NAME.
00332 { 00333 char message[THD_MAX_NAME]; /* error message */ 00334 char * cpt; /* pointer to column suffix */ 00335 char filename[THD_MAX_NAME]; /* time series file name w/o column index */ 00336 char subv[THD_MAX_NAME]; /* string containing column index */ 00337 MRI_IMAGE * im, * flim; /* pointers to image structures 00338 -- used to read 1D ASCII */ 00339 float * far; /* pointer to MRI_IMAGE floating point data */ 00340 int nx; /* number of time points in time series */ 00341 int ny; /* number of columns in time series file */ 00342 int iy; /* time series file column index */ 00343 int ipt; /* time point index */ 00344 float * ts_data = NULL; /* input time series data */ 00345 00346 00347 /*----- First, check for empty filename -----*/ 00348 if (ts_filename == NULL) 00349 FDR_error ("Missing input time series file name"); 00350 00351 00352 /*----- Read the time series file -----*/ 00353 flim = mri_read_1D(ts_filename) ; 00354 if (flim == NULL) 00355 { 00356 sprintf (message, "Unable to read time series file: %s", ts_filename); 00357 FDR_error (message); 00358 } 00359 00360 far = MRI_FLOAT_PTR(flim); 00361 nx = flim->nx; 00362 ny = flim->ny; iy = 0 ; 00363 if( ny > 1 ){ 00364 fprintf(stderr,"WARNING: time series %s has more than 1 column\n",ts_filename); 00365 } 00366 00367 00368 /*----- Save the time series data -----*/ 00369 *ts_length = nx; 00370 ts_data = (float *) malloc (sizeof(float) * nx); 00371 MTEST (ts_data); 00372 for (ipt = 0; ipt < nx; ipt++) 00373 ts_data[ipt] = far[ipt + iy*nx]; 00374 00375 00376 mri_free (flim); flim = NULL; 00377 00378 return (ts_data); 00379 } |
|
Definition at line 738 of file 3dFDR.c. References far, FDR_MAX_LL, FDR_nxyz, voxel::ixyz, voxel::next_voxel, and voxel::pvalue. Referenced by process_volume().
00739 { 00740 int ixyz, ibin; 00741 voxel * voxel_ptr = NULL; /* pointer to voxel */ 00742 00743 00744 /*----- Initialize all voxels to zero -----*/ 00745 for (ixyz = 0; ixyz < FDR_nxyz; ixyz++) 00746 far[ixyz] = 0.0; 00747 00748 00749 for (ibin = 0; ibin < FDR_MAX_LL; ibin++) 00750 { 00751 voxel_ptr = FDR_head_voxel[ibin]; 00752 00753 while (voxel_ptr != NULL) 00754 { 00755 far[voxel_ptr->ixyz] = voxel_ptr->pvalue; 00756 voxel_ptr = voxel_ptr->next_voxel; 00757 } 00758 00759 } 00760 00761 } |
Variable Documentation
|
Definition at line 70 of file 3dFDR.c. Referenced by initialize_program(), and process_dataset(). |
|
Definition at line 57 of file 3dFDR.c. Referenced by initialize_program(), process_volume(), and read_options(). |
|
|
|
|
|
Definition at line 67 of file 3dFDR.c. Referenced by initialize_program(), and process_1ddata(). |
|
Definition at line 62 of file 3dFDR.c. Referenced by initialize_program(), main(), process_volume(), and read_options(). |
|
Definition at line 61 of file 3dFDR.c. Referenced by initialize_program(), and read_options(). |
|
Definition at line 55 of file 3dFDR.c. Referenced by read_options(). |
|
Definition at line 66 of file 3dFDR.c. Referenced by initialize_program(), and process_volume(). |
|
Definition at line 63 of file 3dFDR.c. Referenced by initialize_program(), and read_options(). |
|
Definition at line 56 of file 3dFDR.c. Referenced by initialize_program(), and read_options(). |
|
Definition at line 59 of file 3dFDR.c. Referenced by initialize_program(), and process_volume(). |
|
Definition at line 58 of file 3dFDR.c. Referenced by initialize_program(), process_1ddata(), process_subbrick(), process_volume(), and save_all_voxels(). |
|
Definition at line 64 of file 3dFDR.c. Referenced by initialize_program(), process_1ddata(), process_dataset(), and read_options(). |
|
Definition at line 54 of file 3dFDR.c. Referenced by initialize_program(), output_results(), process_dataset(), process_subbrick(), and read_options(). |