1 #include "mrilib.h"
2
3 /*---------------------------------------------------------------------------*/
4 /** This program is really just a quick hack for the Spanish Inquisition!! **/
5 /** [They need to torture some data, according to Chief Inquisitor Javier] **/
6 /*---------------------------------------------------------------------------*/
7
8 /*---------------------------------------------------------------------------*/
9 /* Inputs:
10 ntin = number of time points in input array tsar
11 tsar = input array to be interpolated
12 tgar = time grid of input array -- monotonic increasing, length = ntin
13 toff = additional time offset for input array (allows for slice timing)
14 ntout = number of time points in output array
15 dtout = time step of output array
16 tzout = time offset of output array -- i-th point is at i*dtout+tzout
17 outar = output array -- length = ntout (duh) -- linearly interpolated
18 Function does not check for bad inputs -- that's supposed to be in the
19 domain of the caller function.
20 *//*-------------------------------------------------------------------------*/
21
THD_resample_timeseries_linear(int ntin,float * tsar,float * tgar,float toff,int ntout,float dtout,float tzout,float * outar)22 void THD_resample_timeseries_linear(
23 int ntin , float *tsar, float *tgar, float toff,
24 int ntout, float dtout, float tzout, float *outar )
25 {
26 int iin,iout , ntin2=ntin-2 ;
27 float tout , f0 ;
28
29 /* loop over output points */
30
31 for( iin=0,iout=0 ; iout < ntout ; iout++ ){
32 tout = tzout + iout*dtout ; /* time of this output point */
33
34 /* find next iin such that tout lies in the input interval iin .. iin+1 */
35
36 for( ; iin < ntin2 && tout > tgar[iin+1]+toff ; iin++ ) ; /*nada */
37
38 /* f0 = linear interpolation factor for tsar[iin] */
39
40 f0 = (tgar[iin+1]+toff-tout) / (tgar[iin+1]-tgar[iin]) ;
41
42 /* if beyond time bounds of input, patch the factor */
43
44 if( f0 < 0.0f ) f0 = 0.0f ; else if( f0 > 1.0f ) f0 = 1.0f ;
45
46 /* and interpolate */
47
48 outar[iout] = f0 * tsar[iin] + (1.0f-f0) * tsar[iin+1] ;
49 }
50
51 return ;
52 }
53
54 /*---------------------------------------------------------------------------*/
55 /* Resample an irregular TR dataset to a fixed TR dataset.
56 inset = input dataset
57 tgin = holds input time grid (monotonic increasing)
58 ntout = number of time points in output
59 dtout = TR of output dataset
60 tzout = time offset of output dataset
61 icode = interpolation method (at present, only MRI_LINEAR works)
62 If dtout <= 0, then it will be computed as the average TR
63 If ntout <= 1, then it will be computed to fit the input dataset
64 *//*-------------------------------------------------------------------------*/
65
THD_resample_irregular_dataset(THD_3dim_dataset * inset,MRI_IMAGE * tgin,char * prefix,int ntout,float dtout,float tzout,int icode)66 THD_3dim_dataset * THD_resample_irregular_dataset(
67 THD_3dim_dataset *inset , MRI_IMAGE *tgin ,
68 char *prefix ,
69 int ntout , float dtout , float tzout , int icode )
70 {
71 float *tgar=NULL , *otar,*itar ;
72 float toff ;
73 int ii , ntin , nvox , nbad ;
74 THD_3dim_dataset *outset=NULL ;
75
76 ENTRY("THD_resample_irregular_dataset") ;
77
78 /* check inputs for sanity */
79
80 if( ! ISVALID_DSET(inset) ) RETURN(NULL) ; /* crazy */
81 ntin = DSET_NVALS(inset) ;
82 if( ntin < 2 ) RETURN(NULL) ; /* whacko */
83 DSET_load(inset) ;
84 if( ! DSET_LOADED(inset) ) RETURN(NULL) ; /* just plain nuts */
85
86 if( tgin == NULL || tgin->nx < ntin ) RETURN(NULL) ; /* gonzo */
87
88 /* check input time grid to see if it is monotonic increasing */
89
90 tgar = MRI_FLOAT_PTR(tgin) ; if( tgar == NULL ) RETURN(NULL) ; /* lunacy */
91 for( nbad=0,ii=1 ; ii < ntin ; ii++ ){
92 if( tgar[ii-1] >= tgar[ii] ) nbad++ ;
93 }
94 if( nbad > 0 ){ /* demented */
95 ERROR_message(
96 "THD_resample_irregular_dataset: input time grid is disordered in %d place%s",
97 nbad , (nbad==1) ? "\0" : "s" ) ;
98 RETURN(NULL) ;
99 }
100
101 /* don't have ntout and/or don't have dtout ==> compute them */
102
103 if( ntout < 2 && dtout <= 0.0f ){
104 dtout = ( tgar[ntin-1] - tgar[0] ) / (ntin-1) ; /* average TR */
105 ntout = ntin ;
106 } else if( ntout < 2 ){ /* just missing ntin? */
107 ntout = 1 + (int)( ( tgar[ntin-1] - tzout ) / dtout + 0.001f ) ;
108 } else if( dtout <= 0.0f ){
109 dtout = ( tgar[ntin-1] - tgar[0] ) / (ntout-1) ; /* average TR */
110 }
111 ININFO_message("Output grid: TR=%.3f nvals=%d",dtout,ntout) ;
112
113 /* create output dataset */
114
115 if( !THD_filename_ok(prefix) ) prefix = "./TRfixed" ; /* pathetic user? */
116
117 outset = EDIT_empty_copy(inset) ;
118 EDIT_dset_items( outset ,
119 ADN_prefix , prefix ,
120 ADN_nvals , ntout ,
121 ADN_ntt , ntout ,
122 ADN_ttdel , dtout , /* set TR */
123 ADN_nsl , 0 , /* time shifting is done */
124 ADN_datum_all , MRI_float ,
125 ADN_brick_fac , NULL ,
126 ADN_none ) ;
127 for( ii=0 ; ii < ntout ; ii++ )
128 EDIT_substitute_brick( outset , ii , MRI_float , NULL ) ;
129
130 /* loop over input voxels */
131
132 nvox = DSET_NVOX(inset) ;
133 otar = (float *)malloc(sizeof(float)*ntout) ; /* output array */
134 itar = (float *)malloc(sizeof(float)*ntin ) ; /* input array */
135
136 ININFO_message("Computing output dataset") ;
137
138 for( ii=0 ; ii < nvox ; ii++ ){
139 /* allow for time shifting */
140 toff = THD_timeof_vox( 0 , ii , inset ) ;
141 /* get input array */
142 (void)THD_extract_array( ii , inset , 0 , itar ) ;
143 /* compute output array */
144 THD_resample_timeseries_linear( ntin , itar , tgar , toff ,
145 ntout, dtout, tzout, otar ) ;
146 /* put it into the output dataset */
147 THD_insert_series( ii , outset , ntout , MRI_float , otar , 1 ) ;
148 }
149
150 free(itar); free(otar); /* take out the trash */
151 RETURN(outset);
152 }
153
154 /*---------------------------------------------------------------------------*/
155
TRfix_help(void)156 void TRfix_help(void)
157 {
158 printf(
159 "Usage: 3dTRfix [options]\n"
160 "\n"
161 "This program will read in a dataset that was sampled on an irregular time\n"
162 "grid and re-sample it via linear interpolation to a regular time grid.\n"
163 "\n"
164 "NOTES:\n"
165 "------\n"
166 "The re-sampling will include the effects of slice time offsets (similarly\n"
167 "to program 3dTshift), if these time offsets are encoded in the input dataset's\n"
168 "header.\n"
169 "\n"
170 "No other processing is performed -- in particular, there is no allowance\n"
171 "(at present) for T1 artifacts resulting from variable TR.\n"
172 "\n"
173 "If the first 1 or 2 time points are abnormally bright due to the NMR\n"
174 "pre-steady-state effect, then their influence might be spread farther\n"
175 "into the output dataset by the interpolation process. You can avoid this\n"
176 "effect by excising these values from the input using the '[2..$]' notation\n"
177 "in the input dataset syntax.\n"
178 "\n"
179 "If the input dataset is catenated from multiple non-contiguous imaging runs,\n"
180 "the program will happily interpolate across the time breaks between the runs.\n"
181 "For this reason, you should not give such a file (e.g., from 3dTcat) to this\n"
182 "program -- you should use 3dTRfix on each run separately, and only later\n"
183 "catenate the runs.\n"
184 "\n"
185 "The output dataset is stored in float format, regardless of the input format.\n"
186 "\n"
187 "** Basically, this program is a hack for the Mad Spaniard.\n"
188 "** When are we going out for tapas y cerveza (sangria es bueno, tambien)?\n"
189 "\n"
190 "OPTIONS:\n"
191 "--------\n"
192 "\n"
193 " -input iii = Input dataset 'iii'. [MANDATORY]\n"
194 "\n"
195 " -TRlist rrr = 1D columnar file of time gaps between sub-bricks in 'iii';\n"
196 " If the input dataset has N time points, this file must\n"
197 " have at least N-1 (positive) values.\n"
198 " * Please note that these time steps (or the time values in\n"
199 " '-TIMElist') should be in seconds, NOT in milliseconds!\n"
200 " * AFNI time units are seconds!!!\n"
201 "\n"
202 " -TIMElist ttt = Alternative to '-TRlist', where you give the N values of\n"
203 " the times at each sub-brick; these values must be monotonic\n"
204 " increasing and non-negative.\n"
205 " * You must give exactly one of '-TIMElist' or '-TRlist'.\n"
206 " * The TR value given in the input dataset header is ignored.\n"
207 "\n"
208 " -prefix ppp = Prefix name for the output dataset.\n"
209 "\n"
210 " -TRout ddd = 'ddd' gives the value for the output dataset's TR (in sec).\n"
211 " If '-TRout' is not given, then the average TR of the input\n"
212 " dataset will be used.\n"
213 "\n"
214 "November 2014 -- Zhark the Fixer\n"
215 "\n"
216 ) ;
217 exit(0) ;
218 }
219
220 /*---------------------------------------------------------------------------*/
221
main(int argc,char * argv[])222 int main( int argc , char *argv[] )
223 {
224 int iarg=1 , ntin=0 , ii , nbad=0 ;
225 THD_3dim_dataset *inset=NULL , *outset=NULL ;
226 char *prefix="./TRfix" ;
227 MRI_IMAGE *TRlist=NULL , *TIMElist=NULL ;
228 float dtout=0.0f , *tgar, *dtar ;
229
230 /*-- help the users if we can, they're feeling down down down --*/
231
232 if( argc < 2 || strcasecmp(argv[1],"-help") == 0 ) TRfix_help() ;
233
234 /*-- the legal formalities --*/
235
236 mainENTRY("3dTRfix"); machdep();
237 AFNI_logger("3dTRfix",argc,argv);
238 PRINT_VERSION("3dTRfix"); AUTHOR("Zhark the Fixer");
239
240 /*-- process args --*/
241
242 while( iarg < argc ){
243
244 /*-- read input dataset --*/
245
246 if( strcasecmp(argv[iarg],"-input") == 0 ){
247 if( inset != NULL ) ERROR_exit("Can't use '-input' twice") ;
248 if( ++iarg >= argc ) ERROR_exit("Need argument after '%s'",argv[iarg-1]) ;
249 inset = THD_open_dataset(argv[iarg]) ;
250 if( inset == NULL ) ERROR_exit("Can't open dataset '%s' :-(",argv[iarg]);
251 ntin = DSET_NVALS(inset) ;
252 if( ntin < 2 ) ERROR_exit("Dataset '%s' doesn't have any time dimension",argv[iarg]) ;
253 DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
254 iarg++ ; continue ;
255 }
256
257 /*-- read input 1D file --*/
258
259 if( strcasecmp(argv[iarg],"-TRlist") == 0 ){
260 if( TRlist != NULL ) ERROR_exit("Can only use '-TRlist' once") ;
261 if( TIMElist != NULL ) ERROR_exit("Can't use '-TRlist' after '-TIMElist'") ;
262 if( ++iarg >= argc ) ERROR_exit("Need argument after '%s'",argv[iarg-1]) ;
263 TRlist = mri_read_1D( argv[iarg] ) ;
264 if( TRlist == NULL ) ERROR_exit("Can't open '-TRlist' file '%s'",argv[iarg]) ;
265 if( TRlist->nx == 1 && TRlist->ny > 1 ){
266 MRI_IMAGE *qim = mri_transpose(TRlist); mri_free(TRlist); TRlist = qim ;
267 }
268 iarg++ ; continue ;
269 }
270
271 /*-- read input 1D file --*/
272
273 if( strcasecmp(argv[iarg],"-TIMElist") == 0 ){
274 if( TIMElist != NULL ) ERROR_exit("Can only use '-TIMElist' once") ;
275 if( TRlist != NULL ) ERROR_exit("Can't use '-TIMElist' after '-TRlist'") ;
276 if( ++iarg >= argc ) ERROR_exit("Need argument after '%s'",argv[iarg-1]) ;
277 TIMElist = mri_read_1D( argv[iarg] ) ;
278 if( TIMElist == NULL ) ERROR_exit("Can't open '-TIMElist' file '%s'",argv[iarg]) ;
279 if( TIMElist->nx == 1 && TIMElist->ny > 1 ){
280 MRI_IMAGE *qim = mri_transpose(TIMElist); mri_free(TIMElist); TIMElist = qim ;
281 }
282 iarg++ ; continue ;
283 }
284
285 /*-- read prefix --*/
286
287 if( strcasecmp(argv[iarg],"-prefix") == 0 ){
288 if( ++iarg >= argc ) ERROR_exit("no argument after '%s' :-(",argv[iarg-1]) ;
289 if( !THD_filename_ok(argv[iarg]) )
290 ERROR_exit("badly formed filename: '%s' '%s' :-(",argv[iarg-1],argv[iarg]) ;
291 if( strcmp(argv[iarg],"NULL") == 0 )
292 WARNING_message("'NULL' prefix will result in an unusual output filename") ;
293 prefix = strdup(argv[iarg]) ; iarg++ ; continue ;
294 }
295
296 /*-- read dtout --*/
297
298 if( strcasecmp(argv[iarg],"-TRout") == 0 || strcasecmp(argv[iarg],"-dt") == 0 ){
299 if( ++iarg >= argc ) ERROR_exit("no argument after '%s' :-(",argv[iarg-1]) ;
300 dtout = (float)strtod(argv[iarg],NULL) ;
301 if( dtout <= 0.0f )
302 WARNING_message("-TRout '%s' means 3dTRfix will choose the output TR",argv[iarg]) ;
303 iarg++ ; continue ;
304 }
305
306 /*-- really, what is going on here? --*/
307
308 ERROR_message("Don't recognize this as an option: '%s'",argv[iarg]) ;
309 suggest_best_prog_option(argv[0], argv[iarg]) ;
310 exit(1) ;
311 }
312
313 /*-- check for user sanity --*/
314
315 if( inset == NULL ){
316 ERROR_message("No input dataset???") ; nbad++ ;
317 }
318 if( TRlist == NULL && TIMElist == NULL ){
319 ERROR_message("Neither -TRlist nor -TIMElist was given???") ; nbad++ ;
320 }
321 if( TRlist != NULL && TRlist->nx < ntin-1 ){
322 ERROR_message("-TRlist file has %d values, but input dataset has %d",TRlist->nx,ntin) ;
323 nbad++ ;
324 } else if( TIMElist != NULL && TIMElist->nx < ntin ){
325 ERROR_message("-TIMElist file has %d values, but input dataset has %d",TIMElist->nx,ntin) ;
326 nbad++ ;
327 }
328
329 /* process TRlist into TIMElist, if needed, by the black art of arithmetic */
330
331 if( TRlist != NULL ){
332 TIMElist = mri_new(ntin,1,MRI_float) ;
333 tgar = MRI_FLOAT_PTR(TIMElist) ;
334 dtar = MRI_FLOAT_PTR(TRlist) ;
335 tgar[0] = DSET_TIMEORIGIN(inset) ; /* usually zero */
336 for( ii=1 ; ii < ntin ; ii++ ){
337 tgar[ii] = tgar[ii-1] + dtar[ii-1] ; /* addition! */
338 if( dtar[ii-1] <= 0.0f ){
339 ERROR_message("-TRlist: #%d=%g -- but must be positive!",ii-1,dtar[ii-1]) ;
340 nbad++ ;
341 }
342 }
343 } else if( TIMElist != NULL ){ /* check TIMElist for monotonicity */
344 tgar = MRI_FLOAT_PTR(TIMElist) ;
345 for( ii=1 ; ii < ntin ; ii++ ){
346 if( tgar[ii] <= tgar[ii-1] ){
347 ERROR_message("-TIMElist: #%d=%g #%d=%g -- these are out of order!" ,
348 ii-1,tgar[ii-1] , ii,tgar[ii] ) ;
349 nbad++ ;
350 }
351 }
352 }
353
354 if( nbad > 0 )
355 ERROR_exit(" !!! )-: 3dTRfix can't continue after such bad-ness :-( !!!") ;
356
357 /*-- compute something that might be useful (or might not be) --*/
358
359 outset = THD_resample_irregular_dataset( inset , TIMElist , prefix ,
360 0 , dtout , 0.0f , MRI_LINEAR ) ;
361
362 if( outset == NULL )
363 ERROR_exit("Can't compute output dataset for some reason!?") ;
364
365 /*-- run and hide, pretending that things are OK --*/
366
367 INFO_message("Writing output dataset %s",DSET_BRIKNAME(outset)) ;
368 DSET_write(outset) ;
369 exit(0) ;
370 }
371