00603 {
00604    float * wt , tau , t , sum ;
00605    int ii , nt2=(ntap-1)/2 ;
00606 
00607    if( ntap < 3 ) return NULL ;
00608 
00609    wt = (float *) calloc(sizeof(float),ntap) ;
00610    tau = nt2 + 1.0 ;
00611 
00612    for( sum=0.0,ii=0 ; ii <= 2*nt2 ; ii++ ){
00613       t = PI*(ii-nt2)/tau ;
00614       wt[ii] = 0.42 + 0.5*cos(t) + 0.08*cos(2*t) ;
00615       sum += wt[ii] ;
00616    }
00617    sum = 1.0 / sum ;
00618    for( ii=0 ; ii < ntap ; ii++ ) wt[ii] *= sum ;
00619 
00620    return wt ;
00621 }

00560 {
00561    MRI_IMAGE * filter_data=NULL;
00562    float * filter_coefficients = NULL;
00563 
00564    filter_data = mri_read_1D(custom_file);
00565    if MRI_IS_1D(filter_data) {
00566      custom_ntaps = filter_data->nx;
00567      if (custom_ntaps % 2) {
00568        filter_coefficients = MRI_FLOAT_PTR(filter_data);
00569      } else {
00570        /*ERROR, EVEN NUMBER OF FILTER TAPS */
00571      }
00572    }
00573 
00574    return filter_coefficients ;
00575 }

00580 {
00581    float * wt , tau , t , sum ;
00582    int ii , nt2=(ntap-1)/2 ;
00583 
00584    if( ntap < 3 ) return NULL ;
00585 
00586    wt = (float *) calloc(sizeof(float),ntap) ;
00587    tau = nt2 + 1.0 ;
00588 
00589    for( sum=0.0,ii=0 ; ii <= 2*nt2 ; ii++ ){
00590       t = PI*(ii-nt2)/tau ;
00591       wt[ii] = 0.54 + 0.46*cos(t) ;
00592       sum += wt[ii] ;
00593    }
00594    sum = 1.0 / sum ;
00595    for( ii=0 ; ii < ntap ; ii++ ) wt[ii] *= sum ;
00596 
00597    return wt ;
00598 }

00459 {
00460    int ii ;
00461    float aa,bb,cc ;
00462 
00463    bb = vec[0] ; cc = vec[1] ; vec[0] = LSUM(bb,bb,cc) ;
00464    for( ii=1 ; ii < num-1 ; ii++ ){
00465       aa = bb ; bb = cc ; cc = vec[ii+1] ;
00466       vec[ii] = LSUM(aa,bb,cc) ;
00467    }
00468    vec[num-1] = LSUM(bb,cc,cc) ;
00469 
00470    return ;
00471 }

00476 {
00477    int ii , nt2=(ntap-1)/2 , jj ;
00478    float sum ;
00479    static int nfar=0 ;
00480    static float *far=NULL ;
00481 
00482    if( npt > nfar ){
00483       if(far != NULL) free(far) ;
00484       far = (float *)malloc(sizeof(float)*npt) ; nfar = npt ;
00485    }
00486 
00487 #undef XX
00488 #define XX(i) ( ((i)<0) ? far[0] : ((i)>npt-1) ? far[npt-1] : far[i] )
00489 
00490    memcpy( far , x , sizeof(float)*npt ) ;
00491 
00492    for( ii=0 ; ii < npt ; ii++ ){
00493       for( sum=0.0,jj=0 ; jj < ntap ; jj++ ) sum += wt[jj] * XX(ii-nt2+jj) ;
00494       x[ii] = sum ;
00495    }
00496 
00497    return ;
00498 }

00530 {
00531    int ii , nt2=(ntap-1)/2 , jj ;
00532    float sum ;
00533    static int nfar=0 ;
00534    static float *far=NULL ;
00535    float a=0.0,b=0.0 ;
00536 
00537    if( npt > nfar ){
00538       if(far != NULL) free(far) ;
00539       far = (float *)malloc(sizeof(float)*npt) ; nfar = npt ;
00540    }
00541 
00542 #undef XX
00543 #define XX(i) ( ((i)<0 || (i)>npt-1) ? (a+b*(i)) : far[i] )
00544 
00545    memcpy( far , x , sizeof(float)*npt ) ;
00546 
00547    get_linear_trend( npt,far , &a,&b ) ; /* cf. thd_detrend.c */
00548 
00549    for( ii=0 ; ii < npt ; ii++ ){
00550       for( sum=0.0,jj=0 ; jj < ntap ; jj++ ) sum += wt[jj] * XX(ii-nt2+jj) ;
00551       x[ii] = sum ;
00552    }
00553 
00554    return ;
00555 }

