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 "thd_maker.h"
8
9 /*-------------------------------------------------------------------------
10 BDW, 24 Feb 1997
11
12 Routine to create a 3D 'fith' dataset from a 3D+time dataset,
13 using a user supplied function to process each time series.
14 "user_func" should like
15
16 void user_func( double tzero , double tdelta ,
17 int npts , float ts[] , double ts_mean , double ts_slope ,
18 void * ud , float * val , float * tval )
19
20 where tzero = time at ts[0]
21 tdelta = time at ts[1] (i.e., ts[k] is at tzero + k*tdelta);
22 tzero and tdelta will be in sec if this is truly "time"
23 npts = number of points in ts array
24 ts = one voxel time series array, ts[0] .. ts[npts-1];
25 note that this will always be a float array, and
26 that ts will start with the "ignore"-th point of
27 the actual voxel time series.
28 ts_mean = mean value of ts array
29 ts_slope = slope of ts array;
30 this will be inversely proportional to tdelta
31 (units of 1/sec);
32 if "detrend" is nonzero, then the mean and slope
33 will been removed from the ts array
34 ud = the user_data pointer passed in here -- this can
35 contain whatever control information the user wants
36 val = pointer to return "intensity" value for this voxel;
37 note that this must be a float
38 tval = pointer to return "threshold" value for this voxel;
39 note that this must be a float between -1 and 1 (inclusive)
40
41 Before the first timeseries is passed, user_func will be called with
42 arguments
43 ( 0.0 , 0.0 , nvox , NULL , 0.0 , 0.0 , user_data , NULL , NULL )
44 where nvox = number of voxels that will be processed.
45 This is to allow for some setup (e.g., malloc, PLUTO_popup_meter, ...).
46
47 After the last timeseries is passed, user_func will be called again with
48 arguments
49 ( 0.0 , 0.0 , 0 , NULL , 0.0 , 0.0 , user_data , NULL , NULL )
50 This is to allow for cleanup (e.g., free of malloc, ...). Note that the
51 only difference between these "notification" calls is the third argument.
52 Note also that the only difference between the user_func for this function
53 and the "to_typed_fim" function is the presence of the last argument.
54
55 The inputs to the present routine are
56 old_dset = pointer to old dataset;
57 note that this dataset must not be warp-on-demand
58 new_prefix = string to use as filename prefix
59 new_datum = type of data to store in output brick;
60 if negative, will use datum from old_dset
61 ignore = number of data points to ignore at the start
62 detrend = if nonzero, this routine will detrend (a+b*t)
63 each time series before passing it to user_func
64 user_func = discussed above
65 user_data = discussed above
66
67 The output consists of pointers to the new dataset. If NULL is returned,
68 some error occurred.
69 ---------------------------------------------------------------------------*/
70 /* BDW, 24 Feb 1997 */
71
72 /*------------------ macro to return workspace at exit -------------------*/
73
74 #undef FREE_WORKSPACE
75 #define FREE_WORKSPACE \
76 do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ; \
77 FREEUP(cptr) ; FREEUP(fxar) ; FREEUP(fac) ; \
78 FREEUP(fout) ; FREEUP(dtr) ; FREEUP(tout) ; \
79 } while(0)
80
81 /*-------------------------------------------------------------------------*/
82 /* BDW, 24 Feb 1997 */
83
MAKER_4D_to_typed_fith(THD_3dim_dataset * old_dset,char * new_prefix,int new_datum,int ignore,int detrend,generic_func * user_func,void * user_data)84 THD_3dim_dataset * MAKER_4D_to_typed_fith( THD_3dim_dataset * old_dset ,
85 char * new_prefix , int new_datum ,
86 int ignore , int detrend ,
87 generic_func * user_func ,
88 void * user_data )
89 {
90 THD_3dim_dataset * new_dset ; /* output dataset */
91
92 byte ** bptr = NULL ; /* one of these will be the array of */
93 short ** sptr = NULL ; /* pointers to input dataset sub-bricks */
94 float ** fptr = NULL ; /* (depending on input datum type) */
95 complex ** cptr = NULL ;
96
97 float * fxar = NULL ; /* array loaded from input dataset */
98 float * fac = NULL ; /* array of brick scaling factors */
99 float * fout = NULL ; /* will be array of output floats (intensity) */
100 float * tout = NULL ; /* will be array of output floats (threshold) */
101 float * dtr = NULL ; /* will be array of detrending coeff */
102
103 float val , d0fac , d1fac , x0,x1;
104 double tzero=0 , tdelta , ts_mean , ts_slope ;
105 int ii , old_datum , nuse , use_fac , iz,izold, nxy,nvox , nbad ;
106 register int kk ;
107
108 void (*ufunc)(double,double,int,float *,double,double,void *,float *,float *)
109 = (void (*)(double,double,int,float *,double,double,void *,float *,float *))user_func ;
110
111 /*----------------------------------------------------------*/
112 /*----- Check inputs to see if they are reasonable-ish -----*/
113
114 if( ! ISVALID_3DIM_DATASET(old_dset) ) return NULL ;
115
116 if( new_datum >= 0 &&
117 new_datum != MRI_short &&
118 new_datum != MRI_float ) return NULL ;
119
120 if( user_func == NULL ) return NULL ;
121
122 if( ignore < 0 ) ignore = 0 ;
123
124 /*--------- set up pointers to each sub-brick in the input dataset ---------*/
125
126 old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* get old dataset datum */
127 nuse = DSET_NUM_TIMES(old_dset) - ignore ; /* # of points on time axis */
128 if( nuse < 2 ) return NULL ;
129
130 if( new_datum < 0 ) new_datum = old_datum ; /* output datum = input */
131 if( new_datum == MRI_complex ) return NULL ; /* but complex = bad news */
132
133 DSET_load( old_dset ) ; /* must be in memory before we get pointers to it */
134
135 kk = THD_count_databricks( old_dset->dblk ) ; /* check if it was */
136 if( kk < DSET_NVALS(old_dset) ){ /* loaded correctly */
137 DSET_unload( old_dset ) ;
138 return NULL ;
139 }
140
141 switch( old_datum ){ /* pointer type depends on input datum type */
142
143 default: /** don't know what to do **/
144 DSET_unload( old_dset ) ;
145 return NULL ;
146
147 /** create array of pointers into old dataset sub-bricks **/
148
149 /*--------- input is bytes ----------*/
150 /* voxel #i at time #k is bptr[k][i] */
151 /* for i=0..nvox-1 and k=0..nuse-1. */
152
153 case MRI_byte:
154 bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;
155 if( bptr == NULL ) return NULL ;
156 for( kk=0 ; kk < nuse ; kk++ )
157 bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk+ignore) ;
158 break ;
159
160 /*--------- input is shorts ---------*/
161 /* voxel #i at time #k is sptr[k][i] */
162 /* for i=0..nvox-1 and k=0..nuse-1. */
163
164 case MRI_short:
165 sptr = (short **) malloc( sizeof(short *) * nuse ) ;
166 if( sptr == NULL ) return NULL ;
167 for( kk=0 ; kk < nuse ; kk++ )
168 sptr[kk] = (short *) DSET_ARRAY(old_dset,kk+ignore) ;
169 break ;
170
171 /*--------- input is floats ---------*/
172 /* voxel #i at time #k is fptr[k][i] */
173 /* for i=0..nvox-1 and k=0..nuse-1. */
174
175 case MRI_float:
176 fptr = (float **) malloc( sizeof(float *) * nuse ) ;
177 if( fptr == NULL ) return NULL ;
178 for( kk=0 ; kk < nuse ; kk++ )
179 fptr[kk] = (float *) DSET_ARRAY(old_dset,kk+ignore) ;
180 break ;
181
182 /*--------- input is complex ---------*/
183 /* voxel #i at time #k is cptr[k][i] */
184 /* for i=0..nvox-1 and k=0..nuse-1. */
185
186 case MRI_complex:
187 cptr = (complex **) malloc( sizeof(complex *) * nuse ) ;
188 if( cptr == NULL ) return NULL ;
189 for( kk=0 ; kk < nuse ; kk++ )
190 cptr[kk] = (complex *) DSET_ARRAY(old_dset,kk+ignore) ;
191 break ;
192
193 } /* end of switch on input type */
194
195 /*---- allocate space for 1 voxel timeseries ----*/
196
197 fxar = (float *) malloc( sizeof(float) * nuse ) ; /* voxel timeseries */
198 if( fxar == NULL ){ FREE_WORKSPACE ; return NULL ; }
199
200 /*--- get scaling factors for sub-bricks ---*/
201
202 fac = (float *) malloc( sizeof(float) * nuse ) ; /* factors */
203 if( fac == NULL ){ FREE_WORKSPACE ; return NULL ; }
204
205 use_fac = 0 ;
206 for( kk=0 ; kk < nuse ; kk++ ){
207 fac[kk] = DSET_BRICK_FACTOR(old_dset,kk+ignore) ;
208 if( fac[kk] != 0.0 ) use_fac++ ;
209 else fac[kk] = 1.0 ;
210 }
211 if( !use_fac ) FREEUP(fac) ;
212
213 /*--- setup for detrending ---*/
214
215 dtr = (float *) malloc( sizeof(float) * nuse ) ;
216 if( dtr == NULL ){ FREE_WORKSPACE ; return NULL ; }
217
218 d0fac = 1.0 / nuse ;
219 d1fac = 12.0 / nuse / (nuse*nuse - 1.0) ;
220 for( kk=0 ; kk < nuse ; kk++ )
221 dtr[kk] = kk - 0.5 * (nuse-1) ; /* linear trend, orthogonal to 1 */
222
223 /*---------------------- make a new dataset ----------------------*/
224
225 new_dset = EDIT_empty_copy( old_dset ) ; /* start with copy of old one */
226
227 /*-- edit some of its internal parameters --*/
228
229 ii = EDIT_dset_items(
230 new_dset ,
231 ADN_prefix , new_prefix , /* filename prefix */
232 ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
233 ADN_datum_all , new_datum , /* atomic datum */
234 ADN_nvals , 2 , /* # sub-bricks */
235 ADN_ntt , 0 , /* # time points */
236 ADN_type , ISHEAD(old_dset) /* dataset type */
237 ? HEAD_FUNC_TYPE
238 : GEN_FUNC_TYPE ,
239 ADN_func_type , FUNC_THR_TYPE , /* function type */
240 ADN_none ) ;
241
242 if( ii != 0 ){
243 ERROR_message("Error creating dataset '%s'",new_prefix) ;
244 THD_delete_3dim_dataset( new_dset , False ) ; /* some error above */
245 FREE_WORKSPACE ; return NULL ;
246 }
247
248 /*------ make floating point output brick
249 (only at the end will scale to byte or shorts) ------*/
250
251 nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;
252
253 fout = (float *) malloc( sizeof(float) * nvox ); /* ptr to intensity brick */
254 tout = (float *) malloc( sizeof(float) * nvox ); /* ptr to threshold brick */
255
256 if( ( fout == NULL ) || (tout == NULL) ){
257 THD_delete_3dim_dataset( new_dset , False ) ;
258 FREE_WORKSPACE ; return NULL ;
259 }
260
261 /*----- set up to find time at each voxel -----*/
262
263 tdelta = old_dset->taxis->ttdel ;
264 if( DSET_TIMEUNITS(old_dset) == UNITS_MSEC_TYPE ) tdelta *= 0.001 ;
265 if( tdelta == 0.0 ) tdelta = 1.0 ;
266
267 izold = -666 ;
268 nxy = old_dset->daxes->nxx * old_dset->daxes->nyy ;
269
270 /*----------------------------------------------------*/
271 /*----- Setup has ended. Now do some real work. -----*/
272
273 /* start notification */
274
275 #if 0
276 user_func( 0.0 , 0.0 , nvox , NULL,0.0,0.0 , user_data , NULL , NULL ) ;
277 #else
278 ufunc( 0.0 , 0.0 , nvox , NULL,0.0,0.0 , user_data , NULL , NULL ) ;
279 #endif
280
281 /***** loop over voxels *****/
282
283 for( ii=0 ; ii < nvox ; ii++ ){ /* 1 time series at a time */
284
285 /*** load data from input dataset, depending on type ***/
286
287 switch( old_datum ){
288
289 /*** input = bytes ***/
290
291 case MRI_byte:
292 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = bptr[kk][ii] ;
293 break ;
294
295 /*** input = shorts ***/
296
297 case MRI_short:
298 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = sptr[kk][ii] ;
299 break ;
300
301 /*** input = floats ***/
302
303 case MRI_float:
304 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = fptr[kk][ii] ;
305 break ;
306
307 /*** input = complex (note we use absolute value) ***/
308
309 case MRI_complex:
310 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = CABS(cptr[kk][ii]) ;
311 break ;
312
313 } /* end of switch over input type */
314
315 /*** scale? ***/
316
317 if( use_fac )
318 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] *= fac[kk] ;
319
320 /** compute mean and slope **/
321
322 x0 = x1 = 0.0 ;
323 for( kk=0 ; kk < nuse ; kk++ ){
324 x0 += fxar[kk] ; x1 += fxar[kk] * dtr[kk] ;
325 }
326
327 x0 *= d0fac ; x1 *= d1fac ; /* factors to remove mean and trend */
328
329 ts_mean = x0 ;
330 ts_slope = x1 / tdelta ;
331
332 /** detrend? **/
333
334 if( detrend )
335 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] -= (x0 + x1 * dtr[kk]) ;
336
337 /** compute start time of this timeseries **/
338
339 iz = ii / nxy ; /* which slice am I in? */
340
341 if( iz != izold ){ /* in a new slice? */
342 tzero = THD_timeof( ignore ,
343 old_dset->daxes->zzorg
344 + iz*old_dset->daxes->zzdel , old_dset->taxis ) ;
345 izold = iz ;
346
347 if( DSET_TIMEUNITS(old_dset) == UNITS_MSEC_TYPE ) tzero *= 0.001 ;
348 }
349
350 /*** compute output ***/
351
352 #if 0
353 user_func( tzero,tdelta , nuse,fxar,ts_mean,ts_slope , user_data ,
354 fout+ii , tout+ii ) ;
355 #else
356 ufunc( tzero,tdelta , nuse,fxar,ts_mean,ts_slope , user_data ,
357 fout+ii , tout+ii ) ;
358 #endif
359
360 /*** limit threshold data to [-1,+1] ***/
361 if (tout[ii] > 1.0) tout[ii] = 1.0;
362 if (tout[ii] < -1.0) tout[ii] = -1.0;
363
364 } /* end of outer loop over 1 voxels at a time */
365
366 DSET_unload( old_dset ) ; /* don't need this no more */
367
368 /* end notification */
369
370 #if 0
371 user_func( 0.0 , 0.0 , 0 , NULL,0.0,0.0 , user_data , NULL , NULL ) ;
372 #else
373 ufunc( 0.0 , 0.0 , 0 , NULL,0.0,0.0 , user_data , NULL , NULL ) ;
374 #endif
375
376 nbad = thd_floatscan(nvox,fout) + thd_floatscan(nvox,tout) ;
377 if( nbad > 0 )
378 fprintf(stderr,
379 "++ Warning: %d bad floats computed in MAKER_4D_to_typed_fith\n\a",
380 nbad ) ;
381
382 /*-------------------------------------------------------------------*/
383 /*------- The output is now in fout[ii] and tout[ii], ii=0..nvox-1.
384 We must now put this into the output dataset --------------*/
385
386 switch( new_datum ){
387
388 /*** output is floats is the simplest:
389 we just have to attach the fout and tout bricks to the dataset ***/
390
391 case MRI_float:
392 EDIT_substitute_brick( new_dset , 0 , MRI_float , fout ) ;
393 EDIT_substitute_brick( new_dset , 1 , MRI_float , tout ) ;
394 fout = NULL ; tout = NULL; /* so it won't be freed later */
395 break ;
396
397 /*** output is shorts:
398 we have to create scaled sub-bricks from fout and tout ***/
399
400 case MRI_short:{
401 short * bout ;
402 short * tbout ;
403 float sfac[2] ;
404
405 /*-- get output sub-bricks --*/
406
407 bout = (short *) malloc( sizeof(short) * nvox ) ;
408 tbout = (short *) malloc( sizeof(short) * nvox ) ;
409 if( (bout == NULL) || (tbout == NULL) ){
410 fprintf(stderr,
411 "\nFinal malloc error in MAKER_4D_to_fith - is memory exhausted?\n\a") ;
412 EXIT(1) ;
413 }
414
415 /*-- find scaling and then scale --*/
416
417 sfac[0] = MCW_vol_amax( nvox,1,1 , MRI_float , fout ) ;
418 if( sfac[0] > 0.0 ){
419 sfac[0] = 32767.0 / sfac[0] ;
420 EDIT_coerce_scale_type( nvox,sfac[0] ,
421 MRI_float,fout , MRI_short,bout ) ;
422 sfac[0] = 1.0 / sfac[0] ;
423 }
424 sfac[1] = FUNC_THR_SCALE_SHORT;
425 EDIT_coerce_scale_type( nvox,sfac[1] ,
426 MRI_float,tout , MRI_short,tbout ) ;
427 sfac[1] = 1.0 / sfac[1];
428
429 /*-- put output bricks into dataset, and store scale factor --*/
430
431 EDIT_substitute_brick( new_dset , 0 , MRI_short , bout ) ;
432 EDIT_substitute_brick( new_dset , 1 , MRI_short , tbout ) ;
433 EDIT_dset_items( new_dset , ADN_brick_fac , sfac , ADN_none ) ;
434 }
435 break ;
436
437
438 } /* end of switch on output data type */
439
440 /*-------------- Cleanup and go home ----------------*/
441
442 FREE_WORKSPACE ;
443 return new_dset ;
444 }
445