00066 {
00067    THD_3dim_dataset * new_dset ;
00068    int ifim , it,iv , nvox , ngood_ref , ntime , it1 , dtyp , nxyz;
00069    float * vval , * tsar , * aval , * rbest , * abest ;
00070    int   * indx ;
00071    short * bar ;
00072    void  * ptr ;
00073    float stataux[MAX_STAT_AUX];
00074    float fthr , topval ;
00075    int nx_ref , ny_ref , ivec , nnow ;
00076    PCOR_references ** pc_ref ;
00077    PCOR_voxel_corr ** pc_vc ;
00078    int save_resam ;
00079 
00080    int fim_nref , nx_ort , ny_ort , internal_ort ;    /* 10 Dec 1996 */
00081    static float * ref_vec = NULL ;
00082    static int    nref_vec = -666 ;
00083 
00084    float * ref_ts_min = NULL, 
00085          * ref_ts_max = NULL, 
00086          * baseline   = NULL;      /* 19 May 1997 */
00087 
00088    int i;
00089    
00090    int nupdt      = 0 ,  /* number of updates done yet */
00091        min_updt   = 5 ;  /* min number needed for display */
00092 
00093 
00094    /*--- check for legal inputs ---*/      /* 14 Jan 1998 */
00095 
00096    if (!DSET_GRAPHABLE(dset_time)) 
00097      {
00098        fprintf (stderr, "Error:  Invalid 3d+time input data file \n");
00099        RETURN (NULL);
00100      }
00101    
00102    if (ref_ts == NULL)
00103      {
00104        fprintf (stderr, "Error:  No ideal time series \n");
00105        RETURN (NULL);
00106      }
00107 
00108    for (i = 0;  i < ref_ts->num;  i++)
00109      if (ref_ts->tsarr[i]->len < DSET_NUM_TIMES(dset_time))
00110        { 
00111          char str[256] ;
00112          sprintf (str,
00113            "Error:  ideal time series is too short: ntime=%d num_ts=%d \n",
00114                   DSET_NUM_TIMES(dset_time), 
00115                   ref_ts->tsarr[i]->len);
00116          fprintf (stderr, str) ;    
00117          RETURN (NULL) ;
00118        }
00119 
00120 
00121    /** 10 Dec 1996: allow for orts **/
00122 
00123    if( ort_ts->num > 0 )      /** 05 Sept 1997 **/
00124      {
00125        internal_ort = 0;
00126        ny_ort = ort_ts->num;
00127        for (i = 0;  i < ny_ort;  i++)
00128          {
00129            nx_ort = ort_ts->tsarr[i]->len ;
00130            if (nx_ort < DSET_NUM_TIMES(dset_time))   /* 14 Jan 1998 */
00131              { 
00132                char str[256] ;
00133                sprintf (str,
00134                  "Error:  ort time series is too short: ntime=%d ort_ts=%d \n",
00135                         DSET_NUM_TIMES(dset_time), 
00136                         ort_ts->tsarr[i]->len);
00137                fprintf (stderr, str) ;    
00138                RETURN (NULL) ;
00139              }     
00140          }
00141      } 
00142    else 
00143      {
00144        internal_ort = 1 ;
00145      }
00146    fim_nref = (internal_ort) ? 3 : (ny_ort+3) ;
00147 
00148    if( nref_vec < fim_nref )
00149      {
00150        ref_vec = (float *) malloc (sizeof(float)*fim_nref) ;
00151        nref_vec = fim_nref;
00152      }
00153 
00154 
00155    /* arrays to store maximum change in the ideal time series */
00156    if (max_percent > 0.0)    /* 19 May 1997 */
00157      {
00158        ref_ts_max = (float *) malloc (sizeof(float) * (ref_ts->num));
00159        ref_ts_min = (float *) malloc (sizeof(float) * (ref_ts->num));
00160      }
00161 
00162 
00163    nx_ref    = ref_ts->tsarr[0]->len;
00164    ny_ref    = ref_ts->num;
00165    ntime     = DSET_NUM_TIMES(dset_time) ;
00166    ngood_ref = 0 ;
00167    it1      = -1 ;
00168    for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00169       tsar = ref_ts->tsarr[ivec]->ts;
00170       ifim = 0 ;
00171 
00172       if (max_percent > 0.0)       /* 19 May 1997 */
00173         {
00174           ref_ts_min[ivec] = (float) SO_BIG;              
00175           ref_ts_max[ivec] = - (float) SO_BIG;
00176         }
00177 
00178       for( it=itbot ; it < ntime ; it++ )
00179         {
00180          if( tsar[it] < SO_BIG )
00181            { 
00182              ifim++ ; 
00183              if( it1 < 0 ) it1 = it ;
00184 
00185              if (max_percent > 0.0)      /* 19 May 1997 */
00186                {
00187                  if (tsar[it] > ref_ts_max[ivec])  ref_ts_max[ivec] = tsar[it];
00188                  if (tsar[it] < ref_ts_min[ivec])  ref_ts_min[ivec] = tsar[it];
00189                }
00190            }
00191         }
00192 
00193       if( ifim < min_updt ){
00194          STATUS("ref_ts has too few good entries!") ;
00195          RETURN(NULL) ;
00196       }
00197 
00198       ngood_ref = MAX( ifim , ngood_ref ) ;
00199    }
00200 
00201    /** at this point, ngood_ref = max number of good reference points,
00202        and                  it1 = index of first point used in first reference **/
00203    
00204    dtyp = DSET_BRICK_TYPE(dset_time,it1) ;
00205    if( ! AFNI_GOOD_FUNC_DTYPE(dtyp) ){
00206       STATUS("illegal input data type!") ;
00207       RETURN(NULL) ;
00208    }
00209 
00210 
00211 #ifdef AFNI_DEBUG
00212 { char str[256] ;
00213   sprintf(str,"new prefix = %s",new_prefix) ; STATUS(str) ; }
00214 #endif
00215 
00216    /*--- FIM: find values above threshold to fim ---*/
00217 
00218    THD_load_datablock( dset_time->dblk ) ;
00219 
00220    nxyz =  dset_time->dblk->diskptr->dimsizes[0]
00221          * dset_time->dblk->diskptr->dimsizes[1]
00222          * dset_time->dblk->diskptr->dimsizes[2] ;
00223 
00224    /** find the mean of the first array,
00225        compute the threshold (fthr) from it,
00226        make indx[i] be the 3D index of the i-th voxel above threshold **/
00227 
00228    switch( dtyp ){
00229 
00230       case MRI_short:{
00231          short * dar = (short *) DSET_ARRAY(dset_time,it1) ;
00232          for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += abs(dar[iv]) ;
00233          fthr = FIM_THR * fthr / nxyz ;
00234          for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
00235             if( abs(dar[iv]) > fthr ) nvox++ ;
00236          indx = (int *) malloc( sizeof(int) * nvox ) ;
00237          if( indx == NULL ){
00238             fprintf(stderr,"\n*** indx malloc failure in fim3d_fimmer_compute\n") ;
00239             RETURN(NULL) ;
00240          }
00241          for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
00242             if( abs(dar[iv]) > fthr ) indx[nvox++] = iv ;
00243       }
00244       break ;
00245 
00246       case MRI_float:{
00247          float * dar = (float *) DSET_ARRAY(dset_time,it1) ;
00248          for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += fabs(dar[iv]) ;
00249          fthr = FIM_THR * fthr / nxyz ;
00250          for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
00251             if( fabs(dar[iv]) > fthr ) nvox++ ;
00252          indx = (int *) malloc( sizeof(int) * nvox ) ;
00253          if( indx == NULL ){
00254             fprintf(stderr,"\n*** indx malloc failure in fim3d_fimmer_compute\n") ;
00255             RETURN(NULL) ;
00256          }
00257          for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
00258             if( fabs(dar[iv]) > fthr ) indx[nvox++] = iv ;
00259       }
00260       break ;
00261 
00262       case MRI_byte:{
00263          byte * dar = (byte *) DSET_ARRAY(dset_time,it1) ;
00264          for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += dar[iv] ;
00265          fthr = FIM_THR * fthr / nxyz ;
00266          for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
00267             if( dar[iv] > fthr ) nvox++ ;
00268          indx = (int *) malloc( sizeof(int) * nvox ) ;
00269          if( indx == NULL ){
00270             fprintf(stderr,"\n*** indx malloc failure in fim3d_fimmer_compute\n") ;
00271             RETURN(NULL) ;
00272          }
00273          for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
00274             if( dar[iv] > fthr ) indx[nvox++] = iv ;
00275       }
00276       break ;
00277    }
00278 
00279    /** allocate space for voxel values **/
00280 
00281    vval = (float *) malloc( sizeof(float) * nvox) ;
00282    if( vval == NULL ){
00283       fprintf(stderr,"\n*** vval malloc failure in fim3d_fimmer_compute\n") ;
00284       free(indx) ; RETURN(NULL) ;
00285    }
00286 
00287   
00288    /*----- allocate space for baseline values -----*/
00289    if (max_percent > 0.0)    /* 19 May 1997 */
00290      {
00291        baseline = (float *) malloc (sizeof(float) * nvox);
00292        if (baseline == NULL)
00293          {
00294            fprintf(stderr,
00295                    "\n*** baseline malloc failure in fim3d_fimmer_compute\n") ;
00296            free(indx) ; free(vval); RETURN(NULL) ;
00297          }
00298        else  /* initialize baseline values to zero */
00299          for (iv = 0;  iv < nvox;  iv++)
00300            baseline[iv] = 0.0;
00301      } 
00302 
00303 
00304    /** allocate extra space for comparing results from multiple ref vectors **/
00305 
00306    if( ny_ref > 1 ){
00307       aval  = (float *) malloc( sizeof(float) * nvox) ;
00308       rbest = (float *) malloc( sizeof(float) * nvox) ;
00309       abest = (float *) malloc( sizeof(float) * nvox) ;
00310       if( aval==NULL || rbest==NULL || abest==NULL ){
00311          fprintf(stderr,"\n*** abest malloc failure in fim3d_fimmer_compute\n") ;
00312          free(vval) ; free(indx) ;
00313          if( aval  != NULL ) free(aval) ;
00314          if( rbest != NULL ) free(rbest) ;
00315          if( abest != NULL ) free(abest) ;
00316          RETURN(NULL) ;
00317       }
00318    } else {
00319       aval = rbest = abest = NULL ;
00320    }
00321 
00322 #ifdef AFNI_DEBUG
00323 { char str[256] ;
00324   sprintf(str,"nxyz = %d  nvox = %d",nxyz,nvox) ; STATUS(str) ; }
00325 #endif
00326 
00327    /*--- FIM: initialize recursive updates ---*/
00328 
00329    pc_ref = (PCOR_references **) malloc( sizeof(PCOR_references *) * ny_ref ) ;
00330    pc_vc  = (PCOR_voxel_corr **) malloc( sizeof(PCOR_voxel_corr *) * ny_ref ) ;
00331 
00332    if( pc_ref == NULL || pc_vc == NULL ){
00333       free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
00334       if( aval  != NULL ) free(aval) ;
00335       if( rbest != NULL ) free(rbest) ;
00336       if( abest != NULL ) free(abest) ;
00337       fprintf(stderr,"\n*** FIM initialization fails in fim3d_fimmer_compute\n") ;
00338       RETURN(NULL) ;
00339    }
00340 
00341    ifim = 0 ;
00342    for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00343       pc_ref[ivec] = new_PCOR_references( fim_nref ) ;
00344       pc_vc[ivec]  = new_PCOR_voxel_corr( nvox , fim_nref ) ;
00345       if( pc_ref[ivec] == NULL || pc_vc[ivec] == NULL ) ifim++ ;
00346    }
00347 
00348    if( ifim > 0 ){
00349       for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00350          free_PCOR_references(pc_ref[ivec]) ;
00351          free_PCOR_voxel_corr(pc_vc[ivec]) ;
00352       }
00353       free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
00354       if( aval  != NULL ) free(aval) ;
00355       if( rbest != NULL ) free(rbest) ;
00356       if( abest != NULL ) free(abest) ;
00357       fprintf(stderr,"\n*** FIM initialization fails in fim3d_fimmer_compute\n") ;
00358       RETURN(NULL) ;
00359    }
00360 
00361    /*--- Make a new dataset to hold the output ---*/
00362 
00363    new_dset = EDIT_empty_copy( dset_time ) ;
00364 
00365    it = EDIT_dset_items( new_dset ,
00366                             ADN_prefix      , new_prefix ,
00367                             ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00368                             ADN_type        , ISHEAD(dset_time)
00369                                               ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
00370                             ADN_func_type   , FUNC_COR_TYPE ,
00371                             ADN_nvals       , FUNC_nvals[FUNC_COR_TYPE] ,
00372                             ADN_datum_all   , MRI_short ,
00373                             ADN_ntt         , 0 ,
00374                          ADN_none ) ;
00375 
00376    if( it > 0 ){
00377       fprintf(stderr,
00378               "\n*** EDIT_dset_items error %d in fim3d_fimmer_compute\n",it) ;
00379       THD_delete_3dim_dataset( new_dset , False ) ;
00380       for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00381          free_PCOR_references(pc_ref[ivec]) ;
00382          free_PCOR_voxel_corr(pc_vc[ivec]) ;
00383       }
00384       free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
00385       if( aval  != NULL ) free(aval) ;
00386       if( rbest != NULL ) free(rbest) ;
00387       if( abest != NULL ) free(abest) ;
00388       RETURN(NULL) ;
00389    }
00390 
00391    for( iv=0 ; iv < new_dset->dblk->nvals ; iv++ ){
00392       ptr = malloc( DSET_BRICK_BYTES(new_dset,iv) ) ;
00393       mri_fix_data_pointer( ptr ,  DSET_BRICK(new_dset,iv) ) ;
00394    }
00395 
00396    if( THD_count_databricks(new_dset->dblk) < new_dset->dblk->nvals ){
00397       fprintf(stderr,
00398               "\n*** failure to malloc new bricks in fim3d_fimmer_compute\n") ;
00399       THD_delete_3dim_dataset( new_dset , False ) ;
00400       for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00401          free_PCOR_references(pc_ref[ivec]) ;
00402          free_PCOR_voxel_corr(pc_vc[ivec]) ;
00403       }
00404       free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
00405       if( aval  != NULL ) free(aval) ;
00406       if( rbest != NULL ) free(rbest) ;
00407       if( abest != NULL ) free(abest) ;
00408       RETURN(NULL) ;
00409    }
00410 
00411 
00412    /*--- FIM: do recursive updates ---*/
00413 
00414    for( it=itbot ; it < ntime ; it++ ){
00415 
00416       nnow = 0 ;
00417       for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00418          tsar = ref_ts->tsarr[ivec]->ts ;
00419          if( tsar[it] >= SO_BIG ) continue ;  /* skip this */
00420 
00421          ref_vec[0] = 1.0 ;         /* we always supply orts */
00422          ref_vec[1] = (float) it ;  /* for mean and linear trend */
00423 
00424          if (internal_ort)          /* 10 Dec 1996 */
00425            {
00426              ref_vec[2] = tsar[it] ;
00427            } 
00428          else 
00429            {
00430              for( iv=0 ; iv < ny_ort ; iv++ )
00431                ref_vec[iv+2] = ort_ts->tsarr[iv]->ts[it];
00432              ref_vec[ny_ort+2] = tsar[it] ;
00433            }
00434 
00435 
00436 #ifdef AFNI_DEBUG
00437 { char str[256] ;
00438   sprintf(str,"time index=%d  ideal[%d]=%f" , it,ivec,tsar[it] ) ;
00439   if (ivec == 0) STATUS(str) ; }
00440 #endif
00441 
00442 
00443          update_PCOR_references( ref_vec , pc_ref[ivec] ) ;
00444 
00445          switch( dtyp ){
00446             case MRI_short:{
00447                short * dar = (short *) DSET_ARRAY(dset_time,it) ;
00448                for( iv=0 ; iv < nvox ; iv++ ) vval[iv] = (float) dar[indx[iv]] ;
00449             }
00450             break ;
00451 
00452             case MRI_float:{
00453                float * dar = (float *) DSET_ARRAY(dset_time,it) ;
00454                for( iv=0 ; iv < nvox ; iv++ ) vval[iv] = (float) dar[indx[iv]] ;
00455             }
00456             break ;
00457 
00458             case MRI_byte:{
00459                byte * dar = (byte *) DSET_ARRAY(dset_time,it) ;
00460                for( iv=0 ; iv < nvox ; iv++ ) vval[iv] = (float) dar[indx[iv]] ;
00461             }
00462             break ;
00463          }
00464 
00465          PCOR_update_float( vval , pc_ref[ivec] , pc_vc[ivec] ) ;
00466          nnow++ ;
00467 
00468          /*----- update baseline value calculation -----*/
00469          if (max_percent > 0.0)    /* 19 May 1997 */
00470            if (ivec == 0)
00471              for (iv = 0;  iv < nvox;  iv++)
00472                baseline[iv] += vval[iv] / ngood_ref;
00473  
00474       }
00475       if( nnow > 0 ) nupdt++ ;
00476 
00477 
00478       /*--- Load results into the dataset and redisplay it ---*/
00479 
00480       if( nupdt == ngood_ref ) 
00481       {
00482          /*--- set the statistical parameters ---*/
00483 
00484          stataux[0] = nupdt ;               /* number of points used */
00485          stataux[1] = (ny_ref==1) ? 1 : 2 ; /* number of references  */
00486          stataux[2] = fim_nref - 1 ;     /* number of orts */  /* 12 Dec 96 */
00487          for( iv=3 ; iv < MAX_STAT_AUX ; iv++ ) stataux[iv] = 0.0 ;
00488 
00489 STATUS("setting statistical parameters") ;
00490 
00491          (void) EDIT_dset_items( new_dset ,
00492                                     ADN_stat_aux , stataux ,
00493                                  ADN_none ) ;
00494 
00495          /*** Compute brick arrays for new dataset ***/
00496 
00497          if( ny_ref == 1 ){
00498 
00499          /*** Just 1 ref vector --> load values directly into dataset ***/
00500 
00501             /*--- get alpha (coef) into vval,
00502                   find max value, scale into brick array ---*/
00503 
00504 STATUS("getting 1 ref alpha") ;
00505 
00506             PCOR_get_coef( pc_ref[0] , pc_vc[0] , vval ) ;
00507 
00508             /*--- replace alpha with percentage change, if so requested ---*/
00509             if (max_percent > 0.0)    /* 19 May 1997 */
00510               {
00511                 for (iv = 0;  iv < nvox;  iv++)
00512                   {
00513                     vval[iv] *= 100.0 * (ref_ts_max[0] - ref_ts_min[0]);
00514                     if (fabs(vval[iv]) < max_percent * fabs(baseline[iv]))
00515                       vval[iv] = fabs( vval[iv] / baseline[iv] );
00516                     else
00517                       vval[iv] = max_percent;
00518                   }
00519                 topval = max_percent;
00520               }
00521             else 
00522               {
00523                 topval = 0.0 ;
00524                 for( iv=0 ; iv < nvox ; iv++ )
00525                   if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
00526               }
00527 
00528             bar = DSET_ARRAY( new_dset , FUNC_ival_fim[FUNC_COR_TYPE] ) ;
00529 
00530 #ifdef DONT_USE_MEMCPY
00531             for( iv=0 ; iv < nxyz ; iv++ ) bar[iv] = 0 ;
00532 #else
00533             memset( bar , 0 , sizeof(short)*nxyz ) ;
00534 #endif
00535 
00536             if( topval > 0.0 ){
00537                topval = MRI_TYPE_maxval[MRI_short] / topval ;
00538                for( iv=0 ; iv < nvox ; iv++ )
00539                   bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
00540 
00541                stataux[0] = 1.0/topval ;
00542             } else {
00543                stataux[0] = 0.0 ;
00544             }
00545 
00546             /*--- get correlation coefficient (pcor) into vval,
00547                   scale into brick array (with fixed scaling factor) ---*/
00548 
00549 STATUS("getting 1 ref pcor") ;
00550 
00551             PCOR_get_pcor( pc_ref[0] , pc_vc[0] , vval ) ;
00552 
00553             bar = DSET_ARRAY( new_dset , FUNC_ival_thr[FUNC_COR_TYPE] ) ;
00554 
00555 #ifdef DONT_USE_MEMCPY
00556             for( iv=0 ; iv < nxyz ; iv++ ) bar[iv] = 0 ;
00557 #else
00558             memset( bar , 0 , sizeof(short)*nxyz ) ;
00559 #endif
00560 
00561             for( iv=0 ; iv < nvox ; iv++ )
00562                bar[indx[iv]] = (short)(FUNC_COR_SCALE_SHORT * vval[iv] + 0.499) ;
00563 
00564             stataux[1] = 1.0 / FUNC_COR_SCALE_SHORT ;
00565 
00566          } else {
00567 
00568          /*** Multiple references --> find best correlation at each voxel ***/
00569 
00570             /*--- get first ref results into abest and rbest (best so far) ---*/
00571 
00572             PCOR_get_coef( pc_ref[0] , pc_vc[0] , abest ) ;
00573 
00574             /*--- modify alpha for percentage change calculation ---*/
00575             if (max_percent > 0.0)    /* 19 May 1997 */
00576               for (iv = 0;  iv < nvox;  iv++)
00577                 abest[iv] *= 100.0 * (ref_ts_max[0] - ref_ts_min[0]);          
00578               
00579             PCOR_get_pcor( pc_ref[0] , pc_vc[0] , rbest ) ;
00580 
00581             /*--- for each succeeding ref vector,
00582                   get results into aval and vval,
00583                   if |vval| > |rbest|, then use that result instead ---*/
00584 
00585             for( ivec=1 ; ivec < ny_ref ; ivec++ ){
00586 
00587                PCOR_get_coef( pc_ref[ivec] , pc_vc[ivec] , aval ) ;
00588 
00589                PCOR_get_pcor( pc_ref[ivec] , pc_vc[ivec] , vval ) ;
00590 
00591                for( iv=0 ; iv < nvox ; iv++ ){
00592                   if( fabs(vval[iv]) > fabs(rbest[iv]) ){
00593                      rbest[iv] = vval[iv] ;
00594                      abest[iv] = aval[iv] ;
00595 
00596                      /*--- modify alpha for percentage change calculation ---*/
00597                      if (max_percent > 0.0)    /* 19 May 1997 */
00598                        abest[iv] *= 100.0 *
00599                          (ref_ts_max[ivec] - ref_ts_min[ivec]);
00600 
00601                   }
00602                }
00603 
00604             }
00605 
00606             /*--- at this point, abest and rbest are the best
00607                   results, so scale them into the dataset bricks ---*/
00608 
00609             /*--- finish percentage change calculation, if so requested ---*/
00610             if (max_percent > 0.0)    /* 19 May 1997 */
00611               {
00612                 for (iv = 0;  iv < nvox;  iv++)
00613                   {
00614                     if (fabs(abest[iv]) < max_percent * fabs(baseline[iv]))
00615                       abest[iv] = fabs( abest[iv] / baseline[iv] );
00616                     else
00617                       abest[iv] = max_percent;
00618                   }
00619                 topval = max_percent;
00620               }
00621             else
00622               {
00623                 topval = 0.0 ;
00624                 for( iv=0 ; iv < nvox ; iv++ )
00625                   if( fabs(abest[iv]) > topval ) topval = fabs(abest[iv]) ;
00626               }
00627 
00628             bar = DSET_ARRAY( new_dset , FUNC_ival_fim[FUNC_COR_TYPE] ) ;
00629 
00630 #ifdef DONT_USE_MEMCPY
00631             for( iv=0 ; iv < nxyz ; iv++ ) bar[iv] = 0 ;
00632 #else
00633             memset( bar , 0 , sizeof(short)*nxyz ) ;
00634 #endif
00635 
00636             if( topval > 0.0 ){
00637                topval = MRI_TYPE_maxval[MRI_short] / topval ;
00638                for( iv=0 ; iv < nvox ; iv++ )
00639                   bar[indx[iv]] = (short)(topval * abest[iv] + 0.499) ;
00640 
00641                stataux[0] = 1.0/topval ;
00642             } else {
00643                stataux[0] = 0.0 ;
00644             }
00645 
00646             bar = DSET_ARRAY( new_dset , FUNC_ival_thr[FUNC_COR_TYPE] ) ;
00647 
00648 #ifdef DONT_USE_MEMCPY
00649             for( iv=0 ; iv < nxyz ; iv++ ) bar[iv] = 0 ;
00650 #else
00651             memset( bar , 0 , sizeof(short)*nxyz ) ;
00652 #endif
00653 
00654             for( iv=0 ; iv < nvox ; iv++ )
00655                bar[indx[iv]] = (short)(FUNC_COR_SCALE_SHORT * rbest[iv] + 0.499) ;
00656 
00657             stataux[1] = 1.0 / FUNC_COR_SCALE_SHORT ;
00658 
00659          }
00660 
00661 STATUS("setting brick_fac") ;
00662 
00663          (void) EDIT_dset_items( new_dset ,
00664                                     ADN_brick_fac , stataux ,
00665                                  ADN_none ) ;
00666 
00667       }
00668    }
00669 
00670  
00671    /*--- End of recursive updates; now free temporary workspaces ---*/
00672 
00673    for( ivec=0 ; ivec < ny_ref ; ivec++ ){
00674       free_PCOR_references(pc_ref[ivec]) ;
00675       free_PCOR_voxel_corr(pc_vc[ivec]) ;
00676    }
00677    free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
00678    if( aval  != NULL ) free(aval) ;
00679    if( rbest != NULL ) free(rbest) ;
00680    if( abest != NULL ) free(abest) ;
00681 
00682    if (ref_ts_min != NULL)  free (ref_ts_min);    /* 19 May 1997 */
00683    if (ref_ts_max != NULL)  free (ref_ts_max);
00684    if (baseline != NULL)    free (baseline);
00685 
00686 
00687    /* --- load the statistics --- */
00688    THD_load_statistics (new_dset);
00689    
00690    /*--- Return new dataset ---*/
00691 
00692    RETURN(new_dset) ;
00693 }

