Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search
uuu3.c

00001 #include "mrilib.h"
00002 
00003 /*------------------------------------------------------------
00004   Set the one-sided tail probability at which we will cutoff
00005   the unusuality test.
00006 --------------------------------------------------------------*/
00007 
00008 static float zstar = 0.0 ;            /* the actual cutoff */
00009 static float pstar = 0.0 ;            /* tail probability  */
00010 
00011 void set_unusuality_tail( float p )
00012 {
00013    if( p > 0.0 && p < 1.0 ){
00014       zstar = qginv(p) ;
00015       pstar = p ;
00016    }
00017    return ;
00018 }
00019 
00020 #define NEAR1 0.999
00021 
00022 /*------------------------------------------------------------
00023   Inputs:  rr[0..nr-1] = array of correlation coefficients.
00024   Outputs: ihi[0..*nhi-1] = indexes of highly correlations
00025              from rr.  ihi must be declared at least length
00026              nr to allow for the worst case.
00027 --------------------------------------------------------------*/
00028 
00029 void find_unusual_correlations( int nrin  , float * rr ,
00030                                 int * nhi , int * ihi   )
00031 {
00032    register int ii,jj,nr ;
00033    register float *zz ;
00034    float zmid,zsig,zmed , rmid,rcut ;
00035 
00036    if( nhi == NULL ) return ;  /* illegal entry */
00037 
00038    *nhi = 0 ;                  /* default return */
00039 
00040    if( nrin < 1000 || rr == NULL || ihi == NULL  ) return ;  /* illegal entry */
00041 
00042    /*-- make workspace --*/
00043 
00044    zz = (float *) malloc(sizeof(float)*nrin) ;
00045 
00046    if( zstar <= 0.0 ){               /* initialize zstar if needed */
00047       char * cp = getenv("PTAIL") ;
00048       float pp = 0.0001 ;
00049       if( cp != NULL ){
00050          float xx = strtod( cp , NULL ) ;
00051          if( xx > 0.0 && xx < 1.0 ) pp = xx ;
00052       }
00053       set_unusuality_tail( pp ) ;
00054    }
00055 
00056    /*-- copy data into workspace (eliding 1's) --*/
00057 
00058    for( ii=jj=0 ; ii < nrin ; ii++ )
00059       if( rr[ii] <= NEAR1 ) zz[jj++] = rr[ii] ;
00060    nr = jj ;
00061    if( nr < 2 ){ free(zz); return; } /* shouldn't happen */
00062 
00063    rmid = qmed_float( nr , zz ) ;  /* median of correlations */
00064    zmid = atanh(rmid) ;            /* median of atanh(rr) */
00065 
00066    /*-- zz <- fabs( tanh( atanh(zz) - atanh(rmid) ) ) --*/
00067 
00068    for( ii=0 ; ii < nr ; ii++ )
00069       zz[ii] = fabs( (zz[ii]-rmid)/(1.0-zz[ii]*rmid) ) ;
00070 
00071    /*-- find MAD of atanh(rr) --*/
00072 
00073    zmed = qmed_float( nr , zz ) ; /* MAD = median absolute deviation of rr */
00074    zmed = atanh(zmed) ;           /* MAD of atanh(rr) */
00075    zsig = 1.4826 * zmed ;         /* estimate standard dev. of atanh(rr) */
00076                                   /* 1.4826 = 1/(sqrt(2)*erfinv(0.5))    */
00077    free(zz) ;
00078 
00079    if( zsig <= 0.0 ) return ; /* shouldn't happen */
00080 
00081    /*-- find values of correlation greater than threshold --*/
00082    /*-- that is, with (atanh(rr)-zmid)/zsig > zstar       --*/
00083 
00084    rcut = tanh( zsig*zstar + zmid ) ;
00085 
00086    for( ii=jj=0 ; ii < nrin ; ii++ )
00087       if( rr[ii] >= rcut && rr[ii] <= NEAR1 ) ihi[jj++] = ii ;
00088 
00089    *nhi = jj ; return ;
00090 }
AFNI/NIfTI Server

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

uuu3.c
