1 /*****************************************************************************
2    Major portions of this software are copyrighted by the Medical College
3    of Wisconsin, 1994-2000, and are released under the Gnu General Public
4    License, Version 2.  See the file README.Copyright for details.
5 ******************************************************************************/
6 
7 /*
8   This program performs the nonparametric Friedman test for randomized
9   complete block design experiments.
10 
11   File:    3dFriedman.c
12   Author:  B. Douglas Ward
13   Date:    23 July 1997
14 
15   Mod:     Added changes for incorporating History notes.
16   Date:    08 September 1999
17 
18   Mod:     Replaced dataset input code with calls to THD_open_dataset,
19            to allow operator selection of individual sub-bricks for input.
20   Date:    03 December 1999
21 
22   Mod:     Added call to AFNI_logger.
23   Date:    14 March 2002
24 
25   Mod:     Modified routine write_afni_fict of NPstats.c so that all output
26            subbricks will now have the scaled short integer format.
27   Date:    14 March 2002
28 
29   Mod:     Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
30   Date:    02 December 2002
31 */
32 
33 
34 /*---------------------------------------------------------------------------*/
35 
36 #define PROGRAM_NAME "3dFriedman"                    /* name of this program */
37 #define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
38 #define PROGRAM_INITIAL "23 July 1997"    /* date of initial program release */
39 #define PROGRAM_LATEST "02 December 2002" /* date of latest program revision */
40 
41 /*---------------------------------------------------------------------------*/
42 
43 #include <stdio.h>
44 #include <math.h>
45 #include "mrilib.h"
46 
47 #define MAX_TREATMENTS 333     /* max. number of treatments */
48 #define MAX_OBSERVATIONS 333   /* max. number of observations per treatment */
49 #define MAX_NAME_LENGTH THD_MAX_NAME   /* max. string length for file names */
50 #define MEGA  1048576          /* one megabyte */
51 
52 #define USE_ARRAY
53 #ifdef  USE_ARRAY
54  static int   ntar = 0 ;     /* 26 Oct 2007 */
55  static float *tar = NULL ;
56 #endif
57 
58 typedef struct NP_options
59 {
60   int   datum;                  /* data type for "intensity" data subbrick */
61   char  session[MAX_NAME_LENGTH];     /* name of output directory */
62 
63 
64   int   nvoxel;                 /* number of voxel for special output */
65 
66   int   s;                      /* number of treatments */
67   int   n[MAX_TREATMENTS];      /* number of observations per treatment */
68 
69   char  *** xname;              /* names of the input data files */
70   char  * first_dataset;        /* name of the first data set */
71 
72   int   nx, ny, nz;             /* data set dimensions */
73   int   nxyz;                   /* number of voxels per image */
74 
75   int workmem;                  /* working memory */
76 
77   char  * outfile;              /* name of output file  */
78 
79 
80 } NP_options;
81 
82 
83 #include "NPstats.c"
84 
85 
86 /*---------------------------------------------------------------------------*/
87 /*
88    Routine to display 3dFriedman help menu.
89 */
90 
display_help_menu()91 void display_help_menu()
92 {
93   printf
94     (
95      "This program performs nonparametric Friedman test for               \n"
96      "randomized complete block design experiments.                     \n\n"
97      "Usage:                                                              \n"
98      "3dFriedman                                                          \n"
99      "-levels s                      s = number of treatments             \n"
100      "-dset 1 filename               data set for treatment #1            \n"
101      " . . .                           . . .                              \n"
102      "-dset 1 filename               data set for treatment #1            \n"
103      " . . .                           . . .                              \n"
104      "-dset s filename               data set for treatment #s            \n"
105      " . . .                           . . .                              \n"
106      "-dset s filename               data set for treatment #s            \n"
107      "                                                                    \n"
108      "[-workmem mega]                number of megabytes of RAM to use    \n"
109      "                                 for statistical workspace          \n"
110      "[-voxel num]                   screen output for voxel # num        \n"
111      "-out prefixname                Friedman statistics are written      \n"
112      "                                 to file prefixname                 \n"
113      "\n");
114 
115   printf
116     (
117      "\n"
118      "N.B.: For this program, the user must specify 1 and only 1 sub-brick  \n"
119      "      with each -dset command. That is, if an input dataset contains  \n"
120      "      more than 1 sub-brick, a sub-brick selector must be used, e.g.: \n"
121      "      -dset 2 'fred+orig[3]'                                          \n"
122      );
123 
124    printf("\n" MASTER_SHORTHELP_STRING ) ;
125 
126   PRINT_COMPILE_DATE ; exit(0);
127 }
128 
129 
130 /*---------------------------------------------------------------------------*/
131 /*
132    Routine to initialize input options.
133 */
134 
initialize_options(NP_options * option_data)135 void initialize_options (NP_options * option_data)
136 {
137   int i;          /* index */
138 
139   option_data->datum = ILLEGAL_TYPE;
140   strcpy (option_data->session, "./");
141 
142 
143   option_data->nvoxel = -1;
144 
145   option_data->s = 0;
146 
147   for (i = 0;  i < MAX_TREATMENTS;  i++)
148     option_data->n[i] = 0;
149 
150   option_data->workmem = 266;
151 
152   /*----- allocate memory for storing data file names -----*/
153   option_data->xname = (char ***) malloc (sizeof(char **) * MAX_TREATMENTS);
154   for (i = 0;  i < MAX_TREATMENTS;  i++)
155     option_data->xname[i]
156       = (char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
157 
158   option_data->first_dataset = NULL;
159 
160   option_data->nx = 0;
161   option_data->ny = 0;
162   option_data->nz = 0;
163   option_data->nxyz = 0;
164 
165   option_data->outfile = NULL;
166 
167 }
168 
169 
170 /*---------------------------------------------------------------------------*/
171 /*
172    Routine to get user specified Friedman options.
173 */
174 
get_options(int argc,char ** argv,NP_options * option_data)175 void get_options (int argc, char ** argv, NP_options * option_data)
176 {
177   int nopt = 1;                  /* input option argument counter */
178   int ival;                      /* integer input */
179   int nijk;                      /* count of data files */
180   float fval;                    /* float input */
181   THD_3dim_dataset * dset=NULL;             /* test whether data set exists */
182   char message[MAX_NAME_LENGTH];            /* error message */
183 
184 
185   /*----- does user request help menu? -----*/
186   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();
187 
188 
189   /*----- add to program log -----*/
190   AFNI_logger (PROGRAM_NAME,argc,argv);
191 
192 
193   /*----- initialize the input options -----*/
194   initialize_options (option_data);
195 
196 
197   /*----- main loop over input options -----*/
198   while (nopt < argc)
199     {
200 
201 
202       /*-----   -datum type   -----*/
203       if( strncmp(argv[nopt],"-datum",6) == 0 ){
204 	if( ++nopt >= argc ) NP_error("need an argument after -datum!") ;
205 
206 	if( strcmp(argv[nopt],"short") == 0 ){
207 	  option_data->datum = MRI_short ;
208 	} else if( strcmp(argv[nopt],"float") == 0 ){
209 	  option_data->datum = MRI_float ;
210 	} else {
211 	  char buf[256] ;
212 	  sprintf(buf,
213 		  "-datum of type '%s' is not supported in 3dFriedman.",
214 		  argv[nopt] ) ;
215 	  NP_error(buf) ;
216 	}
217 	nopt++ ; continue ;  /* go to next arg */
218       }
219 
220 
221       /*-----   -session dirname    -----*/
222       if( strncmp(argv[nopt],"-session",6) == 0 ){
223 	nopt++ ;
224 	if( nopt >= argc ) NP_error("need argument after -session!") ;
225 	strcpy(option_data->session , argv[nopt++]) ;
226 	continue ;
227       }
228 
229 
230       /*-----   -voxel num  -----*/
231       if (strncmp(argv[nopt], "-voxel", 6) == 0)
232 	{
233 	  nopt++;
234 	  if (nopt >= argc)  NP_error ("need argument after -voxel ");
235 	  sscanf (argv[nopt], "%d", &ival);
236 	  if (ival <= 0)
237 	    NP_error ("illegal argument after -voxel ");
238 	  option_data->nvoxel = ival;
239 	  nopt++;
240 	  continue;
241 	}
242 
243 
244       /*-----   -workmem megabytes  -----*/
245 
246       if( strncmp(argv[nopt],"-workmem",6) == 0 ){
247          nopt++ ;
248          if( nopt >= argc ) NP_error ("need argument after -workmem!") ;
249          sscanf (argv[nopt], "%d", &ival);
250          if( ival <= 0 ) NP_error ("illegal argument after -workmem!") ;
251          option_data->workmem = ival ;
252          nopt++ ; continue ;
253       }
254 
255 
256       /*-----   -levels s  -----*/
257       if (strncmp(argv[nopt], "-levels", 7) == 0)
258 	{
259 	  nopt++;
260 	  if (nopt >= argc)  NP_error ("need argument after -levels ");
261 	  sscanf (argv[nopt], "%d", &ival);
262 	  if ((ival <= 0) || (ival > MAX_TREATMENTS))
263 	    NP_error ("illegal argument after -levels ");
264 	  option_data->s = ival;
265 	  nopt++;
266 	  continue;
267 	}
268 
269 
270       /*-----   -dset level filename   -----*/
271       if (strncmp(argv[nopt], "-dset", 5) == 0)
272 	{
273 	  nopt++;
274 	  if (nopt+1 >= argc)  NP_error ("need 2 arguments after -dset ");
275 	  sscanf (argv[nopt], "%d", &ival);
276 	  if ((ival <= 0) || (ival > option_data->s))
277 	    NP_error ("illegal argument after -dset ");
278 
279 	  option_data->n[ival-1] += 1;
280 
281 	  if (option_data->n[ival-1] > MAX_OBSERVATIONS)
282 	    NP_error ("too many data files");
283 	  nijk = option_data->n[ival-1];
284 
285 	  /*--- check whether input files exist ---*/
286 	  nopt++;
287 	  dset = THD_open_dataset( argv[nopt] ) ;
288      CHECK_OPEN_ERROR(dset,argv[nopt]) ;
289 
290 	  /*--- check number of selected sub-bricks ---*/
291 	  if (DSET_NVALS(dset) != 1)
292 	    {
293 	      sprintf(message,"Must specify exactly 1 sub-brick for file %s\n",
294 		      argv[nopt]);
295 	      NP_error (message);
296 	    }
297 
298 	  THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
299 
300 	  option_data->xname[ival-1][nijk-1]
301 	    =  malloc (sizeof(char) * MAX_NAME_LENGTH);
302 	  strcpy (option_data->xname[ival-1][nijk-1], argv[nopt]);
303 	  nopt++;
304 	  continue;
305 	}
306 
307 
308       /*-----   -out filename   -----*/
309       if (strncmp(argv[nopt], "-out", 4) == 0)
310 	{
311 	  nopt++;
312 	  if (nopt >= argc)  NP_error ("need argument after -out ");
313 	  option_data->outfile = malloc (sizeof(char) * MAX_NAME_LENGTH);
314 	  strcpy (option_data->outfile, argv[nopt]);
315 	  nopt++;
316 	  continue;
317 	}
318 
319 
320       /*----- unknown command -----*/
321       NP_error ("unrecognized command line option ");
322     }
323 
324 }
325 
326 
327 /*---------------------------------------------------------------------------*/
328 /*
329   Routine to check for valid inputs.
330 */
331 
check_for_valid_inputs(NP_options * option_data)332 void check_for_valid_inputs (NP_options * option_data)
333 {
334   int i, n;
335   char message[MAX_NAME_LENGTH];            /* error message */
336 
337 
338   n = option_data->n[0];
339   if (n < 1)
340     NP_error ("Sample size is too small");
341 
342   for (i = 1;  i < option_data->s;  i++)
343     if (option_data->n[i] != n)
344       NP_error ("Must have equal sample sizes for all treatments");
345 
346 
347   if (option_data->nvoxel > option_data->nxyz)
348     NP_error ("argument of -voxel is too large");
349 
350 }
351 
352 
353 /*---------------------------------------------------------------------------*/
354 /*
355   Routine to perform all Friedman initialization.
356 */
357 
initialize(int argc,char ** argv,NP_options ** option_data,float ** best,float ** qstat)358 void initialize
359 (
360   int argc,                    /* number of input arguments */
361   char ** argv,                /* array of input arguments */
362   NP_options ** option_data,   /* user input options */
363   float ** best,               /* index of best treatment */
364   float ** qstat               /* Friedman statistic */
365 )
366 
367 {
368 
369 
370   /*----- allocate memory space for input data -----*/
371   *option_data = (NP_options *) malloc(sizeof(NP_options));
372   if (*option_data == NULL)
373     NP_error ("memory allocation error");
374 
375   /*----- get command line inputs -----*/
376   get_options(argc, argv, *option_data);
377 
378   /*----- use first data set to get data set dimensions -----*/
379   (*option_data)->first_dataset = (*option_data)->xname[0][0];
380   get_dimensions (*option_data);
381   printf ("Data set dimensions:  nx = %d  ny = %d  nz = %d  nxyz = %d \n",
382 	  (*option_data)->nx, (*option_data)->ny,
383 	  (*option_data)->nz, (*option_data)->nxyz);
384 
385 
386   /*----- check for valid inputs -----*/
387   check_for_valid_inputs (*option_data);
388 
389   /*----- check whether output files already exist -----*/
390   if( THD_deathcon() ) check_one_output_file (*option_data, (*option_data)->outfile);
391 
392   /*----- allocate memory -----*/
393   *best = (float *) malloc(sizeof(float) * (*option_data)->nxyz);
394   if (*best == NULL)
395     NP_error ("memory allocation error");
396   *qstat = (float *) malloc(sizeof(float) * (*option_data)->nxyz);
397   if (*qstat == NULL)
398     NP_error ("memory allocation error");
399 
400 
401 }
402 
403 
404 /*---------------------------------------------------------------------------*/
405 /*
406   Calculate the Friedman statistic.
407 */
408 
calc_stat(int nvox,int s,int n,float ** xarray,float * best,float * qstat)409 void calc_stat
410 (
411   int nvox,                          /* flag for voxel output */
412   int s,                             /* number of treatments */
413   int n,                             /* number of observations per treatment */
414   float ** xarray,                   /* array of data arrays */
415   float * best,                      /* index of best treatment */
416   float * qstat                      /* Friedman statistic */
417 )
418 
419 {
420   const float EPSILON = 1.0e-10;      /* protection from roundoff error */
421   int i, j;                   /* array indices */
422   node * head = NULL;         /* points to head of list */
423   node * ptr = NULL;          /* points to current position in list */
424   int NN;                     /* total number of sample points */
425   float rsum;                 /* sum of squares of ranks */
426   int d;                      /* count of number of ties */
427   float corr;                 /* correction to account for ties */
428   float rank;                 /* rank of data point */
429   float ranksum;              /* sum of ranks for ith treatment */
430   float qnum;                 /* numerator of Friedman statistic */
431   float qden;                 /* denominator of Friedman statistic */
432   float best_rank;            /* best average rank for a treatment */
433   float ** rank_array;        /* array to store ranks for all observations */
434 
435 #ifdef USE_ARRAY
436   int kk=0 ;
437 #endif
438 
439 
440   /*----- allocate memory for storing ranks -----*/
441   rank_array = (float **) malloc (sizeof(float *) * s);
442   for (i = 0;  i < s;  i++)
443     rank_array[i] = (float *) malloc (sizeof(float) * n);
444 
445 
446   /*----- loop over blocks -----*/
447   corr = 0.0;
448   for (j = 0;  j < n;  j++)
449     {
450 
451       /*----- enter and sort data for each treatment within block j -----*/
452       for (i = 0;  i < s;  i++){
453 #ifdef USE_ARRAY
454         tar[kk++] = xarray[i][j] ;
455 #else
456         node_addvalue (&head, xarray[i][j]);
457 #endif
458       }
459 #ifdef USE_ARRAY
460       node_allatonce( &head , kk , tar ) ;
461 #endif
462 
463 
464       /*----- store the ranks for each treatment within block j -----*/
465       for (i = 0;  i < s;  i++) rank_array[i][j] = node_get_rank (head, xarray[i][j]);
466 
467 
468       /*----- calculate the ties correction factor -----*/
469       ptr = head;
470       while (ptr != NULL)
471 	{
472 	  d = ptr->d;
473 	  corr += d*d*d - d;
474 	  ptr = ptr->next;
475 	}
476 
477       list_delete (&head);
478 #ifdef USE_ARRAY
479       kk = 0 ;
480 #endif
481 
482     } /* j loop */
483 
484 
485   /*----- if display voxel, write the ranks of the input data -----*/
486   if (nvox > 0)
487     {
488       printf ("\n");
489       for (i = 0;  i < s;  i++)
490 	{
491 	  printf ("Y%d ranks: ", i+1);
492 	  for (j = 0;  j < n;  j++)
493 	    {
494 	      rank = rank_array[i][j];
495 	      printf (" %6.1f", rank);
496 	      if (((j+1) % 8 == 0) && (j < n-1))
497 		printf ("\n          ");
498 	    }
499 	  printf ("\n");
500 	  if (n > 8)  printf ("\n");
501 	}
502       printf ("\n");
503       for (i = 0;  i < s;  i++)
504 	{
505 	  printf ("Y%d: ", i+1);
506 	  ranksum = 0.0;
507 	  for (j = 0;  j < n;  j++)
508 	    {
509 	      rank = rank_array[i][j];
510 	      ranksum += rank;
511 	    }
512 	  printf ("   Rank sum = %6.1f    Rank average = %6.1f \n",
513 		  ranksum, ranksum/n);
514 	}
515       printf ("\n");
516     }
517 
518 
519   /*----- calculate the sum of the ranks -----*/
520   rsum = 0.0;
521   *best = 0.0;
522   best_rank = (s + 1.0) / 2.0 + EPSILON;
523   for (i = 0;  i < s;  i++)
524     {
525       ranksum = 0.0;
526       for (j = 0;  j < n;  j++)
527 	ranksum += rank_array[i][j];
528       rsum += ranksum * ranksum;
529 
530       if (ranksum/n > best_rank)
531 	{
532 	  *best = (float) (i+1);
533 	  best_rank = ranksum / n;
534 	}
535     }
536 
537 
538   /*----- numerator of Friedman statistic -----*/
539   qnum = (12.0/(n*s*(s+1)))*rsum - 3.0*n*(s+1);
540 
541   /*----- denominator of Friedman statistic -----*/
542   qden = 1.0 - (corr / (n*s*(s*s-1)));
543 
544   /*----- calculate Friedman statistic -----*/
545   if (qden < EPSILON)
546     *qstat = 0.0;
547   else
548     *qstat = qnum / qden;
549   if (nvox > 0)  printf ("Q = %f \n", *qstat);
550 
551 
552   /*----- deallocate memory -----*/
553   for (i = 0;  i < s;  i++)
554     {
555       free (rank_array[i]);
556       rank_array[i] = NULL;
557     }
558   free (rank_array);
559   rank_array = NULL;
560 
561 }
562 
563 
564 /*---------------------------------------------------------------------------*/
565 /*
566   Calculate results for a single voxel.
567 */
568 
process_voxel(int nvox,int s,int n,float ** xarray,float * best,float * qstat)569 void process_voxel
570 (
571   int nvox,                          /* flag for voxel output */
572   int s,                             /* number of treatments */
573   int n,                             /* number of observations per treatment */
574   float ** xarray,                   /* array of data arrays */
575   float * best,                      /* index of best treatment */
576   float * qstat                      /* Friedman statistic */
577 )
578 
579 {
580   int i;                             /* treatment index */
581   int j;                             /* array index */
582 
583 #ifdef USE_ARRAY
584   if( ntar == 0 ){
585     ntar = s*n+1 ; tar = (float *)malloc(sizeof(float)*ntar) ;
586     if( tar == NULL ) ERROR_exit("Can't malloc 'tar[%d]'!",ntar) ;
587   }
588 #endif
589 
590 
591   /*----- check for voxel output  -----*/
592   if (nvox > 0)
593     {
594       printf ("\nResults for voxel #%d : \n\n", nvox);
595 
596       for (i = 0;  i < s;  i++)
597 	{
598 	  printf ("Y%d data:  ", i+1);
599 	  for (j = 0;  j < n;  j++)
600 	    {
601 	    printf (" %6.1f", xarray[i][j]);
602 	    if (((j+1) % 8 == 0) && (j < n-1))
603 	      printf ("\n          ");
604 	    }
605 	  printf ("\n");
606 	  if (n > 8)  printf ("\n");
607 	}
608       if (n <= 8)  printf ("\n");
609     }
610 
611 
612   /*----- calculate Friedman statistic -----*/
613   calc_stat (nvox, s, n, xarray, best, qstat);
614 
615 }
616 
617 
618 /*---------------------------------------------------------------------------*/
619 /*
620   Calculate the Friedman statistics for all voxels  (by breaking the datasets
621   into sub-volumes, if necessary).
622 */
623 
calculate_results(NP_options * option_data,float * best,float * qstat)624 void calculate_results
625 (
626   NP_options * option_data,    /* user input options */
627   float * best,                /* index of best treatment */
628   float * qstat                /* Friedman statistics */
629 )
630 
631 {
632   int i, j, m;                 /* array indices */
633   int s;                       /* number of treatments */
634   int n;                       /* number of observations per treatment */
635   int nxyz;                    /* number of voxels per dataset */
636   int num_datasets;            /* total number of datasets */
637   int piece_size;              /* number of voxels in dataset sub-volume */
638   int num_pieces;              /* dataset is divided into this many pieces */
639   int piece;                   /* piece index */
640   int piece_len;               /* number of voxels in current sub-volume */
641   int fim_offset;              /* array offset to current sub-volume */
642   int ivox;                    /* index to voxels in current sub-volume */
643   int nvox;                    /* index of voxel within entire volume */
644   float b;                     /* index of best treatment */
645   float q;                     /* Friedman statistic */
646   float ** xfimar;             /* array of sub-volumes of datasets */
647   float ** xarray;             /* array of data arrays */
648 
649 
650   /*----- initialize local variables -----*/
651   s = option_data->s;
652   nxyz = option_data->nxyz;
653   num_datasets = 0;
654   n = option_data->n[0];
655   num_datasets = n * s;
656 
657 
658   /*----- break problem into smaller pieces -----*/
659   piece_size = option_data->workmem * MEGA / (num_datasets * sizeof(float));
660   if (piece_size > nxyz)  piece_size = nxyz;
661   num_pieces = (nxyz + piece_size - 1) / piece_size;
662   printf ("num_pieces = %d    piece_size = %d \n", num_pieces, piece_size);
663 
664 
665   /*----- allocate memory space -----*/
666   xarray = (float **) malloc (sizeof(float *) * s);  MTEST(xarray);
667   for (i = 0;  i < s;  i++)
668     {
669       xarray[i] = (float *) malloc (sizeof(float) * option_data->n[i]);
670       MTEST(xarray[i]);
671     }
672 
673   xfimar = (float **) malloc (sizeof(float *) * num_datasets);  MTEST(xfimar);
674   for (i = 0;  i < num_datasets;  i++)
675     {
676       xfimar[i] = (float *) malloc (sizeof(float) * piece_size);
677       MTEST(xfimar[i]);
678     }
679 
680 
681   /*----- loop over the pieces of the input datasets -----*/
682   nvox = 0;
683   for (piece = 0;  piece < num_pieces;  piece++)
684     {
685       printf ("piece = %d \n", piece);
686       fim_offset = piece * piece_size;
687       if (piece < num_pieces-1)
688 	piece_len = piece_size;
689       else
690 	piece_len = nxyz - fim_offset;
691 
692 
693       /*----- read in sub-volume of data from each dataset -----*/
694       m = 0;
695       for (i = 0;  i < s;  i++)
696 	for (j = 0;  j < option_data->n[i];  j++)
697 	  {
698 	    read_afni_data (option_data, option_data->xname[i][j],
699 			    piece_len, fim_offset, xfimar[m]);
700 	    m++;
701 	  }
702 
703 
704       /*----- loop over voxels in this piece -----*/
705       for (ivox = 0;  ivox < piece_len;  ivox++)
706 	{
707 	  nvox += 1;
708 
709 	  m = 0;
710 	  for (i = 0;  i < s;  i++)
711 	    for (j = 0;  j < option_data->n[i];  j++)
712 	      {
713 		xarray[i][j] = xfimar[m][ivox];
714 		m++;
715 	      }
716 
717 
718 	  /*----- calculate results for this voxel -----*/
719 	  if (nvox == option_data->nvoxel)
720 	    process_voxel (nvox, s, n, xarray, &b, &q);
721 	  else
722 	    process_voxel (-1, s, n, xarray, &b, &q);
723 
724 
725 	  /*----- save results for this voxel -----*/
726 	  best[ivox+fim_offset] = b;
727 	  qstat[ivox+fim_offset] = q;
728 
729 	}
730 
731     }  /* loop over pieces */
732 
733 
734   /*----- deallocate memory -----*/
735   for (i = 0;  i < s;  i++)
736     {
737       free (xarray[i]);   xarray[i] = NULL;
738     }
739   free (xarray);   xarray = NULL;
740 
741   for (i = 0;  i < num_datasets;  i++)
742     {
743       free (xfimar[i]);   xfimar[i] = NULL;
744     }
745   free (xfimar);   xfimar = NULL;
746 }
747 
748 
749 /*---------------------------------------------------------------------------*/
750 /*
751   Generate the requested output.
752 */
753 
output_results(int argc,char ** argv,NP_options * option_data,float * best,float * qstat)754 void output_results
755 (
756   int argc,                         /* number of input arguments */
757   char ** argv,                     /* array of input arguments */
758   NP_options * option_data,         /* user input options */
759   float * best,                     /* index of best treatment */
760   float * qstat                     /* Friedman statistic */
761 )
762 
763 {
764 
765   /*----- write out afni fict data file -----*/
766   write_afni_fict (argc, argv, option_data, option_data->outfile,
767 		   best, qstat, option_data->s - 1);
768 
769 }
770 
771 
772 
773 /*---------------------------------------------------------------------------*/
774 /*
775    Routine to release memory and remove any remaining temporary data files.
776 */
777 
terminate(NP_options ** option_data,float ** best,float ** qstat)778 void terminate
779 (
780   NP_options ** option_data,   /* user input options */
781   float ** best,               /* index of best treatment */
782   float ** qstat               /* Friedman statistics */
783 )
784 
785 {
786    int i, j;                       /* dataset indices */
787 
788 
789    /*----- deallocate memory -----*/
790    for (i = 0;  i < (*option_data)->s;  i++)
791      for (j = 0;  j < (*option_data)->n[i];  j++)
792        {
793 	 free ((*option_data)->xname[i][j]);
794 	 (*option_data)->xname[i][j] = NULL;
795        }
796    for (i = 0;  i < (*option_data)->s;  i++)
797      {
798        free ((*option_data)->xname[i]);
799        (*option_data)->xname[i] = NULL;
800      }
801    free ((*option_data)->xname);
802    (*option_data)->xname = NULL;
803 
804    if ((*option_data)->outfile != NULL)
805    {
806       free ((*option_data)-> outfile);
807       (*option_data)->outfile = NULL;
808    }
809 
810    free (*option_data);   *option_data = NULL;
811 
812    free (*best);          *best = NULL;
813 
814    free (*qstat);         *qstat = NULL;
815 }
816 
817 
818 /*---------------------------------------------------------------------------*/
819 /*
820    Perform nonparametric Friedman test for comparison of multiple
821    treatments.
822 */
823 
main(int argc,char ** argv)824 int main (int argc, char ** argv)
825 {
826   NP_options * option_data = NULL;    /* user input options */
827   float * best;                       /* index of best treatment */
828   float * qstat;                      /* Friedman statistic */
829 
830 
831   /*----- Identify software -----*/
832 #if 0
833   printf ("\n\n");
834   printf ("Program: %s \n", PROGRAM_NAME);
835   printf ("Author:  %s \n", PROGRAM_AUTHOR);
836   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
837   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
838   printf ("\n");
839 #endif
840 
841    PRINT_VERSION("3dFriedman") ; AUTHOR(PROGRAM_AUTHOR); mainENTRY("3dFriedman main") ; machdep() ;
842 
843    /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
844 
845    { int new_argc ; char ** new_argv ;
846      addto_args( argc , argv , &new_argc , &new_argv ) ;
847      if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
848    }
849 
850 
851   /*----- program initialization -----*/
852   initialize (argc, argv, &option_data, &best, &qstat);
853 
854   /*----- calculate nonparameteric Friedman statistics -----*/
855   calculate_results (option_data, best, qstat);
856 
857   /*----- generate requested output -----*/
858   output_results (argc, argv, option_data, best, qstat);
859 
860   /*----- terminate program -----*/
861   terminate (&option_data, &best, &qstat);
862 
863 }
864 
865 
866 
867 
868 
869 
870 
871 
872 
873 
874 
875 
876 
877 
878 
879 
880 
881