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 MRI_IMAGE * UC_ref = NULL ;         /* vector */
00019 
00020 static char UC_prefix[THD_MAX_PREFIX] = "uc" ;
00021 
00022 static int UC_be_quiet = 1 ;
00023 
00024 static byte * UC_mask = NULL ;
00025 static int    UC_mask_nvox = 0 ;
00026 static int    UC_mask_hits = 0 ;
00027 
00028 static int    UC_nvec = 0 ;
00029 static int    UC_vdim = 0 ;
00030 
00031 static float ** UC_vec = NULL ;  /* UC_vec[k][i] is the i-th component  */
00032                                  /* of the k-th vector 0 <= i < UC_vdim */
00033                                  /*                    0 <= k < UC_nvec */
00034 void UC_syntax(char * msg) ;
00035 
00036 /*-------------------------------------------------------------------
00037      detrend: routine to remove unwanted components from time series
00038 ---------------------------------------------------------------------*/
00039 
00040 void detrend( int n , float vec[] )
00041 {
00042    register int ii ;
00043    register float sum0 , sum1 , cf , lf ;
00044    float sum2 , det ;
00045 
00046    static int nold = -1 ;             /* initialization flag */
00047    static float cf0,cf1 , lf0,lf1 ;   /* to be initialized */
00048 
00049  /*** initialize coefficients for detrending ***/
00050 
00051    if( n != nold ){
00052       nold = n ; sum0 = sum1 = sum2 = 0.0 ;
00053       for( ii=0 ; ii < n ; ii++ ){
00054          sum0 += 1.0 ; sum1 += ii ; sum2 += ii*ii ;
00055       }
00056       det = sum0 * sum2 - sum1 * sum1 ;
00057       cf0 =  sum2 / det ;     /* constant factor for sum0 */
00058       cf1 = -sum1 / det ;     /* constant factor for sum1 */
00059       lf0 = cf1 ;             /* linear factor for sum0 */
00060       lf1 =  sum0 / det ;     /* linear factor for sum1 */
00061    }
00062 
00063  /*** remove mean and linear trend ***/
00064 
00065    sum0 = sum1 = 0.0 ;
00066    for( ii=0 ; ii < n ; ii++ ){
00067       sum0 += vec[ii] ; sum1 += vec[ii] * ii ;
00068    }
00069 
00070    cf = cf0 * sum0 + cf1 * sum1 ;
00071    lf = lf0 * sum0 + lf1 * sum1 ;
00072    for( ii=0 ; ii < n ; ii++ ) vec[ii] -= cf + ii*lf ;
00073 }
00074 
00075 /*----------------------------------------------------------
00076    normalize: routine to scale a time series to unit vector
00077 ------------------------------------------------------------*/
00078 
00079 void normalize( int n , float vec[] )
00080 {
00081    register int ii ;
00082    register float sqsum ;
00083 
00084    detrend( n , vec ) ;
00085 
00086    sqsum = 0.0 ;
00087    for( ii=0 ; ii < n ; ii++ ) sqsum += vec[ii] * vec[ii] ;
00088 
00089    if( sqsum < 1.e-10 ){
00090       for( ii=0 ; ii < n ; ii++ ) vec[ii] = 0.0 ;
00091    } else {
00092       sqsum = 1.0 / sqrt(sqsum) ;
00093       for( ii=0 ; ii < n ; ii++ ) vec[ii] *= sqsum ;
00094    }
00095 }
00096 
00097 /*------------------------------------------------------------------*/
00098 
00099 void UC_read_opts( int argc , char * argv[] )
00100 {
00101    int nopt = 1 ;
00102    float val ;
00103    int  kk, nxyz, mm,nn ;
00104    float * vv , * bb ;
00105 
00106    while( nopt < argc && argv[nopt][0] == '-' ){
00107 
00108       /**** -verbose ****/
00109 
00110       if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00111          UC_be_quiet = 0 ;
00112          nopt++ ; continue ;
00113       }
00114 
00115       /**** -ref file.1D ****/
00116 
00117       if( strncmp(argv[nopt],"-ref",4) == 0 ){
00118          MRI_IMAGE * im ;
00119          nopt++ ;
00120          if( nopt >= argc ) UC_syntax("-ref needs an argument!") ;
00121          im = mri_read( argv[nopt] ) ;
00122          if( im == NULL ) UC_syntax("Can't read -ref file!") ;
00123          if( im->kind == MRI_float ){
00124             UC_ref = im ;
00125          } else {
00126             UC_ref = mri_to_float(im) ; mri_free(im) ;
00127          }
00128          im = mri_transpose(UC_ref) ; mri_free(UC_ref) ; UC_ref = im ;
00129          nopt++ ; continue ;
00130       }
00131 
00132       /**** -prefix prefix ****/
00133 
00134       if( strncmp(argv[nopt],"-prefix",6) == 0 ){
00135          nopt++ ;
00136          if( nopt >= argc ) UC_syntax("-prefix needs an argument!") ;
00137          MCW_strncpy( UC_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00138          continue ;
00139       }
00140 
00141       /**** -mask mset ****/
00142 
00143       if( strncmp(argv[nopt],"-mask",5) == 0 ){
00144          THD_3dim_dataset * mset ; int ii ;
00145          nopt++ ;
00146          if( nopt >= argc ) UC_syntax("need arguments after -mask!") ;
00147          mset = THD_open_dataset( argv[nopt] ) ;
00148          if( mset == NULL ) UC_syntax("can't open -mask dataset!") ;
00149          UC_mask = THD_makemask( mset , 0 , 1.0,0.0 ) ;
00150          UC_mask_nvox = DSET_NVOX(mset) ;
00151          DSET_delete(mset) ;
00152          if( UC_mask == NULL ) UC_syntax("can't use -mask dataset!") ;
00153          UC_mask_hits = THD_countmask( UC_mask_nvox , UC_mask ) ;
00154          if( UC_mask_hits == 0 ) UC_syntax("mask is all zeros!") ;
00155          if( !UC_be_quiet ) printf("--- %d voxels in mask\n",UC_mask_hits) ;
00156          nopt++ ; continue ;
00157       }
00158 
00159       /**** unknown switch ****/
00160 
00161       fprintf(stderr,"\n*** unrecognized option %s\n",argv[nopt]) ;
00162       exit(1) ;
00163 
00164    }  /* end of loop over options */
00165 
00166    /*--- a simple consistency check ---*/
00167 
00168    /*--- last input is dataset name ---*/
00169 
00170    if( nopt >= argc ) UC_syntax("no input dataset name?") ;
00171 
00172    UC_dset = THD_open_dataset( argv[nopt] ) ;
00173    if( !ISVALID_3DIM_DATASET(UC_dset) ){
00174       fprintf(stderr,"\n*** can't open dataset file %s\n",argv[nopt]) ;
00175       exit(1) ;
00176    }
00177 
00178    nxyz = DSET_NVOX(UC_dset) ;
00179    if( UC_mask != NULL && nxyz != UC_mask_nvox )
00180       UC_syntax("mask and input dataset size mismatch!") ;
00181 
00182    /*--- load vectors ---*/
00183 
00184    UC_nvec = (UC_mask_hits > 0) ? UC_mask_hits : nxyz ;
00185    UC_vdim = DSET_NVALS(UC_dset) ;
00186    if( UC_vdim < 4 )
00187       UC_syntax("input dataset needs at least 4 sub-bricks!") ;
00188 
00189    if( UC_ref == NULL || UC_ref->nx < UC_vdim )
00190       UC_syntax("input ref not long enough for input dataset!") ;
00191 
00192    vv     = (float *) malloc( sizeof(float) * UC_nvec * UC_vdim ) ;
00193    UC_vec = (float **) malloc( sizeof(float *) * UC_nvec ) ;
00194    for( kk=0 ; kk < UC_nvec ; kk++ ) UC_vec[kk] = vv + (kk*UC_vdim) ;
00195 
00196    if( !UC_be_quiet ) printf("--- reading dataset\n") ;
00197    DSET_load(UC_dset) ;
00198    if( ! DSET_LOADED(UC_dset) )
00199       UC_syntax("Can't load input dataset bricks!") ;
00200 
00201    /* copy brick data into float storage */
00202 
00203    if( !UC_be_quiet ) printf("--- loading vectors\n") ;
00204 
00205    bb = (float *) malloc( sizeof(float) * nxyz ) ;
00206    for( mm=0 ; mm < UC_vdim ; mm++ ){
00207 
00208       EDIT_coerce_type( nxyz , DSET_BRICK_TYPE(UC_dset,mm) ,
00209                                DSET_ARRAY(UC_dset,mm) ,
00210                         MRI_float , bb ) ;
00211 
00212       DSET_unload_one( UC_dset , mm ) ;
00213 
00214       if( UC_mask == NULL ){
00215          for( kk=0 ; kk < nxyz ; kk++ ) UC_vec[kk][mm] = bb[kk] ;
00216       } else {
00217          for( nn=kk=0 ; kk < nxyz ; kk++ )
00218             if( UC_mask[kk] ) UC_vec[nn++][mm] = bb[kk] ;
00219       }
00220    }
00221    free(bb) ; DSET_unload( UC_dset ) ;
00222 
00223    /* detrend and normalize vectors */
00224 
00225    if( !UC_be_quiet ) printf("--- normalizing vectors\n") ;
00226 
00227    for( kk=0 ; kk < UC_nvec ; kk++ )
00228       normalize( UC_vdim , UC_vec[kk] ) ;
00229 
00230    for( kk=0 ; kk < UC_ref->ny ; kk++ )
00231       normalize( UC_vdim , MRI_FLOAT_PTR(UC_ref) + kk*UC_ref->nx ) ;
00232 
00233    return ;
00234 }
00235 
00236 #include "uuu.c"
00237 
00238 /*-----------------------------------------------------------------------
00239   Compute the unusuality index of a reference vector in a set of vectors
00240 -------------------------------------------------------------------------*/
00241 
00242 float UC_unusuality( int ndim , float * ref , int nvec , float ** vec )
00243 {
00244    register int ii , kk ;
00245    float psum,msum , * vv , val ;
00246    float zmid , zmed , zsig , zplus,zminus , uval ;
00247 
00248    static int     nvold=-1   ;
00249    static float * zval =NULL , *aval=NULL ;
00250    static float   pstar , zstar ;
00251 
00252    if( ndim < 4 || nvec < 4 || ref == NULL || vec == NULL ) return 0.0 ;
00253 
00254    /* initialize if number of vectors has changed */
00255 
00256    if( nvold != nvec ){
00257       if( zval != NULL ) free(zval) ;
00258       if( aval != NULL ) free(aval) ;
00259       zval = (float *) malloc(sizeof(float)*nvec) ;
00260       aval = (float *) malloc(sizeof(float)*nvec) ;
00261       nvold = nvec ;
00262       pstar = 10.0 / nvec ;
00263       zstar = qginv(0.5*pstar) ;
00264    }
00265 
00266    /* compute dot products */
00267 
00268    for( kk=0 ; kk < nvec ; kk++ ){
00269       psum = 0.0 ; vv = vec[kk] ;
00270       for( ii=0 ; ii < ndim ; ii++ ) psum += ref[ii] * vv[ii] ;
00271       zval[kk] = psum ;
00272    }
00273 
00274    { char * ps = getenv("PTAIL") ;
00275      float pp=0.0 ;
00276      if( ps != NULL ) pp = strtod(ps,NULL) ;
00277      set_unusuality_tail(pp) ;
00278    }
00279 
00280    psum = unusuality( nvec, zval ) ;
00281 
00282    for( kk=0 ; kk < nvec ; kk++ ) zval[kk] = -zval[kk] ;
00283 
00284    msum = unusuality( nvec, zval ) ;
00285 
00286    uval = psum - msum ;
00287 
00288    printf("psum=%.1f  msum=%.1f  total=%.1f  uval=%.1f\n",
00289           psum,msum,psum+msum,uval) ;
00290    return (psum+msum) ;
00291 }
00292 
00293 /*---------------------------------------------------------------------------*/
00294 
00295 void UC_syntax(char * msg)
00296 {
00297    if( msg != NULL ){ fprintf(stderr,"\n*** %s\n",msg) ; exit(1) ; }
00298 
00299    printf(
00300     "Unusual Component Analysis of 3D Datasets\n"
00301     "Usage: 3duca [options] dataset ...\n"
00302     "\n"
00303     "The input dataset may have a sub-brick selector list.\n"
00304     "Otherwise, all sub-bricks from a dataset will be used.\n"
00305     "\n"
00306     "OPTIONS:\n"
00307     "  -prefix pname \n"
00308     "  -verbose\n"
00309     "  -mask mset \n"
00310     "  -ref file.1D\n"
00311     printf("\n" MASTER_SHORTHELP_STRING ) ;
00312    ) ;
00313 
00314    exit(0) ;
00315 }
00316 
00317 /*------------------------------------------------------------------*/
00318 
00319 int main( int argc , char * argv[] )
00320 {
00321    int kk ;
00322 
00323    /*-- read command line arguments --*/
00324 
00325    if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ) UC_syntax(NULL) ;
00326 
00327    (void) my_getenv("junk") ;
00328 
00329    UC_read_opts( argc , argv ) ;
00330 
00331    for( kk=0 ; kk < UC_ref->ny ; kk++ )
00332       (void) UC_unusuality( UC_vdim,
00333                             MRI_FLOAT_PTR(UC_ref) + kk*UC_ref->nx,
00334                             UC_nvec, UC_vec ) ;
00335    exit(0) ;
00336 }