Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
thd_tshift.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
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 ;
00026 static int TS_mode = MRI_FOURIER ;
00027
00028
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
00041
00042 if( !ISVALID_DSET(TS_dset) || ignore < 0 ) RETURN(1) ;
00043
00044 nxx = DSET_NX(TS_dset) ;
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 ){
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 ;
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 ){
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
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 ,
00094 ADN_ttdel , TS_TR ,
00095 ADN_tunits , TS_tunits ,
00096 ADN_nsl , 0 ,
00097 ADN_ttorg , 0.0 ,
00098 ADN_ttdur , 0.0 ,
00099 ADN_none ) ;
00100
00101
00102
00103 SHIFT_set_method( TS_mode ) ;
00104
00105 nup = csfft_nextup_one35( ntt+4 ) ;
00106
00107 for( kk=0 ; kk < nzz ; kk++ ){
00108
00109 tshift = (TS_tzero - TS_tpat[kk]) / TS_TR ;
00110 #if 1
00111 tshift = -tshift ;
00112 #endif
00113
00114 if( fabs(tshift) < 0.001 ) continue ;
00115
00116 for( ii=0 ; ii < nxy ; ii+=2 ){
00117
00118 flim = THD_extract_series( ii+kk*nxy , TS_dset , 0 );
00119 far = MRI_FLOAT_PTR(flim) ;
00120
00121 if( TS_rlt == 0 ){
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 ) ;
00129
00130 for( fmin=fmax=far[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
00131 if( far[jj] < fmin ) fmin = far[jj] ;
00132 else if( far[jj] > fmax ) fmax = far[jj] ;
00133 }
00134
00135 if( ii < nxy-1 ){
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 ;
00163 else if( far[jj] > fmax ) far[jj] = fmax ;
00164 switch( TS_rlt ){
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
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
00202
00203 mri_free(flim) ; if( gar != NULL ) mri_free(glim) ;
00204
00205 }
00206
00207 }
00208
00209 if( freepat ){ free(TS_tpat) ; TS_tpat = NULL ; }
00210 TS_tzero = -1.0 ; TS_TR = 0.0 ; TS_tunits = UNITS_SEC_TYPE ;
00211
00212 RETURN(0) ;
00213 }