Doxygen Source Code Documentation
thd_maker.h File Reference
#include "mrilib.h"Go to the source code of this file.
Defines | |
| #define | FREEUP(x) do{if((x) != NULL){free((x)); (x)=NULL;}}while(0) |
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) |
| 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) |
| 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
|
|
Definition at line 12 of file thd_maker.h. |
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 }
|
|
||||||||||||||||||||||||||||||||
|
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 }
|
|
||||||||||||||||||||||||||||||||
|
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 }
|
