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