00001
00002
00003
00004
00005
00006
00007 #include "thd_maker.h"
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067 #undef FREE_WORKSPACE
00068 #define FREE_WORKSPACE \
00069 do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ; \
00070 FREEUP(cptr) ; FREEUP(fxar) ; FREEUP(fac) ; \
00071 FREEUP(fout) ; FREEUP(dtr) ; \
00072 } while(0)
00073
00074
00075 THD_3dim_dataset * MAKER_4D_to_typed_fim( THD_3dim_dataset * old_dset ,
00076 char * new_prefix , int new_datum ,
00077 int ignore , int detrend ,
00078 generic_func * user_func ,
00079 void * user_data )
00080 {
00081 THD_3dim_dataset * new_dset ;
00082
00083 byte ** bptr = NULL ;
00084 short ** sptr = NULL ;
00085 float ** fptr = NULL ;
00086 complex ** cptr = NULL ;
00087
00088 float * fxar = NULL ;
00089 float * fac = NULL ;
00090 float * fout = NULL ;
00091 float * dtr = NULL ;
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
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
00116
00117 old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ;
00118 nuse = DSET_NUM_TIMES(old_dset) - ignore ;
00119 if( nuse < 2 ) return NULL ;
00120
00121 if( new_datum < 0 ) new_datum = old_datum ;
00122 if( new_datum == MRI_complex ) return NULL ;
00123
00124 DSET_load( old_dset ) ;
00125
00126 kk = THD_count_databricks( old_dset->dblk ) ;
00127 if( kk < DSET_NVALS(old_dset) ){
00128 DSET_unload( old_dset ) ;
00129 return NULL ;
00130 }
00131
00132 switch( old_datum ){
00133
00134 default:
00135 DSET_unload( old_dset ) ;
00136 return NULL ;
00137
00138
00139
00140
00141
00142
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
00152
00153
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
00163
00164
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
00174
00175
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 }
00185
00186
00187
00188 fxar = (float *) malloc( sizeof(float) * nuse ) ;
00189 if( fxar == NULL ){ FREE_WORKSPACE ; return NULL ; }
00190
00191
00192
00193 fac = (float *) malloc( sizeof(float) * nuse ) ;
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
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) ;
00213
00214
00215
00216 new_dset = EDIT_empty_copy( old_dset ) ;
00217
00218
00219
00220 ii = EDIT_dset_items(
00221 new_dset ,
00222 ADN_prefix , new_prefix ,
00223 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00224 ADN_datum_all , new_datum ,
00225 ADN_nvals , 1 ,
00226 ADN_ntt , 0 ,
00227 ADN_type , ISHEAD(old_dset)
00228 ? HEAD_FUNC_TYPE
00229 : GEN_FUNC_TYPE ,
00230 ADN_func_type , FUNC_FIM_TYPE ,
00231 ADN_none ) ;
00232
00233 if( ii != 0 ){
00234 THD_delete_3dim_dataset( new_dset , False ) ;
00235 FREE_WORKSPACE ; return NULL ;
00236 }
00237
00238
00239
00240
00241 nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;
00242
00243 fout = (float *) malloc( sizeof(float) * nvox ) ;
00244
00245 if( fout == NULL ){
00246 THD_delete_3dim_dataset( new_dset , False ) ;
00247 FREE_WORKSPACE ; return NULL ;
00248 }
00249
00250
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
00261
00262
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
00271
00272 for( ii=0 ; ii < nvox ; ii++ ){
00273
00274
00275
00276 switch( old_datum ){
00277
00278
00279
00280 case MRI_byte:
00281 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = bptr[kk][ii] ;
00282 break ;
00283
00284
00285
00286 case MRI_short:
00287 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = sptr[kk][ii] ;
00288 break ;
00289
00290
00291
00292 case MRI_float:
00293 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = fptr[kk][ii] ;
00294 break ;
00295
00296
00297
00298 case MRI_complex:
00299 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = CABS(cptr[kk][ii]) ;
00300 break ;
00301
00302 }
00303
00304
00305
00306 if( use_fac )
00307 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] *= fac[kk] ;
00308
00309
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 ;
00317
00318 ts_mean = x0 ;
00319 ts_slope = x1 / tdelta ;
00320
00321
00322
00323 if( detrend )
00324 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] -= (x0 + x1 * dtr[kk]) ;
00325
00326
00327
00328 iz = ii / nxy ;
00329
00330 if( iz != izold ){
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
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 }
00348
00349 DSET_unload( old_dset ) ;
00350
00351
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 ) ;
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
00367
00368
00369 switch( new_datum ){
00370
00371
00372
00373
00374 case MRI_float:
00375 EDIT_substitute_brick( new_dset , 0 , MRI_float , fout ) ;
00376 fout = NULL ;
00377 break ;
00378
00379
00380
00381
00382 case MRI_short:{
00383 short * bout ;
00384 float sfac ;
00385
00386
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
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
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
00413
00414
00415 case MRI_byte:{
00416 byte * bout ;
00417 float sfac ;
00418
00419
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
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
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 }
00446
00447
00448
00449 FREE_WORKSPACE ;
00450 return new_dset ;
00451 }