/*-------------------------------------------------------------------
Actually carry out the FIM computations:
im3d
---------------------------------------------------------------------*/
#ifndef FIM_THR
#define FIM_THR 0.0999
#endif
/** Jan 1998: code = combinations of the FIM_*_MASK values **/
/** 30 May 1999: allow polort to be variable (formerly fixed at 1) **/
/** 08 Sep 1999: added PAVE and AVER options **/
/** 03 Jan 2000: added PTOP, TOPL, and SIGM options **/
/** 01 Feb 2000: added ucode, for user written functions **/
THD_3dim_dataset * AFNI_fimmer_compute( Three_D_View * im3d ,
THD_3dim_dataset * dset_time ,
MRI_IMAGE * ref_ts , MRI_IMAGE * ort_ts ,
THD_session * sess , int code, int ucode )
{
THD_3dim_dataset * new_dset=NULL ;
char new_prefix[THD_MAX_PREFIX] ;
char old_prefix[THD_MAX_PREFIX] ;
THD_slist_find fff ;
int ifim , it,iv , nvox , ngood_ref , ntime , it1 , dtyp , nxyz , itbot ;
float * vval, * tsar, * aval, * rbest, * abest, * pbest, * pval, * bbest, * bval;
int * indx ;
short * bar ;
short * ibest ; /* 15 Dec 1997 */
void * ptr ;
float stataux[MAX_STAT_AUX] ;
float fthr , topval ;
int nx_ref , ny_ref , ivec , nnow ;
PCOR_references ** pc_ref ;
PCOR_voxel_corr ** pc_vc ;
int fim_nref , nx_ort , ny_ort , internal_ort ;
float * ortar ;
static float * ref_vec = NULL ;
static int nref_vec = -666 ;
int ibr_best , ibr_perc , ibr_fim , ibr_corr , ibr_base , nbrik ;
int polort = im3d->fimdata->polort , ip ; /* 30 May 1999 */
float top_perc = 0.0 ; /* 30 Aug 1999 */
int ibr_pave , ibr_aver ; /* 08 Sep 1999 */
float * paval , * avval , * pabest , * avbest ; /* 08 Sep 1999 */
int ibr_ptop , ibr_topl , ibr_sigm ; /* 03 Jan 2000 */
float * ptval , * tlval , *sgval ,
* ptbest, * tlbest, *sgbest ;
#ifndef DONT_USE_METER
Widget meter = NULL ;
int meter_perc , meter_pold ;
#endif
int nupdt = 0 , /* number of updates done yet */
min_updt = 5 , /* min number needed for display */
first_updt = 1 ; /* flag to indicate that first update is yet to be displayed */
ENTRY("AFNI_fimmer_compute") ;
/*--- check for legal inputs ---*/
if( ! DSET_GRAPHABLE(dset_time) ||
ref_ts == NULL ||
ref_ts->kind != MRI_float ||
! IM3D_OPEN(im3d) ||
im3d->type != AFNI_3DDATA_VIEW ||
(code == 0 && ucode == 0) || /* Jan 1998 & Feb 2000 */
ref_ts->nx < DSET_NUM_TIMES(dset_time) ){
if(PRINT_TRACING)
{ char str[256] ;
sprintf(str,"illegal inputs: ntime=%d num_ts=%d",
DSET_NUM_TIMES(dset_time), (ref_ts==NULL) ? (0) : (ref_ts->nx) ) ;
STATUS(str) ; }
RETURN(NULL) ;
}
/** 13 Nov 1996: allow for orts **/
if( ort_ts != NULL ){
nx_ort = ort_ts->nx ;
ny_ort = ort_ts->ny ;
ortar = MRI_FLOAT_PTR(ort_ts) ;
internal_ort = (nx_ort < DSET_NUM_TIMES(dset_time)) ;
} else {
internal_ort = 1 ;
}
fim_nref = (internal_ort) ? (polort+2) : (ny_ort+polort+2) ;
if( nref_vec < fim_nref ){
ref_vec = (float *) XtRealloc( (char *)ref_vec, sizeof(float)*fim_nref ) ;
nref_vec = fim_nref ;
}
itbot = im3d->fimdata->init_ignore ;
nx_ref = ref_ts->nx ;
ny_ref = ref_ts->ny ;
ntime = DSET_NUM_TIMES(dset_time) ;
ngood_ref = 0 ;
it1 = -1 ;
for( ivec=0 ; ivec < ny_ref ; ivec++ ){
tsar = MRI_FLOAT_PTR(ref_ts) + (ivec*nx_ref) ;
ifim = 0 ;
for( it=itbot ; it < ntime ; it++ ){
if( tsar[it] < WAY_BIG ){ ifim++ ; if( it1 < 0 ) it1 = it ; }
}
if( ifim < min_updt ){
STATUS("ref_ts has too few good entries!") ;
RETURN(NULL) ;
}
ngood_ref = MAX( ifim , ngood_ref ) ;
}
/** at this point, ngood_ref = max number of good reference points,
and it1 = index of first point used in first reference **/
dtyp = DSET_BRICK_TYPE(dset_time,it1) ;
if( ! AFNI_GOOD_FUNC_DTYPE(dtyp) ){
STATUS("illegal input data type!") ;
RETURN(NULL) ;
}
/*--- Create a new prefix ---*/
MCW_strncpy( old_prefix , DSET_PREFIX(dset_time) , THD_MAX_PREFIX-3 ) ;
if( ! ISVALID_SESSION(sess) ){
sprintf( new_prefix , "%s@%d" , old_prefix , 1 ) ;
} else {
for( ifim=1 ; ifim < 99 ; ifim++ ){
sprintf( new_prefix , "%s@%d" , old_prefix , ifim ) ;
fff = THD_dset_in_session( FIND_PREFIX , new_prefix , sess ) ;
if( fff.dset == NULL ) break ;
}
if( ifim == 99 ){
STATUS("can't create new prefix!") ;
RETURN(NULL) ; /* can't make a new prefix! */
}
}
if(PRINT_TRACING)
{ char str[256] ;
sprintf(str,"new prefix = %s",new_prefix) ; STATUS(str) ; }
/*--- FIM: find values above threshold to fim ---*/
THD_load_datablock( dset_time->dblk , AFNI_purge_unused_dsets ) ;
nxyz = dset_time->dblk->diskptr->dimsizes[0]
* dset_time->dblk->diskptr->dimsizes[1]
* dset_time->dblk->diskptr->dimsizes[2] ;
/** find the mean of the first array,
compute the threshold (fthr) from it,
make indx[i] be the 3D index of the i-th voxel above threshold **/
switch( dtyp ){
case MRI_short:{
short * dar = (short *) DSET_ARRAY(dset_time,it1) ;
for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += abs(dar[iv]) ;
fthr = FIM_THR * fthr / nxyz ;
if(PRINT_TRACING)
{ char str[256] ; sprintf(str,"fthr = %g",fthr) ; STATUS(str) ; }
for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
if( abs(dar[iv]) > fthr ) nvox++ ;
indx = (int *) malloc( sizeof(int) * nvox ) ;
if( indx == NULL ){
fprintf(stderr,"\n*** indx malloc failure in AFNI_fimmer_compute\n") ;
RETURN(NULL) ;
}
for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
if( abs(dar[iv]) > fthr ) indx[nvox++] = iv ;
}
break ;
case MRI_float:{
float * dar = (float *) DSET_ARRAY(dset_time,it1) ;
for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += fabs(dar[iv]) ;
fthr = FIM_THR * fthr / nxyz ;
if(PRINT_TRACING)
{ char str[256] ; sprintf(str,"fthr = %g",fthr) ; STATUS(str) ; }
for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
if( fabs(dar[iv]) > fthr ) nvox++ ;
indx = (int *) malloc( sizeof(int) * nvox ) ;
if( indx == NULL ){
fprintf(stderr,"\n*** indx malloc failure in AFNI_fimmer_compute\n") ;
RETURN(NULL) ;
}
for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
if( fabs(dar[iv]) > fthr ) indx[nvox++] = iv ;
}
break ;
case MRI_byte:{
byte * dar = (byte *) DSET_ARRAY(dset_time,it1) ;
for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += dar[iv] ;
fthr = FIM_THR * fthr / nxyz ;
if(PRINT_TRACING)
{ char str[256] ; sprintf(str,"fthr = %g",fthr) ; STATUS(str) ; }
for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
if( dar[iv] > fthr ) nvox++ ;
indx = (int *) malloc( sizeof(int) * nvox ) ;
if( indx == NULL ){
fprintf(stderr,"\n*** indx malloc failure in AFNI_fimmer_compute\n") ;
RETURN(NULL) ;
}
for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
if( dar[iv] > fthr ) indx[nvox++] = iv ;
}
break ;
}
if(PRINT_TRACING)
{ char str[256] ; sprintf(str,"number of voxels = %d",nvox) ; STATUS(str) ; }
/** allocate space for voxel values **/
vval = (float *) malloc( sizeof(float) * nvox) ;
if( vval == NULL ){
fprintf(stderr,"\n*** vval malloc failure in AFNI_fimmer_compute\n") ;
free(indx) ; RETURN(NULL) ;
}
/** compute number of output bricks **/
ibr_fim=ibr_corr=ibr_best=ibr_perc=ibr_base = -1 ; nbrik = 0 ;
ibr_pave = ibr_aver = -1 ;
ibr_ptop = ibr_topl = ibr_sigm = -1 ;
if( (code & FIM_ALPHA_MASK)!= 0) { ibr_fim = nbrik; nbrik++; }
if( (code & FIM_BEST_MASK) != 0 && ny_ref > 1){ ibr_best = nbrik; nbrik++; }
if( (code & FIM_PERC_MASK) != 0) { ibr_perc = nbrik; nbrik++; }
if( (code & FIM_PAVE_MASK) != 0) { ibr_pave = nbrik; nbrik++; }
if( (code & FIM_BASE_MASK) != 0) { ibr_base = nbrik; nbrik++; }
if( (code & FIM_AVER_MASK) != 0) { ibr_aver = nbrik; nbrik++; }
if( (code & FIM_CORR_MASK) != 0) { ibr_corr = nbrik; nbrik++; }
if( (code & FIM_PTOP_MASK) != 0) { ibr_ptop = nbrik; nbrik++; }
if( (code & FIM_TOPL_MASK) != 0) { ibr_topl = nbrik; nbrik++; }
if( (code & FIM_SIGM_MASK) != 0) { ibr_sigm = nbrik; nbrik++; }
/** 01 Feb 2000: if no normal FIM stuff (code), skip to the ucode stuff **/
if(PRINT_TRACING)
{ char str[256] ; sprintf(str,"number of bricks = %d",nbrik) ; STATUS(str) ; }
if( nbrik == 0 ){
#ifndef DONT_USE_METER
meter = MCW_popup_meter( im3d->vwid->top_shell , METER_TOP_WIDE ) ;
meter_pold = 0 ;
#endif
goto ucode_stuff ; /* way below */
}
/** normal case: do the normal recursive FIMming **/
if(PRINT_TRACING)
{ char str[256] ;
sprintf(str,"code of FIM_MASKs = %d",code) ; STATUS(str) ; }
/** allocate extra space for comparing results from multiple ref vectors **/
if( ny_ref > 1 ){
aval = (float *) malloc(sizeof(float) * nvox) ;
rbest = (float *) malloc(sizeof(float) * nvox) ;
abest = (float *) malloc(sizeof(float) * nvox) ;
ibest = (short *) malloc(sizeof(short) * nvox) ; /* 15 Dec 1997 */
pbest = (float *) malloc(sizeof(float) * nvox) ; /* 16 Jan 1998 */
bbest = (float *) malloc(sizeof(float) * nvox) ; /* 16 Jan 1998 */
pval = (float *) malloc(sizeof(float) * nvox) ; /* 16 Jan 1998 */
bval = (float *) malloc(sizeof(float) * nvox) ; /* 16 Jan 1998 */
paval = (float *) malloc(sizeof(float) * nvox) ; /* 08 Sep 1999 */
avval = (float *) malloc(sizeof(float) * nvox) ; /* 08 Sep 1999 */
pabest = (float *) malloc(sizeof(float) * nvox) ; /* 16 Jan 1998 */
avbest = (float *) malloc(sizeof(float) * nvox) ; /* 16 Jan 1998 */
ptval = (float *) malloc(sizeof(float) * nvox) ; /* 03 Jan 2000 */
tlval = (float *) malloc(sizeof(float) * nvox) ; /* 03 Jan 2000 */
sgval = (float *) malloc(sizeof(float) * nvox) ; /* 03 Jan 2000 */
ptbest = (float *) malloc(sizeof(float) * nvox) ; /* 03 Jan 2000 */
tlbest = (float *) malloc(sizeof(float) * nvox) ; /* 03 Jan 2000 */
sgbest = (float *) malloc(sizeof(float) * nvox) ; /* 03 Jan 2000 */
if( sgbest == NULL ){
fprintf(stderr,"\n*** 'best' malloc failure in AFNI_fimmer_compute\n") ;
free(vval) ; free(indx) ;
if( aval != NULL ) free(aval) ;
if( rbest != NULL ) free(rbest) ;
if( abest != NULL ) free(abest) ;
if( ibest != NULL ) free(ibest) ; /* 15 Dec 1997 */
if( pbest != NULL ) free(pbest) ; /* 16 Jan 1998 */
if( bbest != NULL ) free(bbest) ; /* 16 Jan 1998 */
if( pval != NULL ) free(pval) ; /* 16 Jan 1998 */
if( bval != NULL ) free(bval) ; /* 16 Jan 1998 */
if( paval != NULL ) free(paval) ; /* 08 Sep 1999 */
if( avval != NULL ) free(avval) ; /* 08 Sep 1999 */
if( pabest!= NULL ) free(pabest); /* 08 Sep 1999 */
if( avbest!= NULL ) free(avbest); /* 08 Sep 1999 */
if( ptval != NULL ) free(ptval) ; /* 03 Jan 2000 */
if( tlval != NULL ) free(tlval) ; /* 03 Jan 2000 */
if( sgval != NULL ) free(sgval) ; /* 03 Jan 2000 */
if( ptbest!= NULL ) free(ptbest); /* 03 Jan 2000 */
if( tlbest!= NULL ) free(tlbest); /* 03 Jan 2000 */
if( sgbest!= NULL ) free(sgbest); /* 03 Jan 2000 */
RETURN(NULL) ;
}
} else {
aval = rbest = abest = pbest = bbest = pval = bval = NULL ;
paval = avval = pabest = avbest = NULL ; /* 08 Sep 1999 */
ptval = tlval = ptbest = tlbest = NULL ; /* 03 Jan 2000 */
sgval = sgbest = NULL ; /* 03 Jan 2000 */
ibest = NULL ; /* 15 Dec 1997 */
}
if(PRINT_TRACING)
{ char str[256] ;
sprintf(str,"nxyz = %d nvox = %d",nxyz,nvox) ; STATUS(str) ; }
/*--- FIM: initialize recursive updates ---*/
pc_ref = (PCOR_references **) malloc( sizeof(PCOR_references *) * ny_ref ) ;
pc_vc = (PCOR_voxel_corr **) malloc( sizeof(PCOR_voxel_corr *) * ny_ref ) ;
if( pc_ref == NULL || pc_vc == NULL ){
free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
if( aval != NULL ) free(aval) ;
if( rbest != NULL ) free(rbest) ;
if( abest != NULL ) free(abest) ;
if( ibest != NULL ) free(ibest) ; /* 15 Dec 1997 */
if( pbest != NULL ) free(pbest) ; /* 16 Jan 1998 */
if( bbest != NULL ) free(bbest) ; /* 16 Jan 1998 */
if( pval != NULL ) free(pval) ; /* 16 Jan 1998 */
if( bval != NULL ) free(bval) ; /* 16 Jan 1998 */
if( paval != NULL ) free(paval) ; /* 08 Sep 1999 */
if( avval != NULL ) free(avval) ; /* 08 Sep 1999 */
if( pabest!= NULL ) free(pabest); /* 08 Sep 1999 */
if( avbest!= NULL ) free(avbest); /* 08 Sep 1999 */
if( ptval != NULL ) free(ptval) ; /* 03 Jan 2000 */
if( tlval != NULL ) free(tlval) ; /* 03 Jan 2000 */
if( sgval != NULL ) free(sgval) ; /* 03 Jan 2000 */
if( ptbest!= NULL ) free(ptbest); /* 03 Jan 2000 */
if( tlbest!= NULL ) free(tlbest); /* 03 Jan 2000 */
if( sgbest!= NULL ) free(sgbest); /* 03 Jan 2000 */
fprintf(stderr,"\n*** FIM initialization fails in AFNI_fimmer_compute\n") ;
RETURN(NULL) ;
}
ifim = 0 ;
for( ivec=0 ; ivec < ny_ref ; ivec++ ){
pc_ref[ivec] = new_PCOR_references( fim_nref ) ;
pc_vc[ivec] = new_PCOR_voxel_corr( nvox , fim_nref ) ;
if( pc_ref[ivec] == NULL || pc_vc[ivec] == NULL ) ifim++ ;
}
if( ifim > 0 ){
for( ivec=0 ; ivec < ny_ref ; ivec++ ){
free_PCOR_references(pc_ref[ivec]) ;
free_PCOR_voxel_corr(pc_vc[ivec]) ;
}
free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
if( aval != NULL ) free(aval) ;
if( rbest != NULL ) free(rbest) ;
if( abest != NULL ) free(abest) ;
if( ibest != NULL ) free(ibest) ; /* 15 Dec 1997 */
if( pbest != NULL ) free(pbest) ; /* 16 Jan 1998 */
if( bbest != NULL ) free(bbest) ; /* 16 Jan 1998 */
if( pval != NULL ) free(pval) ; /* 16 Jan 1998 */
if( bval != NULL ) free(bval) ; /* 16 Jan 1998 */
if( paval != NULL ) free(paval) ; /* 08 Sep 1999 */
if( avval != NULL ) free(avval) ; /* 08 Sep 1999 */
if( pabest!= NULL ) free(pabest); /* 08 Sep 1999 */
if( avbest!= NULL ) free(avbest); /* 08 Sep 1999 */
if( ptval != NULL ) free(ptval) ; /* 03 Jan 2000 */
if( tlval != NULL ) free(tlval) ; /* 03 Jan 2000 */
if( sgval != NULL ) free(sgval) ; /* 03 Jan 2000 */
if( ptbest!= NULL ) free(ptbest); /* 03 Jan 2000 */
if( tlbest!= NULL ) free(tlbest); /* 03 Jan 2000 */
if( sgbest!= NULL ) free(sgbest); /* 03 Jan 2000 */
fprintf(stderr,"\n*** FIM initialization fails in AFNI_fimmer_compute\n") ;
RETURN(NULL) ;
}
/*--- Make a new dataset to hold the output ---*/
STATUS("making new dataset") ;
new_dset = EDIT_empty_copy( dset_time ) ;
if( nbrik == 1 && ucode == 0 ){ /* 1 brick out --> a 'fim' dataset */
it = EDIT_dset_items( new_dset ,
ADN_prefix , new_prefix ,
ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
ADN_type , ISHEAD(dset_time)
? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
ADN_func_type , FUNC_FIM_TYPE ,
ADN_nvals , 1 ,
ADN_datum_all , MRI_short ,
ADN_ntt , 0 ,
ADN_none ) ;
/* 2 bricks, 2nd corr --> 'fico' */
} else if( nbrik == 2 && ibr_corr == 1 && ucode == 0 ){
it = EDIT_dset_items( new_dset ,
ADN_prefix , new_prefix ,
ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
ADN_type , ISHEAD(dset_time)
? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
ADN_func_type , FUNC_COR_TYPE ,
ADN_nvals , 2 ,
ADN_datum_all , MRI_short ,
ADN_ntt , 0 ,
ADN_none ) ;
} else if( nbrik > 0 ){ /* otherwise --> 'fbuc' (bucket) */
it = EDIT_dset_items( new_dset ,
ADN_prefix , new_prefix ,
ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
ADN_type , ISHEAD(dset_time)
? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
ADN_func_type , FUNC_BUCK_TYPE ,
ADN_nvals , nbrik ,
ADN_datum_all , MRI_short ,
ADN_ntt , 0 ,
ADN_none ) ;
} else {
it = 999 ;
}
if( it > 0 ){
fprintf(stderr,
"\n*** EDIT_dset_items error %d in AFNI_fimmer_compute\n",it) ;
THD_delete_3dim_dataset( new_dset , False ) ;
for( ivec=0 ; ivec < ny_ref ; ivec++ ){
free_PCOR_references(pc_ref[ivec]) ;
free_PCOR_voxel_corr(pc_vc[ivec]) ;
}
free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
if( aval != NULL ) free(aval) ;
if( rbest != NULL ) free(rbest) ;
if( abest != NULL ) free(abest) ;
if( ibest != NULL ) free(ibest) ; /* 15 Dec 1997 */
if( pbest != NULL ) free(pbest) ; /* 16 Jan 1998 */
if( bbest != NULL ) free(bbest) ; /* 16 Jan 1998 */
if( pval != NULL ) free(pval) ; /* 16 Jan 1998 */
if( bval != NULL ) free(bval) ; /* 16 Jan 1998 */
if( paval != NULL ) free(paval) ; /* 08 Sep 1999 */
if( avval != NULL ) free(avval) ; /* 08 Sep 1999 */
if( pabest!= NULL ) free(pabest); /* 08 Sep 1999 */
if( avbest!= NULL ) free(avbest); /* 08 Sep 1999 */
if( ptval != NULL ) free(ptval) ; /* 03 Jan 2000 */
if( tlval != NULL ) free(tlval) ; /* 03 Jan 2000 */
if( sgval != NULL ) free(sgval) ; /* 03 Jan 2000 */
if( ptbest!= NULL ) free(ptbest); /* 03 Jan 2000 */
if( tlbest!= NULL ) free(tlbest); /* 03 Jan 2000 */
if( sgbest!= NULL ) free(sgbest); /* 03 Jan 2000 */
RETURN(NULL) ;
}
/* modify labels for each brick */
if( ibr_fim >= 0 )
EDIT_BRICK_LABEL( new_dset , ibr_fim , "Fit Coef" ) ;
if( ibr_corr >= 0 )
EDIT_BRICK_LABEL( new_dset , ibr_corr , "Correlation" ) ;
if( ibr_best >= 0 )
EDIT_BRICK_LABEL( new_dset , ibr_best , "Best Index" ) ;
if( ibr_perc >= 0 )
EDIT_BRICK_LABEL( new_dset , ibr_perc , "% Change" ) ;
if( ibr_base >= 0 )
EDIT_BRICK_LABEL( new_dset , ibr_base , "Baseline" ) ;
if( ibr_pave >= 0 )
EDIT_BRICK_LABEL( new_dset , ibr_pave , "% From Ave" ) ; /* 08 Sep 1999 */
if( ibr_aver >= 0 )
EDIT_BRICK_LABEL( new_dset , ibr_aver , "Average" ) ;
if( ibr_ptop >= 0 )
EDIT_BRICK_LABEL( new_dset , ibr_ptop , "% From Top" ) ; /* 03 Jan 2000 */
if( ibr_topl >= 0 )
EDIT_BRICK_LABEL( new_dset , ibr_topl , "Topline" ) ;
if( ibr_sigm >= 0 )
EDIT_BRICK_LABEL( new_dset , ibr_sigm , "Sigma Resid" ) ;
/*-- 30 Aug 1999: set limits on percent change --*/
if( ibr_perc >= 0 || ibr_pave >= 0 || ibr_ptop >= 0 ){
char * cp = my_getenv("AFNI_FIM_PERCENT_LIMIT") ;
if( cp != NULL ){
float tp = strtod(cp,NULL) ;
if( tp > 0.0 ) top_perc = tp ;
}
}
/* create bricks */
STATUS("making output bricks") ;
for( iv=0 ; iv < new_dset->dblk->nvals ; iv++ ){
ptr = malloc( DSET_BRICK_BYTES(new_dset,iv) ) ;
mri_fix_data_pointer( ptr , DSET_BRICK(new_dset,iv) ) ;
}
if( THD_count_databricks(new_dset->dblk) < new_dset->dblk->nvals ){
fprintf(stderr,
"\n*** failure to malloc new bricks in AFNI_fimmer_compute\n") ;
THD_delete_3dim_dataset( new_dset , False ) ;
for( ivec=0 ; ivec < ny_ref ; ivec++ ){
free_PCOR_references(pc_ref[ivec]) ;
free_PCOR_voxel_corr(pc_vc[ivec]) ;
}
free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
if( aval != NULL ) free(aval) ;
if( rbest != NULL ) free(rbest) ;
if( abest != NULL ) free(abest) ;
if( ibest != NULL ) free(ibest) ; /* 15 Dec 1997 */
if( pbest != NULL ) free(pbest) ; /* 16 Jan 1998 */
if( bbest != NULL ) free(bbest) ; /* 16 Jan 1998 */
if( pval != NULL ) free(pval) ; /* 16 Jan 1998 */
if( bval != NULL ) free(bval) ; /* 16 Jan 1998 */
if( paval != NULL ) free(paval) ; /* 08 Sep 1999 */
if( avval != NULL ) free(avval) ; /* 08 Sep 1999 */
if( pabest!= NULL ) free(pabest); /* 08 Sep 1999 */
if( avbest!= NULL ) free(avbest); /* 08 Sep 1999 */
if( ptval != NULL ) free(ptval) ; /* 03 Jan 2000 */
if( tlval != NULL ) free(tlval) ; /* 03 Jan 2000 */
if( sgval != NULL ) free(sgval) ; /* 03 Jan 2000 */
if( ptbest!= NULL ) free(ptbest); /* 03 Jan 2000 */
if( tlbest!= NULL ) free(tlbest); /* 03 Jan 2000 */
if( sgbest!= NULL ) free(sgbest); /* 03 Jan 2000 */
RETURN(NULL) ;
}
/*---------------------------------*/
/*--- FIM: do recursive updates ---*/
#ifndef DONT_USE_METER
meter = MCW_popup_meter( im3d->vwid->top_shell , METER_TOP_WIDE ) ;
meter_pold = 0 ;
#endif
STATUS("starting recursive least squares") ;
for( it=itbot ; it < ntime ; it++ ){ /* loop over time */
nnow = 0 ; /* number of updates done at this time point */
for( ivec=0 ; ivec < ny_ref ; ivec++ ){ /* loop over ref vects */
tsar = MRI_FLOAT_PTR(ref_ts) + (ivec*nx_ref) ; /* ptr to vect */
if( tsar[it] >= WAY_BIG ) continue ; /* skip this */
ref_vec[0] = 1.0 ; /* we always supply ort for constant */
for( ip=1 ; ip <= polort ; ip++ ) /* 30 May 1999: */
ref_vec[ip] = ref_vec[ip-1] * ((float)it) ; /* and polynomials */
if( internal_ort ){ /* no external orts */
ref_vec[ip] = tsar[it] ; /* ref value */
} else {
for( iv=0 ; iv < ny_ort ; iv++ ) /* external */
ref_vec[iv+ip] = ortar[it + iv*nx_ort] ; /* orts */
ref_vec[ny_ort+ip] = tsar[it] ; /* ref value */
}
if(PRINT_TRACING)
{ char str[256] ;
sprintf(str,"time index=%d ideal[%d]=%f" , it,ivec,tsar[it] ) ;
STATUS(str) ; }
/* process the ort+ref update */
update_PCOR_references( ref_vec , pc_ref[ivec] ) ;
/* first time thru: load data from dataset */
if( nnow == 0 ){
switch( dtyp ){
case MRI_short:{
short * dar = (short *) DSET_ARRAY(dset_time,it) ;
for( iv=0; iv < nvox; iv++ ) vval[iv] = (float) dar[indx[iv]];
}
break ;
case MRI_float:{
float * dar = (float *) DSET_ARRAY(dset_time,it) ;
for( iv=0; iv < nvox; iv++ ) vval[iv] = (float) dar[indx[iv]];
}
break ;
case MRI_byte:{
byte * dar = (byte *) DSET_ARRAY(dset_time,it) ;
for( iv=0; iv < nvox; iv++ ) vval[iv] = (float) dar[indx[iv]];
}
break ;
}
}
/* process the data update */
PCOR_update_float( vval , pc_ref[ivec] , pc_vc[ivec] ) ;
nnow++ ; /* one more update at this time point */
} /* end of loop over ref vects */
if( nnow > 0 ) nupdt++ ; /* number of time points that had updates */
#ifndef DONT_USE_METER
meter_perc = (int) ( 100.0 * nupdt / ngood_ref ) ;
if( meter_perc != meter_pold ){
MCW_set_meter( meter , meter_perc ) ;
meter_pold = meter_perc ;
}
#endif
} /* end of loop over time */
/*-------------------------------------------*/
/*--- Load final results into the dataset ---*/
/*--- set the statistical parameters ---*/
stataux[0] = nupdt ; /* number of points used */
stataux[1] = (ny_ref==1) ? 1 : 2 ; /* number of references */
stataux[2] = fim_nref - 1 ; /* number of orts */
for( iv=3 ; iv < MAX_STAT_AUX ; iv++ ) stataux[iv] = 0.0 ;
if( ibr_corr >= 0 ){
if( new_dset->func_type == FUNC_COR_TYPE )
EDIT_dset_items( new_dset, ADN_stat_aux, stataux, ADN_none ) ;
EDIT_BRICK_TO_FICO( new_dset, ibr_corr, stataux[0],stataux[1],stataux[2] ) ;
}
#ifndef DONT_USE_METER
# define METERIZE(ib) do { meter_perc = (int) ( 100.0 * (ib) / nbrik ) ; \
if( meter_perc != meter_pold ){ \
MCW_set_meter( meter , meter_perc ) ; \
meter_pold = meter_perc ; \
} } while(0)
#else
# define METERIZE(ib) /*nada*/
#endif
/*** Compute brick arrays for new dataset ***/
/* [load scale factors into stataux, too] */
if( ny_ref == 1 ){
/*** Just 1 ref vector --> load values directly into dataset ***/
if( ibr_fim >= 0 ){
STATUS("getting 1 ref alpha") ;
PCOR_get_coef( pc_ref[0] , pc_vc[0] , vval ) ;
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_fim ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
stataux[ibr_fim] = 1.0/topval ;
} else {
stataux[ibr_fim] = 0.0 ;
}
METERIZE(ibr_fim) ;
}
if( ibr_corr >= 0 ){
STATUS("getting 1 ref pcor") ;
PCOR_get_pcor( pc_ref[0] , pc_vc[0] , vval ) ;
bar = DSET_ARRAY( new_dset , ibr_corr ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(FUNC_COR_SCALE_SHORT * vval[iv] + 0.499) ;
stataux[ibr_corr] = 1.0 / FUNC_COR_SCALE_SHORT ;
METERIZE(ibr_corr) ;
}
if( ibr_perc >= 0 ){
STATUS("getting 1 ref perc") ;
PCOR_get_perc( pc_ref[0] , pc_vc[0] , vval , NULL , 0 ) ;
if( top_perc > 0.0 ) EDIT_clip_float( top_perc , nvox , vval ) ;
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_perc ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
stataux[ibr_perc] = 1.0/topval ;
} else {
stataux[ibr_perc] = 0.0 ;
}
METERIZE(ibr_perc) ;
}
if( ibr_pave >= 0 ){ /* 08 Sep 1999 */
STATUS("getting 1 ref pave") ;
PCOR_get_perc( pc_ref[0] , pc_vc[0] , vval , NULL , 1 ) ;
if( top_perc > 0.0 ) EDIT_clip_float( top_perc , nvox , vval ) ;
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_pave ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
stataux[ibr_pave] = 1.0/topval ;
} else {
stataux[ibr_pave] = 0.0 ;
}
METERIZE(ibr_pave) ;
}
if( ibr_ptop >= 0 ){ /* 03 Jan 2000 */
STATUS("getting 1 ref ptop") ;
PCOR_get_perc( pc_ref[0] , pc_vc[0] , vval , NULL , 2 ) ;
if( top_perc > 0.0 ) EDIT_clip_float( top_perc , nvox , vval ) ;
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_ptop ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
stataux[ibr_ptop] = 1.0/topval ;
} else {
stataux[ibr_ptop] = 0.0 ;
}
METERIZE(ibr_ptop) ;
}
if( ibr_base >= 0 ){
STATUS("getting 1 ref base") ;
PCOR_get_perc( pc_ref[0] , pc_vc[0] , NULL , vval , 0 ) ;
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_base ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
stataux[ibr_base] = 1.0/topval ;
} else {
stataux[ibr_base] = 0.0 ;
}
METERIZE(ibr_base) ;
}
if( ibr_aver >= 0 ){ /* 08 Sep 1999 */
STATUS("getting 1 ref aver") ;
PCOR_get_perc( pc_ref[0] , pc_vc[0] , NULL , vval , 1 ) ;
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_aver ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
stataux[ibr_aver] = 1.0/topval ;
} else {
stataux[ibr_aver] = 0.0 ;
}
METERIZE(ibr_aver) ;
}
if( ibr_topl >= 0 ){ /* 03 Jan 2000 */
STATUS("getting 1 ref topl") ;
PCOR_get_perc( pc_ref[0] , pc_vc[0] , NULL , vval , 2 ) ;
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_topl ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
stataux[ibr_topl] = 1.0/topval ;
} else {
stataux[ibr_topl] = 0.0 ;
}
METERIZE(ibr_topl) ;
}
if( ibr_sigm >= 0 ){ /* 03 Jan 2000 */
STATUS("getting 1 ref sigm") ;
PCOR_get_stdev( pc_vc[0] , vval ) ;
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_sigm ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
stataux[ibr_sigm] = 1.0/topval ;
} else {
stataux[ibr_sigm] = 0.0 ;
}
METERIZE(ibr_sigm) ;
}
} else {
/*** Multiple references --> find best correlation at each voxel ***/
/*--- get first ref results into abest and rbest (best so far) ---*/
STATUS("getting first ref results") ;
PCOR_get_coef( pc_ref[0] , pc_vc[0] , abest ) ;
PCOR_get_pcor( pc_ref[0] , pc_vc[0] , rbest ) ;
PCOR_get_perc( pc_ref[0] , pc_vc[0] , pbest , bbest , 0 ) ;
PCOR_get_perc( pc_ref[0] , pc_vc[0] , pabest, avbest, 1 ) ;
PCOR_get_perc( pc_ref[0] , pc_vc[0] , ptbest, tlbest, 2 ) ;
PCOR_get_stdev( pc_vc[0] , sgbest ) ;
for( iv=0 ; iv < nvox ; iv++ ) ibest[iv] = 1 ; /* 15 Dec 1997 */
/*--- for each succeeding ref vector,
get results into aval and vval,
if |vval| > |rbest|, then use that result instead ---*/
for( ivec=1 ; ivec < ny_ref ; ivec++ ){
STATUS(" == getting results for next ref") ;
PCOR_get_coef( pc_ref[ivec] , pc_vc[ivec] , aval ) ;
PCOR_get_pcor( pc_ref[ivec] , pc_vc[ivec] , vval ) ;
PCOR_get_perc( pc_ref[ivec] , pc_vc[ivec] , pval , bval , 0 ) ;
PCOR_get_perc( pc_ref[ivec] , pc_vc[ivec] , paval, avval, 1 ) ;
PCOR_get_perc( pc_ref[ivec] , pc_vc[ivec] , ptval, tlval, 2 ) ;
PCOR_get_stdev( pc_vc[ivec] , sgval ) ;
STATUS(" == and finding the best results") ;
for( iv=0 ; iv < nvox ; iv++ ){
if( fabs(vval[iv]) > fabs(rbest[iv]) ){
rbest[iv] = vval[iv] ;
abest[iv] = aval[iv] ;
ibest[iv] = (ivec+1) ; /* 15 Dec 1997 */
pbest[iv] = pval[iv] ; /* Jan 1998 */
bbest[iv] = bval[iv] ;
pabest[iv]= paval[iv] ; /* 08 Sep 1999 */
avbest[iv]= avval[iv] ;
ptbest[iv]= ptval[iv] ; /* 03 Jan 1999 */
tlbest[iv]= tlval[iv] ;
sgbest[iv]= sgval[iv] ;
}
}
}
/*--- at this point, abest and rbest are the best
results, so scale them into the dataset bricks ---*/
/** fim brick **/
if( ibr_fim >= 0 ){
if(PRINT_TRACING)
{ char str[256]; sprintf(str,"getting ibr_fim=%d",ibr_fim); STATUS(str); }
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(abest[iv]) > topval ) topval = fabs(abest[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_fim ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * abest[iv] + 0.499) ;
stataux[ibr_fim] = 1.0/topval ;
} else {
stataux[ibr_fim] = 0.0 ;
}
METERIZE(ibr_fim) ;
}
/** threshold brick **/
if( ibr_corr >= 0 ){
if(PRINT_TRACING)
{ char str[256]; sprintf(str,"getting ibr_corr=%d",ibr_corr); STATUS(str); }
bar = DSET_ARRAY( new_dset , ibr_corr ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(FUNC_COR_SCALE_SHORT * rbest[iv] + 0.499) ;
stataux[ibr_corr] = 1.0 / FUNC_COR_SCALE_SHORT ;
METERIZE(ibr_corr) ;
}
/** best index brick (15 Dec 1997) */
if( ibr_best >= 0 ){
if(PRINT_TRACING)
{ char str[256]; sprintf(str,"getting ibr_best=%d",ibr_best); STATUS(str); }
bar = DSET_ARRAY( new_dset , ibr_best ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
for( iv=0 ; iv < nvox ; iv++ ) bar[indx[iv]] = ibest[iv] ;
stataux[ibr_best] = 0.0 ; /* no scaling */
METERIZE(ibr_best) ;
}
/** perc brick */
if( ibr_perc >= 0 ){
if(PRINT_TRACING)
{ char str[256]; sprintf(str,"getting ibr_perc=%d",ibr_perc); STATUS(str); }
if( top_perc > 0.0 ) EDIT_clip_float( top_perc , nvox , pbest ) ;
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(pbest[iv]) > topval ) topval = fabs(pbest[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_perc ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * pbest[iv] + 0.499) ;
stataux[ibr_perc] = 1.0/topval ;
} else {
stataux[ibr_perc] = 0.0 ;
}
METERIZE(ibr_perc) ;
}
/** pave brick [08 Sep 1999] */
if( ibr_pave >= 0 ){
if(PRINT_TRACING)
{ char str[256]; sprintf(str,"getting ibr_pave=%d",ibr_pave); STATUS(str); }
if( top_perc > 0.0 ) EDIT_clip_float( top_perc , nvox , pabest ) ;
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(pabest[iv]) > topval ) topval = fabs(pabest[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_pave ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * pabest[iv] + 0.499) ;
stataux[ibr_pave] = 1.0/topval ;
} else {
stataux[ibr_pave] = 0.0 ;
}
METERIZE(ibr_pave) ;
}
/** ptop brick [03 Jan 2000] */
if( ibr_ptop >= 0 ){
if(PRINT_TRACING)
{ char str[256]; sprintf(str,"getting ibr_ptop=%d",ibr_ptop); STATUS(str); }
if( top_perc > 0.0 ) EDIT_clip_float( top_perc , nvox , ptbest ) ;
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(ptbest[iv]) > topval ) topval = fabs(ptbest[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_ptop ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * ptbest[iv] + 0.499) ;
stataux[ibr_ptop] = 1.0/topval ;
} else {
stataux[ibr_ptop] = 0.0 ;
}
METERIZE(ibr_ptop) ;
}
/** base brick */
if( ibr_base >= 0 ){
if(PRINT_TRACING)
{ char str[256]; sprintf(str,"getting ibr_base=%d",ibr_base); STATUS(str); }
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(bbest[iv]) > topval ) topval = fabs(bbest[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_base ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * bbest[iv] + 0.499) ;
stataux[ibr_base] = 1.0/topval ;
} else {
stataux[ibr_base] = 0.0 ;
}
METERIZE(ibr_base) ;
}
/** aver brick [08 Sep 1999] */
if( ibr_aver >= 0 ){
if(PRINT_TRACING)
{ char str[256]; sprintf(str,"getting ibr_aver=%d",ibr_aver); STATUS(str); }
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(avbest[iv]) > topval ) topval = fabs(avbest[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_aver ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * avbest[iv] + 0.499) ;
stataux[ibr_aver] = 1.0/topval ;
} else {
stataux[ibr_aver] = 0.0 ;
}
METERIZE(ibr_aver) ;
}
/** topl brick [03 Jan 2000] */
if( ibr_topl >= 0 ){
if(PRINT_TRACING)
{ char str[256]; sprintf(str,"getting ibr_topl=%d",ibr_topl); STATUS(str); }
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(tlbest[iv]) > topval ) topval = fabs(tlbest[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_topl ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * tlbest[iv] + 0.499) ;
stataux[ibr_topl] = 1.0/topval ;
} else {
stataux[ibr_topl] = 0.0 ;
}
METERIZE(ibr_topl) ;
}
/** sigm brick [03 Jan 2000]**/
if( ibr_sigm >= 0 ){
if(PRINT_TRACING)
{ char str[256]; sprintf(str,"getting ibr_sigm=%d",ibr_sigm); STATUS(str); }
topval = 0.0 ;
for( iv=0 ; iv < nvox ; iv++ )
if( fabs(sgbest[iv]) > topval ) topval = fabs(sgbest[iv]) ;
bar = DSET_ARRAY( new_dset , ibr_sigm ) ;
memset( bar , 0 , sizeof(short)*nxyz ) ;
if( topval > 0.0 ){
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( iv=0 ; iv < nvox ; iv++ )
bar[indx[iv]] = (short)(topval * sgbest[iv] + 0.499) ;
stataux[ibr_sigm] = 1.0/topval ;
} else {
stataux[ibr_sigm] = 0.0 ;
}
METERIZE(ibr_sigm) ;
}
} /* end of multiple reference case */
/*** Set the brick factors for the new dataset,
no matter how it was computed above. ***/
STATUS("setting brick_fac") ;
(void) EDIT_dset_items( new_dset , ADN_brick_fac , stataux , ADN_none ) ;
#ifndef DONT_USE_METER
MCW_set_meter( meter , 100 ) ;
#endif
/*--- End of recursive updates; now free temporary workspaces ---*/
for( ivec=0 ; ivec < ny_ref ; ivec++ ){
free_PCOR_references(pc_ref[ivec]) ;
free_PCOR_voxel_corr(pc_vc[ivec]) ;
}
free(pc_ref) ; free(pc_vc) ;
if( aval != NULL ) free(aval) ;
if( rbest != NULL ) free(rbest) ;
if( abest != NULL ) free(abest) ;
if( ibest != NULL ) free(ibest) ; /* 15 Dec 1997 */
if( pbest != NULL ) free(pbest) ; /* 16 Jan 1998 */
if( bbest != NULL ) free(bbest) ; /* 16 Jan 1998 */
if( pval != NULL ) free(pval) ; /* 16 Jan 1998 */
if( bval != NULL ) free(bval) ; /* 16 Jan 1998 */
if( paval != NULL ) free(paval) ; /* 08 Sep 1999 */
if( avval != NULL ) free(avval) ; /* 08 Sep 1999 */
if( pabest!= NULL ) free(pabest); /* 08 Sep 1999 */
if( avbest!= NULL ) free(avbest); /* 08 Sep 1999 */
if( ptval != NULL ) free(ptval) ; /* 03 Jan 2000 */
if( tlval != NULL ) free(tlval) ; /* 03 Jan 2000 */
if( sgval != NULL ) free(sgval) ; /* 03 Jan 2000 */
if( ptbest!= NULL ) free(ptbest); /* 03 Jan 2000 */
if( tlbest!= NULL ) free(tlbest); /* 03 Jan 2000 */
if( sgbest!= NULL ) free(sgbest); /* 03 Jan 2000 */
/*-----------------------------------------------------*/
/*--- 01 Feb 2000: execute user specified functions ---*/
ucode_stuff:
#define MAXUFUN 64 /* should be at least sizeof(int) */
#define MAXTS 32 /* number of timeseries to process at once */
if( ucode != 0 ){
MCW_function_list * rlist = &(GLOBAL_library.registered_fim) ;
int uuse[MAXUFUN] , nbrik[MAXUFUN] , brik1[MAXUFUN] ;
void * udata[MAXUFUN] ;
generic_func * ufunc[MAXUFUN] ;
int nuse , uu , newbrik , oldbrik ;
FIMdata fd ;
MRI_IMAGE * tsim ;
float * tsar , * val , ** vbr ;
short * sar ;
int nts , jts ;
MRI_IMARR * imar ;
/* mark which ones to execute */
for( newbrik=nuse=uu=0 ; uu < rlist->num && nuse < MAXUFUN ; uu++ ){
if( (ucode & (1<<uu)) != 0 ){
uuse [nuse] = uu ;
ufunc[nuse] = rlist->funcs[uu] ; /* user_func for this func */
nbrik[nuse] = rlist->flags[uu] ; /* new bricks for this func */
udata[nuse] = rlist->func_data[uu] ; /* user_data for this func */
brik1[nuse] = newbrik ; /* index of 1st brick */
newbrik += nbrik[nuse] ; /* total number of new bricks */
nuse++ ;
}
}
if( nuse == 0 ) goto final_exit ; /* shouldn't happen */
/* do the initialization calls to the user_func functions */
fd.ref_ts = ref_ts ;
fd.ort_ts = ort_ts ;
fd.nvox = nvox ;
fd.ignore = im3d->fimdata->init_ignore ;
fd.polort = polort ;
#ifndef DONT_USE_METER
MCW_set_meter( meter , 0 ) ; meter_pold = 0.0 ;
#endif
for( uu=0 ; uu < nuse ; uu++ )
#if 0
ufunc[uu]( ntime , NULL , udata[uu] , nbrik[uu] , (void *)(&fd) ) ;
#else
AFNI_CALL_fim_function( ufunc[uu] ,
ntime, NULL, udata[uu], nbrik[uu], (&fd) ) ;
#endif
/* loop over voxels,
assemble time series,
call functions to put results in val[],
store float outputs in vbr[][] */
vbr = (float **) malloc(sizeof(float *)*newbrik) ;
for( iv=0 ; iv < newbrik ; iv++ )
vbr[iv] = (float *) malloc(sizeof(float)*nvox) ;
val = (float *) malloc(sizeof(float)*newbrik) ;
for( iv=0 ; iv < nvox ; iv+=MAXTS ){
nts = MIN( MAXTS , nvox-iv ) ;
imar = THD_extract_many_series( nts,indx+iv , dset_time ) ;
for( jts=0 ; jts < nts ; jts++ ){
tsim = IMARR_SUBIMAGE(imar,jts) ; /* data */
tsar = MRI_FLOAT_PTR(tsim) ;
for( uu=0 ; uu < nuse ; uu++ ){
#if 0
ufunc[uu]( ntime , tsar , /* func */
udata[uu] , nbrik[uu] , (void *) val ) ;
#else
AFNI_CALL_fim_function( ufunc[uu] ,
ntime, tsar, udata[uu], nbrik[uu], val ) ;
#endif
for( it=0 ; it < nbrik[uu] ; it++ ) /* storage */
vbr[it+brik1[uu]][iv+jts] = val[it] ;
}
}
DESTROY_IMARR(imar) ; /* garbage disposal */
#ifndef DONT_USE_METER
meter_perc = (int) ( 100.0 * iv / nvox ) ;
if( meter_perc != meter_pold ){
MCW_set_meter( meter , meter_perc ) ;
meter_pold = meter_perc ;
}
#endif
}
free(val) ; /* no longer needed */
#ifndef DONT_USE_METER
MCW_set_meter( meter , 100 ) ;
#endif
/* if necessary, make the new dataset now */
if( new_dset != NULL ){
oldbrik = DSET_NVALS(new_dset) ; /* number of bricks it has now */
} else {
oldbrik = 0 ;
new_dset = EDIT_empty_copy( dset_time ) ;
EDIT_dset_items( new_dset ,
ADN_prefix , new_prefix ,
ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
ADN_type , ISHEAD(dset_time)
? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
ADN_func_type , FUNC_BUCK_TYPE ,
ADN_nvals , newbrik ,
ADN_datum_all , MRI_short ,
ADN_ntt , 0 ,
ADN_none ) ;
}
/* for each output brick:
make short space for it,
scale and stored floats into this space ,
attach it to the output dataset as a new brick */
for( iv=0 ; iv < newbrik ; iv++ ){
tsar = vbr[iv] ; /* float data */
topval = 0.0 ; /* find range of data */
for( it=0 ; it < nvox ; it++ )
if( fabs(tsar[it]) > topval ) topval = fabs(tsar[it]) ;
sar = (short *) calloc(sizeof(short),nxyz) ; /* new brick */
if( topval > 0.0 ){ /* scale to shorts */
topval = MRI_TYPE_maxval[MRI_short] / topval ;
for( it=0 ; it < nvox ; it++ )
sar[indx[it]] = (short)(topval * tsar[it] + 0.499) ;
topval = 1.0/topval ; /* scale factor */
}
free(tsar) ;
if( oldbrik > 0 ){
EDIT_add_brick( new_dset , MRI_short , topval , sar ) ;
} else {
mri_fix_data_pointer( sar , DSET_BRICK(new_dset,iv) ) ;
EDIT_BRICK_FACTOR( new_dset , iv , topval ) ;
}
}
free(vbr) ;
/* do the ending calls to user_func */
for( uu=0 ; uu < nuse ; uu++ )
#if 0
ufunc[uu]( -(brik1[uu]+oldbrik) , NULL ,
udata[uu] , nbrik[uu] , (void *) new_dset ) ;
#else
AFNI_CALL_fim_function( ufunc[uu] ,
-(brik1[uu]+oldbrik) , NULL ,
udata[uu] , nbrik[uu] , new_dset ) ;
#endif
} /* 01 Feb 2000: end of user_func addition to FIMming */
final_exit:
free(vval) ; free(indx) ; /* can finally free these */
/*--- Return new dataset ---*/
#ifndef DONT_USE_METER
MCW_popdown_meter(meter) ;
#endif
RETURN(new_dset) ;
}
syntax highlighted by Code2HTML, v. 0.9.1