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 "mrilib.h"
00008 
00009 /*-----------------------------------------------------------------------------*/
00010 /*! Time shift a dataset to a common origin.  Stores the data on top of the
00011    input dataset!  Return value is 0 for good processing, nonzero for errors.
00012    This routine is like doing the default actions of 3dTshift (from which
00013    the code is adapted).  -- RWCox -- 15 Feb 2001
00014 -------------------------------------------------------------------------------*/
00015 
00016  /*-- At present, these static variables are not alterable from the
00017       outside; they correspond to the command line options of 3dTshift.
00018       In the future, there may be a way to changes these to affect things. --*/
00019 
00020 static float   TS_TR      = 0.0 ;
00021 static int     TS_tunits  = UNITS_SEC_TYPE ;
00022 static float * TS_tpat    = NULL ;
00023 static float   TS_tzero   = -1.0 ;
00024 static int     TS_slice   = -1 ;
00025 static int     TS_rlt     = 0 ;   /* 0=add both in; 1=add neither; 2=add mean */
00026 static int     TS_mode    = MRI_FOURIER ; /* one of MRI_FOURIER MRI_CUBIC
00027                                                     MRI_QUINTIC MRI_HEPTIC
00028                                                     MRI_LINEAR             */
00029 
00030 int THD_dataset_tshift( THD_3dim_dataset * TS_dset , int ignore )
00031 {
00032    int nzz, ii,jj,kk , ntt,nxx,nyy,nxy , nup , freepat=0 ;
00033    float tomax,tomin , tshift , fmin,fmax , gmin,gmax , f0,f1 , g0,g1 ;
00034    float ffmin,ffmax , ggmin,ggmax ;
00035    MRI_IMAGE * flim , * glim ;
00036    float * far , * gar ;
00037 
00038 ENTRY("THD_dataset_tshift") ;
00039 
00040    /*- extract dataset values, check for errors -*/
00041 
00042    if( !ISVALID_DSET(TS_dset) || ignore < 0 ) RETURN(1) ;
00043 
00044    nxx = DSET_NX(TS_dset) ;                      /* get dimensions */
00045    nyy = DSET_NY(TS_dset) ; nxy = nxx * nyy ;
00046    nzz = DSET_NZ(TS_dset) ;
00047    ntt = DSET_NVALS(TS_dset) ;
00048    if( ignore > ntt-4 ) RETURN(1) ;
00049 
00050    if( DSET_NVALS(TS_dset) < 2 ) RETURN(1) ;
00051    if( TS_slice >= nzz ) RETURN(1) ;
00052 
00053    if( TS_dset->taxis == NULL ){
00054       if( TS_TR == 0.0 || TS_tpat == NULL ) RETURN(1) ;
00055    } else if( TS_tpat == NULL && TS_dset->taxis->toff_sl == NULL ){
00056       RETURN(1) ;
00057    }
00058 
00059    if( TS_TR == 0.0 ){                                    /* set TR from dataset */
00060       TS_TR     = DSET_TIMESTEP(TS_dset) ;
00061       TS_tunits = TS_dset->taxis->units_type ;
00062    }
00063 
00064    if( TS_tpat == NULL ){
00065       if( TS_dset->taxis->nsl < nzz ) RETURN(1) ;
00066       TS_tpat = (float *) malloc( sizeof(float) * nzz ) ;
00067       memcpy( TS_tpat , TS_dset->taxis->toff_sl , sizeof(float)*nzz ) ;
00068       freepat = 1 ;
00069    }
00070 
00071    tomin = WAY_BIG ; tomax = -WAY_BIG ;                      /* check pattern */
00072    for( ii=0 ; ii < nzz ; ii++ ){
00073       if( TS_tpat[ii] > tomax ) tomax = TS_tpat[ii] ;
00074       if( TS_tpat[ii] < tomin ) tomin = TS_tpat[ii] ;
00075    }
00076    if( tomin < 0.0 || tomax > TS_TR ) RETURN(1) ;
00077    else if( tomin >= tomax )          RETURN(1) ;
00078 
00079    if( TS_slice >= 0 && TS_slice < nzz ){                   /* set common time point */
00080       TS_tzero = TS_tpat[TS_slice] ;
00081    } else if( TS_tzero < 0.0 ){
00082       TS_tzero = 0.0 ;
00083       for( ii=0 ; ii < nzz ; ii++ ) TS_tzero += TS_tpat[ii] ;
00084       TS_tzero /= nzz ;
00085    }
00086 
00087    /*- reconfigure the time axis -*/
00088 
00089    DSET_mallocize( TS_dset) ;
00090    DSET_load( TS_dset ) ; if( !DSET_LOADED(TS_dset) ) RETURN(1) ;
00091 
00092    EDIT_dset_items( TS_dset ,
00093                        ADN_ntt    , ntt       ,  /* in case not already set */
00094                        ADN_ttdel  , TS_TR     ,  /* may have changed */
00095                        ADN_tunits , TS_tunits ,  /* may have changed */
00096                        ADN_nsl    , 0         ,  /* will have no offsets when done */
00097                        ADN_ttorg  , 0.0       ,  /* in case not already set */
00098                        ADN_ttdur  , 0.0       ,  /* in case not already set */
00099                     ADN_none ) ;
00100 
00101    /*---- do the temporal shifting! ----*/
00102 
00103    SHIFT_set_method( TS_mode ) ;
00104 
00105    nup = csfft_nextup_one35( ntt+4 ) ;
00106 
00107    for( kk=0 ; kk < nzz ; kk++ ){                            /* loop over slices */
00108 
00109       tshift = (TS_tzero - TS_tpat[kk]) / TS_TR ;  /* rightward fractional shift */
00110 #if 1
00111       tshift = -tshift ;  /* 24 Apr 2003 -- OOG */
00112 #endif
00113 
00114       if( fabs(tshift) < 0.001 ) continue ;                   /* skip this slice */
00115 
00116       for( ii=0 ; ii < nxy ; ii+=2 ){          /* loop over voxel pairs in slice */
00117 
00118          flim = THD_extract_series( ii+kk*nxy , TS_dset , 0 ); /* get this voxel */
00119          far  = MRI_FLOAT_PTR(flim) ;
00120 
00121          if( TS_rlt == 0 ){                             /* range of data: before */
00122             for( ffmin=ffmax=far[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
00123                     if( far[jj] < ffmin ) ffmin = far[jj] ;
00124                else if( far[jj] > ffmax ) ffmax = far[jj] ;
00125             }
00126          }
00127 
00128          THD_linear_detrend( ntt-ignore , far+ignore , &f0,&f1 ) ;             /* remove trend */
00129 
00130          for( fmin=fmax=far[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
00131                  if( far[jj] < fmin ) fmin = far[jj] ;   /* range of data: after */
00132             else if( far[jj] > fmax ) fmax = far[jj] ;
00133          }
00134 
00135          if( ii < nxy-1 ){                                       /* get next voxel */
00136             glim = THD_extract_series( ii+kk*nxy+1 , TS_dset , 0 ) ;
00137             gar  = MRI_FLOAT_PTR(glim) ;
00138             if( TS_rlt == 0 ){
00139                for( ggmin=ggmax=gar[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
00140                        if( gar[jj] < ggmin ) ggmin = gar[jj] ;
00141                   else if( gar[jj] > ggmax ) ggmax = gar[jj] ;
00142                }
00143             }
00144 
00145             THD_linear_detrend( ntt-ignore , gar+ignore , &g0,&g1 ) ;
00146 
00147             for( gmin=gmax=gar[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
00148                     if( gar[jj] < gmin ) gmin = gar[jj] ;
00149                else if( gar[jj] > gmax ) gmax = gar[jj] ;
00150             }
00151 
00152          } else {
00153             gar  = NULL ;
00154          }
00155 
00156          if( gar != NULL )
00157             SHIFT_two_rows( ntt-ignore,nup, tshift,far+ignore , tshift, gar+ignore ) ;
00158          else
00159             SHIFT_two_rows( ntt-ignore,nup, tshift,far+ignore , tshift, NULL ) ;
00160 
00161          for( jj=ignore ; jj < ntt ; jj++ ){
00162                  if( far[jj] < fmin ) far[jj] = fmin ;           /* clip to input range */
00163             else if( far[jj] > fmax ) far[jj] = fmax ;
00164             switch( TS_rlt ){                                    /* restore trend? */
00165                case 0:
00166                   far[jj] += (f0 + (jj-ignore)*f1) ;
00167                        if( far[jj] < ffmin ) far[jj] = ffmin ;
00168                   else if( far[jj] > ffmax ) far[jj] = ffmax ;
00169                break ;
00170 
00171                case 2:
00172                   far[jj] += f0 ;
00173                break ;
00174             }
00175          }
00176 
00177          if( gar != NULL ){
00178             for( jj=ignore ; jj < ntt ; jj++ ){
00179                     if( gar[jj] < gmin ) gar[jj] = gmin ;
00180                else if( gar[jj] > gmax ) gar[jj] = gmax ;
00181                switch( TS_rlt ){
00182                   case 0:
00183                      gar[jj] += (g0 + (jj-ignore)*g1) ;
00184                           if( gar[jj] < ggmin ) gar[jj] = ggmin ;
00185                      else if( gar[jj] > ggmax ) gar[jj] = ggmax ;
00186                   break ;
00187 
00188                   case 2:
00189                      gar[jj] += g0 ;
00190                   break ;
00191                }
00192             }
00193          }
00194 
00195          /* put back into dataset */
00196 
00197          THD_insert_series( ii+kk*nxy , TS_dset , ntt , MRI_float , far , 0 ) ;
00198          if( gar != NULL )
00199             THD_insert_series( ii+kk*nxy+1 , TS_dset , ntt , MRI_float , gar , 0 ) ;
00200 
00201          /* throw out the trash */
00202 
00203          mri_free(flim) ; if( gar != NULL ) mri_free(glim) ;
00204 
00205       } /* end of loop over voxel pairs */
00206 
00207    } /* end of loop over slices */
00208 
00209    if( freepat ){ free(TS_tpat) ; TS_tpat = NULL ; }  /* resets */
00210    TS_tzero = -1.0 ; TS_TR = 0.0 ; TS_tunits = UNITS_SEC_TYPE ;
00211 
00212    RETURN(0) ;
00213 }