Doxygen Source Code Documentation
afni_fimmer_compute.c File Reference
Go to the source code of this file.
Defines | |
#define | FIM_THR 0.0999 |
#define | METERIZE(ib) |
#define | MAXUFUN 64 |
#define | MAXTS 32 |
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) |
Define Documentation
|
Definition at line 9 of file afni_fimmer_compute.c. Referenced by AFNI_fimmer_compute(). |
|
|
|
|
|
Value: 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) |
Function Documentation
|
01 Feb 2000: added ucode, for user written functions * Definition at line 18 of file afni_fimmer_compute.c. References abs, ADN_brick_fac, ADN_datum_all, ADN_func_type, ADN_malloc_type, ADN_none, ADN_ntt, ADN_nvals, ADN_prefix, ADN_stat_aux, ADN_type, AFNI_3DDATA_VIEW, AFNI_CALL_fim_function, AFNI_GOOD_FUNC_DTYPE, calloc, DATABLOCK_MEM_MALLOC, THD_3dim_dataset::dblk, DESTROY_IMARR, THD_diskptr::dimsizes, THD_datablock::diskptr, THD_slist_find::dset, DSET_ARRAY, DSET_BRICK, DSET_BRICK_BYTES, DSET_BRICK_TYPE, DSET_GRAPHABLE, DSET_NUM_TIMES, DSET_NVALS, DSET_PREFIX, EDIT_add_brick(), EDIT_BRICK_FACTOR, EDIT_BRICK_LABEL, EDIT_BRICK_TO_FICO, EDIT_clip_float(), EDIT_dset_items(), EDIT_empty_copy(), ENTRY, fd, FIM_ALPHA_MASK, FIM_AVER_MASK, FIM_BASE_MASK, FIM_BEST_MASK, FIM_CORR_MASK, FIM_PAVE_MASK, FIM_PERC_MASK, FIM_PTOP_MASK, FIM_SIGM_MASK, FIM_THR, FIM_TOPL_MASK, Three_D_View::fimdata, FIND_PREFIX, MCW_function_list::flags, free, free_PCOR_references(), free_PCOR_voxel_corr(), FUNC_BUCK_TYPE, FUNC_COR_SCALE_SHORT, MCW_function_list::func_data, FUNC_FIM_TYPE, THD_3dim_dataset::func_type, MCW_function_list::funcs, GEN_FUNC_TYPE, generic_func, GLOBAL_library, HEAD_FUNC_TYPE, FIMdata::ignore, IM3D_OPEN, IMARR_SUBIMAGE, AFNI_fimmer_type::init_ignore, ISHEAD, ISVALID_SESSION, MRI_IMAGE::kind, malloc, MAX, MAX_STAT_AUX, MCW_popdown_meter(), MCW_popup_meter(), MCW_set_meter(), MCW_strncpy, METER_TOP_WIDE, MIN, mri_fix_data_pointer(), MRI_FLOAT_PTR, my_getenv(), new_PCOR_references(), new_PCOR_voxel_corr(), MCW_function_list::num, THD_datablock::nvals, FIMdata::nvox, MRI_IMAGE::nx, MRI_IMAGE::ny, FIMdata::ort_ts, PCOR_get_coef(), PCOR_get_pcor(), PCOR_get_perc(), PCOR_get_stdev(), PCOR_update_float(), FIMdata::polort, AFNI_fimmer_type::polort, FIMdata::ref_ts, AFNI_library_type::registered_fim, RETURN, STATUS, strtod(), THD_count_databricks(), THD_delete_3dim_dataset(), THD_dset_in_session(), THD_extract_many_series(), THD_load_datablock(), THD_MAX_PREFIX, AFNI_widget_set::top_shell, Three_D_View::type, update_PCOR_references(), Three_D_View::vwid, and XtRealloc. Referenced by AFNI_fimmer_execute().
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 ; /* 15 Dec 1997 */ 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 ; /* 30 May 1999 */ 00047 00048 float top_perc = 0.0 ; /* 30 Aug 1999 */ 00049 00050 int ibr_pave , ibr_aver ; /* 08 Sep 1999 */ 00051 float * paval , * avval , * pabest , * avbest ; /* 08 Sep 1999 */ 00052 00053 int ibr_ptop , ibr_topl , ibr_sigm ; /* 03 Jan 2000 */ 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 , /* number of updates done yet */ 00063 min_updt = 5 , /* min number needed for display */ 00064 first_updt = 1 ; /* flag to indicate that first update is yet to be displayed */ 00065 00066 ENTRY("AFNI_fimmer_compute") ; 00067 00068 /*--- check for legal inputs ---*/ 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) || /* Jan 1998 & Feb 2000 */ 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 /** 13 Nov 1996: allow for orts **/ 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 /** at this point, ngood_ref = max number of good reference points, 00127 and it1 = index of first point used in first reference **/ 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 /*--- Create a new prefix ---*/ 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) ; /* can't make a new prefix! */ 00150 } 00151 } 00152 00153 if(PRINT_TRACING) 00154 { char str[256] ; 00155 sprintf(str,"new prefix = %s",new_prefix) ; STATUS(str) ; } 00156 00157 /*--- FIM: find values above threshold to fim ---*/ 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 /** find the mean of the first array, 00166 compute the threshold (fthr) from it, 00167 make indx[i] be the 3D index of the i-th voxel above threshold **/ 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 /** allocate space for voxel values **/ 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 /** compute number of output bricks **/ 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 /** 01 Feb 2000: if no normal FIM stuff (code), skip to the ucode stuff **/ 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 ; /* way below */ 00273 } 00274 00275 /** normal case: do the normal recursive FIMming **/ 00276 00277 if(PRINT_TRACING) 00278 { char str[256] ; 00279 sprintf(str,"code of FIM_MASKs = %d",code) ; STATUS(str) ; } 00280 00281 /** allocate extra space for comparing results from multiple ref vectors **/ 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) ; /* 15 Dec 1997 */ 00288 pbest = (float *) malloc(sizeof(float) * nvox) ; /* 16 Jan 1998 */ 00289 bbest = (float *) malloc(sizeof(float) * nvox) ; /* 16 Jan 1998 */ 00290 pval = (float *) malloc(sizeof(float) * nvox) ; /* 16 Jan 1998 */ 00291 bval = (float *) malloc(sizeof(float) * nvox) ; /* 16 Jan 1998 */ 00292 00293 paval = (float *) malloc(sizeof(float) * nvox) ; /* 08 Sep 1999 */ 00294 avval = (float *) malloc(sizeof(float) * nvox) ; /* 08 Sep 1999 */ 00295 pabest = (float *) malloc(sizeof(float) * nvox) ; /* 16 Jan 1998 */ 00296 avbest = (float *) malloc(sizeof(float) * nvox) ; /* 16 Jan 1998 */ 00297 00298 ptval = (float *) malloc(sizeof(float) * nvox) ; /* 03 Jan 2000 */ 00299 tlval = (float *) malloc(sizeof(float) * nvox) ; /* 03 Jan 2000 */ 00300 sgval = (float *) malloc(sizeof(float) * nvox) ; /* 03 Jan 2000 */ 00301 ptbest = (float *) malloc(sizeof(float) * nvox) ; /* 03 Jan 2000 */ 00302 tlbest = (float *) malloc(sizeof(float) * nvox) ; /* 03 Jan 2000 */ 00303 sgbest = (float *) malloc(sizeof(float) * nvox) ; /* 03 Jan 2000 */ 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) ; /* 15 Dec 1997 */ 00312 if( pbest != NULL ) free(pbest) ; /* 16 Jan 1998 */ 00313 if( bbest != NULL ) free(bbest) ; /* 16 Jan 1998 */ 00314 if( pval != NULL ) free(pval) ; /* 16 Jan 1998 */ 00315 if( bval != NULL ) free(bval) ; /* 16 Jan 1998 */ 00316 if( paval != NULL ) free(paval) ; /* 08 Sep 1999 */ 00317 if( avval != NULL ) free(avval) ; /* 08 Sep 1999 */ 00318 if( pabest!= NULL ) free(pabest); /* 08 Sep 1999 */ 00319 if( avbest!= NULL ) free(avbest); /* 08 Sep 1999 */ 00320 if( ptval != NULL ) free(ptval) ; /* 03 Jan 2000 */ 00321 if( tlval != NULL ) free(tlval) ; /* 03 Jan 2000 */ 00322 if( sgval != NULL ) free(sgval) ; /* 03 Jan 2000 */ 00323 if( ptbest!= NULL ) free(ptbest); /* 03 Jan 2000 */ 00324 if( tlbest!= NULL ) free(tlbest); /* 03 Jan 2000 */ 00325 if( sgbest!= NULL ) free(sgbest); /* 03 Jan 2000 */ 00326 RETURN(NULL) ; 00327 } 00328 } else { 00329 aval = rbest = abest = pbest = bbest = pval = bval = NULL ; 00330 paval = avval = pabest = avbest = NULL ; /* 08 Sep 1999 */ 00331 ptval = tlval = ptbest = tlbest = NULL ; /* 03 Jan 2000 */ 00332 sgval = sgbest = NULL ; /* 03 Jan 2000 */ 00333 ibest = NULL ; /* 15 Dec 1997 */ 00334 } 00335 00336 if(PRINT_TRACING) 00337 { char str[256] ; 00338 sprintf(str,"nxyz = %d nvox = %d",nxyz,nvox) ; STATUS(str) ; } 00339 00340 /*--- FIM: initialize recursive updates ---*/ 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) ; /* 15 Dec 1997 */ 00351 if( pbest != NULL ) free(pbest) ; /* 16 Jan 1998 */ 00352 if( bbest != NULL ) free(bbest) ; /* 16 Jan 1998 */ 00353 if( pval != NULL ) free(pval) ; /* 16 Jan 1998 */ 00354 if( bval != NULL ) free(bval) ; /* 16 Jan 1998 */ 00355 if( paval != NULL ) free(paval) ; /* 08 Sep 1999 */ 00356 if( avval != NULL ) free(avval) ; /* 08 Sep 1999 */ 00357 if( pabest!= NULL ) free(pabest); /* 08 Sep 1999 */ 00358 if( avbest!= NULL ) free(avbest); /* 08 Sep 1999 */ 00359 if( ptval != NULL ) free(ptval) ; /* 03 Jan 2000 */ 00360 if( tlval != NULL ) free(tlval) ; /* 03 Jan 2000 */ 00361 if( sgval != NULL ) free(sgval) ; /* 03 Jan 2000 */ 00362 if( ptbest!= NULL ) free(ptbest); /* 03 Jan 2000 */ 00363 if( tlbest!= NULL ) free(tlbest); /* 03 Jan 2000 */ 00364 if( sgbest!= NULL ) free(sgbest); /* 03 Jan 2000 */ 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) ; /* 15 Dec 1997 */ 00386 if( pbest != NULL ) free(pbest) ; /* 16 Jan 1998 */ 00387 if( bbest != NULL ) free(bbest) ; /* 16 Jan 1998 */ 00388 if( pval != NULL ) free(pval) ; /* 16 Jan 1998 */ 00389 if( bval != NULL ) free(bval) ; /* 16 Jan 1998 */ 00390 if( paval != NULL ) free(paval) ; /* 08 Sep 1999 */ 00391 if( avval != NULL ) free(avval) ; /* 08 Sep 1999 */ 00392 if( pabest!= NULL ) free(pabest); /* 08 Sep 1999 */ 00393 if( avbest!= NULL ) free(avbest); /* 08 Sep 1999 */ 00394 if( ptval != NULL ) free(ptval) ; /* 03 Jan 2000 */ 00395 if( tlval != NULL ) free(tlval) ; /* 03 Jan 2000 */ 00396 if( sgval != NULL ) free(sgval) ; /* 03 Jan 2000 */ 00397 if( ptbest!= NULL ) free(ptbest); /* 03 Jan 2000 */ 00398 if( tlbest!= NULL ) free(tlbest); /* 03 Jan 2000 */ 00399 if( sgbest!= NULL ) free(sgbest); /* 03 Jan 2000 */ 00400 fprintf(stderr,"\n*** FIM initialization fails in AFNI_fimmer_compute\n") ; 00401 RETURN(NULL) ; 00402 } 00403 00404 /*--- Make a new dataset to hold the output ---*/ 00405 00406 STATUS("making new dataset") ; 00407 00408 new_dset = EDIT_empty_copy( dset_time ) ; 00409 00410 if( nbrik == 1 && ucode == 0 ){ /* 1 brick out --> a 'fim' dataset */ 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 /* 2 bricks, 2nd corr --> 'fico' */ 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 ){ /* otherwise --> 'fbuc' (bucket) */ 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) ; /* 15 Dec 1997 */ 00462 if( pbest != NULL ) free(pbest) ; /* 16 Jan 1998 */ 00463 if( bbest != NULL ) free(bbest) ; /* 16 Jan 1998 */ 00464 if( pval != NULL ) free(pval) ; /* 16 Jan 1998 */ 00465 if( bval != NULL ) free(bval) ; /* 16 Jan 1998 */ 00466 if( paval != NULL ) free(paval) ; /* 08 Sep 1999 */ 00467 if( avval != NULL ) free(avval) ; /* 08 Sep 1999 */ 00468 if( pabest!= NULL ) free(pabest); /* 08 Sep 1999 */ 00469 if( avbest!= NULL ) free(avbest); /* 08 Sep 1999 */ 00470 if( ptval != NULL ) free(ptval) ; /* 03 Jan 2000 */ 00471 if( tlval != NULL ) free(tlval) ; /* 03 Jan 2000 */ 00472 if( sgval != NULL ) free(sgval) ; /* 03 Jan 2000 */ 00473 if( ptbest!= NULL ) free(ptbest); /* 03 Jan 2000 */ 00474 if( tlbest!= NULL ) free(tlbest); /* 03 Jan 2000 */ 00475 if( sgbest!= NULL ) free(sgbest); /* 03 Jan 2000 */ 00476 RETURN(NULL) ; 00477 } 00478 00479 /* modify labels for each brick */ 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" ) ; /* 08 Sep 1999 */ 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" ) ; /* 03 Jan 2000 */ 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 /*-- 30 Aug 1999: set limits on percent change --*/ 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 /* create bricks */ 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) ; /* 15 Dec 1997 */ 00536 if( pbest != NULL ) free(pbest) ; /* 16 Jan 1998 */ 00537 if( bbest != NULL ) free(bbest) ; /* 16 Jan 1998 */ 00538 if( pval != NULL ) free(pval) ; /* 16 Jan 1998 */ 00539 if( bval != NULL ) free(bval) ; /* 16 Jan 1998 */ 00540 if( paval != NULL ) free(paval) ; /* 08 Sep 1999 */ 00541 if( avval != NULL ) free(avval) ; /* 08 Sep 1999 */ 00542 if( pabest!= NULL ) free(pabest); /* 08 Sep 1999 */ 00543 if( avbest!= NULL ) free(avbest); /* 08 Sep 1999 */ 00544 if( ptval != NULL ) free(ptval) ; /* 03 Jan 2000 */ 00545 if( tlval != NULL ) free(tlval) ; /* 03 Jan 2000 */ 00546 if( sgval != NULL ) free(sgval) ; /* 03 Jan 2000 */ 00547 if( ptbest!= NULL ) free(ptbest); /* 03 Jan 2000 */ 00548 if( tlbest!= NULL ) free(tlbest); /* 03 Jan 2000 */ 00549 if( sgbest!= NULL ) free(sgbest); /* 03 Jan 2000 */ 00550 RETURN(NULL) ; 00551 } 00552 00553 /*---------------------------------*/ 00554 /*--- FIM: do recursive updates ---*/ 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++ ){ /* loop over time */ 00564 00565 nnow = 0 ; /* number of updates done at this time point */ 00566 00567 for( ivec=0 ; ivec < ny_ref ; ivec++ ){ /* loop over ref vects */ 00568 00569 tsar = MRI_FLOAT_PTR(ref_ts) + (ivec*nx_ref) ; /* ptr to vect */ 00570 if( tsar[it] >= WAY_BIG ) continue ; /* skip this */ 00571 00572 ref_vec[0] = 1.0 ; /* we always supply ort for constant */ 00573 for( ip=1 ; ip <= polort ; ip++ ) /* 30 May 1999: */ 00574 ref_vec[ip] = ref_vec[ip-1] * ((float)it) ; /* and polynomials */ 00575 00576 if( internal_ort ){ /* no external orts */ 00577 ref_vec[ip] = tsar[it] ; /* ref value */ 00578 } else { 00579 for( iv=0 ; iv < ny_ort ; iv++ ) /* external */ 00580 ref_vec[iv+ip] = ortar[it + iv*nx_ort] ; /* orts */ 00581 00582 ref_vec[ny_ort+ip] = tsar[it] ; /* ref value */ 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 /* process the ort+ref update */ 00591 00592 update_PCOR_references( ref_vec , pc_ref[ivec] ) ; 00593 00594 /* first time thru: load data from dataset */ 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 /* process the data update */ 00619 00620 PCOR_update_float( vval , pc_ref[ivec] , pc_vc[ivec] ) ; 00621 nnow++ ; /* one more update at this time point */ 00622 } /* end of loop over ref vects */ 00623 00624 if( nnow > 0 ) nupdt++ ; /* number of time points that had updates */ 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 } /* end of loop over time */ 00635 00636 /*-------------------------------------------*/ 00637 /*--- Load final results into the dataset ---*/ 00638 00639 /*--- set the statistical parameters ---*/ 00640 00641 stataux[0] = nupdt ; /* number of points used */ 00642 stataux[1] = (ny_ref==1) ? 1 : 2 ; /* number of references */ 00643 stataux[2] = fim_nref - 1 ; /* number of orts */ 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) /*nada*/ 00661 #endif 00662 00663 /*** Compute brick arrays for new dataset ***/ 00664 /* [load scale factors into stataux, too] */ 00665 00666 if( ny_ref == 1 ){ 00667 00668 /*** Just 1 ref vector --> load values directly into dataset ***/ 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 ){ /* 08 Sep 1999 */ 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 ){ /* 03 Jan 2000 */ 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 ){ /* 08 Sep 1999 */ 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 ){ /* 03 Jan 2000 */ 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 ){ /* 03 Jan 2000 */ 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 /*** Multiple references --> find best correlation at each voxel ***/ 00904 00905 /*--- get first ref results into abest and rbest (best so far) ---*/ 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 ; /* 15 Dec 1997 */ 00917 00918 /*--- for each succeeding ref vector, 00919 get results into aval and vval, 00920 if |vval| > |rbest|, then use that result instead ---*/ 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) ; /* 15 Dec 1997 */ 00940 pbest[iv] = pval[iv] ; /* Jan 1998 */ 00941 bbest[iv] = bval[iv] ; 00942 pabest[iv]= paval[iv] ; /* 08 Sep 1999 */ 00943 avbest[iv]= avval[iv] ; 00944 ptbest[iv]= ptval[iv] ; /* 03 Jan 1999 */ 00945 tlbest[iv]= tlval[iv] ; 00946 sgbest[iv]= sgval[iv] ; 00947 } 00948 } 00949 } 00950 00951 /*--- at this point, abest and rbest are the best 00952 results, so scale them into the dataset bricks ---*/ 00953 00954 /** fim brick **/ 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 /** threshold brick **/ 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 /** best index brick (15 Dec 1997) */ 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 ; /* no scaling */ 01010 01011 METERIZE(ibr_best) ; 01012 } 01013 01014 /** perc brick */ 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 /** pave brick [08 Sep 1999] */ 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 /** ptop brick [03 Jan 2000] */ 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 /** base brick */ 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 /** aver brick [08 Sep 1999] */ 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 /** topl brick [03 Jan 2000] */ 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 /** sigm brick [03 Jan 2000]**/ 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 } /* end of multiple reference case */ 01210 01211 /*** Set the brick factors for the new dataset, 01212 no matter how it was computed above. ***/ 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 /*--- End of recursive updates; now free temporary workspaces ---*/ 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) ; /* 15 Dec 1997 */ 01233 if( pbest != NULL ) free(pbest) ; /* 16 Jan 1998 */ 01234 if( bbest != NULL ) free(bbest) ; /* 16 Jan 1998 */ 01235 if( pval != NULL ) free(pval) ; /* 16 Jan 1998 */ 01236 if( bval != NULL ) free(bval) ; /* 16 Jan 1998 */ 01237 if( paval != NULL ) free(paval) ; /* 08 Sep 1999 */ 01238 if( avval != NULL ) free(avval) ; /* 08 Sep 1999 */ 01239 if( pabest!= NULL ) free(pabest); /* 08 Sep 1999 */ 01240 if( avbest!= NULL ) free(avbest); /* 08 Sep 1999 */ 01241 if( ptval != NULL ) free(ptval) ; /* 03 Jan 2000 */ 01242 if( tlval != NULL ) free(tlval) ; /* 03 Jan 2000 */ 01243 if( sgval != NULL ) free(sgval) ; /* 03 Jan 2000 */ 01244 if( ptbest!= NULL ) free(ptbest); /* 03 Jan 2000 */ 01245 if( tlbest!= NULL ) free(tlbest); /* 03 Jan 2000 */ 01246 if( sgbest!= NULL ) free(sgbest); /* 03 Jan 2000 */ 01247 01248 /*-----------------------------------------------------*/ 01249 /*--- 01 Feb 2000: execute user specified functions ---*/ 01250 01251 ucode_stuff: 01252 01253 #define MAXUFUN 64 /* should be at least sizeof(int) */ 01254 #define MAXTS 32 /* number of timeseries to process at once */ 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 /* mark which ones to execute */ 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] ; /* user_func for this func */ 01275 nbrik[nuse] = rlist->flags[uu] ; /* new bricks for this func */ 01276 udata[nuse] = rlist->func_data[uu] ; /* user_data for this func */ 01277 brik1[nuse] = newbrik ; /* index of 1st brick */ 01278 newbrik += nbrik[nuse] ; /* total number of new bricks */ 01279 nuse++ ; 01280 } 01281 } 01282 01283 if( nuse == 0 ) goto final_exit ; /* shouldn't happen */ 01284 01285 /* do the initialization calls to the user_func functions */ 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 /* loop over voxels, 01306 assemble time series, 01307 call functions to put results in val[], 01308 store float outputs in vbr[][] */ 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) ; /* data */ 01322 tsar = MRI_FLOAT_PTR(tsim) ; 01323 01324 for( uu=0 ; uu < nuse ; uu++ ){ 01325 #if 0 01326 ufunc[uu]( ntime , tsar , /* func */ 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++ ) /* storage */ 01334 vbr[it+brik1[uu]][iv+jts] = val[it] ; 01335 } 01336 } 01337 01338 DESTROY_IMARR(imar) ; /* garbage disposal */ 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) ; /* no longer needed */ 01349 #ifndef DONT_USE_METER 01350 MCW_set_meter( meter , 100 ) ; 01351 #endif 01352 01353 /* if necessary, make the new dataset now */ 01354 01355 if( new_dset != NULL ){ 01356 oldbrik = DSET_NVALS(new_dset) ; /* number of bricks it has now */ 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 /* for each output brick: 01375 make short space for it, 01376 scale and stored floats into this space , 01377 attach it to the output dataset as a new brick */ 01378 01379 for( iv=0 ; iv < newbrik ; iv++ ){ 01380 tsar = vbr[iv] ; /* float data */ 01381 topval = 0.0 ; /* find range of data */ 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) ; /* new brick */ 01386 01387 if( topval > 0.0 ){ /* scale to shorts */ 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 ; /* scale factor */ 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 /* do the ending calls to user_func */ 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 } /* 01 Feb 2000: end of user_func addition to FIMming */ 01419 01420 final_exit: 01421 free(vval) ; free(indx) ; /* can finally free these */ 01422 01423 /*--- Return new dataset ---*/ 01424 01425 #ifndef DONT_USE_METER 01426 MCW_popdown_meter(meter) ; 01427 #endif 01428 01429 RETURN(new_dset) ; 01430 } |