00728 {
00729   
00730    int nopt ;
00731    float ftemp ;
00732    MRI_IMAGE *im ;
00733    time_series *ideal;
00734    time_series *ort;
00735    char err_note[128];
00736    Boolean ok;
00737    char input_filename[THD_MAX_NAME];
00738 
00739 
00740    /* --- give help if needed --- */
00741    if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) Syntax(NULL) ;
00742   
00743    /*----- add to program log -----*/
00744    AFNI_logger (PROGRAM_NAME,argc,argv); 
00745 
00746    /* --- setup opt with defaults --- */
00747    strcpy (opt->prefix_name, "");   /* no prefix name yet */
00748    opt->first_im = 0 ;              /* first image to actually use */
00749    opt->max_percent = 0.0;
00750 
00751    /* --- initialize array of time series data --- */
00752    INIT_TSARR(opt->idts) ;
00753 
00754    /* --- initialize array of ort time series data --- */  /* 10 Dec 1996 */
00755    INIT_TSARR (opt->ortts);
00756 
00757    /* --- set defaults in local variables --- */
00758    ideal    = NULL ;                /* time_series of ideal response vector */
00759    strcpy (input_filename, "");     /* no input file name yet */
00760 
00761 
00762    /* --- read command line arguments and process them:
00763       coding technique is to examine each argv, and if it matches
00764       something we expect (e.g., -ideal), process it, then skip to
00765       the next loop through ('continue' statements in each strncmp if).
00766       Note that the 'while' will increment the argument index (nopt)  --- */
00767 
00768    nopt = 1 ;
00769    
00770    do{ 
00771 
00772 #ifdef DEBUG
00773 #  define DBOPT(str) fprintf(stderr,"at option %s: %s\n",argv[nopt],str)
00774 #else
00775 #  define DBOPT(str) /* nothing */
00776 #endif
00777 
00778       /* ---  -im1 #  ==>  index (1...) of first image to actually use  --- */
00779       if( strncmp(argv[nopt],"-im1",4) == 0 )
00780       {
00781          DBOPT("-im1") ;
00782          if( ++nopt >= argc ) Syntax("-im1 needs an argument") ;
00783          ftemp = strtod( argv[nopt] , NULL ) ;
00784          if( ftemp < 1.0 )
00785          {
00786             sprintf( err_note , "illegal -im1 value: %f" , ftemp ) ;
00787             Syntax(err_note) ;
00788          }
00789          opt->first_im = (int)(ftemp+0.499) - 1;
00790          continue ;
00791       }
00792 
00793 
00794       /* --- -ideal file  ==>  ideal response vector time series --- */
00795       if( strncmp(argv[nopt],"-ideal",5) == 0 )
00796       {
00797          DBOPT("-ideal") ;
00798          if( ++nopt >= argc ) Syntax("-ideal needs a filename") ;
00799 
00800          /** July 1995: read a bunch of ideals using [ a b c ... ] format **/
00801 
00802          if( ! IS_OPENER(argv[nopt]) )
00803          {  /* --- one file --- */
00804             ideal = RWC_read_time_series( argv[nopt] ) ;
00805             if( ideal == NULL )  Syntax ("Unable to read ideal time series.");
00806             ADDTO_TSARR( opt->idts , ideal ) ;
00807          } 
00808          else 
00809          {  /* --- many files --- */
00810             for( nopt++ ; !IS_CLOSER(argv[nopt]) && nopt<argc ; nopt++ )
00811             {
00812                ideal = RWC_read_time_series( argv[nopt] ) ;
00813                if( ideal == NULL )  Syntax ("Unable to read ideal time series.");
00814                ADDTO_TSARR( opt->idts , ideal ) ;
00815             }
00816             if( nopt == argc ) Syntax("-ideal never finishes!") ;
00817          }
00818          continue ;
00819       }
00820 
00821 
00822       /*-----  -percent p   Calculate percentage change from baseline -----*/
00823       if( strncmp(argv[nopt],"-percent",8) == 0 )
00824       {
00825          DBOPT("-percent") ;
00826          if( ++nopt >= argc ) Syntax("-percent needs an argument") ;
00827          ftemp = strtod( argv[nopt] , NULL ) ;
00828          if( ftemp <= 0.0 )
00829          {
00830             sprintf( err_note , "illegal -percent value: %f" , ftemp ) ;
00831             Syntax(err_note) ;
00832          }
00833          opt->max_percent = ftemp;
00834          continue ;
00835       }
00836 
00837 
00838       /* --- -ort file  ==>  ort vector time series --- */   /* 10 Dec 1996 */
00839       if( strncmp(argv[nopt],"-ort",4) == 0 )
00840       {
00841          DBOPT("-ort") ;
00842          if( ++nopt >= argc ) Syntax("-ort needs a filename") ;
00843 
00844          /**  read a bunch of orts using [ a b c ... ] format **/
00845 
00846          if( ! IS_OPENER(argv[nopt]) )
00847          {  /* --- one file --- */
00848             ort = RWC_read_time_series( argv[nopt] ) ;
00849             if( ort == NULL )  Syntax ("Unable to read ort time series.");
00850             ADDTO_TSARR( opt->ortts , ort ) ;
00851          } 
00852          else 
00853          {  /* --- many files --- */
00854             for( nopt++ ; !IS_CLOSER(argv[nopt]) && nopt<argc ; nopt++ )
00855             {
00856                ort = RWC_read_time_series( argv[nopt] ) ;
00857                if( ort == NULL )  Syntax ("Unable to read ort time series.");
00858                ADDTO_TSARR( opt->ortts , ort ) ;
00859             }
00860             if( nopt == argc ) Syntax("-ort never finishes!") ;
00861          }
00862          continue ;
00863       }
00864 
00865 
00866       /* ---  -prefix name  ==>  prefix name of output file --- */
00867       if( strncmp(argv[nopt], "-prefix", 5) == 0 ){
00868          DBOPT("-prefix") ;
00869          if( ++nopt >= argc ) Syntax("-prefix needs a name") ;
00870          if( strcmp(opt->prefix_name, "") )  
00871             Syntax("Cannot have two prefix output names!" ) ;
00872          strcpy (opt->prefix_name, argv[nopt]) ;
00873          DBOPT("stored as prefix") ; 
00874          continue ;
00875       }
00876 
00877       /* ---  -input fname  ==>  file name of input 3d+time data --- */
00878 
00879       if( strncmp(argv[nopt], "-input", 5) == 0 )
00880       {
00881          DBOPT("-input") ;
00882          if( ++nopt >= argc ) Syntax("-input needs a name") ;
00883          if( strcmp(input_filename, "") )
00884             Syntax ("Cannot have two input names!" ) ;
00885          strcpy (input_filename, argv[nopt] );
00886          /* open 3d data set */
00887          opt->dset = THD_open_one_dataset( argv[nopt] ) ;
00888          if( opt->dset == NULL ) 
00889          {
00890             sprintf (err_note, "Unable to open 3d+time data file: %s", argv[nopt]);
00891             Syntax (err_note); 
00892          }
00893          continue ;
00894       }
00895 
00896       /* unidentified option */
00897       sprintf( err_note , "Unknown option %s" , argv[nopt] ) ;
00898       Syntax(err_note) ;
00899       
00900    } while( ++nopt < argc ) ;  /* loop until all args are read, or we break */
00901 
00902 
00903    /* --- check for valid user inputs --- */
00904    if (!strcmp(input_filename,""))  Syntax ("No input file name was given.");
00905    if( opt->idts->num == 0 ) Syntax("No ideal response vector was defined!") ; 
00906    if (!strcmp(opt->prefix_name,""))  Syntax ("No prefix name was specified.");
00907 
00908 
00909    /* --- Done! --- */
00910 
00911    return ;
00912 }

