Doxygen Source Code Documentation
thd_makefim.c File Reference
#include "thd_maker.h"Go to the source code of this file.
Defines | |
| #define | FREE_WORKSPACE |
Functions | |
| THD_3dim_dataset * | MAKER_4D_to_typed_fim (THD_3dim_dataset *old_dset, char *new_prefix, int new_datum, int ignore, int detrend, generic_func *user_func, void *user_data) |
Define Documentation
|
|
Value: do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ; \ FREEUP(cptr) ; FREEUP(fxar) ; FREEUP(fac) ; \ FREEUP(fout) ; FREEUP(dtr) ; \ } while(0) Definition at line 68 of file thd_makefim.c. Referenced by MAKER_4D_to_typed_fim(). |
Function Documentation
|
||||||||||||||||||||||||||||||||
|
Definition at line 75 of file thd_makefim.c. References ADN_brick_fac, ADN_datum_all, ADN_func_type, ADN_malloc_type, ADN_none, ADN_ntt, ADN_nvals, ADN_prefix, ADN_type, CABS, DATABLOCK_MEM_MALLOC, THD_3dim_dataset::daxes, THD_3dim_dataset::dblk, DSET_ARRAY, DSET_BRICK_FACTOR, DSET_BRICK_TYPE, DSET_load, DSET_NUM_TIMES, DSET_NVALS, DSET_TIMEUNITS, DSET_unload, EDIT_coerce_scale_type(), EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), EXIT, fout, FREE_WORKSPACE, FREEUP, FUNC_FIM_TYPE, GEN_FUNC_TYPE, generic_func, HEAD_FUNC_TYPE, ISHEAD, ISVALID_3DIM_DATASET, malloc, MCW_vol_amax(), THD_dataxes::nxx, THD_dataxes::nyy, THD_dataxes::nzz, THD_3dim_dataset::taxis, THD_count_databricks(), THD_delete_3dim_dataset(), thd_floatscan(), THD_timeof(), THD_timeaxis::ttdel, UNITS_MSEC_TYPE, user_data, x0, THD_dataxes::zzdel, and THD_dataxes::zzorg. Referenced by PLUTO_4D_to_typed_fim().
00080 {
00081 THD_3dim_dataset * new_dset ; /* output dataset */
00082
00083 byte ** bptr = NULL ; /* one of these will be the array of */
00084 short ** sptr = NULL ; /* pointers to input dataset sub-bricks */
00085 float ** fptr = NULL ; /* (depending on input datum type) */
00086 complex ** cptr = NULL ;
00087
00088 float * fxar = NULL ; /* array loaded from input dataset */
00089 float * fac = NULL ; /* array of brick scaling factors */
00090 float * fout = NULL ; /* will be array of output floats */
00091 float * dtr = NULL ; /* will be array of detrending coeff */
00092
00093 float val , d0fac , d1fac , x0,x1;
00094 double tzero , tdelta , ts_mean , ts_slope ;
00095 int ii , old_datum , nuse , use_fac , iz,izold, nxy,nvox , nbad ;
00096 register int kk ;
00097
00098 void (*ufunc)(double,double,int,float *,double,double,void *,float *)
00099 = (void (*)(double,double,int,float *,double,double,void *,float *)) user_func ;
00100
00101 /*----------------------------------------------------------*/
00102 /*----- Check inputs to see if they are reasonable-ish -----*/
00103
00104 if( ! ISVALID_3DIM_DATASET(old_dset) ) return NULL ;
00105
00106 if( new_datum >= 0 &&
00107 new_datum != MRI_byte &&
00108 new_datum != MRI_short &&
00109 new_datum != MRI_float ) return NULL ;
00110
00111 if( user_func == NULL ) return NULL ;
00112
00113 if( ignore < 0 ) ignore = 0 ;
00114
00115 /*--------- set up pointers to each sub-brick in the input dataset ---------*/
00116
00117 old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* get old dataset datum */
00118 nuse = DSET_NUM_TIMES(old_dset) - ignore ; /* # of points on time axis */
00119 if( nuse < 2 ) return NULL ;
00120
00121 if( new_datum < 0 ) new_datum = old_datum ; /* output datum = input */
00122 if( new_datum == MRI_complex ) return NULL ; /* but complex = bad news */
00123
00124 DSET_load( old_dset ) ; /* must be in memory before we get pointers to it */
00125
00126 kk = THD_count_databricks( old_dset->dblk ) ; /* check if it was */
00127 if( kk < DSET_NVALS(old_dset) ){ /* loaded correctly */
00128 DSET_unload( old_dset ) ;
00129 return NULL ;
00130 }
00131
00132 switch( old_datum ){ /* pointer type depends on input datum type */
00133
00134 default: /** don't know what to do **/
00135 DSET_unload( old_dset ) ;
00136 return NULL ;
00137
00138 /** create array of pointers into old dataset sub-bricks **/
00139
00140 /*--------- input is bytes ----------*/
00141 /* voxel #i at time #k is bptr[k][i] */
00142 /* for i=0..nvox-1 and k=0..nuse-1. */
00143
00144 case MRI_byte:
00145 bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;
00146 if( bptr == NULL ) return NULL ;
00147 for( kk=0 ; kk < nuse ; kk++ )
00148 bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk+ignore) ;
00149 break ;
00150
00151 /*--------- input is shorts ---------*/
00152 /* voxel #i at time #k is sptr[k][i] */
00153 /* for i=0..nvox-1 and k=0..nuse-1. */
00154
00155 case MRI_short:
00156 sptr = (short **) malloc( sizeof(short *) * nuse ) ;
00157 if( sptr == NULL ) return NULL ;
00158 for( kk=0 ; kk < nuse ; kk++ )
00159 sptr[kk] = (short *) DSET_ARRAY(old_dset,kk+ignore) ;
00160 break ;
00161
00162 /*--------- input is floats ---------*/
00163 /* voxel #i at time #k is fptr[k][i] */
00164 /* for i=0..nvox-1 and k=0..nuse-1. */
00165
00166 case MRI_float:
00167 fptr = (float **) malloc( sizeof(float *) * nuse ) ;
00168 if( fptr == NULL ) return NULL ;
00169 for( kk=0 ; kk < nuse ; kk++ )
00170 fptr[kk] = (float *) DSET_ARRAY(old_dset,kk+ignore) ;
00171 break ;
00172
00173 /*--------- input is complex ---------*/
00174 /* voxel #i at time #k is cptr[k][i] */
00175 /* for i=0..nvox-1 and k=0..nuse-1. */
00176
00177 case MRI_complex:
00178 cptr = (complex **) malloc( sizeof(complex *) * nuse ) ;
00179 if( cptr == NULL ) return NULL ;
00180 for( kk=0 ; kk < nuse ; kk++ )
00181 cptr[kk] = (complex *) DSET_ARRAY(old_dset,kk+ignore) ;
00182 break ;
00183
00184 } /* end of switch on input type */
00185
00186 /*---- allocate space for 1 voxel timeseries ----*/
00187
00188 fxar = (float *) malloc( sizeof(float) * nuse ) ; /* voxel timeseries */
00189 if( fxar == NULL ){ FREE_WORKSPACE ; return NULL ; }
00190
00191 /*--- get scaling factors for sub-bricks ---*/
00192
00193 fac = (float *) malloc( sizeof(float) * nuse ) ; /* factors */
00194 if( fac == NULL ){ FREE_WORKSPACE ; return NULL ; }
00195
00196 use_fac = 0 ;
00197 for( kk=0 ; kk < nuse ; kk++ ){
00198 fac[kk] = DSET_BRICK_FACTOR(old_dset,kk+ignore) ;
00199 if( fac[kk] != 0.0 ) use_fac++ ;
00200 else fac[kk] = 1.0 ;
00201 }
00202 if( !use_fac ) FREEUP(fac) ;
00203
00204 /*--- setup for detrending ---*/
00205
00206 dtr = (float *) malloc( sizeof(float) * nuse ) ;
00207 if( dtr == NULL ){ FREE_WORKSPACE ; return NULL ; }
00208
00209 d0fac = 1.0 / nuse ;
00210 d1fac = 12.0 / nuse / (nuse*nuse - 1.0) ;
00211 for( kk=0 ; kk < nuse ; kk++ )
00212 dtr[kk] = kk - 0.5 * (nuse-1) ; /* linear trend, orthogonal to 1 */
00213
00214 /*---------------------- make a new dataset ----------------------*/
00215
00216 new_dset = EDIT_empty_copy( old_dset ) ; /* start with copy of old one */
00217
00218 /*-- edit some of its internal parameters --*/
00219
00220 ii = EDIT_dset_items(
00221 new_dset ,
00222 ADN_prefix , new_prefix , /* filename prefix */
00223 ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
00224 ADN_datum_all , new_datum , /* atomic datum */
00225 ADN_nvals , 1 , /* # sub-bricks */
00226 ADN_ntt , 0 , /* # time points */
00227 ADN_type , ISHEAD(old_dset) /* dataset type */
00228 ? HEAD_FUNC_TYPE
00229 : GEN_FUNC_TYPE ,
00230 ADN_func_type , FUNC_FIM_TYPE , /* function type */
00231 ADN_none ) ;
00232
00233 if( ii != 0 ){
00234 THD_delete_3dim_dataset( new_dset , False ) ; /* some error above */
00235 FREE_WORKSPACE ; return NULL ;
00236 }
00237
00238 /*------ make floating point output brick
00239 (only at the end will scale to byte or shorts) ------*/
00240
00241 nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;
00242
00243 fout = (float *) malloc( sizeof(float) * nvox ) ; /* ptr to brick */
00244
00245 if( fout == NULL ){
00246 THD_delete_3dim_dataset( new_dset , False ) ;
00247 FREE_WORKSPACE ; return NULL ;
00248 }
00249
00250 /*----- set up to find time at each voxel -----*/
00251
00252 tdelta = old_dset->taxis->ttdel ;
00253 if( DSET_TIMEUNITS(old_dset) == UNITS_MSEC_TYPE ) tdelta *= 0.001 ;
00254 if( tdelta == 0.0 ) tdelta = 1.0 ;
00255
00256 izold = -666 ;
00257 nxy = old_dset->daxes->nxx * old_dset->daxes->nyy ;
00258
00259 /*----------------------------------------------------*/
00260 /*----- Setup has ended. Now do some real work. -----*/
00261
00262 /* start notification */
00263
00264 #if 0
00265 user_func( 0.0 , 0.0 , nvox , NULL,0.0,0.0 , user_data , NULL ) ;
00266 #else
00267 ufunc( 0.0 , 0.0 , nvox , NULL,0.0,0.0 , user_data , NULL ) ;
00268 #endif
00269
00270 /***** loop over voxels *****/
00271
00272 for( ii=0 ; ii < nvox ; ii++ ){ /* 1 time series at a time */
00273
00274 /*** load data from input dataset, depending on type ***/
00275
00276 switch( old_datum ){
00277
00278 /*** input = bytes ***/
00279
00280 case MRI_byte:
00281 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = bptr[kk][ii] ;
00282 break ;
00283
00284 /*** input = shorts ***/
00285
00286 case MRI_short:
00287 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = sptr[kk][ii] ;
00288 break ;
00289
00290 /*** input = floats ***/
00291
00292 case MRI_float:
00293 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = fptr[kk][ii] ;
00294 break ;
00295
00296 /*** input = complex (note we use absolute value) ***/
00297
00298 case MRI_complex:
00299 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = CABS(cptr[kk][ii]) ;
00300 break ;
00301
00302 } /* end of switch over input type */
00303
00304 /*** scale? ***/
00305
00306 if( use_fac )
00307 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] *= fac[kk] ;
00308
00309 /** compute mean and slope **/
00310
00311 x0 = x1 = 0.0 ;
00312 for( kk=0 ; kk < nuse ; kk++ ){
00313 x0 += fxar[kk] ; x1 += fxar[kk] * dtr[kk] ;
00314 }
00315
00316 x0 *= d0fac ; x1 *= d1fac ; /* factors to remove mean and trend */
00317
00318 ts_mean = x0 ;
00319 ts_slope = x1 / tdelta ;
00320
00321 /** detrend? **/
00322
00323 if( detrend )
00324 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] -= (x0 + x1 * dtr[kk]) ;
00325
00326 /** compute start time of this timeseries **/
00327
00328 iz = ii / nxy ; /* which slice am I in? */
00329
00330 if( iz != izold ){ /* in a new slice? */
00331 tzero = THD_timeof( ignore ,
00332 old_dset->daxes->zzorg
00333 + iz*old_dset->daxes->zzdel , old_dset->taxis ) ;
00334 izold = iz ;
00335
00336 if( DSET_TIMEUNITS(old_dset) == UNITS_MSEC_TYPE ) tzero *= 0.001 ;
00337 }
00338
00339 /*** compute output ***/
00340
00341 #if 0
00342 user_func( tzero,tdelta , nuse,fxar,ts_mean,ts_slope , user_data , fout+ii ) ;
00343 #else
00344 ufunc( tzero,tdelta , nuse,fxar,ts_mean,ts_slope , user_data , fout+ii ) ;
00345 #endif
00346
00347 } /* end of outer loop over 1 voxels at a time */
00348
00349 DSET_unload( old_dset ) ; /* don't need this no more */
00350
00351 /* end notification */
00352
00353 #if 0
00354 user_func( 0.0 , 0.0 , 0 , NULL,0.0,0.0 , user_data , NULL ) ;
00355 #else
00356 ufunc( 0.0 , 0.0 , 0 , NULL,0.0,0.0 , user_data , NULL ) ;
00357 #endif
00358
00359 nbad = thd_floatscan( nvox , fout ) ; /* 08 Aug 2000 */
00360 if( nbad > 0 )
00361 fprintf(stderr,
00362 "++ Warning: %d bad floats computed in MAKER_4D_to_typed_fim\n\a",
00363 nbad ) ;
00364
00365 /*------------------------------------------------------------*/
00366 /*------- The output is now in fout[ii], ii=0..nvox-1.
00367 We must now put this into the output dataset -------*/
00368
00369 switch( new_datum ){
00370
00371 /*** output is floats is the simplest:
00372 we just have to attach the fout brick to the dataset ***/
00373
00374 case MRI_float:
00375 EDIT_substitute_brick( new_dset , 0 , MRI_float , fout ) ;
00376 fout = NULL ; /* so it won't be freed later */
00377 break ;
00378
00379 /*** output is shorts:
00380 we have to create a scaled sub-brick from fout ***/
00381
00382 case MRI_short:{
00383 short * bout ;
00384 float sfac ;
00385
00386 /*-- get output sub-brick --*/
00387
00388 bout = (short *) malloc( sizeof(short) * nvox ) ;
00389 if( bout == NULL ){
00390 fprintf(stderr,
00391 "\nFinal malloc error in MAKER_4D_to_fim - is memory exhausted?\n\a");
00392 EXIT(1) ;
00393 }
00394
00395 /*-- find scaling and then scale --*/
00396
00397 sfac = MCW_vol_amax( nvox,1,1 , MRI_float , fout ) ;
00398 if( sfac > 0.0 ){
00399 sfac = 32767.0 / sfac ;
00400 EDIT_coerce_scale_type( nvox,sfac ,
00401 MRI_float,fout , MRI_short,bout ) ;
00402 sfac = 1.0 / sfac ;
00403 }
00404
00405 /*-- put output brick into dataset, and store scale factor --*/
00406
00407 EDIT_substitute_brick( new_dset , 0 , MRI_short , bout ) ;
00408 EDIT_dset_items( new_dset , ADN_brick_fac , &sfac , ADN_none ) ;
00409 }
00410 break ;
00411
00412 /*** output is bytes (byte = unsigned char)
00413 we have to create a scaled sub-brick from fout ***/
00414
00415 case MRI_byte:{
00416 byte * bout ;
00417 float sfac ;
00418
00419 /*-- get output sub-brick --*/
00420
00421 bout = (byte *) malloc( sizeof(byte) * nvox ) ;
00422 if( bout == NULL ){
00423 fprintf(stderr,
00424 "\nFinal malloc error in MAKER_4D_to_fim - is memory exhausted?\n\a");
00425 EXIT(1) ;
00426 }
00427
00428 /*-- find scaling and then scale --*/
00429
00430 sfac = MCW_vol_amax( nvox,1,1 , MRI_float , fout ) ;
00431 if( sfac > 0.0 ){
00432 sfac = 255.0 / sfac ;
00433 EDIT_coerce_scale_type( nvox,sfac ,
00434 MRI_float,fout , MRI_byte,bout ) ;
00435 sfac = 1.0 / sfac ;
00436 }
00437
00438 /*-- put output brick into dataset, and store scale factor --*/
00439
00440 EDIT_substitute_brick( new_dset , 0 , MRI_byte , bout ) ;
00441 EDIT_dset_items( new_dset , ADN_brick_fac , &sfac , ADN_none ) ;
00442 }
00443 break ;
00444
00445 } /* end of switch on output data type */
00446
00447 /*-------------- Cleanup and go home ----------------*/
00448
00449 FREE_WORKSPACE ;
00450 return new_dset ;
00451 }
|
