00001
00002
00003
00004
00005
00006
00007
00008 #ifndef FIM_THR
00009 #define FIM_THR 0.0999
00010 #endif
00011
00012
00013
00014
00015
00016
00017
00018 THD_3dim_dataset * AFNI_fimmer_compute( Three_D_View * im3d ,
00019 THD_3dim_dataset * dset_time ,
00020 MRI_IMAGE * ref_ts , MRI_IMAGE * ort_ts ,
00021 THD_session * sess , int code, int ucode )
00022 {
00023 THD_3dim_dataset * new_dset=NULL ;
00024 char new_prefix[THD_MAX_PREFIX] ;
00025 char old_prefix[THD_MAX_PREFIX] ;
00026 THD_slist_find fff ;
00027 int ifim , it,iv , nvox , ngood_ref , ntime , it1 , dtyp , nxyz , itbot ;
00028 float * vval, * tsar, * aval, * rbest, * abest, * pbest, * pval, * bbest, * bval;
00029 int * indx ;
00030 short * bar ;
00031 short * ibest ;
00032 void * ptr ;
00033 float stataux[MAX_STAT_AUX] ;
00034 float fthr , topval ;
00035 int nx_ref , ny_ref , ivec , nnow ;
00036 PCOR_references ** pc_ref ;
00037 PCOR_voxel_corr ** pc_vc ;
00038
00039 int fim_nref , nx_ort , ny_ort , internal_ort ;
00040 float * ortar ;
00041 static float * ref_vec = NULL ;
00042 static int nref_vec = -666 ;
00043
00044 int ibr_best , ibr_perc , ibr_fim , ibr_corr , ibr_base , nbrik ;
00045
00046 int polort = im3d->fimdata->polort , ip ;
00047
00048 float top_perc = 0.0 ;
00049
00050 int ibr_pave , ibr_aver ;
00051 float * paval , * avval , * pabest , * avbest ;
00052
00053 int ibr_ptop , ibr_topl , ibr_sigm ;
00054 float * ptval , * tlval , *sgval ,
00055 * ptbest, * tlbest, *sgbest ;
00056
00057 #ifndef DONT_USE_METER
00058 Widget meter = NULL ;
00059 int meter_perc , meter_pold ;
00060 #endif
00061
00062 int nupdt = 0 ,
00063 min_updt = 5 ,
00064 first_updt = 1 ;
00065
00066 ENTRY("AFNI_fimmer_compute") ;
00067
00068
00069
00070 if( ! DSET_GRAPHABLE(dset_time) ||
00071 ref_ts == NULL ||
00072 ref_ts->kind != MRI_float ||
00073 ! IM3D_OPEN(im3d) ||
00074 im3d->type != AFNI_3DDATA_VIEW ||
00075 (code == 0 && ucode == 0) ||
00076 ref_ts->nx < DSET_NUM_TIMES(dset_time) ){
00077
00078 if(PRINT_TRACING)
00079 { char str[256] ;
00080 sprintf(str,"illegal inputs: ntime=%d num_ts=%d",
00081 DSET_NUM_TIMES(dset_time), (ref_ts==NULL) ? (0) : (ref_ts->nx) ) ;
00082 STATUS(str) ; }
00083
00084 RETURN(NULL) ;
00085 }
00086
00087
00088
00089 if( ort_ts != NULL ){
00090 nx_ort = ort_ts->nx ;
00091 ny_ort = ort_ts->ny ;
00092 ortar = MRI_FLOAT_PTR(ort_ts) ;
00093
00094 internal_ort = (nx_ort < DSET_NUM_TIMES(dset_time)) ;
00095 } else {
00096 internal_ort = 1 ;
00097 }
00098 fim_nref = (internal_ort) ? (polort+2) : (ny_ort+polort+2) ;
00099
00100 if( nref_vec < fim_nref ){
00101 ref_vec = (float *) XtRealloc( (char *)ref_vec, sizeof(float)*fim_nref ) ;
00102 nref_vec = fim_nref ;
00103 }
00104
00105 itbot = im3d->fimdata->init_ignore ;
00106 nx_ref = ref_ts->nx ;
00107 ny_ref = ref_ts->ny ;
00108 ntime = DSET_NUM_TIMES(dset_time) ;
00109 ngood_ref = 0 ;
00110 it1 = -1 ;
00111 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00112 tsar = MRI_FLOAT_PTR(ref_ts) + (ivec*nx_ref) ;
00113 ifim = 0 ;
00114 for( it=itbot ; it < ntime ; it++ ){
00115 if( tsar[it] < WAY_BIG ){ ifim++ ; if( it1 < 0 ) it1 = it ; }
00116 }
00117
00118 if( ifim < min_updt ){
00119 STATUS("ref_ts has too few good entries!") ;
00120 RETURN(NULL) ;
00121 }
00122
00123 ngood_ref = MAX( ifim , ngood_ref ) ;
00124 }
00125
00126
00127
00128
00129 dtyp = DSET_BRICK_TYPE(dset_time,it1) ;
00130 if( ! AFNI_GOOD_FUNC_DTYPE(dtyp) ){
00131 STATUS("illegal input data type!") ;
00132 RETURN(NULL) ;
00133 }
00134
00135
00136
00137 MCW_strncpy( old_prefix , DSET_PREFIX(dset_time) , THD_MAX_PREFIX-3 ) ;
00138
00139 if( ! ISVALID_SESSION(sess) ){
00140 sprintf( new_prefix , "%s@%d" , old_prefix , 1 ) ;
00141 } else {
00142 for( ifim=1 ; ifim < 99 ; ifim++ ){
00143 sprintf( new_prefix , "%s@%d" , old_prefix , ifim ) ;
00144 fff = THD_dset_in_session( FIND_PREFIX , new_prefix , sess ) ;
00145 if( fff.dset == NULL ) break ;
00146 }
00147 if( ifim == 99 ){
00148 STATUS("can't create new prefix!") ;
00149 RETURN(NULL) ;
00150 }
00151 }
00152
00153 if(PRINT_TRACING)
00154 { char str[256] ;
00155 sprintf(str,"new prefix = %s",new_prefix) ; STATUS(str) ; }
00156
00157
00158
00159 THD_load_datablock( dset_time->dblk , AFNI_purge_unused_dsets ) ;
00160
00161 nxyz = dset_time->dblk->diskptr->dimsizes[0]
00162 * dset_time->dblk->diskptr->dimsizes[1]
00163 * dset_time->dblk->diskptr->dimsizes[2] ;
00164
00165
00166
00167
00168
00169 switch( dtyp ){
00170
00171 case MRI_short:{
00172 short * dar = (short *) DSET_ARRAY(dset_time,it1) ;
00173 for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += abs(dar[iv]) ;
00174 fthr = FIM_THR * fthr / nxyz ;
00175
00176 if(PRINT_TRACING)
00177 { char str[256] ; sprintf(str,"fthr = %g",fthr) ; STATUS(str) ; }
00178
00179 for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
00180 if( abs(dar[iv]) > fthr ) nvox++ ;
00181 indx = (int *) malloc( sizeof(int) * nvox ) ;
00182 if( indx == NULL ){
00183 fprintf(stderr,"\n*** indx malloc failure in AFNI_fimmer_compute\n") ;
00184 RETURN(NULL) ;
00185 }
00186 for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
00187 if( abs(dar[iv]) > fthr ) indx[nvox++] = iv ;
00188 }
00189 break ;
00190
00191 case MRI_float:{
00192 float * dar = (float *) DSET_ARRAY(dset_time,it1) ;
00193 for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += fabs(dar[iv]) ;
00194 fthr = FIM_THR * fthr / nxyz ;
00195
00196 if(PRINT_TRACING)
00197 { char str[256] ; sprintf(str,"fthr = %g",fthr) ; STATUS(str) ; }
00198
00199 for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
00200 if( fabs(dar[iv]) > fthr ) nvox++ ;
00201 indx = (int *) malloc( sizeof(int) * nvox ) ;
00202 if( indx == NULL ){
00203 fprintf(stderr,"\n*** indx malloc failure in AFNI_fimmer_compute\n") ;
00204 RETURN(NULL) ;
00205 }
00206 for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
00207 if( fabs(dar[iv]) > fthr ) indx[nvox++] = iv ;
00208 }
00209 break ;
00210
00211 case MRI_byte:{
00212 byte * dar = (byte *) DSET_ARRAY(dset_time,it1) ;
00213 for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += dar[iv] ;
00214 fthr = FIM_THR * fthr / nxyz ;
00215
00216 if(PRINT_TRACING)
00217 { char str[256] ; sprintf(str,"fthr = %g",fthr) ; STATUS(str) ; }
00218
00219 for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
00220 if( dar[iv] > fthr ) nvox++ ;
00221 indx = (int *) malloc( sizeof(int) * nvox ) ;
00222 if( indx == NULL ){
00223 fprintf(stderr,"\n*** indx malloc failure in AFNI_fimmer_compute\n") ;
00224 RETURN(NULL) ;
00225 }
00226 for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
00227 if( dar[iv] > fthr ) indx[nvox++] = iv ;
00228 }
00229 break ;
00230 }
00231
00232 if(PRINT_TRACING)
00233 { char str[256] ; sprintf(str,"number of voxels = %d",nvox) ; STATUS(str) ; }
00234
00235
00236
00237 vval = (float *) malloc( sizeof(float) * nvox) ;
00238 if( vval == NULL ){
00239 fprintf(stderr,"\n*** vval malloc failure in AFNI_fimmer_compute\n") ;
00240 free(indx) ; RETURN(NULL) ;
00241 }
00242
00243
00244
00245 ibr_fim=ibr_corr=ibr_best=ibr_perc=ibr_base = -1 ; nbrik = 0 ;
00246 ibr_pave = ibr_aver = -1 ;
00247 ibr_ptop = ibr_topl = ibr_sigm = -1 ;
00248
00249 if( (code & FIM_ALPHA_MASK)!= 0) { ibr_fim = nbrik; nbrik++; }
00250 if( (code & FIM_BEST_MASK) != 0 && ny_ref > 1){ ibr_best = nbrik; nbrik++; }
00251 if( (code & FIM_PERC_MASK) != 0) { ibr_perc = nbrik; nbrik++; }
00252 if( (code & FIM_PAVE_MASK) != 0) { ibr_pave = nbrik; nbrik++; }
00253 if( (code & FIM_BASE_MASK) != 0) { ibr_base = nbrik; nbrik++; }
00254 if( (code & FIM_AVER_MASK) != 0) { ibr_aver = nbrik; nbrik++; }
00255 if( (code & FIM_CORR_MASK) != 0) { ibr_corr = nbrik; nbrik++; }
00256 if( (code & FIM_PTOP_MASK) != 0) { ibr_ptop = nbrik; nbrik++; }
00257 if( (code & FIM_TOPL_MASK) != 0) { ibr_topl = nbrik; nbrik++; }
00258 if( (code & FIM_SIGM_MASK) != 0) { ibr_sigm = nbrik; nbrik++; }
00259
00260
00261
00262 if(PRINT_TRACING)
00263 { char str[256] ; sprintf(str,"number of bricks = %d",nbrik) ; STATUS(str) ; }
00264
00265 if( nbrik == 0 ){
00266
00267 #ifndef DONT_USE_METER
00268 meter = MCW_popup_meter( im3d->vwid->top_shell , METER_TOP_WIDE ) ;
00269 meter_pold = 0 ;
00270 #endif
00271
00272 goto ucode_stuff ;
00273 }
00274
00275
00276
00277 if(PRINT_TRACING)
00278 { char str[256] ;
00279 sprintf(str,"code of FIM_MASKs = %d",code) ; STATUS(str) ; }
00280
00281
00282
00283 if( ny_ref > 1 ){
00284 aval = (float *) malloc(sizeof(float) * nvox) ;
00285 rbest = (float *) malloc(sizeof(float) * nvox) ;
00286 abest = (float *) malloc(sizeof(float) * nvox) ;
00287 ibest = (short *) malloc(sizeof(short) * nvox) ;
00288 pbest = (float *) malloc(sizeof(float) * nvox) ;
00289 bbest = (float *) malloc(sizeof(float) * nvox) ;
00290 pval = (float *) malloc(sizeof(float) * nvox) ;
00291 bval = (float *) malloc(sizeof(float) * nvox) ;
00292
00293 paval = (float *) malloc(sizeof(float) * nvox) ;
00294 avval = (float *) malloc(sizeof(float) * nvox) ;
00295 pabest = (float *) malloc(sizeof(float) * nvox) ;
00296 avbest = (float *) malloc(sizeof(float) * nvox) ;
00297
00298 ptval = (float *) malloc(sizeof(float) * nvox) ;
00299 tlval = (float *) malloc(sizeof(float) * nvox) ;
00300 sgval = (float *) malloc(sizeof(float) * nvox) ;
00301 ptbest = (float *) malloc(sizeof(float) * nvox) ;
00302 tlbest = (float *) malloc(sizeof(float) * nvox) ;
00303 sgbest = (float *) malloc(sizeof(float) * nvox) ;
00304
00305 if( sgbest == NULL ){
00306 fprintf(stderr,"\n*** 'best' malloc failure in AFNI_fimmer_compute\n") ;
00307 free(vval) ; free(indx) ;
00308 if( aval != NULL ) free(aval) ;
00309 if( rbest != NULL ) free(rbest) ;
00310 if( abest != NULL ) free(abest) ;
00311 if( ibest != NULL ) free(ibest) ;
00312 if( pbest != NULL ) free(pbest) ;
00313 if( bbest != NULL ) free(bbest) ;
00314 if( pval != NULL ) free(pval) ;
00315 if( bval != NULL ) free(bval) ;
00316 if( paval != NULL ) free(paval) ;
00317 if( avval != NULL ) free(avval) ;
00318 if( pabest!= NULL ) free(pabest);
00319 if( avbest!= NULL ) free(avbest);
00320 if( ptval != NULL ) free(ptval) ;
00321 if( tlval != NULL ) free(tlval) ;
00322 if( sgval != NULL ) free(sgval) ;
00323 if( ptbest!= NULL ) free(ptbest);
00324 if( tlbest!= NULL ) free(tlbest);
00325 if( sgbest!= NULL ) free(sgbest);
00326 RETURN(NULL) ;
00327 }
00328 } else {
00329 aval = rbest = abest = pbest = bbest = pval = bval = NULL ;
00330 paval = avval = pabest = avbest = NULL ;
00331 ptval = tlval = ptbest = tlbest = NULL ;
00332 sgval = sgbest = NULL ;
00333 ibest = NULL ;
00334 }
00335
00336 if(PRINT_TRACING)
00337 { char str[256] ;
00338 sprintf(str,"nxyz = %d nvox = %d",nxyz,nvox) ; STATUS(str) ; }
00339
00340
00341
00342 pc_ref = (PCOR_references **) malloc( sizeof(PCOR_references *) * ny_ref ) ;
00343 pc_vc = (PCOR_voxel_corr **) malloc( sizeof(PCOR_voxel_corr *) * ny_ref ) ;
00344
00345 if( pc_ref == NULL || pc_vc == NULL ){
00346 free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
00347 if( aval != NULL ) free(aval) ;
00348 if( rbest != NULL ) free(rbest) ;
00349 if( abest != NULL ) free(abest) ;
00350 if( ibest != NULL ) free(ibest) ;
00351 if( pbest != NULL ) free(pbest) ;
00352 if( bbest != NULL ) free(bbest) ;
00353 if( pval != NULL ) free(pval) ;
00354 if( bval != NULL ) free(bval) ;
00355 if( paval != NULL ) free(paval) ;
00356 if( avval != NULL ) free(avval) ;
00357 if( pabest!= NULL ) free(pabest);
00358 if( avbest!= NULL ) free(avbest);
00359 if( ptval != NULL ) free(ptval) ;
00360 if( tlval != NULL ) free(tlval) ;
00361 if( sgval != NULL ) free(sgval) ;
00362 if( ptbest!= NULL ) free(ptbest);
00363 if( tlbest!= NULL ) free(tlbest);
00364 if( sgbest!= NULL ) free(sgbest);
00365 fprintf(stderr,"\n*** FIM initialization fails in AFNI_fimmer_compute\n") ;
00366 RETURN(NULL) ;
00367 }
00368
00369 ifim = 0 ;
00370 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00371 pc_ref[ivec] = new_PCOR_references( fim_nref ) ;
00372 pc_vc[ivec] = new_PCOR_voxel_corr( nvox , fim_nref ) ;
00373 if( pc_ref[ivec] == NULL || pc_vc[ivec] == NULL ) ifim++ ;
00374 }
00375
00376 if( ifim > 0 ){
00377 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00378 free_PCOR_references(pc_ref[ivec]) ;
00379 free_PCOR_voxel_corr(pc_vc[ivec]) ;
00380 }
00381 free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
00382 if( aval != NULL ) free(aval) ;
00383 if( rbest != NULL ) free(rbest) ;
00384 if( abest != NULL ) free(abest) ;
00385 if( ibest != NULL ) free(ibest) ;
00386 if( pbest != NULL ) free(pbest) ;
00387 if( bbest != NULL ) free(bbest) ;
00388 if( pval != NULL ) free(pval) ;
00389 if( bval != NULL ) free(bval) ;
00390 if( paval != NULL ) free(paval) ;
00391 if( avval != NULL ) free(avval) ;
00392 if( pabest!= NULL ) free(pabest);
00393 if( avbest!= NULL ) free(avbest);
00394 if( ptval != NULL ) free(ptval) ;
00395 if( tlval != NULL ) free(tlval) ;
00396 if( sgval != NULL ) free(sgval) ;
00397 if( ptbest!= NULL ) free(ptbest);
00398 if( tlbest!= NULL ) free(tlbest);
00399 if( sgbest!= NULL ) free(sgbest);
00400 fprintf(stderr,"\n*** FIM initialization fails in AFNI_fimmer_compute\n") ;
00401 RETURN(NULL) ;
00402 }
00403
00404
00405
00406 STATUS("making new dataset") ;
00407
00408 new_dset = EDIT_empty_copy( dset_time ) ;
00409
00410 if( nbrik == 1 && ucode == 0 ){
00411 it = EDIT_dset_items( new_dset ,
00412 ADN_prefix , new_prefix ,
00413 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00414 ADN_type , ISHEAD(dset_time)
00415 ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
00416 ADN_func_type , FUNC_FIM_TYPE ,
00417 ADN_nvals , 1 ,
00418 ADN_datum_all , MRI_short ,
00419 ADN_ntt , 0 ,
00420 ADN_none ) ;
00421
00422 } else if( nbrik == 2 && ibr_corr == 1 && ucode == 0 ){
00423 it = EDIT_dset_items( new_dset ,
00424 ADN_prefix , new_prefix ,
00425 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00426 ADN_type , ISHEAD(dset_time)
00427 ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
00428 ADN_func_type , FUNC_COR_TYPE ,
00429 ADN_nvals , 2 ,
00430 ADN_datum_all , MRI_short ,
00431 ADN_ntt , 0 ,
00432 ADN_none ) ;
00433
00434 } else if( nbrik > 0 ){
00435 it = EDIT_dset_items( new_dset ,
00436 ADN_prefix , new_prefix ,
00437 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00438 ADN_type , ISHEAD(dset_time)
00439 ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
00440 ADN_func_type , FUNC_BUCK_TYPE ,
00441 ADN_nvals , nbrik ,
00442 ADN_datum_all , MRI_short ,
00443 ADN_ntt , 0 ,
00444 ADN_none ) ;
00445 } else {
00446 it = 999 ;
00447 }
00448
00449 if( it > 0 ){
00450 fprintf(stderr,
00451 "\n*** EDIT_dset_items error %d in AFNI_fimmer_compute\n",it) ;
00452 THD_delete_3dim_dataset( new_dset , False ) ;
00453 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00454 free_PCOR_references(pc_ref[ivec]) ;
00455 free_PCOR_voxel_corr(pc_vc[ivec]) ;
00456 }
00457 free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
00458 if( aval != NULL ) free(aval) ;
00459 if( rbest != NULL ) free(rbest) ;
00460 if( abest != NULL ) free(abest) ;
00461 if( ibest != NULL ) free(ibest) ;
00462 if( pbest != NULL ) free(pbest) ;
00463 if( bbest != NULL ) free(bbest) ;
00464 if( pval != NULL ) free(pval) ;
00465 if( bval != NULL ) free(bval) ;
00466 if( paval != NULL ) free(paval) ;
00467 if( avval != NULL ) free(avval) ;
00468 if( pabest!= NULL ) free(pabest);
00469 if( avbest!= NULL ) free(avbest);
00470 if( ptval != NULL ) free(ptval) ;
00471 if( tlval != NULL ) free(tlval) ;
00472 if( sgval != NULL ) free(sgval) ;
00473 if( ptbest!= NULL ) free(ptbest);
00474 if( tlbest!= NULL ) free(tlbest);
00475 if( sgbest!= NULL ) free(sgbest);
00476 RETURN(NULL) ;
00477 }
00478
00479
00480
00481 if( ibr_fim >= 0 )
00482 EDIT_BRICK_LABEL( new_dset , ibr_fim , "Fit Coef" ) ;
00483 if( ibr_corr >= 0 )
00484 EDIT_BRICK_LABEL( new_dset , ibr_corr , "Correlation" ) ;
00485 if( ibr_best >= 0 )
00486 EDIT_BRICK_LABEL( new_dset , ibr_best , "Best Index" ) ;
00487 if( ibr_perc >= 0 )
00488 EDIT_BRICK_LABEL( new_dset , ibr_perc , "% Change" ) ;
00489 if( ibr_base >= 0 )
00490 EDIT_BRICK_LABEL( new_dset , ibr_base , "Baseline" ) ;
00491
00492 if( ibr_pave >= 0 )
00493 EDIT_BRICK_LABEL( new_dset , ibr_pave , "% From Ave" ) ;
00494 if( ibr_aver >= 0 )
00495 EDIT_BRICK_LABEL( new_dset , ibr_aver , "Average" ) ;
00496
00497 if( ibr_ptop >= 0 )
00498 EDIT_BRICK_LABEL( new_dset , ibr_ptop , "% From Top" ) ;
00499 if( ibr_topl >= 0 )
00500 EDIT_BRICK_LABEL( new_dset , ibr_topl , "Topline" ) ;
00501 if( ibr_sigm >= 0 )
00502 EDIT_BRICK_LABEL( new_dset , ibr_sigm , "Sigma Resid" ) ;
00503
00504
00505
00506 if( ibr_perc >= 0 || ibr_pave >= 0 || ibr_ptop >= 0 ){
00507 char * cp = my_getenv("AFNI_FIM_PERCENT_LIMIT") ;
00508 if( cp != NULL ){
00509 float tp = strtod(cp,NULL) ;
00510 if( tp > 0.0 ) top_perc = tp ;
00511 }
00512 }
00513
00514
00515
00516 STATUS("making output bricks") ;
00517
00518 for( iv=0 ; iv < new_dset->dblk->nvals ; iv++ ){
00519 ptr = malloc( DSET_BRICK_BYTES(new_dset,iv) ) ;
00520 mri_fix_data_pointer( ptr , DSET_BRICK(new_dset,iv) ) ;
00521 }
00522
00523 if( THD_count_databricks(new_dset->dblk) < new_dset->dblk->nvals ){
00524 fprintf(stderr,
00525 "\n*** failure to malloc new bricks in AFNI_fimmer_compute\n") ;
00526 THD_delete_3dim_dataset( new_dset , False ) ;
00527 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00528 free_PCOR_references(pc_ref[ivec]) ;
00529 free_PCOR_voxel_corr(pc_vc[ivec]) ;
00530 }
00531 free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
00532 if( aval != NULL ) free(aval) ;
00533 if( rbest != NULL ) free(rbest) ;
00534 if( abest != NULL ) free(abest) ;
00535 if( ibest != NULL ) free(ibest) ;
00536 if( pbest != NULL ) free(pbest) ;
00537 if( bbest != NULL ) free(bbest) ;
00538 if( pval != NULL ) free(pval) ;
00539 if( bval != NULL ) free(bval) ;
00540 if( paval != NULL ) free(paval) ;
00541 if( avval != NULL ) free(avval) ;
00542 if( pabest!= NULL ) free(pabest);
00543 if( avbest!= NULL ) free(avbest);
00544 if( ptval != NULL ) free(ptval) ;
00545 if( tlval != NULL ) free(tlval) ;
00546 if( sgval != NULL ) free(sgval) ;
00547 if( ptbest!= NULL ) free(ptbest);
00548 if( tlbest!= NULL ) free(tlbest);
00549 if( sgbest!= NULL ) free(sgbest);
00550 RETURN(NULL) ;
00551 }
00552
00553
00554
00555
00556 #ifndef DONT_USE_METER
00557 meter = MCW_popup_meter( im3d->vwid->top_shell , METER_TOP_WIDE ) ;
00558 meter_pold = 0 ;
00559 #endif
00560
00561 STATUS("starting recursive least squares") ;
00562
00563 for( it=itbot ; it < ntime ; it++ ){
00564
00565 nnow = 0 ;
00566
00567 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00568
00569 tsar = MRI_FLOAT_PTR(ref_ts) + (ivec*nx_ref) ;
00570 if( tsar[it] >= WAY_BIG ) continue ;
00571
00572 ref_vec[0] = 1.0 ;
00573 for( ip=1 ; ip <= polort ; ip++ )
00574 ref_vec[ip] = ref_vec[ip-1] * ((float)it) ;
00575
00576 if( internal_ort ){
00577 ref_vec[ip] = tsar[it] ;
00578 } else {
00579 for( iv=0 ; iv < ny_ort ; iv++ )
00580 ref_vec[iv+ip] = ortar[it + iv*nx_ort] ;
00581
00582 ref_vec[ny_ort+ip] = tsar[it] ;
00583 }
00584
00585 if(PRINT_TRACING)
00586 { char str[256] ;
00587 sprintf(str,"time index=%d ideal[%d]=%f" , it,ivec,tsar[it] ) ;
00588 STATUS(str) ; }
00589
00590
00591
00592 update_PCOR_references( ref_vec , pc_ref[ivec] ) ;
00593
00594
00595
00596 if( nnow == 0 ){
00597 switch( dtyp ){
00598 case MRI_short:{
00599 short * dar = (short *) DSET_ARRAY(dset_time,it) ;
00600 for( iv=0; iv < nvox; iv++ ) vval[iv] = (float) dar[indx[iv]];
00601 }
00602 break ;
00603
00604 case MRI_float:{
00605 float * dar = (float *) DSET_ARRAY(dset_time,it) ;
00606 for( iv=0; iv < nvox; iv++ ) vval[iv] = (float) dar[indx[iv]];
00607 }
00608 break ;
00609
00610 case MRI_byte:{
00611 byte * dar = (byte *) DSET_ARRAY(dset_time,it) ;
00612 for( iv=0; iv < nvox; iv++ ) vval[iv] = (float) dar[indx[iv]];
00613 }
00614 break ;
00615 }
00616 }
00617
00618
00619
00620 PCOR_update_float( vval , pc_ref[ivec] , pc_vc[ivec] ) ;
00621 nnow++ ;
00622 }
00623
00624 if( nnow > 0 ) nupdt++ ;
00625
00626 #ifndef DONT_USE_METER
00627 meter_perc = (int) ( 100.0 * nupdt / ngood_ref ) ;
00628 if( meter_perc != meter_pold ){
00629 MCW_set_meter( meter , meter_perc ) ;
00630 meter_pold = meter_perc ;
00631 }
00632 #endif
00633
00634 }
00635
00636
00637
00638
00639
00640
00641 stataux[0] = nupdt ;
00642 stataux[1] = (ny_ref==1) ? 1 : 2 ;
00643 stataux[2] = fim_nref - 1 ;
00644 for( iv=3 ; iv < MAX_STAT_AUX ; iv++ ) stataux[iv] = 0.0 ;
00645
00646 if( ibr_corr >= 0 ){
00647 if( new_dset->func_type == FUNC_COR_TYPE )
00648 EDIT_dset_items( new_dset, ADN_stat_aux, stataux, ADN_none ) ;
00649
00650 EDIT_BRICK_TO_FICO( new_dset, ibr_corr, stataux[0],stataux[1],stataux[2] ) ;
00651 }
00652
00653 #ifndef DONT_USE_METER
00654 # define METERIZE(ib) do { meter_perc = (int) ( 100.0 * (ib) / nbrik ) ; \
00655 if( meter_perc != meter_pold ){ \
00656 MCW_set_meter( meter , meter_perc ) ; \
00657 meter_pold = meter_perc ; \
00658 } } while(0)
00659 #else
00660 # define METERIZE(ib)
00661 #endif
00662
00663
00664
00665
00666 if( ny_ref == 1 ){
00667
00668
00669
00670 if( ibr_fim >= 0 ){
00671
00672 STATUS("getting 1 ref alpha") ;
00673
00674 PCOR_get_coef( pc_ref[0] , pc_vc[0] , vval ) ;
00675
00676 topval = 0.0 ;
00677 for( iv=0 ; iv < nvox ; iv++ )
00678 if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
00679
00680 bar = DSET_ARRAY( new_dset , ibr_fim ) ;
00681 memset( bar , 0 , sizeof(short)*nxyz ) ;
00682
00683 if( topval > 0.0 ){
00684 topval = MRI_TYPE_maxval[MRI_short] / topval ;
00685 for( iv=0 ; iv < nvox ; iv++ )
00686 bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
00687
00688 stataux[ibr_fim] = 1.0/topval ;
00689 } else {
00690 stataux[ibr_fim] = 0.0 ;
00691 }
00692
00693 METERIZE(ibr_fim) ;
00694 }
00695
00696 if( ibr_corr >= 0 ){
00697
00698 STATUS("getting 1 ref pcor") ;
00699
00700 PCOR_get_pcor( pc_ref[0] , pc_vc[0] , vval ) ;
00701
00702 bar = DSET_ARRAY( new_dset , ibr_corr ) ;
00703 memset( bar , 0 , sizeof(short)*nxyz ) ;
00704
00705 for( iv=0 ; iv < nvox ; iv++ )
00706 bar[indx[iv]] = (short)(FUNC_COR_SCALE_SHORT * vval[iv] + 0.499) ;
00707
00708 stataux[ibr_corr] = 1.0 / FUNC_COR_SCALE_SHORT ;
00709
00710 METERIZE(ibr_corr) ;
00711 }
00712
00713 if( ibr_perc >= 0 ){
00714
00715 STATUS("getting 1 ref perc") ;
00716
00717 PCOR_get_perc( pc_ref[0] , pc_vc[0] , vval , NULL , 0 ) ;
00718
00719 if( top_perc > 0.0 ) EDIT_clip_float( top_perc , nvox , vval ) ;
00720
00721 topval = 0.0 ;
00722 for( iv=0 ; iv < nvox ; iv++ )
00723 if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
00724
00725 bar = DSET_ARRAY( new_dset , ibr_perc ) ;
00726 memset( bar , 0 , sizeof(short)*nxyz ) ;
00727
00728 if( topval > 0.0 ){
00729 topval = MRI_TYPE_maxval[MRI_short] / topval ;
00730 for( iv=0 ; iv < nvox ; iv++ )
00731 bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
00732
00733 stataux[ibr_perc] = 1.0/topval ;
00734 } else {
00735 stataux[ibr_perc] = 0.0 ;
00736 }
00737
00738 METERIZE(ibr_perc) ;
00739 }
00740
00741 if( ibr_pave >= 0 ){
00742
00743 STATUS("getting 1 ref pave") ;
00744
00745 PCOR_get_perc( pc_ref[0] , pc_vc[0] , vval , NULL , 1 ) ;
00746
00747 if( top_perc > 0.0 ) EDIT_clip_float( top_perc , nvox , vval ) ;
00748
00749 topval = 0.0 ;
00750 for( iv=0 ; iv < nvox ; iv++ )
00751 if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
00752
00753 bar = DSET_ARRAY( new_dset , ibr_pave ) ;
00754 memset( bar , 0 , sizeof(short)*nxyz ) ;
00755
00756 if( topval > 0.0 ){
00757 topval = MRI_TYPE_maxval[MRI_short] / topval ;
00758 for( iv=0 ; iv < nvox ; iv++ )
00759 bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
00760
00761 stataux[ibr_pave] = 1.0/topval ;
00762 } else {
00763 stataux[ibr_pave] = 0.0 ;
00764 }
00765
00766 METERIZE(ibr_pave) ;
00767 }
00768
00769 if( ibr_ptop >= 0 ){
00770
00771 STATUS("getting 1 ref ptop") ;
00772
00773 PCOR_get_perc( pc_ref[0] , pc_vc[0] , vval , NULL , 2 ) ;
00774
00775 if( top_perc > 0.0 ) EDIT_clip_float( top_perc , nvox , vval ) ;
00776
00777 topval = 0.0 ;
00778 for( iv=0 ; iv < nvox ; iv++ )
00779 if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
00780
00781 bar = DSET_ARRAY( new_dset , ibr_ptop ) ;
00782 memset( bar , 0 , sizeof(short)*nxyz ) ;
00783
00784 if( topval > 0.0 ){
00785 topval = MRI_TYPE_maxval[MRI_short] / topval ;
00786 for( iv=0 ; iv < nvox ; iv++ )
00787 bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
00788
00789 stataux[ibr_ptop] = 1.0/topval ;
00790 } else {
00791 stataux[ibr_ptop] = 0.0 ;
00792 }
00793
00794 METERIZE(ibr_ptop) ;
00795 }
00796
00797 if( ibr_base >= 0 ){
00798
00799 STATUS("getting 1 ref base") ;
00800
00801 PCOR_get_perc( pc_ref[0] , pc_vc[0] , NULL , vval , 0 ) ;
00802
00803 topval = 0.0 ;
00804 for( iv=0 ; iv < nvox ; iv++ )
00805 if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
00806
00807 bar = DSET_ARRAY( new_dset , ibr_base ) ;
00808 memset( bar , 0 , sizeof(short)*nxyz ) ;
00809
00810 if( topval > 0.0 ){
00811 topval = MRI_TYPE_maxval[MRI_short] / topval ;
00812 for( iv=0 ; iv < nvox ; iv++ )
00813 bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
00814
00815 stataux[ibr_base] = 1.0/topval ;
00816 } else {
00817 stataux[ibr_base] = 0.0 ;
00818 }
00819
00820 METERIZE(ibr_base) ;
00821 }
00822
00823 if( ibr_aver >= 0 ){
00824
00825 STATUS("getting 1 ref aver") ;
00826
00827 PCOR_get_perc( pc_ref[0] , pc_vc[0] , NULL , vval , 1 ) ;
00828
00829 topval = 0.0 ;
00830 for( iv=0 ; iv < nvox ; iv++ )
00831 if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
00832
00833 bar = DSET_ARRAY( new_dset , ibr_aver ) ;
00834 memset( bar , 0 , sizeof(short)*nxyz ) ;
00835
00836 if( topval > 0.0 ){
00837 topval = MRI_TYPE_maxval[MRI_short] / topval ;
00838 for( iv=0 ; iv < nvox ; iv++ )
00839 bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
00840
00841 stataux[ibr_aver] = 1.0/topval ;
00842 } else {
00843 stataux[ibr_aver] = 0.0 ;
00844 }
00845
00846 METERIZE(ibr_aver) ;
00847 }
00848
00849 if( ibr_topl >= 0 ){
00850
00851 STATUS("getting 1 ref topl") ;
00852
00853 PCOR_get_perc( pc_ref[0] , pc_vc[0] , NULL , vval , 2 ) ;
00854
00855 topval = 0.0 ;
00856 for( iv=0 ; iv < nvox ; iv++ )
00857 if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
00858
00859 bar = DSET_ARRAY( new_dset , ibr_topl ) ;
00860 memset( bar , 0 , sizeof(short)*nxyz ) ;
00861
00862 if( topval > 0.0 ){
00863 topval = MRI_TYPE_maxval[MRI_short] / topval ;
00864 for( iv=0 ; iv < nvox ; iv++ )
00865 bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
00866
00867 stataux[ibr_topl] = 1.0/topval ;
00868 } else {
00869 stataux[ibr_topl] = 0.0 ;
00870 }
00871
00872 METERIZE(ibr_topl) ;
00873 }
00874
00875 if( ibr_sigm >= 0 ){
00876
00877 STATUS("getting 1 ref sigm") ;
00878
00879 PCOR_get_stdev( pc_vc[0] , vval ) ;
00880
00881 topval = 0.0 ;
00882 for( iv=0 ; iv < nvox ; iv++ )
00883 if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
00884
00885 bar = DSET_ARRAY( new_dset , ibr_sigm ) ;
00886 memset( bar , 0 , sizeof(short)*nxyz ) ;
00887
00888 if( topval > 0.0 ){
00889 topval = MRI_TYPE_maxval[MRI_short] / topval ;
00890 for( iv=0 ; iv < nvox ; iv++ )
00891 bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
00892
00893 stataux[ibr_sigm] = 1.0/topval ;
00894 } else {
00895 stataux[ibr_sigm] = 0.0 ;
00896 }
00897
00898 METERIZE(ibr_sigm) ;
00899 }
00900
00901 } else {
00902
00903
00904
00905
00906
00907 STATUS("getting first ref results") ;
00908
00909 PCOR_get_coef( pc_ref[0] , pc_vc[0] , abest ) ;
00910 PCOR_get_pcor( pc_ref[0] , pc_vc[0] , rbest ) ;
00911 PCOR_get_perc( pc_ref[0] , pc_vc[0] , pbest , bbest , 0 ) ;
00912 PCOR_get_perc( pc_ref[0] , pc_vc[0] , pabest, avbest, 1 ) ;
00913 PCOR_get_perc( pc_ref[0] , pc_vc[0] , ptbest, tlbest, 2 ) ;
00914 PCOR_get_stdev( pc_vc[0] , sgbest ) ;
00915
00916 for( iv=0 ; iv < nvox ; iv++ ) ibest[iv] = 1 ;
00917
00918
00919
00920
00921
00922 for( ivec=1 ; ivec < ny_ref ; ivec++ ){
00923
00924 STATUS(" == getting results for next ref") ;
00925
00926 PCOR_get_coef( pc_ref[ivec] , pc_vc[ivec] , aval ) ;
00927 PCOR_get_pcor( pc_ref[ivec] , pc_vc[ivec] , vval ) ;
00928 PCOR_get_perc( pc_ref[ivec] , pc_vc[ivec] , pval , bval , 0 ) ;
00929 PCOR_get_perc( pc_ref[ivec] , pc_vc[ivec] , paval, avval, 1 ) ;
00930 PCOR_get_perc( pc_ref[ivec] , pc_vc[ivec] , ptval, tlval, 2 ) ;
00931 PCOR_get_stdev( pc_vc[ivec] , sgval ) ;
00932
00933 STATUS(" == and finding the best results") ;
00934
00935 for( iv=0 ; iv < nvox ; iv++ ){
00936 if( fabs(vval[iv]) > fabs(rbest[iv]) ){
00937 rbest[iv] = vval[iv] ;
00938 abest[iv] = aval[iv] ;
00939 ibest[iv] = (ivec+1) ;
00940 pbest[iv] = pval[iv] ;
00941 bbest[iv] = bval[iv] ;
00942 pabest[iv]= paval[iv] ;
00943 avbest[iv]= avval[iv] ;
00944 ptbest[iv]= ptval[iv] ;
00945 tlbest[iv]= tlval[iv] ;
00946 sgbest[iv]= sgval[iv] ;
00947 }
00948 }
00949 }
00950
00951
00952
00953
00954
00955
00956 if( ibr_fim >= 0 ){
00957
00958 if(PRINT_TRACING)
00959 { char str[256]; sprintf(str,"getting ibr_fim=%d",ibr_fim); STATUS(str); }
00960
00961 topval = 0.0 ;
00962 for( iv=0 ; iv < nvox ; iv++ )
00963 if( fabs(abest[iv]) > topval ) topval = fabs(abest[iv]) ;
00964
00965 bar = DSET_ARRAY( new_dset , ibr_fim ) ;
00966 memset( bar , 0 , sizeof(short)*nxyz ) ;
00967
00968 if( topval > 0.0 ){
00969 topval = MRI_TYPE_maxval[MRI_short] / topval ;
00970 for( iv=0 ; iv < nvox ; iv++ )
00971 bar[indx[iv]] = (short)(topval * abest[iv] + 0.499) ;
00972
00973 stataux[ibr_fim] = 1.0/topval ;
00974 } else {
00975 stataux[ibr_fim] = 0.0 ;
00976 }
00977
00978 METERIZE(ibr_fim) ;
00979 }
00980
00981
00982
00983 if( ibr_corr >= 0 ){
00984
00985 if(PRINT_TRACING)
00986 { char str[256]; sprintf(str,"getting ibr_corr=%d",ibr_corr); STATUS(str); }
00987
00988 bar = DSET_ARRAY( new_dset , ibr_corr ) ;
00989 memset( bar , 0 , sizeof(short)*nxyz ) ;
00990
00991 for( iv=0 ; iv < nvox ; iv++ )
00992 bar[indx[iv]] = (short)(FUNC_COR_SCALE_SHORT * rbest[iv] + 0.499) ;
00993
00994 stataux[ibr_corr] = 1.0 / FUNC_COR_SCALE_SHORT ;
00995
00996 METERIZE(ibr_corr) ;
00997 }
00998
00999
01000
01001 if( ibr_best >= 0 ){
01002
01003 if(PRINT_TRACING)
01004 { char str[256]; sprintf(str,"getting ibr_best=%d",ibr_best); STATUS(str); }
01005
01006 bar = DSET_ARRAY( new_dset , ibr_best ) ;
01007 memset( bar , 0 , sizeof(short)*nxyz ) ;
01008 for( iv=0 ; iv < nvox ; iv++ ) bar[indx[iv]] = ibest[iv] ;
01009 stataux[ibr_best] = 0.0 ;
01010
01011 METERIZE(ibr_best) ;
01012 }
01013
01014
01015
01016 if( ibr_perc >= 0 ){
01017
01018 if(PRINT_TRACING)
01019 { char str[256]; sprintf(str,"getting ibr_perc=%d",ibr_perc); STATUS(str); }
01020
01021 if( top_perc > 0.0 ) EDIT_clip_float( top_perc , nvox , pbest ) ;
01022
01023 topval = 0.0 ;
01024 for( iv=0 ; iv < nvox ; iv++ )
01025 if( fabs(pbest[iv]) > topval ) topval = fabs(pbest[iv]) ;
01026
01027 bar = DSET_ARRAY( new_dset , ibr_perc ) ;
01028 memset( bar , 0 , sizeof(short)*nxyz ) ;
01029
01030 if( topval > 0.0 ){
01031 topval = MRI_TYPE_maxval[MRI_short] / topval ;
01032 for( iv=0 ; iv < nvox ; iv++ )
01033 bar[indx[iv]] = (short)(topval * pbest[iv] + 0.499) ;
01034
01035 stataux[ibr_perc] = 1.0/topval ;
01036 } else {
01037 stataux[ibr_perc] = 0.0 ;
01038 }
01039
01040 METERIZE(ibr_perc) ;
01041 }
01042
01043
01044
01045 if( ibr_pave >= 0 ){
01046
01047 if(PRINT_TRACING)
01048 { char str[256]; sprintf(str,"getting ibr_pave=%d",ibr_pave); STATUS(str); }
01049
01050 if( top_perc > 0.0 ) EDIT_clip_float( top_perc , nvox , pabest ) ;
01051
01052 topval = 0.0 ;
01053 for( iv=0 ; iv < nvox ; iv++ )
01054 if( fabs(pabest[iv]) > topval ) topval = fabs(pabest[iv]) ;
01055
01056 bar = DSET_ARRAY( new_dset , ibr_pave ) ;
01057 memset( bar , 0 , sizeof(short)*nxyz ) ;
01058
01059 if( topval > 0.0 ){
01060 topval = MRI_TYPE_maxval[MRI_short] / topval ;
01061 for( iv=0 ; iv < nvox ; iv++ )
01062 bar[indx[iv]] = (short)(topval * pabest[iv] + 0.499) ;
01063
01064 stataux[ibr_pave] = 1.0/topval ;
01065 } else {
01066 stataux[ibr_pave] = 0.0 ;
01067 }
01068
01069 METERIZE(ibr_pave) ;
01070 }
01071
01072
01073
01074 if( ibr_ptop >= 0 ){
01075
01076 if(PRINT_TRACING)
01077 { char str[256]; sprintf(str,"getting ibr_ptop=%d",ibr_ptop); STATUS(str); }
01078
01079 if( top_perc > 0.0 ) EDIT_clip_float( top_perc , nvox , ptbest ) ;
01080
01081 topval = 0.0 ;
01082 for( iv=0 ; iv < nvox ; iv++ )
01083 if( fabs(ptbest[iv]) > topval ) topval = fabs(ptbest[iv]) ;
01084
01085 bar = DSET_ARRAY( new_dset , ibr_ptop ) ;
01086 memset( bar , 0 , sizeof(short)*nxyz ) ;
01087
01088 if( topval > 0.0 ){
01089 topval = MRI_TYPE_maxval[MRI_short] / topval ;
01090 for( iv=0 ; iv < nvox ; iv++ )
01091 bar[indx[iv]] = (short)(topval * ptbest[iv] + 0.499) ;
01092
01093 stataux[ibr_ptop] = 1.0/topval ;
01094 } else {
01095 stataux[ibr_ptop] = 0.0 ;
01096 }
01097
01098 METERIZE(ibr_ptop) ;
01099 }
01100
01101
01102
01103 if( ibr_base >= 0 ){
01104
01105 if(PRINT_TRACING)
01106 { char str[256]; sprintf(str,"getting ibr_base=%d",ibr_base); STATUS(str); }
01107
01108 topval = 0.0 ;
01109 for( iv=0 ; iv < nvox ; iv++ )
01110 if( fabs(bbest[iv]) > topval ) topval = fabs(bbest[iv]) ;
01111
01112 bar = DSET_ARRAY( new_dset , ibr_base ) ;
01113 memset( bar , 0 , sizeof(short)*nxyz ) ;
01114
01115 if( topval > 0.0 ){
01116 topval = MRI_TYPE_maxval[MRI_short] / topval ;
01117 for( iv=0 ; iv < nvox ; iv++ )
01118 bar[indx[iv]] = (short)(topval * bbest[iv] + 0.499) ;
01119
01120 stataux[ibr_base] = 1.0/topval ;
01121 } else {
01122 stataux[ibr_base] = 0.0 ;
01123 }
01124
01125 METERIZE(ibr_base) ;
01126 }
01127
01128
01129
01130 if( ibr_aver >= 0 ){
01131
01132 if(PRINT_TRACING)
01133 { char str[256]; sprintf(str,"getting ibr_aver=%d",ibr_aver); STATUS(str); }
01134
01135 topval = 0.0 ;
01136 for( iv=0 ; iv < nvox ; iv++ )
01137 if( fabs(avbest[iv]) > topval ) topval = fabs(avbest[iv]) ;
01138
01139 bar = DSET_ARRAY( new_dset , ibr_aver ) ;
01140 memset( bar , 0 , sizeof(short)*nxyz ) ;
01141
01142 if( topval > 0.0 ){
01143 topval = MRI_TYPE_maxval[MRI_short] / topval ;
01144 for( iv=0 ; iv < nvox ; iv++ )
01145 bar[indx[iv]] = (short)(topval * avbest[iv] + 0.499) ;
01146
01147 stataux[ibr_aver] = 1.0/topval ;
01148 } else {
01149 stataux[ibr_aver] = 0.0 ;
01150 }
01151
01152 METERIZE(ibr_aver) ;
01153 }
01154
01155
01156
01157 if( ibr_topl >= 0 ){
01158
01159 if(PRINT_TRACING)
01160 { char str[256]; sprintf(str,"getting ibr_topl=%d",ibr_topl); STATUS(str); }
01161
01162 topval = 0.0 ;
01163 for( iv=0 ; iv < nvox ; iv++ )
01164 if( fabs(tlbest[iv]) > topval ) topval = fabs(tlbest[iv]) ;
01165
01166 bar = DSET_ARRAY( new_dset , ibr_topl ) ;
01167 memset( bar , 0 , sizeof(short)*nxyz ) ;
01168
01169 if( topval > 0.0 ){
01170 topval = MRI_TYPE_maxval[MRI_short] / topval ;
01171 for( iv=0 ; iv < nvox ; iv++ )
01172 bar[indx[iv]] = (short)(topval * tlbest[iv] + 0.499) ;
01173
01174 stataux[ibr_topl] = 1.0/topval ;
01175 } else {
01176 stataux[ibr_topl] = 0.0 ;
01177 }
01178
01179 METERIZE(ibr_topl) ;
01180 }
01181
01182
01183
01184 if( ibr_sigm >= 0 ){
01185
01186 if(PRINT_TRACING)
01187 { char str[256]; sprintf(str,"getting ibr_sigm=%d",ibr_sigm); STATUS(str); }
01188
01189 topval = 0.0 ;
01190 for( iv=0 ; iv < nvox ; iv++ )
01191 if( fabs(sgbest[iv]) > topval ) topval = fabs(sgbest[iv]) ;
01192
01193 bar = DSET_ARRAY( new_dset , ibr_sigm ) ;
01194 memset( bar , 0 , sizeof(short)*nxyz ) ;
01195
01196 if( topval > 0.0 ){
01197 topval = MRI_TYPE_maxval[MRI_short] / topval ;
01198 for( iv=0 ; iv < nvox ; iv++ )
01199 bar[indx[iv]] = (short)(topval * sgbest[iv] + 0.499) ;
01200
01201 stataux[ibr_sigm] = 1.0/topval ;
01202 } else {
01203 stataux[ibr_sigm] = 0.0 ;
01204 }
01205
01206 METERIZE(ibr_sigm) ;
01207 }
01208
01209 }
01210
01211
01212
01213
01214 STATUS("setting brick_fac") ;
01215
01216 (void) EDIT_dset_items( new_dset , ADN_brick_fac , stataux , ADN_none ) ;
01217
01218 #ifndef DONT_USE_METER
01219 MCW_set_meter( meter , 100 ) ;
01220 #endif
01221
01222
01223
01224 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
01225 free_PCOR_references(pc_ref[ivec]) ;
01226 free_PCOR_voxel_corr(pc_vc[ivec]) ;
01227 }
01228 free(pc_ref) ; free(pc_vc) ;
01229 if( aval != NULL ) free(aval) ;
01230 if( rbest != NULL ) free(rbest) ;
01231 if( abest != NULL ) free(abest) ;
01232 if( ibest != NULL ) free(ibest) ;
01233 if( pbest != NULL ) free(pbest) ;
01234 if( bbest != NULL ) free(bbest) ;
01235 if( pval != NULL ) free(pval) ;
01236 if( bval != NULL ) free(bval) ;
01237 if( paval != NULL ) free(paval) ;
01238 if( avval != NULL ) free(avval) ;
01239 if( pabest!= NULL ) free(pabest);
01240 if( avbest!= NULL ) free(avbest);
01241 if( ptval != NULL ) free(ptval) ;
01242 if( tlval != NULL ) free(tlval) ;
01243 if( sgval != NULL ) free(sgval) ;
01244 if( ptbest!= NULL ) free(ptbest);
01245 if( tlbest!= NULL ) free(tlbest);
01246 if( sgbest!= NULL ) free(sgbest);
01247
01248
01249
01250
01251 ucode_stuff:
01252
01253 #define MAXUFUN 64
01254 #define MAXTS 32
01255
01256 if( ucode != 0 ){
01257 MCW_function_list * rlist = &(GLOBAL_library.registered_fim) ;
01258 int uuse[MAXUFUN] , nbrik[MAXUFUN] , brik1[MAXUFUN] ;
01259 void * udata[MAXUFUN] ;
01260 generic_func * ufunc[MAXUFUN] ;
01261 int nuse , uu , newbrik , oldbrik ;
01262 FIMdata fd ;
01263 MRI_IMAGE * tsim ;
01264 float * tsar , * val , ** vbr ;
01265 short * sar ;
01266 int nts , jts ;
01267 MRI_IMARR * imar ;
01268
01269
01270
01271 for( newbrik=nuse=uu=0 ; uu < rlist->num && nuse < MAXUFUN ; uu++ ){
01272 if( (ucode & (1<<uu)) != 0 ){
01273 uuse [nuse] = uu ;
01274 ufunc[nuse] = rlist->funcs[uu] ;
01275 nbrik[nuse] = rlist->flags[uu] ;
01276 udata[nuse] = rlist->func_data[uu] ;
01277 brik1[nuse] = newbrik ;
01278 newbrik += nbrik[nuse] ;
01279 nuse++ ;
01280 }
01281 }
01282
01283 if( nuse == 0 ) goto final_exit ;
01284
01285
01286
01287 fd.ref_ts = ref_ts ;
01288 fd.ort_ts = ort_ts ;
01289 fd.nvox = nvox ;
01290 fd.ignore = im3d->fimdata->init_ignore ;
01291 fd.polort = polort ;
01292
01293 #ifndef DONT_USE_METER
01294 MCW_set_meter( meter , 0 ) ; meter_pold = 0.0 ;
01295 #endif
01296
01297 for( uu=0 ; uu < nuse ; uu++ )
01298 #if 0
01299 ufunc[uu]( ntime , NULL , udata[uu] , nbrik[uu] , (void *)(&fd) ) ;
01300 #else
01301 AFNI_CALL_fim_function( ufunc[uu] ,
01302 ntime, NULL, udata[uu], nbrik[uu], (&fd) ) ;
01303 #endif
01304
01305
01306
01307
01308
01309
01310 vbr = (float **) malloc(sizeof(float *)*newbrik) ;
01311 for( iv=0 ; iv < newbrik ; iv++ )
01312 vbr[iv] = (float *) malloc(sizeof(float)*nvox) ;
01313
01314 val = (float *) malloc(sizeof(float)*newbrik) ;
01315
01316 for( iv=0 ; iv < nvox ; iv+=MAXTS ){
01317 nts = MIN( MAXTS , nvox-iv ) ;
01318 imar = THD_extract_many_series( nts,indx+iv , dset_time ) ;
01319
01320 for( jts=0 ; jts < nts ; jts++ ){
01321 tsim = IMARR_SUBIMAGE(imar,jts) ;
01322 tsar = MRI_FLOAT_PTR(tsim) ;
01323
01324 for( uu=0 ; uu < nuse ; uu++ ){
01325 #if 0
01326 ufunc[uu]( ntime , tsar ,
01327 udata[uu] , nbrik[uu] , (void *) val ) ;
01328 #else
01329 AFNI_CALL_fim_function( ufunc[uu] ,
01330 ntime, tsar, udata[uu], nbrik[uu], val ) ;
01331 #endif
01332
01333 for( it=0 ; it < nbrik[uu] ; it++ )
01334 vbr[it+brik1[uu]][iv+jts] = val[it] ;
01335 }
01336 }
01337
01338 DESTROY_IMARR(imar) ;
01339
01340 #ifndef DONT_USE_METER
01341 meter_perc = (int) ( 100.0 * iv / nvox ) ;
01342 if( meter_perc != meter_pold ){
01343 MCW_set_meter( meter , meter_perc ) ;
01344 meter_pold = meter_perc ;
01345 }
01346 #endif
01347 }
01348 free(val) ;
01349 #ifndef DONT_USE_METER
01350 MCW_set_meter( meter , 100 ) ;
01351 #endif
01352
01353
01354
01355 if( new_dset != NULL ){
01356 oldbrik = DSET_NVALS(new_dset) ;
01357 } else {
01358 oldbrik = 0 ;
01359
01360 new_dset = EDIT_empty_copy( dset_time ) ;
01361
01362 EDIT_dset_items( new_dset ,
01363 ADN_prefix , new_prefix ,
01364 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
01365 ADN_type , ISHEAD(dset_time)
01366 ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
01367 ADN_func_type , FUNC_BUCK_TYPE ,
01368 ADN_nvals , newbrik ,
01369 ADN_datum_all , MRI_short ,
01370 ADN_ntt , 0 ,
01371 ADN_none ) ;
01372 }
01373
01374
01375
01376
01377
01378
01379 for( iv=0 ; iv < newbrik ; iv++ ){
01380 tsar = vbr[iv] ;
01381 topval = 0.0 ;
01382 for( it=0 ; it < nvox ; it++ )
01383 if( fabs(tsar[it]) > topval ) topval = fabs(tsar[it]) ;
01384
01385 sar = (short *) calloc(sizeof(short),nxyz) ;
01386
01387 if( topval > 0.0 ){
01388 topval = MRI_TYPE_maxval[MRI_short] / topval ;
01389 for( it=0 ; it < nvox ; it++ )
01390 sar[indx[it]] = (short)(topval * tsar[it] + 0.499) ;
01391
01392 topval = 1.0/topval ;
01393 }
01394
01395 free(tsar) ;
01396
01397 if( oldbrik > 0 ){
01398 EDIT_add_brick( new_dset , MRI_short , topval , sar ) ;
01399 } else {
01400 mri_fix_data_pointer( sar , DSET_BRICK(new_dset,iv) ) ;
01401 EDIT_BRICK_FACTOR( new_dset , iv , topval ) ;
01402 }
01403 }
01404 free(vbr) ;
01405
01406
01407
01408 for( uu=0 ; uu < nuse ; uu++ )
01409 #if 0
01410 ufunc[uu]( -(brik1[uu]+oldbrik) , NULL ,
01411 udata[uu] , nbrik[uu] , (void *) new_dset ) ;
01412 #else
01413 AFNI_CALL_fim_function( ufunc[uu] ,
01414 -(brik1[uu]+oldbrik) , NULL ,
01415 udata[uu] , nbrik[uu] , new_dset ) ;
01416 #endif
01417
01418 }
01419
01420 final_exit:
01421 free(vval) ; free(indx) ;
01422
01423
01424
01425 #ifndef DONT_USE_METER
01426 MCW_popdown_meter(meter) ;
01427 #endif
01428
01429 RETURN(new_dset) ;
01430 }