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 #include <stdio.h>
00036 #include <math.h>
00037 #include <stdlib.h>
00038 #if !(defined(DARWIN) || defined(FreeBSD))
00039 # include <malloc.h>
00040 #endif
00041 #include <string.h>
00042 #include <sys/types.h>
00043 #include <sys/stat.h>
00044
00045 #include "mrilib.h"
00046 #include "afni.h"
00047 #include "plug_maskcalc.h"
00048
00049 #ifndef ALLOW_PLUGINS
00050 # error "Plugins not properly set up -- see machdep.h"
00051 #endif
00052
00053 char * MASKCALC_main( PLUGIN_interface * );
00054
00055 static char grMessage[ R_MESSAGE_L ];
00056
00057 static char * gr_help_message =
00058 "maskcalc plugin - rickr";
00059
00060
00061 DEFINE_PLUGIN_PROTOTYPE
00062
00063 PLUGIN_interface * PLUGIN_init( int ncall )
00064 {
00065 PLUGIN_interface * plint;
00066
00067 if ( ncall > 0 )
00068 return NULL;
00069
00070 plint = PLUTO_new_interface( "maskcalc", "masked computations on datasets",
00071 gr_help_message, PLUGIN_CALL_VIA_MENU, MASKCALC_main );
00072
00073 PLUTO_add_hint( plint, "Wouldn't some cookies be right tasty?" );
00074
00075 PLUTO_set_sequence( plint , "z:Reynolds" ) ;
00076
00077
00078
00079 PLUTO_add_option( plint, "Function", "op_st", TRUE );
00080 PLUTO_add_hint ( plint, "function to perform on the data" );
00081 PLUTO_add_string( plint, "Operation", gr_num_ops, gr_op_strings, 0 );
00082 PLUTO_add_hint ( plint, "function to perform on the data" );
00083
00084
00085
00086 PLUTO_add_option ( plint, "Dataset", "mask_st", TRUE );
00087 PLUTO_add_hint ( plint, "dataset to be used as mask" );
00088 PLUTO_add_dataset( plint, "Mask", ANAT_ALL_MASK , FUNC_ALL_MASK,
00089 DIMEN_3D_MASK | DIMEN_4D_MASK | BRICK_SHORT_MASK );
00090 PLUTO_add_hint ( plint, "dataset to be used as mask" );
00091
00092
00093
00094 PLUTO_add_option ( plint, "Dataset", "dset_st", TRUE );
00095 PLUTO_add_hint ( plint, "computational dataset" );
00096 PLUTO_add_dataset( plint, "Dset", ANAT_ALL_MASK, FUNC_ALL_MASK,
00097 DIMEN_ALL_MASK | BRICK_SHORT_MASK );
00098 PLUTO_add_hint ( plint, "dataset to be used for computation" );
00099
00100
00101
00102 PLUTO_add_option( plint, "Output", "ofile_st", FALSE );
00103 PLUTO_add_string( plint, "Outfile", 0, NULL, 0 );
00104 PLUTO_add_hint ( plint, "file for statistical output" );
00105 PLUTO_add_string( plint, "Overwrite", gr_num_yn_strings, gr_yn_strings, 1 );
00106 PLUTO_add_hint ( plint, "option to overwrite output file" );
00107
00108
00109 PLUTO_add_option( plint, "Cutoff", "min_st", FALSE );
00110 PLUTO_add_number( plint, "Min", -10000, 10000, 1, 0, 1 );
00111 PLUTO_add_hint ( plint, "exclude values below this cutoff" );
00112
00113
00114 PLUTO_add_option( plint, "Cutoff", "max_st", FALSE );
00115 PLUTO_add_number( plint, "Max", -10000, 10000, 1, 0, 1 );
00116 PLUTO_add_hint ( plint, "exclude values above this cutoff" );
00117
00118
00119 PLUTO_add_option( plint, "Tails", "tails_st", FALSE );
00120 PLUTO_add_hint ( plint, "apply min and max as tail cutoffs" );
00121
00122
00123 PLUTO_add_option( plint, "Histogram", "bins_st", FALSE );
00124 PLUTO_add_number( plint, "Bins", 1, 1000, 0, 20, 1 );
00125 PLUTO_add_hint ( plint, "number of bins in histogram" );
00126
00127 return plint;
00128 }
00129
00130 char *
00131 MASKCALC_main ( PLUGIN_interface * plint )
00132 {
00133 r_afni_s A;
00134 mask_opt_s M;
00135 char * ret_string;
00136
00137 memset( &A, 0, sizeof( A ) );
00138 memset( &M, 0, sizeof( M ) );
00139
00140 if ( ( ret_string = process_args( &A, &M, plint ) ) != NULL )
00141 return( ret_string );
00142
00143 return( process( &A, &M ) );
00144 }
00145
00146
00147
00148
00149
00150
00151
00152 static char *
00153 process_args(
00154 r_afni_s * A,
00155 mask_opt_s * M,
00156 PLUGIN_interface * plint
00157 )
00158 {
00159 MCW_idcode * idc;
00160 char * tag, * str;
00161 int op_val;
00162
00163 A->must_be_short = 1;
00164 A->want_floats = 0;
00165 A->subs_must_equal = 0;
00166 A->max_subs = 0;
00167
00168 if ( plint == NULL )
00169 return "----------------------\n"
00170 "arguments : NULL input\n"
00171 "----------------------\n";
00172
00173 while ( ( tag = PLUTO_get_optiontag( plint ) ) != NULL )
00174 {
00175 if ( ! strcmp( tag, "op_st" ) )
00176 {
00177 str = PLUTO_get_string( plint );
00178 op_val = 1 + PLUTO_string_index( str, gr_num_ops, gr_op_strings );
00179 if ( ( op_val <= (int)no_op ) || ( op_val >= (int)last_op ) )
00180 {
00181 sprintf( grMessage, "-------------------\n"
00182 "Illegal operation : '%s'\n"
00183 "Value is : %d\n"
00184 "-------------------\n",
00185 str, M->operation );
00186 return grMessage;
00187 }
00188 M->operation = (op_enum)op_val;
00189 continue;
00190 }
00191 else if ( ! strcmp( tag, "mask_st" ) )
00192 {
00193 idc = PLUTO_get_idcode( plint );
00194 A->dset[0] = PLUTO_find_dset( idc );
00195 if ( A->dset[0] == NULL )
00196 return "----------------------\n"
00197 "arg : bad mask dataset\n"
00198 "----------------------";
00199 DSET_load( A->dset[0] );
00200 A->num_dsets++;
00201 continue;
00202 }
00203 else if ( ! strcmp( tag, "dset_st" ) )
00204 {
00205 idc = PLUTO_get_idcode( plint );
00206 A->dset[1] = PLUTO_find_dset( idc );
00207 if ( A->dset[1] == NULL )
00208 return "-----------------------\n"
00209 "arg : bad inupt dataset\n"
00210 "-----------------------";
00211 DSET_load( A->dset[1] );
00212 A->num_dsets++;
00213 continue;
00214 }
00215 else if ( ! strcmp( tag, "ofile_st" ) )
00216 {
00217 M->outfile = PLUTO_get_string( plint );
00218 str = PLUTO_get_string( plint );
00219 if ( ( *str != 'y' ) && file_exists( M->outfile, "" ) )
00220 {
00221 sprintf( grMessage,
00222 "-------------------------------\n"
00223 "output file '%s' already exists\n"
00224 "consider the 'overwrite' option\n"
00225 "-------------------------------", M->outfile );
00226 return( grMessage );
00227 }
00228 continue;
00229 }
00230 else if ( ! strcmp( tag, "min_st" ) )
00231 {
00232 M->min = PLUTO_get_number( plint );
00233 M->use_min = 1;
00234 continue;
00235 }
00236
00237 if ( ! strcmp( tag, "max_st" ) )
00238 {
00239 M->max = PLUTO_get_number( plint );
00240 M->use_max = 1;
00241 continue;
00242 }
00243
00244 if ( ! strcmp( tag, "tails_st" ) )
00245 {
00246 M->use_tails = 1;
00247 continue;
00248 }
00249
00250 if ( ! strcmp( tag, "bins_st" ) )
00251 {
00252 M->num_bins = (int)PLUTO_get_number( plint );
00253 if ( ( M->num_bins <= 0 ) || ( M->num_bins > R_MAX_BINS ) )
00254 {
00255 sprintf( grMessage, "-----------------------------\n"
00256 "Illegal number of bins : %d\n"
00257 "(must be in range [1,%d])\n"
00258 "-----------------------------",
00259 M->num_bins, R_MAX_BINS );
00260 return grMessage;
00261 }
00262
00263 continue;
00264 }
00265
00266
00267
00268 sprintf( grMessage, "-----------------------\n"
00269 "Unknown optiontag : %s\n"
00270 "-----------------------", tag );
00271 return grMessage;
00272 }
00273
00274 if ( M->use_tails && ( ! M->use_min || ! M->use_max ) )
00275 {
00276 sprintf( grMessage, "------------------------------------------\n"
00277 "'tails' option requires min and max values\n"
00278 "------------------------------------------" );
00279 return grMessage;
00280 }
00281
00282 if ( M->num_bins && ( M->operation != hist_op ) )
00283 {
00284 return "----------------------------------------------------\n"
00285 "choosing # bins applies only to the 'hist' operation\n"
00286 "----------------------------------------------------";
00287 }
00288 else if ( ! M->num_bins )
00289 M->num_bins = 20;
00290
00291 if ( ( str = fill_afni_struct( A ) ) != NULL )
00292 return str;
00293
00294 if ( M->outfile && *M->outfile )
00295 {
00296 if ( ( M->outfp = open_file( M->outfile, "w" ) ) == NULL )
00297 {
00298 sprintf( grMessage, "--------------------------------\n"
00299 "Failed to open '%s' for writing.\n"
00300 "--------------------------------",
00301 M->outfile );
00302 return grMessage;
00303 }
00304 }
00305 else
00306 M->outfp = stdout;
00307
00308 return NULL;
00309 }
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321 static char *
00322 process( r_afni_s * A, mask_opt_s * M )
00323 {
00324 char * ret_string;
00325
00326 switch ( M->operation )
00327 {
00328 case hist_op:
00329
00330 ret_string = calc_hist( A, M );
00331
00332 break;
00333
00334 case stats_op:
00335
00336 ret_string = calc_stats( A, M );
00337
00338 break;
00339
00340 default:
00341
00342 ret_string = "--------------------\n"
00343 "Error: maskcalc_p_00\n"
00344 "Invalid operation.\n"
00345 "--------------------";
00346 }
00347
00348 PURGE_DSET( A->dset[0] );
00349 PURGE_DSET( A->dset[1] );
00350
00351 if ( M->outfp != stdout )
00352 fclose( M->outfp );
00353
00354 return ret_string;
00355 }
00356
00357
00358
00359
00360
00361
00362
00363
00364 static long
00365 get_mask_size(
00366 r_afni_s * A,
00367 int index,
00368 int subbrick
00369 )
00370 {
00371 long count, size;
00372 short * ptr;
00373
00374 for ( count = 0, size = 0, ptr = A->simage[ index ][ subbrick ];
00375 count < A->nvox;
00376 count++, ptr++ )
00377 if ( *ptr )
00378 size++;
00379
00380 return size;
00381 }
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391 int short_test(
00392 const void * p1,
00393 const void * p2
00394 )
00395 {
00396 short * s1 = ( short * )p1;
00397 short * s2 = ( short * )p2;
00398
00399 if ( *s1 < *s2 )
00400 return -1;
00401
00402 return ( *s1 > *s2 );
00403 }
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414 static char *
00415 calc_hist( r_afni_s * A, mask_opt_s * M )
00416 {
00417 float * data;
00418 float * ptr;
00419 float bin_size, cum, junk;
00420
00421 long size, new_size = 0;
00422 long count;
00423 int * bins;
00424
00425 int places;
00426
00427
00428 if (( data = (float *)malloc(A->subs[1] * A->nvox * sizeof(float))) == NULL)
00429 {
00430 sprintf( grMessage, "Error: maskcalc_ch_00\n"
00431 "Failed to allocate memory for %d floats.\n",
00432 A->nvox * A->subs[1] );
00433 return grMessage;
00434 }
00435
00436 if ( ( size = mask_all_shorts_to_float( A, 1, 0, data ) ) == 0 )
00437 {
00438 sprintf( grMessage, "Error: 5090\n"
00439 "Masking shorts results in empty array.\n" );
00440 return grMessage;
00441 }
00442
00443 if ( !M->use_min && !M->use_max )
00444 assign_min_max( data, size, &M->min, &M->max );
00445 else if ( !M->use_max )
00446 {
00447 assign_min_max( data, size, &junk, &M->max );
00448
00449 if ( M->min > M->max )
00450 {
00451 sprintf( grMessage, "Error: maskcalc_ch_10\n"
00452 "Min of %f is greater than max of %f\n",
00453 M->min, M->max );
00454 return grMessage;
00455 }
00456 }
00457
00458 junk = ( fabsf( M->max ) > fabsf( M->min ) ) ?
00459 fabsf( M->max ) : fabsf( M->min );
00460
00461 if ( junk == 0 )
00462 places = 2;
00463 else
00464 {
00465 places = 4 - ( int )floor( log10( junk ) );
00466
00467 if ( places > 7 )
00468 places = 7;
00469 else if ( places < 0 )
00470 places = 0;
00471 }
00472
00473 if ( ( bins = (int *)calloc( M->num_bins, sizeof( int ) ) ) == NULL )
00474 {
00475 sprintf( grMessage, "Error: maskcalc_ch_30\n"
00476 "Failed to allocate for %d longs.\n", M->num_bins );
00477 return grMessage;
00478 }
00479
00480 bin_size = ( M->max - M->min ) / M->num_bins;
00481 bin_size += 0.000001 * bin_size;
00482 if ( bin_size == 0.0 )
00483 bin_size = 1.0e-34;
00484
00485 for ( count = 0, ptr = data; count < size; count++, ptr++ )
00486 if ( ( *ptr <= M->max ) && ( *ptr >= M->min ) )
00487 {
00488 bins[ ( int )( ( *ptr - M->min ) / bin_size ) ]++;
00489 new_size++;
00490 }
00491
00492 if ( new_size == 0 )
00493 new_size = 1;
00494
00495
00496 fprintf( M->outfp, "\n range \t #vox \t %% \t cum %%\n");
00497 fprintf( M->outfp, "------------------- \t--------\t------\t-------\n");
00498
00499 cum = 0.0;
00500 for ( count = 0; count < M->num_bins; count++ )
00501 {
00502 cum += 100.0 * bins[ count ] / new_size;
00503
00504 fprintf( M->outfp, "[%8.*f,%8.*f) \t%8d\t%6.3f\t%7.3f\n",
00505 places, M->min + count * bin_size,
00506 places, M->min + (count+1) * bin_size,
00507 bins[ count ],
00508 100.0 * bins[ count ] / new_size,
00509 cum );
00510 }
00511 fputc( '\n', M->outfp );
00512
00513 return NULL;
00514 }
00515
00516
00517
00518
00519
00520
00521
00522
00523 static void
00524 assign_min_max(
00525 float * data,
00526 long size,
00527 float * min,
00528 float * max
00529 )
00530 {
00531 float * ptr = data;
00532 long count;
00533
00534
00535 *min = *data;
00536 *max = *data;
00537
00538 for ( count = 1; count < size; count++, ptr++ )
00539 {
00540 if ( *ptr < *min )
00541 *min = *ptr;
00542
00543 if ( *ptr > *max )
00544 *max = *ptr;
00545 }
00546 }
00547
00548
00549
00550
00551
00552
00553
00554
00555 static char *
00556 calc_stats( r_afni_s * A, mask_opt_s * M )
00557 {
00558 float * data;
00559 float min, max, savemin, savemax;
00560
00561 long size;
00562 int sub;
00563
00564
00565 if ( ( data = ( float * )malloc( A->nvox * sizeof(float) ) )
00566 == NULL )
00567 {
00568 sprintf( grMessage, "Error: 5130\n"
00569 "Failed to allocate memory for %d floats.\n",
00570 A->nvox );
00571 return grMessage;
00572 }
00573
00574 print_stats_header( M->outfp );
00575
00576 for ( sub = 0; sub < A->subs[1]; sub++ )
00577 {
00578 if ( ( size = mask_shorts_to_float( A, data, sub, 0, 1 ) ) == 0 )
00579 {
00580 sprintf( grMessage, "Error: 5140\n"
00581 "Masking shorts results in empty array.\n" );
00582 return grMessage;
00583 }
00584
00585 assign_min_max( data, size, &savemin, &savemax );
00586
00587 if ( ! M->use_min && ! M->use_max )
00588 {
00589 min = savemin;
00590 max = savemax;
00591
00592 do_stats( A, data, size, min, max, sub, M->outfp, NULL, NULL, NULL);
00593 }
00594 else if ( ! M->use_max )
00595 {
00596 min = M->min;
00597 max = savemax;
00598
00599 if ( min <= max )
00600 do_stats( A, data, size, min, max, sub, M->outfp,
00601 NULL, NULL, NULL);
00602 else
00603 print_empty_stats( M->outfp );
00604 }
00605 else
00606 {
00607
00608
00609 min = savemin;
00610 max = M->min;
00611
00612 if ( min <= max )
00613 do_stats( A, data, size, min, max, sub, M->outfp,
00614 NULL, NULL, NULL);
00615 else
00616 print_empty_stats( M->outfp );
00617
00618 min = M->max;
00619 max = savemax;
00620
00621 if ( min <= max )
00622 do_stats( A, data, size, min, max, sub, M->outfp,
00623 NULL, NULL, NULL);
00624 else
00625 print_empty_stats( M->outfp );
00626 }
00627 }
00628
00629 return NULL;
00630 }
00631
00632
00633
00634
00635
00636
00637
00638
00639 static void
00640 print_stats_header ( FILE * fp )
00641 {
00642 fprintf( fp, " Sub\t vol \t%% BRIK\t min \t max "
00643 "\t mean \t SEM \t STD \t 95%% "
00644 "\t thresh # \t thresh %%\n" );
00645 fprintf( fp, "-----\t-------\t------\t------\t------"
00646 "\t------\t------\t------\t-------------"
00647 "\t----------\t---------\n" );
00648 }
00649
00650
00651
00652
00653
00654
00655
00656
00657 static void
00658 print_empty_stats ( FILE * fp )
00659 {
00660 fprintf( fp, " 0 \t 0 \t 0 \t 0 \t 0 "
00661 "\t 0 \t 0 \t 0 \t ( 0, 0 ) "
00662 "\t 0 \t 0 \n" );
00663 }
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673 static void
00674 do_stats(
00675 r_afni_s * A,
00676 float * data,
00677 long size,
00678 float min,
00679 float max,
00680 int sub,
00681 FILE * fp,
00682 long * ret_size,
00683 float * average,
00684 float * ret_var
00685 )
00686 {
00687 double sum_diff2 = 0.0;
00688 double dtmp;
00689
00690 float mean, SEM, STD;
00691 float * ptr = data;
00692 float tmp;
00693 float local_min = max, local_max = min;
00694
00695 long count;
00696 long new_size;
00697
00698
00699
00700
00701 dtmp = 0.0;
00702 new_size = 0;
00703 for ( count = 0, ptr = data; count < size; count++, ptr++ )
00704 if ( ( *ptr >= min ) && ( *ptr <= max ) )
00705 {
00706 new_size++;
00707
00708 dtmp += *ptr;
00709
00710 if ( *ptr < local_min )
00711 local_min = *ptr;
00712
00713 if ( *ptr > local_max )
00714 local_max = *ptr;
00715 }
00716
00717 if ( new_size > 0 )
00718 mean = dtmp / new_size;
00719 else
00720 mean = 0;
00721
00722
00723 sum_diff2 = 0.0;
00724 for ( count = 0, ptr = data; count < size; count++, ptr++ )
00725 if ( ( *ptr >= min ) && ( *ptr <= max ) )
00726 {
00727 tmp = *ptr - mean;
00728 sum_diff2 += tmp * tmp;
00729 }
00730
00731 if ( new_size < 2 )
00732 {
00733 STD = 0;
00734 SEM = 0;
00735 }
00736 else
00737 {
00738 STD = sqrt( sum_diff2 / ( new_size - 1 ) );
00739 SEM = STD / sqrt( new_size );
00740 }
00741
00742 fprintf( fp,
00743 "%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)"
00744 "\t %8ld \t %6.*f\n",
00745 sub, size,
00746 num_places( 100.0*size/A->nvox, 5 ), 100.0*size/A->nvox,
00747 num_places( local_min, 6 ), local_min,
00748 num_places( local_max, 6 ), local_max,
00749 num_places( mean, 6 ), mean,
00750 num_places( SEM, 6 ), SEM,
00751 num_places( STD, 6 ), STD,
00752 num_places( mean-1.96*SEM, 5 ), mean-1.96*SEM,
00753 num_places( mean+1.96*SEM, 5 ), mean+1.96*SEM,
00754 new_size,
00755 num_places( 100.0*new_size/size, 6 ), 100.0*new_size/size );
00756
00757 if ( ret_size != NULL )
00758 *ret_size = new_size;
00759
00760 if ( average != NULL )
00761 *average = mean;
00762
00763 if ( ret_var != NULL )
00764 *ret_var = STD*STD;
00765 }
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775 static int
00776 num_places(
00777 float num,
00778 int size
00779 )
00780 {
00781 float junk;
00782 int places;
00783
00784 junk = fabsf( num );
00785
00786 junk = ( junk == 0 ) ? 1.0 : junk;
00787
00788
00789
00790
00791
00792
00793 places = size - 2 - ( int )floor( log10( junk ) );
00794
00795 if ( places < 0 )
00796 places = 0;
00797 else if ( places >= size )
00798 places = size - 1;
00799
00800 return places;
00801 }
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812 static long
00813 mask_shorts_to_float(
00814 r_afni_s * A,
00815 float * data,
00816 int sub,
00817 int mask_index,
00818 int data_index
00819 )
00820 {
00821 float * fptr;
00822 short * mptr;
00823 short * sptr;
00824
00825 float factor;
00826 long count;
00827
00828
00829 fptr = data;
00830
00831 if ( A->subs[ mask_index ] > 1 )
00832 mptr = A->simage[ mask_index ][ sub ];
00833 else
00834 mptr = A->simage[ mask_index ][ 0 ];
00835
00836 sptr = A->simage[ data_index ][ sub ];
00837
00838
00839
00840
00841
00842
00843 if ( ( factor = A->factor[ data_index ][ sub ] ) == 1 )
00844 {
00845 for( count = 0; count < A->nvox; count++, mptr++, sptr++ )
00846 if ( *mptr )
00847 *fptr++ = *sptr;
00848 }
00849 else
00850 {
00851 for( count = 0; count < A->nvox; count++, mptr++, sptr++ )
00852 if ( *mptr )
00853 *fptr++ = *sptr * factor;
00854 }
00855
00856
00857 return( ( long )( fptr - data ) );
00858 }
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869 static long
00870 mask_shorts_to_short(
00871 r_afni_s * A,
00872 short * data,
00873 int sub,
00874 int mask_index,
00875 int data_index
00876 )
00877 {
00878 short * dptr;
00879 short * mptr;
00880 short * sptr;
00881
00882 long count;
00883
00884
00885 dptr = data;
00886
00887 if ( A->subs[ mask_index ] > 1 )
00888 mptr = A->simage[ mask_index ][ sub ];
00889 else
00890 mptr = A->simage[ mask_index ][ 0 ];
00891
00892 sptr = A->simage[ data_index ][ sub ];
00893
00894
00895
00896
00897
00898
00899 for( count = 0; count < A->nvox; count++, mptr++, sptr++ )
00900 if ( *mptr )
00901 *dptr++ = *sptr;
00902
00903 return( ( long )( dptr - data ) );
00904 }
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915 static long
00916 mask_all_shorts_to_float(
00917 r_afni_s * A,
00918 int data_index,
00919 int mask_index,
00920 float * data
00921 )
00922 {
00923 float * fptr;
00924 short * mptr;
00925 short * sptr;
00926
00927 float factor;
00928 long count;
00929 int sub;
00930
00931
00932 fptr = data;
00933
00934 for ( sub = 0; sub < A->subs[data_index]; sub ++ )
00935 {
00936 if ( A->subs[mask_index] == A->subs[data_index] )
00937 mptr = A->simage[ mask_index ][ sub ];
00938 else
00939 mptr = A->simage[ mask_index ][ 0 ];
00940
00941 sptr = A->simage[ data_index ][ sub ];
00942
00943
00944
00945
00946
00947
00948 if ( ( factor = A->factor[ data_index ][ sub ] ) == 1 )
00949 {
00950 for( count = 0; count < A->nvox; count++, mptr++, sptr++ )
00951 if ( *mptr )
00952 *fptr++ = *sptr;
00953 }
00954 else
00955 {
00956 for( count = 0; count < A->nvox; count++, mptr++, sptr++ )
00957 if ( *mptr )
00958 *fptr++ = *sptr * factor;
00959 }
00960 }
00961
00962 return( ( long )( fptr - data ) );
00963 }
00964
00965
00966
00967
00968
00969
00970
00971
00972 static FILE *
00973 open_file(
00974 char * file,
00975 char * mode
00976 )
00977 {
00978 return fopen( file, mode );
00979 }
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996 static int
00997 file_exists(
00998 char * filename,
00999 char * suffix
01000 )
01001 {
01002 struct stat buf;
01003 char file[ R_FILE_L + 6 ] = "";
01004 char * filep = file;
01005
01006 if ( suffix == NULL )
01007 filep = filename;
01008 else if ( ! strcmp( suffix, filename+strlen(filename)-strlen(suffix) ) )
01009 strcpy( file, filename );
01010 else if ( filename[ strlen( filename ) - 1 ] == '.' )
01011 sprintf( file, "%s%s", filename, suffix );
01012 else
01013 sprintf( file, "%s.%s", filename, suffix );
01014
01015
01016 return ( stat( filep, &buf ) == 0 );
01017 }
01018
01019
01020
01021
01022
01023
01024
01025
01026 static char *
01027 fill_afni_struct( r_afni_s * A )
01028 {
01029 u_short mus;
01030 int sub, brick;
01031
01032 for ( brick = 0; brick < A->num_dsets; brick++ )
01033 {
01034 A->subs[ brick ] = DSET_NVALS( A->dset[ brick ] );
01035
01036 if ( A->max_subs && ( A->subs[brick] > A->max_subs ) )
01037 {
01038 sprintf( grMessage, "------------------------------------\n"
01039 "Error: maskcalc_fas_00\n"
01040 "Brick #%d contains %d sub-bricks.\n"
01041 "The limit is %d.\n"
01042 "------------------------------------",
01043 brick, A->subs[brick], A->max_subs );
01044 return grMessage;
01045 }
01046
01047 if ( A->subs_must_equal && ( A->subs[brick] != A->subs[0] ) )
01048 {
01049 sprintf( grMessage, "------------------------------------\n"
01050 "Error: maskcalc_fas_02\n"
01051 "Brick #%d contains %d sub-bricks.\n"
01052 "Brick #%d contains %d sub-bricks.\n"
01053 "We are requiring them to be equal.\n"
01054 "------------------------------------",
01055 0, A->subs[0],
01056 brick, A->subs[brick] );
01057 return grMessage;
01058 }
01059
01060 if ( ( A->simage[brick] =
01061 (short **)malloc( A->subs[brick]*sizeof(short *)) ) == NULL )
01062 {
01063 return "-------------------------\n"
01064 "Error: maskcalc_fas_05\n"
01065 "memory allocation failure\n"
01066 "-------------------------";
01067 }
01068 if ( ( A->factor[brick] =
01069 (float *)malloc( A->subs[brick]*sizeof(float)) ) == NULL)
01070 {
01071 return "-------------------------\n"
01072 "Error: maskcalc_fas_10\n"
01073 "memory allocation failure\n"
01074 "-------------------------";
01075 }
01076
01077 for ( sub = 0; sub < A->subs[brick]; sub++ )
01078 {
01079 A->simage[brick][sub] = (short *)DSET_ARRAY(A->dset[brick],sub);
01080 A->factor[brick][sub] = DSET_BRICK_FACTOR(A->dset[brick],sub);
01081 if ( A->factor[brick][sub] == 0.0 )
01082 A->factor[brick][sub] = 1.0;
01083 }
01084
01085 if ( brick == 0 )
01086 {
01087 A->nx = A->dset[brick]->daxes->nxx;
01088 A->ny = A->dset[brick]->daxes->nyy;
01089 A->nz = A->dset[brick]->daxes->nzz;
01090 A->nvox = A->nx * A->ny * A->nz;
01091 }
01092 else if ( ( A->nx != A->dset[brick]->daxes->nxx ) ||
01093 ( A->ny != A->dset[brick]->daxes->nyy ) ||
01094 ( A->nz != A->dset[brick]->daxes->nzz ) )
01095 {
01096 sprintf( grMessage,
01097 "--------------------------------\n"
01098 "Error : maskcalc_fas_20\n"
01099 "Unaligned dimensions.\n"
01100 "(%d,%d,%d) != (%d,%d,%d)\n"
01101 "--------------------------------",
01102 A->dset[brick]->daxes->nxx, A->dset[brick]->daxes->nyy,
01103 A->dset[brick]->daxes->nzz, A->nx, A->ny, A->nz );
01104 return grMessage;
01105 }
01106 }
01107
01108 if ( A->want_floats && ! assign_afni_floats( A ) )
01109 return "-----------------------------\n"
01110 "Error: maskcalc_fas_30\n"
01111 "Failed to create afni fimage.\n"
01112 "-----------------------------";
01113
01114 mus = 0;
01115 A->max_u_short = 0;
01116 for ( brick = 1; brick < A->num_dsets; brick++ )
01117 for ( sub = 0; sub < A->subs[brick]; sub++ )
01118 {
01119 mus = r_get_max_u_short( (u_short *)A->simage[brick][sub], A->nvox);
01120 if ( mus > A->max_u_short )
01121 A->max_u_short = mus;
01122 }
01123
01124 return NULL;
01125 }
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135 static int
01136 assign_afni_floats( r_afni_s * A )
01137 {
01138 short * sptr;
01139 float * fptr;
01140 float factor = A->factor[1][0];
01141 int count;
01142
01143
01144 if ( ( A->fimage[1] = (float *)malloc( A->nvox * sizeof(float))) == NULL )
01145 return 0;
01146
01147 for ( count = 0, fptr = A->fimage[1], sptr = A->simage[1][0];
01148 count < A->nvox;
01149 count++, fptr++, sptr++ )
01150 *fptr = factor * *sptr;
01151
01152 return 1;
01153 }
01154
01155
01156
01157
01158
01159
01160
01161
01162 static u_short
01163 r_get_max_u_short( u_short * S, int size )
01164 {
01165 u_short * usptr, max = *S;
01166 int c = 0;
01167
01168 for ( c = 0, usptr = S; c < size; c++, usptr++ )
01169 {
01170 if ( *usptr > max )
01171 max = *usptr;
01172 }
01173
01174 return max;
01175 }
01176
01177
