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