00001
00002
00003
00004
00005
00006
00007 #include "afni.h"
00008
00009 #ifndef ALLOW_PLUGINS
00010 # error "Plugins not properly set up -- see machdep.h"
00011 #endif
00012
00013
00014
00015
00016
00017 char * MASKAVE_main( PLUGIN_interface * ) ;
00018
00019 static char helpstring[] =
00020 " Purpose: Compute statistics over a mask-selected ROI.\n"
00021 "\n"
00022 " Source: Dataset = data to be averaged\n"
00023 " Sub-brick = which one to use\n"
00024 " [if -1 is input, will do all sub-bricks]\n\n"
00025 " Mask: Dataset = masking dataset\n"
00026 " Sub-brick = which one to use\n\n"
00027 " Range: Bottom = min value from mask dataset to use\n"
00028 " Top = max value from mask dataset to use\n"
00029 " [if Bottom > Top, then all nonzero mask voxels will be used; ]\n"
00030 " [if Bottom <= Top, then only nonzero mask voxels in this range]\n"
00031 " [ will be used in computing the statistics. ]\n\n"
00032 " 1D Save: Name = If all input sub-bricks are used (i.e., setting\n"
00033 " 'Source Sub-brick' = -1 above), then this option\n"
00034 " lets you save the resulting average time series\n"
00035 " into the interal list of timeseries available via\n"
00036 " the 'Choose Timeseries', 'Pick Ideal', ... buttons.\n"
00037 " To Disk? = If 'YES' is chosen, then will also write the\n"
00038 " timeseries to disk in the *.1D format.\n\n"
00039 " Author -- RW Cox"
00040 ;
00041
00042 #define NUM_yesno_list 2
00043 static char *yesno_list[] = { "YES" , "NO" } ;
00044
00045
00046
00047
00048
00049
00050
00051 DEFINE_PLUGIN_PROTOTYPE
00052
00053 PLUGIN_interface * PLUGIN_init( int ncall )
00054 {
00055 PLUGIN_interface * plint ;
00056
00057 if( ncall > 0 ) return NULL ;
00058
00059
00060
00061 plint = PLUTO_new_interface( "ROI Average" , "Average Dataset over ROI" , helpstring ,
00062 PLUGIN_CALL_VIA_MENU , MASKAVE_main ) ;
00063
00064 PLUTO_add_hint( plint , "Average Dataset over ROI" ) ;
00065
00066 PLUTO_set_sequence( plint , "A:afniinfo:dset" ) ;
00067
00068
00069
00070 PLUTO_add_option( plint , "Source" , "Source" , TRUE ) ;
00071 PLUTO_add_dataset(plint , "Dataset" ,
00072 ANAT_ALL_MASK , FUNC_ALL_MASK ,
00073 DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
00074
00075 PLUTO_add_number( plint , "Sub-brick" , -1,9999,0 , 0,1 ) ;
00076
00077
00078
00079 PLUTO_add_option( plint , "Mask" , "Mask" , TRUE ) ;
00080 PLUTO_add_dataset( plint , "Dataset" ,
00081 ANAT_ALL_MASK , FUNC_ALL_MASK ,
00082 DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
00083 PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
00084
00085
00086
00087 PLUTO_add_option( plint , "Range" , "Range" , FALSE ) ;
00088 PLUTO_add_number( plint , "Bottom" , -99999,99999, 1, 0,1 ) ;
00089 PLUTO_add_number( plint , "Top" , -99999,99999,-1, 0,1 ) ;
00090
00091
00092
00093 PLUTO_add_option( plint , "1D Save" , "1D Save" , FALSE ) ;
00094 PLUTO_add_string( plint , "Name" , 0,NULL , 12 ) ;
00095 PLUTO_add_string( plint , "To Disk?" , NUM_yesno_list , yesno_list , 1 ) ;
00096
00097 return plint ;
00098 }
00099
00100
00101
00102
00103
00104 char * MASKAVE_main( PLUGIN_interface * plint )
00105 {
00106 MCW_idcode * idc ;
00107 THD_3dim_dataset * input_dset , * mask_dset ;
00108 int iv , mcount , nvox , ii , sigmait , nvals , doall , ivbot,ivtop ;
00109 float mask_bot=666.0 , mask_top=-666.0 ;
00110 double sum , sigma ;
00111 float * sumar=NULL , * sigmar=NULL ;
00112 char * tag , * str , buf[64] , abuf[32],sbuf[32] ;
00113 byte * mmm ;
00114
00115 char * cname=NULL ;
00116 int cdisk=0 ;
00117 int miv=0 ;
00118
00119
00120
00121
00122 if( plint == NULL )
00123 return "*************************\n"
00124 "MASKAVE_main: NULL input\n"
00125 "*************************" ;
00126
00127
00128
00129 PLUTO_next_option(plint) ;
00130 idc = PLUTO_get_idcode(plint) ;
00131 input_dset = PLUTO_find_dset(idc) ;
00132 if( input_dset == NULL )
00133 return "********************************\n"
00134 "MASKAVE_main: bad input dataset\n"
00135 "********************************" ;
00136
00137 iv = (int) PLUTO_get_number(plint) ;
00138 if( iv >= DSET_NVALS(input_dset) )
00139 return "**********************************\n"
00140 "MASKAVE_main: bad input sub-brick\n"
00141 "**********************************" ;
00142 doall = (iv < 0) ;
00143 if( doall ){
00144 nvals = DSET_NVALS(input_dset) ;
00145 ivbot = 0 ; ivtop = nvals-1 ;
00146 } else {
00147 ivbot = ivtop = iv ;
00148 }
00149 DSET_load(input_dset) ;
00150 if( DSET_ARRAY(input_dset,ivbot) == NULL )
00151 return "*********************************\n"
00152 "MASKAVE_main: can't load dataset\n"
00153 "*********************************" ;
00154 nvox = DSET_NVOX(input_dset) ;
00155
00156
00157
00158 PLUTO_next_option(plint) ;
00159 idc = PLUTO_get_idcode(plint) ;
00160 mask_dset = PLUTO_find_dset(idc) ;
00161
00162 if( mask_dset == NULL )
00163 return "*******************************\n"
00164 "MASKAVE_main: bad mask dataset\n"
00165 "*******************************" ;
00166
00167 if( DSET_NVOX(mask_dset) != nvox )
00168 return "*************************************************************\n"
00169 "MASKAVE_main: mask input dataset doesn't match source dataset\n"
00170 "*************************************************************" ;
00171
00172 miv = (int) PLUTO_get_number(plint) ;
00173 if( miv >= DSET_NVALS(mask_dset) )
00174 return "*****************************************************\n"
00175 "MASKAVE_main: mask dataset sub-brick index is too big\n"
00176 "*****************************************************" ;
00177
00178 DSET_load(mask_dset) ;
00179 if( DSET_ARRAY(mask_dset,0) == NULL )
00180 return "**************************************\n"
00181 "MASKAVE_main: can't load mask dataset\n"
00182 "**************************************" ;
00183
00184
00185
00186 while( (tag=PLUTO_get_optiontag(plint)) != NULL ){
00187
00188 if( strcmp(tag,"Range") == 0 ){
00189 mask_bot = PLUTO_get_number(plint) ;
00190 mask_top = PLUTO_get_number(plint) ;
00191 continue ;
00192 }
00193
00194 if( strcmp(tag,"1D Save") == 0 ){
00195 char * yn ;
00196 cname = PLUTO_get_string(plint) ;
00197 yn = PLUTO_get_string(plint) ;
00198 cdisk = (strcmp(yn,yesno_list[0]) == 0) ;
00199 continue ;
00200 }
00201
00202 }
00203
00204
00205
00206
00207
00208
00209 mmm = (byte *) malloc( sizeof(byte) * nvox ) ;
00210 if( mmm == NULL )
00211 return "*** Can't malloc workspace! ***" ;
00212
00213
00214
00215 switch( DSET_BRICK_TYPE(mask_dset,miv) ){
00216 default:
00217 free(mmm) ;
00218 return "*** Can't use mask dataset -- illegal data type! ***" ;
00219
00220 case MRI_short:{
00221 short mbot , mtop ;
00222 short * mar = (short *) DSET_ARRAY(mask_dset,miv) ;
00223 float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ;
00224 if( mfac == 0.0 ) mfac = 1.0 ;
00225 if( mask_bot <= mask_top ){
00226 mbot = SHORTIZE(mask_bot/mfac) ;
00227 mtop = SHORTIZE(mask_top/mfac) ;
00228 } else {
00229 mbot = (short) -MRI_TYPE_maxval[MRI_short] ;
00230 mtop = (short) MRI_TYPE_maxval[MRI_short] ;
00231 }
00232 for( mcount=0,ii=0 ; ii < nvox ; ii++ )
00233 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 ){ mmm[ii] = 1 ; mcount++ ; }
00234 else { mmm[ii] = 0 ; }
00235 }
00236 break ;
00237
00238 case MRI_byte:{
00239 byte mbot , mtop ;
00240 byte * mar = (byte *) DSET_ARRAY(mask_dset,miv) ;
00241 float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ;
00242 if( mfac == 0.0 ) mfac = 1.0 ;
00243 if( mask_bot <= mask_top ){
00244 mbot = BYTEIZE(mask_bot/mfac) ;
00245 mtop = BYTEIZE(mask_top/mfac) ;
00246 if( mtop == 0 ){
00247 free(mmm) ;
00248 return "*** Illegal mask range for mask dataset of bytes. ***" ;
00249 }
00250 } else {
00251 mbot = 0 ;
00252 mtop = (byte) MRI_TYPE_maxval[MRI_short] ;
00253 }
00254 for( mcount=0,ii=0 ; ii < nvox ; ii++ )
00255 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 ){ mmm[ii] = 1 ; mcount++ ; }
00256 else { mmm[ii] = 0 ; }
00257 }
00258 break ;
00259
00260 case MRI_float:{
00261 float mbot , mtop ;
00262 float * mar = (float *) DSET_ARRAY(mask_dset,miv) ;
00263 float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ;
00264 if( mfac == 0.0 ) mfac = 1.0 ;
00265 if( mask_bot <= mask_top ){
00266 mbot = (float) (mask_bot/mfac) ;
00267 mtop = (float) (mask_top/mfac) ;
00268 } else {
00269 mbot = -WAY_BIG ;
00270 mtop = WAY_BIG ;
00271 }
00272 for( mcount=0,ii=0 ; ii < nvox ; ii++ )
00273 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 ){ mmm[ii] = 1 ; mcount++ ; }
00274 else { mmm[ii] = 0 ; }
00275 }
00276 break ;
00277 }
00278
00279 if( mcount == 0 ){
00280 free(mmm) ;
00281 return "*** No voxels survive the masking operations! ***" ;
00282 }
00283 sigmait = (mcount > 1) ;
00284
00285
00286
00287 if( doall ){
00288 sumar = (float *) malloc( sizeof(float) * nvals ) ;
00289 sigmar = (float *) malloc( sizeof(float) * nvals ) ;
00290 }
00291
00292 for( iv=ivbot ; iv <= ivtop ; iv++ ){
00293 sum = sigma = 0.0 ;
00294 switch( DSET_BRICK_TYPE(input_dset,iv) ){
00295
00296 default:
00297 free(mmm) ; if( doall ){ free(sumar) ; free(sigmar) ; }
00298 return "*** Can't use source dataset -- illegal data type! ***" ;
00299
00300 case MRI_short:{
00301 short * bar = (short *) DSET_ARRAY(input_dset,iv) ;
00302 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
00303 if( mfac == 0.0 ) mfac = 1.0 ;
00304
00305 for( ii=0 ; ii < nvox ; ii++ ) if( mmm[ii] ) sum += bar[ii] ;
00306 sum = sum / mcount ;
00307
00308 if( sigmait ){
00309 for( ii=0 ; ii < nvox ; ii++ )
00310 if( mmm[ii] ) sigma += SQR(bar[ii]-sum) ;
00311 sigma = mfac * sqrt( sigma/(mcount-1) ) ;
00312 }
00313 sum = mfac * sum ;
00314 }
00315 break ;
00316
00317 case MRI_byte:{
00318 byte * bar = (byte *) DSET_ARRAY(input_dset,iv) ;
00319 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
00320 if( mfac == 0.0 ) mfac = 1.0 ;
00321
00322 for( ii=0 ; ii < nvox ; ii++ ) if( mmm[ii] ) sum += bar[ii] ;
00323 sum = sum / mcount ;
00324
00325 if( sigmait ){
00326 for( ii=0 ; ii < nvox ; ii++ )
00327 if( mmm[ii] ) sigma += SQR(bar[ii]-sum) ;
00328 sigma = mfac * sqrt( sigma/(mcount-1) ) ;
00329 }
00330 sum = mfac * sum ;
00331 }
00332 break ;
00333
00334 case MRI_float:{
00335 float * bar = (float *) DSET_ARRAY(input_dset,iv) ;
00336 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
00337 if( mfac == 0.0 ) mfac = 1.0 ;
00338
00339 for( ii=0 ; ii < nvox ; ii++ ) if( mmm[ii] ) sum += bar[ii] ;
00340 sum = sum / mcount ;
00341
00342 if( sigmait ){
00343 for( ii=0 ; ii < nvox ; ii++ )
00344 if( mmm[ii] ) sigma += SQR(bar[ii]-sum) ;
00345 sigma = mfac * sqrt( sigma/(mcount-1) ) ;
00346 }
00347 sum = mfac * sum ;
00348 }
00349 break ;
00350 }
00351
00352 if( doall ){ sumar[iv] = sum ; sigmar[iv] = sigma ; }
00353 }
00354
00355 free(mmm) ;
00356
00357
00358
00359 if( doall ){
00360 str = (char *) malloc( 1024 + 64*nvals ) ;
00361 sprintf(str," ****** ROI statistics ****** \n"
00362 " Source = %s [all sub-bricks] \n"
00363 " Mask = %s [%s]" ,
00364 DSET_FILECODE(input_dset) ,
00365 DSET_FILECODE(mask_dset) , DSET_BRICK_LABEL(mask_dset,miv) ) ;
00366 if( mask_bot <= mask_top ){
00367 sprintf(buf," [range %g .. %g]" , mask_bot , mask_top ) ;
00368 strcat(str,buf) ;
00369 }
00370 strcat(str," \n") ;
00371 sprintf(buf," Count = %d voxels\n",mcount) ; strcat(str,buf) ;
00372 for( iv=0 ; iv < nvals ; iv++ ){
00373 AV_fval_to_char( sumar[iv] , abuf ) ;
00374 AV_fval_to_char( sigmar[iv] , sbuf ) ;
00375 sprintf(buf," Average = %9.9s Sigma = %9.9s [%s] \n",
00376 abuf,sbuf , DSET_BRICK_LABEL(input_dset,iv) ) ;
00377 strcat(str,buf) ;
00378 }
00379 PLUTO_popup_textwin( plint , str ) ;
00380
00381
00382
00383 if( cname != NULL && cname[0] != '\0' ){
00384 MRI_IMAGE * qim = mri_new_vol_empty( nvals,1,1 , MRI_float ) ;
00385 mri_fix_data_pointer( sumar , qim ) ;
00386 PLUTO_register_timeseries( cname , qim ) ;
00387
00388 if( cdisk ){
00389 if( PLUTO_prefix_ok(cname) ){
00390 char * cn ;
00391 if( strstr(cname,".1D") == NULL ){
00392 cn = malloc(strlen(cname)+8) ;
00393 strcpy(cn,cname) ; strcat(cn,".1D") ;
00394 } else {
00395 cn = cname ;
00396 }
00397 mri_write_1D( cn , qim ) ;
00398 if( cn != cname ) free(cn) ;
00399 } else {
00400 PLUTO_popup_transient(plint," \n"
00401 "** Illegal filename **\n"
00402 "** in 'To Disk?' !! **\n" ) ;
00403 }
00404 }
00405
00406 mri_fix_data_pointer( NULL , qim ) ; mri_free(qim) ;
00407 }
00408
00409 free(str) ; free(sumar) ; free(sigmar) ;
00410
00411 } else if( mask_bot <= mask_top ){
00412 str = (char *) malloc( 1024 ) ;
00413 sprintf( str , " *** ROI Statistics *** \n"
00414 " Source = %s [%s] \n"
00415 " Mask = %s [%s] [range %g .. %g] \n"
00416 " Count = %d voxels \n"
00417 " Average = %g \n"
00418 " Sigma = %g " ,
00419 DSET_FILECODE(input_dset) , DSET_BRICK_LABEL(input_dset,ivbot) ,
00420 DSET_FILECODE(mask_dset) , DSET_BRICK_LABEL(mask_dset,miv) ,
00421 mask_bot , mask_top , mcount , sum , sigma ) ;
00422 PLUTO_popup_message(plint,str) ;
00423 free(str) ;
00424
00425 } else {
00426 str = (char *) malloc( 1024 ) ;
00427 sprintf( str , " *** ROI Statistics *** \n"
00428 " Source = %s [%s] \n"
00429 " Mask = %s [%s] \n"
00430 " Count = %d voxels \n"
00431 " Average = %g \n"
00432 " Sigma = %g " ,
00433 DSET_FILECODE(input_dset) , DSET_BRICK_LABEL(input_dset,ivbot) ,
00434 DSET_FILECODE(mask_dset) , DSET_BRICK_LABEL(mask_dset,miv) ,
00435 mcount , sum , sigma ) ;
00436 PLUTO_popup_message(plint,str) ;
00437 free(str) ;
00438 }
00439
00440 return NULL ;
00441 }