00998 {
00999    line_opt  opt ;               /* holds data constructed from command line */
01000    THD_3dim_dataset * new_dset;  /* functional data set to be calculated */
01001    Boolean ok;                   /* true if 3d write is successful */
01002    
01003 
01004    /*----- Identify software -----*/
01005    printf ("\n\n");
01006    printf ("Program: %s \n", PROGRAM_NAME);
01007    printf ("Author:  %s \n", PROGRAM_AUTHOR);
01008    printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
01009    printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
01010    printf ("\n");
01011 
01012 
01013    /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
01014 
01015    { int new_argc ; char ** new_argv ;
01016      addto_args( argc , argv , &new_argc , &new_argv ) ;
01017      if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
01018    }
01019 
01020 
01021    /* --- read command line --- */
01022    get_line_opt( argc , argv , &opt ) ;
01023 
01024    /* --- calculate functional image --- */
01025    new_dset = fim3d_fimmer_compute (opt.dset, opt.idts, opt.ortts, 
01026                                     opt.first_im, opt.prefix_name, 
01027                                     opt.max_percent);  /* 19 May 1997 */
01028 
01029    /*----- Record history of dataset -----*/
01030    tross_Copy_History( opt.dset , new_dset ) ;
01031    tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
01032    
01033    /* --- write 3d functional image data --- */
01034    ok = THD_write_3dim_dataset ("", opt.prefix_name, new_dset, TRUE);
01035    if (!ok)  Syntax ("Failure to write functional data set.");
01036    
01037    /* --- cleanup --- */
01038    THD_delete_3dim_dataset (new_dset, FALSE);
01039    
01040    return (0);
01041 }

