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