1 /*********************** 3dBrickStat.c *************************************/
2 /* Author: Daniel Glen, 26 Apr 2005 */
3 #include "mrilib.h"
4 #include "thd_shear3d.h"
5
6 static int datum = MRI_float ;
7 static void Print_Header_MinMax(int Minflag, int Maxflag,
8 THD_3dim_dataset * dset);
9 static void Max_func(int Minflag, int Maxflag, int Meanflag, int Countflag,
10 int Posflag, int Negflag, int Zeroflag, int Absflag,
11 int nan_flag, int Sumflag,
12 int Varflag, int Volflag, THD_3dim_dataset * dset,
13 byte *mmm, int mmvox);
14
usage_3dBrickStat(int detail)15 void usage_3dBrickStat(int detail) {
16 printf(
17 "Usage: 3dBrickStat [options] dataset\n"
18 "Compute maximum and/or minimum voxel values of an input dataset\n"
19 "\n"
20 "The output is a number to the console. The input dataset\n"
21 "may use a sub-brick selection list, as in program 3dcalc.\n"
22 "\n"
23 "Note that this program computes ONE number as the output; e.g.,\n"
24 "the mean over all voxels and time points. If you want (say) the\n"
25 "mean over all voxels but for each time point individually, see\n"
26 "program 3dmaskave.\n"
27 "\n"
28 "Note: If you don't specify one sub-brick, the parameter you get\n"
29 "----- back is computed from all the sub-bricks in dataset.\n"
30 "Options :\n"
31 " -quick = get the information from the header only (default)\n"
32 " -slow = read the whole dataset to find the min and max values\n"
33 " all other options except min and max imply slow\n"
34 " -min = print the minimum value in dataset\n"
35 " -max = print the maximum value in dataset (default)\n"
36 " -mean = print the mean value in dataset \n"
37 " -sum = print the sum of values in the dataset\n"
38 " -var = print the variance in the dataset \n"
39 " -stdev = print the standard deviation in the dataset \n"
40 " -stdev and -var are mutually exclusive\n"
41 " -count = print the number of voxels included\n"
42 " -volume = print the volume of voxels included in microliters\n"
43 " -positive = include only positive voxel values \n"
44 " -negative = include only negative voxel values \n"
45 " -zero = include only zero voxel values \n"
46 " -non-positive = include only voxel values 0 or negative \n"
47 " -non-negative = include only voxel values 0 or greater \n"
48 " -non-zero = include only voxel values not equal to 0 \n"
49 " -absolute = use absolute value of voxel values for all calculations\n"
50 " can be combined with restrictive non-positive, non-negative,\n"
51 " etc. even if not practical. Ignored for percentile and\n"
52 " median computations.\n"
53 " -nan = include only voxel values that are finite numbers, \n"
54 " not NaN, or inf. -nan forces -slow mode.\n"
55 " -nonan =exclude voxel values that are not numbers\n"
56 " -mask dset = use dset as mask to include/exclude voxels\n"
57 " -mrange MIN MAX = Only accept values between MIN and MAX (inclusive)\n"
58 " from the mask. Default it to accept all non-zero\n"
59 " voxels.\n"
60 " -mvalue VAL = Only accept values equal to VAL from the mask.\n"
61 " -automask = automatically compute mask for dataset\n"
62 " Can not be combined with -mask\n"
63 " -percentile p0 ps p1 write the percentile values starting\n"
64 " at p0%% and ending at p1%% at a step of ps%%\n"
65 " Output is of the form p%% value p%% value ...\n"
66 " Percentile values are output first. \n"
67 " Only one sub-brick is accepted as input with this option.\n"
68 " Write the author if you REALLY need this option\n"
69 " to work with multiple sub-bricks.\n"
70 " -perclist NUM_PERC PERC1 PERC2 ...\n"
71 " Like -percentile, but output the given percentiles, rather\n"
72 " than a list on an evenly spaced grid using 'ps'.\n"
73 " -median a shortcut for -percentile 50 1 50 (or -perclist 1 50)\n"
74 " -perc_quiet = only print percentile results, not input percentile cutoffs\n"
75 " -ver = print author and version info\n"
76 " -help = print this help screen\n"
77 ) ;
78 printf("\n" MASTER_SHORTHELP_STRING ) ;
79 PRINT_COMPILE_DATE ;
80 return;
81 }
82 /*
83 static void Max_tsfunc( double tzero , double tdelta ,
84 int npts , float ts[] , double ts_mean ,
85 double ts_slope , void *ud , int nbriks, float *val ) ;
86 static float minvalue=1E10, maxvalue=-1E10;
87 */
88
89 /*! compute the overall minimum and maximum voxel values for a dataset */
main(int argc,char * argv[])90 int main( int argc , char * argv[] )
91 {
92 THD_3dim_dataset * old_dset , * new_dset ; /* input and output datasets */
93 int nopt, nbriks;
94 int slow_flag, quick_flag, min_flag, max_flag, mean_flag,
95 automask,count_flag, sum_flag, var_flag, absolute_flag;
96 int positive_flag, negative_flag, zero_flag, nan_flag, perc_flag, vol_flag;
97
98 byte * mmm=NULL ;
99 int mmvox=0 ;
100 int nxyz, i;
101 float *dvec = NULL, mmin=0.0, mmax=0.0;
102 int N_mp=0, perc_quiet=0;
103 double *mpv=NULL, *perc = NULL;
104 double mp =0.0f, mp0 = 0.0f, mps = 0.0f, mp1 = 0.0f, di =0.0f ;
105 byte *mmf = NULL;
106 MRI_IMAGE *anat_im = NULL;
107 char *mask_dset_name=NULL;
108 void *tmp_vec = NULL;
109
110 /*----- Read command line -----*/
111
112 mainENTRY("3dBrickStat main"); machdep();
113 AFNI_logger("3dBrickStat",argc,argv);
114 nopt = 1 ;
115
116 min_flag = 0;
117 max_flag = -1;
118 mean_flag = 0;
119 sum_flag = 0;
120 var_flag = 0;
121 slow_flag = 0;
122 quick_flag = -1;
123 automask = 0;
124 count_flag = 0;
125 vol_flag = 0;
126 positive_flag = -1;
127 negative_flag = -1;
128 absolute_flag = 0;
129 zero_flag = -1;
130 nan_flag = -1;
131 perc_flag = 0;
132 mmin = 1.0;
133 mmax = -1.0;
134 mask_dset_name = NULL;
135
136 datum = MRI_float;
137 while( nopt < argc && argv[nopt][0] == '-' ){
138 if( strcmp(argv[nopt],"-help") == 0 ||
139 strcmp(argv[nopt],"-h") == 0){
140 usage_3dBrickStat(strlen(argv[nopt])> 3 ? 2:1);
141 exit(0);
142 }
143
144 if( strcmp(argv[nopt],"-ver") == 0 ){
145 PRINT_VERSION("3dBrickStat"); AUTHOR("Daniel Glen");
146 nopt++; continue;
147 }
148
149 if( strcmp(argv[nopt],"-quick") == 0 ){
150 quick_flag = 1;
151 nopt++; continue;
152 }
153
154 if( strcmp(argv[nopt],"-percentile") == 0 ){
155 perc_flag = 1;
156 ++nopt;
157 if (nopt + 2 >= argc) {
158 ERROR_exit( "** Error: Need 3 parameter after -percentile\n");
159 }
160 mp0 = atof(argv[nopt])/100.0f; ++nopt;
161 mps = atof(argv[nopt])/100.0f; ++nopt;
162 mp1 = atof(argv[nopt])/100.0f;
163 if (mps == 0.0f) {
164 ERROR_exit( "** Error: step cannot be 0" );
165 }
166 if (mp0 < 0 || mp0 > 100 || mp1 < 0 || mp1 > 100) {
167 ERROR_exit( "** Error: p0 and p1 must be >=0 and <= 100" );
168 }
169
170 nopt++; continue;
171 }
172
173 if( strcmp(argv[nopt],"-perclist") == 0 ) {
174 /* initialize N_mp and mpv from user list */
175 perc_flag = 1;
176 ++nopt;
177
178 /* init N_mp */
179 if (nopt + 1 >= argc)
180 ERROR_exit( "Need at least 2 params after -perclist\n");
181
182 N_mp = atoi(argv[nopt]);
183 if( N_mp <= 0 )
184 ERROR_exit("need NUM_PERC (and percs) after -perclist\n"
185 " (have NUM_PERC = '%s')\n", argv[nopt]);
186 ++nopt;
187
188 /* allocate mpv */
189 mpv = (double *)malloc(sizeof(double)*N_mp);
190 if (!mpv)
191 ERROR_exit("Failed to allocate %d doubles for mpv", N_mp);
192 if (nopt + N_mp >= argc)
193 ERROR_exit("Need %d percentiles for -perclist\n", N_mp);
194
195 /* and populate mpv */
196 for(i=0; i<N_mp; i++) {
197 /* use mp0 just for a local var */
198 mp0 = atof(argv[nopt])/100.0f;
199 if (mp0 < 0 || mp0 > 100
200 || (mp0 == 0 && ! isdigit(argv[nopt][0])) ) {
201 ERROR_message("** Error: bad -perclist perc #%d of %d: %s\n"
202 " percentiles should be in the range [0,100]\n",
203 i+1, N_mp, argv[nopt]);
204 exit(1);
205 }
206
207 mpv[i] = mp0;
208 ++nopt;
209 }
210
211 continue;
212 }
213
214 if( strcmp(argv[nopt],"-perc_quiet") == 0 ){
215 perc_quiet = 1;
216 nopt++; continue;
217 }
218
219 if( strcmp(argv[nopt],"-median") == 0 ){
220 perc_flag = 1;
221 mp0 = 0.50f;
222 mps = 0.01f;
223 mp1 = 0.50f;
224 nopt++; continue;
225 }
226
227 if( strcmp(argv[nopt],"-slow") == 0 ){
228 slow_flag = 1;
229 nopt++; continue;
230 }
231
232 if( strcmp(argv[nopt],"-min") == 0 ){
233 min_flag = 1;
234 nopt++; continue;
235 }
236
237 if( strcmp(argv[nopt],"-max") == 0 ){
238 max_flag = 1;
239 nopt++; continue;
240 }
241
242 if( strcmp(argv[nopt],"-sum") == 0 ){
243 sum_flag = 1;
244 nopt++; continue;
245 }
246
247 if( strcmp(argv[nopt],"-mean") == 0 ){
248 mean_flag = 1;
249 nopt++; continue;
250 }
251
252 if( strcmp(argv[nopt],"-var") == 0 ){
253 if (var_flag) {
254 ERROR_message("Looks like -stdev is already used.\n"
255 "-var and -stdev are mutually exclusive");
256 exit (1);
257 }
258 var_flag = 1;
259 nopt++; continue;
260 }
261
262 if( strcmp(argv[nopt],"-stdev") == 0 ){
263 if (var_flag) {
264 ERROR_message("Looks like -var is already used.\n"
265 "-var and -stdev are mutually exclusive");
266 exit (1);
267 }
268 var_flag = 2;
269 nopt++; continue;
270 }
271
272 if( strcmp(argv[nopt],"-count") == 0 ){
273 count_flag = 1;
274 nopt++; continue;
275 }
276
277 if( strcmp(argv[nopt],"-volume") == 0 ){
278 vol_flag = 1;
279 nopt++; continue;
280 }
281
282 if( strcmp(argv[nopt],"-positive") == 0 ){
283 if(positive_flag!=-1) {
284 ERROR_exit( "Can not use multiple +/-/0 options");
285
286 }
287 positive_flag = 1;
288 negative_flag = 0;
289 zero_flag = 0;
290 nopt++; continue;
291 }
292
293 if( strcmp(argv[nopt],"-negative") == 0 ){
294 if(positive_flag!=-1) {
295 ERROR_exit( "Can not use multiple +/-/0 options");
296 }
297 positive_flag = 0;
298 negative_flag = 1;
299 zero_flag = 0;
300 nopt++; continue;
301 }
302
303 if( strcmp(argv[nopt],"-zero") == 0 ){
304 if(positive_flag!=-1) {
305 ERROR_exit( "Can not use multiple +/-/0 options");
306 }
307 positive_flag = 0;
308 negative_flag = 0;
309 zero_flag = 1;
310 nopt++; continue;
311 }
312
313 if( strcmp(argv[nopt],"-non-positive") == 0 ){
314 if(positive_flag!=-1) {
315 ERROR_exit( "Can not use multiple +/-/0 options");
316 }
317 positive_flag = 0;
318 negative_flag = 1;
319 zero_flag = 1;
320 nopt++; continue;
321 }
322 if( strcmp(argv[nopt],"-non-negative") == 0 ){
323 if(positive_flag!=-1) {
324 ERROR_exit( "Can not use multiple +/-/0 options");
325 }
326 positive_flag = 1;
327 negative_flag = 0;
328 zero_flag = 1;
329 nopt++; continue;
330 }
331
332 if( strcmp(argv[nopt],"-non-zero") == 0 ){
333 if(positive_flag!=-1) {
334 ERROR_exit( "Can not use multiple +/-/0 options");
335 }
336 positive_flag = 1;
337 negative_flag = 1;
338 zero_flag = 0;
339 nopt++; continue;
340 }
341
342 if( strcmp(argv[nopt],"-absolute") == 0 ){
343 absolute_flag = 1;
344 nopt++; continue;
345 }
346
347 if( strcmp(argv[nopt],"-nan") == 0 ){
348 if(nan_flag!=-1) {
349 ERROR_exit( "Can not use both -nan -nonan options");
350 }
351 nan_flag = 1;
352 nopt++; continue;
353 }
354
355 if( strcmp(argv[nopt],"-nonan") == 0 ){
356 if(nan_flag!=-1) {
357 ERROR_exit( "Can not use both -nan -nonan options");
358
359 }
360 nan_flag = 0;
361 nopt++; continue;
362 }
363
364 if( strcmp(argv[nopt],"-autoclip") == 0 ||
365 strcmp(argv[nopt],"-automask") == 0 ){
366
367 if( mmm != NULL ){
368 ERROR_exit(" ERROR: can't use -autoclip/mask with -mask!");
369
370 }
371 automask = 1 ; nopt++ ; continue ;
372 }
373
374 if( strcmp(argv[nopt],"-mrange") == 0 ){
375 if (nopt+2 >= argc) {
376 ERROR_exit(" ERROR: Need two values after -mrange");
377 }
378 mmin = atof(argv[++nopt]);
379 mmax = atof(argv[++nopt]);
380 if (mmax < mmin) {
381 ERROR_exit(
382 "1st value in -mrange %s %s should be the smallest one",
383 argv[nopt-1], argv[nopt]);
384 }
385 nopt++ ; continue ;
386 }
387
388 if( strcmp(argv[nopt],"-mvalue") == 0 ){
389 if (nopt+1 >= argc) {
390 ERROR_exit(" ERROR: Need 1 value after -mvalue");
391 }
392 mmin = atof(argv[++nopt]);
393 mmax = mmin ;
394 nopt++ ; continue ;
395 }
396
397 if( strcmp(argv[nopt],"-mask") == 0 ){
398 if( mask_dset_name != NULL )
399 ERROR_exit(" ERROR: can't have 2 -mask options!");
400 mask_dset_name = argv[++nopt];
401 nopt++ ; continue ;
402 }
403
404 ERROR_message( " Error - unknown option %s", argv[nopt]);
405 suggest_best_prog_option(argv[0], argv[nopt]);
406 exit(1);
407 }
408
409 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
410 ERROR_message("Too few options");
411 usage_3dBrickStat(0);
412 exit(1) ;
413 }
414
415 if (mask_dset_name) {
416 int ninmask = 0;
417 THD_3dim_dataset * mask_dset ;
418 if( automask ){
419 ERROR_exit(" ERROR: can't use -mask with -automask!");
420 }
421 mask_dset = THD_open_dataset(mask_dset_name) ;
422 CHECK_OPEN_ERROR(mask_dset,mask_dset_name) ;
423
424 mmm = THD_makemask( mask_dset , 0 , mmin, mmax ) ;
425 mmvox = DSET_NVOX( mask_dset ) ;
426 ninmask = THD_countmask (mmvox, mmm);
427 if (!ninmask) {
428 ERROR_exit(" No voxels in mask !");
429 }
430 /* text output program, so avoid extras 26 Dec 2013 [rickr] */
431 /* INFO_message("%d voxels in mask\n", ninmask); */
432 DSET_delete(mask_dset) ;
433 }
434
435 if(((mmm!=NULL) && (quick_flag))||(automask &&quick_flag)) {
436 if(quick_flag==1)
437 WARNING_message( "+++ Warning - can't have quick option with mask");
438 quick_flag = 0;
439 slow_flag = 1;
440 }
441
442 /* if max_flag is not set by user, check if other user options set */
443 if(max_flag==-1) {
444 if(min_flag || mean_flag || count_flag || vol_flag || sum_flag
445 || perc_flag || var_flag)
446 max_flag = 0;
447 else
448 max_flag = 1; /* otherwise check only for max */
449 }
450
451 if((var_flag==1)||(mean_flag==1)||(count_flag==1)||
452 (vol_flag==1)||(absolute_flag==1) ||
453 (positive_flag!=-1)||(nan_flag!=-1)||
454 (sum_flag == 1)||(perc_flag == 1) || (var_flag==2))
455 {
456 /* mean flag or count_flag implies slow */
457 slow_flag = 1;
458 }
459
460 /* check slow and quick options */
461 if((slow_flag) && (quick_flag!=1)) /* if user asked for slow, do so */
462 quick_flag = 0;
463 else
464 quick_flag = 1;
465
466 if((max_flag==0) && (min_flag==0)) /* if the user only asked for mean */
467 quick_flag = 0; /* no need to do quick way */
468
469 if((quick_flag) &&
470 ((absolute_flag==1)||(positive_flag==1)||(negative_flag==1)||(zero_flag==1)))
471 WARNING_message( " Warning - ignoring +/-/0/abs flags for quick computations");
472
473 if(positive_flag==-1) { /* if no +/-/0 options set, allow all voxels */
474 positive_flag = 1;
475 negative_flag = 1;
476 zero_flag = 1;
477 }
478
479 /*----- read input dataset -----*/
480
481 if( nopt >= argc ){
482 ERROR_exit(" No input dataset!?");
483 }
484
485 old_dset = THD_open_dataset( argv[nopt] ) ;
486 CHECK_OPEN_ERROR(old_dset,argv[nopt]) ;
487
488 nxyz = DSET_NVOX(old_dset) ;
489 if( mmm != NULL && mmvox != nxyz ){
490 ERROR_exit(" Mask and input datasets not the same size!") ;
491
492 }
493
494 if(automask && mmm == NULL ){
495 mmm = THD_automask( old_dset ) ;
496 for(i=0;i<nxyz;i++) {
497 if(mmm[i]!=0) ++mmvox;
498 }
499 }
500
501 if(quick_flag)
502 Print_Header_MinMax(min_flag, max_flag, old_dset);
503
504 if(slow_flag!=1)
505 exit(0);
506
507 /* ZSS do some diddlyiddly sorting - DO not affect Daniel's
508 function later on */
509 if (perc_flag == 1) {
510 DSET_mallocize (old_dset);
511 DSET_load (old_dset);
512 if (DSET_NVALS(old_dset) != 1) {
513 ERROR_exit( "-percentile can only be used on one sub-brick only.\n"
514 "Use sub-brick selectors '[.]' to specify "
515 "sub-brick of interest.\n");
516 }
517
518 /* prep for input and output of percentiles */
519 /* if N_mp > 0, it and mpv have already been initialized 15 Mar 2021 */
520 /* (just allocate perc) */
521 if( N_mp > 0 ) {
522 perc = (double *)malloc(sizeof(double)*N_mp);
523 if (!mpv || !perc) {
524 ERROR_message("Failed to allocate mpv or perc");
525 exit(1);
526 }
527 } else {
528 if (mp0 > mp1) {
529 N_mp = 1;
530 } else {
531 /* allocate one above ceiling to prevent truncation error
532 (and crash), N_mp is recomputed anyway 16 Mar 2009 [rickr] */
533 N_mp = (int)((double)(mp1-mp0)/(double)mps) + 2;
534 }
535 mpv = (double *)malloc(sizeof(double)*N_mp);
536 perc = (double *)malloc(sizeof(double)*N_mp);
537 if (!mpv || !perc) {
538 ERROR_message("Failed to allocate for mpv or perc");
539 exit(1);
540 }
541 N_mp = 0;
542 mp = mp0;
543 do {
544 mpv[N_mp] = mp; ++N_mp; mp += mps;
545 } while (mp <= mp1+.00000001);
546 }
547
548 // [PT: March 24, 2021] Squash a bug that affects mean/stdev
549 // calcs when a non-full-FOV mask is used, by setting the arg
550 // 'option' to be 1, not 0 in Percentate(). This way, the input
551 // dset is not sorted/changed (while the mask wasn't). Below,
552 // non-percentile calcs should now be OK.
553 tmp_vec = Percentate (DSET_ARRAY(old_dset, 0), mmm, nxyz,
554 DSET_BRICK_TYPE(old_dset, 0), mpv, N_mp,
555 1, perc,
556 zero_flag, positive_flag, negative_flag );
557
558 if ( !tmp_vec ) {
559 ERROR_message("Failed to compute percentiles.");
560 exit(1);
561 }
562
563 /* take care of brick factor */
564 if (DSET_BRICK_FACTOR(old_dset,0)) {
565 for (i=0; i<N_mp; ++i) {
566 perc[i] = perc[i]*DSET_BRICK_FACTOR(old_dset,0);
567 }
568 }
569
570 for (i=0; i<N_mp; ++i) {
571 if( perc_quiet )
572 fprintf(stdout,"%f ", perc[i]);
573 else
574 fprintf(stdout,"%.1f %f ", mpv[i]*100.0f, perc[i]);
575 }
576
577 free(mpv); mpv = NULL;
578 free(perc); perc = NULL;
579 free(tmp_vec); tmp_vec = NULL;
580 }
581
582 Max_func(min_flag, max_flag, mean_flag,count_flag,
583 positive_flag, negative_flag, zero_flag, absolute_flag,
584 nan_flag, sum_flag, var_flag, vol_flag,old_dset, mmm, mmvox);
585
586
587 if(mmm!=NULL)
588 free(mmm);
589
590 exit(0);
591
592 /* unused code time series method for extracting data */
593 #if 0
594 EDIT_dset_items( old_dset ,
595 ADN_ntt , DSET_NVALS(old_dset) ,
596 ADN_ttorg , 0.0 ,
597 ADN_ttdel , 1.0 ,
598 ADN_tunits , UNITS_SEC_TYPE ,
599 NULL ) ;
600 nbriks = 1;
601
602 /*------------- ready to compute new min, max -----------*/
603 new_dset = MAKER_4D_to_typed_fbuc(
604 old_dset , /* input dataset */
605 "temp" , /* output prefix */
606 datum , /* output datum */
607 0 , /* ignore count */
608 0 , /* can't detrend in maker function KRH 12/02*/
609 nbriks , /* number of briks */
610 Max_tsfunc , /* timeseries processor */
611 NULL, /* data for tsfunc */
612 NULL, /* mask */
613 0 /* Allow auto scaling of output */
614 ) ;
615 if(min_flag)
616 printf("%-13.6g ", minvalue);
617 if(max_flag)
618 printf("%-13.6g", maxvalue);
619 printf("\n");
620 exit(0) ;
621 #endif
622 }
623
624 /*! Print the minimum and maximum values from the header */
625 static void
Print_Header_MinMax(Minflag,Maxflag,dset)626 Print_Header_MinMax(Minflag, Maxflag, dset)
627 int Minflag, Maxflag;
628 THD_3dim_dataset * dset;
629 {
630 int ival, nval_per;
631 float tf=0.0;
632 double scaledmin, scaledmax, internalmin;
633 double internalmax, overallmin, overallmax;
634
635 overallmin = 1E10;
636 overallmax = -1E10;
637
638 ENTRY("Print_Header_MinMax");
639 /* print out stuff for each sub-brick */
640 nval_per = dset->dblk->nvals ;
641 for( ival=0 ; ival < nval_per ; ival++ ){
642 tf = DSET_BRICK_FACTOR(dset,ival) ;
643 if( ISVALID_STATISTIC(dset->stats) ){
644 if( tf != 0.0 ){
645 internalmin = dset->stats->bstat[ival].min/tf;
646 internalmax = dset->stats->bstat[ival].max/tf;
647 }
648 scaledmin = dset->stats->bstat[ival].min;
649 scaledmax = dset->stats->bstat[ival].max;
650 if( tf != 0.0 ){
651 if(internalmin < overallmin)
652 overallmin = scaledmin;
653 if(internalmax > overallmax)
654 overallmax = scaledmax;
655 }
656 else {
657 if(scaledmin < overallmin)
658 overallmin = scaledmin;
659 if(scaledmax > overallmax)
660 overallmax = scaledmax;
661 }
662 }
663 else {
664 WARNING_message("No valid statistics in header. \n"
665 "Use -slow option to generate a new one.") ;
666 EXRETURN;
667 }
668 }
669
670 if(Minflag)
671 printf("%-13.6g ", overallmin);
672 if(Maxflag)
673 printf("%-13.6g", overallmax);
674 if( tf != 0.0)
675 printf(" [*%g]\n",tf) ;
676 else
677 printf("\n");
678
679 EXRETURN;
680 }
681
682
683 /*! search whole dataset for minimum and maximum */
684 /* load all at one time */
Max_func(int Minflag,int Maxflag,int Meanflag,int Countflag,int Posflag,int Negflag,int Zeroflag,int Absflag,int nan_flag,int Sumflag,int Varflag,int Volflag,THD_3dim_dataset * dset,byte * mmm,int mmvox)685 static void Max_func(int Minflag, int Maxflag, int Meanflag, int Countflag,
686 int Posflag, int Negflag, int Zeroflag, int Absflag,
687 int nan_flag,
688 int Sumflag,
689 int Varflag, int Volflag, THD_3dim_dataset * dset,
690 byte *mmm, int mmvox)
691 {
692 double overallmin, overallmax, overallmean;
693 double voxval, fac, sum, sum2, vr;
694 int nvox, npts;
695 int i,k;
696 int test_flag;
697
698 MRI_IMAGE *data_im = NULL;
699
700 ENTRY("Max_func");
701
702 /* maybe the mask came up empty [11 Jun 2019 rickr] */
703 if( mmvox > 0 ) {
704 overallmin = 1E10;
705 overallmax = -1E10;
706 } else {
707 overallmin = 0;
708 overallmax = 0;
709 }
710 sum = 0.0;
711 vr = 0.0;
712 sum2 = 0.0;
713 npts = 0;
714 DSET_mallocize (dset);
715 DSET_load (dset); /* load dataset */
716 npts = 0; /* keep track of number of points */
717 for(i=0;i<dset->dblk->nvals; i++) { /* for each sub-brik in dataset */
718 data_im = DSET_BRICK (dset, i); // set pointer to 0th sub-brik of dset
719 fac = DSET_BRICK_FACTOR(dset, i); /* get scale factor for each sub-brik*/
720 if(fac==0.0) fac=1.0;
721 if( mmm != NULL) /* masked out */
722 nvox = mmvox;
723 else
724 nvox = data_im->nvox; /* number of voxels in the sub-brik */
725
726 for(k=0;k<nvox;k++) {
727 if( mmm != NULL && mmm[k] == 0 ) continue ; /* masked out */
728
729 switch( data_im->kind ){
730 case MRI_short:{
731 short *ar = mri_data_pointer(data_im) ;
732 voxval = ar[k];
733 }
734 break ;
735
736 case MRI_byte:{
737 byte *ar = mri_data_pointer(data_im) ;
738 voxval = ar[k];
739 }
740 break ;
741
742 case MRI_float:{
743 float *ar = mri_data_pointer(data_im) ;
744 voxval = ar[k];
745 }
746 break ;
747
748 case MRI_double:{
749 double *ar = mri_data_pointer(data_im) ;
750 voxval = ar[k];
751 }
752 break ;
753
754 case MRI_int:{
755 int *ar = mri_data_pointer(data_im) ;
756 voxval = ar[k];
757 }
758 break ;
759
760 case MRI_complex:{
761 complex *ar = mri_data_pointer(data_im) ;
762 voxval = CABS(ar[k]);
763 }
764 break ;
765
766 case MRI_rgb:{
767 byte *ar = mri_data_pointer(data_im) ;
768 voxval = 0.299*ar[3*k]+0.587*ar[3*k+1]+0.114*ar[3*k+2];
769 }
770 break ;
771
772 default: /* unknown type */
773 voxval = 0.0; /* ignore this voxel */
774 k = nvox; /* skip to next sub-brik */
775 WARNING_message("Unknown type, %s, in sub-brik %d",
776 MRI_TYPE_name[data_im->kind], i);
777 break;
778 }
779
780 if( mmm == NULL || ((mmm!=NULL) && mmm[k] != 0 )){ // masked in voxel?
781 voxval = voxval * fac; /* apply scale factor */
782 if(nan_flag!=-1) { // check for various not a numbers
783 test_flag = isfinite(voxval);
784 if((nan_flag==1) && (test_flag==1)) /* only looking for NaNs*/
785 continue;
786 if((nan_flag==0) && (test_flag==0)) // only looking for finites
787 continue;
788 if(test_flag==0) { /* not a number */
789 ++npts;
790 continue;
791 }
792 }
793 /* use only absolute values */
794 if(Absflag) voxval = abs(voxval);
795 /* limit data by sign */
796 if(((voxval<0)&&Negflag)||((voxval==0)&&Zeroflag)
797 ||((voxval>0)&&Posflag)) {
798 sum += voxval;
799 if (Varflag) sum2 += voxval*voxval;
800 ++npts;
801 if(voxval<overallmin)
802 overallmin = voxval;
803 if(voxval>overallmax)
804 overallmax = voxval;
805 }
806 }
807 }
808 }
809 if(Minflag)
810 printf("%-13.6g ", overallmin);
811 if(Maxflag)
812 printf("%-13.6g ", overallmax);
813 if(Meanflag)
814 {
815 /* 11 Jun 2019 */
816 if( npts > 0 ) overallmean = sum/npts;
817 else overallmean = 0.0;
818 printf("%-13.6g ", overallmean);
819 }
820 if(Countflag)
821 printf("%-13d", npts);
822
823 if(Volflag)
824 printf("%-13.6f", DSET_VOXVOL(dset) * npts);
825
826
827 if (Sumflag)
828 printf("%-13.6g ", sum);
829 if (Varflag) {
830 /* 11 Jun 2019 */
831 if( npts > 1 ) vr = (sum2-sum*sum/(double)npts)/(double)(npts-1);
832 else vr = 0.0;
833
834 if (Varflag == 2) printf("%-13.6g ", sqrt(vr));
835 else printf("%-13.6g ", vr);
836 }
837 printf("\n");
838
839 mri_free (data_im);
840 /* DSET_unload_one (dset, 0);*/
841 EXRETURN;
842 }
843
844 /* unused code time series method for extracting data */
845 #if 0
846
847 /*! search whole dataset for minimum and maximum */
848 static void Max_tsfunc( double tzero, double tdelta ,
849 int npts, float ts[],
850 double ts_mean, double ts_slope,
851 void * ud, int nbriks, float * val )
852 {
853 static int nvox, ncall;
854 int i;
855
856 ENTRY("Max_tsfunc");
857 /* ts is input vector data */
858
859 /** is this a "notification"? **/
860 if( val == NULL ){
861
862 if( npts > 0 ){ /* the "start notification" */
863
864 nvox = npts ; /* keep track of */
865 ncall = 0 ; /* number of calls */
866
867 } else { /* the "end notification" */
868
869 /* nothing to do here */
870 }
871 return ;
872 }
873 ncall++;
874 for(i=0;i<npts;i++) {
875 if(ts[i]>maxvalue)
876 maxvalue = ts[i];
877 if(ts[i]<minvalue)
878 minvalue = ts[i];
879 }
880
881 EXRETURN;
882 }
883 #endif
884