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
00018
00019 static char * HISTO_main( PLUGIN_interface * ) ;
00020
00021 static PLUGIN_interface * CORREL_init(void) ;
00022
00023 static char helpstring[] =
00024 " Purpose: Plot a histogram of data from a dataset brick.\n"
00025 "\n"
00026 " Source: Dataset = data to be processed\n"
00027 " Sub-brick = which one to use\n\n"
00028 " Values: Bottom = minimum value from dataset to include\n"
00029 " Top = maximum value from dataset to include\n\n"
00030 " Bins: Number = number of bins to use\n"
00031 " Max Count = maximum count per bin\n\n"
00032 " Smooth = +/- bin range to smooth histogram\n"
00033 " Mask: Dataset = masking dataset\n"
00034 " Sub-brick = which one to use\n\n"
00035 " Range: Bottom = min value from mask dataset to use\n"
00036 " Top = max value from mask dataset to use\n"
00037 " [if Bottom > Top, then all nonzero mask voxels will be used; ]\n"
00038 " [if Bottom <= Top, then only nonzero mask voxels in this range]\n"
00039 " [ will be used in computing the statistics. ]\n"
00040 " Aboot: If activated, then only voxels within a distance of Radius mm\n"
00041 " of the current crosshairs will be used in the histogram.\n"
00042 " Output: Name of the ascii file to which histogram values are written.\n"
00043 "\n"
00044 " Author -- RW Cox - 30 September 1999\n"
00045 " Output feature added by V Roopchansingh, Feb 2002\n"
00046 ;
00047
00048
00049
00050
00051
00052
00053 DEFINE_PLUGIN_PROTOTYPE
00054
00055 PLUGIN_interface * PLUGIN_init( int ncall )
00056 {
00057 PLUGIN_interface * plint ;
00058
00059 #if 0
00060 if( ncall > 0 ) return NULL ;
00061 #else
00062 if( ncall == 1 ) return CORREL_init() ;
00063 if( ncall > 1 ) return NULL ;
00064 #endif
00065
00066
00067
00068 plint = PLUTO_new_interface( "Histogram" ,
00069 "Histogram of Dataset Brick" ,
00070 helpstring ,
00071 PLUGIN_CALL_VIA_MENU , HISTO_main ) ;
00072
00073 PLUTO_add_hint( plint , "Histogram of Dataset Brick" ) ;
00074
00075 PLUTO_set_sequence( plint , "A:afniinfo:dsethistog" ) ;
00076
00077 PLUTO_set_runlabels( plint , "Plot+Keep" , "Plot+Close" ) ;
00078
00079
00080
00081 PLUTO_add_option( plint , "Source" , "Source" , TRUE ) ;
00082 PLUTO_add_dataset(plint , "Dataset" ,
00083 ANAT_ALL_MASK , FUNC_ALL_MASK ,
00084 DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
00085
00086 PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
00087
00088
00089
00090 PLUTO_add_option( plint , "Values" , "Values" , FALSE ) ;
00091 PLUTO_add_number( plint , "Bottom" , -99999,99999, 0, 1,1 ) ;
00092 PLUTO_add_number( plint , "Top" , -99999,99999, 0, -1,1 ) ;
00093
00094
00095
00096 PLUTO_add_option( plint , "Bins" , "Bins" , FALSE ) ;
00097 PLUTO_add_number( plint , "Number" , 10,1000,0, 100,1 ) ;
00098 PLUTO_add_number( plint , "Max Count" , 0,999999999,0 , 0,1 ) ;
00099 PLUTO_add_number( plint , "Smooth" , 0,99,0 , 0,1 ) ;
00100
00101
00102
00103 PLUTO_add_option( plint , "Mask" , "Mask" , FALSE ) ;
00104 PLUTO_add_dataset( plint , "Dataset" ,
00105 ANAT_ALL_MASK , FUNC_ALL_MASK ,
00106 DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
00107 PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
00108
00109
00110
00111 PLUTO_add_option( plint , "Range" , "Range" , FALSE ) ;
00112 PLUTO_add_number( plint , "Bottom" , -99999,99999, 0, 1,1 ) ;
00113 PLUTO_add_number( plint , "Top" , -99999,99999, 0, -1,1 ) ;
00114
00115
00116
00117 PLUTO_add_option( plint , "Aboot" , "Aboot" , FALSE ) ;
00118 PLUTO_add_number( plint , "Radius" , 2,100,0,10,1 ) ;
00119
00120
00121
00122 PLUTO_add_option( plint , "Output", "Output", FALSE ) ;
00123 PLUTO_add_string( plint , "Filename", 0 , NULL , 20 ) ;
00124
00125 return plint ;
00126 }
00127
00128
00129
00130
00131
00132 static char * HISTO_main( PLUGIN_interface * plint )
00133 {
00134 MCW_idcode * idc ;
00135 THD_3dim_dataset * input_dset , * mask_dset=NULL ;
00136 int iv , mcount , nvox , ii,jj , nbin=-1 , do_mval=0,mval ;
00137 float mask_bot=666.0 , mask_top=-666.0 , hbot,htop ;
00138 float val_bot=666.0 , val_top=-666.0 ;
00139 char * tag , * str , buf[THD_MAX_NAME+16] ;
00140 byte * mmm ;
00141 MRI_IMAGE * flim ;
00142 float * flar ;
00143 int * hbin ;
00144 int smooth=0 ;
00145
00146 int miv=0 ;
00147
00148 int maxcount=0 ;
00149 float hrad=0.0 ;
00150
00151 char * histout=NULL ;
00152 FILE * HISTOUT=NULL ;
00153 int writehist=0 ;
00154 float dx ;
00155
00156
00157
00158
00159 if( plint == NULL )
00160 return "***********************\n"
00161 "HISTO_main: NULL input\n"
00162 "***********************" ;
00163
00164
00165
00166 PLUTO_next_option(plint) ;
00167 idc = PLUTO_get_idcode(plint) ;
00168 input_dset = PLUTO_find_dset(idc) ;
00169 if( input_dset == NULL )
00170 return "******************************\n"
00171 "HISTO_main: bad input dataset\n"
00172 "******************************" ;
00173
00174 iv = (int) PLUTO_get_number(plint) ;
00175 if( iv >= DSET_NVALS(input_dset) || iv < 0 )
00176 return "********************************\n"
00177 "HISTO_main: bad input sub-brick\n"
00178 "********************************" ;
00179
00180 DSET_load(input_dset) ;
00181 if( DSET_ARRAY(input_dset,iv) == NULL )
00182 return "*******************************\n"
00183 "HISTO_main: can't load dataset\n"
00184 "*******************************" ;
00185 nvox = DSET_NVOX(input_dset) ;
00186
00187
00188
00189 while( (tag=PLUTO_get_optiontag(plint)) != NULL ){
00190
00191
00192
00193 if( strcmp(tag,"Values") == 0 ){
00194 val_bot = PLUTO_get_number(plint) ;
00195 val_top = PLUTO_get_number(plint) ;
00196 do_mval = (val_bot < val_top) ;
00197 continue ;
00198 }
00199
00200
00201
00202 if( strcmp(tag,"Bins") == 0 ){
00203 nbin = PLUTO_get_number(plint) ;
00204 maxcount = PLUTO_get_number(plint) ;
00205 smooth = PLUTO_get_number(plint) ;
00206 continue ;
00207 }
00208
00209
00210
00211 if( strcmp(tag,"Range") == 0 ){
00212 if( mask_dset == NULL )
00213 return "*****************************************\n"
00214 "HISTO_main: Can't use Range without Mask\n"
00215 "*****************************************" ;
00216
00217 mask_bot = PLUTO_get_number(plint) ;
00218 mask_top = PLUTO_get_number(plint) ;
00219 continue ;
00220 }
00221
00222
00223
00224 if( strcmp(tag,"Mask") == 0 ){
00225
00226 idc = PLUTO_get_idcode(plint) ;
00227 mask_dset = PLUTO_find_dset(idc) ;
00228
00229 if( mask_dset == NULL )
00230 return "*****************************\n"
00231 "HISTO_main: bad mask dataset\n"
00232 "*****************************" ;
00233
00234 if( DSET_NVOX(mask_dset) != nvox )
00235 return "***********************************************************\n"
00236 "HISTO_main: mask input dataset doesn't match source dataset\n"
00237 "***********************************************************" ;
00238
00239 miv = (int) PLUTO_get_number(plint) ;
00240 if( miv >= DSET_NVALS(mask_dset) || miv < 0 )
00241 return "***************************************************\n"
00242 "HISTO_main: mask dataset sub-brick index is illegal\n"
00243 "***************************************************" ;
00244
00245 DSET_load(mask_dset) ;
00246 if( DSET_ARRAY(mask_dset,miv) == NULL )
00247 return "************************************\n"
00248 "HISTO_main: can't load mask dataset\n"
00249 "************************************" ;
00250 continue ;
00251 }
00252
00253
00254
00255 if( strcmp(tag,"Aboot") == 0 ){
00256 hrad = PLUTO_get_number(plint) ;
00257 continue ;
00258 }
00259
00260
00261
00262 if( strcmp(tag,"Output") == 0 ){
00263 histout = PLUTO_get_string(plint) ;
00264 writehist = 1 ;
00265 continue ;
00266 }
00267 }
00268
00269
00270
00271
00272
00273
00274 if( mask_dset == NULL ){
00275 mmm = (byte *) malloc( sizeof(byte) * nvox ) ;
00276 if( mmm == NULL )
00277 return " \n*** Can't malloc workspace! ***\n" ;
00278 memset( mmm , 1, nvox ) ; mcount = nvox ;
00279 } else {
00280
00281 mmm = THD_makemask( mask_dset , miv , mask_bot , mask_top ) ;
00282 if( mmm == NULL )
00283 return " \n*** Can't make mask for some reason! ***\n" ;
00284 mcount = THD_countmask( nvox , mmm ) ;
00285
00286 if( !EQUIV_DSETS(mask_dset,input_dset) ) DSET_unload(mask_dset) ;
00287 if( mcount < 3 ){
00288 free(mmm) ;
00289 return " \n*** Less than 3 voxels survive the mask! ***\n" ;
00290 }
00291 sprintf(buf," \n"
00292 " %d voxels in the mask\n"
00293 " out of %d dataset voxels\n ",mcount,nvox) ;
00294 PLUTO_popup_transient(plint,buf) ;
00295 }
00296
00297
00298
00299 if( hrad > 0.0 ){
00300 MCW_cluster * cl ;
00301 short *di,*dj,*dk ;
00302 int nd , xx,yy,zz , dd , nx,ny,nz,nxy, nx1,ny1,nz1 , ip,jp,kp ;
00303
00304 cl = MCW_build_mask( 0,0,0 , fabs(DSET_DX(input_dset)) ,
00305 fabs(DSET_DY(input_dset)) ,
00306 fabs(DSET_DZ(input_dset)) , hrad ) ;
00307
00308 if( cl == NULL || cl->num_pt < 6 ){
00309 KILL_CLUSTER(cl);
00310 PLUTO_popup_transient(plint, " \n"
00311 " Aboot Radius too small\n"
00312 " for this dataset!\n" ) ;
00313 } else {
00314 ADDTO_CLUSTER(cl,0,0,0,0) ;
00315 di = cl->i ; dj = cl->j ; dk = cl->k ; nd = cl->num_pt ;
00316 nx = DSET_NX(input_dset) ; nx1 = nx-1 ;
00317 ny = DSET_NY(input_dset) ; ny1 = ny-1 ; nxy = nx*ny ;
00318 nz = DSET_NZ(input_dset) ; nz1 = nz-1 ;
00319 xx = plint->im3d->vinfo->i1 ;
00320 yy = plint->im3d->vinfo->j2 ;
00321 zz = plint->im3d->vinfo->k3 ;
00322 for( dd=0 ; dd < nd ; dd++ ){
00323 ip = xx+di[dd] ; if( ip < 0 || ip > nx1 ) continue ;
00324 jp = yy+dj[dd] ; if( jp < 0 || jp > ny1 ) continue ;
00325 kp = zz+dk[dd] ; if( kp < 0 || kp > nz1 ) continue ;
00326 mmm[ip+jp*nx+kp*nxy]++ ;
00327 }
00328 KILL_CLUSTER(cl) ;
00329 for( dd=0 ; dd < nvox ; dd++ ) if( mmm[dd] == 1 ) mmm[dd] = 0 ;
00330 mcount = THD_countmask( nvox , mmm ) ;
00331
00332 if( mcount < 3 ){
00333 free(mmm) ;
00334 return " \n*** Less than 3 voxels survive the mask+radius! ***\n" ;
00335 }
00336 sprintf(buf," \n"
00337 " %d voxels in the mask+radius\n"
00338 " out of %d dataset voxels\n ",mcount,nvox) ;
00339 PLUTO_popup_transient(plint,buf) ;
00340 }
00341 }
00342
00343
00344
00345 if ( writehist )
00346 {
00347 static char hbuf[1024] ;
00348 if ( ( histout == NULL ) || ( strlen (histout) == 0 ) ){
00349 sprintf( hbuf , "%s.histog" , DSET_PREFIX(input_dset) ) ;
00350 } else {
00351 strcpy( hbuf , histout ) ;
00352 if( strstr(hbuf,".hist") == NULL ) strcat( hbuf , ".histog" ) ;
00353 }
00354 histout = hbuf ;
00355
00356 if (THD_is_file(histout))
00357 {
00358 free(mmm) ;
00359
00360 return "*******************************\n"
00361 "Outfile exists, won't overwrite\n"
00362 "*******************************\n" ;
00363 }
00364 else {
00365 HISTOUT = fopen (histout, "w") ;
00366 if( HISTOUT == NULL ){
00367 free(mmm) ;
00368 return "**********************************\n"
00369 "Can't open Outfile for some reason\n"
00370 "**********************************\n" ;
00371 }
00372 }
00373 }
00374
00375
00376
00377 flim = mri_new( mcount , 1 , MRI_float ) ;
00378 flar = MRI_FLOAT_PTR(flim) ;
00379
00380
00381
00382 switch( DSET_BRICK_TYPE(input_dset,iv) ){
00383 default:
00384 free(mmm) ; mri_free(flim) ;
00385 return "*** Can't use source dataset -- illegal data type! ***" ;
00386
00387 case MRI_short:{
00388 short * bar = (short *) DSET_ARRAY(input_dset,iv) ;
00389 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
00390 if( mfac == 0.0 ) mfac = 1.0 ;
00391 if( do_mval ){
00392 float val ;
00393 for( ii=jj=0 ; ii < nvox ; ii++ ){
00394 if( mmm[ii] ){
00395 val = mfac*bar[ii] ;
00396 if( val >= val_bot && val <= val_top ) flar[jj++] = val ;
00397 }
00398 }
00399 mval = jj ;
00400 } else {
00401 for( ii=jj=0 ; ii < nvox ; ii++ )
00402 if( mmm[ii] ) flar[jj++] = mfac*bar[ii] ;
00403 }
00404 }
00405 break ;
00406
00407 case MRI_byte:{
00408 byte * bar = (byte *) DSET_ARRAY(input_dset,iv) ;
00409 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
00410 if( mfac == 0.0 ) mfac = 1.0 ;
00411 if( do_mval ){
00412 float val ;
00413 for( ii=jj=0 ; ii < nvox ; ii++ ){
00414 if( mmm[ii] ){
00415 val = mfac*bar[ii] ;
00416 if( val >= val_bot && val <= val_top ) flar[jj++] = val ;
00417 }
00418 }
00419 mval = jj ;
00420 } else {
00421 for( ii=jj=0 ; ii < nvox ; ii++ )
00422 if( mmm[ii] ) flar[jj++] = mfac*bar[ii] ;
00423 }
00424 }
00425 break ;
00426
00427 case MRI_float:{
00428 float * bar = (float *) DSET_ARRAY(input_dset,iv) ;
00429 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
00430 if( mfac == 0.0 ) mfac = 1.0 ;
00431 if( do_mval ){
00432 float val ;
00433 for( ii=jj=0 ; ii < nvox ; ii++ ){
00434 if( mmm[ii] ){
00435 val = mfac*bar[ii] ;
00436 if( val >= val_bot && val <= val_top ) flar[jj++] = val ;
00437 }
00438 }
00439 mval = jj ;
00440 } else {
00441 for( ii=jj=0 ; ii < nvox ; ii++ )
00442 if( mmm[ii] ) flar[jj++] = mfac*bar[ii] ;
00443 }
00444 }
00445 break ;
00446 }
00447
00448 if( do_mval ){
00449 if( mval == 0 ){
00450 free(mmm) ; mri_free(flim) ;
00451 return "*** Can't use source dataset -- no data in Values range ***" ;
00452 }
00453 flim->nx = flim->nvox = mcount = mval ;
00454 }
00455
00456
00457
00458 if( val_bot > val_top ){
00459 hbot = mri_min(flim) ; htop = mri_max(flim) ;
00460 if( hbot >= htop ){
00461 free(mmm) ; mri_free(flim) ;
00462 return "***********************************\n"
00463 "Selected voxels have no data range!\n"
00464 "***********************************" ;
00465 }
00466 } else {
00467 hbot = val_bot ; htop = val_top ;
00468 }
00469
00470 if( nbin < 10 || nbin > 1000 ){
00471 switch( DSET_BRICK_TYPE(input_dset,iv) ){
00472 case MRI_float:
00473 nbin = (int) sqrt((double)mcount) ;
00474 break ;
00475
00476 case MRI_short:
00477 case MRI_byte:{
00478 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
00479 if( mfac == 0.0 || mfac == 1.0 )
00480 nbin = (int)( htop - hbot ) ;
00481 else
00482 nbin = (int) sqrt((double)mcount) ;
00483 }
00484 break ;
00485
00486 }
00487 if( nbin < 10 ) nbin = 10 ; else if( nbin > 1000 ) nbin = 1000 ;
00488 }
00489
00490
00491
00492 hbin = (int *) calloc((nbin+1),sizeof(int)) ;
00493
00494 mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ;
00495
00496 if( smooth > 0 ){
00497 int nwid=smooth , *gbin=(int *)calloc((nbin+1),sizeof(int)) , ibot,itop ;
00498 float ws,wss , *wt ;
00499
00500 ws = 0.0 ;
00501 wt = (float *)malloc(sizeof(float)*(2*nwid+1)) ;
00502 for( ii=0 ; ii <= 2*nwid ; ii++ ){
00503 wt[ii] = nwid-abs(nwid-ii) + 0.5f ;
00504 ws += wt[ii] ;
00505 }
00506 for( ii=0 ; ii <= 2*nwid ; ii++ ) wt[ii] /= ws ;
00507
00508 for( jj=0 ; jj <= nbin ; jj++ ){
00509 ibot = jj-nwid ; if( ibot < 0 ) ibot = 0 ;
00510 itop = jj+nwid ; if( itop > nbin ) itop = nbin ;
00511 ws = wss = 0.0 ;
00512 for( ii=ibot ; ii <= itop ; ii++ ){
00513 ws += wt[nwid-jj+ii] * hbin[ii] ; wss += wt[nwid-jj+ii] ;
00514 }
00515 gbin[jj] = rint(ws/wss) ;
00516 }
00517 memcpy(hbin,gbin,sizeof(int)*(nbin+1)) ;
00518 free((void *)wt) ; free((void *)gbin) ;
00519 }
00520
00521 if( maxcount > 0 ){
00522 for( ii=0 ; ii <= nbin ; ii++ ) hbin[ii] = MIN( hbin[ii] , maxcount ) ;
00523 }
00524 sprintf(buf,"\\noesc %s[%d] %d voxels",DSET_FILECODE(input_dset),iv,mcount);
00525 PLUTO_histoplot( nbin,hbot,htop,hbin , NULL , NULL , buf , 0,NULL ) ;
00526
00527
00528
00529 if ( HISTOUT != NULL )
00530 {
00531 if( hbot >= htop ){ hbot = 0.0 ; htop = nbin ;}
00532
00533 dx = (htop-hbot)/nbin ;
00534
00535 for( ii=0 ; ii <= nbin ; ii++ )
00536 fprintf (HISTOUT, "%12.6f %13d \n", hbot+ii*dx, hbin[ii]) ;
00537
00538 fclose (HISTOUT) ;
00539
00540 fprintf (stderr, "%s written to disk \n", histout) ;
00541 }
00542
00543
00544
00545 free(hbin) ; free(mmm) ; mri_free(flim) ; return NULL ;
00546 }
00547
00548
00549
00550 #define FIT_FISHER
00551 #undef DO_GREEN
00552
00553 static char c_helpstring[] =
00554 "Purpose: Plot a histogram of correlation coefficient of a 3D+time\n"
00555 " dataset with an input time series file.\n"
00556 "\n"
00557 "Source: Dataset = 3D+time dataset to use\n"
00558 " Ignore = number of points at beginning to ignore\n"
00559 "\n"
00560 "Vector: 1D File = time series to correlate the dataset with\n"
00561 " Polort = degree of polynomial detrending to use\n"
00562 "\n"
00563 "Mask: Dataset = optional mask dataset to use, to restrict the voxels\n"
00564 " which will be processed\n"
00565 " Sub-brick = which sub-brick of the mask dataset to use\n"
00566 "\n"
00567 "Range: Bottom = minimum value of mask dataset to use\n"
00568 " Top = maximum value of mask dataset to use\n"
00569 " [default is to use all nonzero values in the mask]\n"
00570 "\n"
00571 "All selected voxel time series from the source are correlated with the\n"
00572 "input 1D vector, using the same algorithms as the AFNI internal FIM.\n"
00573 "The array of correlation coefficients is then put into 100 bins, ranging\n"
00574 "from -1.0 to 1.0, and the histogram graph is popped up to the display.\n"
00575 "\n"
00576 "Overlaid on the histogram are other graphs:\n"
00577 " * The first [red] is a normal fit to the Fisher z-transform of the\n"
00578 " correlation coefficient.\n"
00579 #ifdef DO_GREEN
00580 " * The second [green] is the nominal fit to the number of degrees of\n"
00581 " freedom, assuming the data time series are normal white noise.\n"
00582 #endif
00583 "\n"
00584 "-- Bob Cox - October 1999\n"
00585 ;
00586
00587 static char * CORREL_main( PLUGIN_interface * ) ;
00588
00589 static PLUGIN_interface * CORREL_init(void)
00590 {
00591 PLUGIN_interface * plint ;
00592
00593
00594
00595 plint = PLUTO_new_interface( "Histogram: CC" ,
00596 "Histogram of Correlation Coefficient" ,
00597 c_helpstring ,
00598 PLUGIN_CALL_VIA_MENU , CORREL_main ) ;
00599
00600 PLUTO_add_hint( plint , "Histogram: Correlation Coefficient" ) ;
00601
00602 PLUTO_set_sequence( plint , "A:afniinfo:dset" ) ;
00603
00604
00605
00606 PLUTO_add_option( plint , "Source" , "Source" , TRUE ) ;
00607
00608 PLUTO_add_dataset( plint ,
00609 "Dataset" ,
00610 ANAT_ALL_MASK ,
00611 FUNC_FIM_MASK ,
00612 DIMEN_4D_MASK |
00613 BRICK_ALLREAL_MASK
00614 ) ;
00615
00616 PLUTO_add_number( plint ,
00617 "Ignore" ,
00618 0 ,
00619 999 ,
00620 0 ,
00621 4 ,
00622 TRUE
00623 ) ;
00624
00625
00626
00627 PLUTO_add_option( plint , "Vector" , "Vector" , TRUE ) ;
00628
00629 PLUTO_add_timeseries( plint , "1D File" ) ;
00630
00631 PLUTO_add_number( plint , "Polort" , 0,MAX_POLORT,0,1,FALSE ) ;
00632
00633
00634
00635 PLUTO_add_option( plint , "Mask" , "Mask" , FALSE ) ;
00636 PLUTO_add_dataset( plint , "Dataset" ,
00637 ANAT_ALL_MASK , FUNC_ALL_MASK ,
00638 DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
00639 PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
00640
00641
00642
00643 PLUTO_add_option( plint , "Range" , "Range" , FALSE ) ;
00644 PLUTO_add_number( plint , "Bottom" , -99999,99999, 1, 0,1 ) ;
00645 PLUTO_add_number( plint , "Top" , -99999,99999,-1, 0,1 ) ;
00646
00647
00648 return plint ;
00649 }
00650
00651
00652
00653 #include "uuu.c"
00654
00655 static char * CORREL_main( PLUGIN_interface * plint )
00656 {
00657 MCW_idcode * idc ;
00658 THD_3dim_dataset * input_dset , * mask_dset = NULL ;
00659 MRI_IMAGE * tsim , * flim ;
00660 int ignore , nvox , ntim , polort , miv , it , ip , nupdt , nbin ;
00661 int mcount , ii , jj ;
00662 float * tsar ;
00663 float mask_bot=666.0 , mask_top=-666.0 , hbot,htop ;
00664 byte * mmm ;
00665 char buf[THD_MAX_NAME+16] , * tag ;
00666
00667 PCOR_references * pc_ref ;
00668 PCOR_voxel_corr * pc_vc ;
00669 int fim_nref ;
00670 float * ref_vec , * vval , * zval , * aval ;
00671 int * hbin , * jbin , * kbin , *jist[2] ;
00672 float sum , sumq , dbin , gval,rval,gg , sqp , zmid,zmed,zsig ;
00673 float pstar,zstar,zplus,zminus,psum,msum ;
00674
00675
00676
00677
00678 if( plint == NULL )
00679 return "************************\n"
00680 "CORREL_main: NULL input\n"
00681 "************************" ;
00682
00683
00684
00685 PLUTO_next_option(plint) ;
00686 idc = PLUTO_get_idcode(plint) ;
00687 input_dset = PLUTO_find_dset(idc) ;
00688 if( input_dset == NULL )
00689 return "*******************************\n"
00690 "CORREL_main: bad input dataset\n"
00691 "*******************************" ;
00692
00693 nvox = DSET_NVOX(input_dset) ;
00694 ntim = DSET_NUM_TIMES(input_dset) ;
00695
00696 ignore = (int) PLUTO_get_number(plint) ;
00697 if( ignore >= ntim-5 || ignore < 0 )
00698 return "******************************\n"
00699 "CORREL_main: bad ignore count\n"
00700 "******************************" ;
00701
00702 DSET_load(input_dset) ;
00703 if( DSET_ARRAY(input_dset,0) == NULL )
00704 return "********************************\n"
00705 "CORREL_main: can't load dataset\n"
00706 "********************************" ;
00707
00708
00709
00710 PLUTO_next_option(plint) ;
00711 tsim = PLUTO_get_timeseries(plint) ;
00712 if( tsim == NULL || tsim->nx < ntim )
00713 return "*****************************\n"
00714 "CORREL_main: bad input vector\n"
00715 "*****************************" ;
00716
00717 flim = mri_to_float(tsim) ; tsar = MRI_FLOAT_PTR(flim) ;
00718
00719 polort = (int) PLUTO_get_number(plint) ;
00720
00721
00722
00723 while( (tag=PLUTO_get_optiontag(plint)) != NULL ){
00724
00725
00726
00727 if( strcmp(tag,"Range") == 0 ){
00728 if( mask_dset == NULL ){
00729 mri_free(flim) ;
00730 return "******************************************\n"
00731 "CORREL_main: Can't use Range without Mask\n"
00732 "******************************************" ;
00733 }
00734
00735 mask_bot = PLUTO_get_number(plint) ;
00736 mask_top = PLUTO_get_number(plint) ;
00737 continue ;
00738 }
00739
00740
00741
00742 if( strcmp(tag,"Mask") == 0 ){
00743
00744 idc = PLUTO_get_idcode(plint) ;
00745 mask_dset = PLUTO_find_dset(idc) ;
00746
00747 if( mask_dset == NULL ){
00748 mri_free(flim) ;
00749 return "******************************\n"
00750 "CORREL_main: bad mask dataset\n"
00751 "******************************" ;
00752 }
00753
00754 if( DSET_NVOX(mask_dset) != nvox ){
00755 mri_free(flim) ;
00756 return "************************************************************\n"
00757 "CORREL_main: mask input dataset doesn't match source dataset\n"
00758 "************************************************************" ;
00759 }
00760
00761 miv = (int) PLUTO_get_number(plint) ;
00762 if( miv >= DSET_NVALS(mask_dset) || miv < 0 ){
00763 mri_free(flim) ;
00764 return "****************************************************\n"
00765 "CORREL_main: mask dataset sub-brick index is illegal\n"
00766 "****************************************************" ;
00767 }
00768
00769 DSET_load(mask_dset) ;
00770 if( DSET_ARRAY(mask_dset,miv) == NULL ){
00771 mri_free(flim) ;
00772 return "*************************************\n"
00773 "CORREL_main: can't load mask dataset\n"
00774 "*************************************" ;
00775 }
00776 continue ;
00777 }
00778 }
00779
00780
00781
00782
00783
00784
00785 if( mask_dset == NULL ){
00786 mmm = (byte *) malloc( sizeof(byte) * nvox ) ;
00787 if( mmm == NULL )
00788 return " \n*** Can't malloc workspace! ***\n" ;
00789 memset( mmm , 1, nvox ) ; mcount = nvox ;
00790 } else {
00791
00792 mmm = THD_makemask( mask_dset , miv , mask_bot , mask_top ) ;
00793 if( mmm == NULL )
00794 return " \n*** Can't make mask for some reason! ***\n" ;
00795 mcount = THD_countmask( nvox , mmm ) ;
00796
00797 if( !EQUIV_DSETS(mask_dset,input_dset) ) DSET_unload(mask_dset) ;
00798 if( mcount < 3 ){
00799 free(mmm) ;
00800 return " \n*** Less than 3 voxels survive the mask! ***\n" ;
00801 }
00802 sprintf(buf," \n"
00803 " %d voxels in the mask\n"
00804 " out of %d dataset voxels\n ",mcount,nvox) ;
00805 PLUTO_popup_transient(plint,buf) ;
00806 }
00807
00808
00809
00810 fim_nref = polort+2 ;
00811 ref_vec = (float *) malloc( sizeof(float) * fim_nref ) ;
00812 vval = (float *) malloc( sizeof(float) * mcount ) ;
00813
00814 pc_ref = new_PCOR_references( fim_nref ) ;
00815 pc_vc = new_PCOR_voxel_corr( mcount , fim_nref ) ;
00816
00817
00818
00819 for( nupdt=0,it=ignore ; it < ntim ; it++ ){
00820 if( tsar[it] >= WAY_BIG ) continue ;
00821
00822 ref_vec[0] = 1.0 ;
00823 for( ip=1 ; ip <= polort ; ip++ )
00824 ref_vec[ip] = ref_vec[ip-1] * ((float)it) ;
00825
00826 ref_vec[ip] = tsar[it] ;
00827
00828 update_PCOR_references( ref_vec , pc_ref ) ;
00829
00830
00831
00832 switch( DSET_BRICK_TYPE(input_dset,it) ){
00833
00834 case MRI_short:{
00835 short * bar = (short *) DSET_ARRAY(input_dset,it) ;
00836 float mfac = DSET_BRICK_FACTOR(input_dset,it) ;
00837 if( mfac == 0.0 ) mfac = 1.0 ;
00838 for( ii=jj=0 ; ii < nvox ; ii++ )
00839 if( mmm[ii] ) vval[jj++] = mfac*bar[ii] ;
00840 }
00841 break ;
00842
00843 case MRI_byte:{
00844 byte * bar = (byte *) DSET_ARRAY(input_dset,it) ;
00845 float mfac = DSET_BRICK_FACTOR(input_dset,it) ;
00846 if( mfac == 0.0 ) mfac = 1.0 ;
00847 for( ii=jj=0 ; ii < nvox ; ii++ )
00848 if( mmm[ii] ) vval[jj++] = mfac*bar[ii] ;
00849 }
00850 break ;
00851
00852 case MRI_float:{
00853 float * bar = (float *) DSET_ARRAY(input_dset,it) ;
00854 float mfac = DSET_BRICK_FACTOR(input_dset,it) ;
00855 if( mfac == 0.0 ) mfac = 1.0 ;
00856 for( ii=jj=0 ; ii < nvox ; ii++ )
00857 if( mmm[ii] ) vval[jj++] = mfac*bar[ii] ;
00858 }
00859 break ;
00860 }
00861
00862 PCOR_update_float( vval , pc_ref , pc_vc ) ;
00863 nupdt++ ;
00864 }
00865
00866 free(ref_vec) ; mri_free(flim) ; free(mmm) ;
00867
00868
00869
00870 PCOR_get_pcor( pc_ref , pc_vc , vval ) ;
00871
00872 free_PCOR_references(pc_ref) ; free_PCOR_voxel_corr(pc_vc) ;
00873
00874
00875
00876 sum = 0.0 ;
00877 for( ii=0 ; ii < mcount ; ii++ ) sum += vval[ii] ;
00878 sum /= mcount ; sumq = 0.0 ;
00879 for( ii=0 ; ii < mcount ; ii++ )
00880 sumq += (vval[ii]-sum)*(vval[ii]-sum) ;
00881 sumq = sqrt(sumq/mcount) ;
00882
00883
00884
00885 zval = (float *) malloc( sizeof(float) * mcount ) ;
00886 aval = (float *) malloc( sizeof(float) * mcount ) ;
00887 for( ii=0 ; ii < mcount ; ii++ ) zval[ii] = atanh(vval[ii]) ;
00888 qsort_float( mcount , zval ) ;
00889 if( mcount%2 == 1 )
00890 zmid = zval[mcount/2] ;
00891 else
00892 zmid = 0.5 * ( zval[mcount/2] + zval[mcount/2-1] ) ;
00893
00894 for( ii=0 ; ii < mcount ; ii++ ) aval[ii] = fabs(zval[ii]-zmid) ;
00895 qsort_float( mcount , aval ) ;
00896 if( mcount%2 == 1 )
00897 zmed = aval[mcount/2] ;
00898 else
00899 zmed = 0.5 * ( aval[mcount/2] + aval[mcount/2-1] ) ;
00900 zsig = 1.4826 * zmed ;
00901
00902 free(aval) ;
00903 free(zval) ;
00904
00905
00906
00907 hbot = -1.0 ; htop = 1.0 ; nbin = 100 ; dbin = (htop-hbot)/nbin ;
00908
00909 hbin = (int *) calloc((nbin+1),sizeof(int)) ;
00910 jbin = (int *) calloc((nbin+1),sizeof(int)) ;
00911
00912 #ifndef FIT_FISHER
00913 sqp = 1.0/(sqrt(2.0*PI)*sumq) ;
00914 for( ii=0 ; ii < nbin ; ii++ ){
00915 gval = hbot + (ii+0.5)*dbin - sum ;
00916 gval = sqp * exp( -0.5*gval*gval/(sumq*sumq) ) ;
00917 jbin[ii] = (int)( mcount * dbin * gval + 0.5 ) ;
00918 }
00919 #else
00920 sqp = 1.0/(sqrt(2.0*PI)*zsig) ;
00921 for( ii=0 ; ii < nbin ; ii++ ){
00922 rval = hbot + (ii+0.5)*dbin ;
00923 gval = atanh(rval) - zmid ;
00924 gval = sqp * exp( -0.5*gval*gval/(zsig*zsig) )/sqrt(1.0-rval*rval) ;
00925 jbin[ii] = (int)( mcount * dbin * gval + 0.5 ) ;
00926 }
00927 sum = tanh(zmid) ;
00928 sumq = 0.5*( tanh(zmid+zmed) - tanh(zmid-zmed) ) ;
00929 #endif
00930
00931 kbin = (int *) calloc((nbin+1),sizeof(int)) ;
00932 for( ii=0 ; ii < nbin ; ii++ ){
00933 gval = correl_t2p( fabs(hbot+ii*dbin) ,
00934 (double)nupdt , (double)1 , (double)(polort+1) ) ;
00935 sqp = correl_t2p( fabs(hbot+(ii+1)*dbin) ,
00936 (double)nupdt , (double)1 , (double)(polort+1) ) ;
00937 kbin[ii] = (int)( 0.5*mcount*fabs(gval-sqp) ) ;
00938 }
00939 jist[0] = jbin ; jist[1] = kbin ;
00940
00941 flim = mri_new_vol_empty( mcount,1,1 , MRI_float ) ;
00942 mri_fix_data_pointer( vval , flim ) ;
00943 mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ;
00944
00945 { char * ps = my_getenv("PTAIL") ;
00946 float pp=0.0 ;
00947 if( ps != NULL ) pp = strtod(ps,NULL) ;
00948 set_unusuality_tail(pp) ;
00949 }
00950
00951 for( ii=0 ; ii < mcount ; ii++ ) vval[ii] = -vval[ii] ;
00952 msum = unusuality( mcount , vval ) ;
00953 for( ii=0 ; ii < mcount ; ii++ ) vval[ii] = -vval[ii] ;
00954 psum = unusuality( mcount , vval ) ;
00955
00956 sprintf(buf,"%s^.%s:\\rho_{mid}=%.2f\\pm%.2f,u=%.1f,%.1f",
00957 DSET_FILECODE(input_dset),
00958 (tsim->name != NULL) ? THD_trailname(tsim->name,0) : " " ,
00959 sum,sumq , psum,msum ) ;
00960 #ifdef DO_GREEN
00961 PLUTO_histoplot( nbin,hbot,htop,hbin ,
00962 "Correlation Coefficient",NULL,buf , 2,jist ) ;
00963 #else
00964 PLUTO_histoplot( nbin,hbot,htop,hbin ,
00965 "Correlation Coefficient",NULL,buf , 1,jist ) ;
00966 #endif
00967
00968 mri_clear_data_pointer(flim) ; mri_free(flim) ;
00969 free(vval) ; free(hbin) ; free(jbin) ; free(kbin) ;
00970 return NULL ;
00971 }
