00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   This program calculates a functional image from an ANFI 3D+time data set, 
00010   by calculating the correlation between the image time series and one or
00011   more reference ("ideal") time series.
00012   
00013   File:     3dfim.c
00014   Date:     06 September 1996
00015 
00016   Incorporated "ort" time series into 3dfim program.
00017   BDW       12 December 1996
00018  
00019   Added -percent option for calculating relative effect of reference waveform
00020   upon observed time series.
00021   BDW       19 May 1997
00022 
00023   Corrected reference to "ort" time series data structure.
00024   BDW       05 Sept 1997
00025 
00026   Print a more explicit error message when ideal or "ort" time series are not 
00027   of sufficient length.
00028   BDW       14 January 1998
00029 
00030   Mod:     Added changes for incorporating History notes.
00031   Date:    10 September 1999
00032 
00033   Mod:     Added call to AFNI_logger.
00034   Date:    15 August 2001
00035 
00036 */
00037 
00038 /*---------------------------------------------------------------------------*/
00039 
00040 #define PROGRAM_NAME "3dfim"                         /* name of this program */
00041 #define PROGRAM_AUTHOR "R. W. Cox and B. D. Ward"          /* program author */
00042 #define PROGRAM_INITIAL "06 Sept 1996"    /* date of initial program release */
00043 #define PROGRAM_LATEST  "15 August 2001"  /* date of latest program revision */
00044 
00045 /*---------------------------------------------------------------------------*/
00046 
00047 #define SO_BIG 33333
00048 
00049 #include "afni.h"
00050 #include "mrilib.h"
00051 #include "ts.h"
00052 #include "afni_pcor.c"
00053 
00054 /*-------------------------------------------------------------------
00055    Lots of CPU work in here!
00056 ---------------------------------------------------------------------*/
00057 
00058 #undef  FIM_THR
00059 #define FIM_THR  0.0999
00060 
00061 
00062 THD_3dim_dataset * fim3d_fimmer_compute ( THD_3dim_dataset * dset_time ,
00063    time_series_array * ref_ts , time_series_array * ort_ts , 
00064    int itbot, char * new_prefix, 
00065    float max_percent        /* 19 May 1997 */ ) 
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 }
00694 
00695 /*--------------------------------------------------------------------------*/
00696 
00697 /** structure type to hold results from scanning the command line options **/
00698 
00699 typedef struct 
00700 {
00701    char prefix_name[THD_MAX_NAME];    /* filename root for output */
00702    THD_3dim_dataset * dset;           /* input 3d+time data set */
00703    time_series_array * idts;          /* array of ideal time series */
00704    time_series_array * ortts;         /* array of ort time series */ 
00705    int first_im;                      /* first time series point to be used */ 
00706    float max_percent;                 /* max. percentage change due to ideal
00707                                          time series,  19 May 1997 */
00708 } line_opt ;
00709 
00710 
00711 /*** internal prototype ***/
00712 
00713 void Syntax( char * ) ;
00714 
00715 /*--------------------------------------------------------------------------*/
00716 
00717 /*** read and interpret command line arguments ***/
00718 
00719 
00720 #define OPENERS "[{/%"
00721 #define CLOSERS "]}/%"
00722 
00723 #define IS_OPENER(sss) (strlen((sss))==1 && strstr(OPENERS,(sss))!=NULL)
00724 #define IS_CLOSER(sss) (strlen((sss))==1 && strstr(CLOSERS,(sss))!=NULL)
00725 
00726 
00727 void get_line_opt( int argc , char *argv[] , line_opt *opt )
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 }
00913 
00914 /*----------------------------------------------------------------------------*/
00915 
00916 /* --- give some help and exit --- */
00917 
00918 void Syntax( char *note )
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 }
00993 
00994 
00995 /*----------------------------------------------------------------------------*/
00996 
00997 int main( int argc , char *argv[] )
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 }
01042 
01043 /*---------------------------------------------------------------------------*/