00919 {
00920    static char *mesg[] = {
00921    "Program:   3dfim \n\n"
00922    "Purpose:   Calculate functional image from 3d+time data file. \n"
00923    "Usage:     3dfim  [-im1 num]  -input fname  -prefix name \n"
00924    "              -ideal fname  [-ideal fname] [-ort fname] \n"
00925    " " ,
00926    "options are:",
00927    "-im1 num        num   = index of first image to be used in time series ",
00928    "                        correlation; default is 1  ",
00929    " ",
00930    "-input fname    fname = filename of 3d + time data file for input",
00931    " " ,
00932    "-prefix name    name  = prefix of filename for saving functional data",
00933    " ",
00934    "-ideal fname    fname = filename of a time series to which the image data",
00935    "                        is to be correlated. ",
00936    " ",
00937    "-percent p      Calculate percentage change due to the ideal time series ",
00938    "                p     = maximum allowed percentage change from baseline ",
00939    "                        Note: values greater than p are set equal to p. ",
00940    " ",
00941    "-ort fname      fname = filename of a time series to which the image data",
00942    "                        is to be orthogonalized ",
00943    " ", 
00944    "            N.B.: It is possible to specify more than",
00945    "            one ideal time series file. Each one is separately correlated",
00946    "            with the image time series and the one most highly correlated",
00947    "            is selected for each pixel.  Multiple ideals are specified",
00948    "            using more than one '-ideal fname' option, or by using the",
00949    "            form '-ideal [ fname1 fname2 ... ]' -- this latter method",
00950    "            allows the use of wildcarded ideal filenames.",
00951    "            The '[' character that indicates the start of a group of",
00952    "            ideals can actually be any ONE of these: " OPENERS ,
00953    "            and the ']' that ends the group can be:  " CLOSERS ,
00954    " ",
00955    "            [Format of ideal time series files:",
00956    "            ASCII; one number per line;",
00957    "            Same number of lines as images in the time series;",
00958    "            Value over 33333 --> don't use this image in the analysis]",
00959    " ",
00960    "            N.B.: It is also possible to specify more than",
00961    "            one ort time series file.  The image time series is  ",
00962    "            orthogonalized to each ort time series.  Multiple orts are ",
00963    "            specified by using more than one '-ort fname' option, ",
00964    "            or by using the form '-ort [ fname1 fname2 ... ]'.  This ",
00965    "            latter method allows the use of wildcarded ort filenames.",
00966    "            The '[' character that indicates the start of a group of",
00967    "            ideals can actually be any ONE of these: " OPENERS ,
00968    "            and the ']' that ends the group can be:  " CLOSERS ,
00969    " ",
00970    "            [Format of ort time series files:",
00971    "            ASCII; one number per line;",
00972    "            At least same number of lines as images in the time series]",
00973    " ",
00974    " ",
00975    } ;
00976 
00977    int nsize , ii ;
00978 
00979    if( note != NULL && (nsize=strlen(note)) > 0 ){
00980       fprintf(stderr,"\n") ;
00981       for(ii=0;ii<nsize+9;ii++) fprintf(stderr,"*") ;
00982       fprintf(stderr,"\a\n Error: %s\n", note ) ;
00983       for(ii=0;ii<nsize+9;ii++) fprintf(stderr,"*") ;
00984       fprintf(stderr,"\nTry 3dfim -help for instructions\a\n") ;
00985       exit(-1) ;
00986    }
00987 
00988    for( ii=0 ; ii < sizeof(mesg)/sizeof(char *) ; ii++ ){
00989       printf( " %s\n" , mesg[ii] );
00990    }
00991    exit(0) ;
00992 }