00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #define PROGRAM_NAME "3dFriedman"
00037 #define PROGRAM_AUTHOR "B. Douglas Ward"
00038 #define PROGRAM_INITIAL "23 July 1997"
00039 #define PROGRAM_LATEST "02 December 2002"
00040
00041
00042
00043 #include <stdio.h>
00044 #include <math.h>
00045 #include "mrilib.h"
00046
00047 #define MAX_TREATMENTS 100
00048 #define MAX_OBSERVATIONS 100
00049 #define MAX_NAME_LENGTH THD_MAX_NAME
00050 #define MEGA 1048576
00051
00052
00053 typedef struct NP_options
00054 {
00055 int datum;
00056 char session[MAX_NAME_LENGTH];
00057
00058
00059 int nvoxel;
00060
00061 int s;
00062 int n[MAX_TREATMENTS];
00063
00064 char *** xname;
00065 char * first_dataset;
00066
00067 int nx, ny, nz;
00068 int nxyz;
00069
00070 int workmem;
00071
00072 char * outfile;
00073
00074
00075 } NP_options;
00076
00077
00078 #include "NPstats.c"
00079
00080
00081
00082
00083
00084
00085
00086 void display_help_menu()
00087 {
00088 printf
00089 (
00090 "This program performs nonparametric Friedman test for \n"
00091 "randomized complete block design experiments. \n\n"
00092 "Usage: \n"
00093 "3dFriedman \n"
00094 "-levels s s = number of treatments \n"
00095 "-dset 1 filename data set for treatment #1 \n"
00096 " . . . . . . \n"
00097 "-dset 1 filename data set for treatment #1 \n"
00098 " . . . . . . \n"
00099 "-dset s filename data set for treatment #s \n"
00100 " . . . . . . \n"
00101 "-dset s filename data set for treatment #s \n"
00102 " \n"
00103 "[-workmem mega] number of megabytes of RAM to use \n"
00104 " for statistical workspace \n"
00105 "[-voxel num] screen output for voxel # num \n"
00106 "-out prefixname Friedman statistics are written \n"
00107 " to file prefixname \n"
00108 "\n");
00109
00110 printf
00111 (
00112 "\n"
00113 "N.B.: For this program, the user must specify 1 and only 1 sub-brick \n"
00114 " with each -dset command. That is, if an input dataset contains \n"
00115 " more than 1 sub-brick, a sub-brick selector must be used, e.g.: \n"
00116 " -dset 2 'fred+orig[3]' \n"
00117 );
00118
00119 printf("\n" MASTER_SHORTHELP_STRING ) ;
00120
00121 exit(0);
00122 }
00123
00124
00125
00126
00127
00128
00129
00130 void initialize_options (NP_options * option_data)
00131 {
00132 int i;
00133
00134 option_data->datum = ILLEGAL_TYPE;
00135 strcpy (option_data->session, "./");
00136
00137
00138 option_data->nvoxel = -1;
00139
00140 option_data->s = 0;
00141
00142 for (i = 0; i < MAX_TREATMENTS; i++)
00143 option_data->n[i] = 0;
00144
00145 option_data->workmem = 12;
00146
00147
00148 option_data->xname = (char ***) malloc (sizeof(char **) * MAX_TREATMENTS);
00149 for (i = 0; i < MAX_TREATMENTS; i++)
00150 option_data->xname[i]
00151 = (char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
00152
00153 option_data->first_dataset = NULL;
00154
00155 option_data->nx = 0;
00156 option_data->ny = 0;
00157 option_data->nz = 0;
00158 option_data->nxyz = 0;
00159
00160 option_data->outfile = NULL;
00161
00162 }
00163
00164
00165
00166
00167
00168
00169
00170 void get_options (int argc, char ** argv, NP_options * option_data)
00171 {
00172 int nopt = 1;
00173 int ival;
00174 int nijk;
00175 float fval;
00176 THD_3dim_dataset * dset=NULL;
00177 char message[MAX_NAME_LENGTH];
00178
00179
00180
00181 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
00182
00183
00184
00185 AFNI_logger (PROGRAM_NAME,argc,argv);
00186
00187
00188
00189 initialize_options (option_data);
00190
00191
00192
00193 while (nopt < argc)
00194 {
00195
00196
00197
00198 if( strncmp(argv[nopt],"-datum",6) == 0 ){
00199 if( ++nopt >= argc ) NP_error("need an argument after -datum!") ;
00200
00201 if( strcmp(argv[nopt],"short") == 0 ){
00202 option_data->datum = MRI_short ;
00203 } else if( strcmp(argv[nopt],"float") == 0 ){
00204 option_data->datum = MRI_float ;
00205 } else {
00206 char buf[256] ;
00207 sprintf(buf,
00208 "-datum of type '%s' is not supported in 3dFriedman.",
00209 argv[nopt] ) ;
00210 NP_error(buf) ;
00211 }
00212 nopt++ ; continue ;
00213 }
00214
00215
00216
00217 if( strncmp(argv[nopt],"-session",6) == 0 ){
00218 nopt++ ;
00219 if( nopt >= argc ) NP_error("need argument after -session!") ;
00220 strcpy(option_data->session , argv[nopt++]) ;
00221 continue ;
00222 }
00223
00224
00225
00226 if (strncmp(argv[nopt], "-voxel", 6) == 0)
00227 {
00228 nopt++;
00229 if (nopt >= argc) NP_error ("need argument after -voxel ");
00230 sscanf (argv[nopt], "%d", &ival);
00231 if (ival <= 0)
00232 NP_error ("illegal argument after -voxel ");
00233 option_data->nvoxel = ival;
00234 nopt++;
00235 continue;
00236 }
00237
00238
00239
00240
00241 if( strncmp(argv[nopt],"-workmem",6) == 0 ){
00242 nopt++ ;
00243 if( nopt >= argc ) NP_error ("need argument after -workmem!") ;
00244 sscanf (argv[nopt], "%d", &ival);
00245 if( ival <= 0 ) NP_error ("illegal argument after -workmem!") ;
00246 option_data->workmem = ival ;
00247 nopt++ ; continue ;
00248 }
00249
00250
00251
00252 if (strncmp(argv[nopt], "-levels", 7) == 0)
00253 {
00254 nopt++;
00255 if (nopt >= argc) NP_error ("need argument after -levels ");
00256 sscanf (argv[nopt], "%d", &ival);
00257 if ((ival <= 0) || (ival > MAX_TREATMENTS))
00258 NP_error ("illegal argument after -levels ");
00259 option_data->s = ival;
00260 nopt++;
00261 continue;
00262 }
00263
00264
00265
00266 if (strncmp(argv[nopt], "-dset", 5) == 0)
00267 {
00268 nopt++;
00269 if (nopt+1 >= argc) NP_error ("need 2 arguments after -dset ");
00270 sscanf (argv[nopt], "%d", &ival);
00271 if ((ival <= 0) || (ival > option_data->s))
00272 NP_error ("illegal argument after -dset ");
00273
00274 option_data->n[ival-1] += 1;
00275
00276 if (option_data->n[ival-1] > MAX_OBSERVATIONS)
00277 NP_error ("too many data files");
00278 nijk = option_data->n[ival-1];
00279
00280
00281 nopt++;
00282 dset = THD_open_dataset( argv[nopt] ) ;
00283 if( ! ISVALID_3DIM_DATASET(dset) )
00284 {
00285 sprintf(message,"Unable to open dataset file %s\n", argv[nopt]);
00286 NP_error (message);
00287 }
00288
00289
00290 if (DSET_NVALS(dset) != 1)
00291 {
00292 sprintf(message,"Must specify exactly 1 sub-brick for file %s\n",
00293 argv[nopt]);
00294 NP_error (message);
00295 }
00296
00297 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00298
00299 option_data->xname[ival-1][nijk-1]
00300 = malloc (sizeof(char) * MAX_NAME_LENGTH);
00301 strcpy (option_data->xname[ival-1][nijk-1], argv[nopt]);
00302 nopt++;
00303 continue;
00304 }
00305
00306
00307
00308 if (strncmp(argv[nopt], "-out", 4) == 0)
00309 {
00310 nopt++;
00311 if (nopt >= argc) NP_error ("need argument after -out ");
00312 option_data->outfile = malloc (sizeof(char) * MAX_NAME_LENGTH);
00313 strcpy (option_data->outfile, argv[nopt]);
00314 nopt++;
00315 continue;
00316 }
00317
00318
00319
00320 NP_error ("unrecognized command line option ");
00321 }
00322
00323 }
00324
00325
00326
00327
00328
00329
00330
00331 void check_for_valid_inputs (NP_options * option_data)
00332 {
00333 int i, n;
00334 char message[MAX_NAME_LENGTH];
00335
00336
00337 n = option_data->n[0];
00338 if (n < 1)
00339 NP_error ("Sample size is too small");
00340
00341 for (i = 1; i < option_data->s; i++)
00342 if (option_data->n[i] != n)
00343 NP_error ("Must have equal sample sizes for all treatments");
00344
00345
00346 if (option_data->nvoxel > option_data->nxyz)
00347 NP_error ("argument of -voxel is too large");
00348
00349 }
00350
00351
00352
00353
00354
00355
00356
00357 void initialize
00358 (
00359 int argc,
00360 char ** argv,
00361 NP_options ** option_data,
00362 float ** best,
00363 float ** qstat
00364 )
00365
00366 {
00367
00368
00369
00370 *option_data = (NP_options *) malloc(sizeof(NP_options));
00371 if (*option_data == NULL)
00372 NP_error ("memory allocation error");
00373
00374
00375 get_options(argc, argv, *option_data);
00376
00377
00378 (*option_data)->first_dataset = (*option_data)->xname[0][0];
00379 get_dimensions (*option_data);
00380 printf ("Data set dimensions: nx = %d ny = %d nz = %d nxyz = %d \n",
00381 (*option_data)->nx, (*option_data)->ny,
00382 (*option_data)->nz, (*option_data)->nxyz);
00383
00384
00385
00386 check_for_valid_inputs (*option_data);
00387
00388
00389 check_one_output_file (*option_data, (*option_data)->outfile);
00390
00391
00392 *best = (float *) malloc(sizeof(float) * (*option_data)->nxyz);
00393 if (*best == NULL)
00394 NP_error ("memory allocation error");
00395 *qstat = (float *) malloc(sizeof(float) * (*option_data)->nxyz);
00396 if (*qstat == NULL)
00397 NP_error ("memory allocation error");
00398
00399
00400 }
00401
00402
00403
00404
00405
00406
00407
00408 void calc_stat
00409 (
00410 int nvox,
00411 int s,
00412 int n,
00413 float ** xarray,
00414 float * best,
00415 float * qstat
00416 )
00417
00418 {
00419 const float EPSILON = 1.0e-10;
00420 int i, j;
00421 node * head = NULL;
00422 node * ptr = NULL;
00423 int NN;
00424 float rsum;
00425 int d;
00426 float corr;
00427 float rank;
00428 float ranksum;
00429 float qnum;
00430 float qden;
00431 float best_rank;
00432 float ** rank_array;
00433
00434
00435
00436
00437 rank_array = (float **) malloc (sizeof(float *) * s);
00438 for (i = 0; i < s; i++)
00439 rank_array[i] = (float *) malloc (sizeof(float) * n);
00440
00441
00442
00443 corr = 0.0;
00444 for (j = 0; j < n; j++)
00445 {
00446
00447
00448 for (i = 0; i < s; i++)
00449 node_addvalue (&head, xarray[i][j]);
00450
00451
00452
00453 for (i = 0; i < s; i++)
00454 rank_array[i][j] = node_get_rank (head, xarray[i][j]);
00455
00456
00457
00458 ptr = head;
00459 while (ptr != NULL)
00460 {
00461 d = ptr->d;
00462 corr += d*d*d - d;
00463 ptr = ptr->next;
00464 }
00465
00466 list_delete (&head);
00467
00468 }
00469
00470
00471
00472 if (nvox > 0)
00473 {
00474 printf ("\n");
00475 for (i = 0; i < s; i++)
00476 {
00477 printf ("Y%d ranks: ", i+1);
00478 for (j = 0; j < n; j++)
00479 {
00480 rank = rank_array[i][j];
00481 printf (" %6.1f", rank);
00482 if (((j+1) % 8 == 0) && (j < n-1))
00483 printf ("\n ");
00484 }
00485 printf ("\n");
00486 if (n > 8) printf ("\n");
00487 }
00488 printf ("\n");
00489 for (i = 0; i < s; i++)
00490 {
00491 printf ("Y%d: ", i+1);
00492 ranksum = 0.0;
00493 for (j = 0; j < n; j++)
00494 {
00495 rank = rank_array[i][j];
00496 ranksum += rank;
00497 }
00498 printf (" Rank sum = %6.1f Rank average = %6.1f \n",
00499 ranksum, ranksum/n);
00500 }
00501 printf ("\n");
00502 }
00503
00504
00505
00506 rsum = 0.0;
00507 *best = 0.0;
00508 best_rank = (s + 1.0) / 2.0 + EPSILON;
00509 for (i = 0; i < s; i++)
00510 {
00511 ranksum = 0.0;
00512 for (j = 0; j < n; j++)
00513 ranksum += rank_array[i][j];
00514 rsum += ranksum * ranksum;
00515
00516 if (ranksum/n > best_rank)
00517 {
00518 *best = (float) (i+1);
00519 best_rank = ranksum / n;
00520 }
00521 }
00522
00523
00524
00525 qnum = (12.0/(n*s*(s+1)))*rsum - 3.0*n*(s+1);
00526
00527
00528 qden = 1.0 - (corr / (n*s*(s*s-1)));
00529
00530
00531 if (qden < EPSILON)
00532 *qstat = 0.0;
00533 else
00534 *qstat = qnum / qden;
00535 if (nvox > 0) printf ("Q = %f \n", *qstat);
00536
00537
00538
00539 for (i = 0; i < s; i++)
00540 {
00541 free (rank_array[i]);
00542 rank_array[i] = NULL;
00543 }
00544 free (rank_array);
00545 rank_array = NULL;
00546
00547 }
00548
00549
00550
00551
00552
00553
00554
00555 void process_voxel
00556 (
00557 int nvox,
00558 int s,
00559 int n,
00560 float ** xarray,
00561 float * best,
00562 float * qstat
00563 )
00564
00565 {
00566 int i;
00567 int j;
00568
00569
00570
00571 if (nvox > 0)
00572 {
00573 printf ("\nResults for voxel #%d : \n\n", nvox);
00574
00575 for (i = 0; i < s; i++)
00576 {
00577 printf ("Y%d data: ", i+1);
00578 for (j = 0; j < n; j++)
00579 {
00580 printf (" %6.1f", xarray[i][j]);
00581 if (((j+1) % 8 == 0) && (j < n-1))
00582 printf ("\n ");
00583 }
00584 printf ("\n");
00585 if (n > 8) printf ("\n");
00586 }
00587 if (n <= 8) printf ("\n");
00588 }
00589
00590
00591
00592 calc_stat (nvox, s, n, xarray, best, qstat);
00593
00594 }
00595
00596
00597
00598
00599
00600
00601
00602
00603 void calculate_results
00604 (
00605 NP_options * option_data,
00606 float * best,
00607 float * qstat
00608 )
00609
00610 {
00611 int i, j, m;
00612 int s;
00613 int n;
00614 int nxyz;
00615 int num_datasets;
00616 int piece_size;
00617 int num_pieces;
00618 int piece;
00619 int piece_len;
00620 int fim_offset;
00621 int ivox;
00622 int nvox;
00623 float b;
00624 float q;
00625 float ** xfimar;
00626 float ** xarray;
00627
00628
00629
00630 s = option_data->s;
00631 nxyz = option_data->nxyz;
00632 num_datasets = 0;
00633 n = option_data->n[0];
00634 num_datasets = n * s;
00635
00636
00637
00638 piece_size = option_data->workmem * MEGA / (num_datasets * sizeof(float));
00639 if (piece_size > nxyz) piece_size = nxyz;
00640 num_pieces = (nxyz + piece_size - 1) / piece_size;
00641 printf ("num_pieces = %d piece_size = %d \n", num_pieces, piece_size);
00642
00643
00644
00645 xarray = (float **) malloc (sizeof(float *) * s); MTEST(xarray);
00646 for (i = 0; i < s; i++)
00647 {
00648 xarray[i] = (float *) malloc (sizeof(float) * option_data->n[i]);
00649 MTEST(xarray[i]);
00650 }
00651
00652 xfimar = (float **) malloc (sizeof(float *) * num_datasets); MTEST(xfimar);
00653 for (i = 0; i < num_datasets; i++)
00654 {
00655 xfimar[i] = (float *) malloc (sizeof(float) * piece_size);
00656 MTEST(xfimar[i]);
00657 }
00658
00659
00660
00661 nvox = 0;
00662 for (piece = 0; piece < num_pieces; piece++)
00663 {
00664 printf ("piece = %d \n", piece);
00665 fim_offset = piece * piece_size;
00666 if (piece < num_pieces-1)
00667 piece_len = piece_size;
00668 else
00669 piece_len = nxyz - fim_offset;
00670
00671
00672
00673 m = 0;
00674 for (i = 0; i < s; i++)
00675 for (j = 0; j < option_data->n[i]; j++)
00676 {
00677 read_afni_data (option_data, option_data->xname[i][j],
00678 piece_len, fim_offset, xfimar[m]);
00679 m++;
00680 }
00681
00682
00683
00684 for (ivox = 0; ivox < piece_len; ivox++)
00685 {
00686 nvox += 1;
00687
00688 m = 0;
00689 for (i = 0; i < s; i++)
00690 for (j = 0; j < option_data->n[i]; j++)
00691 {
00692 xarray[i][j] = xfimar[m][ivox];
00693 m++;
00694 }
00695
00696
00697
00698 if (nvox == option_data->nvoxel)
00699 process_voxel (nvox, s, n, xarray, &b, &q);
00700 else
00701 process_voxel (-1, s, n, xarray, &b, &q);
00702
00703
00704
00705 best[ivox+fim_offset] = b;
00706 qstat[ivox+fim_offset] = q;
00707
00708 }
00709
00710 }
00711
00712
00713
00714 for (i = 0; i < s; i++)
00715 {
00716 free (xarray[i]); xarray[i] = NULL;
00717 }
00718 free (xarray); xarray = NULL;
00719
00720 for (i = 0; i < num_datasets; i++)
00721 {
00722 free (xfimar[i]); xfimar[i] = NULL;
00723 }
00724 free (xfimar); xfimar = NULL;
00725 }
00726
00727
00728
00729
00730
00731
00732
00733 void output_results
00734 (
00735 int argc,
00736 char ** argv,
00737 NP_options * option_data,
00738 float * best,
00739 float * qstat
00740 )
00741
00742 {
00743
00744
00745 write_afni_fict (argc, argv, option_data, option_data->outfile,
00746 best, qstat, option_data->s - 1);
00747
00748 }
00749
00750
00751
00752
00753
00754
00755
00756
00757 void terminate
00758 (
00759 NP_options ** option_data,
00760 float ** best,
00761 float ** qstat
00762 )
00763
00764 {
00765 int i, j;
00766
00767
00768
00769 for (i = 0; i < (*option_data)->s; i++)
00770 for (j = 0; j < (*option_data)->n[i]; j++)
00771 {
00772 free ((*option_data)->xname[i][j]);
00773 (*option_data)->xname[i][j] = NULL;
00774 }
00775 for (i = 0; i < (*option_data)->s; i++)
00776 {
00777 free ((*option_data)->xname[i]);
00778 (*option_data)->xname[i] = NULL;
00779 }
00780 free ((*option_data)->xname);
00781 (*option_data)->xname = NULL;
00782
00783 if ((*option_data)->outfile != NULL)
00784 {
00785 free ((*option_data)-> outfile);
00786 (*option_data)->outfile = NULL;
00787 }
00788
00789 free (*option_data); *option_data = NULL;
00790
00791 free (*best); *best = NULL;
00792
00793 free (*qstat); *qstat = NULL;
00794 }
00795
00796
00797
00798
00799
00800
00801
00802
00803 int main (int argc, char ** argv)
00804 {
00805 NP_options * option_data = NULL;
00806 float * best;
00807 float * qstat;
00808
00809
00810
00811 printf ("\n\n");
00812 printf ("Program: %s \n", PROGRAM_NAME);
00813 printf ("Author: %s \n", PROGRAM_AUTHOR);
00814 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
00815 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
00816 printf ("\n");
00817
00818
00819
00820
00821 machdep() ;
00822 { int new_argc ; char ** new_argv ;
00823 addto_args( argc , argv , &new_argc , &new_argv ) ;
00824 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00825 }
00826
00827
00828
00829 initialize (argc, argv, &option_data, &best, &qstat);
00830
00831
00832 calculate_results (option_data, best, qstat);
00833
00834
00835 output_results (argc, argv, option_data, best, qstat);
00836
00837
00838 terminate (&option_data, &best, &qstat);
00839
00840 }
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857