00503 {
00504    int ii , nt2=(ntap-1)/2 , jj ;
00505    float sum ;
00506    static int nfar=0 ;
00507    static float *far=NULL ;
00508 
00509    if( npt > nfar ){
00510       if(far != NULL) free(far) ;
00511       far = (float *)malloc(sizeof(float)*npt) ; nfar = npt ;
00512    }
00513 
00514 #undef XX
00515 #define XX(i) ( ((i)<0 || (i)>npt-1) ? 0 : far[i] )
00516 
00517    memcpy( far , x , sizeof(float)*npt ) ;
00518 
00519    for( ii=0 ; ii < npt ; ii++ ){
00520       for( sum=0.0,jj=0 ; jj < ntap ; jj++ ) sum += wt[jj] * XX(ii-nt2+jj) ;
00521       x[ii] = sum ;
00522    }
00523 
00524    return ;
00525 }

00027 {
00028    void (*smth)(int,float *) = linear3_func ;     /* default filter */
00029 
00030    char prefix[256] = "smooth" ;
00031    int new_datum = ILLEGAL_TYPE , old_datum ;
00032    int nopt ;
00033    THD_3dim_dataset * old_dset , * new_dset ;
00034    int ii,kk , nxyz , ntime , use_fac , ityp , nbytes ;
00035    void * new_brick ;
00036 
00037    byte  ** bptr = NULL ;  /* one of these will be the array of */
00038    short ** sptr = NULL ;  /* pointers to input dataset sub-bricks */
00039    float ** fptr = NULL ;  /* (depending on input datum type) */
00040 
00041    byte  ** new_bptr = NULL ;  /* one of these will be the array of */
00042    short ** new_sptr = NULL ;  /* pointers to output dataset sub-bricks */
00043    float ** new_fptr = NULL ;  /* (depending on output datum type) */
00044 
00045    float * fxar = NULL ;  /* array loaded from input dataset */
00046    float * fac  = NULL ;  /* array of brick scaling factors */
00047    float * faci = NULL ;
00048 
00049 #define BLACKMAN 1
00050 #define HAMMING  2
00051 #define CUSTOM   3
00052 
00053 #define EXTEND   77
00054 #define ZERO     78
00055 #define TREND    79
00056 
00057    int ntap=0 ;      /* 01 Mar 2001 */
00058    float *ftap=NULL ;
00059    int nwin=0,nfil=EXTEND ;      /* 03 Mar 2001 */
00060 
00061    void (*lfil)(int,float *,int,float *) = linear_filter_extend ;
00062    float * (*lwin)(int) = NULL ;
00063 
00064    /* start of code */
00065 
00066    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00067       printf("Usage: 3dTsmooth [options] dataset\n"
00068              "Smooths each voxel time series in a 3D+time dataset and produces\n"
00069              "as output a new 3D+time dataset (e.g., lowpass filter in time).\n"
00070              "\n"
00071              "General Options:\n"
00072              "  -prefix ppp  = Sets the prefix of the output dataset to be 'ppp'.\n"
00073              "                   [default = 'smooth']\n"
00074              "  -datum type  = Coerce output dataset to be stored as the given type.\n"
00075              "                   [default = input data type]\n"
00076              "\n"
00077              "Three Point Filtering Options [07 July 1999]\n"
00078              "--------------------------------------------\n"
00079              "The following options define the smoothing filter to be used.\n"
00080              "All these filters  use 3 input points to compute one output point:\n"
00081              "  Let a = input value before the current point\n"
00082              "      b = input value at the current point\n"
00083              "      c = input value after the current point\n"
00084              "           [at the left end, a=b; at the right end, c=b]\n"
00085              "\n"
00086              "  -lin = 3 point linear filter: 0.15*a + 0.70*b + 0.15*c\n"
00087              "           [This is the default smoother]\n"
00088              "  -med = 3 point median filter: median(a,b,c)\n"
00089              "  -osf = 3 point order statistics filter:\n"
00090              "           0.15*min(a,b,c) + 0.70*median(a,b,c) + 0.15*max(a,b,c)\n"
00091              "\n"
00092              "  -3lin m = 3 point linear filter: 0.5*(1-m)*a + m*b + 0.5*(1-m)*c\n"
00093              "              Here, 'm' is a number strictly between 0 and 1.\n"
00094              "\n"
00095              "General Linear Filtering Options [03 Mar 2001]\n"
00096              "----------------------------------------------\n"
00097              "  -hamming N  = Use N point Hamming or Blackman windows.\n"
00098              "  -blackman N     (N must be odd and bigger than 1.)\n"
00099              "  -custom coeff_filename.1D (odd # of coefficients must be in a \n"
00100              "                             single column in ASCII file)\n"
00101              "   (-custom added Jan 2003)\n"
00102              "    WARNING: If you use long filters, you do NOT want to include the\n"
00103              "             large early images in the program.  Do something like\n"
00104              "                3dTsmooth -hamming 13 'fred+orig[4..$]'\n"
00105              "             to eliminate the first 4 images (say).\n"
00106              " The following options determing how the general filters treat\n"
00107              " time points before the beginning and after the end:\n"
00108              "  -EXTEND = BEFORE: use the first value; AFTER: use the last value\n"
00109              "  -ZERO   = BEFORE and AFTER: use zero\n"
00110              "  -TREND  = compute a linear trend, and extrapolate BEFORE and AFTER\n"
00111              " The default is -EXTEND.  These options do NOT affect the operation\n"
00112              " of the 3 point filters described above, which always use -EXTEND.\n"
00113            ) ;
00114       printf("\n" MASTER_SHORTHELP_STRING ) ;
00115       exit(0) ;
00116    }
00117 
00118    mainENTRY("3dTsmooth main"); machdep(); AFNI_logger("3dTsmooth",argc,argv);
00119    PRINT_VERSION("3dTsmooth") ;
00120 
00121    /* parse options */
00122 
00123    nopt = 1 ;
00124 
00125    while( nopt < argc && argv[nopt][0] == '-' ){
00126 
00127       if( strcmp(argv[nopt],"-EXTEND") == 0 ){         /* 03 Mar 2001 */
00128          nfil = EXTEND ; lfil = linear_filter_extend ;
00129          nopt++ ; continue ;
00130       }
00131 
00132       if( strcmp(argv[nopt],"-ZERO") == 0 ){           /* 03 Mar 2001 */
00133          nfil = ZERO ; lfil = linear_filter_zero ;
00134          nopt++ ; continue ;
00135       }
00136 
00137       if( strcmp(argv[nopt],"-TREND") == 0 ){          /* 03 Mar 2001 */
00138          nfil = TREND ; lfil = linear_filter_trend ;
00139          nopt++ ; continue ;
00140       }
00141 
00142       if( strcmp(argv[nopt],"-hamming") == 0 ){
00143          if( ++nopt >= argc ){fprintf(stderr,"*** Illegal -hamming!\n");exit(1);}
00144          ntap = (int) strtod(argv[nopt],NULL) ;
00145          if( ntap < 3 || ntap%2 != 1 ){fprintf(stderr,"*** Illegal -hamming!\n");exit(1);}
00146          nwin = HAMMING ; lwin = hamming_window ;
00147          nopt++ ; continue ;
00148       }
00149 
00150       if( strcmp(argv[nopt],"-blackman") == 0 ){
00151          if( ++nopt >= argc ){fprintf(stderr,"*** Illegal -blackman!\n");exit(1);}
00152          ntap = (int) strtod(argv[nopt],NULL) ;
00153          if( ntap < 3 || ntap%2 != 1 ){fprintf(stderr,"*** Illegal -blackman!\n");exit(1);}
00154          nwin = BLACKMAN ; lwin = blackman_window ;
00155          nopt++ ; continue ;
00156       }
00157 
00158       if( strcmp(argv[nopt],"-custom") == 0 ){
00159          if( ++nopt >= argc ){fprintf(stderr,"*** Illegal -custom!\n");exit(1);}
00160          strcpy(custom_file, argv[nopt]) ;
00161          nwin = CUSTOM ; lwin = custom_filter ;
00162          ntap = 1;
00163          nopt++ ; continue ;
00164       }
00165 
00166       if( strcmp(argv[nopt],"-prefix") == 0 ){
00167          if( ++nopt >= argc ){fprintf(stderr,"*** Illegal -prefix!\n");exit(1);}
00168          strcpy(prefix,argv[nopt]) ;
00169          if( !THD_filename_ok(prefix) ){fprintf(stderr,"*** Illegal -prefix!\n");exit(1);}
00170          nopt++ ; continue ;
00171       }
00172 
00173       if( strcmp(argv[nopt],"-datum") == 0 ){
00174          if( ++nopt >= argc ){fprintf(stderr,"*** Illegal -datum!\n");exit(1);}
00175          if( strcmp(argv[nopt],"short") == 0 ){
00176             new_datum = MRI_short ;
00177          } else if( strcmp(argv[nopt],"float") == 0 ){
00178             new_datum = MRI_float ;
00179          } else if( strcmp(argv[nopt],"byte") == 0 ){
00180             new_datum = MRI_byte ;
00181          } else {
00182             fprintf(stderr,"*** Illegal -datum!\n");exit(1);
00183          }
00184          nopt++ ; continue ;
00185       }
00186 
00187       if( strcmp(argv[nopt],"-lin") == 0 ){
00188          bf = 0.70 ; af = cf = 0.15 ;
00189          smth = linear3_func ;
00190          nopt++ ; continue ;
00191       }
00192 
00193       if( strcmp(argv[nopt],"-med") == 0 ){
00194          smth = median3_func ;
00195          nopt++ ; continue ;
00196       }
00197 
00198       if( strcmp(argv[nopt],"-osf") == 0 ){
00199          smth = osfilt3_func ;
00200          nopt++ ; continue ;
00201       }
00202 
00203       if( strcmp(argv[nopt],"-3lin") == 0 ){
00204          if( ++nopt >= argc ){fprintf(stderr,"*** Illegal -3lin!\n");exit(1);}
00205          bf = strtod( argv[nopt] , NULL ) ;
00206          if( bf <= 0.0 || bf >= 1.0 ){fprintf(stderr,"*** Illegal -3lin!\n");exit(1);}
00207          af = cf = 0.5*(1.0-bf) ;
00208          smth = linear3_func ;
00209          nopt++ ; continue ;
00210       }
00211 
00212       fprintf(stderr,"*** Unknown option: %s\n",argv[nopt]) ; exit(1) ;
00213 
00214    }  /* end of loop over options */
00215 
00216    if( nopt >= argc ){
00217       fprintf(stderr,"*** No input dataset?!\n") ; exit(1) ;
00218    }
00219 
00220    /* open dataset */
00221 
00222    old_dset = THD_open_dataset( argv[nopt] ) ;
00223    if( old_dset == NULL ){
00224       fprintf(stderr,"*** Can't open dataset %s\n",argv[nopt]) ; exit(1) ;
00225    }
00226 
00227    ntime = DSET_NUM_TIMES(old_dset) ;
00228    nxyz  = DSET_NVOX(old_dset) ;
00229 
00230    if( ntime < 4 ){
00231       fprintf(stderr,"*** Can't smooth dataset with less than 4 time points!\n") ;
00232       exit(1) ;
00233    }
00234 
00235    DSET_load(old_dset) ;
00236    if( !DSET_LOADED(old_dset) ){
00237       fprintf(stderr,"*** Can't load datsaet from disk!\n") ; exit(1) ;
00238    }
00239 
00240    old_datum = DSET_BRICK_TYPE(old_dset,0) ;
00241    if( new_datum < 0 ) new_datum = old_datum ;
00242 
00243    switch( old_datum ){  /* pointer type depends on input datum type */
00244 
00245       /** create array of pointers into old dataset sub-bricks **/
00246 
00247       /*--------- input is bytes ----------*/
00248       /* voxel #i at time #k is bptr[k][i] */
00249       /* for i=0..nxyz-1 and k=0..ntime-1. */
00250 
00251       case MRI_byte:
00252          bptr = (byte **) malloc( sizeof(byte *) * ntime ) ;
00253          for( kk=0 ; kk < ntime ; kk++ )
00254             bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk) ;
00255       break ;
00256 
00257       /*--------- input is shorts ---------*/
00258       /* voxel #i at time #k is sptr[k][i] */
00259       /* for i=0..nxyz-1 and k=0..ntime-1. */
00260 
00261       case MRI_short:
00262          sptr = (short **) malloc( sizeof(short *) * ntime ) ;
00263          for( kk=0 ; kk < ntime ; kk++ )
00264             sptr[kk] = (short *) DSET_ARRAY(old_dset,kk) ;
00265       break ;
00266 
00267       /*--------- input is floats ---------*/
00268       /* voxel #i at time #k is fptr[k][i] */
00269       /* for i=0..nxyz-1 and k=0..ntime-1. */
00270 
00271       case MRI_float:
00272          fptr = (float **) malloc( sizeof(float *) * ntime ) ;
00273          for( kk=0 ; kk < ntime ; kk++ )
00274             fptr[kk] = (float *) DSET_ARRAY(old_dset,kk) ;
00275       break ;
00276 
00277    } /* end of switch on input type */
00278 
00279    /*---- allocate space for 1 voxel timeseries ----*/
00280 
00281    fxar = (float *) malloc( sizeof(float) * ntime ) ;   /* voxel timeseries */
00282 
00283    /*--- get scaling factors for sub-bricks ---*/
00284 
00285    fac = (float *) malloc( sizeof(float) * ntime ) ;   /* factors */
00286 
00287    use_fac = 0 ;
00288    for( kk=0 ; kk < ntime ; kk++ ){
00289       fac[kk] = DSET_BRICK_FACTOR(old_dset,kk) ;
00290       if( fac[kk] != 0.0 ) use_fac++ ;
00291       else                 fac[kk] = 1.0 ;
00292    }
00293    if( !use_fac ){
00294       free(fac) ; fac == NULL ;
00295    } else {
00296       faci = (float *) malloc( sizeof(float) * ntime ) ;
00297       for( kk=0 ; kk < ntime ; kk++ ) faci[kk] = 1.0 / fac[kk] ;
00298    }
00299 
00300    /*---------------------- make a new dataset ----------------------*/
00301 
00302    new_dset = EDIT_empty_copy( old_dset ) ;
00303 
00304    tross_Copy_History( old_dset , new_dset ) ;
00305    tross_Make_History( "3dTsmooth" , argc,argv , new_dset ) ;
00306 
00307    /*-- edit some of its internal parameters --*/
00308 
00309    EDIT_dset_items( new_dset ,
00310                       ADN_prefix      , prefix ,
00311                       ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00312                       ADN_datum_all   , new_datum ,
00313                     ADN_none ) ;
00314 
00315    /*-- make brick(s) for this dataset --*/
00316 
00317    switch( new_datum ){
00318       case MRI_byte:
00319          new_bptr = (byte **) malloc( sizeof(byte *) * ntime ) ;
00320       break ;
00321 
00322       case MRI_short:
00323          new_sptr = (short **) malloc( sizeof(short *) * ntime ) ;
00324       break ;
00325 
00326       case MRI_float:
00327          new_fptr = (float **) malloc( sizeof(float *) * ntime ) ;
00328       break ;
00329    }
00330 
00331    for( kk=0 ; kk < ntime ; kk++ ){
00332       ityp      = DSET_BRICK_TYPE(new_dset,kk) ;   /* type of data */
00333       nbytes    = DSET_BRICK_BYTES(new_dset,kk) ;  /* how much data */
00334       new_brick = malloc( nbytes ) ;               /* make room */
00335 
00336       if( new_brick == NULL ){
00337         fprintf(stderr,"*** Can't get memory for output dataset!\n") ; exit(1) ;
00338       }
00339 
00340       EDIT_substitute_brick( new_dset , kk , ityp , new_brick ) ;
00341 
00342       switch( new_datum ){
00343          case MRI_byte:  new_bptr[kk] = (byte * ) new_brick ; break ;
00344          case MRI_short: new_sptr[kk] = (short *) new_brick ; break ;
00345          case MRI_float: new_fptr[kk] = (float *) new_brick ; break ;
00346       }
00347    }
00348 
00349    if( lwin != NULL && ntap > 0 ){        /* 03 Mar 2001 */
00350       ftap = lwin(ntap) ;
00351       if( lfil == NULL ) lfil = linear_filter_extend ;
00352       if( nwin == CUSTOM ) ntap = custom_ntaps ;
00353    }
00354 
00355    /*----------------------------------------------------*/
00356    /*----- Setup has ended.  Now do some real work. -----*/
00357 
00358    /***** loop over voxels *****/
00359 
00360    for( ii=0 ; ii < nxyz ; ii++  ){  /* 1 time series at a time */
00361 
00362       /*** load data from input dataset, depending on type ***/
00363 
00364       switch( old_datum ){
00365          case MRI_byte:
00366             for( kk=0 ; kk < ntime ; kk++ ) fxar[kk] = bptr[kk][ii] ;
00367          break ;
00368 
00369          case MRI_short:
00370             for( kk=0 ; kk < ntime ; kk++ ) fxar[kk] = sptr[kk][ii] ;
00371          break ;
00372 
00373          case MRI_float:
00374             for( kk=0 ; kk < ntime ; kk++ ) fxar[kk] = fptr[kk][ii] ;
00375          break ;
00376       } /* end of switch over input type */
00377 
00378       if( use_fac )
00379          for( kk=0 ; kk < ntime ; kk++ ) fxar[kk] *= fac[kk] ;
00380 
00381       /* do smoothing */
00382 
00383       if( ftap != NULL )
00384          lfil( ntap,ftap , ntime,fxar ) ;  /* 01 Mar 2001 */
00385       else
00386          smth( ntime , fxar ) ;            /* 3 point smoother */
00387 
00388       /*** put data into output dataset ***/
00389 
00390       switch( new_datum ){
00391 
00392          case MRI_byte:
00393             if( use_fac )
00394                for( kk=0 ; kk < ntime ; kk++ ) new_bptr[kk][ii] = (byte)(fxar[kk] * faci[kk]) ;
00395             else
00396                for( kk=0 ; kk < ntime ; kk++ ) new_bptr[kk][ii] = (byte) fxar[kk] ;
00397          break ;
00398 
00399          case MRI_short:
00400             if( use_fac )
00401                for( kk=0 ; kk < ntime ; kk++ ) new_sptr[kk][ii] = (short)(fxar[kk] * faci[kk]) ;
00402             else
00403                for( kk=0 ; kk < ntime ; kk++ ) new_sptr[kk][ii] = (short) fxar[kk] ;
00404          break ;
00405 
00406          case MRI_float:
00407             if( use_fac )
00408                for( kk=0 ; kk < ntime ; kk++ ) new_fptr[kk][ii] = (float)(fxar[kk] * faci[kk]) ;
00409             else
00410                for( kk=0 ; kk < ntime ; kk++ ) new_fptr[kk][ii] = (float) fxar[kk] ;
00411          break ;
00412       }
00413    }  /* end of loop over voxels */
00414 
00415    DSET_unload(old_dset) ; free(ftap) ;
00416    DSET_write(new_dset) ;
00417    fprintf(stderr,"++ output dataset: %s\n",DSET_BRIKNAME(new_dset)) ;
00418    exit(0) ;
00419 }

00441 {
00442    int ii ;
00443    float aa,bb,cc ;
00444 
00445    bb = vec[0] ; cc = vec[1] ;
00446    for( ii=1 ; ii < num-1 ; ii++ ){
00447       aa = bb ; bb = cc ; cc = vec[ii+1] ;
00448       vec[ii] = MEDIAN(aa,bb,cc) ;         /* see mrilib.h */
00449    }
00450 
00451    return ;
00452 }

00424 {
00425    int ii ;
00426    float aa,bb,cc ;
00427 
00428    bb = vec[0] ; cc = vec[1] ; vec[0] = OSFILT(bb,bb,cc) ;
00429    for( ii=1 ; ii < num-1 ; ii++ ){
00430       aa = bb ; bb = cc ; cc = vec[ii+1] ;
00431       vec[ii] = OSFILT(aa,bb,cc) ;         /* see mrilib.h */
00432    }
00433    vec[num-1] = OSFILT(bb,cc,cc) ;
00434 
00435    return ;
00436 }