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 *
9 * plug_maskcalc.c - plugin to do mask-based computations
10 *
11 * $Log$
12 * Revision 1.7 2005/04/19 21:07:17 rwcox
13 * Cput
14 *
15 * Revision 1.6 2004/01/07 19:50:37 rwcox
16 * Cput
17 *
18 * Revision 1.5 2003/07/15 13:28:30 rwcox
19 * Cput
20 *
21 * Revision 1.4 2003/06/25 20:45:07 rwcox
22 * Cput
23 *
24 * Revision 1.3 2000/12/21 16:10:54 cox
25 * AFNI
26 *
27 * Revision 1.1 1999/08/06 19:10:48 cox
28 * AFNI
29 *
30 * Revision 1.1 1998/11/12 18:11:06 rickr
31 * Initial revision
32 *
33 ***********************************************************************
34 */
35 #include <stdio.h>
36 #include <math.h>
37 #include <stdlib.h>
38 #if !(defined(DARWIN) || defined(FreeBSD))
39 # include <malloc.h>
40 #endif
41 #include <string.h>
42 #include <sys/types.h>
43 #include <sys/stat.h>
44
45 #include "mrilib.h"
46 #include "afni.h"
47 #include "plug_maskcalc.h"
48
49 #ifndef ALLOW_PLUGINS
50 # error "Plugins not properly set up -- see machdep.h"
51 #endif
52
53 char * MASKCALC_main( PLUGIN_interface * );
54
55 static char grMessage[ R_MESSAGE_L ];
56
57 static char * gr_help_message =
58 "maskcalc plugin - rickr";
59
60
61 DEFINE_PLUGIN_PROTOTYPE
62
PLUGIN_init(int ncall)63 PLUGIN_interface * PLUGIN_init( int ncall )
64 {
65 PLUGIN_interface * plint;
66
67 if ( ncall > 0 ) return NULL; /* only one interface */
68
69 CHECK_IF_ALLOWED("MASKCALC","maskcalc") ; /* 30 Sep 2016 */
70
71 plint = PLUTO_new_interface( "maskcalc", "masked computations on datasets",
72 gr_help_message, PLUGIN_CALL_VIA_MENU, MASKCALC_main );
73
74 PLUTO_add_hint( plint, "Wouldn't some cookies be right tasty?" );
75
76 PLUTO_set_sequence( plint , "z:Reynolds" ) ;
77
78 /* first input : the operation */
79
80 PLUTO_add_option( plint, "Function", "op_st", TRUE );
81 PLUTO_add_hint ( plint, "function to perform on the data" );
82 PLUTO_add_string( plint, "Operation", gr_num_ops, gr_op_strings, 0 );
83 PLUTO_add_hint ( plint, "function to perform on the data" );
84
85 /* second input : the mask */
86
87 PLUTO_add_option ( plint, "Dataset", "mask_st", TRUE );
88 PLUTO_add_hint ( plint, "dataset to be used as mask" );
89 PLUTO_add_dataset( plint, "Mask", ANAT_ALL_MASK , FUNC_ALL_MASK,
90 DIMEN_3D_MASK | DIMEN_4D_MASK | BRICK_SHORT_MASK );
91 PLUTO_add_hint ( plint, "dataset to be used as mask" );
92
93 /* third input : the computational dataset */
94
95 PLUTO_add_option ( plint, "Dataset", "dset_st", TRUE );
96 PLUTO_add_hint ( plint, "computational dataset" );
97 PLUTO_add_dataset( plint, "Dset", ANAT_ALL_MASK, FUNC_ALL_MASK,
98 DIMEN_ALL_MASK | BRICK_SHORT_MASK );
99 PLUTO_add_hint ( plint, "dataset to be used for computation" );
100
101 /* fourth input : optional output file */
102
103 PLUTO_add_option( plint, "Output", "ofile_st", FALSE );
104 PLUTO_add_string( plint, "Outfile", 0, NULL, 0 );
105 PLUTO_add_hint ( plint, "file for statistical output" );
106 PLUTO_add_string( plint, "Overwrite", gr_num_yn_strings, gr_yn_strings, 1 );
107 PLUTO_add_hint ( plint, "option to overwrite output file" );
108
109 /* fifth input : minimum cutoff */
110 PLUTO_add_option( plint, "Cutoff", "min_st", FALSE );
111 PLUTO_add_number( plint, "Min", -10000, 10000, 1, 0, 1 );
112 PLUTO_add_hint ( plint, "exclude values below this cutoff" );
113
114 /* sixth input : maximum cutoff */
115 PLUTO_add_option( plint, "Cutoff", "max_st", FALSE );
116 PLUTO_add_number( plint, "Max", -10000, 10000, 1, 0, 1 );
117 PLUTO_add_hint ( plint, "exclude values above this cutoff" );
118
119 /* seventh input : tails */
120 PLUTO_add_option( plint, "Tails", "tails_st", FALSE );
121 PLUTO_add_hint ( plint, "apply min and max as tail cutoffs" );
122
123 /* eighth input : number of bins for histogram */
124 PLUTO_add_option( plint, "Histogram", "bins_st", FALSE );
125 PLUTO_add_number( plint, "Bins", 1, 1000, 0, 20, 1 );
126 PLUTO_add_hint ( plint, "number of bins in histogram" );
127
128 return plint;
129 }
130
131 char *
MASKCALC_main(PLUGIN_interface * plint)132 MASKCALC_main ( PLUGIN_interface * plint )
133 {
134 r_afni_s A;
135 mask_opt_s M;
136 char * ret_string;
137
138 memset( &A, 0, sizeof( A ) );
139 memset( &M, 0, sizeof( M ) );
140
141 if ( ( ret_string = process_args( &A, &M, plint ) ) != NULL )
142 return( ret_string );
143
144 return( process( &A, &M ) );
145 }
146
147
148
149 /*--------------------------------------------------------------------
150 Read the arguments, load the afni and mask structures.
151 --------------------------------------------------------------------
152 */
153 static char *
process_args(r_afni_s * A,mask_opt_s * M,PLUGIN_interface * plint)154 process_args(
155 r_afni_s * A,
156 mask_opt_s * M,
157 PLUGIN_interface * plint
158 )
159 {
160 MCW_idcode * idc;
161 char * tag, * str;
162 int op_val;
163
164 A->must_be_short = 1;
165 A->want_floats = 0;
166 A->subs_must_equal = 0;
167 A->max_subs = 0;
168
169 if ( plint == NULL )
170 return "----------------------\n"
171 "arguments : NULL input\n"
172 "----------------------\n";
173
174 while ( ( tag = PLUTO_get_optiontag( plint ) ) != NULL )
175 {
176 if ( ! strcmp( tag, "op_st" ) )
177 {
178 str = PLUTO_get_string( plint );
179 op_val = 1 + PLUTO_string_index( str, gr_num_ops, gr_op_strings );
180 if ( ( op_val <= (int)no_op ) || ( op_val >= (int)last_op ) )
181 {
182 sprintf( grMessage, "-------------------\n"
183 "Illegal operation : '%s'\n"
184 "Value is : %d\n"
185 "-------------------\n",
186 str, M->operation );
187 return grMessage;
188 }
189 M->operation = (op_enum)op_val;
190 continue;
191 }
192 else if ( ! strcmp( tag, "mask_st" ) )
193 {
194 idc = PLUTO_get_idcode( plint );
195 A->dset[0] = PLUTO_find_dset( idc );
196 if ( A->dset[0] == NULL )
197 return "----------------------\n"
198 "arg : bad mask dataset\n"
199 "----------------------";
200 DSET_load( A->dset[0] );
201 A->num_dsets++;
202 continue;
203 }
204 else if ( ! strcmp( tag, "dset_st" ) )
205 {
206 idc = PLUTO_get_idcode( plint );
207 A->dset[1] = PLUTO_find_dset( idc );
208 if ( A->dset[1] == NULL )
209 return "-----------------------\n"
210 "arg : bad inupt dataset\n"
211 "-----------------------";
212 DSET_load( A->dset[1] );
213 A->num_dsets++;
214 continue;
215 }
216 else if ( ! strcmp( tag, "ofile_st" ) )
217 {
218 M->outfile = PLUTO_get_string( plint );
219 str = PLUTO_get_string( plint );
220 if ( ( *str != 'y' ) && file_exists( M->outfile, "" ) )
221 {
222 sprintf( grMessage,
223 "-------------------------------\n"
224 "output file '%s' already exists\n"
225 "consider the 'overwrite' option\n"
226 "-------------------------------", M->outfile );
227 return( grMessage );
228 }
229 continue;
230 }
231 else if ( ! strcmp( tag, "min_st" ) )
232 {
233 M->min = PLUTO_get_number( plint );
234 M->use_min = 1;
235 continue;
236 }
237
238 if ( ! strcmp( tag, "max_st" ) )
239 {
240 M->max = PLUTO_get_number( plint );
241 M->use_max = 1;
242 continue;
243 }
244
245 if ( ! strcmp( tag, "tails_st" ) )
246 {
247 M->use_tails = 1; /* need to check with min/max */
248 continue;
249 }
250
251 if ( ! strcmp( tag, "bins_st" ) )
252 {
253 M->num_bins = (int)PLUTO_get_number( plint );
254 if ( ( M->num_bins <= 0 ) || ( M->num_bins > R_MAX_BINS ) )
255 {
256 sprintf( grMessage, "-----------------------------\n"
257 "Illegal number of bins : %d\n"
258 "(must be in range [1,%d])\n"
259 "-----------------------------",
260 M->num_bins, R_MAX_BINS );
261 return grMessage;
262 }
263
264 continue;
265 }
266
267 /* we should not get to this point */
268
269 sprintf( grMessage, "-----------------------\n"
270 "Unknown optiontag : %s\n"
271 "-----------------------", tag );
272 return grMessage;
273 }
274
275 if ( M->use_tails && ( ! M->use_min || ! M->use_max ) )
276 {
277 sprintf( grMessage, "------------------------------------------\n"
278 "'tails' option requires min and max values\n"
279 "------------------------------------------" );
280 return grMessage;
281 }
282
283 if ( M->num_bins && ( M->operation != hist_op ) )
284 {
285 return "----------------------------------------------------\n"
286 "choosing # bins applies only to the 'hist' operation\n"
287 "----------------------------------------------------";
288 }
289 else if ( ! M->num_bins )
290 M->num_bins = 20;
291
292 if ( ( str = fill_afni_struct( A ) ) != NULL )
293 return str;
294
295 if ( M->outfile && *M->outfile )
296 {
297 if ( ( M->outfp = open_file( M->outfile, "w" ) ) == NULL )
298 {
299 sprintf( grMessage, "--------------------------------\n"
300 "Failed to open '%s' for writing.\n"
301 "--------------------------------",
302 M->outfile );
303 return grMessage;
304 }
305 }
306 else
307 M->outfp = stdout;
308
309 return NULL;
310 }
311
312
313 /***********************************************************************
314 **
315 ** Perform computation, storing the result back in image1 (space).
316 **
317 ** return 1 - successful completion
318 ** 0 - failure
319 **
320 ************************************************************************
321 */
322 static char *
process(r_afni_s * A,mask_opt_s * M)323 process( r_afni_s * A, mask_opt_s * M )
324 {
325 char * ret_string;
326
327 switch ( M->operation )
328 {
329 case hist_op:
330
331 ret_string = calc_hist( A, M );
332
333 break;
334
335 case stats_op:
336
337 ret_string = calc_stats( A, M );
338
339 break;
340
341 default:
342
343 ret_string = "--------------------\n"
344 "Error: maskcalc_p_00\n"
345 "Invalid operation.\n"
346 "--------------------";
347 } /* end switch */
348
349 PURGE_DSET( A->dset[0] );
350 PURGE_DSET( A->dset[1] );
351
352 if ( M->outfp != stdout )
353 fclose( M->outfp );
354
355 return ret_string;
356 }
357
358
359 /***********************************************************************
360 **
361 ** Mann-Whitney U-test for comparison of two arbitrary datasets.
362 **
363 ************************************************************************
364 */
365 static long
get_mask_size(r_afni_s * A,int index,int subbrick)366 get_mask_size(
367 r_afni_s * A,
368 int index, /* index into simage array */
369 int subbrick
370 )
371 {
372 long count, size;
373 short * ptr;
374
375 for ( count = 0, size = 0, ptr = A->simage[ index ][ subbrick ];
376 count < A->nvox;
377 count++, ptr++ )
378 if ( *ptr )
379 size++;
380
381 return size;
382 }
383
384
385
386 /***********************************************************************
387 **
388 ** Comparison function for two shorts, for use in qsort.
389 **
390 ************************************************************************
391 */
short_test(const void * p1,const void * p2)392 int short_test(
393 const void * p1,
394 const void * p2
395 )
396 {
397 short * s1 = ( short * )p1;
398 short * s2 = ( short * )p2;
399
400 if ( *s1 < *s2 )
401 return -1;
402
403 return ( *s1 > *s2 ); /* if >, return 1, else ==, so 0 */
404 }
405
406
407
408
409 /***********************************************************************
410 **
411 ** Output histogram from masked image.
412 **
413 ************************************************************************
414 */
415 static char *
calc_hist(r_afni_s * A,mask_opt_s * M)416 calc_hist( r_afni_s * A, mask_opt_s * M )
417 {
418 float * data;
419 float * ptr;
420 float bin_size, cum, junk;
421
422 long size, new_size = 0;
423 long count;
424 int * bins;
425
426 int places; /* decimal places in output */
427
428
429 if (( data = (float *)malloc(A->subs[1] * A->nvox * sizeof(float))) == NULL)
430 {
431 sprintf( grMessage, "Error: maskcalc_ch_00\n"
432 "Failed to allocate memory for %d floats.\n",
433 A->nvox * A->subs[1] );
434 return grMessage;
435 }
436
437 if ( ( size = mask_all_shorts_to_float( A, 1, 0, data ) ) == 0 )
438 {
439 sprintf( grMessage, "Error: 5090\n"
440 "Masking shorts results in empty array.\n" );
441 return grMessage;
442 }
443
444 if ( !M->use_min && !M->use_max )
445 assign_min_max( data, size, &M->min, &M->max );
446 else if ( !M->use_max )
447 {
448 assign_min_max( data, size, &junk, &M->max );
449
450 if ( M->min > M->max )
451 {
452 sprintf( grMessage, "Error: maskcalc_ch_10\n"
453 "Min of %f is greater than max of %f\n",
454 M->min, M->max );
455 return grMessage;
456 }
457 }
458
459 junk = ( fabsf( M->max ) > fabsf( M->min ) ) ?
460 fabsf( M->max ) : fabsf( M->min );
461
462 if ( junk == 0 )
463 places = 2;
464 else
465 {
466 places = 4 - ( int )floor( log10( junk ) );
467
468 if ( places > 7 )
469 places = 7;
470 else if ( places < 0 )
471 places = 0;
472 }
473
474 if ( ( bins = (int *)calloc( M->num_bins, sizeof( int ) ) ) == NULL )
475 {
476 sprintf( grMessage, "Error: maskcalc_ch_30\n"
477 "Failed to allocate for %d longs.\n", M->num_bins );
478 return grMessage;
479 }
480
481 bin_size = ( M->max - M->min ) / M->num_bins;
482 bin_size += 0.000001 * bin_size;
483 if ( bin_size == 0.0 )
484 bin_size = 1.0e-34;
485
486 for ( count = 0, ptr = data; count < size; count++, ptr++ )
487 if ( ( *ptr <= M->max ) && ( *ptr >= M->min ) )
488 {
489 bins[ ( int )( ( *ptr - M->min ) / bin_size ) ]++;
490 new_size++;
491 }
492
493 if ( new_size == 0 )
494 new_size = 1;
495
496
497 fprintf( M->outfp, "\n range \t #vox \t %% \t cum %%\n");
498 fprintf( M->outfp, "------------------- \t--------\t------\t-------\n");
499
500 cum = 0.0;
501 for ( count = 0; count < M->num_bins; count++ )
502 {
503 cum += 100.0 * bins[ count ] / new_size;
504
505 fprintf( M->outfp, "[%8.*f,%8.*f) \t%8d\t%6.3f\t%7.3f\n",
506 places, M->min + count * bin_size,
507 places, M->min + (count+1) * bin_size,
508 bins[ count ],
509 100.0 * bins[ count ] / new_size,
510 cum );
511 }
512 fputc( '\n', M->outfp );
513
514 return NULL;
515 }
516
517
518 /***********************************************************************
519 **
520 ** Assign min and max values to global variables.
521 **
522 ************************************************************************
523 */
524 static void
assign_min_max(float * data,long size,float * min,float * max)525 assign_min_max(
526 float * data,
527 long size,
528 float * min,
529 float * max
530 )
531 {
532 float * ptr = data;
533 long count;
534
535
536 *min = *data;
537 *max = *data;
538
539 for ( count = 1; count < size; count++, ptr++ )
540 {
541 if ( *ptr < *min )
542 *min = *ptr;
543
544 if ( *ptr > *max )
545 *max = *ptr;
546 }
547 }
548
549
550 /***********************************************************************
551 **
552 ** Output statistics from a masked image.
553 **
554 ************************************************************************
555 */
556 static char *
calc_stats(r_afni_s * A,mask_opt_s * M)557 calc_stats( r_afni_s * A, mask_opt_s * M )
558 {
559 float * data;
560 float min, max, savemin, savemax;
561
562 long size;
563 int sub;
564
565
566 if ( ( data = ( float * )malloc( A->nvox * sizeof(float) ) )
567 == NULL )
568 {
569 sprintf( grMessage, "Error: 5130\n"
570 "Failed to allocate memory for %d floats.\n",
571 A->nvox );
572 return grMessage;
573 }
574
575 print_stats_header( M->outfp );
576
577 for ( sub = 0; sub < A->subs[1]; sub++ )
578 {
579 if ( ( size = mask_shorts_to_float( A, data, sub, 0, 1 ) ) == 0 )
580 {
581 sprintf( grMessage, "Error: 5140\n"
582 "Masking shorts results in empty array.\n" );
583 return grMessage;
584 }
585
586 assign_min_max( data, size, &savemin, &savemax );
587
588 if ( ! M->use_min && ! M->use_max )
589 {
590 min = savemin; /* use actual min/max for data */
591 max = savemax;
592
593 do_stats( A, data, size, min, max, sub, M->outfp, NULL, NULL, NULL);
594 }
595 else if ( ! M->use_max ) /* so use_min is set */
596 {
597 min = M->min; /* use user input cutoff */
598 max = savemax;
599
600 if ( min <= max )
601 do_stats( A, data, size, min, max, sub, M->outfp,
602 NULL, NULL, NULL);
603 else
604 print_empty_stats( M->outfp );
605 }
606 else /* use_min AND use_max are set */
607 {
608 /* NOTE : we are using the tails here */
609
610 min = savemin;
611 max = M->min;
612
613 if ( min <= max )
614 do_stats( A, data, size, min, max, sub, M->outfp,
615 NULL, NULL, NULL);
616 else
617 print_empty_stats( M->outfp );
618
619 min = M->max;
620 max = savemax;
621
622 if ( min <= max )
623 do_stats( A, data, size, min, max, sub, M->outfp,
624 NULL, NULL, NULL);
625 else
626 print_empty_stats( M->outfp );
627 }
628 }
629
630 return NULL;
631 }
632
633
634 /***********************************************************************
635 **
636 ** Print stats header.
637 **
638 ************************************************************************
639 */
640 static void
print_stats_header(FILE * fp)641 print_stats_header ( FILE * fp )
642 {
643 fprintf( fp, " Sub\t vol \t%% BRIK\t min \t max "
644 "\t mean \t SEM \t STD \t 95%% "
645 "\t thresh # \t thresh %%\n" );
646 fprintf( fp, "-----\t-------\t------\t------\t------"
647 "\t------\t------\t------\t-------------"
648 "\t----------\t---------\n" );
649 }
650
651
652 /***********************************************************************
653 **
654 ** Print empty stats.
655 **
656 ************************************************************************
657 */
658 static void
print_empty_stats(FILE * fp)659 print_empty_stats ( FILE * fp )
660 {
661 fprintf( fp, " 0 \t 0 \t 0 \t 0 \t 0 "
662 "\t 0 \t 0 \t 0 \t ( 0, 0 ) "
663 "\t 0 \t 0 \n" );
664 }
665
666
667 /***********************************************************************
668 **
669 ** Output statistics for each masked subbrick. We are given
670 ** a min and max value to further restrict our brick.
671 **
672 ************************************************************************
673 */
674 static void
do_stats(r_afni_s * A,float * data,long size,float min,float max,int sub,FILE * fp,long * ret_size,float * average,float * ret_var)675 do_stats(
676 r_afni_s * A,
677 float * data, /* IN */
678 long size, /* IN */
679 float min, /* IN */
680 float max, /* IN */
681 int sub, /* IN */
682 FILE * fp, /* IN */
683 long * ret_size, /* OUT */
684 float * average, /* OUT */
685 float * ret_var /* OUT */
686 )
687 {
688 double sum_diff2 = 0.0;
689 double dtmp;
690
691 float mean, SEM, STD;
692 float * ptr = data;
693 float tmp;
694 float local_min = max, local_max = min;
695
696 long count;
697 long new_size;
698
699
700 /* start by getting the mean, min and max (across mask and in range) */
701
702 dtmp = 0.0;
703 new_size = 0;
704 for ( count = 0, ptr = data; count < size; count++, ptr++ )
705 if ( ( *ptr >= min ) && ( *ptr <= max ) )
706 {
707 new_size++;
708
709 dtmp += *ptr;
710
711 if ( *ptr < local_min )
712 local_min = *ptr;
713
714 if ( *ptr > local_max )
715 local_max = *ptr;
716 }
717
718 if ( new_size > 0 )
719 mean = dtmp / new_size;
720 else
721 mean = 0;
722
723 /* now get variance */
724 sum_diff2 = 0.0;
725 for ( count = 0, ptr = data; count < size; count++, ptr++ )
726 if ( ( *ptr >= min ) && ( *ptr <= max ) )
727 {
728 tmp = *ptr - mean;
729 sum_diff2 += tmp * tmp;
730 }
731
732 if ( new_size < 2 )
733 {
734 STD = 0;
735 SEM = 0;
736 }
737 else
738 {
739 STD = sqrt( sum_diff2 / ( new_size - 1 ) );
740 SEM = STD / sqrt( new_size );
741 }
742
743 fprintf( fp,
744 "%5d\t%7ld\t %5.*f\t%6.*f\t%6.*f\t%6.*f\t%6.*f\t%6.*f\t(%-5.*f, %5.*f)"
745 "\t %8ld \t %6.*f\n",
746 sub, size,
747 num_places( 100.0*size/A->nvox, 5 ), 100.0*size/A->nvox,
748 num_places( local_min, 6 ), local_min,
749 num_places( local_max, 6 ), local_max,
750 num_places( mean, 6 ), mean,
751 num_places( SEM, 6 ), SEM,
752 num_places( STD, 6 ), STD,
753 num_places( mean-1.96*SEM, 5 ), mean-1.96*SEM,
754 num_places( mean+1.96*SEM, 5 ), mean+1.96*SEM,
755 new_size,
756 num_places( 100.0*new_size/size, 6 ), 100.0*new_size/size );
757
758 if ( ret_size != NULL )
759 *ret_size = new_size;
760
761 if ( average != NULL )
762 *average = mean;
763
764 if ( ret_var != NULL )
765 *ret_var = STD*STD;
766 }
767
768
769 /***********************************************************************
770 **
771 ** Calculate the number of decimal places needed for a float
772 ** given the magnitude and the total space.
773 **
774 ************************************************************************
775 */
776 static int
num_places(float num,int size)777 num_places(
778 float num,
779 int size
780 )
781 {
782 float junk;
783 int places;
784
785 junk = fabsf( num );
786
787 junk = ( junk == 0 ) ? 1.0 : junk;
788
789 /*
790 ** Allow for at least one place to the left of the decimal,
791 ** and the decimal itself.
792 */
793
794 places = size - 2 - ( int )floor( log10( junk ) );
795
796 if ( places < 0 )
797 places = 0;
798 else if ( places >= size )
799 places = size - 1;
800
801 return places;
802 }
803
804
805 /***********************************************************************
806 **
807 ** Copy masked and scaled shorts into float array.
808 **
809 ** Return the number of floats.
810 **
811 ************************************************************************
812 */
813 static long
mask_shorts_to_float(r_afni_s * A,float * data,int sub,int mask_index,int data_index)814 mask_shorts_to_float(
815 r_afni_s * A,
816 float * data, /* output data location */
817 int sub, /* current subbrick */
818 int mask_index, /* index of mask */
819 int data_index /* index of data */
820 )
821 {
822 float * fptr; /* floats */
823 short * mptr; /* mask */
824 short * sptr; /* shorts */
825
826 float factor;
827 long count; /* voxel counter */
828
829
830 fptr = data; /* floats */
831
832 if ( A->subs[ mask_index ] > 1 )
833 mptr = A->simage[ mask_index ][ sub ]; /* mask */
834 else
835 mptr = A->simage[ mask_index ][ 0 ]; /* mask */
836
837 sptr = A->simage[ data_index ][ sub ]; /* shorts */
838
839 /*
840 ** Note that mptr and sptr get incremented continuously, where
841 ** fptr gets incremented only when we have a good mask value;
842 */
843
844 if ( ( factor = A->factor[ data_index ][ sub ] ) == 1 )
845 {
846 for( count = 0; count < A->nvox; count++, mptr++, sptr++ )
847 if ( *mptr )
848 *fptr++ = *sptr;
849 }
850 else
851 {
852 for( count = 0; count < A->nvox; count++, mptr++, sptr++ )
853 if ( *mptr )
854 *fptr++ = *sptr * factor;
855 }
856
857
858 return( ( long )( fptr - data ) ); /* return # of floats */
859 }
860
861
862 /***********************************************************************
863 **
864 ** Copy masked shorts into separate array.
865 **
866 ** Return the number copied.
867 **
868 ************************************************************************
869 */
870 static long
mask_shorts_to_short(r_afni_s * A,short * data,int sub,int mask_index,int data_index)871 mask_shorts_to_short(
872 r_afni_s * A,
873 short * data, /* output data location */
874 int sub, /* current subbrick */
875 int mask_index, /* index of mask */
876 int data_index /* index of data */
877 )
878 {
879 short * dptr; /* destination */
880 short * mptr; /* mask */
881 short * sptr; /* shorts */
882
883 long count; /* voxel counter */
884
885
886 dptr = data; /* destination pointer */
887
888 if ( A->subs[ mask_index ] > 1 )
889 mptr = A->simage[ mask_index ][ sub ]; /* mask */
890 else
891 mptr = A->simage[ mask_index ][ 0 ]; /* mask */
892
893 sptr = A->simage[ data_index ][ sub ]; /* shorts */
894
895 /*
896 ** Note that mptr and sptr get incremented continuously, where
897 ** dptr gets incremented only when we have a good mask value;
898 */
899
900 for( count = 0; count < A->nvox; count++, mptr++, sptr++ )
901 if ( *mptr )
902 *dptr++ = *sptr;
903
904 return( ( long )( dptr - data ) );
905 }
906
907
908 /***********************************************************************
909 **
910 ** Copy masked and scaled shorts into float array.
911 **
912 ** Return the number of floats.
913 **
914 ************************************************************************
915 */
916 static long
mask_all_shorts_to_float(r_afni_s * A,int data_index,int mask_index,float * data)917 mask_all_shorts_to_float(
918 r_afni_s * A,
919 int data_index,
920 int mask_index,
921 float * data
922 )
923 {
924 float * fptr; /* floats */
925 short * mptr; /* mask */
926 short * sptr; /* shorts */
927
928 float factor;
929 long count; /* voxel counter */
930 int sub; /* sub-brick */
931
932
933 fptr = data; /* floats */
934
935 for ( sub = 0; sub < A->subs[data_index]; sub ++ )
936 {
937 if ( A->subs[mask_index] == A->subs[data_index] )
938 mptr = A->simage[ mask_index ][ sub ]; /* mask */
939 else
940 mptr = A->simage[ mask_index ][ 0 ];
941
942 sptr = A->simage[ data_index ][ sub ]; /* shorts */
943
944 /*
945 ** Note that mptr and sptr get incremented continuously, where
946 ** fptr gets incremented only when we have a good mask value;
947 */
948
949 if ( ( factor = A->factor[ data_index ][ sub ] ) == 1 )
950 {
951 for( count = 0; count < A->nvox; count++, mptr++, sptr++ )
952 if ( *mptr )
953 *fptr++ = *sptr;
954 }
955 else
956 {
957 for( count = 0; count < A->nvox; count++, mptr++, sptr++ )
958 if ( *mptr )
959 *fptr++ = *sptr * factor;
960 }
961 }
962
963 return( ( long )( fptr - data ) ); /* return # of floats */
964 }
965
966
967 /*----------------------------------------------------------------------
968 **
969 ** Open a file in the given mode (just to shorten code).
970 **
971 **----------------------------------------------------------------------
972 */
973 static FILE *
open_file(char * file,char * mode)974 open_file(
975 char * file,
976 char * mode
977 )
978 {
979 return fopen( file, mode );
980 }
981
982
983 /***********************************************************************
984 **
985 ** Check for existence of output file.
986 **
987 ** If suffix is NULL, only consider filename.
988 ** If filename already ends in the suffix, only consider filname.
989 ** If the filename ends in a period, only add the suffix.
990 ** Otherwise add a period and the suffix before checking its existence.
991 **
992 ** return 1 - exists
993 ** 0 - does not exist
994 **
995 ************************************************************************
996 */
997 static int
file_exists(char * filename,char * suffix)998 file_exists(
999 char * filename,
1000 char * suffix
1001 )
1002 {
1003 struct stat buf;
1004 char file[ R_FILE_L + 6 ] = "";
1005 char * filep = file;
1006
1007 if ( suffix == NULL ) /* we don't need to worry about memory */
1008 filep = filename;
1009 else if ( ! strcmp( suffix, filename+strlen(filename)-strlen(suffix) ) )
1010 strcpy( file, filename );
1011 else if ( filename[ strlen( filename ) - 1 ] == '.' )
1012 sprintf( file, "%s%s", filename, suffix );
1013 else
1014 sprintf( file, "%s.%s", filename, suffix );
1015
1016 /* stat returns 0 on existence */
1017 return ( stat( filep, &buf ) == 0 );
1018 }
1019
1020
1021 /***********************************************************************
1022 **
1023 ** Use the dset variable to fill the rest of the structure.
1024 **
1025 ************************************************************************
1026 */
1027 static char *
fill_afni_struct(r_afni_s * A)1028 fill_afni_struct( r_afni_s * A )
1029 {
1030 u_short mus;
1031 int sub, brick;
1032
1033 for ( brick = 0; brick < A->num_dsets; brick++ )
1034 {
1035 A->subs[ brick ] = DSET_NVALS( A->dset[ brick ] );
1036
1037 if ( A->max_subs && ( A->subs[brick] > A->max_subs ) )
1038 {
1039 sprintf( grMessage, "------------------------------------\n"
1040 "Error: maskcalc_fas_00\n"
1041 "Brick #%d contains %d sub-bricks.\n"
1042 "The limit is %d.\n"
1043 "------------------------------------",
1044 brick, A->subs[brick], A->max_subs );
1045 return grMessage;
1046 }
1047
1048 if ( A->subs_must_equal && ( A->subs[brick] != A->subs[0] ) )
1049 {
1050 sprintf( grMessage, "------------------------------------\n"
1051 "Error: maskcalc_fas_02\n"
1052 "Brick #%d contains %d sub-bricks.\n"
1053 "Brick #%d contains %d sub-bricks.\n"
1054 "We are requiring them to be equal.\n"
1055 "------------------------------------",
1056 0, A->subs[0],
1057 brick, A->subs[brick] );
1058 return grMessage;
1059 }
1060
1061 if ( ( A->simage[brick] =
1062 (short **)malloc( A->subs[brick]*sizeof(short *)) ) == NULL )
1063 {
1064 return "-------------------------\n"
1065 "Error: maskcalc_fas_05\n"
1066 "memory allocation failure\n"
1067 "-------------------------";
1068 }
1069 if ( ( A->factor[brick] =
1070 (float *)malloc( A->subs[brick]*sizeof(float)) ) == NULL)
1071 {
1072 return "-------------------------\n"
1073 "Error: maskcalc_fas_10\n"
1074 "memory allocation failure\n"
1075 "-------------------------";
1076 }
1077
1078 for ( sub = 0; sub < A->subs[brick]; sub++ )
1079 {
1080 A->simage[brick][sub] = (short *)DSET_ARRAY(A->dset[brick],sub);
1081 A->factor[brick][sub] = DSET_BRICK_FACTOR(A->dset[brick],sub);
1082 if ( A->factor[brick][sub] == 0.0 )
1083 A->factor[brick][sub] = 1.0;
1084 }
1085
1086 if ( brick == 0 )
1087 {
1088 A->nx = A->dset[brick]->daxes->nxx;
1089 A->ny = A->dset[brick]->daxes->nyy;
1090 A->nz = A->dset[brick]->daxes->nzz;
1091 A->nvox = A->nx * A->ny * A->nz;
1092 }
1093 else if ( ( A->nx != A->dset[brick]->daxes->nxx ) ||
1094 ( A->ny != A->dset[brick]->daxes->nyy ) ||
1095 ( A->nz != A->dset[brick]->daxes->nzz ) )
1096 {
1097 sprintf( grMessage,
1098 "--------------------------------\n"
1099 "Error : maskcalc_fas_20\n"
1100 "Unaligned dimensions.\n"
1101 "(%d,%d,%d) != (%d,%d,%d)\n"
1102 "--------------------------------",
1103 A->dset[brick]->daxes->nxx, A->dset[brick]->daxes->nyy,
1104 A->dset[brick]->daxes->nzz, A->nx, A->ny, A->nz );
1105 return grMessage;
1106 }
1107 }
1108
1109 if ( A->want_floats && ! assign_afni_floats( A ) )
1110 return "-----------------------------\n"
1111 "Error: maskcalc_fas_30\n"
1112 "Failed to create afni fimage.\n"
1113 "-----------------------------";
1114
1115 mus = 0;
1116 A->max_u_short = 0;
1117 for ( brick = 1; brick < A->num_dsets; brick++ )
1118 for ( sub = 0; sub < A->subs[brick]; sub++ )
1119 {
1120 mus = r_get_max_u_short( (u_short *)A->simage[brick][sub], A->nvox);
1121 if ( mus > A->max_u_short )
1122 A->max_u_short = mus;
1123 }
1124
1125 return NULL;
1126 }
1127
1128
1129 /***********************************************************************
1130 **
1131 ** Create a float brick corresponding to the first subbrick of
1132 ** every non-mask brick.
1133 **
1134 ************************************************************************
1135 */
1136 static int
assign_afni_floats(r_afni_s * A)1137 assign_afni_floats( r_afni_s * A )
1138 {
1139 short * sptr;
1140 float * fptr;
1141 float factor = A->factor[1][0];
1142 int count;
1143
1144 /* at this point, only brick 1 is a non-mask brick */
1145 if ( ( A->fimage[1] = (float *)malloc( A->nvox * sizeof(float))) == NULL )
1146 return 0;
1147
1148 for ( count = 0, fptr = A->fimage[1], sptr = A->simage[1][0];
1149 count < A->nvox;
1150 count++, fptr++, sptr++ )
1151 *fptr = factor * *sptr;
1152
1153 return 1;
1154 }
1155
1156
1157 /***********************************************************************
1158 **
1159 ** Calculate the maximum unsigned short in the array.
1160 **
1161 ************************************************************************
1162 */
1163 static u_short
r_get_max_u_short(u_short * S,int size)1164 r_get_max_u_short( u_short * S, int size )
1165 {
1166 u_short * usptr, max = *S;
1167 int c = 0;
1168
1169 for ( c = 0, usptr = S; c < size; c++, usptr++ )
1170 {
1171 if ( *usptr > max )
1172 max = *usptr;
1173 }
1174
1175 return max;
1176 }
1177
1178
1179