00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 #include <string.h>
00008 #include "mrilib.h"
00009 #include <stdlib.h>
00010 #include <ctype.h>
00011 
00012 /*-------------------------- global data --------------------------*/
00013 
00014 /** inputs **/
00015 
00016 static THD_3dim_dataset * UC_dset = NULL ; /* dataset */
00017 
00018 static char UC_prefix[THD_MAX_PREFIX] = "uuu2" ;
00019 
00020 static int UC_be_quiet = 1 ;
00021 
00022 static byte * UC_mask = NULL ;
00023 static int    UC_mask_nvox = 0 ;
00024 static int    UC_mask_hits = 0 ;
00025 
00026 static int    UC_nvec = 0 ;      /* # of vectors to use from dataset */
00027 static int    UC_vdim = 0 ;      /* length of each vector            */
00028 
00029 static float ** UC_vec = NULL ;  /* UC_vec[k][i] is the i-th component  */
00030                                  /* of the k-th vector 0 <= i < UC_vdim */
00031                                  /*                    0 <= k < UC_nvec */
00032 
00033 static int    * UC_iv  = NULL ;  /* UC_vec[k] comes from voxel # UC_iv[k] */
00034 
00035 static float UC_ptail = 0.0001 ;
00036 static int   UC_mtail = 2 ;
00037 
00038 void UC_syntax(char * msg) ;
00039 
00040 #include "uuu3.c"
00041 
00042 #include "graph_compon.c"
00043 
00044 /*-------------------------------------------------------------------
00045      detrend: routine to remove unwanted components from time series
00046 ---------------------------------------------------------------------*/
00047 
00048 void detrend( int n , float vec[] )
00049 {
00050    register int ii ;
00051    register float sum0 , sum1 , cf , lf ;
00052    float sum2 , det ;
00053 
00054    static int nold = -1 ;             /* initialization flag */
00055    static float cf0,cf1 , lf0,lf1 ;   /* to be initialized */
00056 
00057  /*** initialize coefficients for detrending ***/
00058 
00059    if( n != nold ){
00060       nold = n ; sum0 = sum1 = sum2 = 0.0 ;
00061       for( ii=0 ; ii < n ; ii++ ){
00062          sum0 += 1.0 ; sum1 += ii ; sum2 += ii*ii ;
00063       }
00064       det = sum0 * sum2 - sum1 * sum1 ;
00065       cf0 =  sum2 / det ;     /* constant factor for sum0 */
00066       cf1 = -sum1 / det ;     /* constant factor for sum1 */
00067       lf0 = cf1 ;             /* linear factor for sum0 */
00068       lf1 =  sum0 / det ;     /* linear factor for sum1 */
00069    }
00070 
00071  /*** remove mean and linear trend ***/
00072 
00073    sum0 = sum1 = 0.0 ;
00074    for( ii=0 ; ii < n ; ii++ ){
00075       sum0 += vec[ii] ; sum1 += vec[ii] * ii ;
00076    }
00077 
00078    cf = cf0 * sum0 + cf1 * sum1 ;
00079    lf = lf0 * sum0 + lf1 * sum1 ;
00080    for( ii=0 ; ii < n ; ii++ ) vec[ii] -= cf + ii*lf ;
00081 }
00082 
00083 /*----------------------------------------------------------
00084    normalize: routine to scale a time series to unit vector
00085 ------------------------------------------------------------*/
00086 
00087 void normalize( int n , float vec[] )
00088 {
00089    register int ii ;
00090    register float sqsum ;
00091 
00092    detrend( n , vec ) ;
00093 
00094    sqsum = 0.0 ;
00095    for( ii=0 ; ii < n ; ii++ ) sqsum += vec[ii] * vec[ii] ;
00096 
00097    if( sqsum < 1.e-10 ){
00098       for( ii=0 ; ii < n ; ii++ ) vec[ii] = 0.0 ;
00099    } else {
00100       sqsum = 1.0 / sqrt(sqsum) ;
00101       for( ii=0 ; ii < n ; ii++ ) vec[ii] *= sqsum ;
00102    }
00103 }
00104 
00105 /*------------------------------------------------------------------*/
00106 
00107 void UC_read_opts( int argc , char * argv[] )
00108 {
00109    int nopt = 1 ;
00110    float val ;
00111    int  kk, nxyz, mm,nn ;
00112    float * vv , * bb ;
00113 
00114    while( nopt < argc && argv[nopt][0] == '-' ){
00115 
00116       /**** -verbose ****/
00117 
00118       if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00119          UC_be_quiet = 0 ;
00120          nopt++ ; continue ;
00121       }
00122 
00123       /**** -prefix prefix ****/
00124 
00125       if( strncmp(argv[nopt],"-prefix",6) == 0 ){
00126          nopt++ ;
00127          if( nopt >= argc ) UC_syntax("-prefix needs an argument!") ;
00128          MCW_strncpy( UC_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00129          continue ;
00130       }
00131 
00132       /**** -mask mset ****/
00133 
00134       if( strncmp(argv[nopt],"-mask",5) == 0 ){
00135          THD_3dim_dataset * mset ; int ii,nn ;
00136          nopt++ ;
00137          if( nopt >= argc ) UC_syntax("need arguments after -mask!") ;
00138          mset = THD_open_dataset( argv[nopt] ) ;
00139          if( mset == NULL ) UC_syntax("can't open -mask dataset!") ;
00140          UC_mask = THD_makemask( mset , 0 , 1.0,0.0 ) ;
00141          UC_mask_nvox = DSET_NVOX(mset) ;
00142          DSET_delete(mset) ;
00143          if( UC_mask == NULL ) UC_syntax("can't use -mask dataset!") ;
00144          UC_mask_hits = THD_countmask( UC_mask_nvox , UC_mask ) ;
00145          if( UC_mask_hits == 0 ) UC_syntax("mask is all zeros!") ;
00146          if( !UC_be_quiet ) printf("--- %d voxels in mask\n",UC_mask_hits) ;
00147 
00148          UC_iv = (int *) malloc( sizeof(int) * UC_mask_hits ) ;
00149          for( nn=ii=0 ; ii < UC_mask_nvox ; ii++ )
00150             if( UC_mask[ii] ) UC_iv[nn++] = ii ;
00151 
00152          nopt++ ; continue ;
00153       }
00154 
00155       /**** -ptail p ****/
00156 
00157       if( strcmp(argv[nopt],"-ptail") == 0 ){
00158          if( ++nopt >= argc ) UC_syntax("-ptail needs an argument!") ;
00159          UC_ptail = strtod( argv[nopt] , NULL ) ;
00160          if( UC_ptail <= 0.0 || UC_ptail >= 0.499 )
00161             UC_syntax("value after -ptail is illegal!") ;
00162          nopt++ ; continue ;
00163       }
00164 
00165       /**** -mtail m ****/
00166 
00167       if( strcmp(argv[nopt],"-mtail") == 0 ){
00168          if( ++nopt >= argc ) UC_syntax("-mtail needs an argument!") ;
00169          UC_mtail = strtod( argv[nopt] , NULL ) ;
00170          nopt++ ; continue ;
00171       }
00172 
00173       /**** unknown switch ****/
00174 
00175       fprintf(stderr,"\n*** unrecognized option %s\n",argv[nopt]) ;
00176       exit(1) ;
00177 
00178    }  /* end of loop over options */
00179 
00180    /*--- a simple consistency check ---*/
00181 
00182    /*--- last input is dataset name ---*/
00183 
00184    if( nopt >= argc ) UC_syntax("no input dataset name?") ;
00185 
00186    UC_dset = THD_open_dataset( argv[nopt] ) ;
00187    if( !ISVALID_3DIM_DATASET(UC_dset) ){
00188       fprintf(stderr,"\n*** can't open dataset file %s\n",argv[nopt]) ;
00189       exit(1) ;
00190    }
00191 
00192    nxyz = DSET_NVOX(UC_dset) ;
00193    if( UC_mask != NULL && nxyz != UC_mask_nvox )
00194       UC_syntax("mask and input dataset size mismatch!") ;
00195 
00196    /*--- load vectors ---*/
00197 
00198    UC_nvec = (UC_mask_hits > 0) ? UC_mask_hits : nxyz ;
00199    UC_vdim = DSET_NVALS(UC_dset) ;
00200    if( UC_vdim < 4 )
00201       UC_syntax("input dataset needs at least 4 sub-bricks!") ;
00202 
00203    vv     = (float *) malloc( sizeof(float) * UC_nvec * UC_vdim ) ;
00204    UC_vec = (float **) malloc( sizeof(float *) * UC_nvec ) ;
00205    for( kk=0 ; kk < UC_nvec ; kk++ ) UC_vec[kk] = vv + (kk*UC_vdim) ;
00206 
00207    if( !UC_be_quiet ) printf("--- reading input dataset\n") ;
00208    DSET_load(UC_dset) ;
00209    if( ! DSET_LOADED(UC_dset) )
00210       UC_syntax("Can't load input dataset bricks!") ;
00211 
00212    /* copy brick data into float storage */
00213 
00214    if( !UC_be_quiet ) printf("--- loading vectors\n") ;
00215 
00216    bb = (float *) malloc( sizeof(float) * nxyz ) ;
00217    for( mm=0 ; mm < UC_vdim ; mm++ ){
00218 
00219       EDIT_coerce_type( nxyz ,
00220                         DSET_BRICK_TYPE(UC_dset,mm) , DSET_ARRAY(UC_dset,mm) ,
00221                         MRI_float , bb ) ;
00222 
00223       DSET_unload_one( UC_dset , mm ) ;
00224 
00225       if( UC_mask == NULL ){
00226          for( kk=0 ; kk < nxyz ; kk++ ) UC_vec[kk][mm] = bb[kk] ;
00227       } else {
00228          for( nn=kk=0 ; kk < nxyz ; kk++ )
00229             if( UC_mask[kk] ) UC_vec[nn++][mm] = bb[kk] ;
00230       }
00231    }
00232    free(bb) ; DSET_unload( UC_dset ) ;
00233 
00234    /* detrend and normalize vectors */
00235 
00236    if( !UC_be_quiet ) printf("--- normalizing vectors\n") ;
00237 
00238    for( kk=0 ; kk < UC_nvec ; kk++ )
00239       normalize( UC_vdim , UC_vec[kk] ) ;
00240 
00241    return ;
00242 }
00243 
00244 /*-----------------------------------------------------------------------
00245   Compute the unusuality index of a reference vector in a set of vectors
00246 -------------------------------------------------------------------------*/
00247 
00248 static int * UC_ihi = NULL ;  /* cheap, but it works */
00249 
00250 int UC_unusuality( int ndim , float * ref , int nvec , float ** vec )
00251 {
00252    register int ii , kk ;
00253    register float psum , * vv ;
00254    int nhi ;
00255 
00256    static int     nvold = -1   ;
00257    static float * zval  = NULL ;
00258 
00259    if( ndim < 4 || nvec < 4 || ref == NULL || vec == NULL ) return 0 ;
00260 
00261    /* initialize if number of vectors has changed */
00262 
00263    if( nvold != nvec ){
00264       if( zval != NULL ) free(zval) ;
00265       if( UC_ihi  != NULL ) free(UC_ihi) ;
00266       zval = (float *) malloc(sizeof(float)*nvec) ;
00267       UC_ihi  = (int *)   malloc(sizeof(int)  *nvec) ;
00268       nvold = nvec ;
00269    }
00270 
00271    /* compute dot products */
00272 
00273    for( kk=0 ; kk < nvec ; kk++ ){
00274       psum = 0.0 ; vv = vec[kk] ;
00275       for( ii=0 ; ii < ndim ; ii++ ) psum += ref[ii] * vv[ii] ;
00276       zval[kk] = psum ;
00277    }
00278 
00279    find_unusual_correlations( nvec , zval , &nhi , UC_ihi ) ;
00280 
00281    return (short) nhi ;
00282 }
00283 
00284 /*---------------------------------------------------------------------------*/
00285 
00286 void UC_syntax(char * msg)
00287 {
00288    if( msg != NULL ){ fprintf(stderr,"\n*** %s\n",msg) ; exit(1) ; }
00289 
00290    printf(
00291     "Usage: 3duuu2 [options] dataset ...\n"
00292     "\n"
00293     "The input dataset may have a sub-brick selector list.\n"
00294     "Otherwise, all sub-bricks from a dataset will be used.\n"
00295     "\n"
00296     "OPTIONS:\n"
00297     "  -prefix pname \n"
00298     "  -verbose\n"
00299     "  -mask mset\n"
00300     "  -ptail p\n"
00301     "  -mtail m\n"
00302    ) ;
00303 
00304    exit(0) ;
00305 }
00306 
00307 /*------------------------------------------------------------------*/
00308 
00309 int main( int argc , char * argv[] )
00310 {
00311    int kk , nvox , ii , jj , uval,ncom , aa ;
00312    THD_3dim_dataset * oset ;
00313    short * sar ;
00314    int ** gmat , ** cmat ;
00315 
00316    /*-- read command line arguments --*/
00317 
00318    if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ) UC_syntax(NULL) ;
00319 
00320    (void) my_getenv("junk") ;
00321 
00322    UC_read_opts( argc , argv ) ;
00323    set_unusuality_tail( UC_ptail ) ;
00324 
00325    oset = EDIT_empty_copy( UC_dset ) ;
00326    EDIT_dset_items( oset ,
00327                        ADN_prefix      , UC_prefix ,
00328                        ADN_ntt         , 0 ,
00329                        ADN_nvals       , 1 ,
00330                        ADN_func_type   , ANAT_BUCK_TYPE ,
00331                        ADN_datum_all   , MRI_short ,
00332                        ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00333                     ADN_none ) ;
00334 
00335    nvox = DSET_NVOX(oset) ;
00336    sar = (short *) calloc( nvox , sizeof(short) ) ;
00337    EDIT_substitute_brick( oset , 0 , MRI_short , sar ) ;
00338 
00339    gmat = (int **) malloc( sizeof(int *) * UC_nvec ) ;
00340 
00341    if( !UC_be_quiet ){ printf("--- computing u") ; fflush(stdout) ; }
00342 
00343    for( kk=0 ; kk < UC_nvec ; kk++ ){
00344       ii = (UC_iv == NULL) ? kk : UC_iv[kk] ;
00345       uval = UC_unusuality( UC_vdim, UC_vec[kk], UC_nvec, UC_vec ) ;
00346 
00347       if( uval < UC_mtail ) uval = 0 ;
00348 
00349       sar[ii] = uval ;
00350 
00351       /* make graph data */
00352 
00353       gmat[kk]    = (int *) malloc( sizeof(int) * (uval+1) ) ;
00354       gmat[kk][0] = uval ;
00355       for( jj=0 ; jj < uval ; jj++ ) gmat[kk][jj+1] = UC_ihi[jj] ;
00356 
00357       if( !UC_be_quiet && kk%1000==999 ){
00358          printf("%d",(kk/1000)%10);fflush(stdout);
00359       }
00360    }
00361    if( !UC_be_quiet ) printf("\n") ;
00362 
00363    if( !UC_be_quiet ) printf("--- fixing graph\n") ;
00364 
00365 #undef ADDTHEM
00366 
00367    for( kk=0 ; kk < UC_nvec ; kk++ ){   /* loop over pts */
00368       uval = gmat[kk][0] ;              /* # pts connected to kk */
00369       for( jj=0 ; jj < uval ; jj++ ){   /* loop over pts connected to kk */
00370          ii = gmat[kk][jj+1] ;          /* jj-th pt connected to kk */
00371 
00372          for( aa=1 ; aa <= gmat[ii][0] ; aa++ )  /* see if kk is in ii's list */
00373             if( gmat[ii][aa] == kk ) break ;
00374 
00375          if( aa > gmat[ii][0] ){        /* wasn't in list */
00376 #ifdef ADDTHEM
00377            /* add kk to ii's list */
00378 
00379             gmat[ii] = (int *) realloc( sizeof(int)*(gmat[ii][0]+2) ) ;
00380             gmat[ii][++(gmat[ii][0])] = kk ;
00381 #else
00382            /* remove ii from kk's list */
00383 
00384            gmat[kk][jj+1] = -1 ;  /* flag as a bad connection */
00385 #endif
00386          }
00387 
00388       } /* end of loop over pts connected to kk */
00389    } /* end of loop over pts */
00390 
00391    if( !UC_be_quiet ) printf("--- finding components\n") ;
00392 
00393    GRAPH_find_components( UC_nvec , gmat , &ncom , &cmat ) ;
00394 
00395    if( !UC_be_quiet ) printf("--- found %d components\n",ncom) ;
00396 
00397    sar = (short *) calloc( nvox , sizeof(short) ) ;
00398    EDIT_add_brick( oset , MRI_short , 0.0 , sar ) ;
00399 
00400    for( kk=0 ; kk < ncom ; kk++ ){
00401       if( !UC_be_quiet )
00402          printf("--- component %d has %d voxels\n",kk,cmat[kk][0]) ;
00403 
00404       if( cmat[kk][0] < 2 ) break ;
00405 
00406       for( ii=1 ; ii <= cmat[kk][0] ; ii++ ){
00407          jj = (UC_iv == NULL) ? cmat[kk][ii] : UC_iv[cmat[kk][ii]] ;
00408          sar[jj] = (kk+1) ;
00409       }
00410    }
00411 
00412    if( !UC_be_quiet ) printf("--- writing output\n") ;
00413 
00414    DSET_write(oset) ;
00415    exit(0) ;
00416 }