00001 
00002 
00003 /*-------------------------------------------------------------------
00004   Actually carry out the FIM computations:
00005     im3d
00006 ---------------------------------------------------------------------*/
00007 
00008 #ifndef FIM_THR
00009 #define FIM_THR  0.0999
00010 #endif
00011 
00012 /** Jan 1998:    code = combinations of the FIM_*_MASK values      **/
00013 /** 30 May 1999: allow polort to be variable (formerly fixed at 1) **/
00014 /** 08 Sep 1999: added PAVE and AVER options                       **/
00015 /** 03 Jan 2000: added PTOP, TOPL, and SIGM options                **/
00016 /** 01 Feb 2000: added ucode, for user written functions           **/
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 ;  /* 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 }