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
00068
00069
00070
00071
00072
00073
00074
00075
00076 #undef FREE_FOUT
00077 #define FREE_FOUT \
00078 do{ int jv ; \
00079 if( fout != NULL ){ \
00080 for( jv=0 ; jv < nbrik ; jv++ ) FREEUP(fout[jv]) ; \
00081 free(fout) ; \
00082 } } while(0)
00083
00084 #undef FREE_WORKSPACE
00085 #define FREE_WORKSPACE \
00086 do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ; \
00087 FREEUP(cptr) ; FREEUP(fxar) ; FREEUP(fac) ; \
00088 FREEUP(dtr) ; FREEUP(val) ; FREE_FOUT ; \
00089 } while(0)
00090
00091
00092
00093 THD_3dim_dataset * MAKER_4D_to_typed_fbuc( THD_3dim_dataset * old_dset ,
00094 char * new_prefix , int new_datum ,
00095 int ignore , int detrend ,
00096 int nbrik , generic_func * user_func ,
00097 void * user_data )
00098 {
00099 THD_3dim_dataset * new_dset ;
00100
00101 byte ** bptr = NULL ;
00102 short ** sptr = NULL ;
00103 float ** fptr = NULL ;
00104 complex ** cptr = NULL ;
00105
00106 float * fxar = NULL ;
00107 float * fac = NULL ;
00108 float ** fout = NULL ;
00109 float * dtr = NULL ;
00110 float * val = NULL ;
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 ;
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
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
00138
00139 old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ;
00140 nuse = DSET_NUM_TIMES(old_dset) - ignore ;
00141 if( nuse < 2 ) return NULL ;
00142
00143 if( new_datum < 0 ) new_datum = old_datum ;
00144 if( new_datum == MRI_complex ) return NULL ;
00145
00146 DSET_load( old_dset ) ;
00147
00148 kk = THD_count_databricks( old_dset->dblk ) ;
00149 if( kk < DSET_NVALS(old_dset) ){
00150 DSET_unload( old_dset ) ;
00151 return NULL ;
00152 }
00153
00154 switch( old_datum ){
00155
00156 default:
00157 DSET_unload( old_dset ) ;
00158 return NULL ;
00159
00160
00161
00162
00163
00164
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
00174
00175
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
00185
00186
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
00196
00197
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 }
00207
00208
00209
00210 fxar = (float *) malloc( sizeof(float) * nuse ) ;
00211 if( fxar == NULL ){ FREE_WORKSPACE ; return NULL ; }
00212
00213
00214
00215 fac = (float *) malloc( sizeof(float) * nuse ) ;
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
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) ;
00235
00236
00237
00238 new_dset = EDIT_empty_copy( old_dset ) ;
00239
00240
00241
00242 ii = EDIT_dset_items(
00243 new_dset ,
00244 ADN_prefix , new_prefix ,
00245 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00246 ADN_datum_all , new_datum ,
00247 ADN_nvals , nbrik ,
00248 ADN_ntt , 0 ,
00249 ADN_type , ISHEAD(old_dset)
00250 ? HEAD_FUNC_TYPE
00251 : GEN_FUNC_TYPE ,
00252 ADN_func_type , FUNC_BUCK_TYPE ,
00253 ADN_none ) ;
00254
00255 if( ii != 0 ){
00256 THD_delete_3dim_dataset( new_dset , False ) ;
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
00265
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
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
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
00305
00306
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
00315
00316 for( ii=0 ; ii < nvox ; ii++ ){
00317
00318
00319
00320 switch( old_datum ){
00321
00322
00323
00324 case MRI_byte:
00325 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = bptr[kk][ii] ;
00326 break ;
00327
00328
00329
00330 case MRI_short:
00331 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = sptr[kk][ii] ;
00332 break ;
00333
00334
00335
00336 case MRI_float:
00337 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = fptr[kk][ii] ;
00338 break ;
00339
00340
00341
00342 case MRI_complex:
00343 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = CABS(cptr[kk][ii]) ;
00344 break ;
00345
00346 }
00347
00348
00349
00350 if( use_fac )
00351 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] *= fac[kk] ;
00352
00353
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 ;
00361
00362 ts_mean = x0 ;
00363 ts_slope = x1 / tdelta ;
00364
00365
00366
00367 if( detrend )
00368 for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] -= (x0 + x1 * dtr[kk]) ;
00369
00370
00371
00372 iz = ii / nxy ;
00373
00374 if( iz != izold ){
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
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 }
00396
00397 DSET_unload( old_dset ) ;
00398
00399
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
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
00419
00420
00421 switch( new_datum ){
00422
00423
00424
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 ;
00430 }
00431 break ;
00432
00433
00434
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
00462
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 }
00490
00491
00492
00493 FREE_WORKSPACE ;
00494 return new_dset ;
00495 }