1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 /*-----------------------------------------------------------------------------*/
10 /*! Time shift a dataset to a common origin. Stores the data on top of the
11 input dataset! Return value is 0 for good processing, nonzero for errors.
12 This routine is like doing the default actions of 3dTshift (from which
13 the code is adapted). -- RWCox -- 15 Feb 2001
14 -------------------------------------------------------------------------------*/
15
16 /*-- At present, these static variables are not alterable from the
17 outside; they correspond to the command line options of 3dTshift.
18 In the future, there may be a way to changes these to affect things. --*/
19
20 static float TS_TR = 0.0 ;
21 static int TS_tunits = UNITS_SEC_TYPE ;
22 static float * TS_tpat = NULL ;
23 static float TS_tzero = -1.0 ;
24 static int TS_slice = -1 ;
25 static int TS_rlt = 0 ; /* 0=add both in; 1=add neither; 2=add mean */
26 static int TS_mode = MRI_FOURIER ; /* one of MRI_FOURIER MRI_CUBIC
27 MRI_QUINTIC MRI_HEPTIC
28 MRI_LINEAR */
29
THD_dataset_tshift(THD_3dim_dataset * TS_dset,int ignore)30 int THD_dataset_tshift( THD_3dim_dataset * TS_dset , int ignore )
31 {
32 int nzz, ii,jj,kk , ntt,nxx,nyy,nxy , nup , freepat=0 ;
33 float tomax,tomin , tshift , fmin,fmax , gmin=0,gmax=0 , f0,f1 , g0,g1 ;
34 float ffmin=0,ffmax=0 , ggmin=0,ggmax=0 ;
35 MRI_IMAGE *flim , *glim=NULL ;
36 float *far , *gar ;
37
38 ENTRY("THD_dataset_tshift") ;
39
40 /*- extract dataset values, check for errors -*/
41
42 if( !ISVALID_DSET(TS_dset) || ignore < 0 ) RETURN(1) ;
43
44 nxx = DSET_NX(TS_dset) ; /* get dimensions */
45 nyy = DSET_NY(TS_dset) ; nxy = nxx * nyy ;
46 nzz = DSET_NZ(TS_dset) ;
47 ntt = DSET_NVALS(TS_dset) ;
48 if( ignore > ntt-4 ) RETURN(1) ;
49
50 if( DSET_NVALS(TS_dset) < 2 ) RETURN(1) ;
51 if( TS_slice >= nzz ) RETURN(1) ;
52
53 if( TS_dset->taxis == NULL ){
54 if( TS_TR == 0.0 || TS_tpat == NULL ) RETURN(1) ;
55 } else if( TS_tpat == NULL && TS_dset->taxis->toff_sl == NULL ){
56 RETURN(1) ;
57 }
58
59 if( TS_TR == 0.0 ){ /* set TR from dataset */
60 TS_TR = DSET_TIMESTEP(TS_dset) ;
61 TS_tunits = TS_dset->taxis->units_type ;
62 }
63
64 if( TS_tpat == NULL ){
65 if( TS_dset->taxis->nsl < nzz ) RETURN(1) ;
66 TS_tpat = (float *) malloc( sizeof(float) * nzz ) ;
67 memcpy( TS_tpat , TS_dset->taxis->toff_sl , sizeof(float)*nzz ) ;
68 freepat = 1 ;
69 }
70
71 tomin = WAY_BIG ; tomax = -WAY_BIG ; /* check pattern */
72 for( ii=0 ; ii < nzz ; ii++ ){
73 if( TS_tpat[ii] > tomax ) tomax = TS_tpat[ii] ;
74 if( TS_tpat[ii] < tomin ) tomin = TS_tpat[ii] ;
75 }
76 if( tomin < 0.0 || tomax > TS_TR ) RETURN(1) ;
77 else if( tomin >= tomax ) RETURN(1) ;
78
79 if( TS_slice >= 0 && TS_slice < nzz ){ /* set common time point */
80 TS_tzero = TS_tpat[TS_slice] ;
81 } else if( TS_tzero < 0.0 ){
82 TS_tzero = 0.0 ;
83 for( ii=0 ; ii < nzz ; ii++ ) TS_tzero += TS_tpat[ii] ;
84 TS_tzero /= nzz ;
85 }
86
87 /*- reconfigure the time axis -*/
88
89 DSET_mallocize( TS_dset) ;
90 DSET_load( TS_dset ) ; if( !DSET_LOADED(TS_dset) ) RETURN(1) ;
91
92 EDIT_dset_items( TS_dset ,
93 ADN_ntt , ntt , /* in case not already set */
94 ADN_ttdel , TS_TR , /* may have changed */
95 ADN_tunits , TS_tunits , /* may have changed */
96 ADN_nsl , 0 , /* will have no offsets when done */
97 ADN_ttorg , 0.0 , /* in case not already set */
98 ADN_ttdur , 0.0 , /* in case not already set */
99 ADN_none ) ;
100
101 /*---- do the temporal shifting! ----*/
102
103 SHIFT_set_method( TS_mode ) ;
104
105 nup = csfft_nextup_one35( ntt+4 ) ;
106
107 for( kk=0 ; kk < nzz ; kk++ ){ /* loop over slices */
108
109 tshift = (TS_tzero - TS_tpat[kk]) / TS_TR ; /* rightward fractional shift */
110 #if 1
111 tshift = -tshift ; /* 24 Apr 2003 -- OOG */
112 #endif
113
114 if( fabs(tshift) < 0.001 ) continue ; /* skip this slice */
115
116 for( ii=0 ; ii < nxy ; ii+=2 ){ /* loop over voxel pairs in slice */
117
118 flim = THD_extract_series( ii+kk*nxy , TS_dset , 0 ); /* get this voxel */
119 far = MRI_FLOAT_PTR(flim) ;
120
121 if( TS_rlt == 0 ){ /* range of data: before */
122 for( ffmin=ffmax=far[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
123 if( far[jj] < ffmin ) ffmin = far[jj] ;
124 else if( far[jj] > ffmax ) ffmax = far[jj] ;
125 }
126 }
127
128 THD_linear_detrend( ntt-ignore , far+ignore , &f0,&f1 ) ; /* remove trend */
129
130 for( fmin=fmax=far[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
131 if( far[jj] < fmin ) fmin = far[jj] ; /* range of data: after */
132 else if( far[jj] > fmax ) fmax = far[jj] ;
133 }
134
135 if( ii < nxy-1 ){ /* get next voxel */
136 glim = THD_extract_series( ii+kk*nxy+1 , TS_dset , 0 ) ;
137 gar = MRI_FLOAT_PTR(glim) ;
138 if( TS_rlt == 0 ){
139 for( ggmin=ggmax=gar[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
140 if( gar[jj] < ggmin ) ggmin = gar[jj] ;
141 else if( gar[jj] > ggmax ) ggmax = gar[jj] ;
142 }
143 }
144
145 THD_linear_detrend( ntt-ignore , gar+ignore , &g0,&g1 ) ;
146
147 for( gmin=gmax=gar[ignore],jj=ignore+1 ; jj < ntt ; jj++ ){
148 if( gar[jj] < gmin ) gmin = gar[jj] ;
149 else if( gar[jj] > gmax ) gmax = gar[jj] ;
150 }
151
152 } else {
153 gar = NULL ;
154 }
155
156 if( gar != NULL )
157 SHIFT_two_rows( ntt-ignore,nup, tshift,far+ignore , tshift, gar+ignore ) ;
158 else
159 SHIFT_two_rows( ntt-ignore,nup, tshift,far+ignore , tshift, NULL ) ;
160
161 for( jj=ignore ; jj < ntt ; jj++ ){
162 if( far[jj] < fmin ) far[jj] = fmin ; /* clip to input range */
163 else if( far[jj] > fmax ) far[jj] = fmax ;
164 switch( TS_rlt ){ /* restore trend? */
165 case 0:
166 far[jj] += (f0 + (jj-ignore)*f1) ;
167 if( far[jj] < ffmin ) far[jj] = ffmin ;
168 else if( far[jj] > ffmax ) far[jj] = ffmax ;
169 break ;
170
171 case 2:
172 far[jj] += f0 ;
173 break ;
174 }
175 }
176
177 if( gar != NULL ){
178 for( jj=ignore ; jj < ntt ; jj++ ){
179 if( gar[jj] < gmin ) gar[jj] = gmin ;
180 else if( gar[jj] > gmax ) gar[jj] = gmax ;
181 switch( TS_rlt ){
182 case 0:
183 gar[jj] += (g0 + (jj-ignore)*g1) ;
184 if( gar[jj] < ggmin ) gar[jj] = ggmin ;
185 else if( gar[jj] > ggmax ) gar[jj] = ggmax ;
186 break ;
187
188 case 2:
189 gar[jj] += g0 ;
190 break ;
191 }
192 }
193 }
194
195 /* put back into dataset */
196
197 THD_insert_series( ii+kk*nxy , TS_dset , ntt , MRI_float , far , 0 ) ;
198 if( gar != NULL )
199 THD_insert_series( ii+kk*nxy+1 , TS_dset , ntt , MRI_float , gar , 0 ) ;
200
201 /* throw out the trash */
202
203 mri_free(flim) ; if( gar != NULL ) mri_free(glim) ;
204
205 } /* end of loop over voxel pairs */
206
207 } /* end of loop over slices */
208
209 if( freepat ){ free(TS_tpat) ; TS_tpat = NULL ; } /* resets */
210 TS_tzero = -1.0 ; TS_TR = 0.0 ; TS_tunits = UNITS_SEC_TYPE ;
211
212 RETURN(0) ;
213 }
214