00001 #include "mrilib.h"
00002 
00003 /************************************************************************/
00004 /******* Hack to remove large spikes from time series data (oog). *******/
00005 /******* 30 Aug 2002 - RWCox                                      *******/
00006 /************************************************************************/
00007 
00008 #define TFAC  0.1    /* scale factor for -ssave dataset */
00009 #define ITFAC 10.0   /* inverse of above */
00010 
00011 /*----------------------------------------------------------------------*/
00012 
00013 static INLINE float mytanh( float x )
00014 {
00015   register float ex , exi ;
00016        if( x >  7.0 ) return  1.0 ;  /* 03 Sep: check for stupid inputs */
00017   else if( x < -7.0 ) return -1.0 ;
00018   ex = exp(x) ; exi = 1.0/ex ;
00019   return (ex-exi)/(ex+exi) ;
00020 }
00021 
00022 /*----------------------------------------------------------------------*/
00023 
00024 int main( int argc , char * argv[] )
00025 {
00026    THD_3dim_dataset *dset , *oset=NULL , *tset=NULL ;
00027    int nvals , iv , nxyz , ii,jj,kk , iarg , kz,kzold ;
00028    float cut1=2.5,cut2=4.0 , sq2p , fq , val , fsig ;
00029    MRI_IMAGE *flim ;
00030    float *far, *dar , *var , *fitar ;
00031    char *prefix="despike" , *tprefix=NULL ;
00032 
00033    int corder=-1 , nref , ignore=0 , polort=2 , nuse , nomask=0 ;
00034    int nspike, nbig, nproc ;
00035    float **ref ;
00036    float  *fit , *ssp , snew , c21,ic21 , pspike,pbig , cls ;
00037    short  *sar , *qar ;
00038    byte   *tar , *mask=NULL ;
00039    float  *zar , *yar ;
00040    int     datum ;
00041 
00042    /*----- Read command line -----*/
00043 
00044    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00045       printf("Usage: 3dDespike [options] dataset\n"
00046              "Removes 'spikes' from the 3D+time input dataset and writes\n"
00047              "a new dataset with the spike values replaced by something\n"
00048              "more pleasing.\n"
00049              "\n"
00050              "Method:\n"
00051              " * L1 fit a smooth-ish curve to each voxel time series\n"
00052              "    [see -corder option for description of the curve].\n"
00053              " * Compute the MAD of the difference between the curve and\n"
00054              "    the data time series (the residuals).\n"
00055              " * Estimate the standard deviation 'sigma' of the residuals\n"
00056              "    as sqrt(PI/2)*MAD.\n"
00057              " * For each voxel value, define s = (value-curve)/sigma.\n"
00058              " * Values with s > c1 are replaced with a value that yields\n"
00059              "    a modified s' = c1+(c2-c1)*tanh((s-c1)/(c2-c1)).\n"
00060              " * c1 is the threshold value of s for a 'spike' [default c1=2.5].\n"
00061              " * c2 is the upper range of the allowed deviation from the curve:\n"
00062              "    s=[c1..infinity) is mapped to s'=[c1..c2)   [default c2=4].\n"
00063              "\n"
00064              "Options:\n"
00065              " -ignore I  = Ignore the first I points in the time series:\n"
00066              "               these values will just be copied to the\n"
00067              "               output dataset [default I=0].\n"
00068              " -corder L  = Set the curve fit order to L:\n"
00069              "               the curve that is fit to voxel data v(t) is\n"
00070              "\n"
00071              "                       k=L [        (2*PI*k*t)          (2*PI*k*t) ]\n"
00072              " f(t) = a+b*t+c*t*t + SUM  [ d * sin(--------) + e * cos(--------) ]\n"
00073              "                       k=1 [  k     (    T   )    k     (    T   ) ]\n"
00074              "\n"
00075              "               where T = duration of time series;\n"
00076              "               the a,b,c,d,e parameters are chosen to minimize\n"
00077              "               the sum over t of |v(t)-f(t)| (L1 regression);\n"
00078              "               this type of fitting is is insensitive to large\n"
00079              "               spikes in the data.  The default value of L is\n"
00080              "               NT/30, where NT = number of time points.\n"
00081              "\n"
00082              " -cut c1 c2 = Alter default values for the spike cut values\n"
00083              "               [default c1=2.5, c2=4.0].\n"
00084              " -prefix pp = Save de-spiked dataset with prefix 'pp'\n"
00085              "               [default pp='despike']\n"
00086              " -ssave ttt = Save 'spikiness' measure s for each voxel into a\n"
00087              "               3D+time dataset with prefix 'ttt' [default=no save]\n"
00088              " -nomask    = Process all voxels\n"
00089              "               [default=use a mask of high-intensity voxels, ]\n"
00090              "               [as created via '3dAutomask -dilate 4 dataset'].\n"
00091              "\n"
00092              "Caveats:\n"
00093              "* Despiking may interfere with image registration, since head\n"
00094              "   movement may produce 'spikes' at the edge of the brain, and\n"
00095              "   this information would be used in the registration process.\n"
00096              "   This possibility has not been explored.\n"
00097              "* Check your data visually before and after despiking and\n"
00098              "   registration!\n"
00099              "   [Hint: open 2 AFNI controllers, and turn Time Lock on.]\n"
00100             ) ;
00101       exit(0) ;
00102    }
00103 
00104    /** AFNI package setup and logging **/
00105 
00106    mainENTRY("3dDespike main"); machdep(); AFNI_logger("3dDespike",argc,argv);
00107    PRINT_VERSION("3dDespike") ;
00108 
00109    /** parse options **/
00110 
00111    iarg = 1 ;
00112    while( iarg < argc && argv[iarg][0] == '-' ){
00113 
00114       /** don't use masking **/
00115 
00116       if( strcmp(argv[iarg],"-nomask") == 0 ){
00117         nomask = 1 ;
00118         iarg++ ; continue ;
00119       }
00120 
00121       /** output dataset prefix **/
00122 
00123       if( strcmp(argv[iarg],"-prefix") == 0 ){
00124         prefix = argv[++iarg] ;
00125         if( !THD_filename_ok(prefix) ){
00126           fprintf(stderr,"** ERROR: -prefix is not good!\n"); exit(1);
00127         }
00128         iarg++ ; continue ;
00129       }
00130 
00131       /** ratio dataset prefix **/
00132 
00133       if( strcmp(argv[iarg],"-ssave") == 0 ){
00134         tprefix = argv[++iarg] ;
00135         if( !THD_filename_ok(tprefix) ){
00136           fprintf(stderr,"** ERROR: -ssave prefix is not good!\n"); exit(1);
00137         }
00138         iarg++ ; continue ;
00139       }
00140 
00141       /** trigonometric polynomial order **/
00142 
00143       if( strcmp(argv[iarg],"-corder") == 0 ){
00144         corder = strtol( argv[++iarg] , NULL , 10 ) ;
00145         if( corder < 0 ){
00146           fprintf(stderr,"** Illegal value of -corder!\n"); exit(1);
00147         }
00148         iarg++ ; continue ;
00149       }
00150 
00151       /** how much to ignore at start **/
00152 
00153       if( strcmp(argv[iarg],"-ignore") == 0 ){
00154         ignore = strtol( argv[++iarg] , NULL , 10 ) ;
00155         if( ignore < 0 ){
00156           fprintf(stderr,"** Illegal value of -ignore!\n"); exit(1);
00157         }
00158         iarg++ ; continue ;
00159       }
00160 
00161       /** thresholds for s ratio **/
00162 
00163       if( strcmp(argv[iarg],"-cut") == 0 ){
00164         cut1 = strtod( argv[++iarg] , NULL ) ;
00165         cut2 = strtod( argv[++iarg] , NULL ) ;
00166         if( cut1 < 1.0 || cut2 < cut1+0.5 ){
00167           fprintf(stderr,"** Illegal values after -cut!\n"); exit(1);
00168         }
00169         iarg++ ; continue ;
00170       }
00171 
00172       fprintf(stderr,"** Unknown option: %s\n",argv[iarg]) ; exit(1) ;
00173    }
00174 
00175    c21 = cut2-cut1 ; ic21 = 1.0/c21 ;
00176 
00177    /*----- read input dataset -----*/
00178 
00179    if( iarg >= argc ){
00180      fprintf(stderr,"** No input dataset!?\n"); exit(1);
00181    }
00182 
00183    dset = THD_open_dataset( argv[iarg] ) ;
00184    if( !ISVALID_DSET(dset) ){
00185      fprintf(stderr,"** Can't open dataset %s\n",argv[iarg]); exit(1);
00186    }
00187    datum = DSET_BRICK_TYPE(dset,0) ;
00188    if( (datum != MRI_short && datum != MRI_float) || !DSET_datum_constant(dset) ){
00189      fprintf(stderr,"** Can't process non-short, non-float dataset!\n") ; exit(1) ;
00190    }
00191    fprintf(stderr,"Input data type = %s\n",MRI_TYPE_name[datum]) ;
00192    nvals = DSET_NUM_TIMES(dset) ; nuse = nvals - ignore ;
00193    if( nuse < 15 ){
00194      fprintf(stderr,"** Can't use dataset with < 15 time points per voxel!\n") ;
00195      exit(1) ;
00196    }
00197    fprintf(stderr,"++ ignoring first %d time points, using last %d\n",ignore,nuse);
00198    if( corder > 0 && 4*corder+2 > nuse ){
00199      fprintf(stderr,"** -corder %d is too big for NT=%d!\n",corder,nvals) ;
00200      exit(1) ;
00201    } else if( corder < 0 ){
00202      corder = rint(nuse/30.0) ;
00203      fprintf(stderr,"++ using %d time points => -corder %d\n",nuse,corder) ;
00204    } else {
00205      fprintf(stderr,"++ -corder %d set from command line\n",corder) ;
00206    }
00207    nxyz = DSET_NVOX(dset) ;
00208    fprintf(stderr,"++ Loading dataset %s\n",argv[iarg]) ;
00209    DSET_load(dset) ;
00210    if( !DSET_LOADED(dset) ){
00211      fprintf(stderr,"** Can't read dataset from disk!\n") ; exit(1) ;
00212    }
00213 
00214    /*-- create automask --*/
00215 
00216    if( !nomask ){
00217      mask = THD_automask( dset ) ;
00218      for( ii=0 ; ii < 4 ; ii++ )
00219        THD_mask_dilate( DSET_NX(dset), DSET_NY(dset), DSET_NZ(dset), mask, 3 ) ;
00220      ii = THD_countmask( DSET_NVOX(dset) , mask ) ;
00221      fprintf(stderr,"++ %d voxels in the mask [out of %d in dataset]\n",ii,DSET_NVOX(dset)) ;
00222    } else {
00223      fprintf(stderr,"++ processing all %d voxels in dataset\n",DSET_NVOX(dset)) ;
00224    }
00225 
00226    /*-- create empty despiked dataset --*/
00227 
00228    oset = EDIT_empty_copy( dset ) ;
00229    EDIT_dset_items( oset ,
00230                       ADN_prefix    , prefix ,
00231                       ADN_brick_fac , NULL ,
00232                       ADN_datum_all , datum ,
00233                     ADN_none ) ;
00234 
00235    if( THD_is_file(DSET_HEADNAME(oset)) ){
00236      fprintf(stderr,"** ERROR: output dataset already exists: %s\n",DSET_HEADNAME(oset));
00237      exit(1);
00238    }
00239 
00240    tross_Copy_History( oset , dset ) ;
00241    tross_Make_History( "3dDespike" , argc , argv , oset ) ;
00242 
00243    /* create bricks (will be filled with zeros) */
00244 
00245    for( iv=0 ; iv < nvals ; iv++ )
00246      EDIT_substitute_brick( oset , iv , datum , NULL ) ;
00247 
00248    /* copy the ignored bricks */
00249 
00250    switch( datum ){
00251      case MRI_short:
00252        for( iv=0 ; iv < ignore ; iv++ ){
00253          sar = DSET_ARRAY(oset,iv) ;
00254          qar = DSET_ARRAY(dset,iv) ;
00255          memcpy( sar , qar , DSET_BRICK_BYTES(dset,iv) ) ;
00256          DSET_unload_one(dset,iv) ;
00257        }
00258      break ;
00259      case MRI_float:
00260        for( iv=0 ; iv < ignore ; iv++ ){
00261          zar = DSET_ARRAY(oset,iv) ;
00262          yar = DSET_ARRAY(dset,iv) ;
00263          memcpy( zar , yar , DSET_BRICK_BYTES(dset,iv) ) ;
00264          DSET_unload_one(dset,iv) ;
00265        }
00266      break ;
00267    }
00268 
00269    /*-- setup to save a threshold statistic dataset, if desired --*/
00270 
00271    if( tprefix != NULL ){
00272      float *fac ;
00273      tset = EDIT_empty_copy( dset ) ;
00274      fac  = (float *) malloc( sizeof(float) * nvals ) ;
00275      for( ii=0 ; ii < nvals ; ii++ ) fac[ii] = TFAC ;
00276      EDIT_dset_items( tset ,
00277                         ADN_prefix    , tprefix ,
00278                         ADN_brick_fac , fac ,
00279                         ADN_datum_all , MRI_byte ,
00280                         ADN_func_type , FUNC_FIM_TYPE ,
00281                       ADN_none ) ;
00282      free(fac) ;
00283 
00284      tross_Copy_History( tset , dset ) ;
00285      tross_Make_History( "3dDespike" , argc , argv , tset ) ;
00286 
00287      if( THD_is_file(DSET_HEADNAME(tset)) ){
00288        fprintf(stderr,"** ERROR: -ssave dataset already exists!\n"); exit(1);
00289      }
00290 
00291      tross_Copy_History( tset , dset ) ;
00292      tross_Make_History( "3dDespike" , argc , argv , tset ) ;
00293 
00294      for( iv=0 ; iv < nvals ; iv++ )
00295        EDIT_substitute_brick( tset , iv , MRI_byte , NULL ) ;
00296    }
00297 
00298    /*-- setup to find spikes --*/
00299 
00300    sq2p  = sqrt(0.5*PI) ;
00301    var   = (float *) malloc( sizeof(float) * nvals ) ;
00302    far   = (float *) malloc( sizeof(float) * nvals ) ;
00303    dar   = (float *) malloc( sizeof(float) * nvals ) ;
00304    fitar = (float *) malloc( sizeof(float) * nvals ) ;
00305    ssp   = (float *) malloc( sizeof(float) * nvals ) ;
00306 
00307    /* make ref functions */
00308 
00309    nref = 2*corder+3 ;
00310    ref  = (float **) malloc( sizeof(float *) * nref ) ;
00311    for( jj=0 ; jj < nref ; jj++ )
00312      ref[jj] = (float *) malloc( sizeof(float) * nuse ) ;
00313 
00314    fit = (float *) malloc( sizeof(float) * nref ) ;
00315 
00316    /* r(t) = 1 */
00317 
00318    for( iv=0 ; iv < nuse ; iv++ ) ref[0][iv] = 1.0 ;
00319    jj = 1 ;
00320 
00321    /* r(t) = t - tmid */
00322 
00323    { float tm = 0.5 * (nuse-1.0) ; float fac = 2.0 / nuse ;
00324      for( iv=0 ; iv < nuse ; iv++ ) ref[1][iv] = (iv-tm)*fac ;
00325      jj = 2 ;
00326 
00327      /* r(t) = (t-tmid)**jj */
00328 
00329      for( ; jj <= polort ; jj++ )
00330        for( iv=0 ; iv < nuse ; iv++ )
00331          ref[jj][iv] = pow( (iv-tm)*fac , (double)jj ) ;
00332    }
00333 
00334    for( kk=1 ; kk <= corder ; kk++ ){
00335      fq = (2.0*PI*kk)/nuse ;
00336 
00337      /* r(t) = sin(2*PI*k*t/N) */
00338 
00339      for( iv=0 ; iv < nuse ; iv++ )
00340        ref[jj][iv] = sin(fq*iv) ;
00341      jj++ ;
00342 
00343      /* r(t) = cos(2*PI*k*t/N) */
00344 
00345      for( iv=0 ; iv < nuse ; iv++ )
00346        ref[jj][iv] = cos(fq*iv) ;
00347      jj++ ;
00348    }
00349 
00350    /*--- loop over voxels and do work ---*/
00351 
00352    fprintf(stderr,"++ edit thresholds: %.1f .. %.1f standard deviations\n",cut1,cut2) ;
00353    fprintf(stderr,"++                [ %.4f%% .. %.4f%% of normal distribution]\n",
00354                   200.0*qg(cut1) , 200.0*qg(cut2) ) ;
00355    fprintf(stderr,"++ %d slices to process\n",DSET_NZ(dset)) ;
00356    kzold  = -1 ;
00357    nspike =  0 ; nbig = 0 ; nproc = 0 ;
00358    for( ii=0 ; ii < nxyz ; ii++ ){   /* ii = voxel index */
00359 
00360       if( mask != NULL && mask[ii] == 0 ) continue ;   /* skip this voxel */
00361 
00362       kz = DSET_index_to_kz(dset,ii) ;       /* starting a new slice */
00363       if( kz != kzold ){
00364         fprintf(stderr, "++ start slice %2d",kz ) ;
00365         if( nproc > 0 ){
00366           pspike = (100.0*nspike)/nproc ;
00367           pbig   = (100.0*nbig  )/nproc ;
00368           fprintf(stderr,"; done %d data points, %d edits [%.3f%%], %d big edits [%.3f%%]",
00369                   nproc,nspike,pspike,nbig,pbig ) ;
00370         }
00371         fprintf(stderr,"\n") ;
00372         kzold = kz ;
00373       }
00374 
00375       /*** extract ii-th time series into far[] ***/
00376 
00377       switch( datum ){
00378         case MRI_short:
00379           for( iv=0 ; iv < nuse ; iv++ ){
00380             qar = DSET_ARRAY(dset,iv+ignore) ;   /* skip ignored data */
00381             far[iv] = (float) qar[ii] ;
00382           }
00383         break ;
00384         case MRI_float:
00385           for( iv=0 ; iv < nuse ; iv++ ){
00386             zar = DSET_ARRAY(dset,iv+ignore) ;
00387             far[iv] = zar[ii] ;
00388           }
00389         break ;
00390       }
00391       memcpy(dar,far,sizeof(float)*nuse) ;   /* copy time series into dar[] */
00392 
00393       /*** solve for L1 fit ***/
00394 
00395       cls = cl1_solve( nuse , nref , far , ref , fit,0 );
00396 
00397       if( cls < 0.0 ){                       /* fit failed! */
00398 #if 0
00399         fprintf(stderr,"curve fit fails at voxel %d %d %d\n",
00400                 DSET_index_to_ix(dset,ii) ,
00401                 DSET_index_to_jy(dset,ii) ,
00402                 DSET_index_to_kz(dset,ii)  ) ;
00403 #endif
00404         continue ;                           /* skip this voxel */
00405       }
00406 
00407       for( iv=0 ; iv < nuse ; iv++ ){        /* detrend */
00408         val =  fit[0]
00409              + fit[1]*ref[1][iv]             /* quadratic part of curve fit */
00410              + fit[2]*ref[2][iv] ;
00411         for( jj=3 ; jj < nref ; jj++ )       /* rest of curve fit */
00412           val += fit[jj] * ref[jj][iv] ;
00413 
00414         fitar[iv] = val ;                    /* save curve value */
00415         var[iv]   = dar[iv]-val ;            /* remove fitted value = resid */
00416         far[iv]   = fabs(var[iv]) ;          /* abs value of resid */
00417       }
00418 
00419       /*** compute estimate standard deviation of detrended data ***/
00420 
00421       fsig = sq2p * qmed_float(nuse,far) ;   /* also mangles far array */
00422 
00423       /*** process time series for spikes, editing data in dar[] ***/
00424 
00425       if( fsig > 0.0 ){                      /* data wasn't fit perfectly */
00426 
00427         /* find spikiness for each point in time */
00428 
00429         fq = 1.0 / fsig ;
00430         for( iv=0 ; iv < nuse ; iv++ ){
00431           ssp[iv] = fq * var[iv] ;           /* spikiness s = how many sigma out */
00432         }
00433 
00434         /* save spikiness in -ssave datset */
00435 
00436         if( tset != NULL ){
00437           for( iv=0 ; iv < nuse ; iv++ ){
00438             tar   = DSET_ARRAY(tset,iv+ignore) ;
00439             snew   = ITFAC*fabs(ssp[iv]) ;   /* scale for byte storage */
00440             tar[ii] = BYTEIZE(snew) ;        /* cf. mrilib.h */
00441           }
00442         }
00443 
00444         /* process values of |s| > cut1, editing dar[] */
00445 
00446         for( iv=0 ; iv < nuse ; iv++ ){
00447           if( ssp[iv] > cut1 ){
00448             snew = cut1 + c21*mytanh((ssp[iv]-cut1)*ic21) ;   /* edit s down */
00449             dar[iv] = fitar[iv] + snew*fsig ;
00450             nspike++ ; if( ssp[iv] > cut2 ) nbig++ ;          /* count edits */
00451           } else if( ssp[iv] < -cut1 ){
00452             snew = -cut1 + c21*mytanh((ssp[iv]+cut1)*ic21) ;  /* edit s up */
00453             dar[iv] = fitar[iv] + snew*fsig ;
00454             nspike++ ; if( ssp[iv] < -cut2 ) nbig++ ;
00455           }
00456         }
00457         nproc += nuse ;  /* number data points processed */
00458 
00459       } /* end of processing time series when fsig is positive */
00460 
00461       /* put dar[] time series (possibly edited above) into output bricks */
00462 
00463       switch( datum ){
00464         case MRI_short:
00465           for( iv=0 ; iv < nuse ; iv++ ){
00466             sar = DSET_ARRAY(oset,iv+ignore) ; /* output brick */
00467             sar[ii] = (short) dar[iv] ;        /* original or mutated data */
00468           }
00469         break ;
00470         case MRI_float:
00471           for( iv=0 ; iv < nuse ; iv++ ){
00472             zar = DSET_ARRAY(oset,iv+ignore) ; /* output brick */
00473             zar[ii] = dar[iv] ;                /* original or mutated data */
00474           }
00475         break ;
00476       }
00477 
00478    } /* end of loop over voxels #ii */
00479 
00480    /*--- finish up ---*/
00481 
00482    DSET_delete(dset) ; /* delete input dataset */
00483 
00484    if( nproc > 0 ){
00485      pspike = (100.0*nspike)/nproc ;
00486      pbig   = (100.0*nbig  )/nproc ;
00487      fprintf(stderr,"++ FINAL: %d data points, %d edits [%.3f%%], %d big edits [%.3f%%]\n",
00488              nproc,nspike,pspike,nbig,pbig ) ;
00489    } else {
00490      fprintf(stderr,"++ FINAL: no good voxels found to process!!!!\n") ;
00491    }
00492 
00493    /* write results */
00494 
00495    DSET_write(oset) ;
00496    fprintf(stderr,"++ Wrote output dataset %s\n",DSET_BRIKNAME(oset)) ;
00497    DSET_delete(oset) ;
00498 
00499    if( tset != NULL ){
00500      DSET_write(tset) ;
00501      fprintf(stderr,"++ Wrote -ssave dataset %s\n",DSET_BRIKNAME(tset)) ;
00502    }
00503 
00504    exit(0) ;
00505 }